0

I am using the lme4 package in R to undertake linear mixed effect models (LMM). Essentially all participants received two interventions (an intervention treatment and a placebo (control)) and were separated by a washout period. However, the order or sequence they received the interventions differed.

enter image description here

An interaction term of intervention and visit was included in the LMM with eight levels including all combinations of intervention (2 levels: control and intervention) and visit (4 levels: visit 1=baseline 1, visit 2, visit 3=post-randomization baseline 2, visit 4).

My question is how do I determine the intervention effect by a post-hoc t-test as the average differences of the differences between interventions, hence between visits 1 and 2 and between visits 3 and 4. I also want to determine the effects of the intervention and control compared to baseline.

Please see code below:

  model1<- lmer(X ~ treatment_type:visit_code + (1|SID) + (1|SID:period), na.action= na.omit, data = data.x)
emm <- emmeans(model1 , ~treatment_type:visit_code)

My results of model 1 is:

emm
 treatment_type visit_code  emmean    SE   df lower.CL upper.CL
 Control        T0         -0.2915 0.167 26.0   -0.635   0.0520
 Intervention   T0         -0.1424 0.167 26.0   -0.486   0.2011
 Control        T1         -0.2335 0.167 26.0   -0.577   0.1100
 Intervention   T1          0.0884 0.167 26.0   -0.255   0.4319
 Control        T2          0.0441 0.167 26.0   -0.299   0.3876
 Intervention   T2         -0.2708 0.168 26.8   -0.616   0.0748
 Control        T3          0.1272 0.167 26.0   -0.216   0.4708
 Intervention   T3          0.0530 0.168 26.8   -0.293   0.3987

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95

I first created a matrix/ vectors: #name vectors

Control.B1<- c(1,0,0,0,0,0,0,0) #control baseline 1 (visit 1)
Intervention.B1<- c(0,1,0,0,0,0,0,0)  #intervention baseline 1 (visit 1)
Control.A2<- c(0,0,1,0,0,0,0,0) #post control 1 (visit 2)
Intervention.A2<- c(0,0,0,1,0,0,0,0) #post intervention 1 (visit 2)
ControlB3<- c(0,0,0,0,1,0,0,0) #control baseline 2 (visit 3)
Intervention.B3<- c(0,0,0,0,0,1,0,0) #intervention baseline 2 (visit 3)
Control.A4<- c(0,0,0,0,0,0,1,0) #post control 2 (visit 4)
Intervention.A4<- c(0,0,0,0,0,0,0,1) #post intervention 2 (visit 4)

Contbaseline = (Control.B1 + Control.B3)/2 # average of control baseline visits
Intbaseline = (Intervention. B1 + Intervention.B3)/2 # average of intervention baseline visits
ControlAfter= (Control.A2 + Control.A4)/2 # average of after control visits
IntervAfter= (Intervention.A2 + Intervention.A4)/2 # average of after intervention visits

Control.vs.Baseline = (ControlAfter-Contbaseline) 
Intervention.vs.Baseline = (IntervAfter-Intbaseline) 
Control.vs.Intervention = ((Control.vs.Baseline)-(Intervention.vs.Baseline))

the output of these are as follows:

> Control.vs.Baseline
[1] -0.5  0.0  0.5  0.0 -0.5  0.0  0.5  0.0
> Intervention.vs.Baseline
[1]  0.0 -0.5  0.0  0.5  0.0 -0.5  0.0  0.5
> Control.vs.Intervention
[1] -0.5  0.5  0.5 -0.5 -0.5  0.5  0.5 -0.5

Is this correct to the average differences of the differences between baseline and treatment periods?

Many thanks in advance!

