|
|
|
||||
|
|
|||||
This paper can be downloaded from Professor Singer's web site at
http://gseweb.harvard.edu/%7Efaculty/singer/Papers/sasprocmixed.pdf
Most of the examples use the nlme package which comes with R and can be made available by typing
> library(nlme)
First set of examples using the hsb12 data file.
To input the hsb12.csv data file in R, use the following code. Change d:/data/ to the where you stored the data files on your computer (and do use forward slashes to separate directory names)
> hsb12 <- read.csv("d:/data/hsb12.csv", header=T)
Page 329
In order to use the lme function the data object must already have the
multilevel structure which is imposed by the groupedData function.
> hsb12.int <- groupedData( mathach ~ 1 | school, data=hsb12)
> pg329 <- lme(mathach ~ 1, data=hsb12.int, random = ~ 1 | school)
> summary(pg329)
Linear mixed-effects model fit by REML
Data: hsb12.int
AIC BIC logLik
47122.79 47143.43 -23558.4
Random effects:
Formula: ~ 1 | school
(Intercept) Residual
StdDev: 2.934966 6.256862
Fixed effects: mathach ~ 1
Value Std.Error DF t-value p-value
(Intercept) 12.63697 0.2443936 7025 51.70747 <.0001
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-3.063125 -0.753874 0.02670132 0.7606217 2.742626
Number of Observations: 7185
Number of Groups: 160
Page 331 Table 2, random intercept, predictor = meanses
> hsb12.meanses <- groupedData(mathach ~ 1 | school, data=hsb12, outer= ~ meanses)
> pg331 <- lme(mathach ~ meanses, data=hsb12.meanses, random = ~ 1 | school)
> summary(pg331)
Linear mixed-effects model fit by REML
Data: hsb12.meanses
AIC BIC logLik
46969.28 46996.8 -23480.64
Random effects:
Formula: ~ 1 | school
(Intercept) Residual
StdDev: 1.62441 6.257562
Fixed effects: mathach ~ meanses
Value Std.Error DF t-value p-value
(Intercept) 12.64944 0.1492801 7025 84.73622 <.0001
meanses 5.86354 0.3614580 158 16.22191 <.0001
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-3.134798 -0.7525634 0.02408657 0.7677326 2.785014
Number of Observations: 7185
Number of Groups: 160
Page 335
# compute cses
> hsb12$cses <- (hsb12$ses - hsb12$meanses)
> hsb12.cses <- groupedData( mathach ~ cses | school, data=hsb12)
# use unstructured correlation (default in lme function, corr = NULL)
> pg335 <- lme(mathach ~ cses, data=hsb12.cses, random = ~ 1 + cses | school )
> summary(pg335)
Linear mixed-effects model fit by REML
Data: hsb12.cses
AIC BIC logLik
46726.24 46767.52 -23357.12
Random effects:
Formula: ~ cses | school
Structure: General positive-definite
StdDev Corr
(Intercept) 2.9462163 (Inter
cses 0.8328983 0.011
Residual 6.0580812
Fixed effects: mathach ~ cses
Value Std.Error DF t-value p-value
(Intercept) 12.64925 0.2444948 7024 51.73627 <.0001
cses 2.19291 0.1282521 7024 17.09842 <.0001
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-3.096421 -0.7313269 0.01767674 0.7537266 2.90052
Number of Observations: 7185
Number of Groups: 160
Pages 337-338,
> hsb12.big <- groupedData( mathach ~ cses | school, data=hsb12)
> pg337 <- lme(mathach ~ cses + meanses + sector + cses:meanses + cses:sector,
data=hsb12.big, random = ~ 1 + cses | school, method="REML", corr=NULL)
> summary(pg337)
Linear mixed-effects model fit by REML
Data: hsb12.big
AIC BIC logLik
46524.79 46593.58 -23252.40
Random effects:
Formula: ~1 + cses | school
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 1.54128855 (Intr)
cses 0.01829966 0.006
Residual 6.06348796
Fixed effects: mathach ~ cses + meanses + sector + cses:meanses + cses:sector
Value Std.Error DF t-value p-value
(Intercept) 12.113821 0.1986555 7022 60.97905 0e+00
cses 2.935840 0.1507202 7022 19.47874 0e+00
meanses 5.342946 0.3690010 157 14.47949 0e+00
sector 1.214628 0.3061363 157 3.96761 1e-04
cses:meanses 1.044054 0.2910687 7022 3.58697 3e-04
cses:sector -1.642061 0.2331191 7022 -7.04387 0e+00
Correlation:
(Intr) cses meanss sector css:mn
cses 0.005
meanses 0.245 0.001
sector -0.697 -0.003 -0.356
cses:meanses 0.001 0.284 0.005 -0.002
cses:sector -0.003 -0.694 -0.002 0.005 -0.351
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-3.17005035 -0.72487781 0.01487615 0.75426738 2.96552094
Number of Observations: 7185
Number of Groups: 160
> pg339 <- lme(mathach ~ cses + meanses + sector + cses:meanses + cses:sector,
data=hsb12.big, random = ~ 1 | school, method="REML", corr=NULL)
> summary(pg339)
Linear mixed-effects model fit by REML
Data: hsb12.big
AIC BIC logLik
46520.79 46575.83 -23252.40
Random effects:
Formula: ~1 | school
(Intercept) Residual
StdDev: 1.541212 6.063506
Fixed effects: mathach ~ cses + meanses + sector + cses:meanses + cses:sector
Value Std.Error DF t-value p-value
(Intercept) 12.113821 0.1986485 7022 60.98117 0e+00
cses 2.935841 0.1507053 7022 19.48068 0e+00
meanses 5.342945 0.3689877 157 14.48001 0e+00
sector 1.214627 0.3061252 157 3.96775 1e-04
cses:meanses 1.044086 0.2910422 7022 3.58741 3e-04
cses:sector -1.642070 0.2330966 7022 -7.04459 0e+00
Correlation:
(Intr) cses meanss sector css:mn
cses 0.005
meanses 0.245 0.001
sector -0.697 -0.003 -0.356
cses:meanses 0.001 0.284 0.005 -0.002
cses:sector -0.003 -0.694 -0.002 0.005 -0.351
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-3.17006274 -0.72487339 0.01488267 0.75425043 2.96552444
Number of Observations: 7185
Number of Groups: 160
> anova(pg337, pg339)
Model df AIC BIC logLik Test L.Ratio p-value pg337 1 10 46524.79 46593.58 -23252.40 pg339 2 8 46520.79 46575.83 -23252.40 1 vs 2 0.003256624 0.9984
Second set of examples using willet data file.
To input the willett.csv data file in R, use the following code. Change d:/data/ to the where you stored the data files on your computer (and do use forward slashes to separate directory names)
> willett <- read.csv("d:/data/willett.csv", header=T)
Page 342
Let's describe the grouping structure.
> willettg <- groupedData( y ~ time | id, data=willett)
Page 342
Including both the intercept and the slope of time as random components.
Note: The intercept is included in the random argument even though it is not
explicitly listed in the random argument. In order to have only the
slope of time included as a random component the random argument would look
like this: random=~ time -1 and it would be the -1 which would indicate to R
that the intercept is not random.
> wfit1 <- lme(y ~ time, data=willettg, random = ~ 1 + time | id)
> summary(wfit1)
Linear mixed-effects model fit by REML
Data: willettg
AIC BIC logLik
1278.823 1296.386 -633.4114
Random effects:
Formula: ~ 1 + time | id
Structure: General positive-definite
StdDev Corr
(Intercept) 34.62336 (Inter
time 11.50654 -0.45
Residual 12.62843
Fixed effects: y ~ time
Value Std.Error DF t-value p-value
(Intercept) 164.3743 6.118849 104 26.86360 <.0001
time 26.9600 2.166604 104 12.44344 <.0001
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.272063 -0.6003579 0.02631387 0.6073068 1.584396
Number of Observations: 140
Number of Groups: 35
Page 344
#centering the covar variable
> willett$ccovar <- willett$covar - mean(willett$covar)
#grouping the data
> willettg1 <- groupedData( y ~ time | id, data=willett, outer=~ccovar)
> wfit2 <- lme(y ~ time*ccovar, data=willettg1, random = ~ time | id)
> summary(wfit2)
Linear mixed-effects model fit by REML
Data: willettg1
AIC BIC logLik
1276.285 1299.586 -630.1424
Random effects:
Formula: ~time | id
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 35.16281 (Intr)
time 10.35610 -0.489
Residual 12.62844
Fixed effects: y ~ time * ccovar
Value Std.Error DF t-value p-value
(Intercept) 164.37429 6.206120 103 26.485836 0.0000
time 26.96000 1.993878 103 13.521386 0.0000
ccovar -0.11355 0.504014 33 -0.225297 0.8231
time:ccovar 0.43286 0.161928 103 2.673155 0.0087
Correlation:
(Intr) time ccovar
time -0.522
ccovar 0.000 0.000
time:ccovar 0.000 0.000 -0.522
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.248168979 -0.618725671 0.004285166 0.614719013 1.556883000
Number of Observations: 140
Number of Groups: 35
Bottom of page 346 and top of page 347
> #compound symmetry > fit.gls.cs <- gls(y~time, data=willett, corr=corCompSymm(, form=~time | id)) > summary(fit.gls.cs)
Generalized least squares fit by REML
Model: y ~ time
Data: willett
AIC BIC logLik
1308.340 1320.049 -650.1698
Correlation Structure: Compound symmetry
Formula: ~time | id
Parameter estimate(s):
Rho
0.7061649
Coefficients:
Value Std.Error t-value p-value
(Intercept) 164.3743 5.774154 28.46725 0
time 26.9600 1.465864 18.39189 0
Correlation:
(Intr)
time -0.381
Standardized residuals:
Min Q1 Med Q3 Max
-2.47038764 -0.71517526 0.03314509 0.80858053 1.80988181
Residual standard error: 35.77345
Degrees of freedom: 140 total; 138 residual
> #autoregressive > fit.gls.ar1 <- gls(y~time, data=willett, corr=corAR1(, form=~time | id)) > summary(fit.gls.ar1)
Generalized least squares fit by REML
Model: y ~ time
Data: willett
AIC BIC logLik
1281.465 1293.174 -636.7327
Correlation Structure: AR(1)
Formula: ~time | id
Parameter estimate(s):
Phi
0.824912
Coefficients:
Value Std.Error t-value p-value
(Intercept) 164.33842 6.136374 26.78103 0
time 27.19786 1.919857 14.16661 0
Correlation:
(Intr)
time -0.469
Standardized residuals:
Min Q1 Med Q3 Max
-2.42825412 -0.71561366 0.03192972 0.78792580 1.76110724
Residual standard error: 36.37940
Degrees of freedom: 140 total; 138 residual
Table on page 347. Note that unstructured is omitted since this covariance structure is not provided. Note you need to multiply logLik by -2 to get the -2RLL value in the table.
> anova(fit.gls.cs, fit.gls.ar1)
Model df AIC BIC logLik fit.gls.cs 1 4 1308.340 1320.049 -650.1698 fit.gls.ar1 2 4 1281.465 1293.174 -636.7327
p. 348
> wfit2.ar1 <- lme( y ~ time*ccovar, data=willettg1, random = ~ 1 + time | id,
corr=corAR1(, form =~ time | id))
> summary(wfit2.ar1)
Linear mixed-effects model fit by REML
Data: willettg1
AIC BIC logLik
1278.045 1304.259 -630.0223
Random effects:
Formula: ~ 1 + time | id
Structure: General positive-definite
StdDev Corr
(Intercept) 35.46865 (Inter
time 10.53770 -0.488
Residual 11.88619
Correlation Structure: AR(1)
Formula: ~ 1 + time | id
Parameter estimate(s):
Phi
-0.1374558
Fixed effects: y ~ time * ccovar
Value Std.Error DF t-value p-value
(Intercept) 164.4229 6.198601 103 26.52581 <.0001
time 26.9080 1.978001 103 13.60363 <.0001
ccovar -0.1234 0.503403 33 -0.24517 0.8078
time:ccovar 0.4357 0.160638 103 2.71256 0.0078
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.318034 -0.6248881 -0.01432131 0.6548167 1.567537
Number of Observations: 140
Number of Groups: 35
UCLA Researchers are invited to our Statistical Consulting Services
We recommend others to our list of Other Resources for Statistical Computing Help
These pages are Copyrighted (c) by UCLA Academic Technology Services