# Repeated measure clinical trial Linear mixed model

I am working on a clinical trial testing an innovative rehabilitation therapy on patients and I would like some suggestions on how to analyse the data.

The study design is: 2-groups: conventional (n=17) vs innovative (n=15) treatment; 4-time points (pre-therapy, T0; halfway through the therapy period, T1; end of therapy, T2; 2 months follow-up, T3). As output, we record a continuous variable: time (in seconds) to walk from point A to point B. We record multiple values for each subject at each time point (T0,T1,T2,T3).

I have seen articles where the authors used ANOVA to evaluate the change in performance within each group. However, I would also like to evaluate the difference between groups, preferably at each time point as well. I have noticed other studies using Linear Mixed Models.

I thought of setting as Fixed effect: Group, Time-point; Random effect: intercept per subject, Group per subject.

However, I am not an expert of this technique and I do not really know how to apply it, or how to write it in R/Statsmodels. Can you please help me understand how to set it up and how to write it on a coding platform?

• How is this part described in the statistical methodology part in the trial protocol and/or the statistical analysis plan? Jul 12, 2020 at 12:04
• Hi Michael. It is a very good point. This study is carried by several institutes, and the one where I work is not in charge of the statistical analysis. Unfortunately, I have no access to the info you asked. However, I want to learn how to do this analysis. Jul 12, 2020 at 12:10
• "I thought of setting as Fixed effect: Group, Time-point; Random effect: intercept per subject, treatment per subject.". You mention Group as a variable and also treatment, but aren't they the same ? Jul 12, 2020 at 12:17
• Yes, sorry. I change the question. Jul 12, 2020 at 13:05

outcome ~ group * time-point + (group | subject)


This will estimate fixed effects for group, time-point and the interaction between them. If you have sufficient statistical power, this will enable you to answer the research questions. The main effect for group will estimate the difference in outcome between the two groups at time-point 0, the main effect for time-point (assuming this is a factor and not continuous) will have 3 estimates, each being the estimated difference in the outcome between each time-point estimate and time-point 0, in the conventional group (that is, if convetional is the reference level for group variable). The interactions will estimate the difference in the outcome between the two groups for each time-point relative to time-poiht 0.

The model also fits random slopes for group which will allow the "effect" of group to be different for each subject (ie as an offset to the main effect for group)

Often in longitudinal studies you would also want the effect of time (and the interactions) to vary by subject, that is, to fit random slopes, but in this case with relatively few subjects and a time variable with 4 levels it is possible that such a model will not be supported by the data. An alternative would be to code time as continuous but then you would need to know the actual time points of each measurement and possibly allow for non linearity.

