8
$\begingroup$

I'm using the lme4 package in R to do some logistic mixed-effects modeling.
My understanding was that sum of each random effects should be zero.

When I make toy linear mixed-models using lmer, the random effects are usually < $10^{-10}$ confirming my belief that the colSums(ranef(model)$groups) ~ 0 But in toy binomial models (and in models of my real binomial data) some of the random effect sum to ~0.9.

Should I be concerned? How do I interpret this?

Here is a linear toy example

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
}

In this case the colMeans always come out $<10^{-6}$

Here is a binomial toy example (I would share my actual data, but it is being submitted for publication and I am not sure what the journal policy is on posting beforehand):


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 }

Now the colMeans sometimes come out above 0.3, and definitely higher, on average than the linear example.

| cite | improve this question | | | | |
$\endgroup$
  • 3
    $\begingroup$ Can you include code to reproduce those toy examples here? It would help in exploring this interesting behavior. $\endgroup$ – Aaron left Stack Overflow Feb 24 '13 at 5:14
  • $\begingroup$ I have seen this same behaviour with my experiments as well. In gaussian case there is a sum to zero constraint, but in non-gausian cases not. I'm not sure if sum-to-zero is necessary condition, as long as the expected value of random effects is zero. It might be helpful in some cases, and apparently it is easy to code in gaussian case so it's there... Hopefully somebody with better understanding chimes in. $\endgroup$ – Jouni Mar 18 '14 at 5:17
3
+50
$\begingroup$

Since @Hemmo's code got slightly mangled in the "Bounty" box, I'm adding this reformatted version as "community wiki". If this is not an appropriate use of the wiki, I apologize in advance. Feel free to remove it.

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")
) 
$\endgroup$
  • 1
    $\begingroup$ Well, it seems that question still didn't receive enough attention. My own conclusion is that there really isn't anything wrong here, there really is no sum-to-zero condition, but it just happens in gaussian cases where everything is linear. Expectation of random effects must be 0 is the real assumption, not that the actual sums of estimated effects are zero. I'll have to award somebody, so you're welcome. :) $\endgroup$ – Jouni Mar 26 '14 at 6:07
  • 2
    $\begingroup$ @Hemmo Yikes, now I feel like I should actually contribute something. You're right that nothing is actually wrong. The short answer (which I was hoping to write up but didn't find the time to), is that the mean will be zero iff the likelihood surface is Gaussian. Informally, we can "prove" this by noting that Gaussian errors * Gaussian random effects leads to another Gaussian. When you have a glmm with a non-Gaussian error function (e.g. Poisson in your case), then the likelihood surface can become non-Gaussian, and all bets are off. Hope this helps. $\endgroup$ – David J. Harris Mar 26 '14 at 7:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.