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.506716 2.638933
## 
## 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.827454 3.461032 
## 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.24169379 -0.62481084 -0.06273466  0.63549118  3.06610136 
## 
## 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.323922 5.506716
## Residual     6.963968 2.638933

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.32392 30.32392 30.32392
## [2,] 30.32392 30.32392 30.32392
## [3,] 30.32392 30.32392 30.32392

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.963967  3.499969 17.09015
## [2,]  3.499969 23.259812 31.93835
## [3,] 17.090154 31.938350 83.43394

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.82389  47.41408
## [2,] 33.82389 53.58373  62.26227
## [3,] 47.41408 62.26227 113.75786

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

meanvector - meanvector_timecatunst

##              [,1]
## [1,] 1.421085e-13
## [2,] 2.629008e-13
## [3,] 2.557954e-13

varcovmatrix - varcovmatrix_timecatunst

##               verb2        verb4        verb6
## verb2 -4.663967e-05 -0.004320914 -0.009193447
## verb4 -4.320914e-03 -0.003030440 -0.007383449
## verb6 -9.193447e-03 -0.007383449 -0.014544496

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