如何修改R代码复现Stata面板数据Bootstrap抽样回归

我需要把Stata中针对面板数据的Bootstrap抽样回归逻辑复现到R语言中。核心要求是从个体标识idcode而非idcode-year层面进行Bootstrap抽样，之后基于抽中的样本估计回归参数，并将结果保存为Stata数据文件。Stata代码里的cluster(idcode)和idcluster(newid)选项正是用来实现按个体而非年度观测抽样的逻辑。

已有@PBulls提供的抽样部分R代码，我需要修改这段代码，完整复现下方Stata代码的全部功能。

现有R抽样代码 #

data <- haven::read_dta("http://www.stata-press.com/data/r9/nlswork.dta") |>
  dplyr::filter(!is.na(ttl_exp) & !is.na(hours))

panelboot <- function(df) {
  ids <- unique(df[["idcode"]])
  
  data.frame(
    idcode = sample(ids, replace=TRUE),
    nidcode = seq_along(ids)
  ) |> dplyr::left_join(df, by = "idcode", relationship = "many-to-many")
}

set.seed(1)
boot1 <- panelboot(data)

## 示例：ID 1130被选中3次，每次贡献8条观测
fid <- which(boot1[,1] == boot1[1,1])

## 重复抽样的ID 1130对应新ID：1, 933, 1671
boot1[fid,2]

Stata原代码及输出 #

Stata代码 #

use "nslwork.dta", clear 
* 重置时间序列设置，确保按个体抽样而非年度观测
tsset, clear
generate long newid = idcode // newid为长格式变量，用于标识抽样后的个体
tsset newid year

capture program drop savemargins
program savemargins, rclass
reg ttl_exp hours
end 
bootstrap _b,  saving(boot_output.dta, replace ) reps(10) cluster(idcode) idcluster(newid) : savemargins 

Stata输出结果 #

Bootstrap replications (10): .........10 done

Linear regression                                       Number of obs = 28,467
                                                        Replications  =     10
                                                        Wald chi2(1)  = 194.08
                                                        Prob > chi2   = 0.0000
                                                        R-squared     = 0.0118
                                                        Adj R-squared = 0.0118
                                                        Root MSE      = 4.6267

                              (Replications based on 4,710 clusters in idcode)
------------------------------------------------------------------------------
             |   Observed   Bootstrap                         Normal-based
     ttl_exp | coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       hours |   .0512427   .0036782    13.93   0.000     .0440335    .0584519
       _cons |   4.343879   .1302096    33.36   0.000     4.088673    4.599085
------------------------------------------------------------------------------

修改后的完整R代码 #

以下代码实现了和Stata完全一致的功能：按个体抽样、重复回归、计算Bootstrap统计量、保存结果到Stata文件。

library(dplyr)
library(haven)
library(broom)

# 加载并清理数据
data <- haven::read_dta("http://www.stata-press.com/data/r9/nlswork.dta") |>
  dplyr::filter(!is.na(ttl_exp) & !is.na(hours))

# 定义单次Bootstrap抽样+回归的函数
panelboot_reg <- function(df) {
  # 获取所有唯一个体ID
  ids <- unique(df[["idcode"]])
  # 有放回抽样个体ID
  sampled_ids <- sample(ids, replace = TRUE)
  # 拼接抽样后的面板数据集
  boot_data <- data.frame(idcode = sampled_ids) |>
    left_join(df, by = "idcode", relationship = "many-to-many")
  # 拟合线性回归并提取系数
  model <- lm(ttl_exp ~ hours, data = boot_data)
  tidy(model) |> pull(estimate)
}

# 设置随机种子，指定Bootstrap重复次数（对应Stata的reps(10)）
set.seed(1)
reps <- 10
# 执行重复抽样回归
boot_results <- replicate(reps, panelboot_reg(data))

# 整理Bootstrap系数结果
boot_coefs <- t(boot_results) |>
  as.data.frame() |>
  rename(hours = V1, `_cons` = V2)

# 计算观测模型的系数及Bootstrap统计量
observed_model <- lm(ttl_exp ~ hours, data = data)
observed_coefs <- tidy(observed_model) |> pull(estimate)

boot_summary <- data.frame(
  variable = c("hours", "_cons"),
  observed = observed_coefs,
  se_bootstrap = apply(boot_coefs, 2, sd),
  z_score = observed_coefs / apply(boot_coefs, 2, sd),
  p_value = 2 * pnorm(-abs(observed_coefs / apply(boot_coefs, 2, sd))),
  ci_lower = observed_coefs - 1.96 * apply(boot_coefs, 2, sd),
  ci_upper = observed_coefs + 1.96 * apply(boot_coefs, 2, sd)
)

# 打印类似Stata格式的结果
cat("Linear regression                                       Number of obs =", nrow(data), "\n")
cat("                                                        Replications  =", reps, "\n\n")
print(boot_summary, row.names = FALSE)

# 将Bootstrap系数保存为Stata数据文件（对应Stata的saving(boot_output.dta)）
haven::write_dta(boot_coefs, "boot_output.dta")

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