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