[R] glmmADMB: Generalized Linear Mixed Models using AD Model Builder

[Roel de Jong]

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")
>>
>>
>>
>>
>>