• thanks for the clear explanation! You said that the main effect for group will estimate the difference in outcome between the two groups at time-point 0. What if I also what to know the difference between the two groups at other time points? Should I perform any PostHoc test to evaluate the difference among groups at other time points? Jul 12, 2020 at 19:22
• The interaction terms will inform you about that but they will be contrasts vs the reference level (which I have assumed will me 0. If you want contrasts vs the other levels then you can change the reference level. This introduces problems with multiple testing but I believe there are some additional packages that will handle that (lsmeans is one I am aware of though there are probably others). I would recommend some research around these methods first. Jul 12, 2020 at 20:40

The model suggested by @RobertLong is a great start (+1). It has one quirk though that may be unintended: there is a treatment effect at each time point, including at time T0 (pre-therapy). \begin{aligned} \operatorname{E}\left\{Y | \text{time} = T_0, \text{group} = C\right\} &= \beta_0 \\ \operatorname{E}\left\{Y | \text{time} = T_0, \text{group} = T\right\} &= \beta_0 + \beta_1 \end{aligned}

However, the treatment can't have an effect on the outcome before the treatment is administered and if the assignment of subjects to a treatment group is randomized, the expected difference between the two groups at baseline is 0. \begin{aligned} \operatorname{E}\left\{Y | \text{time} = T_0, \text{group} = C\right\} &= \beta_0 \\ \operatorname{E}\left\{Y | \text{time} = T_0, \text{group} = T\right\} &= \beta_0 \qquad \end{aligned} For the other time points (T1, T2, T3) the model remains unchanged.

Here is how to specify this model in the R formula syntax. It's a bit wordy but gets the job done.

  outcome ~
I(time == "T1") + I(time == "T2") + I(time == "T3")
+ I(group == "T" & time == "T1")
+ I(group == "T" & time == "T2")
+ I(group == "T" & time == "T3")
+ (group | subject)


The terms I(group == "T" & time == "Ti") are the contrasts between treatment and control for each time point Ti = {T1, T2, T3}. The benefit of this formulation is that we estimate the fixed effects more precisely, ie. with smaller std. error.

I ran a simulation to illustrate (full R code below). First the fixed effects from Model 1, which includes a treatment effect at T0:

#>               Estimate Std. Error t value
#> (Intercept)    0.09103    0.05012   1.816
#> groupT        -0.11822    0.06687  -1.768
#> timeT1        -0.06673    0.04887  -1.366
#> timeT2         0.05765    0.04887   1.180
#> timeT3         0.17911    0.04887   3.665
#> groupT:timeT1  0.21681    0.07137   3.038
#> groupT:timeT2  0.20880    0.07137   2.925
#> groupT:timeT3  0.39799    0.07137   5.576


Of primary interest are the contrasts between treatment and control at time points T1 to T3:

#> time = T1:
#>  contrast estimate     SE   df t.ratio p.value
#>  T - C      0.0986 0.0669 85.2   1.474  0.1441
#> time = T2:
#>  contrast estimate     SE   df t.ratio p.value
#>  T - C      0.0906 0.0669 85.2   1.355  0.1791
#> time = T3:
#>  contrast estimate     SE   df t.ratio p.value
#>  T - C      0.2798 0.0669 85.2   4.184  0.0001


And next the fixed effects from Model 2, which doesn't have a treatment effect at T0. The bottom three rows correspond to the contrasts T - C at time points T1, T2 and T3; the std. errors are smaller.

#>                                    Estimate Std. Error t value
#> (Intercept)                         0.02416    0.03357   0.720
#> I(time == "T1")TRUE                -0.03590    0.04561  -0.787
#> I(time == "T2")TRUE                 0.08849    0.04561   1.940
#> I(time == "T3")TRUE                 0.20995    0.04561   4.603
#> I(group == "T" & time == "T1")TRUE  0.15102    0.06069   2.488
#> I(group == "T" & time == "T2")TRUE  0.14301    0.06069   2.356
#> I(group == "T" & time == "T3")TRUE  0.33221    0.06069   5.474


Finally one more point. The comment

This study is carried by several institutes (...)

suggests this is a multi-center study. If this is the case, the model should account for "center effects" as well; whether as fixed or random depends on the details.

The simulation R code:

library("lme4")

set.seed(123)

# T0: T - C = 0.0 - 0.0 = 0.0
# T1: T - C = 0.2 - 0.0 = 0.2
# T2: T - C = 0.3 - 0.1 = 0.2
# T3: T - C = 0.6 - 0.2 = 0.3

sim_fixef <- data.frame(
group = rep(c("T", "C"), each = 4),
time = rep(c("T0", "T1", "T2", "T3"), times = 2),
fixef = c(0.0, 0.2, 0.3, 0.6, 0.0, 0, 0.1, 0.2)
)

k <- 4 # time points
nC <- 17 # subjects in C (control group)
nT <- 15 # subjects in T (treatment group)
n <- nC + nT # subjects in study
r <- 4 # measurements taken per subject per time point
tau <- 0.1
sigma <- 0.3

sim_ranef <- data.frame(
group = rep(c("C", "T"), c(nC, nT)),
subject = seq(nC + nT),
ranef = rnorm(nC + nT, sd = tau)
)

d <- merge(sim_fixef, sim_ranef)
sim_data <- d[rep(seq(n * k), r), ] # replicated each row r times
sim_data\$outcome <- with(
sim_data,
rnorm(n * k * r, mean = fixef + ranef, sd = sigma)
)

fit1 <- lmer(
outcome ~ group * time + (group | subject),
data = sim_data
)

fit2 <- lmer(
outcome ~
I(time == "T1") + I(time == "T2") + I(time == "T3")
+ I(group == "T" & time == "T1")
+ I(group == "T" & time == "T2")
+ I(group == "T" & time == "T3")
+ (group | subject),
data = sim_data
)

summary(fit1)
summary(fit2)

library("emmeans")
pairs(emmeans(fit1, ~group, by = "time"), reverse = TRUE)