Search code examples
lme4emmeanslsmeans

Contrast emmeans: post-hoc t-test as the average differences of the differences between baseline and treatment periods

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!

Solution

  • 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.