While lme4 seems to be accurate and reliable over a wide range of cases, but there is a discrepancy between its results and those of other implementations (glmmML and glmmADMB from R and SAS PROC MIXED) for the (in)famous toenail data analyzed in a variety of places (see references below and Ch 10 of MLMUS). The toenail data are known to be somewhat pathological, and this discrepancy is not large (about 0.4% of the parameter values, and much smaller than the standard errors of the estimates). Nevertheless, it is puzzling and worrying; it would be nice to know what’s going on here.

library("foreign") ## for read.dta

Get data:

toenail <- read.dta("http://www.stata-press.com/data/r12/toenail.dta")
## summary(toenail)
toenail <- transform(toenail,
                     ## glmmML, glmmADMB prefer numeric 0/1 responses
t_glmer <- system.time(
    m_glmer <- glmer(outcome~treatment+visit+(1|patient),toenail,
t_glmmML <- system.time(
    m_glmmML <- glmmML(n.outcome~treatment+visit,
t_glmmADMB <- system.time(
    m_glmmADMB <- glmmadmb(n.outcome~treatment+visit+(1|patient),toenail,
##             RE_sd (Intercept)  treatment      visit    logLik elapsed
## glmer    4.687377   -1.051835 -0.6968895 -0.9115262 -626.1709   1.076
## glmmML   4.708440   -1.070209 -0.7005596 -0.9126392 -626.1627   0.101
## glmmADMB 4.723452   -1.070187 -0.7005568 -0.9126386 -626.1630 112.825
## SAS      4.708440   -1.070200 -0.7005000 -0.9126000 -626.1650   1.180

The glmmML and SAS PROC MIXED estimates are particularly close. However, for all methods other than glmer, and especially for the fixed effect estimates, the estimates are still much closer to each other than they are to glmer’s estimates.

Scaled RMSE of parameter estimates:

rmse <- function(x) sd(x)/mean(x)
(rmse_avg <- mean(abs(apply(ctab[,1:4],2,rmse))))
## [1] 0.003744686

These differences are not practically important: they are much smaller than the uncertainty in the parameters, e.g. for the fixed effects:

##          (Intercept) treatmentTerbinafine                visit 
##                0.800                0.687                0.074

Check to see if the issue is one of non-convergence/failure of the nonlinear optimizer to find the lowest-deviance parameters …

glmer_devfun <- update(m_glmer,devFunOnly=TRUE)
## 'log Lik.' -626.1709 (df=4)
par0 <- unlist(getME(m_glmer,c("theta","beta")))
par1 <- with(m_glmmML,c(sigma,coefficients))
## [1] 0.0008063085

lme4 thinks it’s correct, or at least that its solution is better than glmmML’s …

  1. De Backer, M., C. De Vroey, E. Lesaffre, I. Scheys, and P. De Keyser. 1998. “Twelve Weeks of Continuous Oral Therapy for Toenail Onychomycosis Caused by Dermatophytes: A Double-Blind Comparative Trial of Terbinafine 250 Mg/day versus Itraconazole 200 Mg/day.” Journal of the American Academy of Dermatology 38 (5, Supplement 2): S57–63. doi:10.1016/S0190-9622(98)70486-4.
  2. Lesaffre, Emmanuel, and Bart Spiessens. 2001. “On the Effect of the Number of Quadrature Points in a Logistic Random Effects Model: An Example.” Journal of the Royal Statistical Society: Series C (Applied Statistics) 50 (3): 325–35. doi:10.1111/1467-9876.00237.
  3. Rabe-Hesketh, Sophia, and Anders Skrondal. 2012. Multilevel and Longitudinal Modeling Using Stata, Volumes I and II, Third Edition. 3 edition. College Station, Tex: Stata Press.
  4. Vock, David M., Marie Davidian, and Anastasios A. Tsiatis. 2014. “SNP_NLMM: A SAS Macro to Implement a Flexible Random Effects Density for Generalized Linear and Nonlinear Mixed Models.” Journal of Statistical Software 56: 2. http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3969790/.