带诊断的固定效应模型IV回归实现问题求助

一、ivreg：确认固定效应属性与诊断补充 #

你用ivreg写的代码通过加入City和Year虚拟变量实现固定效应，这属于虚拟变量法固定效应，和plm的within模型（去组内均值）在平衡面板中是等价的，但估计形式不同。

• 确认固定效应：可以对比不加City/Year的OLS结果，看固定效应是否显著；或者用plm先对数据做双向within变换，再用ivreg回归。
• 补充诊断：summary(fit_ivreg, diagnostics = T)会输出弱工具变量检验（Cragg-Donald F统计量），但Wu-Hausman检验需要手动对比固定效应OLS和IV的系数差异，用hausman()函数（需先拟合固定效应OLS模型）。

修正后的ivreg固定效应IV写法（更清晰）：

# 拟合固定效应OLS用于对比
fit_fe_ols = lm(y ~ x1 + c1 + c2 + c3 + factor(City) + factor(Year), data = data)
# IV回归（虚拟变量法固定效应）
fit_ivreg = ivreg(y ~ x1 + c1 + c2 + c3 + factor(City) + factor(Year) | inst + c1 + c2 + c3 + factor(City) + factor(Year), data = data)
# 查看诊断结果
summary(fit_ivreg, diagnostics = T)
# 执行Wu-Hausman检验
hausman(fit_ivreg, fit_fe_ols)

二、plm：手动补充诊断统计量 #

plm的双向固定效应IV模型（model="within", effect="twoways"）默认summary不输出诊断，需手动计算：

1. 弱工具变量检验：提取第一阶段回归的F统计量

# 提取第一阶段回归结果
first_stage = plm(x1 ~ inst + c1 + c2 + c3, data = data, effect = "twoways", model = "within", index = c("City", "Year"))
summary(first_stage) # 一般F统计量>10说明工具变量不弱

2. Wu-Hausman检验：对比双向固定效应OLS和IV的系数

# 拟合双向固定效应OLS
fit_plm_ols = plm(y ~ x1 + c1 + c2 + c3, data = data, effect = "twoways", model = "within", index = c("City", "Year"))
# 拟合IV模型（你的原代码）
fit_plm = plm(y ~ x1 + c1 + c2 + c3 | inst + c1 + c2 + c3, effect = "twoways", model = "within", index = c("City", "Year"), inst.method = 'baltagi')
# 执行Wu-Hausman检验
hausman(fit_plm, fit_plm_ols)

3. 注意：plm的inst.method='baltagi'仅适用于平衡面板；如果是非平衡面板，建议改用inst.method='iv'。

三、fixest：结果差异原因与正确设定 #

你的fixest代码是虚拟变量法双向固定效应IV，而plm的model="within", effect="twoways"是去均值法双向固定效应，两者在不平衡面板中系数会有差异，这是核心原因：

• 虚拟变量法：直接控制每个City和Year的虚拟变量，保留所有有效观测
• 去均值法：对每个City去组内均值、每个Year去时间均值，不平衡面板中会改变变量中心化方式，导致系数不同

如果要让fixest和plm结果一致，需让fixest使用去均值法：

# 用demean参数实现去均值法双向固定效应IV
fit_fixest_demean = feols(y ~ c1 + c2 + c3 | City + Year, data = data, endog = ~x1, inst = ~inst, demean = "twoways")
summary(fit_fixest_demean)

fixest默认会输出完整诊断：包括第一阶段F统计量（弱工具变量检验）、Wu-Hausman检验（Hausman test of endogeneity），无需额外代码。

另外，fixest默认聚类到第一个固定效应维度（City），如果要和plm的默认标准误一致，可以设置vcov = "iid"：

fit_fixest = feols(y ~ c1 + c2 + c3 | City + Year | x1 ~ inst, data = data, vcov = "iid")
summary(fit_fixest)

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