[R-sig-ME] Reproducing results from an old lmer fit

[Douglas Bates]

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
>