[R] Using quasibinomial family in lmer

[Steven McKinney]

Hi Russell,


You might join in on the discussion on the R-sig-ME
mixed effects listserv

(see e.g. Doug Bates' message today (Sept 16 2008) titled 
Re: [R-sig-ME] glmer and overdispersed Poisson models)

There Prof. Bates suggests that older versions of lmer() may be doing
things more appropriately (though maybe not yet optimally).
(He also asks for input on how this might be
addressed - I'm not yet quasi-specialized enough to
figure it out but maybe someone reading this will
be.)

Here's what I get with a slightly older version of lmer - note that
my quasi-binomial t-value is not as extreme as yours.

Note also that for the output below, the ratio of the 
standard deviation estimates for the intercept estimate
> 0.1219/0.07441
[1] 1.638221
equals the residual standard deviation, so it appears
that for the quasibinomial case lmer or some lower
level routine estimates the scale parameter phi using the
residual standard deviation, though I have not delved
into the source code to verify this.  So some attempt
to modify scale appears to be in place for this version
of lmer.


library(lme4)
set.seed(12)
eta=rnorm(50)
p=exp(eta)/(1+exp(eta))
y=rbinom(50,20,p)/20 #IID overdispersed binomial-normal proportions
#y=rbinom(50,20,0.5)/20 #IID binomial(20,0.5)
Group=rep(c("A","B","C","D","E"),c(10,10,10,10,10))
lmer(y~1+(1|Group),weights=rep(20,50),family="binomial")
lmer(y~1+(1|Group),weights=rep(20,50),family="quasibinomial")
sessionInfo()

> library(lme4)
> set.seed(12)
> eta=rnorm(50)
> p=exp(eta)/(1+exp(eta))
> y=rbinom(50,20,p)/20 #IID overdispersed binomial-normal proportions
> #y=rbinom(50,20,0.5)/20 #IID binomial(20,0.5)
> Group=rep(c("A","B","C","D","E"),c(10,10,10,10,10))
> lmer(y~1+(1|Group),weights=rep(20,50),family="binomial")
Generalized linear mixed model fit using Laplace 
Formula: y ~ 1 + (1 | Group) 
 Family: binomial(logit link)
   AIC   BIC logLik deviance
 149.7 153.6 -72.87    145.7
Random effects:
 Groups Name        Variance  Std.Dev.
 Group  (Intercept) 0.0075745 0.087032
number of obs: 50, groups: Group, 5

Estimated scale (compare to  1 )  1.638542 

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept) -0.13739    0.07441  -1.847   0.0648 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
> lmer(y~1+(1|Group),weights=rep(20,50),family="quasibinomial")
Generalized linear mixed model fit using Laplace 
Formula: y ~ 1 + (1 | Group) 
 Family: quasibinomial(logit link)
   AIC   BIC logLik deviance
 149.7 153.6 -72.87    145.7
Random effects:
 Groups   Name        Variance Std.Dev.
 Group    (Intercept) 0.020336 0.14261 
 Residual             2.684818 1.63854 
number of obs: 50, groups: Group, 5

Fixed effects:
            Estimate Std. Error t value
(Intercept)  -0.1374     0.1219  -1.127
> sessionInfo()
R version 2.7.1 Patched (2008-06-25 r45988) 
powerpc-apple-darwin8.11.1 

locale:
C

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] RMySQL_0.6-0      DBI_0.2-4         RODBC_1.2-3       lme4_0.99875-9   
[5] Matrix_0.999375-9 lattice_0.17-10  

loaded via a namespace (and not attached):
[1] grid_2.7.1
> 



Steven McKinney

Statistician
Molecular Oncology and Breast Cancer Program
British Columbia Cancer Research Centre

email: smckinney +at+ bccrc +dot+ ca

tel: 604-675-8000 x7561

BCCRC
Molecular Oncology
675 West 10th Ave, Floor 4
Vancouver B.C. 
V5Z 1L3
Canada




-----Original Message-----
From: r-help-bounces at r-project.org on behalf of r.millar at auckland.ac.nz
Sent: Tue 9/16/2008 3:18 PM
To: r-help at r-project.org
Subject: [R] Using quasibinomial family in lmer
 
Dear R-Users,

I can't understand the behaviour of quasibinomial in lmer. It doesn't
appear to be calculating a scaling parameter, and looks to be reducing the
standard errors of fixed effects estimates when overdispersion is present
(and when it is not present also)! A simple demo of what I'm seeing is
given below. Comments appreciated?

Thanks,

Russell Millar
Dept of Stat
U. Auckland

PS. I'm using the latest version of lme4 (0.999375-26) with R 2.7.2.

> eta=rnorm(50)
> p=exp(eta)/(1+exp(eta))
> y=rbinom(50,20,p)/20 #IID overdispersed binomial-normal proportions
> #y=rbinom(50,20,0.5)/20 #IID binomial(20,0.5)
>
> Group=rep(c("A","B","C","D","E"),c(10,10,10,10,10))
>
> #library(lme4)
>
> lmer(y~1+(1|Group),weights=rep(20,50),family="binomial")
Generalized linear mixed model fit by the Laplace approximation
Formula: y ~ 1 + (1 | Group)
 AIC   BIC logLik deviance
 211 214.8 -103.5      207
Random effects:
 Groups Name        Variance Std.Dev.
 Group  (Intercept) 0.072891 0.26998
Number of obs: 50, groups: Group, 5

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)   0.2194     0.1367   1.605    0.108
>
> lmer(y~1+(1|Group),weights=rep(20,50),family="quasibinomial")
Generalized linear mixed model fit by the Laplace approximation
Formula: y ~ 1 + (1 | Group)
 AIC   BIC logLik deviance
 213 218.7 -103.5      207
Random effects:
 Groups   Name        Variance  Std.Dev.
 Group    (Intercept) 0.0032632 0.057125
 Residual             0.0447685 0.211586
Number of obs: 50, groups: Group, 5

Fixed effects:
            Estimate Std. Error t value
(Intercept)  0.21940    0.02892   7.586
>

______________________________________________
R-help at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.