[R] glmmADMB: Generalized Linear Mixed Models using AD Model Builder
Roel de Jong
dejongroel at gmail.com
Mon Dec 19 23:57:14 CET 2005
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
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
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)
>  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")
>>optim or nlminb returned message ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
>> in: LMEopt(x = mer, value = cv)
>>PART OF OUTPUT:
>> 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:
>> 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.
>># Load data
>># Run lmer
>>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
>>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")
More information about the R-help