[R-sig-ME] Reproducing results from an old lmer fit
Douglas Bates
bates at stat.wisc.edu
Thu Feb 26 19:17:49 CET 2009
On Thu, Feb 26, 2009 at 10:31 AM, Afshartous, David
<DAfshartous at med.miami.edu> wrote:
>
> All,
>
> For a paper revision I'm trying to reproduce some results from an old lmer
> fit with Rv2.7.1 prior to 5/28/08. However, when I currently load Rv2.7.1
> and lmer, the variance component estimates are slightly different than the
> original fit; the sessionInfo() is as follows:
>
>> sessionInfo()
> R version 2.7.1 (2008-06-23)
> i386-apple-darwin8.10.1
>
> locale:
> en_US.UTF-8/en_US.UTF-8/C/C/en_US.UTF-8/en_US.UTF-8
>
> attached base packages:
> [1] stats graphics grDevices utils datasets methods base
>
> other attached packages:
> [1] lme4_0.999375-24 Matrix_0.999375-11 lattice_0.17-8
>
> loaded via a namespace (and not attached):
> [1] grid_2.7.1 nlme_3.1-89
>
> Thus, I assume that I need to use the same older version of lme4 and/or
> Matrix which might be responsible for the difference in the results? If
> this is possible, how is this done?
>
> Cheers,
> David
>
> PS - for whatever it's worth, if I do the fit with lme (nlme_3.1-89) under
> Rv2.7.1 the results are closer to the original lmer results.
>
>
>
> ___________________________________________________
>
> Original lmer fit from 5/08:
> Model 2:
> AIC BIC logLik MLdeviance REMLdeviance
> 2813 2843 -1397 2829 2795
> Random effects:
> Groups Name Variance Std.Dev. Corr
> subject (Intercept) 2226.3 47.183
> Drug 2132.9 46.184 -0.865
> Residual 13673.6 116.934
>
> Current lmer fit:
> AIC BIC logLik deviance REMLdev
> 2815 2849 -1397 2830 2795
> Random effects:
> Groups Name Variance Std.Dev. Corr
> Patient_no (Intercept) 2165.1 46.531
> Drug.full.reverseC 1386.3 37.233 -1.000
> Residual 13947.5 118.100
Notice the large change in the estimated correlation with very little
change in the log-likelihood or deviance. This is an indication that
the model is over-specified.
Are you able to install R packages from the sources? If so, you could
try the branches/allcoef version from the SVN archive. On an
optimization problem like this it may be more successful in converging
to the global optimum instead of the local optimum.
This, by the way, is why I am always looking for better optimization
code to incorporate in R. The code in the nlme and lme4 packages just
evaluates the log-likelihood or the REML criterion for the model at
the observed data and a proposed value of the parameters. The actual
optimization is done by the nlminb optimizer which is based on very
old Fortran code written by David Gay. Even though the code is old
this optimizer is, in my experience, more reliable than the optimizers
used by optim and by nlm.
It is surprisingly difficult to find good optimization code that is
covered by an open source license. There is not a strong tradition of
open source code in the numerical analysis world. Many users are
enthused about the ipopt library (projects.coin-or.org/Ipopt) but even
though that code is open source it depends on other software, some of
which is commercial.
The optimization in lme4 is minimization of a real-valued function of
real parameters, some of which are subject to non-negativity
constraints. It is not an unconstrained optimization problem but the
constraints are very simple. The objective function can be evaluated
and, in theory, the gradient can also be evaluated. However, for
models with non-nested random effects evaluation of the gradient is
much, much more difficult and time consuming than is evaluation of the
objective function. Thus the ideal optimizer would allow for simple
"box constraints" on the parameters and would be derivative-free or at
least allow for numeric evaluation of the gradient. If anyone knows
of such code covered by a valid open-source license I would be
delighted to hear of it.
> Current lme fit:
> AIC BIC logLik
> 2814.611 2848.638 -1397.305
>
> StdDev Corr
> (Intercept) 47.21031 (Intr)
> Drug.full.reverseC 46.14014 -0.866
> Residual 116.93541
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
More information about the R-sig-mixed-models
mailing list