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.