Jess H
  • 39
  • 7
  • Assuming you subsequently gather these into a named list `cons` and run `contrast(emm, cons)` these look right. There is a disconnect between variable names used in the different steps -- e.g., `Control.B1` vs. `ContB1`. – Russ Lenth May 12 '22 at 13:15
  • In the context of a crossover design that you showed in an earlier version of this Q, , am I correct that visits 1 and 2 are period 1 and visits 3 and 4 are period 2? I think there may be another way to model this that makes this less complicated. – Russ Lenth May 12 '22 at 13:19
  • Hi Russ, thank you so much for your reply. That is correct. So essentially visits 1 is the baseline of period one, and visits 3 is the baseline of period two if you like. However, participants are either receiving the control or intervention treatment in period one, then they swap to the alternative treatment in period two. So still two different treatment groups in each period. I have edited my question and re-added the figure incase its useful. – Jess H May 12 '22 at 22:10
  • Regarding your first question, yes, then I each as: `contvstreat<-data.frame(contrast(emm, method = list(Control.vs.Intervention), adjust ="bonferroni")) # contrasts control vs treatment` – Jess H May 12 '22 at 22:30
  • No this will not work. The method has to be a NAMED list, e.g. `list(CvsB = Control.vs.Baseline, IvS = Intervention.vs.Baseline, CvI = Control.vs.Intervention)` – Russ Lenth May 12 '22 at 23:34
  • Thank you very much Russ, so you're saying I should run all these post hocs together, rather than separately, and include a named list as per your suggestion above? – Jess H May 12 '22 at 23:46

1 Answers1

1

A two-period crossover is the same as a repeated 2x2 Latin square. My suggestion for future such experiments is to structure the data accordingly, using variables for sequence (rows), period (columns), and treatment (assigned in the pattern (A,B) first sequence and (B,A) second sequence. The subjects are randomized to which sequence they are in.

So with your data, you would need to add a variable sequence that has the level AB for those subjects who receive the treatment sequence A, A, B, B, and level BA for those who receive B, B, A, A (though I guess the 1st and 3rd are really baseline for everybody).

Since there are 4 visits, it helps keep things sorted if you recode that as two factors trial and period, as follows:

visit   trial   period
  1      base      1
  2      test      1
  3      base      2
  4      test      2

Then fit the model with formula

model2 <- lmer(X ~ (sequence + period + treatment_type) * trial + 
                   (1|SID:sequence), ...etc...)

The parenthesized part is the standard model for a Latin square. Then the analysis can be done without custom contrasts as follows:

RG <- ref_grid(model2)   # same really as emmeans() for all 4 factors
CHG <- contrast(RG, "consec", simple = "trial")
CHG <- update(CHG, by = NULL, infer = c(TRUE, FALSE))

CHG contains the differences from baseline (trial differences for each combination of the other three factors. The update() step removes the by variables saved from contrast(). Now, we can get the marginal means and comparisons for each factor:

emmeans(CHG, consec ~ treatment_type)
emmeans(CHG, consec ~ period)
emmeans(CHG, consec ~ sequence)

These will be the same results you got the other way via custom contrasts. The one that was a difference of differences before is now handled by sequence. This works because in a 2x2 Latin square, the main effect of each factor is confounded with the two-way interaction of the other two factors.

Russ Lenth
  • 5,922
  • 2
  • 13
  • 21
  • Thank you so much Russ, this is the best explanation I've ever received! I ran your code, however, I got different results when using the means. I am assuming it is due to differences in interactions in the original model. Was my original model incorrect, or this is just a better option? Essentially, the significant treatment_type: trial interaction means there are differences between control and intervention at pre-post measures? I have ~1000 outcomes. I am wondering the best way to adjust p values? Should I only run post hoc tests on those that have a significant treatment:trial interaction – Jess H May 13 '22 at 12:07
  • The trial factor just covers the distinction between baseline and experimental conditions in each period. It shouldn't be surprising if it interacts with things -- it's pretty much supposed to. The first thing we do is compute the differences between experimental and baseline in each arm. Then the trial effect is taken care of and needs no further analysis. – Russ Lenth May 13 '22 at 12:54
  • I think my confusion is related to the output of 'CHG'. It has eight lines, the first two of which includes: ```contrast sequence period treatment_type estimate SE df lower.CL upper.CL treat - base AB Period.1 Control 0.0580 0.0787 64 -0.1583 0.274 treat - base BA Period.1 Control 0.0241 0.0796 64 -0.1947 0.243' ``` However, none of the participants from sequence AB (meaning those who received the intervention first) had the control in period one? Perhaps I misunderstood how to set up the data? – Jess H May 14 '22 at 03:48
  • I did a detailed comparison with a similar simulated example [here](https://www.dropbox.com/s/w0gvn2sa4hq3ca9/xover.html?dl=0) and show that you should get the same results. It is confusing because the reference grid for my model contains cases that don't actually occur; but those cases can still be estimated correctly with this model even though they don't occur. – Russ Lenth May 14 '22 at 04:31
  • I can see now- When I re-run it I also get the same effect change! I cant thank you enough – Jess H May 14 '22 at 06:28