配对Wilcoxon符号秩检验功效分析：R函数使用验证及推荐

一、你当前代码的核心问题 #

你定义的rxy函数只生成了单组数据，完全没用到Group2的均值(1)和标准差(0.01)，也没有体现配对设计中同一受试者两组数据的对应关系，这直接导致模拟结果不可靠，对n=16时的功效值存疑是合理的。

正确的做法是定义一个能同时生成配对的两组数据的函数，示例代码如下：

library(MKpower)

# 生成配对数据的函数：对应同一受试者的两种处理结果
rx_pair <- function(n) {
  group1 <- rnorm(n, mean = 0.9, sd = 1.4)
  group2 <- rnorm(n, mean = 1, sd = 0.01)
  cbind(group1, group2)
}

# 若考虑同一受试者两次测量的相关性，可使用mvtnorm生成多元正态数据
library(mvtnorm)
rx_pair_cor <- function(n) {
  mu <- c(0.9, 1)
  # 假设两组相关系数为0.5，可根据研究实际调整
  sigma <- matrix(c(1.4^2, 0.5*1.4*0.01, 0.5*1.4*0.01, 0.01^2), nrow=2)
  rmvnorm(n, mean=mu, sigma=sigma)
}

# 重新运行功效模拟
sim.ssize.wilcox.test(rx = rx_pair_cor, n.min=8, n.max=30, step.size=2, iter=10000, type="paired")

二、可计算配对Wilcoxon检验功效的R方法/函数 #

1. 自定义模拟函数（最灵活可控） #

直接编写循环模拟逻辑，完全自定义数据生成和检验流程：

wilcox_paired_power <- function(n, mu1, sd1, mu2, sd2, cor=0, alpha=0.05, iter=10000) {
  library(mvtnorm)
  sig_count <- 0
  mu <- c(mu1, mu2)
  sigma <- matrix(c(sd1^2, cor*sd1*sd2, cor*sd1*sd2, sd2^2), nrow=2)
  
  for (i in 1:iter) {
    data <- rmvnorm(n, mean=mu, sigma=sigma)
    res <- wilcox.test(data[,1], data[,2], paired=TRUE)
    if (res$p.value < alpha) sig_count <- sig_count + 1
  }
  sig_count / iter
}

# 计算n=16时的功效（示例）
wilcox_paired_power(n=16, mu1=0.9, sd1=1.4, mu2=1, sd2=0.01, cor=0.5)

2. exactRankTests包结合模拟 #

该包提供精确Wilcoxon检验，适合小样本场景的功效模拟：

library(exactRankTests)

wilcox_exact_power <- function(n, mu1, sd1, mu2, sd2, cor=0, alpha=0.05, iter=10000) {
  library(mvtnorm)
  sig_count <- 0
  mu <- c(mu1, mu2)
  sigma <- matrix(c(sd1^2, cor*sd1*sd2, cor*sd1*sd2, sd2^2), nrow=2)
  
  for (i in 1:iter) {
    data <- rmvnorm(n, mean=mu, sigma=sigma)
    res <- wilcox.exact(data[,1], data[,2], paired=TRUE)
    if (res$p.value < alpha) sig_count <- sig_count + 1
  }
  sig_count / iter
}

# 计算n=16的功效
wilcox_exact_power(n=16, mu1=0.9, sd1=1.4, mu2=1, sd2=0.01, cor=0.5)

3. coin包置换检验结合模拟 #

coin包支持灵活的秩检验框架，可用于置换检验的功效模拟：

library(coin)

wilcox_coin_power <- function(n, mu1, sd1, mu2, sd2, cor=0, alpha=0.05, iter=10000) {
  library(mvtnorm)
  sig_count <- 0
  mu <- c(mu1, mu2)
  sigma <- matrix(c(sd1^2, cor*sd1*sd2, cor*sd1*sd2, sd2^2), nrow=2)
  
  for (i in 1:iter) {
    data <- rmvnorm(n, mean=mu, sigma=sigma)
    df <- data.frame(x=data[,1], y=data[,2], id=1:n)
    res <- wilcoxsign_test(x ~ y | id, data=df)
    if (pvalue(res) < alpha) sig_count <- sig_count + 1
  }
  sig_count / iter
}

# 计算功效
wilcox_coin_power(n=16, mu1=0.9, sd1=1.4, mu2=1, sd2=0.01, cor=0.5)

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