[R] problem with lme in nlme package

[John Fox]

Dear Doug,

First, thanks for getting back to me.

At 08:18 AM 5/2/2002 -0500, Douglas Bates wrote:
>John,
>
>I'm travelling today and won't be able to check this discrepancy.  I
>will try to look at it tomorrow.

There's no rush.

>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


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?)

>  - I believe the data are already in the nlme package as the
>    MathAchieve and MathAchSchool data sets.

Yes, they are -- I didn't know that. Thanks.

>  - In Bryk's analysis the meanses is calculated incorrectly.

I was aware of this -- it's why I recalculated the means rather than using 
Bryk and Raudenbush's -- but the results don't change much.

>  - Some analysis of these data is given in the slides at
>    http://www.stat.wisc.edu/~st850-1/MEcourse6_4up.pdf

Again, thanks. I'll take a look. I'm working up the example for use in my 
appendix on mixed models, a draft of which is at 
<http://socserv.socsci.mcmaster.ca/jfox/Books/Companion/appendix-mixed-models.pdf>.

>  - It is much easier (and faster) to calculate the meanses gsummary.
>    I think it would be gsummary(MathAchieve, which = "ses")

The versions of gsummary that I have don't seem to take a which argument, 
but tapply works (and is faster than my previous code, though in a data set 
this size both are essentially instantaneous):

 >    mses <- tapply(ses, school, mean)
 >    meanses <- mses[as.character(school)]

I expect that the simple answer to my basic question is that the problem is 
ill-conditioned. Interesting, since Bryk and Raudenbush's data have been 
used a lot.

Thanks for the suggestions and for your help.

John





-----------------------------------------------------
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario, Canada L8S 4M4
email: jfox at mcmaster.ca
phone: 905-525-9140x23604
web: www.socsci.mcmaster.ca/jfox
-----------------------------------------------------

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