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