test-power.roc.test.R
In pROC: Display and Analyze ROC Curves

library(pROC)
data(aSAH)

context("power.roc.test")

# define variables shared among multiple tests here

test_that("power.roc.test basic function", {
  res <- power.roc.test(r.s100b)
  expect_equal(as.numeric(res$auc), as.numeric(r.s100b$auc))
  expect_equal(res$ncases, length(r.s100b$cases))
  expect_equal(res$ncontrols, length(r.s100b$controls))
  expect_equal(res$sig.level, 0.05)
  expect_equal(res$power, 0.9904833, tolerance = 0.000001)
})

test_that("power.roc.test with percent works", {
  res <- power.roc.test(r.s100b.percent)
  expect_equal(as.numeric(res$auc), as.numeric(r.s100b$auc))
  expect_equal(res$ncases, length(r.s100b$cases))
  expect_equal(res$ncontrols, length(r.s100b$controls))
  expect_equal(res$sig.level, 0.05)
  expect_equal(res$power, 0.9904833, tolerance = 0.000001)
})

test_that("power.roc.test with given auc function", {
  res <- power.roc.test(ncases = 41, ncontrols = 72, auc = 0.73, sig.level = 0.05)
  expect_equal(as.numeric(res$auc), 0.73)
  expect_equal(res$ncases, 41)
  expect_equal(res$ncontrols, 72)
  expect_equal(res$sig.level, 0.05)
  expect_equal(res$power, 0.9897453, tolerance = 0.000001)
})


test_that("power.roc.test sig.level can be omitted", {
  res <- power.roc.test(ncases = 41, ncontrols = 72, auc = 0.73)
  expect_equal(res$sig.level, 0.05)
  expect_equal(res$power, 0.9897453, tolerance = 0.000001)
})

test_that("power.roc.test can determine ncases & ncontrols", {
  res <- power.roc.test(auc = r.s100b$auc, sig.level = 0.05, power = 0.95, kappa = 1.7)
  expect_equal(as.numeric(res$auc), as.numeric(r.s100b$auc))
  expect_equal(res$sig.level, 0.05)
  expect_equal(res$power, 0.95)
  expect_equal(res$ncases, 29.29764, tolerance = 0.000001)
  expect_equal(res$ncontrols, 49.806, tolerance = 0.000001)
})

test_that("power.roc.test can determine sig.level", {
  res <- power.roc.test(ncases = 41, ncontrols = 72, auc = 0.73, power = 0.95, sig.level = NULL)
  expect_equal(as.numeric(res$auc), 0.73)
  expect_equal(res$ncases, 41)
  expect_equal(res$ncontrols, 72)
  expect_equal(res$power, 0.95)
  expect_equal(res$sig.level, 0.009238584, tolerance = 0.000001)
})

test_that("power.roc.test can determine AUC", {
  res <- power.roc.test(ncases = 41, ncontrols = 72, sig.level = 0.05, power = 0.95)
  expect_equal(res$ncases, 41)
  expect_equal(res$ncontrols, 72)
  expect_equal(res$power, 0.95)
  expect_equal(res$sig.level, 0.05)
  expect_equal(as.numeric(res$auc), 0.6961054, tolerance = 0.000001)
})

test_that("power.roc.test can take 2 ROC curves with DeLong variance", {
  res <- power.roc.test(r.ndka, r.wfns)
  expect_equal(res$ncases, 41)
  expect_equal(res$ncontrols, 72)
  expect_equal(as.numeric(res$auc1), as.numeric(r.ndka$auc))
  expect_equal(as.numeric(res$auc2), as.numeric(r.wfns$auc))
  expect_equal(res$power, 0.7131594, tolerance = 0.000001)
  expect_equal(res$sig.level, 0.05)
  expect_equal(res$alternative, "two.sided")
})


test_that("power.roc.test can take 2 percent ROC curves with DeLong variance", {
  res <- power.roc.test(r.ndka.percent, r.wfns.percent)
  expect_equal(res$ncases, 41)
  expect_equal(res$ncontrols, 72)
  expect_equal(as.numeric(res$auc1), as.numeric(r.ndka$auc))
  expect_equal(as.numeric(res$auc2), as.numeric(r.wfns$auc))
  expect_equal(res$power, 0.7131594, tolerance = 0.000001)
  expect_equal(res$sig.level, 0.05)
  expect_equal(res$alternative, "two.sided")
})

