[R-sig-ME] Population fit with glm works fine: totally off with glmer

[Dieter Menne]

Dear Mixed-Modelers,


I have data from subjects asked repeatedly if they felt satiated (0/1) with
different volumes of meal taken. There are also two treatment groups A and
B.

When I fit the data with glm, discarding subject information, the results
are fine and plotted as first plot below.

When I try the equivalent model with lmer, the results look fine at the
first site, but when plotted the results are totally nonsense. I do not see
any indication why the fit went wrong.

Any help? What did I do wrong? Is there a way to see that this result is
bogus without plotting the prediction (or looking at the coeffs, which is
already evidence enough).

Code and full data from Internet below the results.

Dieter


--------------------------
glm(formula = Satiated ~ MealVol * Group, family = "binomial",   data =
data)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-1.954  -0.757  -0.612   0.991   1.950  

Coefficients:
                    Estimate Std. Error z value Pr(>|z|)  
(Intercept)         -1.84661    1.03213   -1.79    0.074 .
MealVol              0.00738    0.00344    2.15    0.032 *
GroupTreatB         -0.51958    1.49443   -0.35    0.728  
MealVol:GroupTreatB -0.00385    0.00505   -0.76    0.446  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 319.17  on 236  degrees of freedom
Residual deviance: 275.71  on 233  degrees of freedom
AIC: 283.7

------------------------------
Generalized linear mixed model fit by the Laplace approximation Formula: 
Satiated ~ MealVol * Group + (1 | Subject) 
   Data: data 
 AIC BIC logLik deviance
 145 162  -67.6      135
Random effects:
 Groups  Name        Variance Std.Dev.
 Subject (Intercept) 232      15.2    
Number of obs: 237, groups: Subject, 31

Fixed effects:
                     Estimate Std. Error z value Pr(>|z|)   
(Intercept)           2.73101    5.61885    0.49   0.6269   
MealVol               0.02107    0.00944    2.23   0.0256 * 
GroupTreatB         -38.22735   11.66904   -3.28   0.0011 **
MealVol:GroupTreatB   0.05210    0.03004    1.73   0.0829 . 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Correlation of Fixed Effects:
            (Intr) MealVl GrpTrB
MealVol     -0.527              
GroupTreatB -0.482  0.254       
MlVl:GrpTrB  0.165 -0.314 -0.743
>
--------------------------------------------


## Code here; data from Internet
library(lme4)
library(lattice)
options(digits=2)
datafile = "Fullness.rdata"

# Get data (8kB) from Internet or from local cache
if (file.exists(datafile)) 
  load(datafile) else
{
  data= read.csv("http://www.menne-biomed.de/uni/fullness.csv")
#  data= read.csv("fullness.csv")
  data$Subject=as.factor(data$Subject)
  save(data,file=datafile)
}

table(data$Group,data$Satiated)
str(data)

# Brute force predict and plot (lmer has no predict)
plotPrediction = function(fit){
  cf  = coef(fit)[,1]
  sv = seq(50,400,by=10)
  
  score = exp(cf[1]+sv*cf[2])
  pA = score/(1+score)
  score = exp(cf[1]+cf[3]+sv*(cf[2]+cf[4]))
  pV = score/(1+score)

  pd = data.frame(
    group=rep(c("TreatA","TreatB"), each=length(pV)), 
    vol=rep(sv,2),p=c(pA,pV))
  
  xyplot(p~vol,groups=group,data=pd,ylim=c(0,1),type="l",lwd=2,
             auto.key=list(points=FALSE,lines=TRUE,columns=2),
             ylab="p(satiated)", xlab="volume",aspect=1.4)
  
}

## Using a global glm fit, discarding subject correlation
# The results are fine
(fitglm= summary(glm(Satiated~MealVol*Group,family="binomial",data=data)))
p1 = plotPrediction(fitglm)

# Take into account Subjects
# The results are totally off
(fitlmer=   summary(lmer(Satiated~MealVol*Group+(1|Subject), 
                         family=binomial, data=data)))
p2 = plotPrediction(fitlmer)
print(p1,split=c(1,1,2,1),more=TRUE)
print(p2,split=c(2,1,2,1))