[R] problem with lme in nlme package
John Fox
jfox at mcmaster.ca
Thu May 2 16:50:47 CEST 2002
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
More information about the R-help
mailing list