12.6 Repeated Measures MANOVA

Now let’s adjust the RM ANOVA error structure to get a RM MANOVA. The multivariate part of the MANOVA, has to do with relaxing the assumption on the error structure - this means more flexibility than compound symmetry.

Here is our RM ANOVA model with time as within-person factor. The error structure is compound symmetry with heterogeneous variances.

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

summary(timecathet_fit)

## Linear mixed-effects model fit by REML
##   Data: verblong 
##        AIC      BIC    logLik
##   3930.253 3961.136 -1958.127
## 
## Random effects:
##  Formula: ~1 | id
##         (Intercept) Residual
## StdDev:    6.086682 2.871709
## 
## Variance function:
##  Structure: Different standard deviations per stratum
##  Formula: ~1 | time0 
##  Parameter estimates:
##        0        1        2 
## 1.000000 1.254186 2.341456 
## Fixed effects:  verb ~ 1 + time0 
##                 Value Std.Error  DF  t-value p-value
## (Intercept) 25.415343 0.4712021 406 53.93724       0
## time01       7.192402 0.3225105 406 22.30130       0
## time02      18.334559 0.5119102 406 35.81597       0
##  Correlation: 
##        (Intr) time01
## time01 -0.266       
## time02 -0.168  0.245
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.24995336 -0.55415412 -0.05694858  0.51324627  2.99087738 
## 
## Number of Observations: 612
## Number of Groups: 204

12.6.1 Model-Implied Mean Vector

So what is the implied representation of the basic information Making the implied mean vector.

fixef(timecathet_fit)

## (Intercept)      time01      time02 
##   25.415343    7.192402   18.334559

beta <- matrix(
  c(
    fixef(timecathet_fit)[1], 
    fixef(timecathet_fit)[2], 
    fixef(timecathet_fit)[3]
    ), nrow=3, ncol=1)
beta

##           [,1]
## [1,] 25.415343
## [2,]  7.192402
## [3,] 18.334559

Create the model design matrix, X, for the fixed effects. In this model this is a matrix of order \(3 \times 3\).

X <- matrix(c(1,1,1,0,1,0,0,0,1), nrow=3,ncol=3)
X

##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    1    1    0
## [3,]    1    0    1

Creating the model implied mean vector through multiplication

\[ \mathbf{Y} = \mathbf{X}\mathbf{\beta} + 0 + 0 \]

meanvector_timecathet  <- X %*% beta
meanvector_timecathet

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

See the differences in the means across levels of time0. This is exactly the same as in the last model

12.6.2 Model-Implied Covariance Matrix

#parse the between-person variances
Psi <- matrix(c(as.numeric(VarCorr(timecathet_fit)[1])),nrow=1,ncol=1)  
Psi

##         [,1]
## [1,] 37.0477

#create the model design matrix
Z <- matrix(c(1,1,1), nrow=3,ncol=1)
Z

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

So the implied variance covariance of the between-person random effects for the three repeated measures is:

Z %*% Psi %*% t(Z)

##         [,1]    [,2]    [,3]
## [1,] 37.0477 37.0477 37.0477
## [2,] 37.0477 37.0477 37.0477
## [3,] 37.0477 37.0477 37.0477

Now for the within-person residual var-cov.

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

## [1] 2.871709

Note: Now we have another step due to the heterogeneous variances. We have to get the heterogeneous weights for the residual standard deviation from the summary(timecathet_fit).

# Variance function:
#   Structure: Different standard deviations per stratum
# Formula: ~1 | time0 
# Parameter estimates:
#        0        1        2 
# 1.000000 1.254186 2.341456 

Pulling from the output.

sigmahet <- sigma * (diag(c(1.000000, 1.254186, 2.341456), nrow=3,ncol=3))
#Calculate the implied residual error variance-covariance
Lambda   <- sigmahet^2 

Finally calculating the model implied between- + within var-cov structure.

varcovmatrix_timecathet <- Z %*% Psi %*% t(Z) + Lambda
varcovmatrix_timecathet

##          [,1]     [,2]     [,3]
## [1,] 45.29441 37.04770 37.04770
## [2,] 37.04770 50.01964 37.04770
## [3,] 37.04770 37.04770 82.25961

Note the heterogeneity we now have along the diagonal.

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

meanvector - meanvector_timecathet

##               [,1]
## [1,] -5.329071e-14
## [2,] -4.973799e-14
## [3,]  1.136868e-13

varcovmatrix - varcovmatrix_timecathet

##           verb2     verb4    verb6
## verb2 -8.006572 -3.228132 10.35718
## verb4 -3.228132  3.561067 25.20719
## verb6 10.357180 25.207186 31.48370

Getting better. Note, we can formally compare the models.

anova(timecat_fit, timecathet_fit)

##                Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## timecat_fit        1  5 4013.176 4035.235 -2001.588                        
## timecathet_fit     2  7 3930.253 3961.136 -1958.127 1 vs 2 86.92304  <.0001

Remember, the null hypothesis of the LRT states that the more constrained model provides as good a fit for the data as the less constrained model. If the null hypothesis is rejected, then the alternative, unconstrained model provides a significant improvement in fit over the smaller model. Thus, allowing for the additional heterogeneity improved our model.