Verossimilhança para um modelo de efeitos aleatórios

September 25, 2011
By

This post was kindly contributed by Ridículas - go there to comment and to read the full post.

Contornos perfilhados de deviance para regiões de 0,9; 0,95; e 0,99 de confiança. O ponto + no gráfico representa os valores paramétricos usados para simulação dos dados. Imperfeições nos contornos apontam problemas numéricos associados ao passo de integração.

Modelos com termos de efeito aleatórios estão sendo cada vez mais usados na modelagem de dados. Nestes modelos, parte da variabilidade observada na resposta pode ser explicada por um conjunto de fatores de efeito fixo, outra parte pode estar associada a realização de variáveis latentes (ou não observáveis)  a qual é atribuída o efeito aleatório.

A decisão sobre o que é fixo ou aleatório nem sempre é bem clara e depende dos objetivos da pesquisa/inferência. Um caso bem típico é a questão dos blocos nos experimentos agronômicos. Uma maneira simples de pensar é questionar: os blocos presentes no meu experimento foram sorteados de um universo de blocos possíveis? Se sim, considere eles como de efeito aleatório. Há casos em que eles não são de efeito aleatório como, por exemplo, quando o bloco identifica o nível topográfico das parcelas, a classe de peso dos animais. Há casos que é genuinamente de efeito aleatório como, por exemplo, o lote da sementes em teste. Mas como eu ressaltei, a decisão sobre fixo/aleatório depende dos objetivos e difere de pessoa para pessoa.

Quando você declara um fator com efeito aleatório, você assume que os efeitos que ele exerce sobre a resposta são realizações de uma distribuição de probabilidades. Imagine um experimento para verificar a germinação de sementes onde temos 10 lotes de sementes e 4 repetições por lote. Os lotes estudados foram sorteados de uma grande lista de lotes. Ao realizar o experimento, verifica-se variabilidade entre repetições de um mesmo lote e entre lotes. Se lote fosse de efeito fixo, teríamos que estimar 10 parâmetros, um para cada lote. Se lote é de efeito aleatório temos apenas que estimar os parâmetros da distribuição de probabilidades associada a realização dos efeitos. Normalmente é usada a distribuição normal, e nesse caso é estimado a variância dos efeitos aleatórios.

Os efeitos aleatórios não são parâmetros, são realizações de uma distribuição de probabilidades. A verossimilhança para um modelo com efeitos aleatórios é baseado nessa afirmativa. A distribuição marginal da resposta é a integral no domínio dos efeitos aleatórios do produto entre a distribuição condicional da resposta ao valor realizado dos efeitos e a densidade dos efeitos aleatórios, ou seja

f(y) = \int f(y|b) \times f(b)\,\, \text{d}b

em que y é a resposta observada e b o efeito aleatório. Para o nosso experimento com germinação de sementes, temos uma resposta binomial (sementes germinadas em ensaios de 50 sementes), 10 níveis de lote e 4 repetições por lote. A verossimilhança seria

f(\sigma^2_l;y) = \prod_i \int \prod_j \text{Binomial}(y_{ij}|b_i) \times \text{Normal}(b_i, \sigma^2_l)\,\, \text{d}b_i.

No script abaixo eu fiz a simulação de dados para esse modelo, a função de verossimilhança e a otimização usando a optim() bem como a obtenção das regiões de confiança baseados na deviance perfilhada. As imperfeições observadas nos contornos de deviance são causadas por problemas numéricos na etapa de integração. Para checagem fiz o ajuste pela função lme4::glmer(), porém, esta função usa aproximação de Laplace (por isso estimativas levemente diferentes) e possivelmente usa uma função de verossimilhança sem as constantes padronizadoras (por isso máximo da verossimilhança tão diferente).

Para o caso de dados com distribuição Normal, o passo de integração é evitado porque a distribuição marginal possui forma fechada. Além do mais, diversos métodos numéricos podem ser empregados no passo de integração, como aproximação de Laplace, quadratura Gaussiana, integração Monte Carlo.  Outra coisa que é importante são as parametrizações usadas. Em geral, recomenda-se que os parâmetros tenham domínio nos reais. Para o caso da variância é usual otimizar o log do desvio-padrão.

