marginaleffects固定效应OLS回归分组预测取值来源与计算方法

我希望使用marginaleffects包计算固定效应OLS回归的分组预测。当我为预测生成新数据但未指定要使用的固定效应变量具体水平时，预测函数会基于am==0、gear==3和carb==4进行计算（见下方代码）。请问这些取值来自何处？在该示例中，计算固定效应下平均预测的最合理方式是什么？

library(tidyverse)
library(fixest)
library(marginaleffects)

dat <- mtcars

m1 <- dat %>% 
  mutate(cyl = as_factor(cyl)) %>% 
  feols(mpg ~ hp*cyl | am + gear + carb)

etable(m1)
#>                               m1
#> Dependent Var.:              mpg
#>                                 
#> hp              -0.1316 (0.0504)
#> cyl6              -47.87 (27.42)
#> cyl8              -21.31 (14.80)
#> hp x cyl6        0.4434 (0.2619)
#> hp x cyl8        0.1471 (0.0665)
#> Fixed-Effects:  ----------------
#> am                           Yes
#> gear                         Yes
#> carb                         Yes
#> _______________ ________________
#> S.E.: Clustered           by: am
#> Observations                  32
#> R2                       0.88892
#> Within R2                0.38565
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

pred1 <- predictions(m1, 
                     newdata = datagrid(
                       hp = c(min(dat$hp):max(dat$hp)), 
                       cyl = levels(dat$cyl))) 
pred1
#> 
#>  Estimate Std. Error     z Pr(>|z|)   S  2.5 % 97.5 % cyl am gear carb  hp
#>      9.76       8.72 1.119    0.263 1.9  -7.33   26.9   8  0    3    4  52
#>      9.78       8.60 1.136    0.256 2.0  -7.09   26.6   8  0    3    4  53
#>      9.79       8.49 1.154    0.249 2.0  -6.84   26.4   8  0    3    4  54
#>      9.81       8.37 1.172    0.241 2.1  -6.60   26.2   8  0    3    4  55
#>      9.82       8.25 1.190    0.234 2.1  -6.35   26.0   8  0    3    4  56
#> --- 274 rows omitted. See ?avg_predictions and ?print.marginaleffects --- 
#>     14.09      23.89 0.590    0.555 0.8 -32.73   60.9   8  0    3    4 331
#>     14.11      24.01 0.588    0.557 0.8 -32.94   61.2   8  0    3    4 332
#>     14.12      24.12 0.586    0.558 0.8 -33.15   61.4   8  0    3    4 333
#>     14.14      24.24 0.583    0.560 0.8 -33.37   61.6   8  0    3    4 334
#>     14.16      24.36 0.581    0.561 0.8 -33.58   61.9   8  0    3    4 335
#> Columns: rowid, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, mpg, cyl, am, gear, carb, hp

1. 固定效应默认取值的来源 #

当使用datagrid()生成新数据且未指定固定效应变量的水平时，marginaleffects会自动为未指定变量选择数据中最常见的水平（众数）：

• am的最常见水平是0（mtcars中19个观测为0）
• gear的最常见水平是3（15个观测为3）
• carb的最常见水平是4（10个观测为4）

背后的逻辑是datagrid()调用typical_value()函数生成“典型值”：分类/二分类变量取出现频率最高的水平，数值变量取中位数，以此作为未指定变量的基准取值。

2. 固定效应下平均预测的合理方式 #

根据研究目的，有两种常用的合理方式：

方式一：基于原始数据分布的平均预测 #

使用avg_predictions()函数，直接对原始数据中所有观测的预测值取平均，或按指定分组（如cyl）计算边际平均预测：

# 按cyl分组计算平均预测，自动整合所有固定效应的样本分布
avg_pred <- avg_predictions(m1, by = "cyl")
print(avg_pred)

这种方式得到的结果贴近真实数据的分布，反映协变量在实际样本中的平均预期结果。

方式二：覆盖所有固定效应水平的平均预测 #

如果想忽略原始样本分布，对固定效应的每个水平都计算预测后取平均，可以用datagrid()生成包含所有固定效应水平的新数据，再计算并平均预测值：

# 生成包含cyl、am、gear、carb所有水平的新数据，hp固定为中位数
new_dat <- datagrid(
  model = m1,
  hp = median(dat$hp),
  cyl = levels(dat$cyl),
  am = unique(dat$am),
  gear = unique(dat$gear),
  carb = unique(dat$carb)
)

# 计算预测并按cyl分组平均
pred_avg <- predictions(m1, newdata = new_dat) %>%
  group_by(cyl) %>%
  summarize(avg_estimate = mean(estimate))
print(pred_avg)

这种方式适合展示协变量在所有固定效应水平下的平均预期，不受原始样本分布的影响。

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