Constrainted longitudinal data analysis (cLDA) in R. Example with the `fev_data` from the `mmrm` package.

cLDA.R

library(tidyverse)
library(nlme)
library(rms)
library(lmerTest)
library(emmeans)

data("fev_data", package = "mmrm")
data <- fev_data %>%
  as_tibble() %>%
  mutate(ARMCD = as.character(ARMCD)) %>%
  select(-VISITN, -VISITN2) |> 
  filter(!is.na(FEV1))


# Approach 1: custom model matrix -----------------------------------------

# Build alternative model matrix without VIS1
mod_mat <- model.matrix(~ AVISIT + AVISIT:ARMCD, data)
cols <- colnames(mod_mat)
alt_cols <- cols[!str_detect(cols, "VIS1|Intercept")]
alt_mod_mat <- mod_mat[,alt_cols]

# Fit the linear mixed model using this matrix
fit <- Gls(
  FEV1 ~ alt_mod_mat,
  data = data,
  weights = varIdent(form = ~ 1 | AVISIT), 
  correlation=corSymm (form = ~ 1 | USUBJID), 
)

# Predict values with 95% CI
predictions <- cbind(data, predict(fit, data, conf.int=0.95)) |> 
  distinct(AVISIT, ARMCD, .keep_all = TRUE) |> 
  mutate(ARMCD2 = ifelse(AVISIT == "VIS1", "All", ARMCD))

# Plot estimated marginal means over visits
ggplot(predictions, aes(x = AVISIT, y = linear.predictors)) +
  geom_line(aes(group = ARMCD)) +
  geom_pointrange(aes(ymin = lower, ymax = upper, color = ARMCD2)) +
  ylab("Estimated mean with corresponding 95% CI") +
  see::theme_modern()

# Approach 2: using dummy variables ---------------------------------------

# Build a dummy variable for the interaction between each visit (except the first) and the treatment 
dummy_visits <- levels(data$AVISIT)[-1]
for (v in dummy_visits) {
  data[[v]] <- as.numeric(data$ARMCD == "TRT" & data$AVISIT == v)
}

# Build the formula for the model
form <- paste0("FEV1 ~ AVISIT + ", paste(dummy_visits, collapse = " + "), " + (1 | USUBJID)") |> 
  as.formula()

# Fit the model using standard LMM functions
fit2 <- lmer(
  formula = form,
  data = data
)

# Extract the marginal means using emmeans
emm <- emmeans(
  object = fit2,
  specs = c("AVISIT", dummy_visits), 
) |> 
  confint() |> 
  broom::tidy(conf.int = TRUE)

# Build a matrix of only the contrasts of interest
contr <- data.frame(AVISIT = c("VIS1", rep(dummy_visits, each = 2)))
for (i in seq_len(length(dummy_visits))) {
  name <- dummy_visits[i]
  contr[[name]] <- c(rep(0, i*2), 1, rep(0, (3-i)*2))
}

# Select only these contrast in the estimate MM
emm_contr <- emm |> 
  right_join(contr, by = join_by(AVISIT, VIS2, VIS3, VIS4)) |> 
  rowwise() |> 
  mutate(
    ARMCD = ifelse(sum(c_across(dummy_visits)) == 1, "TRT", "PBO"),
    ARMCD2 = ifelse(AVISIT == "VIS1", "All", ARMCD)
  ) |> 
  ungroup()

# Duplicate first row to have lines between points
emm_contr_full <- rbind(
  emm_contr |> slice_head() |> mutate(ARMCD = "TRT"),
  emm_contr
)

# Plot estimated marginal means over visits
ggplot(emm_contr_full, aes(x = AVISIT, y = estimate)) +
  geom_line(aes(group = ARMCD)) +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high, color = ARMCD2)) +
  ylab("Estimated mean with corresponding 95% CI") +
  see::theme_modern()