Applied Longitudinal Data Analysis: Modeling Change and Event Occurrence by Judith D. Singer and John B. Willett Chapter 7: Examining the multilevel model’s error covariance structure
The examples on the page were written in R version 2.4.1.
Reading in the opposites data set and the nlme library. You may need to first install this library.
opposites <- read.table("https://stats.idre.ucla.edu/stat/r/examples/alda/data/opposites_pp.txt",header=TRUE,sep=",")
library(nlme)
Table 7.2, p. 246.
opp.reml <- lme(opp~time*ccog, opposites, random= ~time | id)
summary(opp.reml)
Linear mixed-effects model fit by REML
Data: opposites
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.16282 (Intr)
time 10.35609 -0.489
Residual 12.62843
Fixed effects: opp ~ time * ccog
Value Std.Error DF t-value p-value
(Intercept) 164.37429 6.206122 103 26.485828 0.0000
time 26.95998 1.993878 103 13.521383 0.0000
ccog -0.11355 0.504014 33 -0.225297 0.8231
time:ccog 0.43286 0.161928 103 2.673156 0.0087
Correlation:
(Intr) time ccog
time -0.522
ccog 0.000 0.000
time:ccog 0.000 0.000 -0.522
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.248169084 -0.618725724 0.004284978 0.614719613 1.556883051
Number of Observations: 140
Number of Groups: 35
Table 7.3, p. 258-259. Code contributed by Daniel B. Wright from the University of Sussex.
attach(opposites)
corandcov <- function(glsob,cov=T,...){
corm <- corMatrix(glsob$modelStruct$corStruct)[[5]]
print(corm)
varstruct <- print(glsob$modelStruct$varStruct)
varests <- coef(varstruct, uncons=F, allCoef=T)
covm <- corm*glsob$sigma^2*t(t(varests))%*%t(varests)
return(covm)}
unstruct <- gls(opp~time*ccog,opposites, correlation=corSymm(form = ~ 1 |id), weights=varIdent(form = ~ 1|wave),method="REML")
corandcov(unstruct)
[,1] [,2] [,3] [,4]
[1,] 1.0000000 0.8085045 0.7338888 0.4578986
[2,] 0.8085045 1.0000000 0.8626074 0.7187155
[3,] 0.7338888 0.8626074 1.0000000 0.7939959
[4,] 0.4578986 0.7187155 0.7939959 1.0000000
Variance function structure of class varIdent representing
1 2 3 4
1.0000000 0.9248170 0.9584917 0.9468611
1 2 3 4
1 1345.1224 1005.7731 946.1944 583.1998
2 1005.7731 1150.4649 1028.5350 846.5661
3 946.1944 1028.5350 1235.7723 969.2921
4 583.1998 846.5661 969.2921 1205.9640
comsym <- gls(opp~time*ccog,opposites, correlation=corCompSymm(,form = ~ 1 |id), method="REML") cc <- corMatrix(comsym$modelStruct$corStruct)[[5]] print(cc)
[,1] [,2] [,3] [,4] [1,] 1.0000000 0.7309599 0.7309599 0.7309599 [2,] 0.7309599 1.0000000 0.7309599 0.7309599 [3,] 0.7309599 0.7309599 1.0000000 0.7309599 [4,] 0.7309599 0.7309599 0.7309599 1.0000000
cc * comsym$sigma^2
[,1] [,2] [,3] [,4]
[1,] 1231.3559 900.0718 900.0718 900.0718
[2,] 900.0718 1231.3559 900.0718 900.0718
[3,] 900.0718 900.0718 1231.3559 900.0718
[4,] 900.0718 900.0718 900.0718 1231.3559
hetercom <- gls(opp~time*ccog,opposites, correlation=corCompSymm(,form = ~ 1 |id),weights=varIdent(form = ~1|wave), method="REML")
corandcov(hetercom)
[,1] [,2] [,3] [,4]
[1,] 1.0000000 0.7367232 0.7367232 0.7367232
[2,] 0.7367232 1.0000000 0.7367232 0.7367232
[3,] 0.7367232 0.7367232 1.0000000 0.7367232
[4,] 0.7367232 0.7367232 0.7367232 1.0000000
Variance function structure of class varIdent representing
1 2 3 4
1.0000000 0.8616618 0.8934397 0.9528415
1 2 3 4
1 1438.1404 912.9405 946.6096 1009.5465
2 912.9405 1067.7632 815.6573 869.8876
3 946.6096 815.6573 1147.9734 901.9689
4 1009.5465 869.8876 901.9689 1305.6977
auto1 <- gls(opp~time*ccog,opposites, correlation=corAR1(,form = ~ 1 |id), method="REML")
cc <- corMatrix(auto1$modelStruct$corStruct)[[5]]
print(cc)
[,1] [,2] [,3] [,4]
[1,] 1.0000000 0.8253423 0.6811899 0.5622148
[2,] 0.8253423 1.0000000 0.8253423 0.6811899
[3,] 0.6811899 0.8253423 1.0000000 0.8253423
[4,] 0.5622148 0.6811899 0.8253423 1.0000000
cc * auto1$sigma^2
[,1] [,2] [,3] [,4]
[1,] 1256.6859 1037.1960 856.0417 706.5274
[2,] 1037.1960 1256.6859 1037.1960 856.0417
[3,] 856.0417 1037.1960 1256.6859 1037.1960
[4,] 706.5274 856.0417 1037.1960 1256.6859
hauto1 <- gls(opp~time*ccog,opposites, correlation=corAR1(,form = ~ 1 |id), weights=varIdent(form = ~1|wave), method="REML")
corandcov(hauto1)
[,1] [,2] [,3] [,4]
[1,] 1.0000000 0.8198784 0.6722005 0.5511227
[2,] 0.8198784 1.0000000 0.8198784 0.6722005
[3,] 0.6722005 0.8198784 1.0000000 0.8198784
[4,] 0.5511227 0.6722005 0.8198784 1.0000000
Variance function structure of class varIdent representing
1 2 3 4
1.0000000 0.9104126 0.9512871 0.9593613
1 2 3 4
1 1340.7078 1000.7413 857.3232 708.8668
2 1000.7413 1111.2471 951.9922 787.1427
3 857.3232 951.9922 1213.2696 1003.1765
4 708.8668 787.1427 1003.1765 1233.9528
toep <- gls(opp~time*ccog,opposites, correlation=corARMA(,form = ~ 1 |id,p=3,q=0), method="REML")
cc <- corMatrix(toep$modelStruct$corStruct)[[5]]
print(cc)
[,1] [,2] [,3] [,4]
[1,] 1.0000000 0.8255072 0.7190556 0.5004861
[2,] 0.8255072 1.0000000 0.8255072 0.7190556
[3,] 0.7190556 0.8255072 1.0000000 0.8255072
[4,] 0.5004861 0.7190556 0.8255072 1.0000000
cc * toep$sigma^2
[,1] [,2] [,3] [,4]
[1,] 1246.8832 1029.3111 896.5783 624.0477
[2,] 1029.3111 1246.8832 1029.3111 896.5783
[3,] 896.5783 1029.3111 1246.8832 1029.3111
[4,] 624.0477 896.5783 1029.3111 1246.8832
anova(unstruct,comsym,hetercom,auto1,hauto1,toep)
Model df AIC BIC logLik Test L.Ratio p-value unstruct 1 14 1283.789 1324.566 -627.8944 comsym 2 6 1299.048 1316.524 -643.5238 1 vs 2 31.258855 0.0001 hetercom 3 9 1302.954 1329.168 -642.4770 2 vs 3 2.093731 0.5532 auto1 4 6 1277.876 1295.352 -632.9382 3 vs 4 19.077497 0.0003 hauto1 5 9 1282.840 1309.054 -632.4199 4 vs 5 1.036723 0.7924 toep 6 8 1274.081 1297.382 -629.0404 5 vs 6 6.759007 0.0093
Table 7.4, p. 265.
#Standard error covariance structure
summary(lme(opp ~ time * ccog, opposites, random = ~ time | id))
Linear mixed-effects model fit by REML
Data: opposites
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.16282 (Intr)
time 10.35609 -0.489
Residual 12.62843
Fixed effects: opp ~ time * ccog
Value Std.Error DF t-value p-value
(Intercept) 164.37429 6.206122 103 26.485828 0.0000
time 26.95998 1.993878 103 13.521383 0.0000
ccog -0.11355 0.504014 33 -0.225297 0.8231
time:ccog 0.43286 0.161928 103 2.673156 0.0087
Correlation:
(Intr) time ccog
time -0.522
ccog 0.000 0.000
time:ccog 0.000 0.000 -0.522
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.248169084 -0.618725724 0.004284978 0.614719613 1.556883051
Number of Observations: 140
Number of Groups: 35
#Unstructured error covariance structure
summary(gls(opp~time*ccog,opposites, correlation=corSymm(form = ~ 1 |id), weights=varIdent(form = ~ 1|wave),method="REML"))
Generalized least squares fit by REML
Model: opp ~ time * ccog
Data: opposites
AIC BIC logLik
1283.789 1324.566 -627.8944
Correlation Structure: General
Formula: ~1 | id
Parameter estimate(s):
Correlation:
1 2 3
2 0.809
3 0.734 0.863
4 0.458 0.719 0.794
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | wave
Parameter estimates:
1 2 3 4
1.0000000 0.9248170 0.9584917 0.9468610
Coefficients:
Value Std.Error t-value p-value
(Intercept) 165.83211 5.953009 27.856855 0.0000
time 26.58438 1.925888 13.803698 0.0000
ccog -0.07409 0.483458 -0.153250 0.8784
time:ccog 0.45829 0.156406 2.930157 0.0040
Correlation:
(Intr) time ccog
time -0.508
ccog 0.000 0.000
time:ccog 0.000 0.000 -0.508
Standardized residuals:
Min Q1 Med Q3 Max
-2.45633281 -0.74687031 0.06638661 0.72058752 2.16323956
Residual standard error: 36.67591
Degrees of freedom: 140 total; 136 residual