R语言 | 回归中自变量的交互效应及R语言计算示例
Posted 生物空间站
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了R语言 | 回归中自变量的交互效应及R语言计算示例相关的知识，希望对你有一定的参考价值。
例如前述就是一种典型的没有交互效应的回归，预测变量和响应变量间关系可简单概括为预测变量的加权和形式。以二元线性回归为例，响应变量（Y）和两种预测变量（X1和X2）就可以通过下式描述，其中β1和β2分别为X1和X2的回归系数（斜率），β0是截距：
此时若在上述回归中引入X1和X2的交互效应，则交互效应表示为X1和X2的乘积。此时X1和X2之间的交互效应称为双因素交互（two-way interaction），因为它来自于两个预测变量之间的交互效应：
交互效应也可以为高阶的，但分析时通常避免引入，因为它们非常难以解释。例如三因素交互（three-way interaction）这种，过于复杂了：
注意，当在回归中引入交互效应时，回归都不再表示预测变量和响应变量间的线性关系。
下文通过两个简单示例，展示R语言执行带交互项回归的方法。
示例数据、R代码等，可见网盘附件（提取码，3vk7）：
https://pan.baidu.com/s/1k47Uf7fhM8UxdTafMiUJhA
连续型预测变量间的交互效应
先来看一种最简单类型，仅涉及两组连续变量的交互效应。
示例数据
网盘附件“beetles.txt”来自某项观察研究数据，一区域中某种类型甲虫对猎物的平均捕食频率（attack rates），以及甲虫附近的平均猎物数量密度（prey number）和气温环境（temperature），并期望查看后两种环境特征对甲虫捕食频率的影响。
一个简单的二元线性回归
如前文中所提到的，尽管实际中变量间关系并非线性这样简单，但由于线性关系比较直观且易于理解，大多数分析中首选将问题简化为线性模型去解释。因此，不妨首先拟合甲虫捕食频率与猎物数量和气温的线性回归。
#读取示例数据
dat_beetles <- read.delim('beetles.txt', sep = ' ')
#拟合甲虫捕食频率与周围猎物数量和气温的关系
#不考虑猎物数量和气温的交互效应时，不放首先通过二元线性回归来简单描述三者间关系
fit1 <- lm(attack_rate~prey_num+temperature, data = dat_beetles)
summary(fit1) #不带交互项回归的简单统计
可看到全模型p<<0.05是非常显著的，R2=0.918非常高，并且各预测变量也很显著，表明甲虫捕食频率与猎物数量和气温之间存在很强的线性关系。通过各预测变量的斜率均大于0，很容易推断甲虫捕食频率随猎物数量和气温的增加而增加。
将斜率各截距项提取出，组合回归式为：
气温升高时，甲虫捕食频率增加，可归因于气温回暖调动了昆虫这种冷血动物的活动；猎物的数量增加后，甲虫捕食次数也观察到随之增加的一致性，这也是可以预期的。
添加预测变量的交互效应
然后思考一个问题，甲虫猎物通常是更小的昆虫或其它节肢动物，它们的活动也可能会受到气温的影响：气温的回升通常有助于它们活动频率和范围的增加，以及种群密度的提升，因此温度也可能会影响甲虫猎物数量对甲虫捕食频率的效应。
出于这种考虑，不妨假设甲虫猎物数量和温度之间存在非独立效应，尝试在上述回归中添加气温与猎物数量的交互效应，并观察交互效应是否也能在一定程度上解释甲虫捕食频率。
lm()中，交互效应直接表示为两个预测变量的相乘项。需要注意的是，添加交互项后由于引入了预测变量的乘积，因此不再是线性方程。
#如果考虑猎物数量和气温的交互效应，公式中使用“*”表示二者交互，此时不再是线性回归
fit2 <- lm(attack_rate~prey_num*temperature, data = dat_beetles)
summary(fit2) #带交互项回归的简单统计
与上文的二元线性回归相比，引入交互效应后会极大改变所有系数的原有解释，包括斜率、截距项以及显著性等。
然而，观察到气温与猎物数量的交互效应并不显著，交互项的添加也没有助于模型精度的提升（带交互项R2=0.918，先前不带交互项为R2=0.916）。
据此可知，原来的假设“甲虫猎物数量和温度之间存在交互效应”是不成立的。
若假设上述交互效应是显著的，那么回归方程就可以写作为：
这里只是举例说一下带交互项回归方程式的写法，实际上是不显著的，所以该式不可用，因此还是使用上文的二元线性回归就可以了。
作图查看交互效应
不妨作图查看一下交互效应，可以帮助我们理解：
（1）预测变量间交互作用显著时，怎样体现的；
（2）预测变量间交互作用不明显时，又是怎样的分布状态。
#作图查看交互效应
#effects 包提供了强大而便捷的方法
library(effects)
plot(effect('prey_num*temperature', fit2), multiline=TRUE)
这是借助effects包中的方法，查看的甲虫猎物数量和温度之间的交互效应。该图度量了一个因素不同水平的效应变化依赖于另一个因素的水平的程度。如果交互作用不显著，图中趋势线之间会相互接近平行；如果交互作用越显著，则相交的程度就越明显。
由于上述已知甲虫猎物数量和温度之间的交互效应不明显，也就是甲虫猎物数量的不同水平对甲虫捕食频率的效应变化并不依赖于温度的水平，所以图中的线都是相互平行的。
带类别型预测变量的交互效应
上述示例中交互的两预测变量均为连续型变量。现在再来看另一个示例，当预测变量中存在类别型变量时，如何考虑交互效应。
示例数据
网盘附件“plant.txt”来自某项调查研究数据，测量了生长于两个不同温度（低、高）环境中的某种植物的高度，以及土壤中的氮浓度，期望了解土壤氮浓度及温度对植物生长的影响。
数据中，响应变量plant代表植物高度，预测变量N代表土壤中氮浓度，它们均为连续变量。另一预测变量temperature（温度）并非具体的测量数值，而是以Low（低温）和High（高温）表示，它是一列类别变量。
考虑类别预测变量的交互效应
该示例数据与前文“”中使用的数据一致，前文已初步简介了如何将类别预测变量添加到回归中，并建立线性回归探讨了植物生长高度与其所处环境中土壤氮浓度及温度的关系，提示土壤氮含量高时更有利于植物生长，而低温将抑制植物生长。
不过，当观察植物生长与土壤氮浓度的关系时，其实也能发现这种关系在高温和低温环境下是有区别的：与低温相比，较高温度下植物似乎对土壤氮利用率更高（下图中，红色代表高温关系，蓝色代表低温关系）。因此推测，土壤氮浓度对植物生长的效应随温度的改变而改变，二者存在交互效应。
考虑这些情况，接下来尝试在回归中添加温度与土壤氮浓度的交互效应。
#读取示例数据
dat_plant <- read.delim('plant.txt', sep = ' ')
#探讨植物生长高度与其所处环境中土壤氮浓度及温度的关系
#存在类别型预测变量时，lm() 中推荐使用 factor(类别型预测变量) 写入公式
#不考虑土壤氮浓度及温度的交互效应时，一个简单的带类别型预测变量的二元线性回归
fit1 <- lm(plant~N+factor(temperature), data = dat_plant)
summary(fit1) #不带交互项回归的简单统计
#如果考虑土壤氮浓度及温度的交互效应，同样地，公式中使用“*”表示二者交互，此时不再是线性回归
fit2 <- lm(plant~N*factor(temperature), data = dat_plant)
summary(fit2) #带交互项回归的简单统计
显示全模型p<<0.05是非常显著的，同时各预测变量的效应也很显著，包括交互效应在内，即表明温度与土壤氮浓度交互效应明显。可以看到，在添加了交互项后回归整体R2有所提高，达到很理想的数值。
在结果中，连续型预测变量的斜率理解较为简单，无需多提，可参考前文“”中的描述。
对于类别型预测变量的理解，前文“”中已经提到了，当回归中包含一个类别（或因子）预测变量时，类别变量被转换成取值为0或者1的名义变量（dummy variables）。例如，类别变量temperature有两个等级，Low（低温）和High（高温），R因此构造了一个名义变量，在上述结果中变量名为factor(temperature)Low，当temperature=Low时其值为1，temperature=High时其值为0。factor(temperature)Low的斜率代表了Low和High的斜率差值，在结果中它是显著的。
再结合结果中求出的斜率和截距，这个带交互效应的回归可以表示为：
其中plant是植物生长高度；N是土壤氮含量，一组连续变量；temperature为类别变量，temperature=Low时其值为1，temperature=High时其值为0。因此，不妨将上式按类别变量拆解为两个回归：
作图查看交互效应
最终，该示例中温度与土壤氮浓度的交互效应可以简化为两个线性方程来表示，可分别绘制散点图表示它们。
#作图查看交互效应
library(ggplot2)
p <- ggplot(dat_plant, aes(x = N, y = plant, color = factor(temperature))) +
geom_point() +
scale_color_manual(values = c('red', 'blue'), limit = c('High', 'Low')) +
stat_function(fun = function(x) 2.99*x - 1.17, color = 'blue') + #添加自定义回归方程式
stat_function(fun = function(x) 4.50*x - 1.10, color = 'red') +
theme(panel.grid = element_blank(), legend.key = element_blank(),
panel.background = element_rect(color = 'black', fill = 'transparent')) +
labs(x = 'N', y = 'plant', color = '')
#将回归方程添加在图中
label_fit <- data.frame(
formula_high = sprintf('italic(Y) == %.2f*italic(X) - %.2f', 4.50, 1.10),
formula_low = sprintf('italic(Y) == %.2f*italic(X) - %.2f', 2.99, 1.17))
p +
geom_text(x = 1, y = 40, aes(label = formula_high), data = label_fit, parse = TRUE,
hjust = 0, color = 'red', show.legend = FALSE) +
geom_text(x = 1, y = 35, aes(label = formula_low), data = label_fit, parse = TRUE,
hjust = 0, color = 'blue', show.legend = FALSE)
关于是否考虑保留交互效应
这是一个有待思考的问题。一方面，要保证回归的解释度足够充分，R2不能过低；另一方面，还要使回归易于解释，即简约性原则。上文两个示例中，仅第二个示例的交互效应显著，就继续以它为例作一些简单概括。
尽管添加温度与土壤氮浓度的交互效应后提升了回归的精度（R2=0.998），但事实上即使不考虑交互项，仅温度与土壤氮浓度的二元线性回归精度也已经很出色了（R2=0.96）。如果您不是很确定要不要对交互效应进行取舍，可参考前文中提到的方法，不妨检验比较上述两个嵌套模型的拟合优度，或者AIC（Akaike Information Criterion，赤池信息准则）等值综合衡量回归复杂性与拟合优度的关系，这些均可给出一些提示。
#anova() 检验两个嵌套模型的拟合优度
#比较不含交互项与含交互项的预测效果是否一样好
anova(fit1, fit2)
#AIC 评估两个回归复杂性与拟合优度的关系
#AIC 值较小的回归优先选择，表明较少的预测变量已经获得了足够的拟合度
AIC(fit1, fit2)
anova()检验后p<<0.05，表明交互效应的添加有助于精度的显著提升。类似地，AIC值也显示加入交互效应后更佳。
当然，这只是机器算法自动评估的，它认为考虑温度与土壤氮浓度的交互效应更好。实际情况中还需人为选择，您可以认同机器评估；但如果觉得交互效应思考过多复杂因素（尽管它们可能重要），并且仅温度与土壤氮浓度的二元线性回归的R2=0.96已经很充分了，不考虑交互效应也未尝不可。
此外，以上示例均展示的双变量间的交互效应，它们并不是很难理解。然而更高阶的情况，例如三因素或以上的交互，由于过于复杂难以解释问题，实际应用中尽量避免。
参考资料
以上是关于R语言 | 回归中自变量的交互效应及R语言计算示例的主要内容，如果未能解决你的问题，请参考以下文章
R语言生存分析之COX比例风险模型构建及C-index计算示例
R语言使用lm函数拟合多元线性回归模型假定预测变量之间有交互作用R语言使用effects包的effect函数查看交互作用对于回归模型预测响应变量的影响
R语言使用lm函数拟合多元线性回归模型假定预测变量之间有交互作用R语言使用effects包的effect函数查看交互作用对于回归模型预测响应变量的影响
R语言使用pwr包的pwr.f2.test函数对线性回归模型进行效用分析（power analysis）在已知效应量(effect size)显著性水平效用值的情况下计算需要的样本量