R中带州级聚类标准误的Logistic回归工具变量分析方法问询

我有个体层面数据，用于分析州级教育支出对学生个体成绩的影响。学生成绩是二元变量（0=未通过测试，1=通过测试），已运行带州层面聚类标准误的Logistic回归：

library(miceadds)
df_logit <- data.frame(performance = c(0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0),
                       state = c("MA", "MA", "MB", "MC", "MB", "MD", "MA", "MC", "MB", "MD", "MB", "MC", "MA", "MA", "MA", "MA", "MD", "MA","MB","MA","MA","MD","MC","MA","MA","MC","MB","MB","MD", "MB"),
                       expenditure = c(123000, 123000,654000, 785000, 654000, 468000, 123000,  785000, 654000, 468000, 654000, 785000,123000,123000,123000,123000, 468000,123000, 654000, 123000, 123000, 468000,785000,123000, 123000, 785000, 654000, 654000, 468000,654000),
                       population = c(0.25, 0.25, 0.12, 0.45, 0.12, 0.31, 0.25, 0.45, 0.12, 0.31, 0.12, 0.45, 0.25, 0.25, 0.25, 0.25, 0.31, 0.25, 0.12, 0.25, 0.25, 0.31, 0.45, 0.25, 0.25, 0.45, 0.12, 0.12, 0.31, 0.1),
                       left_wing = c(0.10, 0.10, 0.12, 0.18, 0.12, 0.36, 0.10, 0.18, 0.12, 0.36, 0.12, 0.18, 0.10, 0.10, 0.10, 0.10, 0.36, 0.10, 0.12, 0.10, 0.10, 0.36, 0.18, 0.10, 0.10,0.18, 0.12, 0.12, 0.36, 0.12))

df_logit$performance <- as.factor(df_logit$performance)

glm_clust_1 <- miceadds::glm.cluster(data=df_logit, formula=performance ~ expenditure + population,
                                     cluster="state", family=binomial(link = "logit")) 
summary(glm_clust_1) 

因无法排除支出的内生性，希望用州级左翼政党占比left_wing作为教育支出的工具变量，但未找到能在带州层面聚类标准误的Logistic回归中运行含控制变量的工具变量法的命令，同时希望输出Wu-Hausman检验、弱工具变量检验等常用诊断（类似OLS中ivreg的功能）。理想是将以下OLS-IV命令适配到二元因变量并实现州层面聚类标准误：

iv_1 <- ivreg(performance ~ population + expenditure | left_wing + population, data=df_logit)
summary(iv_1, cluster="state", diagnostics = TRUE)

针对二元因变量的内生性问题，**两阶段残差纳入法（TSRI）**是可靠的一致估计方法，优于直接对Logit模型用线性2SLS（后者系数估计不一致）。以下是完整实现，包含聚类标准误和所需诊断检验：

具体代码实现 #

# 加载所需包
library(lmtest)
library(sandwich)
library(miceadds)

# ----------------------
# 第一阶段：内生变量回归（用于提取残差+弱工具检验）
# ----------------------
first_stage <- lm(expenditure ~ population + left_wing, data = df_logit)
# 计算聚类调整的弱工具变量F统计量
first_stage_clust <- coeftest(first_stage, vcov = vcovCL, cluster = ~state)
f_stat <- waldtest(first_stage, test = "F", vcov = vcovCL, cluster = ~state)$F[2]

# 保存第一阶段残差，用于第二阶段回归
df_logit$resid_exp <- residuals(first_stage)

# ----------------------
# 第二阶段：Logistic回归（含残差）+ 聚类标准误
# ----------------------
tsri_model <- miceadds::glm.cluster(
  data = df_logit,
  formula = performance ~ expenditure + population + resid_exp,
  cluster = "state",
  family = binomial(link = "logit")
)

# 对比模型：不含残差的原Logistic回归
original_model <- miceadds::glm.cluster(
  data = df_logit,
  formula = performance ~ expenditure + population,
  cluster = "state",
  family = binomial(link = "logit")
)

# ----------------------
# 诊断检验输出
# ----------------------
# 1. 弱工具变量检验（聚类调整F统计量）
cat("=== 弱工具变量检验 ===\n")
cat("聚类调整F统计量：", round(f_stat, 3), "\n")
cat("注：F统计量>10通常认为工具变量不存在弱识别问题\n\n")

# 2. Wu-Hausman检验（检验内生性）
coef_original <- coef(original_model)
vcov_original <- vcov(original_model)
coef_tsri <- coef(tsri_model)[names(coef_original)] # 匹配变量
vcov_tsri <- vcov(tsri_model)[names(coef_original), names(coef_original)]

# 构造Hausman检验统计量
hausman_stat <- t(coef_tsri - coef_original) %*% solve(vcov_tsri - vcov_original) %*% (coef_tsri - coef_original)
p_value <- pchisq(hausman_stat, df = length(coef_tsri))

cat("=== Wu-Hausman内生性检验 ===\n")
cat("卡方统计量：", round(hausman_stat, 3), "\n")
cat("P值：", round(p_value, 3), "\n")
cat("注：P值<0.05则拒绝原假设，认为expenditure存在内生性\n\n")

# 输出第二阶段TSRI模型结果
cat("=== 第二阶段TSRI Logistic回归结果（聚类标准误） ===\n")
summary(tsri_model)

补充说明 #

1. TSRI方法优势：专门针对非线性模型（如Logit）的内生性问题，估计量一致，且易结合聚类标准误。
2. 替代包实现：若想使用专门的IV Logit工具，可尝试ivregress包的ivlogit函数，但需手动计算聚类标准误：

library(ivregress)
iv_logit_model <- ivlogit(performance ~ expenditure + population | left_wing + population, data = df_logit)
# 手动添加聚类标准误
vcov_clust <- vcovCL(iv_logit_model, cluster = ~state)
coeftest(iv_logit_model, vcov = vcov_clust)

不过该方法的诊断检验需额外手动计算，TSRI方法更直观且易实现所有所需检验。

内容的提问来源于stack exchange，提问作者R-User