UCLA Academic Technology Services HomeServicesClassesContactJobs
Stat Computing > R > Textbook Examples > Applied Longitudinal Data Analysis

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

How to cite this page

Report an error on this page

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

The content of this web site should not be construed as an endorsement of any particular web site, book, or software product by the University of California.