作者们开发了一个R包，名为{InteractionPoweR}，可用于计算交互回归模型的样本量和统计功效。如果一个回归模型含有一个交互项，它可以使用下方的公式表示：其中Y为因变量，X1和X2为自变量，X1X2为交互项，βs为回归系数。这个R包希望通过统计功效分析来解决以下问题：在预设的显著性水平 (α，即p值) 下，该分析方法具有多少统计功效来检测交互项回归系数 β3 是否存在显著差异。install.packages("InteractionPoweR")
library(InteractionPoweR)
先创建一个模拟的数据，即包含两个自变量和一个交互项的回归模型的数据集，代码如下：
set.seed(123)
mydata <- generate_interaction(N = 200,
r.x1.y = 0.18, # x1和y的相关系数
r.x2.y = 0.11,
r.x1.x2 = 0.19,
r.x1x2.y = 0.2)
head(mydata)
数据比较抽象，可以直接使用一个函数进行作图，代码如下：
plot_interaction(data = mydata,
q = 2) # 简单斜率的数量，默认为2
test_interaction(data = mydata,
q = 2)
假设样本量以及其他的相关参数全部已知，α值设在0.05（这个R包的默认情况），如何计算统计功效？ 代码如下：power_interaction_r2(N = 200,
r.x1.y = 0.18,
r.x2.y = 0.11,
r.x1.x2 = 0.19,
r.x1x2.y = 0.2,
alpha = 0.05)
或者也可以使用这个R包的另外一个函数，它使用蒙特卡罗模拟来计算统计功效，代码如下： set.seed(321) # 为了让结果可以重复
power_interaction(n.iter = 1500, # 模拟的次数；这里是为了节省时间，原文作者推荐10000
alpha = 0.05,
N = 200,
r.x1.y = 0.18,
r.x2.y = 0.11,
r.x1.x2 = 0.19,
r.x1x2.y = 0.2)
下一步，稍微提升一下难度。现在假设样本量未知，需要解决的问题是找到那个最佳的样本量从而让统计功效达到一定的值，比如0.8或0.9，那该怎么做？代码如下: power_test0 <- power_interaction_r2(N = seq(100, 360, by = 10),
r.x1.y = 0.18,
r.x2.y = 0.11,
r.x1.x2 = 0.19,
r.x1x2.y = 0.2,
alpha = 0.05)
power_test0
从代码中可知，我们给N设置了一个范围，再分别计算出所对应的统计功效。从上面的结果可以找到样本量所对应的统计功效。
可以选择将上面的数据进行作图，查看统计功效为0.8所对应的样本量，代码如下： plot_power_curve(power_data = power_test0,
power_target = 0.8) # 默认
图中的那条水平线指的是统计功效=0.8，和曲线交叉点所对于x轴上的值就是要找的样本量，接近200。
如果将统计功效的水平调高到0.9，同样可以画出曲线，代码如下：
plot_power_curve(power_data = power_test0,
power_target = 0.9)
通过肉眼观察，交叉点所对应的样本量大概在260左右。因为上面的数值范围是间隔10个样本，假设现在希望更加的精确，比如精确到1个样本，那该怎么做？代码如下：
power_test1 <- power_interaction_r2(N = seq(250, 270, by = 1),
r.x1.y = 0.18,
r.x2.y = 0.11,
r.x1.x2 = 0.19,
r.x1x2.y = 0.2,
alpha = 0.05)
power_test1
从上面的数据中可知，如果希望统计功效达到0.9，所需的样本量为258。下一步，继续增加难度，假设现在有两个参数未知，以样本量和[交互项与因变量的效应量（相关系数）]为例。power_test2 <- power_interaction_r2(N = seq(100, 360, by = 10),
r.x1.y = 0.18,
r.x2.y = 0.11,
r.x1.x2 = 0.19,
r.x1x2.y = seq(0.16, 0.21, by = 0.005),
alpha = 0.05)
plot_power_curve(power_data = power_test2,
power_target = 0.8,
x = "N",
group = "r.x1x2.y")
上图相当于是将交互项与因变量的效应量（相关系数）按照数值进行分组，再画出对应的样本量和统计功效的关系图。
最后，再举一个例子，这个R包还可以将变量的信度(reliability)考虑在内。因为作者认为，如果变量的测量不可靠，那显然会影响统计功效的计算。下面的代码可以让3个变量具备不同的信度，从而可以观察信度对统计功效分析的影响，代码如下：power_test3 <- power_interaction_r2(N = 200,
r.x1.y = 0.18,
r.x2.y = 0.11,
r.x1.x2 = 0.19,
r.x1x2.y = 0.2,
alpha = 0.05,
rel.x1 = seq(0.3,1,by= 0.1),
rel.x2 = seq(0.3,1,by= 0.1),
rel.y = seq(0.3,1,by= 0.2))
plot_power_curve(power_data = power_test3,
power_target = 0.8,
x = "rel.x1",
group = "rel.x2",
facets= "rel.y")
对于希望继续深入了解这个内容的读者，可以查看下方的参考文献，或者参考原文中提供的文献，如下表： ![]()
文献[1] Table 2好啦，今天的内容就到这里。如果有帮助，记得分享给需要的人![]()
！[1] Power Analyses for Interaction Effects in Cross-Sectional Regressions
公众号核心成员担任SCI杂志Associate Editor!▌课程相关咨询可添加R师妹微信: kefu_rstats
