Mixed model for clinical trials
Enwu Liu
2024-04-19
1. Repeated measurement and clinical tiral
Repeated measures are quite common in clinical trials. A repeated measures study involves measuring outcomes at different times; it does not necessarily mean that covariates or predictor variables were also measured repeatedly. For the design of repeated measures studies, general linear mixed models, generalized linear mixed models, or GEE models are commonly used to analyze the data.
Mixed linear regression models or GEE regression models are true multivariate regression models because they analyze more than one outcome variable simultaneously. This should not be confused with multiple or multivariable regression, where only one outcome and several predictor variables are analyzed.
For example, suppose we conduct a three-arm clinical trial among 15 patients, with treatments 1 and 2 involving drug 1 and drug 2, respectively, and treatment 3 involving a placebo. We randomize 5 patients to each of the three arms and measure the primary outcome variable at baseline (time 0, before taking any drug/placebo) and follow up at 30 minutes, 60 minutes, and 90 minutes after taking the drug/placebo.
The wide form of the first 9 patients’ data looks like:
If we use linear mixed model to analyse the data, we should make the data as long form which look like
In this format, each row represents a single measurement for a patient at a specific time point, under a specific treatment arm, with the corresponding outcome value
2 Use R and SAS to analyse the data
I will use a the dataset phlebitis.csv from the PennState STATA 510 to show how to analyse a three-arm clinical trial with repeated measures.
An experiment is done to examine ways to detect phlebitis during the intravenous administration of a particular drug. Phlebitis is an inflammation of a blood vein. Three intravenous treatments were administered. 15 test animals were randomly divided into three groups of n = 5. Each group is given a different treatment. The treatments are:
Treatment 1: The drug in a solution designed to carry the drug
Treatment 2:* The carrier solution only (no drug)
Treatment 3: Saline solution
The treatments were administered to one ear of the test animal. The y variable = difference in temperature between the treated ear and the untreated ear. This is measured at times 0, 30 minutes, 60 minutes, and 90 minutes after the treatment is administered. It is believed that increased temperature in the treated ear may be an early sign of phlebitis.
We use the following R code to show what the data looks like
library(nlme)
phlebitisdata<- read.csv("https://raw.githubusercontent.com/enwuliu/datasets/refs/heads/main/phlebitis.csv", header=T, sep=",")
interaction.plot (phlebitisdata$Time, factor(phlebitisdata$Treatment), phlebitisdata$Y, lty=c(1:3),lwd=2,ylab="mean of Y", xlab="time", trace.label="Treatment")
The STATA 510 uses Generalized Least Squares (GLS) to analyze the data. Since the dataset is balanced (i.e., with equal sample sizes for each group and no missing data), both the GLS method and the mixed method will yield the same results.
For this study design, the outcome variable is a continuous variable. Many people believe that the baseline value of the outcome variable should be treated as a covariate variable in the mixed model. However, in STATA 510, no adjustment was made for the baseline value of the outcome variable. For the linear mixed model we will adjust the baseline outcome values in the mixed model, therefore, the covariates are “Treatment,” “Time,”, “Treatment*Time” interaction and baseline outcome value. We will perform a random-intercept model that allows the intercept to vary between each animal (subject).
We use the following R code to merge outcome baseline values to the dataset.
library(tidyverse)
baseline<-phlebitisdata |> filter(Time==0) |> select("Animal","Y") |>rename(base_Y=Y)
phlebitisdata2<-phlebitisdata |> left_join(baseline,by="Animal") |> filter(Time!=0)
head(phlebitisdata2,10)## Animal Treatment Time Y base_Y
## 1 1 1 30 -0.2 -0.3
## 2 1 1 60 1.2 -0.3
## 3 1 1 90 3.1 -0.3
## 4 2 1 30 2.2 -0.5
## 5 2 1 60 3.3 -0.5
## 6 2 1 90 3.7 -0.5
## 7 3 1 30 2.4 -1.1
## 8 3 1 60 2.2 -1.1
## 9 3 1 90 2.7 -1.1
## 10 4 1 30 1.7 1.0For a repeated measures trial to be analysed using a random intercept mixed model, we can specify the linear mixed model as follow:
\[y_{ijk}=\mu + \alpha_j+d_{ij}+\tau_k +(\alpha\tau)_{jk}+\beta*\text{baseline} +e_{ijk}\tag{1}\] where \(y_{ijk}\) is the \(i^{th}\) subject’s response for treatment \(j\) at time \(k\), for the above experiment, \(i=1,2,...15\), since there are 15 animals in total, \(j=1,2,3\) since there are three treatment groups, and there are four time points \(k=30,60,90\). note, we treat baseline time as a covariate.
\(\mu\) is the intercept stands for an overall mean.
\(\alpha_j\) is a fixed effect of treatment \(j\)
\(d_{ij}\) is a random effect of subject \(i\) in treatment group \(j\), note we can put \((\mu+d_{ij})\) together, then we think the intercept has some randomness.
\(\tau_k\) is a fixed effect of time \(k\).
\((\alpha \tau)_{jk}\) is a fixed interaction effect of treatment \(j\) with time \(k\).
In term of matrix form of the general linear mixed model, the vector \(\beta\) contains the fixed effect parameters \(\mu, \alpha_j,\tau_k\) and \((\alpha \tau)_{jk}\). The random vector \(U\) contains the between-subject residual errors \(d_{jk}\) and \(e\) contains the within-subject residual \(e_{ijk}\)
For the above study design, there were 3 treatment groups and four time points. We treat both treatment and time point as categorical variables, so there will be two dummy variables for treatment and three dummy variables for time point. In R, we use the factor() function to create these dummy variables. Therefore, our linear mixed model should produce 13 coefficients: one for the intercept, 2 for the two dummy treatment groups, 2 dummy variables for the 3 follow up time points (30,60,90), 4 (2 × 4) coefficients for the treatment and time interaction, and one extra coefficient for the baseline outcome score.
We use the following R code to perform the linear mixed model analysis:
library(lme4)
library(lmerTest)
library(lsmeans)
m1 <- lmer(Y ~ factor(Treatment)*factor(Time) + base_Y +(1|Animal),data=phlebitisdata2)
summary(m1)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Y ~ factor(Treatment) * factor(Time) + base_Y + (1 | Animal)
## Data: phlebitisdata2
##
## REML criterion at convergence: 102.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.05994 -0.49673 0.04788 0.64105 1.60227
##
## Random effects:
## Groups Name Variance Std.Dev.
## Animal (Intercept) 0.2179 0.4668
## Residual 0.5159 0.7183
## Number of obs: 45, groups: Animal, 15
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.37139 0.38758 28.20696 3.538 0.00142
## factor(Treatment)2 -1.91219 0.54815 28.20520 -3.488 0.00161
## factor(Treatment)3 -1.15578 0.55240 27.86610 -2.092 0.04565
## factor(Time)60 0.50000 0.45426 24.00000 1.101 0.28195
## factor(Time)90 1.20000 0.45426 24.00000 2.642 0.01429
## base_Y -0.03586 0.24513 11.00000 -0.146 0.88634
## factor(Treatment)2:factor(Time)60 -0.04000 0.64243 24.00000 -0.062 0.95087
## factor(Treatment)3:factor(Time)60 -1.16000 0.64243 24.00000 -1.806 0.08353
## factor(Treatment)2:factor(Time)90 -0.68000 0.64243 24.00000 -1.058 0.30037
## factor(Treatment)3:factor(Time)90 -1.22000 0.64243 24.00000 -1.899 0.06964
##
## (Intercept) **
## factor(Treatment)2 **
## factor(Treatment)3 *
## factor(Time)60
## factor(Time)90 *
## base_Y
## factor(Treatment)2:factor(Time)60
## factor(Treatment)3:factor(Time)60 .
## factor(Treatment)2:factor(Time)90
## factor(Treatment)3:factor(Time)90 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) fc(T)2 fc(T)3 f(T)60 f(T)90 base_Y f(T)2:(T)6 f(T)3:(T)6
## fctr(Trtm)2 -0.668
## fctr(Trtm)3 -0.656 0.514
## factr(Tm)60 -0.586 0.414 0.411
## factr(Tm)90 -0.586 0.414 0.411 0.500
## base_Y 0.152 0.152 0.195 0.000 0.000
## f(T)2:(T)60 0.414 -0.586 -0.291 -0.707 -0.354 0.000
## f(T)3:(T)60 0.414 -0.293 -0.581 -0.707 -0.354 0.000 0.500
## f(T)2:(T)90 0.414 -0.586 -0.291 -0.354 -0.707 0.000 0.500 0.250
## f(T)3:(T)90 0.414 -0.293 -0.581 -0.354 -0.707 0.000 0.250 0.500
## f(T)2:(T)9
## fctr(Trtm)2
## fctr(Trtm)3
## factr(Tm)60
## factr(Tm)90
## base_Y
## f(T)2:(T)60
## f(T)3:(T)60
## f(T)2:(T)90
## f(T)3:(T)90 0.500lsmeans(m1,pairwise~Treatment|Time, ddf="Satterthwaite",adjust = "none")## $lsmeans
## Time = 30:
## Treatment lsmean SE df lower.CL upper.CL
## 1 1.38932 0.388 28.1 0.594 2.185
## 2 -0.52287 0.384 28.7 -1.308 0.262
## 3 0.23355 0.386 28.4 -0.556 1.023
##
## Time = 60:
## Treatment lsmean SE df lower.CL upper.CL
## 1 1.88932 0.388 28.1 1.094 2.685
## 2 -0.06287 0.384 28.7 -0.848 0.722
## 3 -0.42645 0.386 28.4 -1.216 0.363
##
## Time = 90:
## Treatment lsmean SE df lower.CL upper.CL
## 1 2.58932 0.388 28.1 1.794 3.385
## 2 -0.00287 0.384 28.7 -0.788 0.782
## 3 0.21355 0.386 28.4 -0.576 1.003
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## Time = 30:
## contrast estimate SE df t.ratio p.value
## Treatment1 - Treatment2 1.912 0.548 28.2 3.488 0.0016
## Treatment1 - Treatment3 1.156 0.552 27.9 2.092 0.0456
## Treatment2 - Treatment3 -0.756 0.542 28.7 -1.395 0.1738
##
## Time = 60:
## contrast estimate SE df t.ratio p.value
## Treatment1 - Treatment2 1.952 0.548 28.2 3.561 0.0013
## Treatment1 - Treatment3 2.316 0.552 27.9 4.192 0.0003
## Treatment2 - Treatment3 0.364 0.542 28.7 0.670 0.5080
##
## Time = 90:
## contrast estimate SE df t.ratio p.value
## Treatment1 - Treatment2 2.592 0.548 28.2 4.729 <0.0001
## Treatment1 - Treatment3 2.376 0.552 27.9 4.301 0.0002
## Treatment2 - Treatment3 -0.216 0.542 28.7 -0.399 0.6928
##
## Degrees-of-freedom method: kenward-rogerThe following SAS code will produce the same results as the R code above
filename phleb url 'https://raw.githubusercontent.com/enwuliu/datasets/refs/heads/main/phlebitis.csv';
proc import datafile=phleb dbms=csv out=phleb2 replace;
run;
data base(keep=animal y rename=(y=base_score));
set phleb2;
where time=0;
run;
proc sort data=phleb2;
by animal;
run;
proc sort data=base;
by animal;
run;
data pleb_with_base;
merge phleb2 base;
by animal;
run;
data pleb_with_base2;
set pleb_with_base;
where time^=0;
run;
proc mixed data=pleb_with_base2 method=REML ;
class animal treatment(ref='1') time(ref='30');
model Y = treatment time treatment*time base_score/solution DDFM=SATTERTHWAITE;
random intercept / subject=animal; /*we use default variance structure which is Variance Components*/
lsmeans treatment*time/diff /*ADJUST=TUKEY, not to use*/;
run;We obtained the following SAS results, which are the same as the R results above.
Based on above results we can write the linear mixed effect model as follow:
\[ \begin{aligned} Y_{ik} &=\beta_0+\beta_1*\text{Treatment}_{2i}+\beta_2*\text{Treatment}_{3i}+\beta_3*\text{Time}_{60}+\beta_4*\text{Time}_{90}\\&+ \beta_5*(\text{Treatment}_{2i}*\text{Time}_{60})+\beta_6*(\text{Treatment}_3*\text{Time}_{60})+\beta_7*(\text{Treatment}_2*\text{Time}_{90})\\ &+\beta_8*(\text{Treatment}_3*\text{Time}_{90})+\beta_{9}*\text{BaselineY}+b_i+\epsilon_{ik} \end{aligned}\] Consequently:
\(\beta_0\) is the expected outcome at 30 minutes for Treatment 1 when base_Y = 0.
\(\beta_1\) is the adjusted difference between Treatments 2 and 1 at 30 minutes.
\(\beta_2\) is the adjusted difference between Treatments 3 and 1 at 30 minutes.
\(\beta_3\) is the change from 30 to 60 minutes in Treatment 1.
\(\beta_4\) is the change from 30 to 90 minutes in Treatment 1.
\(\beta_5\)–\(\beta_8\) are difference-in-differences relative to 30 minutes.
\(\beta_9\) is the expected difference in a post-baseline outcome associated with a one-unit difference in baseline outcome, holding treatment and time constant.
where \[b_i\sim N(0,\sigma_b^2), \text{ }\epsilon_{ik}\sim N(0,\sigma^2)\]
Notes for the mixed model
For the lsmeans, we set the degree of freedom (df) be estimate by “Satterthwaite” method, there are other methods to estimate the df in the mixed model, such as ““Kenward-Roger”. In SAS we use DDF=“SATTERTHWAITE” option.
Since there are three arms, R will do multiple comparisons adjustment using Tukey method, we set adjust =‘none’, not to adjust multiple comparison. The default setting for SAS is not to adjust the multiple comparison.