How to test whether the effects of predictors on one DV explains their effects on another in a mixed ANOVA in R

Question

I have data from a psychology experiment in which instructional method was manipulated between subjects via two factors: pretraining (3 levels) and training (2 levels). I assessed participants' performance via two tests, a pretest and a posttest, so I have test accuracy as my dv with test section (pre/post) as a within-subjects factor. I also have some other factors of secondary interest. I analyzed this data using a mixed ANOVA and found a significant (predicted) 3-way interaction of pretraining, training, and test section, indicating that the relative advantages of the two training conditions, as measured by improvement from pretest to posttest, was different depending on pretraining condition.

Now I have noticed that accuracy score in answering questions DURING the training session, which I'll call "study", is also affected by the pretraining factor, i.e. one of the pretraining conditions increases study accuracy relative to the others. I'd like to test whether this effect on study accuracy can explain the subsequent effect of pretraining on the training * test section interaction. In fact I suspect it cannot, but I'd like to test my belief. How can I do this?

My first thought was to use ANCOVA, i.e. add study accuracy as a covariate to my original ANOVA and then see whether the pretrainingtrainingsection interaction is still significant. However, in another thread (here) I was advised not to use ANCOVA when my experimentally manipulated factors are correlated with the covariate, as certainly is the case here (i.e. there is an effect of pretraining on study accuracy). So, what should I do instead?

# Sample data:
nsubj = 215; nsec  = 2; nprob = 6
D = data.frame(
    subjid   = rep( 1:nsubj, each=nsec*nprob ),
    pretrain = rep( sample( c('a','b','c'), nsubj, replace=TRUE ), each=nsec*nprob ),
    training = rep( sample( c('j','k'), nsubj, replace=TRUE ), each=nsec*nprob ),
    study    = rep( sample( 1:6, nsubj, replace=TRUE ), each=nsec*nprob ),
    section  = rep( rep( c( 'pretest', 'posttest' ), each=nprob ), nsubj ),
    probtype = rep( c( 'v', 'w', 'x', 'y', 'z', 'z' ), nsec*nsubj ),
    accuracy = sample( c( 0.0, 0.5, 1.0 ), nsubj * nsec * nprob, replace=TRUE ) )

# Model:
library( afex )
ez.glm( "subjid", "accuracy", D, within=c("section","probtype"), between= c("pretrain","training"), type=3 )

# Model with covariate:
ez.glm( "subjid", "accuracy", D, within=c("section","probtype"), between= c("pretrain","training"), covariate="study", type=3 )

Answer 1

As the Miller & Chapman (2001) paper explains in great detail, using a covariate that is correlated with the independent variables is usually a bad idea. Hence simply adding study to the model will lead to wrong conclusions.

Note however, that in your case only pretraining is correlated with the dv, so you could still use it as a covariate for the training sessions in the following steps:

1. Compute a regression with the posttest scores as dv and only study as iv.
2. Take the residuals of the regression and add the value of the intercept (both from the model computed in step 1).
3. Use the values from step 2 as the new posttest scores (those should be now on the same scale as the pretest scores, but the influence of study should be removed) and run the same ANOVA as your first one.
4. Note that your df for posttest are now one off, so you need to take this manually into account by subtracting 1 from all effects involving posttest (note that this most likely is not trivial as you also need to correctly apply the df-corrections, so for all effects involving probtype you cannot simply subtract one, but you need to subtract one from the uncorrected df and then need to obtain the correct Greenhouse-Geisser epsilon and multiply the new df with it).

As the saying goes, the implementation of this is left as an exercise to the reader.