Slope Model Introduction

Slope model, or random slope and intercept model can be expressed as follows in the context of a clinical trial with two treatment arms and one placebo arm:

\[Y_{ij} = \alpha_0 + \alpha_1*I_{ARM_i = trt1} + \alpha_2*I_{ARM_i = trt2} + a_i + \beta_0 * T_{ij} + \beta_1 * I_{ARM_i = trt1} * T_{ij}+ \\\beta_2 * I_{ARM = trt2} * T_{ij} + \beta_1 * I_{ARM_i = trt1} * T_{ij} + b_i * T_{ij} + \epsilon_{ij}\] where

  • \(Y_{ij}\) is the response variable of \(i\) th patient at \(j\)th timepoint
  • \(ARM_i\) is the treatment arm assignment of \(i\) th patient and \(I_{\{\}}\) is the indicator function
  • \(\alpha_0, \alpha_0 + \alpha_1, \alpha_0 + \alpha_2\) are fixed-effects intercepts for placebo, treatment arm 1 and treatment arm 2 respectively
  • \(\beta_0, \beta_0+\beta_1, \beta_0 + \beta_2\) are the fixed-effects slopes for placebo, treatment arm 1 and treatment arm 2 respectively
  • \(a_i\) and \(b_i\) represent the random intercept and slope for \(i\)th patient where \((a_i, b_i) \overset{\mathrm{iid}}{\sim} N(\mathbf{0}, \Omega)\)
  • \(\epsilon_{ij}\) is the residual error of \(i\) th patient at \(j\) th timepoint and \(\epsilon_{ij} \overset{\mathrm{iid}}{\sim} N(0, \sigma^2)\)

Test Data Generation

expit <- function(x) {
  exp(x)/(1+exp(x))
}

set.seed(1234)
n <- 450
library(MASS)
corr <- matrix(1, nrow=4, ncol=4) + diag(0.5, nrow=4)
data <- mvrnorm(n, mu=c(0,0,0,0), Sigma=corr)

random_num <- runif(n)
trt <- ifelse(random_num<(1/3), "pbo", ifelse(random_num>(2/3), "trt2", "trt1"))
y0 <- exp(data[,1])
y1 <- exp(data[,2])
y2 <- exp(data[,3])
y3 <- exp(data[,4])

#add in effect of treatment
y1 <- y1+(trt=="trt2")*0.5
y2 <- y2+(trt=="trt2")*1
y3 <- y3+(trt=="trt2")*1.5

wideData <- data.frame(id=1:n, trt=trt, y0=y0, y1=y1, y2=y2, y3=y3)

longData <- reshape(wideData, varying=c("y0", "y1", "y2", "y3"),
                                   direction="long", sep="", idvar="id")
longData <- longData[order(longData$trt,longData$id,longData$time),]
longData$time <- longData$time + runif(length(longData$time), -0.1, 0.1)
longData <- longData %>% 
   mutate(log_y = log(y))

#export data to CSV
write.csv(longData, file="slope_model_dat.csv", row.names = FALSE, na=".")

The specs of the test data set are

  • id is the patient ID

  • trt is the treatment group

  • time is the time of the assessment in years

  • y is the response variable

  • log_y is the response variable after the transformation

Slope Model Inference

Firstly, we can use lmer function to fit a slope model. Because the data is right skewed, a log transformation was conducted before we run the model. The results should be interpreted in the log scale. No back transformation is applied.

result1 <- longData%>%
  lmerTest::lmer(log_y ~  time + trt*time + (time|id), REML = T, data = .)

Secondly, we are interested in the annualized rate of change of each arm and the treatment effect between treatment and placebo arm. Hereby, we are reporting the least-square means for the progression rate (i.e., \(\hat{\beta_0}, \hat{\beta_0} + \hat{\beta_1}, \hat{\beta_0} + \hat{\beta_2}\)) with its 95% CI of each arm and the least-square means of the difference in progression rate (i.e., \(\hat{\beta_1}, \hat{\beta_2}\)) with its 95% CI between the treatment groups and placebo.

result2 <- emmeans::emtrends(result1, trt.vs.ctrl~trt, var="time", adjust = "none")
result2
## $emtrends
##  trt  time.trend     SE  df lower.CL upper.CL
##  pbo     -0.0495 0.0257 447  -0.1000 0.000901
##  trt1     0.0208 0.0246 448  -0.0275 0.069086
##  trt2     0.3654 0.0249 444   0.3165 0.414428
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast   estimate     SE  df t.ratio p.value
##  trt1 - pbo   0.0703 0.0355 447   1.979  0.0484
##  trt2 - pbo   0.4150 0.0358 445  11.600  <.0001
## 
## Degrees-of-freedom method: kenward-roger

In addition, the relative difference in progression rate (i.e., \(\hat{\beta_1}/\hat{\beta_0}, \hat{\beta_2}/\hat{\beta_0}\)) could be calculated as follows:

pbo_lsmean <- summary(result2)[["emtrends"]]["time.trend"][[1]][1]
trt1_lsmean <- summary(result2)[["emtrends"]]["time.trend"][[1]][2]
trt2_lsmean <- summary(result2)[["emtrends"]]["time.trend"][[1]][3]

(trt1_lsmean - pbo_lsmean) / pbo_lsmean
## [1] -1.419747
(trt2_lsmean - pbo_lsmean) / pbo_lsmean
## [1] -8.378397

A equivalent SAS code is also provided below with identical results.

# SAS code 
# proc MIXED data = WORK.IMPORT order = internal covtest ;
#  CLASS id trt (REF="pbo") ;
# MODEL log_y = trt time trt*time /solution CL ddfm = KR;
# RANDOM intercept time/ type=un subject=id;
# ESTIMATE "Annual Rate of Change in Placebo" time 1 trt*time 0 0 1 / e cl;
# ESTIMATE "Annual Rate of Change in Treatment 1"    time 1 trt*time 1 0 0 / e cl;
# estimate "Annual Rate of Change in Treatment 2"    time 1 trt*time 0 1 0 / e cl;
# ESTIMATE "Difference in Annual Rates of Change in Treatment 1" trt*time 1 0 -1 / alpha=0.05 e CL;
# ESTIMATE "Difference in Annual Rates of Change in Treatment 2" trt*time 0 1 -1 / alpha=0.05 e CL;
# RUN;