离散时间模型中基于R公式实现cLDA约束的方法问询

我从事临床试验重复测量数据的纵向模型研究，患者被随机分配到不同治疗组，在多个预设时间点接受测量。示例数据集采用mmrm包中的FEV数据集，其中ARMCD代表治疗组，AVISIT代表离散时间点：

library(tidyverse)
data("fev_data", package = "mmrm")
data <- fev_data %>%
  as_tibble() %>%
  mutate(ARMCD = as.character(ARMCD)) %>%
  select(-VISITN, -VISITN2)
data
#> # A tibble: 800 × 8
#>    USUBJID AVISIT ARMCD RACE                      SEX    FEV1_BL  FEV1 WEIGHT
#>    <fct>   <fct>  <chr> <fct>                     <fct>    <dbl> <dbl>  <dbl>
#>  1 PT1     VIS1   TRT   Black or African American Female    25.3  NA    0.677
#>  2 PT1     VIS2   TRT   Black or African American Female    25.3  40.0  0.801
#>  3 PT1     VIS3   TRT   Black or African American Female    25.3  NA    0.709
#>  4 PT1     VIS4   TRT   Black or African American Female    25.3  20.5  0.809
#>  5 PT2     VIS1   PBO   Asian                     Male      45.0  NA    0.465
#>  6 PT2     VIS2   PBO   Asian                     Male      45.0  31.5  0.233
#>  7 PT2     VIS3   PBO   Asian                     Male      45.0  36.9  0.360
#>  8 PT2     VIS4   PBO   Asian                     Male      45.0  48.8  0.507
#>  9 PT3     VIS1   PBO   Black or African American Female    43.5  NA    0.682
#> 10 PT3     VIS2   PBO   Black or African American Female    43.5  36.0  0.892
#> # ℹ 790 more rows
#> # ℹ Use `print(n = ...)` to see more rows

临床研究中常采用约束纵向数据分析（cLDA），要求基线时合并所有治疗组。如果时间是连续型且基线为AVISIT == 0，可以通过以下公式轻松实现cLDA：

formula <- FEV1 ~ AVISIT + AVISIT:ARMCD

但这种方法不适用于离散时间，因为生成的模型矩阵会包含AVISITVIS1:ARMCDTRT项，导致基线时PBO和TRT组被区分开：

colnames(model.matrix(formula, data = data))
#> [1] "(Intercept)"         "AVISITVIS2"          "AVISITVIS3"          "AVISITVIS4"         
#> [5] "AVISITVIS1:ARMCDTRT" "AVISITVIS2:ARMCDTRT" "AVISITVIS3:ARMCDTRT" "AVISITVIS4:ARMCDTRT"

我见过有人用手动修改数据的方式实现cLDA：

count(data, ARMCD, AVISIT)
#> # A tibble: 8 × 3
#>   ARMCD AVISIT     n
#>   <chr> <fct>  <int>
#> 1 PBO   VIS1     105
#> 2 PBO   VIS2     105
#> 3 PBO   VIS3     105
#> 4 PBO   VIS4     105
#> 5 TRT   VIS1      95
#> 6 TRT   VIS2      95
#> 7 TRT   VIS3      95
#> 8 TRT   VIS4      95
clda <- mutate(data, ARMCD = ifelse(AVISIT == "VIS1", "PBO", ARMCD))
count(clda, ARMCD, AVISIT)
#> # A tibble: 7 × 3
#>   ARMCD AVISIT     n
#>   <chr> <fct>  <int>
#> 1 PBO   VIS1     200
#> 2 PBO   VIS2     105
#> 3 PBO   VIS3     105
#> 4 PBO   VIS4     105
#> 5 TRT   VIS2      95
#> 6 TRT   VIS3      95
#> 7 TRT   VIS4      95

但这种方式并没有改变模型的项：

colnames(model.matrix(formula, data = clda))
#> [1] "(Intercept)"         "AVISITVIS2"          "AVISITVIS3"          "AVISITVIS4"         
#> [5] "AVISITVIS1:ARMCDTRT" "AVISITVIS2:ARMCDTRT" "AVISITVIS3:ARMCDTRT" "AVISITVIS4:ARMCDTRT"

更糟的是，模型矩阵不再是满秩：

