[R] problem with lme in nlme package

[Douglas Bates]

John Fox <jfox at mcmaster.ca> writes:

> At 08:18 AM 5/2/2002 -0500, Douglas Bates wrote:
> 
> >A couple of things to consider:
> >
> >  - the underlying optimization software is different in S-PLUS and R.
> >    It may be that the REML estimates are poorly defined and
> >    convergence is being declared in different places by different software.
> >    Have you checked the intervals on those variance components using
> >    intervals(bryk.lme.1)
> 
> It occurred to me that the problem may be ill-conditioned, but I
> didn't check the confidence intervals for the variance components. I
> did look at the log-likelihoods, though. To more decimal places than
> in my previous note, they are:
> 
> 
> R: -23252.3919
> 
> S-PLUS: -23253.2177

Apparently the optimizer in S-PLUS declared convergence prematurely.

> Getting the confidence intervals on the variance components produces:
> 
> R:
> 
>      > intervals(bryk.lme.1)
>      Approximate 95% confidence intervals
> 
>      Fixed effects:
>                              lower      est.     upper
>      (Intercept)         11.7377216 12.128207 12.518692
>      meanses              4.6078631  5.336665  6.065467
>      cses                 2.6457000  2.942145  3.238589
>      sectorCatholic       0.6198986  1.224531  1.829164
>      meanses:cses         0.4738118  1.044406  1.615000
>      cses:sectorCatholic -2.0991261 -1.642148 -1.185170
> 
>      Random Effects:
>      Level: school
>                                  lower        est.        upper
>      sd((Intercept))        1.321746e+00 1.541176850 1.797037e+00
>      sd(cses)               6.686558e-09 0.018173637 4.939478e+04
>      cor((Intercept),cses) -8.630434e-02 0.005809486 9.782482e-02
> 
>      Within-group standard error:
>      lower     est.    upper
>      5.964008 6.063492 6.164636
> 
> S-PLUS
> 
>      > intervals(bryk.lme.1)
>      Problem in intervals.lme(bryk.lme.1): Cannot get confidence inter
>      vals on var-cov components: Non-positive definite approximate var
>      iance-covariance
>      Use traceback() to see the call stack

> Obviously, there's a problem here, as you suspected: the interval for
> sd(cses) in the R solution is very broad, and the covariance matrix in
> the S-PLUS solution isn't positive definite. (Does the S-PLUS version
> of the software use the log-Cholesky parametrization of the covariance
> matrix for the random effects?)

The default parameterization in the S-PLUS version of lme is the
Matrix logarithm.  The difference in parameterization should not have
a major effect on the convergence as both of these are unconstrained,
nonredundant parameterizations.  I think the difference in the results
from different optimization methods is an indication that the optimum
is poorly defined, as you can see from the intervals.

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._