10.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.

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.506836 2.638686
## 
## 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.827573 3.461330 
## Fixed effects:  verb ~ 1 + time0 
##                 Value Std.Error  DF  t-value p-value
## (Intercept) 25.415343 0.4275323 406 59.44660       0
## time01       7.192402 0.3374451 406 21.31428       0
## time02      18.334559 0.5249722 406 34.92482       0
##  Correlation: 
##        (Intr) time01
## time01 -0.118       
## time02  0.221  0.507
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.24169075 -0.62480430 -0.06272641  0.63543916  3.06613520 
## 
## Number of Observations: 612
## Number of Groups: 204

10.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

10.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.325248 5.506836
## Residual     6.962664 2.638686

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.32525 30.32525 30.32525
## [2,] 30.32525 30.32525 30.32525
## [3,] 30.32525 30.32525 30.32525

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.962664  3.499314 17.08695
## [2,]  3.499314 23.255458 31.93237
## [3,] 17.086955 31.932371 83.41832

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

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

##          [,1]     [,2]      [,3]
## [1,] 37.28791 33.82456  47.41220
## [2,] 33.82456 53.58071  62.25762
## [3,] 47.41220 62.25762 113.74357

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

meanvector - meanvector_timecatunst

##              [,1]
## [1,] 1.598721e-13
## [2,] 7.105427e-14
## [3,] 1.563194e-13

varcovmatrix - varcovmatrix_timecatunst

##               verb2         verb4         verb6
## verb2 -6.906778e-05 -4.991762e-03 -0.0073203738
## verb4 -4.991762e-03 -2.480221e-06 -0.0027309696
## verb6 -7.320374e-03 -2.730970e-03 -0.0002526546

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