使用fixest与marginaleffects做2SLS+FE回归时delta法置信区间出NA的问题

问题原因 #

当你的2SLS模型仅包含内生自变量x1、工具变量z和固定效应，且公式第一部分用y ~ 0（无外生变量/常数项）时，marginaleffects的avg_comparisons函数在计算delta法置信区间时，会因为模型参数矩阵的特殊结构，无法正确推导边际效应的标准误，最终返回NA。而添加额外外生变量x2后，模型结构完整，函数可以正常运行。

解决方案 #

方案1：给模型添加常数项（最简方法） #

修改feols的公式，将第一部分的y ~ 0改为y ~ 1（加入常数项），这样模型会包含一个外生的截距项，让avg_comparisons可以正确识别参数并计算标准误，同时不会改变x1的系数估计结果。

代码示例：

library(fixest)
library(marginaleffects)

# 生成可复现数据
i <- c("a","a","a","a","b","b","b","b","c","c","c","c","d","d","d","d")
t <- c(1,2,3,4,1,2,3,4,1,2,3,4,1,2,3,4)
y <- rnorm(16)
x1 <- rnorm(16)
x2 <- rnorm(16)
z <- rnorm(16)
n <- c(2,2,2,2,1,1,1,1,6,6,6,6,8,8,8,8)
df <- data.frame(i,t,y,x1,x2,z,n)

# 加入常数项的2SLS模型
TwoSLS <- feols(y ~ 1 | i + t | x1 ~ z, df)

# 调用avg_comparisons
avg_comparisons(TwoSLS)

运行后会返回包含正确标准误和置信区间的结果。

方案2：手动计算平均边际效应与置信区间 #

对于线性2SLS模型，x1的平均边际效应等于其回归系数（因为预测值随x1线性变化，单位变化的效应就是系数本身）。因此可以直接提取模型系数和方差-协方差矩阵，手动计算置信区间，结果与avg_comparisons完全等价。

代码示例：

# 原模型（无常数项）
TwoSLS <- feols(y ~ 0 | i + t | x1 ~ z, df)

# 提取x1的系数和标准误
coef_x1 <- coef(TwoSLS)["x1"]
se_x1 <- sqrt(vcov(TwoSLS)["x1", "x1"])

# 计算95%置信区间
ci_low <- coef_x1 - 1.96 * se_x1
ci_high <- coef_x1 + 1.96 * se_x1

# 整理结果
result <- data.frame(
  Term = "x1",
  Contrast = "mean(+1)",
  Estimate = round(coef_x1, 3),
  Std.Error = round(se_x1, 3),
  z = round(coef_x1 / se_x1, 3),
  Pr_z = round(2 * pnorm(-abs(coef_x1 / se_x1)), 3),
  "2.5 %" = round(ci_low, 3),
  "97.5 %" = round(ci_high, 3)
)

print(result)

补充说明 #

添加x2后函数正常运行的原因是，x2作为外生自变量存在，让模型的参数矩阵结构更完整，marginaleffects可以正常处理delta法所需的梯度计算和方差传播，因此不会出现NA值。

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