Nos próximos posts vou abordar o perfil de verossimilhança para um parâmetro e a predição dos efeitos aleatórios. Até a próxima ridícula.

#-----------------------------------------------------------------------------
# dados gerados de um experimento de germinação de sementes

sigma <- 1 # desvio-padrão do efeito de lote
mu <- 0    # exp(mu)/(1+exp(mu)) = probabilidade de germinação
i <- 10    # níveis de lote
j <- 4     # número de repetições por lote
m <- 50    # número de sementes por ensaio
b <- rnorm(i, 0, sigma) # efeitos dos lotes
z <- gl(i, j)            # vetor de níveis de lote
y <- rbinom(i*j, size=m, prob=1/(1+exp(-mu-rep(b, e=j)))) # germinação observada

require(lattice)
xyplot(y~z, jitter.x=TRUE, col=z)

#-----------------------------------------------------------------------------
# função de log-verossimilhança

L <- function(par, y, z, m=50){
  ## par: vetor intercepto e parâmetro de variância (numeric)
  ## y: vetor valores observados para o número de sucessos (numeric)
  ## z: vetor de níveis de lote (factor)
  ## m: escalar número de sementes no ensaio (integer)
  ##-------------------------------------------------------
  ## lista em que cada slot são dos dados de um lote
  data <- split(y, z)
  ##-------------------------------------------------------
  ## verossimilhança para cada lote: f(y|b)*f(b)
  ## tem que função vetorial para usar na integrate()
  Li <- function(b, par, yi, m){
    sapply(b,
           function(bj){
             nu <- 1/(1+exp(-(par[1]+bj)))
             lj <- sum(dbinom(yi, prob=nu, size=m, log=TRUE))+
               dnorm(bj, mean=0, sd=exp(par[2]), log=TRUE)
             return(exp(lj))
           })
  }
  ##-------------------------------------------------------
  ## valor das i integrais
  int <-
    lapply(data,
           function(i){
             integrate(Li, -Inf, Inf, par=par,
                       yi=i, m=m)$value
           })
  ##-------------------------------------------------------
  ## valor da log-verossimilhança
  ll <- sum(log(unlist(int)))
  print(c(ll, par))
  return(ll)
  ##-------------------------------------------------------
}

#-----------------------------------------------------------------------------
# usa optim() para encontrar as estimativas de máxima verossimilhança

op <- optim(c(0,0), L, y=y, z=z, m=50,
            method="BFGS", control=list(fnscale=-1))
op$par   # estimativas de máxima verossimilhança para b0 e log(sigma)
op$value # valor do máximo da verossmilhança

#-----------------------------------------------------------------------------
# perfil de verossimilhança conjunto para b0 e log(sigma)

mu.grid <- seq(-2, 1, l=50)
lsigma.grid <- seq(-1, 1, l=50)
grid <- expand.grid(mu=mu.grid, lsigma=lsigma.grid)
grid$ll <- apply(grid, 1, L, y=y, z=z, m=m)

qch <- qchisq(c(0.9,0.95,0.99), df=2)
grid$dev <- -2*(grid$ll-op$value)

wireframe(dev~lsigma+mu, data=grid)

png("f022.png", w=500, h=500)
levelplot(dev~lsigma+mu, data=grid,
          xlab=expression(log(sigma)), ylab=expression(mu),
          panel=function(..., at, contour, region){
            panel.levelplot(..., at=at)#, at=94:195, contour=FALSE, region=TRUE)
            panel.contourplot(..., at=qch, contour=TRUE, region=FALSE)
            panel.abline(v=op$par[2], h=op$par[1], lty=3)
            panel.points(log(sigma), mu, pch=3, col=1)
          })
dev.off()

#-----------------------------------------------------------------------------
# usando a glmer()

require(lme4)

mm0 <- glmer(cbind(y, 50-y)~(1|z), family=binomial)
summary(mm0)
c(op$par[1], exp(op$par[2]))

#-----------------------------------------------------------------------------

Tags: , , , , ,

Comments are closed.