From: Roel de Jong <dejongroel_at_gmail.com>

Date: Tue 20 Dec 2005 - 09:57:14 EST

>>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).
*

>>sum(with(lesaffre, tapply(y, subject, mean)) == 0)

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

*>>
*

*>>
*

*>>
*

*>>
*

*>>
*

R-help@stat.math.ethz.ch mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Tue Dec 20 10:03:06 2005

Date: Tue 20 Dec 2005 - 09:57:14 EST

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

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@mi.uib.no> wrote: >

>>Douglas Bates wrote:

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

R-help@stat.math.ethz.ch mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Tue Dec 20 10:03:06 2005

*
This archive was generated by hypermail 2.1.8
: Fri 03 Mar 2006 - 03:41:40 EST
*