如何在Stata或R中实现Bootstrapping计算Z分数均值及置信上下限

给定如下X和Y的数据集：

X   Y
1   2
3   9
7   16
10  17
14  23

需要通过Bootstrapping有放回重采样Y共250次，生成类似如下的结果结构：

X    Y    y1   y2  z_y z_y1 z_y2  mean_z lower_bound  upper_bound 
1    2    2    9    .    .    .    .        .              . 
3    9    16   23
7   16    16   17
10  17    23   2
14  23    9    16

具体需求：

• 对每个X值，计算原始Y、每次Bootstrap重采样得到的y（如y1、y2...y250）的Z分数
• 计算每个X对应的所有Z分数（原始+250次重采样）的均值mean_z
• 计算该均值的置信区间上下限lower_bound和upper_bound

Stata实现步骤 #

1. 导入原始数据 #

clear
input X Y
1   2
3   9
7   16
10  17
14  23
end

2. 生成250次Bootstrap重采样Y列 #

循环生成250个有放回采样的Y变量：

forvalues i=1/250 {
    gen y`i' = Y[runiformint(1, _N)]
}

3. 计算所有Y相关变量的Z分数 #

先获取原始Y的整体均值和标准差，再批量计算Z分数：

sum Y, meanonly
local mu = r(mean)
local sd = r(sd)

gen z_y = (Y - `mu')/`sd'
forvalues i=1/250 {
    gen z_y`i' = (y`i' - `mu')/`sd'
}

4. 计算Z分数均值及置信区间 #

• 先计算每个X对应的所有Z分数均值：

egen mean_z = rowmean(z_y z_y1-z_y250)

• 再通过Bootstrap法计算均值的95%置信区间（可调整reps参数改变重采样次数）：

bootstrap lower_bound=r(lb) upper_bound=r(ub), reps(1000): by X: ci mean mean_z

运行后会输出每个X对应的置信区间，可手动合并到原数据集中。

R实现步骤 #

1. 构造原始数据框 #

df <- data.frame(
  X = c(1, 3, 7, 10, 14),
  Y = c(2, 9, 16, 17, 23)
)

2. 生成250次Bootstrap重采样Y列 #

设置随机种子保证结果可重复，批量生成重采样列：

set.seed(123)
bootstrap_y <- replicate(250, sample(df$Y, size = nrow(df), replace = TRUE))
colnames(bootstrap_y) <- paste0("y", 1:250)
df <- cbind(df, bootstrap_y)

3. 计算所有Y相关变量的Z分数 #

基于原始Y的均值和标准差，批量计算Z分数：

mu <- mean(df$Y)
sd_val <- sd(df$Y)

# 原始Y的Z分数
df$z_y <- (df$Y - mu)/sd_val

# 所有重采样Y的Z分数
z_col_names <- paste0("z_y", 1:250)
df[z_col_names] <- lapply(df[paste0("y", 1:250)], function(x) (x - mu)/sd_val)

4. 计算Z分数均值及置信区间 #

• 先计算每个X的Z分数均值：

df$mean_z <- rowMeans(df[, c("z_y", z_col_names)])

• 自定义函数计算每个X对应mean_z的95%置信区间：

bootstrap_ci <- function(data, n_reps = 1000) {
  ci_results <- data.frame(X = unique(data$X), lower_bound = NA, upper_bound = NA)
  for (x_val in unique(data$X)) {
    x_subset <- data[data$X == x_val, ]
    boot_samples <- replicate(n_reps, {
      sample_z <- sample(c(x_subset$z_y, unlist(x_subset[z_col_names])), size = 251, replace = TRUE)
      mean(sample_z)
    })
    ci_results[ci_results$X == x_val, ] <- c(x_val, quantile(boot_samples, c(0.025, 0.975)))
  }
  return(ci_results)
}

# 合并置信区间到原数据框
ci_df <- bootstrap_ci(df)
df <- merge(df, ci_df, by = "X")

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