含双交互项的回归模型边际效应可视化问题咨询

我用estimatr包的lm_robust函数运行回归模型，纳入多个控制变量、gender_f分别与SDOsummary、RWAsummary的两个交互项，同时设置国家固定效应并在个体层面聚类标准误。因为多数工具包不兼容稳健标准误，我手动提取系数与协方差矩阵计算边际效应，但得到的两条边际效应曲线y轴起始点完全相同，这明显不合理。

回归模型核心代码（完整代码见下文）：

mod_sdo_rwa <- lm_robust(Vote ~ gender_f + occ_f + age_f + experience_f + children_f + 
    ideology_f + origin_f + SDOsummary + gender_f*SDOsummary + RWAsummary + gender_f*RWAsummary, 
    data=FullLaunch_New, weights = weight, fixed_effects = country, clusters = Individ_f, se_type="stata")

其中gender_f是二元变量，SDOsummary和RWAsummary均为0-1区间的连续变量。

完整代码：

mod_sdo_rwa <- lm_robust(Vote ~ gender_f + occ_f + age_f + experience_f + children_f +    
    ideology_f + origin_f + SDOsummary + gender_f*SDOsummary + RWAsummary + gender_f*RWAsummary, 
    data=FullLaunch_New, weights = weight, fixed_effects = country, clusters = Individ_f,   
    se_type="stata")

beta.hatRWA <- coef(mod_sdo_rwa) 
covRWA <- vcov(mod_sdo_rwa)
z0RWA <- seq(min(FullLaunch_New$RWAsummary, na.rm=TRUE), max(FullLaunch_New$RWAsummary,            
    na.rm=TRUE), length.out = 1000)
dy.dxRWA <- beta.hatRWA["gender_fFemale"] + beta.hatRWA["gender_fFemale:RWAsummary"]*z0RWA
se.dy.dxRWA <- sqrt(covRWA["gender_fFemale", "gender_fFemale"] +      
    z0RWA^2*covRWA["gender_fFemale:RWAsummary", "gender_fFemale:RWAsummary"] +    
    2*z0RWA*covRWA["gender_fFemale", "gender_fFemale:RWAsummary"])
uprRWA <- dy.dxRWA + 1.96*se.dy.dxRWA
lwrRWA <- dy.dxRWA - 1.96*se.dy.dxRWA

beta.hatSDO <- coef(mod_sdo_rwa) 
covSDO <- vcov(mod_sdo_rwa)
z0SDO <- seq(min(FullLaunch_New$SDOsummary, na.rm=TRUE), max(FullLaunch_New$SDOsummary,    
    na.rm=TRUE), length.out = 1000)
dy.dxSDO <- beta.hatSDO["gender_fFemale"] + beta.hatSDO["gender_fFemale:SDOsummary"]*z0SDO
se.dy.dxSDO <- sqrt(covSDO["gender_fFemale", "gender_fFemale"] +     
    z0SDO^2*covSDO["gender_fFemale:SDOsummary", "gender_fFemale:SDOsummary"] + 2*z0SDO*covSDO["gender_fFemale", "gender_fFemale:SDOsummary"])
uprSDO <- dy.dxSDO + 1.96*se.dy.dxSDO
lwrSDO <- dy.dxSDO - 1.96*se.dy.dxSDO

Line <- ggplot(data=NULL,aes(x=z0RWA, y=dy.dxRWA)) +
labs(x="Right Wing Authoritarianism/SDO",y="Estimated Coefficient of Gender") +
geom_line(aes(z0RWA, dy.dxRWA, color="red"), size = 1) +
geom_line(aes(z0RWA, lwrRWA, color="red"), size = 1, 
        linetype = 2) +
geom_line(aes(z0RWA, uprRWA, color="red"), size = 1, 
        linetype = 2) +
geom_hline(yintercept=0) + 
geom_line(data=NULL,aes(x=z0SDO, y=dy.dxSDO, color="blue"), size = 1) + 
geom_line(aes(z0SDO, lwrSDO, color="blue"), size = 1, 
        linetype = 2) +
geom_line(aes(z0SDO, uprSDO, color="blue"), size = 1, 
        linetype = 2)

PanelB <- Line + theme(legend.position="right") +
scale_color_discrete(name="", labels = c("SDO", "RWA")) +
ylim(-.04, .1) +
theme_bw() +
theme(axis.title.x=element_text(size=22),           
    axis.title.y=element_text(size=20),legend.title=element_text(size=20), legend.text=element_text(size=18), legend.position="bottom")

生成的图表：

我的两个问题：

1. 为什么两条曲线的y轴起始点完全相同？
2. 我想可视化性别对Vote的边际效应在SDO和RWA各自取值范围内的变化，但现有工具好像只能单独展示其中一个调节变量的分组效应，是不是我漏了什么方法？

1. 曲线起始点相同的原因 #

你计算边际效应的逻辑有误。在同时包含两个交互项的模型中，性别对Vote的边际效应不能只单独计算单个交互项的线性组合，必须考虑另一个调节变量的取值。

