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

[Ben Bolker]

Dieter Menne  <dieter.menne at ...> writes:

> 
> 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


  Have just begun to investigate this, but the obvious (?)
hallmark of something going wrong is the extremely large
magnitude of the estimate, and the standard error, for GroupTreatB,
which suggests something Hauck-Donnerish and/or complete separation
that is induced once you add subjects to the treatment. (I now
see that you already mentioned looking at the coefficients
as a way to see that something had gone wrong.) The estimated among-subject
variance is also gigantic.  Some possible
solutions (not all necessarily feasible) are (1) bias-reduced
estimation a la brglm (probably not feasible in this context),
(2) a Bayesian solution with a bit more of a prior, via
MCMCglmm or blme (a package by Vincent Dorie: see below)

Looking at pictures makes the source of the problem
seem pretty clear:

library(ggplot2)
sdata <- transform(sdata,Subject=reorder(Subject,Satiated,mean))
g1 <- ggplot(sdata,aes(x=MealVol,y=Satiated,colour=Group))+geom_point()

g1 + geom_smooth(method="glm",family=binomial,formula=y~poly(x,2))

g1 + facet_wrap(~Subject)

 There are 13 individuals who are never satiated regardless of
meal volume (3 in Treatment A, 10 in Treatment B); there are 8 who
are always satiated (7 in treatment A, 1 in Treatment B).

  This (using the aforementioned Dorie package)

library(blme)
(fitbglmer <-   summary(g3 <- bglmer(Satiated~MealVol*Group+(1|Subject), 
                         family=binomial, data=sdata)))

## requires LATEST version of coefplot2 from r-forge:
##  packages won't be rebuilt until tomorrow, probably
library(coefplot2)
coefplot2(list(gf1,gf2,gf3),col=c(1,2,4))
coefplot2(list(gf1,gf2,gf3),xlim=c(-0.05,0.15))

  I was actually very encouraged, because the couple of times
I have tried bglmer before, on more complicated things, it has
broken ...