as.integer(Matrix::rankMatrix(model.matrix(formula, data = clda)))
#> [1] 7

有没有无需重编码数据或手动生成模型矩阵，仅通过公式和/或对比就能实现cLDA的方法？手动生成自定义模型矩阵后调用lm()或brms::brm()的方式太底层，还会丢失brms自带的便捷后处理工具。

方法1：通过自定义对比矩阵实现cLDA #

可以为AVISIT:ARMCD交互项设置自定义对比，强制基线时（VIS1）的治疗组差异为0。

# 生成交互项的所有水平
int_levels <- expand.grid(AVISIT = levels(data$AVISIT), ARMCD = unique(data$ARMCD))
int_levels <- paste(int_levels$AVISIT, int_levels$ARMCD, sep = ":")

# 构造对比矩阵：仅保留VIS2-VIS4的组间差异，VIS1差异约束为0
contrasts_mat <- matrix(0, nrow = length(int_levels), ncol = 3)
colnames(contrasts_mat) <- c("VIS2:TRT-PBO", "VIS3:TRT-PBO", "VIS4:TRT-PBO")

# 为各时间点组间差异赋值
contrasts_mat[int_levels == "VIS2:TRT", 1] <- 1
contrasts_mat[int_levels == "VIS2:PBO", 1] <- -1

contrasts_mat[int_levels == "VIS3:TRT", 2] <- 1
contrasts_mat[int_levels == "VIS3:PBO", 2] <- -1

contrasts_mat[int_levels == "VIS4:TRT", 3] <- 1
contrasts_mat[int_levels == "VIS4:PBO", 3] <- -1

# 将对比矩阵应用到交互项
contrasts(interaction(data$AVISIT, data$ARMCD)) <- contrasts_mat

然后使用包含主效应和交互项的公式拟合模型：

formula_clda <- FEV1 ~ AVISIT + AVISIT:ARMCD
model <- lm(formula_clda, data = data)
summary(model)

这种方式不修改原始数据，模型矩阵保持满秩，同时实现了基线合并的cLDA约束。

方法2：调整公式结构，显式约束基线 #

重新构造公式，将基线作为共享截距，仅对后续时间点设置组间差异：

# 生成时间点指示变量和治疗组指示变量
data <- data %>%
  mutate(
    is_vis2 = as.integer(AVISIT == "VIS2"),
    is_vis3 = as.integer(AVISIT == "VIS3"),
    is_vis4 = as.integer(AVISIT == "VIS4"),
    trt = as.integer(ARMCD == "TRT")
  )

# cLDA公式：截距代表所有组的基线均值，后续时间点分别设置组间差异
formula_clda <- FEV1 ~ 1 + is_vis2 + is_vis3 + is_vis4 + is_vis2:trt + is_vis3:trt + is_vis4:trt
model <- lm(formula_clda, data = data)
summary(model)

公式含义：

• 截距：所有治疗组在VIS1（基线）的均值
• is_vis2/is_vis3/is_vis4：PBO组对应时间点相对于基线的变化
• is_vis2:trt等：TRT组相对于PBO组在对应时间点的额外变化

该方式兼容brms，可直接使用其后续处理工具：

library(brms)
brm_model <- brm(formula_clda, data = data, family = gaussian())
summary(brm_model)

方法3：用emmeans施加事后约束 #

如果已拟合包含完整交互项的模型，可通过emmeans设置自定义对比，强制基线组间差异为0：

library(emmeans)
# 拟合完整交互模型
full_model <- lm(FEV1 ~ AVISIT * ARMCD, data = data)

# 定义对比：基线差异约束为0，保留后续时间点组间差异
clda_contrasts <- list(
  vis2_diff = c(0, 1, 0, 0, 0, -1, 0, 0),
  vis3_diff = c(0, 0, 1, 0, 0, 0, -1, 0),
  vis4_diff = c(0, 0, 0, 1, 0, 0, 0, -1),
  vis1_constraint = c(1, 0, 0, 0, -1, 0, 0, 0)
)

# 应用对比并输出结果
emmeans(full_model, ~ AVISIT * ARMCD) %>%
  contrast(clda_contrasts)

此方法适合已拟合模型后需要施加cLDA约束的场景，无需修改原始数据。

内容的提问来源于stack exchange，提问作者landau