如何在R中为面板数据个体固定效应模型提取预测概率与平均边际效应

针对你提到的百万级样本下lm加个体虚拟变量内存不足的问题，以下三个方案都可以直接基于你已经跑通的model_plm或model_felm输出结果，完全不需要生成个体虚拟变量，内存占用远低于lm方案，和你常用的分析包兼容性也很高：

方案1：直接使用ggeffects包（和你原有操作习惯完全一致） #

最新版ggeffects已经原生支持plm和felm类模型，你只需要把已经拟合好的模型对象直接传入ggpredict即可，参数和你之前用lm时完全一致：

# 基于plm模型生成预测值
pred_plm <- ggpredict(model_plm, c("residence","union"), 
                      vcov.fun = "vcovCL", 
                      vcov.type = "HC1",
                      vcov.args = list(cluster = Males$nr))

# 基于felm模型生成预测值（更适配大样本场景，运行速度更快）
pred_felm <- ggpredict(model_felm, c("residence","union"), 
                       vcov.fun = "vcovCL", 
                       vcov.type = "HC1",
                       vcov.args = list(cluster = Males$nr))

输出结果和你之前用model_lm得到的pred_ggeffects完全一致。

方案2：使用emmeans包提取平均边际效应 #

如果你需要更灵活的边际效应对比，emmeans也支持plm和felm对象：

library(emmeans)

# 以felm模型为例，plm模型用法相同
# 生成residence和union交互项的各组拟合均值
emm_res <- emmeans(model_felm, ~ residence * union,
                   vcov. = vcovCL(model_felm, cluster = ~nr, type = "HC1"))

# 查看各组预测值
emm_res

# 生成组间边际效应对比
pairs(emm_res, by = "residence")

方案3：手动计算（极端大样本下性能最优） #

如果样本量特别大，需要极致运行速度，可以直接提取模型系数和协方差矩阵手动计算预测值，不需要依赖任何额外的边际效应包：

# 提取felm模型的系数和聚类稳健标准误对应的协方差矩阵
coef_vec <- coef(model_felm)
vcov_mat <- vcovCL(model_felm, cluster = ~nr, type = "HC1")

# 取控制变量的样本均值（exper、health控制在平均水平）
mean_exper <- mean(Males$exper)
mean_health <- mean(Males$health)

# 计算residence和union四个组合的预测值（假设两个变量都是二分类0/1，可根据你的变量实际取值调整）
pred_r0u0 <- coef_vec["(Intercept)"] + coef_vec["exper"]*mean_exper + coef_vec["health"]*mean_health
pred_r0u1 <- pred_r0u0 + coef_vec["union"]
pred_r1u0 <- pred_r0u0 + coef_vec["residence"]
pred_r1u1 <- pred_r0u0 + coef_vec["residence"] + coef_vec["union"] + coef_vec["residence:union"]

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