如何修复R代码中的置换检验（Permutation Test）部分？

问题重现 #

原代码调用permTS时触发以下错误：

Error in permTS.default(iris_subset$Sepal.Length, group = iris_subset$Species,  : argument "y" is missing, with no default

添加y=iris_subset$Species后，又出现类型错误：x和y需为数值型，更换其他包也未解决问题。

错误原因 #

1. perm::permTS参数逻辑错误：该函数默认需要两个独立的数值向量（对应两组样本），而非传入单个连续变量加分组参数；自定义检验统计量的格式也不符合要求，它需要接收两组数据而非整个数据框。
2. 直接传入因子类型的Species作为y，不符合函数对数值型输入的要求，导致类型不匹配报错。

修复方案 #

方案1：手动实现置换检验（最直观） #

无需依赖特定包，手动模拟置换过程：

library(dplyr)
data(iris)

# 筛选目标数据集
iris_subset <- iris %>%
  filter(Species %in% c("setosa", "virginica"))

# 计算原始观测的均值差
obs_diff <- mean(iris_subset$Sepal.Length[iris_subset$Species == "setosa"]) - 
            mean(iris_subset$Sepal.Length[iris_subset$Species == "virginica"])

# 模拟10000次置换
set.seed(13)
perm_diffs <- replicate(10000, {
  # 随机打乱分组标签
  perm_groups <- sample(iris_subset$Species)
  # 计算置换后的均值差
  mean(iris_subset$Sepal.Length[perm_groups == "setosa"]) - 
  mean(iris_subset$Sepal.Length[perm_groups == "virginica"])
})

# 计算双侧检验的p值
p_value <- sum(abs(perm_diffs) >= abs(obs_diff)) / length(perm_diffs)

# 输出结果
cat("P-value:", p_value, "\n")

方案2：使用perm包的正确调用方式 #

调整参数为两组数值向量，定义符合要求的检验统计量：

library(dplyr)
library(perm)
data(iris)

# 拆分两组独立数据
setosa_sl <- iris$Sepal.Length[iris$Species == "setosa"]
virginica_sl <- iris$Sepal.Length[iris$Species == "virginica"]

# 定义检验统计量（接收两组数值向量）
mean_diff <- function(x, y) {
  mean(x) - mean(y)
}

# 执行置换检验
set.seed(13)
perm_results <- permTS(x = setosa_sl, y = virginica_sl, 
                       testStatistic = mean_diff, nrepl = 10000)

# 输出p值
cat("P-value:", perm_results$p.value, "\n")

方案3：使用coin包实现置换检验 #

coin包支持公式接口，更适配分组数据场景：

library(dplyr)
library(coin)
data(iris)

iris_subset <- iris %>%
  filter(Species %in% c("setosa", "virginica"))

# 确保分组为因子类型
iris_subset$Species <- factor(iris_subset$Species)

# 执行置换检验（近似分布，10000次重抽样）
set.seed(13)
perm_test <- independence_test(Sepal.Length ~ Species, 
                               data = iris_subset,
                               teststat = "mean",
                               distribution = approximate(nresample = 10000))

# 提取并输出p值
p_value <- pvalue(perm_test)
cat("P-value:", p_value, "\n")

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