glmmADMB: Generalized Linear Mixed Models using AD Model Builder
Well, the dataset which I send in my previous message did without any
doubt come from a multilevel model (generated and fitted under the
binomial probit link), and gave the earlier posted error message while
fitting it with the latest version of lmer:
Warning: IRLS iterations for PQL did not converge
Error in objective(.par, ...) : Unable to invert singular factor of
downdated X'X
When data is generated from a specified model with reasonable parameter
values, it should be possible to fit such a model successful, or is this
me being stupid? We could try other parameter values, link functions
and/or more cases in a class if these given below are somehow implausible.
500 samples are drawn with the model specification (binomial probit):
y = (intercept*f1+pred2*f2+pred3*f3)+(intercept*ri+pred2*rs)
where pred2 and pred3 are predictors distributed N(0,1)
f1..f3 are fixed effects, f1=-1, f2=1.5, f3=0.5
ri is random intercept with associated variance var_ri=0.2
rs is random slope with associated variance var_rs=0.4
the covariance between ri and rs "covr"=0.1
1500 units/dataset, class size=30
Best regards,
Roel de Jong
Douglas Bates wrote:
On 12/19/05, Hans Julius Skaug <Hans.Skaug at mi.uib.no> wrote:
Douglas Bates wrote:
The "Laplace" method in lmer and the default method in glmm.admb, which according to the documentation is the Laplace approximation, produce essentially the same model fit. One difference is the reported value of the log-likelihood, which we should cross-check, and another difference is in the execution time
Yes, glmmADMB has sqrt(2*pi) constants missing. Thanks for pointing that out. Execution time: As pointed out by Roel de Jong, the underlying software AD Model Builder does not use hand-coded derivatives for the Hessian involved in the Laplace approximation, but calculates these by automatic differentiation. There is a cost in terms of execution speed, but on the other hand it is very quick to develop new models, as you do not have to worry about derivatives. I hope to exploit this beyond standard GLMMs, and provide other R packages. Comparison of glmmADMB with lmer: I find that the two packages do not give the same result on one of the standard datasets in the literature (Lesaffre et. al., Appl. Statist. (2001) 50, Part3, pp 325-335).
Ah yes, that example. It is also given as the 'toenail' data set in the 'mlmus' package of data sets from the book "Multilevel and Longitudinal Modeling Using Stata" by Sophia Rabe-Hesketh and Anders Skrondal (Stata Press, 2005). It is not surprising that it is difficult to fit such a model to these data because the data do not look like they come from such a model. You did not include the estimates of the variance of the random effects in your output. It is very large and very poorly determined. If you check the distribution of the posterior modes of the random effects (for linear mixed models these are called the BLUPs - Best Linear Unbiased Predictors - and you could call them BLUPs here too except for the fact that they are not linear and they are not unbiased and there isn't a clear sense in which they are "best") it is clearly not a Gaussian distribution with mean zero. I enclose a density plot. You can see that it is bimodal and the larger of the two peaks is for a negative value. These are the random effects for those subjects that had no positive responses - 163 out of the 294 subjects.
sum(with(lesaffre, tapply(y, subject, mean)) == 0)
[1] 163 There is no information to estimate the random effects for these subjects other than "make it as large and negative as possible". It is pointless to estimate the fixed effects for such a clearly inappropriate model. lattice package.
The full set of R commands used to download data and fit the model is given at the end of this email.
fit_lmer_lapl <- lmer(y~ treat + time + (1|subject),data=lesaffre,family=binomial,method="Laplace")
Warning message:
optim or nlminb returned message ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
in: LMEopt(x = mer, value = cv)
PART OF OUTPUT:
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.626214 0.264996 -2.3631 0.01812 *
treat -0.304660 0.360866 -0.8442 0.39853
time -0.346605 0.026666 -12.9979 < 2e-16 ***
The corresponding estimates with glmmADMB is:
fit_glmmADMB <- glmm.admb(y~ treat + time,random=~1,group="subject",data=lesaffre,family="binomial",link="logit")
PART OF OUTPUT:
Fixed effects:
Log-likelihood: -359.649
Formula: y ~ treat + time
(Intercept) treat time
-2.33210 -0.68795 -0.46134
So, the estimates of the fixed effects differ. lmer() does infact produces a warning, and it appears that
it method="Laplace" and method="PQL" produce the same results.
Best regards,
hans
# Load data
source("http://www.mi.uib.no/~skaug/cash/lesaffre_dat.s")
# Run lmer
library(lme4)
fit_lmer <- lmer(y~ treat + time + (1|subject),data=lesaffre,family=binomial)
fit_lmer_lapl <- lmer(y~ treat + time + (1|subject),data=lesaffre,family=binomial,method="Laplace")
# Run glmmADMB
library(glmmADMB)
example(glmm.admb) # Must be run once in each new directory (this feature will be removed in future version of glmmADMB).
fit_glmmADMB <- glmm.admb(y~ treat + time,random=~1,group="subject",data=lesaffre,family="binomial",link="logit")