我有一个随机对照试验的纵向数据集,其中有两个组 (
arm
),以及在三个时间点 (Y
) 评估的结果 (period
)。第一个时期是 0
,是随机化之前记录的基线分数。
我模拟了一个数据集:
library(simstudy)
set.seed(420)
data <- defData(varname = "arm", dist = "binary", formula = 0.5)
data <- defData(data, varname = "Y0", dist = "normal", formula = 10, variance = 2)
data <- defData(data, varname = "Y1", dist = "normal", formula = "Y0 + 5 + 5 * arm", variance = 2)
data <- defData(data, varname = "Y2", dist = "normal", formula = "Y0 + 10 + 5 * arm", variance = 2)
data <- genData(500, data)
data <- addPeriods(data, nPeriods = 3, idvars = "id", timevars = c("Y0", "Y1", "Y2"), timevarName = "Y")
data <- data %>%
mutate(
arm = as.factor(arm),
period = as.factor(period)
)
每个时间点和每个臂内的分数汇总统计数据为:
data %>%
group_by(arm, period) %>%
summarise(mean = mean(Y))
# A tibble: 6 × 3
# Groups: arm [2]
arm period mean
<fct> <fct> <dbl>
1 0 0 10.1
2 0 1 14.9
3 0 2 20.2
4 1 0 10.2
5 1 1 20.2
6 1 2 25.2
我希望使用约束基线模型对每个时间点的平均差异进行建模(https://bmjopen.bmj.com/content/6/12/e013096)。简而言之,这适合一个具有时间段固定效应以及时间段和治疗组之间相互作用的回归模型(治疗没有主效应)。这具有“约束”基线的效果,或者假设它来自相同的分布(大多数随机对照试验中的合理假设)。包括主效应会导致各组之间的基线有所不同,在这种情况下这是不希望的。
我尝试使用
lme4
包来拟合这个模型:
model1 <- lmer(Y ~ period + period:arm + (1 | id), data = data)
但是,在拟合
emmeans
时,我注意到基线得分(应该相同)反而等于组特定平均值:
emmeans(model1, ~ arm | period)
period = 0:
arm emmean SE df lower.CL upper.CL
0 10.1 0.113 943 9.83 10.3
1 10.2 0.113 943 9.97 10.4
period = 1:
arm emmean SE df lower.CL upper.CL
0 14.9 0.113 943 14.71 15.1
1 20.2 0.113 943 20.01 20.5
period = 2:
arm emmean SE df lower.CL upper.CL
0 20.2 0.113 943 19.94 20.4
1 25.2 0.113 943 25.01 25.5
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
发生这种情况是因为模型中指定的公式已扩展到每个时间点(包括基线)的交互:
colnames(model.matrix(model1))
[1] "(Intercept)" "period1" "period2" "period0:arm1" "period1:arm1" "period2:arm1"
我期待看到:
[1] "(Intercept)" "period1" "period2" "period1:arm1" "period2:arm1"
(没有
period0:arm1
术语)。
我的问题是:如何修改这个模型,使其实际上是一个约束基线模型?
您的
period
变量是一个因子,因此当您在 lme4
中对其进行建模时,它是虚拟编码的,因此除了 period0
和 period1
之外,还包含 period2
的交互作用。据我所知,您不能要求 lme4
仅包含变量的某些值。您必须剔除数据以排除 period0 并仅针对基线模型运行该数据。
baseline_data <- data[data$period != 0, ] # filter out period0
baseline_model <- lmer(Y ~ period + period:arm + (1 | id), data = baseline_data)
在这种情况下,当您运行
colnames(model.matrix(baseline_model))
时,列表中不会有 period1
,因为截距为 period1
(截距之前为 period0
)。