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
library("lme4")
library("glmmML")
library("glmmADMB")
```

Get data:

```
toenail <- read.dta("http://www.stata-press.com/data/r12/toenail.dta")
## summary(toenail)
toenail <- transform(toenail,
patient=factor(patient),
## glmmML, glmmADMB prefer numeric 0/1 responses
n.outcome=as.numeric(outcome)-1)
```

```
t_glmer <- system.time(
m_glmer <- glmer(outcome~treatment+visit+(1|patient),toenail,
family=binomial))
t_glmmML <- system.time(
m_glmmML <- glmmML(n.outcome~treatment+visit,
cluster=toenail$patient,toenail,
family="binomial"))
t_glmmADMB <- system.time(
m_glmmADMB <- glmmadmb(n.outcome~treatment+visit+(1|patient),toenail,
family="binomial"))
```

```
## 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)
logLik(m_glmer)
```

`## 'log Lik.' -626.1709 (df=4)`

```
par0 <- unlist(getME(m_glmer,c("theta","beta")))
par1 <- with(m_glmmML,c(sigma,coefficients))
glmer_devfun(par1)-glmer_devfun(par0)
```

`## [1] 0.0008063085`

`lme4`

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

’s …

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