R语言高效生成0/1向量排列及计算随机试验精确p值优化方法

前置说明 #

445样本中取185个处理组的总排列组合量级为$C_{445}^{185} \approx 10^{131}$，属于完全无法生成和计算的规模，因此随机采样置换是唯一可行的精确p值计算方案，无需追求生成完整排列集合。

核心优化逻辑 #

原实现性能瓶颈主要来自冗余数据框合并、循环内重复子集索引与均值计算，优化核心是用向量化运算完全替代循环：

• 直接提取结果变量re78为向量，无需和置换矩阵合并为数据框
• 利用矩阵乘法一次性计算所有置换下的处理组结果总和，避免循环内重复计算
• 通过数学公式转换批量生成所有置换的效应量，完全消除显式循环

优化后完整代码 #

# 读入数据与基础参数计算
lalonde <- read.csv("lalonde.csv")
re78 <- lalonde$re78
n_total <- length(re78)
n1 <- sum(lalonde$treat)
n0 <- n_total - n1
sum_re78 <- sum(re78)
# 计算原始观测的效应量
tau_hat <- abs(mean(re78[lalonde$treat == 1]) - mean(re78[lalonde$treat == 0]))

# 生成10万次置换矩阵（每列对应1次置换的处理分组）
set.seed(123) # 固定种子保证结果可复现
perms <- replicate(100000, sample(c(rep(1, n1), rep(0, n0))))

# 向量化计算所有置换的效应量
sum_treat <- crossprod(perms, re78)[, 1] # 批量得到每个置换的处理组re78总和
diff_vec <- abs(sum_treat / n1 - (sum_re78 - sum_treat) / n0)

# 计算精确p值
p_val <- mean(diff_vec >= tau_hat)

额外性能提升可选方案 #

如果需要采样超过100万次置换，可进一步做如下优化：

• 用parallel::parReplicate并行生成置换和计算效应量
• 用Rcpp实现核心计算逻辑，性能可再提升1-2个数量级

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