复制Chaisemartin论文时将Stata命令bsample, cluster(nr)转R的需求

一、Stata中bsample, cluster(nr)的作用 #

bsample, cluster(nr)是聚类放回bootstrap抽样命令，核心作用：

• 以nr作为聚类单位，从所有聚类组中随机放回抽取与原样本聚类数量一致的组；
• 每个被抽中的聚类组，会完整保留组内所有观测值；
• 这种抽样方式能保留原数据的组内相关性结构，专门用于聚类稳健的bootstrap推断，适配面板数据的特性。

二、转换后的R代码实现 #

首先安装依赖包：

install.packages(c("fixest", "didmultiplegt", "dplyr"))

核心代码（替换df为你的实际数据框名）：

library(fixest)
library(didmultiplegt)
library(dplyr)

# 设置随机种子，与Stata保持一致
set.seed(1)

# 初始化结果矩阵，存储200次循环的差值
results_mat <- matrix(NA, nrow = 200, ncol = 3)
colnames(results_mat) <- c("fe-fd", "fe-did_m", "fd-did_m")

# 获取所有聚类ID及数量
all_clusters <- unique(df$nr)
n_clust <- length(all_clusters)

# 循环执行200次bootstrap
for(i in 1:200){
  # 1. 聚类放回抽样：模拟Stata的bsample, cluster(nr)
  sampled_clusts <- sample(all_clusters, size = n_clust, replace = TRUE)
  # 构建bootstrap样本，处理重复抽中的聚类
  boot_data <- df %>%
    filter(nr %in% sampled_clusts) %>%
    group_by(nr) %>%
    # 统计每个聚类被抽中的次数，重复对应观测
    mutate(reps = sum(sampled_clusts == nr)) %>%
    uncount(reps) %>%
    ungroup()
  
  # 2. 双向固定效应回归（对应xtreg fe robust）
  fe_reg <- feols(lwage ~ d81 + d82 + d83 + d84 + d85 + d86 + d87 + union | nr,
                  data = boot_data, vcov = ~nr)
  beta_fe <- coef(fe_reg)["union"]
  
  # 3. 一阶差分回归（对应reg ... cluster(nr)）
  fd_reg <- feols(diff_lwage ~ diff_union + d82 + d83 + d84 + d85 + d86 + d87,
                  data = boot_data, vcov = ~nr)
  beta_fd <- coef(fd_reg)["diff_union"]
  
  # 4. did_multiplegt估计（对应Stata命令）
  did_m_reg <- did_multiplegt(y = boot_data$lwage,
                              g = boot_data$nr,
                              t = boot_data$year,
                              treat = boot_data$union)
  did_m <- did_m_reg$ATT[1] # 对应Stata的e(effect_0)
  
  # 5. 存储三个差值结果
  results_mat[i, ] <- c(beta_fe - beta_fd, beta_fe - did_m, beta_fd - did_m)
}

# 查看前5次循环的结果
head(results_mat)

代码说明 #

• fixest包的feols函数实现固定效应和聚类稳健标准误，对应Stata的xtreg fe robust和reg ... cluster(nr)；
• 抽样部分通过sample+uncount精确模拟Stata的bsample, cluster(nr)行为，确保聚类放回抽样的一致性；
• did_multiplegt包的输出中，$ATT[1]对应Stata返回的e(effect_0)。

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