解释回归模型中的交互

Question

我希望是一个简单的问题。

我有一个实验设计，我在其中测量两组的一些反应（比如说血压）：对照组和受影响组，两组都接受了三种治疗：t1、t2、t3。数据在任何意义上都没有配对。

这是一个示例数据：

set.seed(1)
df <- data.frame(response = c(rnorm(5,10,1),rnorm(5,10,1),rnorm(5,10,1),
                         rnorm(5,7,1),rnorm(5,5,1),rnorm(5,10,1)),
                 group = as.factor(c(rep("control",15),rep("affected",15))),
                 treatment = as.factor(rep(c(rep("t1",5),rep("t2",5),rep("t3",5)),2)))

我感兴趣的是量化每种治疗相对于对照组对受影响组的影响。我该如何建模，比如使用线性模型（例如 R 中的 lm）？

我的想法是不是错了：

lm(response ~ 0 + treatment * group, data = df)

相当于：

lm(response ~ 0 + treatment + group + treatment:group, data = df)

不是我需要的吗？我认为在这个模型中，治疗：组交互项是相对于所有基线组和基线治疗测量的平均值。

因此我认为这个模型：

lm(response ~ 0 + treatment:group, data = df)

是我需要的，但它正在量化治疗和群体互动术语的每个组合：治疗t1：群体控制治疗t1：受群体影响的治疗t2：群体控制治疗t2：受群体影响的治疗t3：群体控制治疗t3：受群体影响

所以也许这个模型：

lm(response ~ 0 + treatment + treatment:group, data = df)

是正确的吗？

虽然除了量化每种治疗组合和受群体影响的交互项之外，它还量化了每种治疗的效果。我不确定在此模型中与每个治疗和受组影响的交互项进行比较的基线是什么。

我们将不胜感激。

另外，假设我进行了第四次治疗，这实际上是两种治疗的组合，比如 t1+t3，我不知道它们的组合效果的预期是什么：加法/减法或协同作用。有什么办法可以合并吗？

Answer 1

您的第一个规范很好。

lm(response ~  0 + treatment * group, data = df)

Call:
lm(formula = response ~ 0 + treatment * group, data = df)

Coefficients:
         treatmentt1               treatmentt2               treatmentt3  
               7.460                     5.081                     9.651  
        groupcontrol  treatmentt2:groupcontrol  treatmentt3:groupcontrol  
               2.670                     2.384                    -2.283 

第一个系数 7.460 表示参与者同时接受 t1 治疗和受影响时发生的效果。从左到右，第二个系数 5.081 表示参与者同时受到 t2 治疗和受影响等...

因此，例如，当参与者被 t2 处理并且在控制中，效果是 5.081 + 2.384。

如果我在做这个分析，我会保留拦截。

Call:
lm(formula = response ~ treatment * group, data = df)

Coefficients:
         (Intercept)               treatmentt2               treatmentt3  
               7.460                    -2.378                     2.192  
        groupcontrol  treatmentt2:groupcontrol  treatmentt3:groupcontrol  
               2.670                     2.384                    -2.283  

现在第二个系数从左到右表示接受 t2 治疗并受影响的参与者相对于接受 t1 治疗并受影响的参与者的效果。看到这个通知，7.460 - 2.378 = 5.081（第一个规范中的第二个系数）。我喜欢这种方法，因为它更容易解释相关效果。

大家都说@MrFlick 是对的。这是交叉验证的问题。

Answer 2

交互项告诉您组间的差异取决于治疗，即 t1、t2 和 t3 的受影响和控制之间的差异不同。

我会模拟拦截。

lm(response ~ group + treatment + group:treatment, data=df)

在获得重要的交互项后，我会使用 t.tests 进一步调查并帮助解释。

可以看出，相互作用是由 t2 相对于其他因素的更大影响驱动的。

library(data.table)
library(dplyr)
library(ggplot2)

set.seed(1)
df <- data.frame(response = c(rnorm(5,10,1),rnorm(5,10,1),rnorm(5,10,1),rnorm(5,7,1),rnorm(5,5,1),rnorm(5,10,1)),
             group = as.factor(c(rep("control",15),rep("affected",15))),
             treatment = as.factor(rep(c(rep("t1",5),rep("t2",5),rep("t3",5)),2)))

# t tests of the desired comparisons to see if there is a difference and get 95% confidence intervals
t.test(df$response[df$treatment=="t1"] ~ df$group[df$treatment=="t1"])
t.test(df$response[df$treatment=="t2"] ~ df$group[df$treatment=="t2"])
t.test(df$response[df$treatment=="t3"] ~ df$group[df$treatment=="t3"])

# plot 95% C.I.
ci_plot <- matrix(nrow=3, ncol=3)
ci_plot <- as.data.frame(ci_plot)
colnames(ci_plot) <- c("treatment", "lci", "uci")

ci_plot[,1] <- c("t1", "t2", "t3")
ci_plot[,3] <- c(t.test(df$response[df$treatment=="t1"] ~ df$group[df$treatment=="t1"])$conf.int[1],
             t.test(df$response[df$treatment=="t2"] ~ df$group[df$treatment=="t2"])$conf.int[1],
             t.test(df$response[df$treatment=="t3"] ~ df$group[df$treatment=="t3"])$conf.int[1])
ci_plot[,4] <- c(t.test(df$response[df$treatment=="t1"] ~ df$group[df$treatment=="t1"])$conf.int[2],
             t.test(df$response[df$treatment=="t2"] ~ df$group[df$treatment=="t2"])$conf.int[2],
             t.test(df$response[df$treatment=="t3"] ~ df$group[df$treatment=="t3"])$conf.int[2])

ggplot(ci_plot, aes(x=treatment, y=uci)) +
    geom_errorbar(aes(ymin=uci, ymax=lci), width=0.5, position=position_dodge(0.9), weight=0.5) +
    xlab("Treatment") +
    ylab("Change in mean relative to control (95% C.I.)") +
    theme_bw() +
    theme(panel.border = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          axis.line = element_line(colour = "black"),
          axis.text.x = element_text(angle = 90, hjust = 1))