第 29 章交互作用

线性回归部分目前为止我们讨论过如何用多元回归模型来控制 (或调整) 特定的预测变量 \((X)\) 之外的变量。多元回归的目的之一是为了估计预测变量和因变量之间的回归系数的同时，保持其他 (想要被调整的) 变量不变。如此一来，其实等于是假定了无论其余的调整变量取值如何，\(X,Y\) 之间的回归系数总是相同 (绘制的回归线是一组平行线)。本章讨论的交互作用，就是探讨其中某个变量改变了 \(X,Y\) 之间的关系 (modification effect) 的情况。放宽了之前强制所有直线都平行的限制，探讨两个变量之间的线性关系是否因为某个变量而发生了质的改变。这样的关系，在流行病学中被定义为 交互作用 interaction。

本站会探讨如何利用线性回归模型分析交互作用，如何理解并解释统计忍者包输出的报告结果的意义。具体涉及的例子为：两个连续型变量，两个分类型变量，以及一个连续型，一个分类型变量之间的关系的交互作用。

29.1 两个预测变量之间的线性模型交互作用

29.1.1 交互作用线性模型的一般表达式

假如准备拟合的模型是一个因变量 \(Y\)，两个预测变量 \(X_1, X_2\)。同时模型考虑根据 \(X_2\) 的值，\(X_1, Y\) 之间关系的回归系数可以不相等 (直线的斜率不同，即会出现两条相交的回归直线)。这样的模型其实只要在原有的两个预测变量的回归线性模型中增加一个新的预测变量，新的预测变量是 \(X_1, X_2\) 的乘积即可。很简单，不是么？

\[ \begin{aligned} y_i & = \alpha + \beta_1 x_{1i} + \beta_2 x_{2i} + \beta_3 x_{3i} + \varepsilon_i \\ \text{Where, } & \varepsilon_i \sim \text{NID}(0,\sigma^2) \\ y_i & = \text{value of the dependent variable} \\ x_{1i} & = \text{value of the first predictor variable} \\ x_{2i} & = \text{value of the second predictor variable} \\ x_{3i} & = x_{1i} \times x_{2i} \\ \end{aligned} \tag{29.1} \]

为什么增加一个 \(X_1\times X_2\) 就能够分析交互作用呢 (不同直线的斜率)？要理解其中的奥妙，我们可以这样来理解：当像普通的线性回归模型那样调整了\(X_2\) 之后，也就是当\(X_2\) 固定不变时\((X_2=k)\)，回归方程@ref (eq:lm7-1)，就变成了：

\[ \begin{equation} y_i = (\alpha + \beta_2 k) + (\beta_1 + \beta_3 k)x_{1i} + \varepsilon_i \end{equation} \tag{29.2} \]

此时的 \(X_1\) 的斜率从 \(\beta_1\) 变成了 \((\beta_1 + \beta_3 k)\)，截距从 \(\alpha\) 变成了 \((\alpha + \beta_2 k)\)。

29.1.2 连续型变量和二分类变量之间的交互作用

一个连续型变脸一个二分类变量的交互作用回归方程十分容易理解 (利用哑变量建立模型)：

\[ \begin{array}{ll} y_i = \alpha + \beta_1 x_1i + \varepsilon_i & \text{ when } X_2 = 0 \\ y_i = (\alpha + \beta_2) + (\beta_1+\beta_3)x_{1i} + \varepsilon_i & \text{ when } X_2 =1 \end{array} \tag{29.3} \]

所以，\(X_2\) 取零或者取 \(1\) 代表了不同的分组，上面的回归方程就可以拟合 \(Y, X_1\) 在 \(X_2\) 的两组中不同截距，不同斜率的两条直线。其中各个参数估计，用人话来解释就是：

\(\alpha\) 是当 \(X_2 = 0\) 时的截距；
\(\alpha + \beta_2\) 是当 \(X_2 = 1\) 时的截距，所以 \(\beta_2\) 就是二分类预测变量 \(X_2\) 的两组之间截距的差；
\(\beta_1\) 是当 \(X_2 = 0\) 时的斜率；
\((\beta_1+\beta_3)\) 是当 \(X_2 = 1\) 时的截距，所以 \(\beta_3\) 就是二分类预测变量 \(X_2\) 的两组之间斜率的差。

