Quão perto de zero deve ser a soma dos efeitos aleatórios no GLMM (com lme4)

8

Estou usando o lme4pacote no R para fazer alguma modelagem logística de efeitos mistos.
Meu entendimento era que a soma de cada efeito aleatório deveria ser zero.

Quando eu uso modelos mistos lineares de brinquedo lmer, os efeitos aleatórios geralmente são < confirmando minha crença de que alguns dos modelos binomiais de brinquedos But (e nos modelos dos meus dados binomiais reais) alguns dos efeitos aleatórios somam ~ 0,9.1010colSums(ranef(model)$groups) ~ 0

Devo me preocupar? Como eu interpreto isso?

Aqui está um exemplo de brinquedo linear

toylin<-function(n=30,gn=10,doplot=FALSE){
 require(lme4)
 x=runif(n,0,1000)
 y1=matrix(0,gn,n)
 y2=y1
 for (gx in 1:gn)
 {
   y1[gx,]=2*x*(1+(gx-5.5)/10) + gx-5.5  + rnorm(n,sd=10)
   y2[gx,]=3*x*(1+(gx-5.5)/10) * runif(1,1,10)  + rnorm(n,sd=20)
 }
 c1=y1*0;
 c2=y2*0+1;
 y=c(t(y1[c(1:gn),]),t(y2[c(1:gn),]))
 g=rep(1:gn,each=n,times=2)
 x=rep(x,times=gn*2)
 c=c(c1,c2)
 df=data.frame(list(x=x,y=y,c=factor(c),g=factor(g)))
 (m=lmer(y~x*c + (x*c|g),data=df))
 if (doplot==TRUE)
  {require(lattice)
   df$fit=fitted(m)
   plot1=xyplot(fit ~ x|g,data=df,group=c,pch=19,cex=.1)
   plot2=xyplot(y ~ x|g,data=df,group=c)
   print(plot1+plot2)
  }
 print(colMeans(ranef(m)$g))
 m
}

Nesse caso, o colMeans sempre sai<106

Aqui está um exemplo de brinquedo binomial (eu compartilharia meus dados reais, mas eles estão sendo enviados para publicação e não tenho certeza de qual é a política do diário antes de postar):


toybin<-function(n=100,gn=4,doplot=FALSE){
  require(lme4)
x=runif(n,-16,16) y1=matrix(0,gn,n) y2=y1 for (gx in 1:gn) { com=runif(1,1,5) ucom=runif(1,1,5) y1[gx,]=tanh(x/(com+ucom) + rnorm(1)) > runif(x,-1,1) y2[gx,]=tanh(2*(x+2)/com + rnorm(1)) > runif(x,-1,1) } c1=y1*0; c2=y2*0+1; y=c(t(y1[c(1:gn),]),t(y2[c(1:gn),])) g=rep(1:gn,each=n,times=2) x=rep(x,times=gn*2) c=c(c1,c2) df=data.frame(list(x=x,y=y,c=factor(c),g=factor(g))) (m=lmer(y~x*c + (x*c|g),data=df,family=binomial)) if (doplot==TRUE) {require(lattice) df$fit=fitted(m) print(xyplot(fit ~ x|g,data=df,group=c,pch=19,cex=.1)) } print(colMeans(ranef(m)$g)) m }

Agora, o colMeans às vezes sai acima de 0,3 e, definitivamente, mais alto, em média do que o exemplo linear.

jerlich
fonte
3
Você pode incluir código para reproduzir esses exemplos de brinquedos aqui? Ajudaria a explorar esse comportamento interessante.
Aaron deixou Stack Overflow em
Vi esse mesmo comportamento com meus experimentos também. No caso gaussiano, existe uma soma à restrição zero, mas nos casos não gausianos, não. Não tenho certeza se a soma-zero é condição necessária, desde que o valor esperado dos efeitos aleatórios seja zero. Pode ser útil em alguns casos, e aparentemente é fácil codificar em caso gaussiano, por isso está lá ... Espero que alguém com melhor entendimento entre em contato.
Jouni

Respostas:

3

Como o código do @ Hemmo ficou um pouco distorcido na caixa "Bounty", estou adicionando esta versão reformatada como "wiki da comunidade". Se este não for um uso apropriado do wiki, peço desculpas antecipadamente. Sinta-se livre para removê-lo.

library(mvabund)
library(lme4) 
data(spider) 
Y <- as.matrix(spider$abund)
X <- spider$x 
X <- X[ ,c(1, 4, 5, 6)] 
X <- rbind(X, X, X, X, X, X, X, X, X, X, X, X) 
site <- rep(seq(1, 28), 12) 
dataspider <- data.frame(c(Y), X, site) 
names(dataspider) <- c("Y","soil.dry", "moss", "herb.layer", "reflection", "site") 
fit <- glmer(
  Y ~ soil.dry + moss + herb.layer + reflection + (1|site), 
  family = poisson(link = log), 
  data = dataspider,
  control = glmerControl(optimizer = "bobyqa")
) 
David J. Harris
fonte
1
Bem, parece que a pergunta ainda não recebeu atenção suficiente. Minha própria conclusão é que realmente não há nada errado aqui, realmente não existe uma condição de soma para zero, mas isso acontece apenas nos casos gaussianos em que tudo é linear. A expectativa de efeitos aleatórios deve ser 0 é a suposição real, não que as somas reais dos efeitos estimados sejam zero. Vou ter que premiar alguém, então você é bem-vindo. :)
Jouni
2
@ Hemmo Yikes, agora sinto que devo realmente contribuir com algo. Você está certo que nada está realmente errado. A resposta curta (que eu esperava escrever, mas não encontrou tempo para isso) é que a média será zero se a superfície da probabilidade for gaussiana. Informalmente, podemos "provar" isso observando que erros gaussianos * Os efeitos aleatórios gaussianos levam a outro gaussiano. Quando você tem um glmm com uma função de erro não gaussiana (por exemplo, Poisson no seu caso), a superfície de probabilidade pode se tornar não gaussiana e todas as apostas estão desativadas. Espero que isto ajude.
David J. Harris