R help - Nov 2009 - different results across versions for glmer/lmer with the quasi-poisson or quasi-binomial families: the lattest version might not be accurate...

Dear R-helpers,

this mail is intended to mention a rather trange result and generate potential
useful comments on it. I am not aware of another posts on this issue (
RSiteSearch("quasipoisson lmer version dispersion")).

MUsing the exemple in the reference of the lmer function (in lme4 library) and
turning it into a quasi-poisson or quasi-binomial analysis, we get different
results, especially for the dispersion parameter (section 1 in the following;
cf. bold passages in the following). For the quasi-poisson, we found even
stronger differences with another, bigger database on which we worked (residual
variance from approx. 0.7 to 1.2). We also tested the differences between
versions for plain likelihood families (poisson and binomial): the differences
were slighter, and were actually very small in the case of the poisson family.

My first reaction was that the later version should have been better. However,
the quasibinomial results of this version are rather strange. Therefore, I
simulated a bigger data set corresponding to Poisson and binomial mixed models
(section 2 in the following). I actually found out that the later version of
lme4 was suspect but not the old one. My temptative conclusion is therefore to
use the old versions of lme4 (here:  lme4 version 0.99875-9) when using
quasi-binomial and quasi-poisson methods.

Any comment/insight appreciated.

 
All the best,
 
Frédéric Gosselin
Engineer & Researcher (PhD) in Forest Ecology
Cemagref
Domaine des Barres
F-45290 Nogent sur Vernisson
France


############### 1- EXAMPLE FROM LMER HELP

############# here are the commands for the quasipoisson:
library(lme4)
gm1 <- lmer(incidence ~ period + (1 | herd), family = quasipoisson, data =
cbpp)
summary(gm1)


########## here is the result under R2.5.1 and Package lme4 version 0.99875-9:
Generalized linear mixed model fit using Laplace
Formula: incidence ~ period + (1 | herd)
   Data: cbpp
 Family: quasipoisson(log link)
   AIC   BIC logLik deviance
 112.2 122.3 -51.11    102.2
Random effects:
 Groups   Name        Variance Std.Dev.
 herd     (Intercept) 0.35085  0.59233
 Residual             1.40470  1.18520
number of obs: 56, groups: herd, 15

Fixed effects:
            Estimate Std. Error t value
(Intercept)   1.2812     0.2200   5.824
period2      -1.1240     0.3315  -3.391
period3      -1.3203     0.3579  -3.689
period4      -1.9477     0.4808  -4.051

Correlation of Fixed Effects:
        (Intr) perid2 perid3
period2 -0.339             
period3 -0.314  0.219      
period4 -0.233  0.163  0.151


#############" here is the result on R2.7.1 and Package lme4 version
0.999375-26

Generalized linear mixed model fit by the Laplace approximation
Formula: incidence ~ period + (1 | herd)
   Data: cbpp
   AIC   BIC logLik deviance
 114.2 126.4  -51.1    102.2
Random effects:
 Groups   Name        Variance Std.Dev.
 herd     (Intercept) 0.32421  0.56939
 Residual             1.29474  1.13786
Number of obs: 56, groups: herd, 15

Fixed effects:
            Estimate Std. Error t value
(Intercept)   1.2762     0.2115   6.035
period2      -1.1249     0.3187  -3.530
period3      -1.3190     0.3438  -3.837
period4      -1.9450     0.4615  -4.215

Correlation of Fixed Effects:
        (Intr) perid2 perid3
period2 -0.339             
period3 -0.314  0.219      
period4 -0.233  0.163  0.151



############# now the commands for the quasibinomial:
library(lme4)
toto<-as.double(cbpp$incidence>0)
gm2 <- lmer(toto ~ period + (1 | herd), family = quasibinomial, data = cbpp)
summary(gm2)
 
########## here is the result under R2.5.1 and Package lme4 version 0.99875-9:

Generalized linear mixed model fit using Laplace 
Formula: toto ~ period + (1 | herd) 
   Data: cbpp 
 Family: quasibinomial(logit link)
   AIC   BIC logLik deviance
 72.04 82.17 -31.02    62.04
Random effects:
 Groups   Name        Variance   Std.Dev.  
 herd     (Intercept) 5.0259e-10 2.2418e-05
 Residual             1.0052e+00 1.0026e+00
number of obs: 56, groups: herd, 15
 
Fixed effects:
            Estimate Std. Error t value
(Intercept)    2.662      1.048   2.540
period2       -2.078      1.188  -1.750
period3       -2.950      1.180  -2.501
period4       -3.135      1.194  -2.626
 
Correlation of Fixed Effects:
        (Intr) perid2 perid3
period2 -0.882              
period3 -0.888  0.784       
period4 -0.878  0.775  0.780

 
 
############# here is the result on R2.7.1 and Package lme4 version 0.999375-26

Generalized linear mixed model fit by the Laplace approximation 
Formula: toto ~ period + (1 | herd) 
   Data: cbpp 
   AIC  BIC logLik deviance
 74.04 86.2 -31.02    62.04
Random effects:
 Groups   Name        Variance Std.Dev.
 herd     (Intercept) 0.00000  0.00000 
 Residual             0.13522  0.36772 