test_that("power.roc.test can take 2 ROC curves with Obuchowski variance", {
  res <- power.roc.test(r.ndka, r.wfns, method = "obuchowski")
  expect_equal(res$ncases, 41)
  expect_equal(res$ncontrols, 72)
  expect_equal(as.numeric(res$auc1), as.numeric(r.ndka$auc))
  expect_equal(as.numeric(res$auc2), as.numeric(r.wfns$auc))
  expect_equal(res$power, 0.8061004, tolerance = 0.000001)
  expect_equal(res$sig.level, 0.05)
  expect_equal(res$alternative, "two.sided")
})

test_that("power.roc.test ncases/ncontrols can take 2 ROC curves with DeLong variance", {
  res <- power.roc.test(r.ndka, r.wfns, power = 0.9)
  expect_equal(res$ncases, 64.77777, tolerance = 0.000001)
  expect_equal(res$ncontrols, 113.7561, tolerance = 0.000001)
  expect_equal(as.numeric(res$auc1), as.numeric(r.ndka$auc))
  expect_equal(as.numeric(res$auc2), as.numeric(r.wfns$auc))
  expect_equal(res$power, 0.9)
  expect_equal(res$sig.level, 0.05)
  expect_equal(res$alternative, "two.sided")
})

test_that("power.roc.test ncases/ncontrols can take 2 ROC curves with Obuchowski variance", {
  res <- power.roc.test(r.ndka, r.wfns, power = 0.9, method = "obuchowski")
  expect_equal(res$ncases, 53.23685, tolerance = 0.000001)
  expect_equal(res$ncontrols, 93.48911, tolerance = 0.000001)
  expect_equal(as.numeric(res$auc1), as.numeric(r.ndka$auc))
  expect_equal(as.numeric(res$auc2), as.numeric(r.wfns$auc))
  expect_equal(res$power, 0.9)
  expect_equal(res$sig.level, 0.05)
  expect_equal(res$alternative, "two.sided")
})

test_that("power.roc.test sig.level can take 2 ROC curves with DeLong variance", {
  res <- power.roc.test(r.ndka, r.wfns, power = 0.9, sig.level = NULL)
  expect_equal(res$ncases, 41)
  expect_equal(res$ncontrols, 72)
  expect_equal(as.numeric(res$auc1), as.numeric(r.ndka$auc))
  expect_equal(as.numeric(res$auc2), as.numeric(r.wfns$auc))
  expect_equal(res$power, 0.9)
  expect_equal(res$sig.level, 0.1836639, tolerance = 0.000001)
  expect_equal(res$alternative, "two.sided")
})

test_that("power.roc.test sig.level can take 2 ROC curves with Obuchowski variance", {
  res <- power.roc.test(r.ndka, r.wfns, power = 0.9, sig.level = NULL, method = "obuchowski")
  expect_equal(res$ncases, 41)
  expect_equal(res$ncontrols, 72)
  expect_equal(as.numeric(res$auc1), as.numeric(r.ndka$auc))
  expect_equal(as.numeric(res$auc2), as.numeric(r.wfns$auc))
  expect_equal(res$power, 0.9)
  expect_equal(res$sig.level, 0.1150686, tolerance = 0.000001)
  expect_equal(res$alternative, "two.sided")
})

test_that("power.roc.test works with partial AUC", {
  skip_slow()
  skip("Bootstrap cannot be tested yet")
  r.wfns.partial <- roc(aSAH$outcome, aSAH$wfns, quiet = TRUE, partial.auc = c(1, 0.9))
  r.ndka.partial <- roc(aSAH$outcome, aSAH$ndka, quiet = TRUE, partial.auc = c(1, 0.9))
  power.roc.test(r.wfns.partial, r.ndka.partial, power = 0.9)
})


test_that("power.roc.test works with partial AUC", {
  r.wfns.partial <- roc(aSAH$outcome, aSAH$wfns, quiet = TRUE, partial.auc = c(1, 0.9))
  r.ndka.partial <- roc(aSAH$outcome, aSAH$ndka, quiet = TRUE, partial.auc = c(1, 0.9))
  res <- power.roc.test(r.wfns.partial, r.ndka.partial, power = 0.9, method = "obuchowski")

  expect_equal(res$ncases, 0.5061498, tolerance = 0.000001)
  expect_equal(res$ncontrols, 0.8888484, tolerance = 0.000001)
  expect_equal(as.numeric(res$auc1), as.numeric(r.wfns.partial$auc))
  expect_equal(as.numeric(res$auc2), as.numeric(r.ndka.partial$auc))
  expect_equal(res$power, 0.9)
  expect_equal(res$sig.level, 0.05)
  expect_equal(res$alternative, "two.sided")
})

