基于R语言mtcars数据集用虚拟变量复现交互水平下主效应计算

我们以R内置mtcars数据集为例，完整操作流程如下： #

核心原理 #

含cyl和gear交互项的回归模型表达式为：

mpg = β0 + β_cyl * cyl + β_gear * gear_dummies + β_inter * cyl * gear_dummies
当我们将目标水平gear=4设为分类变量gear的参照组时，其余档位的虚拟变量取值均为0，此时cyl的主效应系数就刚好等于gear=4时cyl对mpg的边际效应，对应的p值就是该效应的显著性检验结果。

具体操作步骤 #

• 第一步：预处理数据集，将gear转为因子并指定gear=4为参照水平

# 载入数据集
car_data <- mtcars
# 重编码gear为因子，指定4为第一水平（参照组）
car_data$gear <- factor(car_data$gear, levels = c(4, 3, 5))

如果不想用R的因子自动生成虚拟变量，也可以手动构造虚拟变量：

# 手动生成非目标水平的虚拟变量
car_data$gear3 <- ifelse(car_data$gear == 3, 1, 0)
car_data$gear5 <- ifelse(car_data$gear == 5, 1, 0)
# 手动生成交互项
car_data$cyl_gear3 <- car_data$cyl * car_data$gear3
car_data$cyl_gear5 <- car_data$cyl * car_data$gear5

• 第二步：构建交互项回归模型，直接读取目标效应
  如果用了因子重编码的方式，运行以下代码即可：

m_dummy <- lm(mpg ~ cyl * gear, data = car_data)
summary(m_dummy)$coefficients["cyl", ]

如果是手动生成虚拟变量的方式，运行以下代码：

m_dummy_manual <- lm(mpg ~ cyl + gear3 + gear5 + cyl_gear3 + cyl_gear5, data = car_data)
summary(m_dummy_manual)$coefficients["cyl", ]

以上两种方式输出的cyl系数、标准误、p值，和multcomp::glht()计算的gear=4时cyl的效应结果完全一致。 #

内容的提问来源于stack exchange，提问作者Bob