12.7 Repeated Measures MANOVA (Unstructured)

Let’s adjust the RM MANOVA error structure to have no constraints. To do this in lme() we can use the correlation argument.

We can now use the weights argument to indicate that the residual variances will be a function of time.

Option correlation = corSymm(form=~1|id) specifies that the correlation structure is unstructured.

verblong$time0 <- factor(verblong$time0, ordered=FALSE)
timecatunst_fit <- lme(
  fixed= verb ~ 1 + time0, 
  random= ~ 1|id,
  weights=varIdent(form=~1|time0),
  correlation=corSymm(form=~1|id),
  data=verblong,
  na.action = na.exclude
)
summary(timecatunst_fit)

## Linear mixed-effects model fit by REML
##   Data: verblong 
##       AIC      BIC   logLik
##   3869.06 3913.178 -1924.53
## 
## Random effects:
##  Formula: ~1 | id
##         (Intercept) Residual
## StdDev:    5.506713 2.638939
## 
## Correlation Structure: General
##  Formula: ~1 | id 
##  Parameter estimate(s):
##  Correlation: 
##   1     2    
## 2 0.275      
## 3 0.709 0.725
## Variance function:
##  Structure: Different standard deviations per stratum
##  Formula: ~1 | time0 
##  Parameter estimates:
##        0        1        2 
## 1.000000 1.827452 3.461025 
## Fixed effects:  verb ~ 1 + time0 
##                 Value Std.Error  DF  t-value p-value
## (Intercept) 25.415343 0.4275322 406 59.44662       0
## time01       7.192402 0.3374453 406 21.31428       0
## time02      18.334559 0.5249720 406 34.92483       0
##  Correlation: 
##        (Intr) time01
## time01 -0.118       
## time02  0.221  0.507
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.24169387 -0.62481099 -0.06273486  0.63549242  3.06610057 
## 
## Number of Observations: 612
## Number of Groups: 204

12.7.1 Model-Implied Mean Vector

There will be no differences in the model implied mean vector in this model.

beta <- matrix(fixef(timecatunst_fit), nrow=3, ncol=1)
X <- matrix(c(1,1,1,0,1,0,0,0,1), nrow=3,ncol=3)
meanvector_timecatunst  <- X %*% beta
meanvector_timecatunst

##          [,1]
## [1,] 25.41534
## [2,] 32.60775
## [3,] 43.74990

12.7.2 Model-Implied Covariance Matrix

Again, we start with the between-part of our model.

#parsing the model variances & covariances
VarCorr(timecatunst_fit)

## id = pdLogChol(1) 
##             Variance  StdDev  
## (Intercept) 30.323891 5.506713
## Residual     6.963999 2.638939

Psi <- matrix(c(as.numeric(VarCorr(timecatunst_fit)[1])),nrow=1,ncol=1)  
Z <- matrix(c(1,1,1),  nrow=3,ncol=1)
Z %*% Psi %*% t(Z)

##          [,1]     [,2]     [,3]
## [1,] 30.32389 30.32389 30.32389
## [2,] 30.32389 30.32389 30.32389
## [3,] 30.32389 30.32389 30.32389

Now it gets a bit more complicated, because we have a fully unstructured matrix for the within-person residual variance-covariance.

#parse the residual/"error" variance-covariance
sigma <- as.numeric(VarCorr(timecatunst_fit)[4]) #note this is standard deviation

#From the summary(timecatunst_fit) above
#  Variance function:
#   Structure: Different standard deviations per stratum
# Formula: ~1 | time0 
# Parameter estimates:
#        0        1        2 
# 1.000000 1.827573 3.461330 
sigmaunst <- sigma * (diag(c(1.000000, 1.827573, 3.461330), nrow=3,ncol=3))

#From the summary(timecatunst_fit) above
# Correlation Structure: General
# Formula: ~1 | id 
# Parameter estimate(s):
#   Correlation: 
#   1     2    
# 2 0.275      
# 3 0.709 0.725
cormatrixunst <- matrix(c(1.000, 0.275, 0.709,
                          0.275, 1.000, 0.725,
                          0.709, 0.725, 1.000),
                        nrow=3,ncol=3)

#Pre and post multiply by SDs to convert Correlation matrix into
#Covariance matrix
covresidunst <- sigmaunst %*% cormatrixunst %*% t(sigmaunst)
covresidunst

##           [,1]      [,2]     [,3]
## [1,]  6.963999  3.499985 17.09023
## [2,]  3.499985 23.259917 31.93850
## [3,] 17.090231 31.938495 83.43432

Finally, calculate the implied between- + within-person variance-covariances

varcovmatrix_timecatunst <- Z %*% Psi %*% t(Z) + covresidunst
varcovmatrix_timecatunst

##          [,1]     [,2]      [,3]
## [1,] 37.28789 33.82388  47.41412
## [2,] 33.82388 53.58381  62.26239
## [3,] 47.41412 62.26239 113.75821

Let’s look at the misfit (real matrix - model implied).

meanvector - meanvector_timecatunst

##               [,1]
## [1,] -1.563194e-13
## [2,] -2.131628e-13
## [3,] -1.421085e-13

varcovmatrix - varcovmatrix_timecatunst

##               verb2        verb4        verb6
## verb2 -0.0000473069 -0.004305830 -0.009240161
## verb4 -0.0043058296 -0.003105210 -0.007497682
## verb6 -0.0092401614 -0.007497682 -0.014892895

It appears we have fully reproduced our observed data. We can formally test the improvement of this model using the anove() function.

anova(timecathet_fit, timecatunst_fit)

##                 Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## timecathet_fit      1  7 3930.253 3961.136 -1958.127                        
## timecatunst_fit     2 10 3869.060 3913.178 -1924.530 1 vs 2 67.19366  <.0001