test_that("power.roc.test works with binormal parameters", {
  ob.params <- list(
    A1 = 2.6, B1 = 1, A2 = 1.9, B2 = 1, rn = 0.6, ra = 0.6, FPR11 = 0,
    FPR12 = 0.2, FPR21 = 0, FPR22 = 0.2, delta = 0.037
  )

  res1 <- power.roc.test(ob.params, power = 0.8, sig.level = 0.05)
  expect_equal(res1$ncases, 107.0238, tolerance = 0.000001)
  expect_equal(res1$ncontrols, 107.0238, tolerance = 0.000001)
  expect_equal(res1$power, 0.8)
  expect_equal(res1$sig.level, 0.05)

  res2 <- power.roc.test(ob.params, power = 0.8, sig.level = NULL, ncases = 107)
  expect_equal(res2$ncases, 107)
  expect_equal(res2$ncontrols, 107)
  expect_equal(res2$power, 0.8)
  expect_equal(res2$sig.level, 0.05004012, tolerance = 0.000001)

  res3 <- power.roc.test(ob.params, power = NULL, sig.level = 0.05, ncases = 107)
  expect_equal(res3$ncases, 107)
  expect_equal(res3$ncontrols, 107)
  expect_equal(res3$sig.level, 0.05)
  expect_equal(res3$power, 0.7999286, tolerance = 0.000001)
})

## With only binormal parameters given
# From example 2 of Obuchowski and McClish, 1997.


test_that("power.roc.test returns correct results from litterature", {
  context("Check results in Obuchowski 2004 Table 4")
  # Note: the table reports at least 10 in each cell, and adapts
  # the complement value to match kappa. So in 0.25/0.95 we have 10/40
  # although both values are < 10.
  # Note2: some values don't match exactly, specifically
  # expected.ncases[0.5, 0.6] and expected.ncontrols[4, 0.6]
  # are off by 1 (< 1%).
  kappas <- c(0.25, 0.5, 1, 2, 4)
  thetas <- c(0.6, 0.7, 0.8, 0.9, 0.95)
  expected.ncontrols <- matrix(
    c(
      84,     20,     10,     10,     10,
      101,    25,     10,     10,     10,
      135,    33,     14,     10,     10,
      203,    50,     21,     20,     20,
      # 339,	84,	40,	40,	40
      340,    84,     40,     40,     40 # Fixed
    ),
    nrow = 5, byrow = TRUE,
    dimnames = list(kappas, thetas)
  )
  expected.ncases <- matrix(
    c(
      334,    80,     40,     40,     40,
      # 201,	49,	20,	20,	20,
      202,    49,     20,     20,     20, # Fixed
      135,    33,     14,     10,     10,
      102,    25,     11,     10,     10,
      85,     21,     10,     10,     10
    ),
    nrow = 5, byrow = TRUE,
    dimnames = list(kappas, thetas)
  )

  for (kappa in kappas) {
    for (theta in thetas) {
      context(sprintf("kappa: %s, theta: %s", kappa, theta))
      pr <- power.roc.test(auc = theta, sig.level = 0.05, power = 0.9, kappa = kappa, alternative = "one.sided")
      expect_equal(max(10, ifelse(ceiling(pr$ncases) < 10, 10, 0) * kappa, ceiling(pr$ncontrols)), expected.ncontrols[as.character(kappa), as.character(theta)])
      expect_equal(max(10, ifelse(ceiling(pr$ncontrols) < 10, 10, 0) / kappa, ceiling(pr$ncases)), expected.ncases[as.character(kappa), as.character(theta)])
    }
  }
})


test_that("kappa works with a single ROC curve", {
  # kappa from data
  res <- power.roc.test(r.s100b, sig.level = 0.05, power = 0.9)
  expect_equal(res$ncases, 23.5598674)
  expect_equal(res$ncontrols, 41.3734257)
  expect_equal(res$ncases / res$ncontrols, length(r.s100b$cases) / length(r.s100b$controls))
  # set kappa
  res <- power.roc.test(r.s100b, sig.level = 0.05, power = 0.9, kappa = 1)
  expect_equal(res$ncases, 29.5697422)
  expect_equal(res$ncontrols, 29.5697422)
})


test_that("kappa works with two ROC curves", {
  # kappa from data
  res <- power.roc.test(r.s100b, r.ndka, sig.level = 0.05, power = 0.9)
  expect_equal(res$ncases, 213.117677)
  expect_equal(res$ncontrols, 374.255432)
  expect_equal(res$ncases / res$ncontrols, length(r.s100b$cases) / length(r.s100b$controls))
  # set kappa
  res <- power.roc.test(r.s100b, r.ndka, sig.level = 0.05, power = 0.9, kappa = 1)
  expect_equal(res$ncases, 213.117677)
  expect_equal(res$ncases, 213.117677)
  # ...
})