\documentclass{article}
% ----------------------------------------------------------------
\usepackage[brazil,brazilian]{babel}
\usepackage[latin1]{inputenc}
\usepackage{graphicx,Rd}
\usepackage[a4paper, text={16cm, 23.7cm}, top=2cm, left=2.5cm,
right=2.5cm, includeheadfoot, headheight=0cm]{geometry}
\title{Infer\^encia bayesiana para modelos com efeitos aleat\'orios em \R}
\author{Elias, Edson, Ana, PJ}
\begin{document}
\maketitle
\section*{Resumo}
Neste trabalho estudamos a implementa\c{c}~ao em \R{}
da infer\^encia em modelos de efeitos aleat\'orios.
A infer\^encia \'e sempre baseada na fun\c{c}\~ao de
verossimilhan\c{c}a dos dados. Fazemos estima\c{c}\~ao
pontual e intervalar por m\'axima verossimilhan\c{c}a
e bayesiana. No primeiro caso utizamos a fun\c{c}\~ao
\code{optim()} para encontrar as estimativas e no
segundo, utilizamos a fun\c{c}\~ao \code{MCMCmetrop1R()}
do pacote \pkg{MCMCpack} para simular amostras da
distribui\c{c}\~ao \`a \textit{posteriori}.
\section{Experimentos em blocos completos casualizados}
Vamos considerar nesta se\c{c}\~ao o caso em que temos um
experimento em blocos completos casualizados. Neste caso
temos que a observa\c{c}~ao $y_{i,j}$ da vari\'avel
resposta $Y$ \'e referente \`a repeti\c{c}\~ao $j$ no
bloco $i$, que $i=1,2,...,I$ e $j=1,2,...,n_i$.
Podemos dizer tamb\'em que temos $I$ indiv\'iduos sendo
observados e em cada um foram feitas $n_i$ observa\c{c}\~oes.
\subsection{Especifica\c{c}\~ao}
Vamos considerar o efeito aleat\'orio $U$ de bloco tem
distribui\c{c}\~ao Normal:
$$U \sim Normal(\mu, \tau^2)\;,$$
em que $\mu$ \'e a m\'edia do efeito de bloco e $\tau^2$
a vari\^ancia. Tamb\'em consideramos que a
$$Y_{ij}\sim Normal(U_j, \sigma^2)\;,$$
A distribui\c{c}\~ao da vari\'avel resposta condicionada
no efeito de bloco \'e obtida integrando em rela\c{c}~ao
ao efeito aleat\'orio
$$[Y] = \int [Y, U] dU$$
Al\'em disso, como os blocos s\~ao independentes, podemos
fatorar a distribui\c{c}\~ao de $Y$ em produtos de $k$
distribui\c{c}\~oes condicionais
$$[Y] = [Y_{1}][Y_{2}]...[Y_{k}]\;,$$
e a distribui\c{c}\~ao da resposta em cada bloco \'e
$$[Y_{j}] = \int [Y_{j}, u_j] du_j$$
que \'e igual \`a
$$[Y_{j}] = \int [Y_{j}|u_j][u_j] du_j$$
Definindo os par\^ametros em \R:
<>=
I <- 10 ## number of blocks
n.sample <- 20 ## number of samples in a block
mu <- 10 ## mean of data
sigma <- 2 ## sd of data
tau <- 0.5 ## sd of the random effect
@
Simulando uma amostra:
<>=
set.seed(1)
block <- rep(1:I, each=n.sample)
rand.eff <- rep(rnorm(I, mu, tau), each=n.sample)
y <- rnorm(I*n.sample, rand.eff, sigma)
@
\subsection{Fun\c{c}\~ao de verossimilhan\c{c}a}
A distribui\c{c}\~ao conjunta $[Y,U]$ \'e um produto de
duas densidades normais. Podemoas escrever a verossimilhan\c{c}a
desse modelo e obter as estimativas de m\'axima verossimilhan\c{c}a
em \R{} utilizando a fun\c{c}\~ao \code{optim()}.
Inicialmente definimos a sub-fun\c{c}\~ao conjunta $[Y,U]$
na escala logaritma:
<>=
dconj <- function(ran, obs, mu, s.ran, s.obs)
dnorm(obs, ran, s.obs)*dnorm(ran, mu, s.ran)
@
Vamos desenhar o gr\'afico da verossimilhan\c{c}a em
fun\c{c}\~ao do efeito aleat\'orio para um valor de
$\theta$ e alguns cada valor de $y$:
A contribui\c{c}\~ao de cada observa\c{c}\~ao para a
verossimilhan\c{c}a \'e a \'area abaixo sua respectiva curva.
Fazendo para as cinco primeiras observa\c{c}\~oes:
<>=
plot(function(x,...)
dconj(x,obs=y[1], mu=10,
s.ran=.5, s.obs=2), 8, 12, ylim=c(0,.16),
xlab="Efeito aleatorio",
ylab="Contribuicao para a verossimilhanca")
for (i in 2:5)
plot(function(x,...)
dconj(x,obs=y[i], mu=10,
s.ran=.5, s.obs=2),8,12, col=i, add=T)
legend(8, .16,
paste('y = ', format(y[1:5], dig=2), sep=''),
col=1:5, lty=1)
@
\begin{figure}\centering
\caption{Contribui\c{c}\~ao individual na verossimilhan\c{c}a}
<>=
<>
@
\end{figure}
A fun\c{c}\~ao de verossimilhan\c{c}a de $Y$ \'e obtida
integrando em rela\c{c}\~ao a $U$. Para isso, definimos
uma malha de valores razo\'aveis para o efeito aleta\'orio,
avaliamos a fun\c{c}\~ao nesses valores.
<>=
lvero <- function(theta, y, id, ran)
sum(log(colSums(outer(ran, y, dconj, mu=theta[1],
s.ran=theta[2],
s.obs=theta[3]))*diff(ran[1:2])))
@
Encontrando estimativas de m\'axima verossimilhan\c{c}a,
utilizando algoritmo de Nelder-Mead:
<<>>=
gran <- seq(5, 15, len=100)
opt1 <- optim(c(10,.5,2), lvero, y=y, id=block, ran=gran,
hessian=TRUE, control=list(fnscale=-1))
opt1$conv
opt1$par
sum(exp(opt1$par[2:3])^2)
var(y)
sqrt(diag(solve(-opt1$hess)))
@
Utilizando algoritmo de BFGS:
<<>>=
opt1b <- optim(c(10,.5,2), lvero, y=y, id=block, ran=gran,
hessian=TRUE, method='B',
control=list(fnscale=-1))
opt1b$conv
opt1b$par
sum(exp(opt1b$par[2:3])^2)
var(y)
sqrt(diag(solve(-opt1b$hess)))
@
A verossimilhan\c{a}a pode ser escrita de outra maneira.
Podemos considerar as observa\c{c}\~oes $y_i$ do
indiv\'iduo $i$ tem distribui\c{c}\~ao Normal $n_i$
variada, com vetor de m\'edias $u_i$ e matriz de covari\^ancia
com os elementos da diagonal igual a $\sigma^2+\tau^2$ e
os elementos fora da diagonal igual a $\tau^2$.
<<>>=
require(mvtnorm)
lvero2 <- function(par,y,id) {
n <- length(y)
ni <- table(id)
sum(sapply(unique(id), function(i) {
mc <- matrix(par[2]^2,ni[i],ni[i])+diag(ni[i])*par[3]^2
dmvnorm(y[id==i], rep(par[1],ni[i]), mc, log=TRUE)
}))
}
opt2 <- optim(c(10,.5,2), lvero2, y=y, id=block, hessian=TRUE,
control=list(fnscale=-1))
opt2$conv
opt2$par
sum(opt2$par[2:3]^2)
var(y)
sqrt(diag(solve(-opt2$hess)))
@
Utilizando a fun\c{c}\~ao \code{glmmPQL()} do pacote \pkg{MASS}.
<>=
require(MASS)
summary(aj.glmm <- glmmPQL(y~1,~1|block,gaussian))
@
E notamos que tem resultados similares \'a abordagem por blocos.
\section{Infer\^encia bayesiana}
Agora vamos tratar $\theta$ como vari\'avel aleat\'oria.
Vamos supor que \'a \textit{priori} $\mu$, $\tau^2$ e
$\sigma^2$ s\~ao independentes.
\'E razo\'avel considerar uma \textit{priori} Normal para
$\mu$ e \textit{prioris} Gamma para $\tau^2$ e $\sigma^2$.
Sem obter express\~oes ana\'iticas para a distribui\c{c}\~oes
condicionais \'a \textit{posteriori} de $\mu$, $\tau^2$ e
$\sigma^2$, vamos simular amostras dessas distribui\c{c}\~oes.
Utilizaremos o algor\'itmo de Metropolis.
O Metropolis random walk esta implementado em C++ no pacote
\pkg{MCMCpack} e pode ser utilizado via \R{} com a fun\c{c}\~ao
\code{MCMCmetrop1R()}. Nesse algoritmo, os valores
propostos para todos os par\^ametros s\~ao simulados de uma
distribui\c{c}\~ao Normal multivariada. O vetor de m\'edias
utilizado a cada itera\c{c}\~ao \'e o valor atual dos par\^ametros.
A matriz de vari\^ancia pode ser informada pelo usu\'ario ou ent\~ao
\'e obtida do hessiano num\'erico calculado internamente na
fun\c{c}\~ao. Neste caso, pode ser passado um vetor do tamanho do
n\'umero de par\^ametros para ajustar a vari\^ancia de acordo
com a taxa de aceita\c{c}\~ao.
Os par\^ametros de vari\^ancia s\~ao positivos e podemos simular
de uma distribui\c{c}\~ao log-normal fazendo uma altera\c{c}\~ao
na func\c{a}\~ao de verossimilhan\c{c}a:
<>=
lvero2p <- function(par,y,id) {
n <- length(y)
ni <- table(id)
sum(sapply(unique(id), function(i) {
mc <- matrix(exp(par[2])^2,ni[i],ni[i])+diag(ni[i])*exp(par[3])^2
dmvnorm(y[id==i], rep(par[1],ni[i]), mc, log=TRUE)
}))
}
@
E podemos entao utilizar esta fun\c{c}\~ao para simular da
distribui\c{c}\~ao \'a \textit{posteriori}.
N\~ao dispomos de muito conhecimento \'a \textit{priori} acerca
de $\theta$, apenas sabemos que $\mu \in \Re$, $\tau^2$ e
$\sigma^2$ s\^ao estritamente positivos. Escolhemos as
seguintes distribui\c{c}\~oes \'a \textit{priori}:
$$\mu \sim Normal(0, 10000)\;,$$
$$log(\tau^2) \sim Normal(0, 9)\; e$$
$$log(\sigma^2) \sim Normal(0, 9)\;.$$
Vizualizando em \R{}:
\begin{figure}\centering
\caption{Distribui\c{c}\~oes \'a \textit{priori}.
O terceiro gr\'afico \'e um close so segundo}
<>=
par(mfrow=c(1,3), mar=c(3,3,1,1), mgp=c(1.7,.5,0))
plot(function(x) dnorm(x,0,100), -300, 300,
xlab=expression(mu), ylab=expression(p(mu)))
plot(function(x) dnorm(log(x),0,3), 0, 300,
xlab=expression(tau), ylab=expression(p(tau)))
plot(function(x) dnorm(log(x),0,3), 0, 5,
xlab=expression(tau), ylab=expression(p(tau)))
@
\end{figure}
O logaritmo da distribui\c{c}\~ao \'a priori em \R{}
pode ser escrito como:
<<>>=
p.theta <- function(mu, tau2, sigma2)
dnorm(mu, 0, 100, TRUE) +
dnorm(log(tau2), 0, 1, TRUE) +
dnorm(log(sigma2), 0, 3, TRUE)
@
O logaritmo da distribui\c{c}\~ao \'a posteriori ent\~ao \'e:
<<>>=
dpost <- function(theta, y, id)
p.theta(theta[1], exp(theta[2])^2, exp(theta[3])^2) +
lvero2p(theta, y, id)
@
Vizualiza\c{c}\~ao das distribui\c{c}\~oes \'a
\textit{posteriori} condicionais.
<>=
par(mfrow=c(1,3), mar=c(3,3,1,1), mgp=c(1.7,.5,0))
plot(function(x)
sapply(x,function(a)
dpost(c(a,log(.5),log(2)), y, block)), 9, 11, n=31,
xlab=expression(mu), ylab=expression(L(mu/y,tau,sigma)))
plot(function(x)
sapply(x,function(a)
dpost(c(10,a,log(2)), y, block)), -2, 0, n=31,
xlab=expression(exp(tau)), ylab=expression(L(exp(tau)/y,mu,sigma)))
plot(function(x)
sapply(x,function(a)
dpost(c(10,log(.5),a), y, block)), 0, 2, n=31,
xlab=expression(exp(sigma)), ylab=expression(L(exp(sigma)/y,mu,tau)))
@
\begin{figure}\centering
\caption{Vizualizando a distribui\c{c}\~oes
\'a posteriori condicionais}
<>=
<>
@
\end{figure}
Uma id\'eia para a matriz de vari\^ancia da
distribui\c{c}\~ao da proposta \'e considerar a
a matriz de vari\^ancia estimada pela inversa da
matriz hessiana, que pode ser obtida numericamente
com a fun\c{c}\~ao \code{optim()}.
<<>>=
opt2b <- optim(c(10,-1,1), dpost, y=y, id=block, hessian=TRUE,
control=list(fnscale=-1))
str(opt2b)
round(mcov <- solve(-opt2b$hess),4)
round(v.metrop <- diag(c(1.5,1.7,1))%*%mcov%*%diag(c(1.5,1.7,1)),4)
@
Neste caso, utilizar o argumento \code{V=v.metrop}, \'e o
mesmo que utilizer \code{V=NULL} e \code{tune=c(1.5,1.5,0.7)}.
Simulando da distribui\c{c}\~ao \'a posteriori:
<>=
library(MCMCpack)
post <- MCMCmetrop1R(dpost, theta.init=c(10,log(.5), log(2)),
100, 3000, 10, verbose=500,
V=v.metrop, y=y, id=block)
@
Visualizando as amostras da distribui\c{c}\~ao \'a \textit{posteriori}.
\begin{figure}
\centering
\caption{Tra\c{c}o e densidades das amostras da
distribui\c{c}\~ao \'a posteriori}
<>=
postb <- post
postb[,2] <- exp(post[,2])
postb[,3] <- exp(post[,3])
par(mar=c(3,3,1,1), mgp=c(1.7,.5,0))
plot(postb)
@
\end{figure}
Calculando quantidades de interesse:
<<>>=
summary(postb)
@
HPD intervalo
<<>>=
ic <- HPDinterval(postb)
ic
apply(postb, 2, median)
@
\subsection{Utilizando o princ\'ipio de aumento de dados}
Uma alternativa que pode ser utilizada em lugar de integrar
os efeitos aleat\'orios, \'e utilizar o princ\'ipio de aumento
de dados e simul\'a-los. Neste caso temos que escrever uma
fun\c{c}\~ao para $[y/u,\tau^2,\sigma^2]$, que recebe como
argumentos o vetor $u$, $log(\tau)$ e $log(\sigma)$.
Esta fun\c{c}\~ao deve retornar
$[\theta][u/tau][y/u,sigma]$.
<>=
post.aug <- function(pars, y, id) {
nb <- length(unique(id))
mu <- mean(pars[1:nb])
p.theta(mu, exp(pars[nb+1])^2, exp(pars[nb+2])^2) +
sum(dnorm(pars[1:nb], mu, exp(pars[nb+1]), log=TRUE)) +
sum(dnorm(y, pars[id], exp(pars[nb+2]), log=TRUE))
}
opt3 <- optim(c(rep(10,10),1,1), post.aug, y=y, id=block,
hess=T, control=list(fnscale=-1))
opt3$conv
c(mean(opt3$par[1:10]), exp(opt3$par[11:12]))
@
Matriz de variancia para algoritmo de Metropolis:
<<>>=
v3.metrop <- diag(c(rep(.6,10),1,1))%*%solve(-opt3$hess)%*%
diag(c(rep(.6,10),1,1))
100*round(v3.metrop,3)
@
<>=
post3 <- MCMCmetrop1R(post.aug, c(rep(10,10), exp(.5), exp(2)),
1000, 10000, 10, verbose=1000,
V=v3.metrop, y=y, id=block)
post3b <- post3[,10:12]
post2b <- colMeans(post3[,1:10])
post3b[,2] <- exp(post3[,11])
post3b[,3] <- exp(post3[,12])
@
\begin{figure}
\centering
\caption{Efeitos aleat\'orios simulados e par\^ametros de
vari\^ancia}
<>=
par(mfrow=c(3,4), mar=c(2,2,1,.5), mgp=c(1.2,.2,0))
for (i in 1:ncol(post3))
plot.ts(post3[,i])
@
\end{figure}
\begin{figure}\centering
\caption{Tra\c{o} e densidade das amostras \'a posteriori
dos par\^ametros de interesse}
<>=
par(mar=c(3,3,1,1), mgp=c(1.7,.5,0))
plot(post3b)
@
\end{figure}
Estat\'isticas dos par\^ametros de interesse:
<>=
summary(post3b)
HPDinterval(post3b)
apply(post3b, 2, median)
@
\end{document}