[R] problem with lme in nlme package
John Fox
jfox at mcmaster.ca
Thu May 2 04:22:28 CEST 2002
Dear R list members,
I've turned up a strange discrepancy between results obtained from the lme
function in the nlme package in R and results obtained with lme in S-PLUS.
I'm using version 3.1-24 of nlme in R 1.4.1 under Windows 2000, and both
S-PLUS 2000 and 6.0, again under Windows 2000.
I've noticed discrepancies in a couple of instances. Here's one, using data
from Bryk and Raudenbush's Hierarchical Linear models:
From R:
> attach(Bryk)
> cses <- meanses <- numeric(length(ses)) # initialize
> for (sc in unique(school)) {
+ meanses[school==sc] <- mean(ses[school==sc])
+ cses[school==sc]<-ses[school==sc]- meanses[school==sc]
+ }
> Bryk$cses <- cses
> Bryk$meanses <- meanses
> Bryk[1:10,]
school ses mathach sector cses meanses
1 1224 0.022 4.583 public 0.456383 -0.43438
2 1224 0.332 20.349 public 0.766383 -0.43438
3 1224 0.372 6.714 public 0.806383 -0.43438
4 1224 -0.078 16.405 public 0.356383 -0.43438
5 1224 -0.158 17.898 public 0.276383 -0.43438
6 1224 -0.298 19.338 public 0.136383 -0.43438
7 1224 -0.458 21.521 public -0.023617 -0.43438
8 1224 -0.468 3.154 public -0.033617 -0.43438
9 1224 -0.468 21.178 public -0.033617 -0.43438
10 1224 -0.528 20.349 public -0.093617 -0.43438
> remove(cses, meanses)
> Bryk$sector <- factor(Bryk$sector, levels=c('public', 'Catholic'))
> contrasts(Bryk$sector)
Catholic
public 0
Catholic 1
>
> bryk.lme.1<-lme(mathach ~ meanses*cses + sector*cses,
+ random = ~ cses | school,
+ data=Bryk)
> summary(bryk.lme.1)
Linear mixed-effects model fit by REML
Data: Bryk
AIC BIC logLik
46525 46594 -23252
Random effects:
Formula: ~cses | school
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 1.541177 (Intr)
cses 0.018174 0.006
Residual 6.063492
Fixed effects: mathach ~ meanses * cses + sector * cses
Value Std.Error DF t-value p-value
(Intercept) 12.1282 0.19920 7022 60.886 <.0001
meanses 5.3367 0.36898 157 14.463 <.0001
cses 2.9421 0.15122 7022 19.456 <.0001
sectorCatholic 1.2245 0.30611 157 4.000 1e-04
meanses:cses 1.0444 0.29107 7022 3.588 3e-04
cses:sectorCatholic -1.6421 0.23312 7022 -7.044 <.0001
Correlation:
(Intr) meanss cses sctrCt mnss:c
meanses 0.256
cses 0.000 0.000
sectorCatholic -0.699 -0.356 0.000
meanses:cses 0.000 0.000 0.295 0.000
cses:sectorCatholic 0.000 0.000 -0.696 0.000 -0.351
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-3.170106 -0.724877 0.014892 0.754263 2.965498
Number of Observations: 7185
Number of Groups: 160
>
From S-PLUS 6.0 (identical results from S-PLUS 2000):
> attach(Bryk)
> cses <- meanses <- numeric(length(ses)) # initialize
> for (sc in unique(school)) {
+ meanses[school==sc] <- mean(ses[school==sc])
+ cses[school==sc]<-ses[school==sc]- meanses[school==sc]
+ }
> Bryk$cses <- cses
> Bryk$meanses <- meanses
> Bryk[1:10,]
school ses mathach sector cses meanses
1 1224 0.022 4.583 public 0.456383 -0.43438
2 1224 0.332 20.349 public 0.766383 -0.43438
3 1224 0.372 6.714 public 0.806383 -0.43438
4 1224 -0.078 16.405 public 0.356383 -0.43438
5 1224 -0.158 17.898 public 0.276383 -0.43438
6 1224 -0.298 19.338 public 0.136383 -0.43438
7 1224 -0.458 21.521 public -0.023617 -0.43438
8 1224 -0.468 3.154 public -0.033617 -0.43438
9 1224 -0.468 21.178 public -0.033617 -0.43438
10 1224 -0.528 20.349 public -0.093617 -0.43438
> remove(c('cses', 'meanses'))
> options(contrasts=c('contr.treatment', 'contr.poly'))
>
> Bryk$sector <- factor(Bryk$sector, levels=c('public', 'Catholic'))
> contrasts(Bryk$sector)
Catholic
public 0
Catholic 1
>
> bryk.lme.1<-lme(mathach ~ meanses*cses + sector*cses,
+ random = ~ cses | school,
+ data=Bryk)
> summary(bryk.lme.1)
Linear mixed-effects model fit by REML
Data: Bryk
AIC BIC logLik
46524 46592 -23252
Random effects:
Formula: ~ cses | school
Structure: General positive-definite
StdDev Corr
(Intercept) 1.54278 (Inter
cses 0.32004 0.395
Residual 6.05976
Fixed effects: mathach ~ meanses * cses + sector * cses
Value Std.Error DF t-value p-value
(Intercept) 12.128 0.19931 7022 60.851 <.0001
meanses 5.333 0.36920 157 14.444 <.0001
cses 2.945 0.15565 7022 18.921 <.0001
sector 1.227 0.30630 157 4.005 0.0001
meanses:cses 1.039 0.29899 7022 3.476 0.0005
sector:cses -1.643 0.23986 7022 -6.849 <.0001
Correlation:
(Intr) meanss cses sector mnss:c
meanses 0.256
cses 0.076 0.020
sector -0.699 -0.356 -0.053
meanses:cses 0.019 0.076 0.293 -0.027
sector:cses -0.053 -0.028 -0.696 0.078 -0.351
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-3.1591 -0.72314 0.01723 0.75461 2.9581
Number of Observations: 7185
Number of Groups: 160
Notice that the results are very close, except for the standard deviation
of the random cses effect and the correlation between the random effects
for cses and the intercept. Published results and SAS PROC MIXED agree with
the S-PLUS output.
I don't think that I've done anything wrong in setting up the problem in R.
In any event, any suggestions would be greatly appreciated. I could make
the data available, if that would help.
Thanks,
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