20.26 Wilcoxon 秩和检验 wilcox.test

单样本 Wilcoxon 秩和检验,两样本 Wilcoxon 符号秩检验,也叫 Mann-Whitney 检验

Wilcoxon Rank Sum and Signed Rank Tests

Performs one- and two-sample Wilcoxon tests on vectors of data; the latter is also known as ‘Mann-Whitney’ test.

usage(wilcox.test)
wilcox.test(x, ...)
usage("wilcox.test.default")
## Default S3 method:
wilcox.test(x, y = NULL, alternative = c("two.sided", "less", "greater"),
    mu = 0, paired = FALSE, exact = NULL, correct = TRUE, conf.int = FALSE,
    conf.level = 0.95, tol.root = 1e-04, digits.rank = Inf, ...)
usage("wilcox.test.formula")
## S3 method for class 'formula'
wilcox.test(formula, data, subset, na.action, ...)
  • coin::wilcox_test for exact, asymptotic and Monte Carlo conditional p-values, including in the presence of ties.

coin 包 (Hothorn et al. 2008) 提供大量基于秩的检验

20.26.1 ROC 曲线和 wilcox.test 检验的关系

https://github.com/xrobin/pROC/wiki/FAQ---Frequently-asked-questions#can-i-test-if-a-single-roc-curve-is-significantly-different-from-05

ROC 曲线越往左上角拱越好, AUC 是 ROC 曲线下的面积,所以 AUC 指标越接近 1 越好。

对每个标签的预测概率指定服从均匀分布,相当于随机猜测,所以最后 ROC 会接近对角线,而且样本量越大越接近,AUC 会越来越接近 0.5

再往深一点就是研究一下 R 内置的排序算法,因为计算 AUC 最核心的步骤是排序。 order 函数默认的排序方法是 auto 即当数据量较小的时候,自动选择 radix 排序,当数据量比较大的时候,自动选择 shell 排序36

# 模拟一些数据
set.seed(2019) # 设置随机数种子
N <- 10^5 # 样本量
sim_dat <- cbind.data.frame(
  pred = runif(N),
  label = rbinom(N, size = 1, prob = 0.95)
)

# 计算 auc 的函数
# dat is a data.frame as input return AUC value
comp_auc <- function(dat, show_roc = TRUE) {
  # order label by predicted probability
  dat <- dat[order(dat$pred, dat$label, decreasing = TRUE), ]

  # total samples
  n_total <- length(dat$label)

  # number of positive label 1
  n_pos <- sum(dat$label)

  # number of negative label 0
  n_neg <- n_total - n_pos

  # calculate TPR and FPR
  tpr <- cumsum(dat$label) / n_pos
  fpr <- (1:n_total - cumsum(dat$label)) / n_neg

  # calculate auc
  auc <- 0
  for (i in 1:(n_total - 1)) {
    auc <- auc + (fpr[i + 1] - fpr[i]) * tpr[i]
  }
  # show ROC curve or not?
  if (show_roc) {
    plot(fpr, tpr, type = "l")
  }
  auc
}

comp_auc(dat = sim_dat, show_roc = FALSE)
## [1] 0.5015558

模拟一个逻辑回归模型测试自编 AUC 计算程序和 R 包 pROC 计算结果

set.seed(2018)
N <- 10^4 # 样本量
x <- rnorm(N)
beta_0 <- 0.5
beta_1 <- 0.3
eta <- beta_0 + beta_1 * x
# 模拟数据集
dat <- data.frame(x = x, y = rbinom(N, 1, prob = exp(eta) / (1 + exp(eta))))

# 数据集分隔
is_train <- sample(1:nrow(dat), N * 0.7)
train <- dat[is_train, ]
test <- dat[-is_train, ]
# 模型拟合
fit <- glm(y ~ x, data = train, family = binomial(link = "logit"))
# 预测
y_pred <- predict(fit, newdata = test, type = "response")

dat2 <- data.frame(pred = y_pred, label = test$y)
# 计算 auc
comp_auc(dat = dat2, show_roc = FALSE)
## [1] 0.5850287

对比 R 包 pROC 的计算结果是一致的

pROC::auc(test$y, y_pred)

计算一下运行时间

# 100 万样本
system.time(comp_auc(dat = dat2, show_roc = FALSE))
##    user  system elapsed 
##   0.002   0.000   0.002

更多关于 auc 计算的讨论见统计之都论坛帖 https://d.cosx.org/d/419436,我感觉这个问题最后会归结到排序问题。

# https://stat.ethz.ch/pipermail/r-help/2005-April/069217.html
trap.rule <- function(x, y) sum(diff(x) * (y[-1] + y[-length(y)])) / 2

参考文献

Hothorn, Torsten, Kurt Hornik, Mark A. van de Wiel, and Achim Zeileis. 2008. “Implementing a Class of Permutation Tests: The coin Package.” Journal of Statistical Software 28 (8): 1–23. https://doi.org/10.18637/jss.v028.i08.

  1. radix 排序翻译过来叫桶排序或基数排序,详细描述见 ?sort↩︎