你当前的计算逻辑是：

• RWA组边际效应：dy.dxRWA = gender_fFemale + gender_fFemale:RWAsummary * z0RWA
• SDO组边际效应：dy.dxSDO = gender_fFemale + gender_fFemale:SDOsummary * z0SDO

这相当于默认把另一个调节变量的取值设为0，而如果SDOsummary和RWAsummary的最小值都是0，两条曲线在x=0时的截距都会等于gender_fFemale的系数，所以起始点完全重合。

正确的计算应该固定另一个调节变量的参考值（比如均值、中位数）：

• 计算RWA取值下的边际效应（固定SDO为均值mean_sdo）：

  dy.dxRWA <- beta.hatRWA["gender_fFemale"] + 
              beta.hatRWA["gender_fFemale:RWAsummary"]*z0RWA +
              beta.hatRWA["gender_fFemale:SDOsummary"]*mean_sdo
  

• 计算SDO取值下的边际效应（固定RWA为均值mean_rwa）：

  dy.dxSDO <- beta.hatSDO["gender_fFemale"] + 
              beta.hatSDO["gender_fFemale:SDOsummary"]*z0SDO +
              beta.hatSDO["gender_fFemale:RWAsummary"]*mean_rwa
  

同时，标准误的计算也要加入交叉项的协方差，不能只考虑单个交互项的方差。

2. 同时可视化两个调节变量的边际效应 #

推荐两种实现方式：

方式1：分面板展示（更清晰） #

用facet_wrap把两个调节变量的边际效应分成子图，避免混淆：

# 先计算参考值
mean_sdo <- mean(FullLaunch_New$SDOsummary, na.rm=TRUE)
mean_rwa <- mean(FullLaunch_New$RWAsummary, na.rm=TRUE)

# 整理RWA边际效应数据
df_rwa <- data.frame(
  x = z0RWA,
  marginal_effect = beta.hatRWA["gender_fFemale"] + 
                    beta.hatRWA["gender_fFemale:RWAsummary"]*z0RWA +
                    beta.hatRWA["gender_fFemale:SDOsummary"]*mean_sdo,
  upper = dy.dxRWA + 1.96*sqrt(covRWA["gender_fFemale","gender_fFemale"] + 
                               z0RWA^2*covRWA["gender_fFemale:RWAsummary","gender_fFemale:RWAsummary"] +
                               mean_sdo^2*covRWA["gender_fFemale:SDOsummary","gender_fFemale:SDOsummary"] +
                               2*z0RWA*covRWA["gender_fFemale","gender_fFemale:RWAsummary"] +
                               2*mean_sdo*covRWA["gender_fFemale","gender_fFemale:SDOsummary"] +
                               2*z0RWA*mean_sdo*covRWA["gender_fFemale:RWAsummary","gender_fFemale:SDOsummary"]),
  lower = dy.dxRWA - 1.96*sqrt(...), # 同上标准误计算
  moderator = "RWA"
)

# 整理SDO边际效应数据
df_sdo <- data.frame(
  x = z0SDO,
  marginal_effect = beta.hatSDO["gender_fFemale"] + 
                    beta.hatSDO["gender_fFemale:SDOsummary"]*z0SDO +
                    beta.hatSDO["gender_fFemale:RWAsummary"]*mean_rwa,
  upper = ..., # 对应修正后的标准误计算
  lower = ...,
  moderator = "SDO"
)

# 合并绘图
df_combined <- rbind(df_rwa, df_sdo)
ggplot(df_combined, aes(x=x, y=marginal_effect, color=moderator)) +
  geom_line(size=1) +
  geom_ribbon(aes(ymin=lower, ymax=upper, fill=moderator), alpha=0.2, color=NA) +
  geom_hline(yintercept=0, linetype=2) +
  facet_wrap(~moderator, scales="free_x") +
  labs(x="调节变量取值", y="性别对Vote的边际效应") +
  theme_bw()

方式2：用专业工具包简化操作 #

marginaleffects包支持estimatr的lm_robust模型，能直接计算并可视化条件边际效应，不用手动计算：

library(marginaleffects)

# 计算RWA调节下的性别边际效应（固定SDO为均值）
me_rwa <- plot_cme(mod_sdo_rwa,
                   variables = "gender_f",
                   condition = list(RWAsummary = seq(min(FullLaunch_New$RWAsummary, na.rm=T), max(...), length.out=100),
                                    SDOsummary = mean(FullLaunch_New$SDOsummary, na.rm=T)))

# 计算SDO调节下的性别边际效应（固定RWA为均值）
me_sdo <- plot_cme(mod_sdo_rwa,
                   variables = "gender_f",
                   condition = list(SDOsummary = seq(min(FullLaunch_New$SDOsummary, na.rm=T), max(...), length.out=100),
                                    RWAsummary = mean(FullLaunch_New$RWAsummary, na.rm=T)))

# 合并或分图展示即可

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