# [R] Using quasibinomial family in lmer

Steven McKinney smckinney at bccrc.ca
Wed Sep 17 03:00:01 CEST 2008

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

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

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