如何在R的lmer线性混合模型中计算Complier Average Causal Effect？

针对你有个体内嵌套观测的实验数据，lme4的lmer确实不支持直接纳入工具变量，以下是几种可行的R实现方案：

方案1：带聚类标准误的两阶段最小二乘(2SLS) #

使用estimatr包的iv_robust函数，支持个体层面的聚类标准误，适配嵌套数据结构，同时完成IV分析：

library(estimatr)

# 变量定义：Z=随机分配的干预（工具变量），D=complier变量（实际接受干预/操纵检查通过）
iv_model <- iv_robust(
  Y ~ D * time + factor(id) | Z * time + factor(id),
  data = df,
  clusters = id,
  se_type = "stata"
)

summary(iv_model)

• 第一阶段会用Z（分配的干预）预测D（complier状态），同时控制时间和个体固定效应
• 输出中D:time的系数即为CACE，代表complier群体在干预前后的平均因果效应
• 可通过第一阶段回归的F值检验工具变量的相关性（F值>10通常认为有效）

方案2：多层结构的广义矩估计(GMM) #

若需保留原模型的个体随机截距，可使用gmm包拟合多层IV模型，贴合你原lmer的混合效应设定：

library(gmm)
library(lme4)

# 定义矩条件函数，适配Y = β0 + β1D + β2time + β3D:time + u_id + e_ij的结构
moments_fun <- function(theta, data) {
  # 提取参数：固定效应系数、随机截距方差、残差方差
  beta <- theta[1:4]
  sigma_u <- theta[5]
  sigma_e <- theta[6]
  
  # 计算模型预测值与残差
  y_pred <- beta[1] + beta[2]*data$D + beta[3]*data$time + beta[4]*data$D*data$time
  resid_total <- data$Y - y_pred
  
  # 矩条件：工具变量Z与总残差不相关，同时控制时间交互
  m1 <- data$Z * resid_total
  m2 <- (data$Z * data$time) * resid_total
  # 随机截距的矩条件
  m3 <- tapply(resid_total, data$id, mean)^2 - sigma_u
  
  cbind(m1, m2, rep(m3, table(data$id)))
}

# 用原lmer模型的结果作为GMM初始参数
init_lmer <- lmer(Y ~ treatment*time + (1 | id), data = df)
init_params <- c(fixef(init_lmer), as.data.frame(VarCorr(init_lmer))$vcov, sigma(init_lmer)^2)

# 拟合GMM模型
gmm_fit <- gmm(moments_fun, data = df, t0 = init_params)
summary(gmm_fit)

• 该模型明确纳入了个体随机截距，和你原分析的结构一致
• 输出中对应D*time的系数就是CACE

方案3：固定效应IV模型 #

如果你的面板数据属于短面板（每个个体观测数较少），可使用固定效应控制个体异质性，结合IV分析：

library(AER)
library(sandwich)
library(lmtest)

# 拟合固定效应IV模型，工具变量为Z
fe_iv <- ivreg(Y ~ D*time | Z*time + factor(id), data = df)
# 加入个体层面的聚类标准误
coeftest(fe_iv, vcov = vcovCL, cluster = ~id)

• factor(id)控制个体固定效应，消除不随时间变化的个体混淆因素
• 聚类标准误修正了个体内观测的相关性

关键注意事项 #

• 工具变量有效性：必须确保Z（分配的干预）与D（complier状态）显著相关（第一阶段F值），且Z与结果Y的误差项无关
• CACE的解释：仅适用于complier群体（即操纵检查通过、实际接受干预的个体），这正是你需要的目标群体效应

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