[R-sig-ME] Binomial GLMM with different versions of lme4
Sundar Dorai-Raj
sundar.dorai-raj at pdf.com
Tue Oct 21 23:57:50 CEST 2008
Hi, Prof. Bates (and all other interested parties),
I haven't been using the lme4 package recently for binomial modeling but
have recently returned for some analysis and found some puzzling
problems. The current package version (0.999375-27) incorrectly
estimates the dispersion parameter when using a quasibinomial. I think
this is referred to indirectly in the following thread:
https://stat.ethz.ch/pipermail/r-sig-mixed-models/2008q4/001501.html
However, I have attached a working example that shows (I believe) the
last known package where the dispersion is believable. For now, I'm
using the older package to do the analysis but would be grateful for any
corrections or workarounds.
<code>
# library(lme4, lib = "D:/R/archive") ## version 0.99875-9
library(lme4) ## version 0.999375-27
set.seed(2)
## simulate data
A <- factor(rep(1:10, each = 10))
n <- length(A)
rA <- rnorm(nlevels(A))
y <- rbinom(n, 1, binomial()$linkinv(0.5 + rA[A]))
(f <- lmer(y ~ (1 | A), family = quasibinomial))
extractVar.lmer <- function (object, ...) {
vc <- VarCorr(object)
vdc <- unlist(lapply(lapply(vc, as.matrix), diag))
names(vdc) <- sub("\\.\\(.*\\)", "", names(vdc))
Var <- c(vdc, attr(vc, "sc")^2)
names(Var)[length(Var)] <- "residual"
Var
}
v <- extractVar.lmer(f)
c(v[1]/v[2], v[2])
sessionInfo()
</code>
<output version="0.99875-9">
> (f <- lmer(y ~ (1 | A), family = quasibinomial))
Generalized linear mixed model fit using Laplace
Formula: y ~ (1 | A)
Family: quasibinomial(logit link)
AIC BIC logLik deviance
130 135.2 -63 126
Random effects:
Groups Name Variance Std.Dev.
A (Intercept) 0.22075 0.46984
Residual 0.95095 0.97517
number of obs: 100, groups: A, 10
Fixed effects:
Estimate Std. Error t value
(Intercept) 0.7443 0.2572 2.893
>
> v <- extractVar.lmer(f)
> c(v[1]/v[2], v[1])
A residual
0.2321391 0.9509488
> sessionInfo()
R version 2.7.2 (2008-08-25)
i386-pc-mingw32
locale:
LC_COLLATE=English_United States.1252;LC_CTYPE=English_United
States.1252;LC_MONETARY=English_United
States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] lme4_0.99875-9 Matrix_0.999375-15 lattice_0.17-15
loaded via a namespace (and not attached):
[1] grid_2.7.2
</output>
<output version="0.999375-27">
> (f <- lmer(y ~ (1 | A), family = quasibinomial))
Generalized linear mixed model fit by the Laplace approximation
Formula: y ~ (1 | A)
AIC BIC logLik deviance
132 139.8 -63 126
Random effects:
Groups Name Variance Std.Dev.
A (Intercept) 0.046833 0.21641
Residual 0.205057 0.45283
Number of obs: 100, groups: A, 10
Fixed effects:
Estimate Std. Error t value
(Intercept) 0.7463 0.1191 6.264
>
> v <- extractVar.lmer(f)
> c(v[1]/v[2], v[1])
A residual
0.2283881 0.2050567
>
> sessionInfo()
R version 2.7.2 (2008-08-25)
i386-pc-mingw32
locale:
LC_COLLATE=English_United States.1252;LC_CTYPE=English_United
States.1252;LC_MONETARY=English_United
States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] lme4_0.999375-27 Matrix_0.999375-15 lattice_0.17-15
loaded via a namespace (and not attached):
[1] grid_2.7.2
</output>
Thanks,
--sundar
More information about the R-sig-mixed-models
mailing list