Julia运行面板固定效应模型后如何提取平均边际效应与预测值

在Julia环境下使用FixedEffectModels.jl开展高维固定效应模型估计时，完成回归后无法直接提取交互项对应的平均边际效应、不同场景下的分组预测基准值。
参考Stata实现逻辑，使用如下代码即可在固定效应模型后输出指定变量的边际效应：

xtreg chronic_illness age country_birth social_class#macro_unemployment, fe
margins crisis, dydx(social_class)

当前运行的Julia回归代码如下：

m = reg(df1, @formula(chronic_illness ~ status + age + social_class*crisis + fe(id) + fe(year) + fe(country)), contrasts=contr, save=true)

变量说明：

• 被解释变量chronic_illness为二元变量，0代表无慢性病
• 核心解释变量crisis为二元变量，0代表无金融危机
• 研究目标：对比金融危机发生前后，不同社会阶层群体的慢性病患病平均水平，当前模型仅输出交互项系数，无法直接得到分组基准预测值。

模型回归输出结果：

Fixed Effect Model                                           
========================================================================================================
Number of obs:                              1468882  Degrees of freedom:                          459252
R2:                                           0.703  R2 Adjusted:                                  0.567
F-Stat:                                     62.8378  p-value:                                      0.000
R2 within:                                    0.001  Iterations:                                      18
========================================================================================================
cillness                           |    Estimate  Std.Error     t value Pr(>|t|)   Lower 95%   Upper 95%
-------------------------------------------------------------------------------------------------------
status: Unemployed                 |   0.0145335 0.00157535     9.22556    0.000   0.0114459   0.0176212
status: missing                    | -0.00702545  0.0136504    -0.51467    0.607  -0.0337797   0.0197288
age                                |  0.00178437    2.79058 0.000639427    0.999    -5.46766     5.47123
class: Lower-middle class          |  0.00458331    250.251  1.83149e-5    1.000    -490.478     490.487
class: Working class               |   0.0286466    163.324 0.000175398    1.000    -320.081     320.138
crisis                             | -0.00600744 0.00156138    -3.84753    0.000 -0.00906768 -0.00294719
class: Lower-middle class & crisis | -0.00189866 0.00192896   -0.984289    0.325 -0.00567936  0.00188205
class: Working class & crisis      | -0.00332881 0.00170221    -1.95558    0.051  -0.0066651  7.46994e-6

你跑的是线性概率模型，两种方法都可以实现和Stata margins一致的输出结果，不需要更换模型框架。

方法1：手动计算（无额外依赖） #

因为你在reg时已经加了save=true参数，模型对象存储了所有估计分量，直接调用内置的predict函数生成反事实预测值即可，逻辑和Stata margins“其他变量取观测原值、核心变量取指定值”的计算规则完全一致：

1. 构造两个反事实数据集，分别对应无危机、有危机场景，除了crisis变量统一赋值外，其余所有变量（包括固定效应对应的id、year、country列）保持和原数据完全一致：

   df_crisis0 = copy(df1)
   df_crisis0.crisis .= 0
   df_crisis1 = copy(df1)
   df_crisis1.crisis .= 1
   

2. 分别调用predict生成两个场景下的预测值，函数会自动处理交互项、固定效应的贡献，不需要手动拼接系数：

   pred_nocrisis = predict(m, df_crisis0)
   pred_crisis = predict(m, df_crisis1)
   

3. 按社会阶层分组求预测值的均值，就是你需要的不同场景下各阶层慢性病患病平均基准值；同阶层两个场景的均值差，就是该阶层对应的金融危机平均边际效应：

   using DataFrames, Statistics
   # 合并预测值到原数据表
   df1[!, :pred_nocrisis] = pred_nocrisis
   df1[!, :pred_crisis] = pred_crisis
   # 分组计算统计量
   res = combine(groupby(df1, :social_class),
       :pred_nocrisis => mean => :avg_rate_nocrisis,
       :pred_crisis => mean => :avg_rate_crisis,
       [:pred_nocrisis, :pred_crisis] => ((x,y) -> mean(y.-x)) => :me_crisis
   )
   

如果需要边际效应的标准误，直接对分组差值做bootstrap抽样即可，大样本下运行速度很快。

方法2：调用边际效应专用包 #

如果需要自动输出带标准误、置信区间的边际效应结果，直接用MarginalEffects.jl即可，该包原生兼容FixedEffectModels.jl的回归对象，语法和Stata margins高度对齐：

using MarginalEffects
# 输出不同社会阶层下，crisis的平均边际效应及标准误
me_res = marginal_effects(m, :crisis, across=:social_class)
# 输出不同crisis取值下，各社会阶层的调整预测均值（即你要的基准值）
pred_res = estimated_margins(m, :crisis, across=:social_class)

补充说明 #

你当前回归结果里age、社会阶层主效应的标准误异常偏大，是因为同时控制了个体、年份、国家三层高维固定效应，这类几乎不随时间变化的变量变异被固定效应完全吸收，属于正常现象，不是代码错误，不影响交互项、边际效应和预测值的计算。
如果需要预测值严格落在[0,1]的概率区间，可以换用包内置的nlreg函数估计固定效应logit模型，后续预测、边际效应计算的逻辑和上述流程完全一致。

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