Number of obs: 56, groups: herd, 15
 
Fixed effects:
            Estimate Std. Error t value
(Intercept)   2.6391     0.3806   6.933
period2      -2.0513     0.4324  -4.744
period3      -2.9267     0.4293  -6.817
period4      -3.1091     0.4345  -7.155
 
Correlation of Fixed Effects:
        (Intr) perid2 perid3
period2 -0.880              
period3 -0.887  0.780       
period4 -0.876  0.771  0.777

 
 
 
 
############### 2- SIMULATED EXAMPLE

############# here are the commands both for the quasibinomial and 
quasipoisson:


library(lme4)

set.seed(1)

period<-rnorm(1000)

herd<-rep(1:50,eac=20)

herd.effect<-rnorm(50)

toto<-rbinom(1000,1,exp(period+herd.effect[herd])/(1+exp(period+herd.effect[herd])))

gm2s <- lmer(toto ~ period + (1 | herd), family = quasibinomial)

summary(gm2s)

 

gm2sL <- lmer(toto ~ period + (1 | herd), family = binomial)

summary(gm2sL)

 

#now poisson with similar data:

#essai avec jeu de données simulé:

set.seed(2)

toto<-rpois(1000,exp(period+herd.effect[herd]))

gm1s <- lmer(toto ~ period + (1 | herd), family = quasipoisson)

summary(gm1s)

 

gm1sL <- lmer(toto ~ period + (1 | herd), family = poisson)

summary(gm1sL)

 

 

########## here is the quasi-poisson result under R2.5.1 and Package lme4
version 0.99875-9: (all is OK)


Generalized linear mixed model fit using Laplace 
Formula: toto ~ period + (1 | herd) 
 Family: quasipoisson(log link)
  AIC  BIC logLik deviance
 1151 1166 -572.7     1145
Random effects:
 Groups   Name        Variance Std.Dev.
 herd     (Intercept) 0.83605  0.91436 
 Residual             0.97980  0.98985 
number of obs: 1000, groups: herd, 50
 
Fixed effects:
            Estimate Std. Error t value
(Intercept)  0.02775    0.13419    0.21
period       0.97424    0.02291   42.52
 
Correlation of Fixed Effects:
       (Intr)
period -0.157

 
############# here is the quasipoisson result on R2.7.1 and Package lme4 version
0.999375-26: the variance estimates are suspect
Generalized linear mixed model fit by the Laplace approximation 
Formula: toto ~ period + (1 | herd) 
  AIC  BIC logLik deviance
 1153 1173 -572.7     1145
Random effects:
 Groups   Name        Variance Std.Dev.
 herd     (Intercept) 0.35305  0.59418 
 Residual             0.41370  0.64320 
Number of obs: 1000, groups: herd, 50
 
Fixed effects:
            Estimate Std. Error t value
(Intercept)  0.02757    0.08720    0.32
period       0.97423    0.01489   65.44
 
Correlation of Fixed Effects:
       (Intr)
period -0.157

 
 

########## here is the quasi-binomial result under R2.5.1 and Package lme4
version 0.99875-9: (all is OK)


Generalized linear mixed model fit using Laplace 
Formula: toto ~ period + (1 | herd) 
 Family: quasibinomial(logit link)
  AIC  BIC logLik deviance
 1145 1160 -569.6     1139
Random effects:
 Groups   Name        Variance Std.Dev.
 herd     (Intercept) 0.85697  0.92573 
 Residual             0.94584  0.97254 
number of obs: 1000, groups: herd, 50
 
Fixed effects:
            Estimate Std. Error t value
(Intercept) -0.05493    0.14949  -0.367
period       1.02879    0.08360  12.306
 
Correlation of Fixed Effects:
       (Intr)
period -0.001

 
############# here is the quasi-binomial result on R2.7.1 and Package lme4
version 0.999375-26: the variance estimates are very suspect
Generalized linear mixed model fit by the Laplace approximation 
Formula: toto ~ period + (1 | herd) 
  AIC  BIC logLik deviance
 1147 1167 -569.6     1139
Random effects:
 Groups   Name        Variance Std.Dev.
 herd     (Intercept) 0.12520  0.35384 
 Residual             0.13819  0.37173 
Number of obs: 1000, groups: herd, 50
 
Fixed effects:
            Estimate Std. Error t value
(Intercept) -0.05503    0.05714   -0.96
period       1.03117    0.03199   32.24
 
Correlation of Fixed Effects:
       (Intr)


 

	[[alternative HTML version deleted]]

Gosselin Frederic <frederic.gosselin <at> cemagref.fr> writes:
> 
> Dear R-helpers,
> 
> this mail is intended to mention a rather trange result and 
> generate potential useful comments on it. I am
> not aware of another posts on this issue 
> ( RSiteSearch("quasipoisson lmer version dispersion")).
> 
[snip snip]

 I would definitely forward this to 
r-sig-mixed-models at lists.r-project.org, if I were you.
At one point I reported what may be the same bug;

http://tinyurl.com/ye3pkve

 I don't know when it appeared and when it was
fixed, so I don't know whether it explains your
results or not.