29.1.3 两个二分类变量之间的交互作用

当两个预测变量都是二分类变量时，可以用两个哑变量来编码各自的分组，拟合下面的回归模型：

\[ \begin{array}{lll} y_i = \alpha + \varepsilon_i & \text{ when } X_1 = 0 \& X_2 = 0 & \mu_{00} \\ y_i = \alpha + \beta_1 + \varepsilon_i & \text{ when } X_1 = 1 \& X_2 =0 & \mu_{10} \\ y_i = \alpha + \beta_2 + \varepsilon_i & \text{ when } X_1 = 0 \& X_2 =1 & \mu_{01} \\ y_i = \alpha + \beta_1 + \beta_2+ \beta_3 + \varepsilon_i & \text{ when } X_1 = 1 \& X_2 = 1 & \mu_{11} \end{array} \tag{29.4} \]

如果用 \(\mu_{ij}\) 表示 \(X_1 = i, X_2 = j\) 时的总体均值 (population mean)，那么模型 (29.4) 各个参数估计及其意义为：

\(\alpha\) 是当 \(X_1 = 0, \& X_2 = 0\) 时 \(Y\) 的均值估计 \((\mu_{00})\)；
\(\alpha+\beta_1\) 是当\(X_1 = 1 \& X_2 = 0\) 时\(Y\) 的均值估计\(\mu_{10}\)，所以\(\beta_1\) 是\(X_2 = 0\) 时\(X_1\)的两组之间\(Y\) 的均值差，\(\mu_{10}-\mu_{00}\)；
\(\alpha+\beta_2\) 是当\(X_1 = 0 \& X_2 = 1\) 时\(Y\) 的均值估计\(\mu_{01}\)，所以\(\beta_2\) 是\(X_1 = 0\) 时\(X_2\)的两组之间\(Y\) 的均值差，\(\mu_{01}-\mu_{00}\)；
\(\alpha + \beta_1 + \beta_2 + \beta_3\) 是当\(X_1 = 1 \& X_2 = 1\) 时的均值估计\(\mu_{11}\)，所以\(\beta_3\) 是\(X_1 = 1\)时，\(X_2\) 的两组之间\(Y\) 的均值差\(\mu_{11}-\mu_{10}\) 减去 \(X_2 = 0\) 时，\(X_1\) 的两组之间的均值差\(\mu_{01}-\mu_{00}\)：\((\mu_{11}-\mu_{10}) - (\mu_{01}-\mu_{00})\)。

当 \(X_1\) 是连续型变量时，交互作用项的回归系数 \(\beta_3\) 的几何意义是两个回归直线斜率的差。但是本例中，两个预测变量都是二分类变量的情况下，\(\beta_3\) 的实际意义就变成了，被\(X_2\) 定义的两组\(X_1\) 之间因变量差的差，\(( \mu_{11}-\mu_{10}) - (\mu_{01}-\mu_{00})\)。

29.1.4 两个连续变量之间的交互作用

前面(Section 29.1.2) 已经讨论过，一个是连续型变量\(X_1\)，另一个是分类变量时\(X_2\)，线性回归的交互作用项回归系数的含义是因变量\(Y\) 和连续性变量\(X_1\) 在不同的\(X_2\) 组中的回归系数之差(斜率之差)。但是，当两个预测变量 \(X_1, X_2\) 都是连续型变量时，交互作用项的回归系数该如何解释呢？直观的说，此时的交互作用项回归系数应该被理解为：预测变量 \(X_2\) 每增加一个单位时，\(Y, X_1\) 之间关系的回归方程的斜率变化。为了更好地解释这个概念，我们沿用前面儿童年龄和身长预测其身高的模型 (Section 27.3.3)，加入年龄和身高的交互作用项结果如下：

growgam1 <- read_dta("../Datas/growgam1.dta")
growgam1$sex <- as.factor(growgam1$sex)

Model1 <- lm(wt ~ age + len + age*len, data=growgam1)
# or equivalently use lm(wt ~  age*len, data=growgam1)

summary(Model1)

