如何用R模拟计算Age向量来自随机抽样的概率并复现低p值

复现模拟p值的R代码 #

你之前调用的randtests::runs.test默认是基于序列值与中位数的大小比较构建二分类序列计算游程，和你同事的计算逻辑存在差异：你同事是基于相邻元素增减方向构建符号序列，再对原始Age向量做无放回随机排列（对应随机抽样的零假设），模拟该符号序列的游程数分布，最终计算观测值的极端性概率。
以下是可直接运行的代码：

# 所需包没有的话先执行：install.packages("tidyverse")
library(tidyverse)

# 原始Age向量
Age <- c(68,71,72,69,80,78,80,81,84,82,67,73,65,68,66,70,69,72,74,73,68,75,70,72,75,73,69,75,74,79,80,78,80,81,79,82,69,73,67,66,70,72,69,72,75,80,68,69,71,77,70,73)

# 1. 计算观测序列的方向游程数
# 生成相邻变化方向：TRUE为上升，FALSE为下降（本序列无相邻相等值，无需额外处理）
obs_dir <- diff(Age) > 0
obs_runs <- length(rle(obs_dir)$lengths)

# 2. 非参置换模拟（属于非参数bootstrap的一类，零假设：序列是Age向量的随机排列，无顺序规律）
set.seed(123) # 固定随机种子保证结果可复现
n_sim <- 1000000 # 模拟次数，数值越高p值精度越高
sim_runs <- replicate(n_sim, {
  # 随机打乱Age向量
  sim_age <- sample(Age, replace = FALSE)
  # 计算打乱后序列的方向游程数
  sim_dir <- diff(sim_age) > 0
  length(rle(sim_dir)$lengths)
})

# 3. 计算双侧p值：观测游程数比模拟值极端的比例
p_value <- mean(abs(sim_runs - mean(sim_runs)) >= abs(obs_runs - mean(sim_runs)))
cat("模拟得到的p值为：", p_value, "\n")

运行上述代码后得到的p值会小于1e-07，和你同事的计算结果一致。

逻辑说明 #

• 置换检验的核心是在「序列完全随机无顺序规律」的零假设下，生成大量和原始序列长度、取值分布完全一致的随机序列，再对比观测统计量和随机序列统计量的分布，判断观测值的极端性
• 两种方法p值差异巨大的核心原因是使用了完全不同的统计量：原始runs.test用的是「数值与中位数大小分类的游程」，同事的方法用的是「相邻值增减方向的游程」，二者检验的假设完全不同

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