R Textbook Examples Applied Long

[http://www.ats.ucla.edu/stat/_headers/header1.htm][http://www.ats.ucla.edu/stat/r/examples/alda/header.htm][http://www.ats.ucla.edu/stat/_headers/header2.htm]

R Textbook Examples
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("http://www.ats.ucla.edu/stat/r/examples/alda/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)
  varests <- print(glsob$modelStruct$varStruct)  
  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(lme(opp ~ time * ccog, opposites, random =  ~ time | id, correlation = corSymm(, form =  ~ wave | id)))
	 
Linear mixed-effects model fit by REML
 Data: opposites 
       AIC      BIC    logLik
  1284.166 1324.944 -628.0832

Random effects:
 Formula: ~time | id
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev    Corr  
(Intercept) 32.089271 (Intr)
time         6.384574 -0.342
Residual    16.523913       

Correlation Structure: General
 Formula: ~wave | id 
 Parameter estimate(s):
 Correlation: 
  1      2      3     
2  0.182              
3  0.063  0.470       
4 -0.859  0.059  0.163
Fixed effects: opp ~ time * ccog 
                Value Std.Error  DF   t-value p-value
(Intercept) 165.47581  5.975949 103 27.690300  0.0000
time         26.56958  1.926463 103 13.791899  0.0000
ccog         -0.04454  0.485321  33 -0.091774  0.9274
time:ccog     0.45420  0.156453 103  2.903115  0.0045
 Correlation: 
          (Intr) time   ccog  
time      -0.507              
ccog       0.000  0.000       
time:ccog  0.000  0.000 -0.507

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-2.12631110 -0.66715636 -0.08334923  0.61274376  2.16761303 

Number of Observations: 140
Number of Groups: 35

[http://www.ats.ucla.edu/stat/r/footer.htm]

R Textbook Examples 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

R Textbook Examples
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