R中plm固定效应模型交互项的结果展示、边际效应计算与可视化问询

背景与模型代码 #

我正在使用R语言的plm包处理面板数据，估计包含交互项的固定效应模型。目前已通过plm得到模型结果并使用stargazer进行展示，但希望以更精细化的方式呈现结果（重点聚焦交互项）、计算交互项的边际效应及其p值并纳入结果表格，同时绘制边际效应图。

模型示例代码：

fe_model_interaction <- plm(life_satisfaction ~ 
                                  employment_level_Full_Time*care_low +
                                  employment_level_Full_Time*care_high +
                                  employment_level_Part_Time*care_low +
                                  employment_level_Part_Time*care_high+
                                  other_control_variables,
                                  data = data_analyse_mother, 
                                  index = c("pid", "syear"), 
                                  model = "within")

问题1：替代stargazer的精细化模型结果展示 #

推荐使用modelsummary包，它比stargazer更灵活，支持自定义系数标签、高亮交互项、灵活配置统计量，输出格式覆盖HTML/LaTeX/Word，表格整洁专业：

1. 安装并加载包：

install.packages("modelsummary")
library(modelsummary)

2. 核心调用示例（聚焦交互项）：

modelsummary(fe_model_interaction, 
             # 配置要展示的统计量
             statistic = c("标准误 = {std.error}", "p值 = {p.value}"),
             # 自定义系数标签，让交互项更易读
             coef_map = c(
               "employment_level_Full_Time:care_low" = "全职 × 低照料",
               "employment_level_Full_Time:care_high" = "全职 × 高照料",
               "employment_level_Part_Time:care_low" = "兼职 × 低照料",
               "employment_level_Part_Time:care_high" = "兼职 × 高照料",
               # 可补充控制变量的标签
               "other_control_variables1" = "控制变量1"
             ),
             # 高亮交互项行，突出重点
             highlight = list(rows = grep("×", names(coef_map)), color = "lightgray"),
             # 表格标题
             title = "固定效应模型结果（含交互项）")

若偏好LaTeX输出，也可使用texreg包，支持更细粒度的格式定制：

install.packages("texreg")
library(texreg)

texreg(fe_model_interaction, 
       custom.coef.names = c("控制变量1", "全职 × 低照料", "全职 × 高照料", "兼职 × 低照料", "兼职 × 高照料"),
       caption = "固定效应模型结果（含交互项）",
       caption.above = TRUE,
       include.ci = FALSE,
       digits = 3)

问题2：计算交互项边际效应并纳入表格 #

你的模型直接估计了分组交互项（全职/兼职 × 低/高照料），这些系数本身就是对应分组下的边际效应。若需计算特定变量的条件边际效应（比如全职状态下，照料水平对生活满意度的边际影响），可使用margins包：

1. 安装并加载包：

install.packages("margins")
library(margins)

2. 计算并提取边际效应：

# 计算全职/兼职状态下，低/高照料的边际效应
marg_eff <- margins(fe_model_interaction, 
                    variables = c("employment_level_Full_Time", "employment_level_Part_Time"),
                    at = list(care_low = c(0, 1), care_high = c(0, 1)))
# 转为数据框方便处理
marg_eff_df <- summary(marg_eff)

3. 将边际效应合并到原模型结果表格：

# 整理边际效应为表格行
marg_rows <- data.frame(
  term = c("全职（低照料）边际效应", "全职（高照料）边际效应", 
           "兼职（低照料）边际效应", "兼职（高照料）边际效应"),
  estimate = marg_eff_df$estimate,
  `标准误` = marg_eff_df$std.error,
  `p值` = marg_eff_df$p.value,
  check.names = FALSE
)

# 用modelsummary合并原模型结果与边际效应
modelsummary(fe_model_interaction,
             statistic = c("标准误 = {std.error}", "p值 = {p.value}"),
             coef_map = coef_map, # 复用问题1中的系数标签
             add_rows = marg_rows,
             title = "固定效应模型结果与边际效应")

问题3：绘制边际效应图的方法与包 #

1. margins + ggplot2（灵活定制） #

基于margins的计算结果，用ggplot2绘制带置信区间的边际效应图：

library(ggplot2)

ggplot(marg_eff_df, aes(x = factor(at), y = estimate, 
                        ymin = estimate - 1.96*std.error, 
                        ymax = estimate + 1.96*std.error)) +
  geom_pointrange(color = "#2c3e50", fatten = 1.2) +
  facet_wrap(~variable, labeller = labeller(variable = c(
    employment_level_Full_Time = "全职状态",
    employment_level_Part_Time = "兼职状态"
  ))) +
  labs(x = "照料水平（0=无/1=有）", y = "边际效应（生活满意度变化）", 
       title = "不同就业状态下的照料边际效应") +
  theme_minimal() +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red")

2. interactions包（交互项可视化专用） #

专门针对交互项设计，支持plm模型，可直接绘制交互效应的预测值图：

install.packages("interactions")
library(interactions)

# 绘制全职状态 × 低照料的交互效应图
interact_plot(fe_model_interaction, 
              pred = employment_level_Full_Time, 
              modx = care_low,
              plot.points = TRUE, # 展示原始数据点
              x.label = "全职就业状态（0=否/1=是）",
              y.label = "生活满意度预测值",
              legend.title = "低照料水平",
              main.title = "全职状态与低照料的交互效应")

# 批量绘制所有交互组合，可结合facet参数

3. emmeans包（边际均值视角） #

通过估计不同分组的边际均值，间接展示边际效应差异：

install.packages("emmeans")
library(emmeans)

# 估计就业状态×照料水平的边际均值
emm <- emmeans(fe_model_interaction, ~ employment_level | care_low + care_high)

# 绘制边际均值对比图
plot(emm, comparisons = TRUE, adjust = "bonferroni") +
  theme_minimal() +
  labs(title = "不同就业-照料组合下的生活满意度边际均值")

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