在含交互项的plm模型中用avg_slopes获取非工作转兼职的边际效应

问题 #

我正在使用R语言的plm包分析面板数据，在解读和计算avg_slopes函数的输出时遇到了问题。我的模型将life_satisfaction作为因变量，自变量包含employment_level（分为non_working、part_time、full_time三类）、其与current_children（分为one_child、two_children、more_than_two_children三类）的交互项，以及其他控制变量。

模型简化形式如下：

mod_2 <- plm(life_satisfaction ~ employment_level * current_children + ...,
             data = data_analyse_mother, 
             index = c("pid", "syear"), 
             model = "within")

我使用avg_slopes()来分析不同current_children水平下，就业状态转换的边际效应，尤其是non_working转full_time、part_time转full_time、non_working转part_time这几种转换。

我的avg_slopes调用代码为：

avg <- avg_slopes(mod_2, variables = "employment_level", by="current_children")

该输出提供了non_working与full_time、part_time与full_time的对比，但未直接给出non_working与part_time的对比，输出表格如下：

| term             | contrast                               | current_children         | estimate    |
|------------------|----------------------------------------|--------------------------|-------------|
| employment_level | mean(Not_Working) - mean(Full_Time)    | more_than_two_children   | -0.33049722 |
| employment_level | mean(Not_Working) - mean(Full_Time)    | one_child                | -0.21408351 |
| employment_level | mean(Not_Working) - mean(Full_Time)    | two_children             | -0.25985500 |
| employment_level | mean(Part_Time) - mean(Full_Time)      | more_than_two_children   | -0.22848561 |
| employment_level | mean(Part_Time) - mean(Full_Time)      | one_child                | -0.09376425 |
| employment_level | mean(Part_Time) - mean(Full_Time)      | two_children             | -0.06093251 |

我需要获取non_working转part_time的效应对比，尝试了成对假设但无效，请问如何用avg_slopes实现这一需求？

解决方案 #

方法1：基于现有输出手动推导计算 #

从输出的对比逻辑可以直接推导non_working与part_time的差值：
mean(Not_Working) - mean(Part_Time) = [mean(Not_Working) - mean(Full_Time)] - [mean(Part_Time) - mean(Full_Time)]

比如针对more_than_two_children组，计算就是：-0.33049722 - (-0.22848561) = -0.10201161

你可以用dplyr直接在现有结果数据框上批量计算：

library(dplyr)

# 计算non_working与part_time的对比
avg_non_part <- avg %>%
  group_by(current_children) %>%
  summarise(
    term = "employment_level",
    contrast = "mean(Not_Working) - mean(Part_Time)",
    estimate = estimate[contrast == "mean(Not_Working) - mean(Full_Time)"] - estimate[contrast == "mean(Part_Time) - mean(Full_Time)"]
  ) %>%
  ungroup()

# 合并到原结果中
avg_full <- bind_rows(avg, avg_non_part)

方法2：直接用avg_slopes指定自定义对比 #

通过variables参数的自定义对比设置，直接生成你需要的三组转换效应，不需要依赖默认的基准组对比：

avg_custom <- avg_slopes(
  mod_2,
  variables = list(employment_level = c("Not_Working - Part_Time", "Not_Working - Full_Time", "Part_Time - Full_Time")),
  by = "current_children"
)

注意：如果你的原始数据中employment_level的水平名称是non_working（而非输出中的Not_Working），要把代码里的名称替换成实际的变量水平值。

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