## 
## Call:
## lm(formula = wt ~ age + len + age * len, data = growgam1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.1998 -0.6506  0.0033  0.5216  2.8945 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.519559   1.909179  -2.367  0.01895 *  
## age         -0.284900   0.104955  -2.714  0.00726 ** 
## len          0.187770   0.026808   7.004  4.4e-11 ***
## age:len      0.003398   0.001287   2.640  0.00899 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9397 on 186 degrees of freedom
## Multiple R-squared:  0.7526, Adjusted R-squared:  0.7487 
## F-statistic: 188.7 on 3 and 186 DF,  p-value: < 2.2e-16

confint(Model1)

##                     2.5 %       97.5 %
## (Intercept) -8.2859886987 -0.753130290
## age         -0.4919561735 -0.077844188
## len          0.1348837205  0.240657157
## age:len      0.0008587879  0.005937418

这里的模型中，儿童的身长和年龄相乘的部分构成了一个交互作用项。我们用这个拟合的回归方程写出当儿童年龄为 12 个月时，身长预测体重的方程：

\[ \begin{aligned} E(\text{weight|age}=12) & = (-4.5196 - 0.2849\times12) \\ & \;\;\;\;\;\;\;\;+ (0.1878 + 0.0034\times12)\times\text{Length} \\ & = -7.9384 + 0.2286\times \text{Length} \end{aligned} \]

类似地，年龄 13 个月时，预测体重的方程是：

\[ \begin{aligned} E(\text{weight|age}=13) & = (-4.5196 - 0.2849\times13) \\ & \;\;\;\;\;\;\;\;+ (0.1878 + 0.0034\times13)\times\text{Length} \\ & = -8.2233 + 0.2320\times \text{Length} \end{aligned} \]

年龄 13 个月时，预测体重的方程是：

\[ \begin{aligned} E(\text{weight|age}=14) & = (-4.5196 - 0.2849\times14) \\ & \;\;\;\;\;\;\;\;+ (0.1878 + 0.0034\times14)\times\text{Length} \\ & = -8.5082 + 0.2354\times \text{Length} \end{aligned} \]

所以你会看到每个给定的儿童年龄时的方程身长预测体重的方程都是线性方程，截距和斜率都在变化。儿童的年龄每增加 \(1\) 个月，身长和体重的相关系数增加 \(0.0034 \text{kg/cm}\)。除了回归系数，其余的数字都是不能用正常的数据来理解的 (没有儿童身长或者体重会等于零)。如果非要解释，那么需要把数据全部中心化 (Section 24.3.1)。

epiDisplay::summ(growgam1$age, graph=FALSE); epiDisplay::summ(growgam1$len, graph=FALSE)

##  obs. mean   median  s.d.   min.   max.  
##  190  16.979 16      8.337  5      36

##  obs. mean   median  s.d.   min.   max.  
##  190  76.697 76.05   7.156  60.1   95.5

growgam1$age_c <- growgam1$age-mean(growgam1$age)
growgam1$len_c <- growgam1$len-mean(growgam1$len)

Model2 <- lm(wt ~ age_c + len_c + age_c*len_c, data=growgam1)
# or equivalently use lm(wt ~  age_c*len_c, data=growgam1)
summary(Model2)

## 
## Call:
## lm(formula = wt ~ age_c + len_c + age_c * len_c, data = growgam1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.1998 -0.6506  0.0033  0.5216  2.8945 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  9.469781   0.095075   99.60  < 2e-16 ***
## age_c       -0.024275   0.017211   -1.41  0.16009    
## len_c        0.245467   0.019470   12.61  < 2e-16 ***
## age_c:len_c  0.003398   0.001287    2.64  0.00899 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9397 on 186 degrees of freedom
## Multiple R-squared:  0.7526, Adjusted R-squared:  0.7487 
## F-statistic: 188.7 on 3 and 186 DF,  p-value: < 2.2e-16

confint(Model2)

##                     2.5 %      97.5 %
## (Intercept)  9.2822176582 9.657344804
## age_c       -0.0582289785 0.009679706
## len_c        0.2070560442 0.283877252
## age_c:len_c  0.0008587879 0.005937418

你会解释上面中心化数据以后拟合的回归方程的结果，和各个参数估计的意义吗？

第 29 章 交互作用