R语言实现Kolmogorov-Smirnov随机排序检验及功效曲线绘制求助

常见错误点 #

• 未指定ks.test的单边检验方向：默认是双边检验，与你的单边备择假设不匹配
• 置换检验逻辑错误：合并样本后随机分配的方式不符合独立样本置换的要求
• 功效计算时未正确重复模拟：样本量小或重复次数不足会导致结果偏差
• 检验统计量与备择假设不匹配：单边检验应使用F_Y - F_X的最大值而非双边的绝对差

修正后的R代码 #

# 设置参数
set.seed(123) # 固定随机种子保证可复现
m <- 60
n <- 70
mu_x <- 5
sigma2 <- 1
sigma <- sqrt(sigma2)
mu_y_seq <- seq(4, 6, by = 0.02)
alpha <- 0.05
n_sim <- 1000 # 模拟次数，可根据需要调整

# 定义功效计算函数
compute_power <- function(mu_y) {
  reject_count <- 0
  for (i in 1:n_sim) {
    # 生成样本
    X <- rnorm(m, mean = mu_x, sd = sigma)
    Y <- rnorm(n, mean = mu_y, sd = sigma)
    
    # 单边K-S检验（使用ks.test的单边选项）
    ks_result <- ks.test(Y, X, alternative = "greater")
    if (ks_result$p.value <= alpha) {
      reject_count <- reject_count + 1
    }
    
    # 如果你想用置换检验手动实现K-S检验（而非直接用ks.test），可以用以下代码：
    # # 合并样本
    # combined <- c(X, Y)
    # # 原样本的单边K-S统计量：sup(F_Y - F_X)
    # ecdf_X <- ecdf(X)
    # ecdf_Y <- ecdf(Y)
    # all_points <- sort(c(X, Y))
    # stat_obs <- max(ecdf_Y(all_points) - ecdf_X(all_points))
    # 
    # # 置换检验计算p值
    # perm_stats <- replicate(999, {
    #   perm <- sample(combined)
    #   perm_X <- perm[1:m]
    #   perm_Y <- perm[(m+1):(m+n)]
    #   perm_ecdf_X <- ecdf(perm_X)
    #   perm_ecdf_Y <- ecdf(perm_Y)
    #   max(perm_ecdf_Y(all_points) - perm_ecdf_X(all_points))
    # })
    # p_val <- (sum(perm_stats >= stat_obs) + 1) / (length(perm_stats) + 1)
    # if (p_val <= alpha) {
    #   reject_count <- reject_count + 1
    # }
  }
  return(reject_count / n_sim)
}

# 计算每个mu_y对应的功效
power_vals <- sapply(mu_y_seq, compute_power)

# 绘制功效曲线
plot(mu_y_seq, power_vals, type = "l", lwd = 2, col = "blue",
     xlab = expression(mu[y]), ylab = "经验功效",
     main = "Kolmogorov-Smirnov检验的经验功效曲线",
     ylim = c(0, 1))
abline(v = mu_x, lty = 2, col = "red") # 原假设边界
abline(h = alpha, lty = 3, col = "gray") # 显著性水平

代码说明 #

1. 单边检验指定：在ks.test中明确alternative = "greater"，对应备择假设H₁: mᵧ > mₓ，此时检验统计量为F_Y(x) - F_X(x)的最大值，p值基于该统计量的分布计算。
2. 置换检验实现（可选）：如果需要手动实现随机排序的K-S检验，注释部分展示了正确逻辑：合并样本后随机分配为新的X和Y，计算置换后的统计量，通过统计量大于等于观测值的比例得到p值。
3. 功效计算：通过重复模拟生成样本，统计拒绝原假设的比例，得到经验功效。
4. 可复现性：设置set.seed保证每次运行结果一致。

结果匹配说明 #

如果之前的代码未指定alternative = "greater"，会默认使用双边检验，导致功效曲线差异极大；另外，手动实现置换检验时若错误使用双边统计量（最大绝对差），也会和ks.test的单边结果不符。修正后，两种方式（直接用ks.test单边/手动置换单边）的结果应一致。

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