[R-sig-ME] Binomial GLMM with different versions of lme4

[Sundar Dorai-Raj]

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