From: dave fournier <otter_at_otter-rsch.com>

Date: Fri 14 Oct 2005 - 06:19:24 EST

Date: Fri 14 Oct 2005 - 06:19:24 EST

Do Users of Nonlinear Mixed Effects Models Know

Whether Their Software Really Works?

Lesaffre et. al. (Appl. Statist. (2001) 50, Part3, pp 325-335)
analyzed

some simple clinical trials data using a logistic random effects
model. Several packages and methods MIXOR, SAS NLMIXED were employed.
They reported obtaining very different parameter estimates and
P values for the log-likelihood with the different packages and
methods. We thought it would be interesting to revisit this example
using the AD Model Builder random effects module which we feel is
the most stable software available for this problem at this time.

http://otter-rsch.com/admodel.htm

You can get Table 2 from Lesaffre et al at

http://otter-rsch.com/downloads/other/lesaffre.pdf

The data and more information are available from the publisher.

http://www.blackwellpublishers.co.uk/rss/Volumes/Cv50p3.htm

We considered three questions:

1.) What are the estimates using the Laplace approximation for integrating out the random effects. 2.) What are the exact MLE's. 3.) How well does hypothesis testing (likelihood-ratio) using the Laplace approximation compare with the exact MLE's.

We first fit the data using ADMB-RE's Laplace approximation option.

Laplace approximation estimates:

# Number of parameters = 4 log-likelihood = -629.817 value std dev P value b_1 -2.3321e+00 7.6973e-01 < 0.0024 b_2 -6.8795e-01 6.6185e-01 0.298 b_3 -4.6134e-01 4.0000e-02 < 0.001 sigma 4.5738e+00 7.0970e-01

The parameter of interest here the treatment effect b_2 which is the parameter reported in Lesaffre et. al.

To calculate the exact MLE we fit the model using 100 point adaptive Gaussian integration. The ADMB-RE results were:

Gaussian integration estimates:

# Number of parameters = 4 log-likelihood = -627.481

name value std dev P value b_1 -1.4463e+00 4.2465e-01 < 0.001 b_2 -5.2225e-01 5.5716e-01 0.348 b_3 -4.5150e-01 3.6663e-02 < 0.001 sigma 4.0137e+00 3.8083e-01

Of the estimates reported in Lesaffre et al. in table 2 only the 50 point quadrature for the program MIXOR appear to be correct for both the log-likelihood value and the parameter estimates while the authors concluded that the SAS NLMIXED parameter estimates they obtained were correct. So even though these authors were looking for pathological behaviour and were presumably very careful, and their paper was presumably peer-reviewed, they came to the wrong conclusion using SAS NLMIXED.

How do we know that our exact MLE's are correct? To confirm our results we used our parameter estimates as initial values in the SAS NLMIXED procedure using 100 point adaptive quadrature. The procedure returned our values, that is it agreed that these are the maximum likelihood estimates. However we verified that to get these estimates from the SAS NLMIXED procedure one must begin with fairly good starting values. In contrast the ADMB-RE procedure is very insensitive to the starting values used. Our conclusion is that while SAS NLMIXED might work for this very simple problem it probably begins to break down when the problem is a bit more difficult.

The ADMB-RE software is more stable because it calculates exact higher oreder derivatives by automatic differentiation for use in its optimization procedure and calculations while other packages do not.

Gauss-Hermite integration for the random effects can be used for this model because the Hessian for the random effects is diagonal which permits one dimensional integration over the random effects to great accuracy. However this procedure does not scale well to problems where the Hessian is not diagonal. Suppose that it takes a 20 point quadrature to obtain reliable parameter estimates with a diagonal Hessian. Then with a block diagonal Hessian where the blocks are of size 4x4 it would take 160,000 points.

Results using R

We fit the model using what appear to be the currently available procedures in R. The two routines lmer (lme4 package) and glmmPQL (MASS library) were tried.

The call

>> lmer(y ~ treat + time + (1|subject),data=lesaffre,family=binomial)

resulted in a warning message from lme4() but both routines produced the same results.

Generalized linear mixed model fit using PQL Formula: y ~ treat + time + (1 | subject) Data: lesaffre Family: binomial(logit link) AIC BIC logLik deviance 1305.859 1333.628 -647.9295 1295.859 Random effects: Groups Name Variance Std.Dev. subject (Intercept) 6.8059 2.6088 # of obs: 1908, groups: subject, 294 Estimated scale (compare to 1) 0.9091945 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 *** sigma 2.608 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Warning message: optim or nlminb returned message ERROR: ABNORMAL_TERMINATION_IN_LNSRCH in: LMEopt(x = mer, value = cv)

The R routine correctly identifies the treatment effect as not significant. However the parameter estimates are poor.

Likelihood Ratio Testing

Accurate calculation of the log-likelihood value is desirable so that hypothesis testing can be carried out using likelihood ratio tests.

However as noted above the use of Gaussian integration is not practical for many nonlinear mixed models. We were interested in seeing how well the use of the approximate log-likelihood values produced by ADMB-RE's Laplace approximation option would perform.

We consider the alternative model with an extra interaction term (b_4) from Lesaffre et al.

Here are the results for the laplace approximation:

# Number of parameters = 5 log-likelihood = -627.809 name value std dev P vlaue b_1 -2.5233e+00 7.8829e-01 < 0.002 b_2 -3.0702e-01 6.8996e-01 0.655 b_3 -4.0009e-01 4.7059e-02 < 0.001 b_4 -1.3726e-01 6.9586e-02 0.044 sigma 4.5783e+00 7.2100e-01

and the exact parameter estimates by 100 point Gaussian integration.

# Number of parameters = 5 log-likelihood = -625.398 name value std dev P value b_1 -1.6183e+00 4.3427e-01 < 0.001 b_2 -1.6077e-01 5.8394e-01 0.783 b_3 -3.9100e-01 4.4380e-02 < 0.001 b_4 -1.3679e-01 6.8013e-02 0.044 sigma 4.0131e+00 3.8044e-01

The log-likelihood differences are 2.01 for the Laplace approximation and 2.08 for Gaussian integration. Since the 95% point for hypothesis testing is 1.92 use of either model results in acceptance of the interaction term.

Conclusions

With the exception of AD Model Builder random effect module none of the packages tested appear to function reliably for this problem. SAS NLMIXED was beginning to exhibit symptoms of instability which would probably render it unreliable on more difficult problems. We can see no reason for using "quasi-likelihoods" to fit nonlinear mixed models when ADMB-RE can fit the models by maximum likelihood with all the advantages that ensue.

Note

We realize that there are many other packages out there. We would welcome results for other packages. If we can find a serious competitor to AD Model Builder then we could move on to comparing the relative performance on more difficult models.

Cheers,

Dave Fournier

-- David A. Fournier P.O. Box 2040, Sidney, B.C. V8l 3s3 Canada http://otter-rsch.com -- Internal Virus Database is out-of-date. ______________________________________________ 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.htmlReceived on Fri Oct 14 06:20:26 2005

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