如何修复 R 代码的排列部分？

Question

我已附上我的代码，其中包含我遇到的问题：

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

# Create subset to only see setosa and virginica species
iris_subset <- iris %>%
  filter(Species %in% c("setosa", "virginica"))

# Define test statistic
difference <- diff(mean(iris_subset$Sepal.Length[iris_subset$Species == "setosa"]),
                   mean(iris_subset$Sepal.Length[iris_subset$Species == "virginica"]))

# Define function
permutation_test <- function(data) {
  permutated_species <- sample(data$Species)
  permutated_diff <- diff(mean(data$Sepal.Length[permutated_species == "setosa"]),
                          mean(data$Sepal.Length[permutated_species == "virginica"]))
  return(permutated_diff)
}

# Permutation test
set.seed(13)
permutation_results <- permTS(iris_subset$Sepal.Length, group = iris_subset$Species, testStatistic = perm_test, nrepl = 10000)

# Calculate p-value
p_value <- permutation_results$p.value

# Print p-value
cat("P-value:", p_value, "\n")

我的代码中不断出现各种错误，所有这些错误都围绕排列测试部分，我不确定如何修复它，以便获得我的 p 值。

我将我的代码放入 ChatGPT 中，看看它是否可以修复我的错误，第一个是：

"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 都需要是数字。我还尝试使用库（coin）和库（exactRankTest）并重新编码我的原始代码以适应包，但仍然存在错误。

Answer 1

你让软件包已经变得简单的地方变得复杂了。两个样本检验是对均值差异的排列检验，无需定义函数来执行

permTS

已经执行的操作。来自

help("permTS")

：

对于 permTS，检验统计量等于一组平均值减去另一组平均值。

library(perm)
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union

data(iris)

# Create subset to only see setosa and virginica species
iris_subset <- iris %>%
  filter(Species %in% c("setosa", "virginica"))

# set conf level, number of Monte Carlo replications, etc
permCntl <- permControl(cm = NULL, nmc = 10^4, seed = 1234321, digits = 12,
                        p.conf.level = 0.95, setSEED = TRUE, tsmethod = "central")

# test default is "two.sided"
permutation_results <- permTS(Sepal.Length ~ Species, data = iris_subset, control = permCntl)

# Calculate p-value
p_value <- permutation_results$p.value
# Print p-value
cat("P-value:", p_value, "\n")
#> P-value: 5.882895e-17

^{创建于 2024-04-09，使用 reprex v2.1.0}