I have conducted a randomized 2x2 cross-over trial of 8 participants measuring the effect of a specific diet (intervention) vs normal diet (control) on the number of sleep hours.
The study design includes one week periods: run-in, intervention/control, washout, intervention/control, and then washout. Thus in total 5 weeks. Thus:
Period 1 | Period 2 | |
---|---|---|
Sequence AB | Treatment A | Treatment B |
Sequence BA | Treatment B | Treatment A |
The data looks like this:
subject | sleep_hours | sequence | period | treatment |
---|---|---|---|---|
1 | 4.3 | AB | runin | 0 |
2 | 6.5 | AB | runin | 0 |
3 | 5.2 | AB | runin | 0 |
4 | 4.4 | AB | runin | 0 |
5 | 4.2 | BA | runin | 0 |
6 | 6.5 | BA | runin | 0 |
7 | 5.2 | BA | runin | 0 |
8 | 4.6 | BA | runin | 0 |
1 | 5.2 | AB | 1 | A |
2 | 4.1 | AB | 1 | A |
3 | 6.5 | AB | 1 | A |
4 | 4.4 | AB | 1 | A |
1 | 7.1 | AB | 2 | B |
2 | 8.7 | AB | 2 | B |
3 | 6.5 | AB | 2 | B |
4 | 7.4 | AB | 2 | B |
5 | 7.2 | BA | 1 | B |
6 | 8.3 | BA | 1 | B |
7 | 6.9 | BA | 1 | B |
8 | 7.4 | BA | 1 | B |
5 | 4.8 | BA | 2 | A |
6 | 5.1 | BA | 2 | A |
7 | 4.2 | BA | 2 | A |
8 | 6.6 | BA | 2 | A |
Here is a reproducible code of the above data:
db <- structure(list(subject = c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L,
2L, 3L, 4L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 5L, 6L, 7L, 8L),
sleep_hours = c(4.3, 6.5, 5.2, 4.4, 4.2, 6.5, 5.2, 4.6, 5.2,
4.1, 6.5, 4.4, 7.1, 8.7, 6.5, 7.4, 7.2, 8.3, 6.9, 7.4, 4.8,
5.1, 4.2, 6.6), sequence = structure(c(1L, 1L, 1L, 1L, 2L,
2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L), levels = c("AB", "BA"), class = "factor"),
period = structure(c(3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 1L,
1L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L
), levels = c("1", "2", "runin"), class = "factor"),
treatment = structure(c(1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 2L, 2L, 2L, 2L), levels = c("0", "A", "B"
), class = "factor")), row.names = c(NA, -24L),
class = "data.frame")
I used a Linear mixed-effects model like to see if the intervention (treatment A) affects sleep_hours compared with the control (treatment B):
install.packages(lme4)
install.packages(lmerTest)
model <- lmer(sleep_hours ~ treatment * period + sequence +
(1|subject), data = db)
summary(model)
This gave these results:
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: sleep_hours ~ treatment * period + sequence + (1 | subject)
Data: db
REML criterion at convergence: 57.7
Scaled residuals:
Min 1Q Median 3Q Max
-1.0220 -0.6944 -0.1703 0.3407 1.5199
Random effects:
Groups Name Variance Std.Dev.
subject (Intercept) 0.0000 0.000
Residual 0.9101 0.954
Number of obs: 24, groups: subject, 8
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 5.100000000000009415 0.477005998098789574 17.999999999553889296 10.692 0.00000000316 ***
treatmentA -0.050000000000007615 0.674588351844622736 17.999999999553857322 -0.074 0.94173
treatmentB 2.325000000000000178 0.674588351844622736 17.999999999553857322 3.447 0.00288 **
period2 -0.000000000000009721 0.954011996197578593 17.999999999553843111 0.000 1.00000
sequenceBA 0.024999999999988885 0.674588351844622736 17.999999999553867980 0.037 0.97085
treatmentA:period2 0.100000000000017311 1.652397248444412492 17.999999999616743906 0.061 0.95241
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) trtmnA trtmnB perid2 sqncBA
treatmentA -0.707
treatmentB 0.000 0.000
period2 -0.500 0.354 -0.707
sequenceBA -0.707 0.500 -0.500 0.707
trtmntA:pr2 0.577 -0.612 0.612 -0.866 -0.816
fit warnings:
fixed-effect model matrix is rank deficient so dropping 4 columns / coefficients
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
I also tried to extract the P value and 95% confidence intervals for the effects of treatment A, inspired by the code of @BenBolker here:
model_p <- pvalue(coef(summary(as(model,
"lmerModLmerTest")))[2,5])
model_ci <- paste0(ciformat1(coef(summary(as(model,
"lmerModLmerTest")))[2,1]), " (", ciformat1(confint(model,
method="Wald")[4,1]), " to ", ciformat1(confint(model,
method="Wald")[4,2]), ")")
Which gave these results:
Intervention period | P | |
---|---|---|
Sleep hours (hours) | -0.1 (-1.4 to 1.3) | .942 |
My interpretation: treatment A did not affect sleep differently than treatment B, mean difference -0.1 hours (95% CI -1.4 to 1.3 hours), P = .942.
Is my code and interpretation correct?