Different ways to include pre-test performance as a covariate in a linear mixed-effect regression. Which is correct?

Question

I have a pre-post experimental design, where I have measured participants' performance in three courses (tasks; A, B, C) at both pre and post-test. The study is a learning experiment. 60 elite alpine skiers have practised 12 runs in each of the three courses (A, B, C). One group have practised these courses in 'random' order, e.g., 'ABC, whereas the other group have practised these courses in a blocked order, e.g. 'AAA, BBB, CCC. I want to build a model to test whether the 'random' group performed better on post-test when controlling for pre-test performance. So the pre-test is just an earlier version of their performance. The data structure looks like this:

subjectID	Course	Group	pre-test.score	post-test.score
1	A	ii	###	#
1	B	ii	###	#
1	C	ii	###	#
2	A	b	###	#
2	B	b	###	#
2	C	b	###	#

I have analysed these data using a linear mixed-effect regression model where I predict post-test performance, controlling for pre-test performance with the + sign:

# I fit these models with lmer in R
CI_post <- lmer(
  post ~ 
    pre +
    group * course 
  + (1|subjectID) , 
  data = dat, 
  REML = FALSE)

Using Satterthwaite's method from the emmeans package I get:

CI_post_interaction_coursegroup <- emmeans(CI_post, specs = c("course", "group"),lmer.df = "satterthwaite")

course group       emmean    SE   df lower.CL upper.CL
 A      blocked      0.311 0.191 6.65  -0.1452    0.768
 B      blocked      0.649 0.180 5.38   0.1954    1.102
 C      blocked      1.141 0.195 7.28   0.6847    1.598
 A      interleaved  0.189 0.194 7.15  -0.2666    0.645
 B      interleaved  0.497 0.179 5.31   0.0451    0.949
 C      interleaved  1.046 0.191 6.72   0.5907    1.502

But I could perhaps also perform the same model adding pre-test with as an interaction term to the model, so that that the model becomes pre.test * course * group

CI_post <- lmer(
  post ~ 
    pre *
    group * course 
  + (1|subjectID),
  data = dat, 
  REML = FALSE)

, which gives very different estimates:

 course group        emmean    SE    df lower.CL upper.CL
 A      blocked     -0.0669 0.188 11.10   -0.481    0.347
 B      blocked      0.6466 0.161  6.09    0.255    1.038
 C      blocked      1.1980 0.194 12.65    0.778    1.618
 A      interleaved -0.1520 0.211 16.76   -0.597    0.293
 B      interleaved  0.4872 0.160  6.12    0.098    0.876
 C      interleaved  1.0593 0.181  9.82    0.654    1.464

I am trying to understand the exact differences between these two models, and which is them is "correct"?. Long (see comments below) gave a helpful comment that "if you want the group:course interaction to vary depending on the value of pre then you fit the 2nd model with the 3-way interaction". But since my pre-test score is a continuous variable that measured each participant's performance on course A, course B and course C, will the model use this data structure when I add an interaction term between the covariate and factors in the model?

Answer 1

A wide variety of interaction terms could be "correct" here, as comments from Robert Long and COOLSerdash suggest. To decide which if any interaction to use involving the pre value, you have to think very carefully about how you want to "control for pre-test performance."

Try pre-post difference as outcome

The first question is whether you might be better off just using the time difference pre-post as the outcome measure of improvement following training. That's one perfectly acceptable way to control for pre-test performance. Or, if you expect proportional improvements rather than linear improvements in time with training, use the difference of the logs of pre and post to get the log of the pre/post ratio. (You shouldn't then make the mistake* of continuing to include the pre value in your model.) That's the idea behind a paired t-test, just expanded in your situation to include the types of courses and types of training. Expressing the difference in this direction has an advantage in that an increasingly positive outcome value is an improvement in performance, which seems to be of most interest.

That would change the interpretation of the fixed and random effects in your model, but maybe in a way more useful in this application. The course*group term would give you how the improvement itself depends on each combination of group and course. The (1|subjectID) would represent estimated differences in improvement among individuals at some baseline combination of group and course.

That approach, if it's reasonable in your application, would greatly simplify both the modeling and explaining results to others.

Using pre as a predictor for post, without further interactions.

Although that's another acceptable way to "control for pre-test performance," it does so in a different way. Your model then is something like:

post ~ b0 + b1 * pre + (terms from course, group, their interactions, with associated coefficients)

That implies that the post value is based around a common intercept b0 plus an additive contribution proportional to the pre value (the same regardless of the particular combination of group and course in question) and separate additive contributions from the combination of group and course. Does that make sense for your application?

Including pre in interactions

That would take the prior model and expand it to allow the proportional contribution of pre to post to differ depending on course and group (and vice-versa). You don't need to include all potential interactions with a pre*course*group term. That full interaction allows each of the course*group interactions to depend linearly on the value of pre. Is that what you need? Is such a linear interaction with the pre value reasonable? Might simply allowing for interactions of pre with course and/or group individually (e.g., with pre:course and/or pre:group terms in R syntax) be enough? Again, there is no "correct" answer in general here, just what makes sense based on your understanding of the subject matter.

Recommendation

In general, it can help to write out the full equations for the models including all coefficients and terms rather than just looking at the shorthand of R (or other modeling) syntax. If you do that, you will see just which contributions to outcome are being modified by including particular interaction terms. You can then choose accordingly.

If these are elite Alpine skiers, you presumably are looking at post-pre differences on the order of a few seconds, around full pre times on the order of many seconds to minutes. In that case, my guess is that the simple t-test-like approach of just using pre-post as the outcome while omitting pre as a predictor otherwise will work quite well and be much easier to explain to others.

Note that nothing above allows for wipe-outs in pre or post values or for considering changes in snow-surface or other environmental conditions. That's true for all of these approaches, however.

---

*An answer from Robert Long on the linked page shows a way to model properly if you think that the change in performance depends on the pre value. In outline, you model the pre and post values separately, with an indicator variable predictor specifying which of pre or post is involved.