suppressMessages(library(dplyr))

test_that("kB factors are correct for normal distribution", {
  cmh17_factors <- matrix(
    c(2, 20.581, 36, 1.725, 70, 1.582, 104, 1.522,
      3, 6.157, 37, 1.718, 71, 1.579, 105, 1.521,
      4, 4.163, 38, 1.711, 72, 1.577, 106, 1.519,
      5, 3.408, 39, 1.704, 73, 1.575, 107, 1.518,
      6, 3.007, 40, 1.698, 74, 1.572, 108, 1.517,
      7, 2.756, 41, 1.692, 75, 1.570, 109, 1.516,
      8, 2.583, 42, 1.686, 76, 1.568, 110, 1.515,
      9, 2.454, 43, 1.680, 77, 1.566, 111, 1.513,
      10, 2.355, 44, 1.675, 78, 1.564, 112, 1.512,
      11, 2.276, 45, 1.669, 79, 1.562, 113, 1.511,
      12, 2.211, 46, 1.664, 80, 1.560, 114, 1.510,
      13, 2.156, 47, 1.660, 81, 1.558, 115, 1.509,
      14, 2.109, 48, 1.655, 82, 1.556, 116, 1.508,
      15, 2.069, 49, 1.650, 83, 1.554, 117, 1.507,
      16, 2.034, 50, 1.646, 84, 1.552, 118, 1.506,
      17, 2.002, 51, 1.642, 85, 1.551, 119, 1.505,
      18, 1.974, 52, 1.638, 86, 1.549, 120, 1.504,
      19, 1.949, 53, 1.634, 87, 1.547, 121, 1.503,
      20, 1.927, 54, 1.630, 88, 1.545, 122, 1.502,
      21, 1.906, 55, 1.626, 89, 1.544, 123, 1.501,
      22, 1.887, 56, 1.623, 90, 1.542, 124, 1.500,
      23, 1.870, 57, 1.619, 91, 1.540, 125, 1.499,
      24, 1.854, 58, 1.616, 92, 1.539, 126, 1.498,
      25, 1.839, 59, 1.613, 93, 1.537, 127, 1.497,
      26, 1.825, 60, 1.609, 94, 1.536, 128, 1.496,
      27, 1.812, 61, 1.606, 95, 1.534, 129, 1.495,
      28, 1.800, 62, 1.603, 96, 1.533, 130, 1.494,
      29, 1.789, 63, 1.600, 97, 1.531, 131, 1.493,
      30, 1.778, 64, 1.597, 98, 1.530, 132, 1.492,
      31, 1.768, 65, 1.595, 99, 1.529, 133, 1.492,
      32, 1.758, 66, 1.592, 100, 1.527, 134, 1.491,
      33, 1.749, 67, 1.589, 101, 1.526, 135, 1.490,
      34, 1.741, 68, 1.587, 102, 1.525, 136, 1.489,
      35, 1.733, 69, 1.584, 103, 1.523, 137, 1.488),
    ncol = 2, byrow = TRUE
  ) %>%
    as.data.frame() %>%
    rename(n = V1) %>%
    rename(kb = V2) %>%
    filter(n <= 95) %>%
    rowwise() %>%
    mutate(calc_kb = k_factor_normal(n, p = 0.90, conf = 0.95)) %>%
    mutate(check = expect_lte(abs(calc_kb - kb), expected = 0.002,
                              label = paste0("Validation failure for ", n, ".",
                                             "CMH-17 gives kB=", kb, ",",
                                             "library gives kB=", calc_kb)))
})

test_that("kA factors are correct for normal distribution", {
  cmh17_factors <- matrix(
    c(2, 37.094, 36, 2.983, 70, 2.765, 104, 2.676,
      3, 10.553, 37, 2.972, 71, 2.762, 105, 2.674,
      4, 7.042, 38, 2.961, 72, 2.758, 106, 2.672,
      5, 5.741, 39, 2.951, 73, 2.755, 107, 2.671,
      6, 5.062, 40, 2.941, 74, 2.751, 108, 2.669,
      7, 4.642, 41, 2.932, 75, 2.748, 109, 2.667,
      8, 4.354, 42, 2.923, 76, 2.745, 110, 2.665,
      9, 4.143, 43, 2.914, 77, 2.742, 111, 2.663,
      10, 3.981, 44, 2.906, 78, 2.739, 112, 2.662,
      11, 3.852, 45, 2.898, 79, 2.736, 113, 2.660,
      12, 3.747, 46, 2.890, 80, 2.733, 114, 2.658,
      13, 3.659, 47, 2.883, 81, 2.730, 115, 2.657,
      14, 3.585, 48, 2.876, 82, 2.727, 116, 2.655,
      15, 3.520, 49, 2.869, 83, 2.724, 117, 2.654,
      16, 3.464, 50, 2.862, 84, 2.721, 118, 2.652,
      17, 3.414, 51, 2.856, 85, 2.719, 119, 2.651,
      18, 3.370, 52, 2.850, 86, 2.716, 120, 2.649,
      19, 3.331, 53, 2.844, 87, 2.714, 121, 2.648,
      20, 3.295, 54, 2.838, 88, 2.711, 122, 2.646,
      21, 3.263, 55, 2.833, 89, 2.709, 123, 2.645,
      22, 3.233, 56, 2.827, 90, 2.706, 124, 2.643,
      23, 3.206, 57, 2.822, 91, 2.704, 125, 2.642,
      24, 3.181, 58, 2.817, 92, 2.701, 126, 2.640,
      25, 3.158, 59, 2.812, 93, 2.699, 127, 2.639,
      26, 3.136, 60, 2.807, 94, 2.697, 128, 2.638,
      27, 3.116, 61, 2.802, 95, 2.695, 129, 2.636,
      28, 3.098, 62, 2.798, 96, 2.692, 130, 2.635,
      29, 3.080, 63, 2.793, 97, 2.690, 131, 2.634,
      30, 3.064, 64, 2.789, 98, 2.688, 132, 2.632,
      31, 3.048, 65, 2.785, 99, 2.686, 133, 2.631,
      32, 3.034, 66, 2.781, 100, 2.684, 134, 2.630,
      33, 3.020, 67, 2.777, 101, 2.682, 135, 2.628,
      34, 3.007, 68, 2.773, 102, 2.680, 136, 2.627,
      35, 2.995, 69, 2.769, 103, 2.678, 137, 2.626),
    ncol = 2, byrow = TRUE
  ) %>%
    as.data.frame() %>%
    rename(n = V1) %>%
    rename(ka = V2) %>%
    filter(n <= 75) %>%
    rowwise() %>%
    mutate(calc_ka = k_factor_normal(n, p = 0.99, conf = 0.95)) %>%
    mutate(check = expect_lte(abs(calc_ka - ka), expected = 0.002,
                              label = paste0("Validation failure for ", n, ".",
                                             "CMH-17 gives kA=", ka, ",",
                                             "library gives kA=", calc_ka)))
})

test_that("(normal) basis value equals mean when sd = 0", {
  expect_equal(
    (data.frame(x = rep(100, 10)) %>%
      basis_normal(x, p = 0.9, conf = 0.95,
                   override = c("anderson_darling_normal",
                                "outliers_within_batch",
                                "between_batch_variability")))$basis,
    100
  )
})

test_that("(normal) basis value approx equals percentile for large samples", {
  m <- 100
  s <- 5
  set.seed(100)  # make sure that this doesn't fail by pure chance
  q <- qnorm(0.10, m, s, lower.tail = TRUE)
  basis <- (data.frame(x = rnorm(50, m, s)) %>%
    basis_normal(x, p = 0.90, conf = 0.95,
                 override = c("outliers_within_batch",
                              "between_batch_variability")))$basis
  expect_lt(abs(basis - q), 0.1)
})

test_that("printing of basis objects works as expected", {
  set.seed(100)
  x <- c(runif(25))

  expect_output(
    print(basis_normal(x = x, p = 0.9, conf = 0.95,
                       override = c("outliers_within_batch",
                                    "between_batch_variability"))),
    "B-Basis"
  )

  expect_output(
    print(basis_normal(x = x, p = 0.99, conf = 0.95,
                       override = c("outliers_within_batch",
                                    "between_batch_variability"))),
    "A-Basis"
  )

  expect_output(
    print(basis_normal(x = x, p = 0.9, conf = 0.9,
                       override = c("outliers_within_batch",
                                    "between_batch_variability"))),
    "[^AB-]Basis"
  )

  expect_error(
    print.basis(x)  # give it the wrong type
  )
})

test_that("normal basis value matches STAT17/ASAP result", {
  data <- c(
    137.4438,
    139.5395,
    150.89,
    141.4474,
    141.8203,
    151.8821,
    143.9245,
    132.9732,
    136.6419,
    138.1723,
    148.7668,
    143.283,
    143.5429,
    141.7023,
    137.4732,
    152.338,
    144.1589,
    128.5218
  )

  res <- basis_normal(x = data, p = 0.9, conf = 0.95,
                      override = c("outliers_within_batch",
                                   "between_batch_variability"))
  expect_equal(res$basis, 129.287, tolerance = 0.0005)
  expect_output(print(res), "b-basis.*129\\.2", ignore.case = TRUE)
  expect_output(print(res), "normal", ignore.case = TRUE)

  res <- basis_normal(x = data, p = 0.99, conf = 0.95,
                      override = c("outliers_within_batch",
                                   "between_batch_variability"))
  expect_equal(res$basis, 120.336, tolerance = 0.0005)
  expect_output(print(res), "a-basis.*120\\.3", ignore.case = TRUE)
  expect_output(print(res), "normal", ignore.case = TRUE)

  expect_match(res$distribution, "normal", ignore.case = TRUE)
})

test_that("normal basis values produce expected diagnostic failures", {
  set.seed(100)
  x <- c(runif(25), runif(25, max = 2), 200)
  batch <- c(rep("A", 25), rep("B", 26))

  expect_warning(
    expect_warning(
      expect_warning(
        expect_warning(
          res <- basis_normal(x = x, batch = batch),
          "outliers_within_batch"
        ),
        "between_batch_variability"
      ),
      "outliers"
    ),
    "anderson_darling_normal"
  )

  # Check that res$... contains the correct value
  expect_equal(res$batch, batch)
  expect_equal(res$diagnostic_failures,
               c("outliers_within_batch",
                 "between_batch_variability",
                 "outliers",
                 "anderson_darling_normal"))
  expect_length(res$override, 0)

  expect_output(print(res),
                regexp = paste("failed.+",
                               "outliers_within_batch",
                               "between_batch_variability",
                               "outliers",
                               "anderson_darling_normal",
                               sep = ".+"))
  output <- capture_output(print(res))
  expect_false(grepl("overridden", output, ignore.case = TRUE))

  # overriding the diagnostics should eliminate the warnings
  res <- basis_normal(x = x, batch = batch,
                      override = c("outliers_within_batch",
                                   "between_batch_variability",
                                   "outliers",
                                   "anderson_darling_normal"))

  expect_equal(res$override,
               c("outliers_within_batch",
                 "between_batch_variability",
                 "outliers",
                 "anderson_darling_normal"))
  expect_length(res$diagnostic_failures, 0)

  expect_output(print(res),
                regexp = paste("overridden.+",
                               "outliers_within_batch",
                               "between_batch_variability",
                               "outliers",
                               "anderson_darling_normal",
                               sep = ".+"))
  output <- capture_output(print(res))
  expect_false(grepl("failed", output, ignore.case = TRUE))

  # overriding the diagnostics using "all" should do the same thing
  res <- basis_normal(x = x, batch = batch,
                      override = "all")
  expect_equal(res$override,
               c("outliers_within_batch",
                 "between_batch_variability",
                 "outliers",
                 "anderson_darling_normal"))
  expect_length(res$diagnostic_failures, 0)

  # call basis_normal without batch
  expect_warning(
    expect_warning(
      expect_message(
        expect_message(
          res <- basis_normal(x = x),
          "outliers_within_batch"
        ),
        "between_batch_variability"
      ),
      "outliers"
    ),
    "anderson_darling_normal"
  )

  # Check that res$... contains the correct value
  expect_equal(res$diagnostic_failures,
               c("outliers",
                 "anderson_darling_normal"))
  expect_length(res$override, 0)

  # overriding the diagnostics should eliminate the warnings
  res <- basis_normal(x = x,
                      override = c("outliers",
                                   "anderson_darling_normal",
                                   "outliers_within_batch",
                                   "between_batch_variability"))

  expect_equal(res$override,
               c("outliers",
                 "anderson_darling_normal",
                 "outliers_within_batch",
                 "between_batch_variability"))
  expect_length(res$diagnostic_failures, 0)
})

test_that("log-normal basis value matches STAT17 result", {
  data <- c(
    137.4438,
    139.5395,
    150.89,
    141.4474,
    141.8203,
    151.8821,
    143.9245,
    132.9732,
    136.6419,
    138.1723,
    148.7668,
    143.283,
    143.5429,
    141.7023,
    137.4732,
    152.338,
    144.1589,
    128.5218
  )

  res <- basis_lognormal(x = data, p = 0.9, conf = 0.95, override = "all")
  expect_equal(res$basis, 129.664, tolerance = 0.0005)
  expect_output(print(res), "b-basis.*129.6", ignore.case = TRUE)
  expect_output(print(res), "normal", ignore.case = TRUE)
  expect_output(print(res), "log", ignore.case = TRUE)

  res <- basis_lognormal(x = data, p = 0.99, conf = 0.95, override = "all")
  expect_equal(res$basis, 121.710, tolerance = 0.0005)
  expect_output(print(res), "a-basis.*121.7", ignore.case = TRUE)
  expect_output(print(res), "normal", ignore.case = TRUE)
  expect_output(print(res), "log", ignore.case = TRUE)

  expect_match(res$distribution, "log", ignore.case = TRUE)
  expect_match(res$distribution, "normal", ignore.case = TRUE)
})

test_that("lognormal basis values produce expected diagnostic failures", {
  set.seed(100)
  x <- c(runif(25), runif(25, max = 2), 200)
  batch <- c(rep("A", 25), rep("B", 26))

  expect_warning(
    expect_warning(
      expect_warning(
        expect_warning(
          res <- basis_lognormal(x = x, batch = batch),
          "outliers_within_batch"
        ),
        "between_batch_variability"
      ),
      "outliers"
    ),
    "anderson_darling_lognormal"
  )

  # Check that res$... contains the correct value
  expect_equal(res$batch, batch)
  expect_equal(res$diagnostic_failures,
               c("outliers_within_batch",
                 "between_batch_variability",
                 "outliers",
                 "anderson_darling_lognormal"))
  expect_length(res$override, 0)

  # overriding the diagnostics should eliminate the warnings
  res <- basis_lognormal(x = x, batch = batch,
                      override = c("outliers_within_batch",
                                   "between_batch_variability",
                                   "outliers",
                                   "anderson_darling_lognormal"))

  expect_equal(res$override,
               c("outliers_within_batch",
                 "between_batch_variability",
                 "outliers",
                 "anderson_darling_lognormal"))
  expect_length(res$diagnostic_failures, 0)

  # overriding the diagnostics with "all" should do the same thing
  res <- basis_lognormal(x = x, batch = batch,
                         override = "all")

  expect_equal(res$override,
               c("outliers_within_batch",
                 "between_batch_variability",
                 "outliers",
                 "anderson_darling_lognormal"))
  expect_length(res$diagnostic_failures, 0)

  # call basis_normal without batch
  expect_warning(
    expect_warning(
      expect_message(
        expect_message(
          res <- basis_lognormal(x = x),
          "outliers_within_batch"
        ),
        "between_batch_variability"
      ),
      "outliers"
    ),
    "anderson_darling_lognormal"
  )

  # Check that res$... contains the correct value
  expect_equal(res$diagnostic_failures,
               c("outliers",
                 "anderson_darling_lognormal"))
  expect_length(res$override, 0)

  # overriding the diagnostics should eliminate the warnings
  res <- basis_lognormal(x = x,
                      override = c("outliers",
                                   "anderson_darling_lognormal",
                                   "outliers_within_batch",
                                   "between_batch_variability"))

  expect_equal(res$override,
               c("outliers",
                 "anderson_darling_lognormal",
                 "outliers_within_batch",
                 "between_batch_variability"))
  expect_length(res$diagnostic_failures, 0)
})

test_that("Weibull basis value matches STAT17 result", {
  data <- c(
    137.4438,
    139.5395,
    150.89,
    141.4474,
    141.8203,
    151.8821,
    143.9245,
    132.9732,
    136.6419,
    138.1723,
    148.7668,
    143.283,
    143.5429,
    141.7023,
    137.4732,
    152.338,
    144.1589,
    128.5218
  )

  # stat17 B-Basis: 125.441
  # stat17 A-Basis: 109.150

  res <- basis_weibull(x = data, p = 0.9, conf = 0.95,
                       override = c("outliers_within_batch",
                                    "between_batch_variability"))
  expect_equal(res$basis, 125.441, tolerance = 0.3)
  expect_output(print(res), "b-basis.*125", ignore.case = TRUE)
  expect_output(print(res), "weibull", ignore.case = TRUE)

  res <- basis_weibull(x = data, p = 0.99, conf = 0.95,
                       override = c("outliers_within_batch",
                                    "between_batch_variability"))
  expect_equal(res$basis, 109.150, tolerance = 0.6)
  expect_output(print(res), "a-basis.*109", ignore.case = TRUE)
  expect_output(print(res), "weibull", ignore.case = TRUE)

  expect_match(res$distribution, "weibull", ignore.case = TRUE)
})

test_that("weibull basis values produce expected diagnostic failures", {
  set.seed(100)
  x <- c(rnorm(10, 100, 2), rnorm(10, 103, 2), 120)
  batch <- c(rep("A", 10), rep("B", 11))

  expect_warning(
    expect_warning(
      expect_warning(
        expect_warning(
          res <- basis_weibull(x = x, batch = batch),
          "outliers_within_batch"
        ),
        "between_batch_variability"
      ),
      "outliers"
    ),
    "anderson_darling_weibull"
  )

  # Check that res$... contains the correct value
  expect_equal(res$batch, batch)
  expect_equal(res$diagnostic_failures,
               c("outliers_within_batch",
                 "between_batch_variability",
                 "outliers",
                 "anderson_darling_weibull"))
  expect_length(res$override, 0)

  # overriding the diagnostics should eliminate the warnings
  res <- basis_weibull(x = x, batch = batch,
                         override = c("outliers_within_batch",
                                      "between_batch_variability",
                                      "outliers",
                                      "anderson_darling_weibull"))

  expect_equal(res$override,
               c("outliers_within_batch",
                 "between_batch_variability",
                 "outliers",
                 "anderson_darling_weibull"))
  expect_length(res$diagnostic_failures, 0)

  # overriding the diagnostics with "all" should do the same thing
  res <- basis_weibull(x = x, batch = batch,
                       override = "all")

  expect_equal(res$override,
               c("outliers_within_batch",
                 "between_batch_variability",
                 "outliers",
                 "anderson_darling_weibull"))
  expect_length(res$diagnostic_failures, 0)

  # call basis_normal without batch
  expect_warning(
    expect_warning(
      expect_message(
        expect_message(
          res <- basis_weibull(x = x),
          "outliers_within_batch"
        ),
        "between_batch_variability"
      ),
      "outliers"
    ),
    "anderson_darling_weibull"
  )

  # Check that res$... contains the correct value
  expect_equal(res$diagnostic_failures,
               c("outliers",
                 "anderson_darling_weibull"))
  expect_length(res$override, 0)

  # overriding the diagnostics should eliminate the warnings
  res <- basis_weibull(x = x,
                         override = c("outliers",
                                      "anderson_darling_weibull",
                                      "outliers_within_batch",
                                      "between_batch_variability"))

  expect_equal(res$override,
               c("outliers",
                 "anderson_darling_weibull",
                 "outliers_within_batch",
                 "between_batch_variability"))
  expect_length(res$diagnostic_failures, 0)
})

test_that("Non-parametric (small sample) basis value matches STAT17 result", {
  data <- c(
    137.4438,
    139.5395,
    150.89,
    141.4474,
    141.8203,
    151.8821,
    143.9245,
    132.9732,
    136.6419,
    138.1723,
    148.7668,
    143.283,
    143.5429,
    141.7023,
    137.4732,
    152.338,
    144.1589,
    128.5218
  )

  res <- basis_hk_ext(x = data, p = 0.9, conf = 0.95,
                      method = "optimum-order",
                      override = c("outliers_within_batch",
                                   "between_batch_variability"))
  expect_equal(res$basis, 124.156, tolerance = 0.002)
  expect_output(print(res), "b-basis.*124", ignore.case = TRUE)
  expect_output(print(res), "nonparametric", ignore.case = TRUE)
  expect_match(res$distribution, "nonparametric.*optimum", ignore.case = TRUE)

  res <- basis_hk_ext(x = data, p = 0.99, conf = 0.95,
                      method = "woodward-frawley",
                      override = c("outliers_within_batch",
                                   "between_batch_variability"))
  expect_equal(res$basis, 99.651, tolerance = 0.002)
  expect_output(print(res), "a-basis.*99", ignore.case = TRUE)
  expect_output(print(res), "nonparametric", ignore.case = TRUE)
  expect_match(res$distribution,
               "nonparametric.*Woodward-Frawley", ignore.case = TRUE)

  expect_error(basis_hk_ext(x = data, method = "something invalid",
                            override = c("outliers_within_batch",
                                         "between_batch_variability")))
})

test_that("non-para (small) basis values produce expected diag failures", {
  set.seed(100)
  x_small <- c(rnorm(10, 100, 2), rnorm(10, 103, 2), 120)
  batch_small <- c(rep("A", 10), rep("B", 11))
  x_large <- c(rnorm(200, 100, 2), rnorm(100, 103, 2), 120)
  batch_large <- c(rep("A", 200), rep("B", 101))

  # woodward-frawley is only for A-Basis. Should fail if we calculate B-Basis
  expect_warning(
    expect_warning(
      expect_warning(
        expect_warning(
          expect_warning(
            res <- basis_hk_ext(
              x = x_large, batch = batch_large, method = "woodward-frawley"),
            "outliers_within_batch"
          ),
          "between_batch_variability"
        ),
        "outliers"
      ),
      "correct_method_used"
    ),
    "sample_size"
  )

  # Check that res$... contains the correct value
  expect_equal(res$batch, batch_large)
  expect_equal(res$diagnostic_failures,
               c("outliers_within_batch",
                 "between_batch_variability",
                 "outliers",
                 "correct_method_used",
                 "sample_size"))
  expect_length(res$override, 0)

  # overriding the diagnostics should eliminate the warnings
  res <- basis_hk_ext(x = x_large, batch = batch_large,
                      method = "woodward-frawley",
                      override = c("outliers_within_batch",
                                   "between_batch_variability",
                                   "outliers",
                                   "correct_method_used",
                                   "sample_size"))

  expect_equal(res$override,
               c("outliers_within_batch",
                 "between_batch_variability",
                 "outliers",
                 "correct_method_used",
                 "sample_size"))
  expect_length(res$diagnostic_failures, 0)

  # overriding the diagnostics with "all" should do the same thing
  res <- basis_hk_ext(x = x_large, batch = batch_large,
                      method = "woodward-frawley",
                      override = "all")
  expect_equal(res$override,
               c("outliers_within_batch",
                 "between_batch_variability",
                 "outliers",
                 "correct_method_used",
                 "sample_size"))
  expect_length(res$diagnostic_failures, 0)

  # optimum-order is only for B-Basis. Should fail if we calculate A-Basis
  expect_warning(
    expect_warning(
      expect_warning(
        expect_warning(
          expect_warning(
            res <- basis_hk_ext(
              x = x_large, batch = batch_large, method = "optimum-order",
              p = 0.99, conf = 0.95),
            "outliers_within_batch"
          ),
          "between_batch_variability"
        ),
        "outliers"
      ),
      "correct_method_used"
    ),
    "sample_size"
  )

  # call basis_normal without batch
  expect_warning(
    expect_message(
      expect_message(
        res <- basis_hk_ext(x = x_small, method = "optimum-order"),
        "outliers_within_batch"
      ),
      "between_batch_variability"
    ),
    "outliers"
  )

  # Check that res$... contains the correct value
  expect_equal(res$diagnostic_failures,
               c("outliers"))
  expect_length(res$override, 0)

  # overriding the diagnostics should eliminate the warnings
  res <- basis_hk_ext(x = x_small, method = "optimum-order",
                      override = c("outliers",
                                   "outliers_within_batch",
                                   "between_batch_variability"))

  expect_equal(res$override,
               c("outliers",
                 "outliers_within_batch",
                 "between_batch_variability"))
  expect_length(res$diagnostic_failures, 0)
})

test_that("Non-parametric (large sample) basis value matches STAT17 result", {
  data <- c(
    137.3603, 135.6665, 136.6914, 154.7919, 159.2037, 137.3277, 128.821,
    138.6304, 138.9004, 147.4598, 148.6622, 144.4948, 131.0851, 149.0203,
    131.8232, 146.4471, 123.8124, 126.3105, 140.7609, 134.4875, 128.7508,
    117.1854, 129.3088, 141.6789, 138.4073, 136.0295, 128.4164, 141.7733,
    134.455,  122.7383, 136.9171, 136.9232, 138.8402, 152.8294, 135.0633,
    121.052,  131.035,  138.3248, 131.1379, 147.3771, 130.0681, 132.7467,
    137.1444, 141.662,  146.9363, 160.7448, 138.5511, 129.1628, 140.2939,
    144.8167, 156.5918, 132.0099, 129.3551, 136.6066, 134.5095, 128.2081,
    144.0896, 141.8029, 130.0149, 140.8813, 137.7864
  )

  res <- basis_nonpara_large_sample(x = data, p = 0.9, conf = 0.95,
                                    override = "all")
  expect_equal(res$basis, 122.738297, tolerance = 0.005)
  expect_output(print(res), "b-basis.*122", ignore.case = TRUE)
  expect_output(print(res), "nonparametric", ignore.case = TRUE)
  expect_match(res$distribution, "nonparametric.*large", ignore.case = TRUE)
})

test_that("non-para (large) basis values produce expected diag failures", {
  set.seed(100)
  x_small <- c(rnorm(13, 100, 2), rnorm(13, 103, 2), 120)
  batch_small <- c(rep("A", 13), rep("B", 14))
  x_large <- c(rnorm(200, 100, 2), rnorm(100, 103, 2), 120)
  batch_large <- c(rep("A", 200), rep("B", 101))

  expect_warning(
    expect_warning(
      expect_warning(
        res <- basis_nonpara_large_sample(
          x = x_large, batch = batch_large),
        "outliers_within_batch"
      ),
      "between_batch_variability"
    ),
    "outliers"
  )

  # Check that res$... contains the correct value
  expect_equal(res$batch, batch_large)
  expect_equal(res$diagnostic_failures,
               c("outliers_within_batch",
                 "between_batch_variability",
                 "outliers"))
  expect_length(res$override, 0)

  # overriding the diagnostics should eliminate the warnings
  res <- basis_nonpara_large_sample(x = x_large, batch = batch_large,
                      override = c("outliers_within_batch",
                                   "between_batch_variability",
                                   "outliers"))

  expect_equal(res$override,
               c("outliers_within_batch",
                 "between_batch_variability",
                 "outliers"))
  expect_length(res$diagnostic_failures, 0)

  # overriding the diagnostics with "all" should do the same thing
  res <- basis_nonpara_large_sample(x = x_large, batch = batch_large,
                                    override = "all")
  expect_equal(res$override,
               c("outliers_within_batch",
                 "between_batch_variability",
                 "outliers",
                 "sample_size"))
  expect_length(res$diagnostic_failures, 0)

  expect_warning(
    expect_warning(
      expect_warning(
        res <- basis_nonpara_large_sample(
          x = x_large, batch = batch_large,
          p = 0.99, conf = 0.95),
        "outliers_within_batch"
      ),
      "between_batch_variability"
    ),
    "outliers"
  )

  # call basis_normal without batch
  expect_warning(
    expect_message(
      expect_message(
        res <- basis_nonpara_large_sample(x = x_large),
        "outliers_within_batch"
      ),
      "between_batch_variability"
    ),
    "outliers"
  )

  # Check that res$... contains the correct value
  expect_equal(res$diagnostic_failures,
               c("outliers"))
  expect_length(res$override, 0)

  # overriding the diagnostics should eliminate the warnings
  res <- basis_nonpara_large_sample(x = x_large,
                      override = c("outliers",
                                   "outliers_within_batch",
                                   "between_batch_variability"))

  expect_equal(res$override,
               c("outliers", "outliers_within_batch",
                 "between_batch_variability"))
  expect_length(res$diagnostic_failures, 0)
})

# data from CMH-17-1G Section 8.3.11.1.1
cmh_17_8_3_11_1_1 <- tribble(
  ~batch, ~strength, ~condition,
  1, 118.3774604, "CTD", 1, 84.9581364, "RTD", 1, 83.7436035, "ETD",
  1, 123.6035612, "CTD", 1, 92.4891822, "RTD", 1, 84.3831677, "ETD",
  1, 115.2238092, "CTD", 1, 96.8212659, "RTD", 1, 94.8030433, "ETD",
  1, 112.6379744, "CTD", 1, 109.030325, "RTD", 1, 94.3931537, "ETD",
  1, 116.5564277, "CTD", 1, 97.8212659, "RTD", 1, 101.702222, "ETD",
  1, 123.1649896, "CTD", 1, 100.921519, "RTD", 1, 86.5372121, "ETD",
  2, 128.5589027, "CTD", 1, 103.699444, "RTD", 1, 92.3772684, "ETD",
  2, 113.1462103, "CTD", 2, 93.7908212, "RTD", 2, 89.2084024, "ETD",
  2, 121.4248107, "CTD", 2, 107.526709, "RTD", 2, 100.686001, "ETD",
  2, 134.3241906, "CTD", 2, 94.5769704, "RTD", 2, 81.0444192, "ETD",
  2, 129.6405117, "CTD", 2, 93.8831373, "RTD", 2, 91.3398070, "ETD",
  2, 117.9818658, "CTD", 2, 98.2296605, "RTD", 2, 93.1441939, "ETD",
  3, 115.4505226, "CTD", 2, 111.346590, "RTD", 2, 85.8204168, "ETD",
  3, 120.0369467, "CTD", 2, 100.817538, "RTD", 3, 94.8966273, "ETD",
  3, 117.1631088, "CTD", 3, 100.382203, "RTD", 3, 95.8068520, "ETD",
  3, 112.9302797, "CTD", 3, 91.5037811, "RTD", 3, 86.7842252, "ETD",
  3, 117.9114501, "CTD", 3, 100.083233, "RTD", 3, 94.4011973, "ETD",
  3, 120.1900159, "CTD", 3, 95.6393615, "RTD", 3, 96.7231171, "ETD",
  3, 110.7295966, "CTD", 3, 109.304779, "RTD", 3, 89.9010384, "ETD",
  3, 100.078562, "RTD", 3, 99.1205847, "RTD", 3, 89.3672306, "ETD",
  1, 106.357525, "ETW", 1, 99.0239966, "ETW2",
  1, 105.898733, "ETW", 1, 103.341238, "ETW2",
  1, 88.4640082, "ETW", 1, 100.302130, "ETW2",
  1, 103.901744, "ETW", 1, 98.4634133, "ETW2",
  1, 80.2058219, "ETW", 1, 92.2647280, "ETW2",
  1, 109.199597, "ETW", 1, 103.487693, "ETW2",
  1, 61.0139431, "ETW", 1, 113.734763, "ETW2",
  2, 99.3207107, "ETW", 2, 108.172659, "ETW2",
  2, 115.861770, "ETW", 2, 108.426732, "ETW2",
  2, 82.6133082, "ETW", 2, 116.260375, "ETW2",
  2, 85.3690411, "ETW", 2, 121.049610, "ETW2",
  2, 115.801622, "ETW", 2, 111.223082, "ETW2",
  2, 44.3217741, "ETW", 2, 104.574843, "ETW2",
  2, 117.328077, "ETW", 2, 103.222552, "ETW2",
  2, 88.6782903, "ETW", 3, 99.3918538, "ETW2",
  3, 107.676986, "ETW", 3, 87.3421658, "ETW2",
  3, 108.960241, "ETW", 3, 102.730741, "ETW2",
  3, 116.122640, "ETW", 3, 96.3694916, "ETW2",
  3, 80.2334815, "ETW", 3, 99.5946088, "ETW2",
  3, 106.145570, "ETW", 3, 97.0712407, "ETW2",
  3, 104.667866, "ETW",
  3, 104.234953, "ETW"
)

test_that("expected diagnostic failures are noted for pooling methods", {
  # This test follows CMH-17-1G Section
  # This section in CMH-17-1G shows the removal of one condition
  # before running Levene's test on the pooled data, so this test
  # will be skipped in this test.
  expect_warning(
    expect_warning(
      expect_warning(
        expect_warning(
          expect_warning(
            res <- basis_pooled_sd(cmh_17_8_3_11_1_1, strength, condition,
                                   batch),
            "outliers_within_batch"
          ),
          "between_group_variability"
        ),
        "outliers_within_group"
      ),
      "pooled_data_normal"
    ),
    "pooled_variance_equal"
  )

  expect_warning(
    expect_warning(
      expect_warning(
        expect_warning(
          res <- cmh_17_8_3_11_1_1 %>%
            filter(condition != "ETW2") %>%
            basis_pooled_sd(strength, condition, batch),
          "outliers_within_batch"
        ),
        "outliers_within_group"
      ),
      "pooled_data_normal"
    ),
    "pooled_variance_equal"
  )

  # removing both ETW and ETW2 should remove all diagnostic failures
  res <- cmh_17_8_3_11_1_1 %>%
    filter(condition != "ETW2" & condition != "ETW") %>%
    basis_pooled_sd(strength, condition, batch)

  expect_equal(res$basis$value[res$basis$group == "CTD"], 108.70,
               tolerance = 0.02)
  expect_equal(res$basis$value[res$basis$group == "RTD"], 88.52,
               tolerance = 0.02)
  expect_equal(res$basis$value[res$basis$group == "ETD"], 80.68,
               tolerance = 0.02)

  expect_warning(
    expect_warning(
      expect_warning(
        expect_warning(
          expect_warning(
            res <- basis_pooled_cv(cmh_17_8_3_11_1_1, strength,
                                   condition, batch),
            "outliers_within_batch"
          ),
          "between_group_variability"
        ),
        "outliers_within_group"
      ),
      "pooled_data_normal"
    ),
    "normalized_variance_equal"
  )

  expect_warning(
    expect_warning(
      expect_warning(
        expect_warning(
          res <- cmh_17_8_3_11_1_1 %>%
            filter(condition != "ETW2") %>%
            basis_pooled_cv(strength, condition, batch),
          "outliers_within_batch"
        ),
        "outliers_within_group"
      ),
      "pooled_data_normal"
    ),
    "normalized_variance_equal"
  )

  # removing both ETW and ETW2 should remove all diagnostic failures
  res <- cmh_17_8_3_11_1_1 %>%
    filter(condition != "ETW2" & condition != "ETW") %>%
    basis_pooled_cv(strength, condition, batch)
})

# data from CMH-17-1G Section 8.3.11.1.2
poolable_data <- tribble(
  ~batch, ~strength, ~condition,
  1, 79.04517, "CTD",
  1, 102.6014, "CTD",
  1, 97.79372, "CTD",
  1, 92.86423, "CTD",
  1, 117.218,  "CTD",
  1, 108.7168, "CTD",
  1, 112.2773, "CTD",
  1, 114.0129, "CTD",
  2, 106.8452, "CTD",
  2, 112.3911, "CTD",
  2, 115.5658, "CTD",
  2, 87.40657, "CTD",
  2, 102.2785, "CTD",
  2, 110.6073, "CTD",
  3, 105.2762, "CTD",
  3, 110.8924, "CTD",
  3, 108.7638, "CTD",
  3, 110.9833, "CTD",
  3, 101.3417, "CTD",
  3, 100.0251, "CTD",
  1, 103.2006, "RTD",
  1, 105.1034, "RTD",
  1, 105.1893, "RTD",
  1, 100.4189, "RTD",
  2, 85.32319, "RTD",
  2, 92.69923, "RTD",
  2, 98.45242, "RTD",
  2, 104.1014, "RTD",
  2, 91.51841, "RTD",
  2, 101.3746, "RTD",
  2, 101.5828, "RTD",
  2, 99.57384, "RTD",
  2, 88.84826, "RTD",
  3, 92.18703, "RTD",
  3, 101.8234, "RTD",
  3, 97.68909, "RTD",
  3, 101.5172, "RTD",
  3, 100.0481, "RTD",
  3, 102.0544, "RTD",
  1, 63.22764, "ETW",
  1, 70.84454, "ETW",
  1, 66.43223, "ETW",
  1, 75.37771, "ETW",
  1, 72.43773, "ETW",
  1, 68.43073, "ETW",
  1, 69.72524, "ETW",
  2, 66.20343, "ETW",
  2, 60.51251, "ETW",
  2, 65.69334, "ETW",
  2, 62.73595, "ETW",
  2, 59.00798, "ETW",
  2, 62.37761, "ETW",
  3, 64.3947,  "ETW",
  3, 72.8491,  "ETW",
  3, 66.56226, "ETW",
  3, 66.56779, "ETW",
  3, 66.00123, "ETW",
  3, 59.62108, "ETW",
  3, 60.61167, "ETW",
  3, 57.65487, "ETW",
  3, 66.51241, "ETW",
  3, 64.89347, "ETW",
  3, 57.73054, "ETW",
  3, 68.94086, "ETW",
  3, 61.63177, "ETW",
  1, 54.09806, "ETW2",
  1, 58.87615, "ETW2",
  1, 61.60167, "ETW2",
  1, 60.23973, "ETW2",
  1, 61.4808,  "ETW2",
  1, 64.55832, "ETW2",
  2, 57.76131, "ETW2",
  2, 49.91463, "ETW2",
  2, 61.49271, "ETW2",
  2, 57.7281,  "ETW2",
  2, 62.11653, "ETW2",
  2, 62.69353, "ETW2",
  3, 61.38523, "ETW2",
  3, 60.39053, "ETW2",
  3, 59.17616, "ETW2",
  3, 60.17616, "ETW2",
  3, 46.47396, "ETW2",
  3, 51.16616, "ETW2"
)

test_that("Pooled SD results match ASAP results", {
  # This data fails the anderson-darling test for normality for the
  # transformed data
  expect_warning(
    expect_message(
      expect_message(
        res_b <- basis_pooled_sd(poolable_data, strength, condition,
                                 override = c("pooled_variance_equal")),
        "outliers_within_batch"
      ),
      "between_group_variability"
    ),
    "pooled_data_normal"
  )

  expect_equal(res_b$basis$value[res_b$basis$group == "CTD"],
               93.64, tolerance = 0.01)
  expect_equal(res_b$basis$value[res_b$basis$group == "RTD"],
               87.30, tolerance = 0.01)
  expect_equal(res_b$basis$value[res_b$basis$group == "ETW"],
               54.33, tolerance = 0.01)
  expect_equal(res_b$basis$value[res_b$basis$group == "ETW2"],
               47.12, tolerance = 0.01)

  expect_equal(res_b$n, 83)
  expect_equal(res_b$r, 4)
  expect_output(print(res_b), "b-basis", ignore.case = TRUE)
  expect_output(print(res_b), "pooled standard deviation", ignore.case = TRUE)
  expect_output(print(res_b), "CTD.*93\\.6", ignore.case = TRUE)
  expect_output(print(res_b), "RTD.*87\\.29", ignore.case = TRUE)
  expect_output(print(res_b), "ETW.*54\\.3", ignore.case = TRUE)
  expect_output(print(res_b), "ETW2.*47\\.07", ignore.case = TRUE)

  res_a <- basis_pooled_sd(poolable_data, strength, condition,
                           p = 0.99, conf = 0.95,
                           override = c("pooled_data_normal",
                                        "pooled_variance_equal",
                                        "outliers_within_batch",
                                        "between_group_variability"))
  expect_equal(res_a$basis$value[res_a$basis$group == "CTD"],
               86.19, tolerance = 0.01)
  expect_equal(res_a$basis$value[res_a$basis$group == "RTD"],
               79.86, tolerance = 0.01)
  expect_equal(res_a$basis$value[res_a$basis$group == "ETW"],
               46.84, tolerance = 0.01)
  expect_equal(res_a$basis$value[res_a$basis$group == "ETW2"],
               39.69, tolerance = 0.01)
})

test_that("Pooled CV results match CMH17STATS", {
  # This data fails the anderson-darling test for normality for the
  # transformed data
  expect_warning(
    expect_message(
      expect_message(
        res_b <- basis_pooled_cv(poolable_data, strength, condition),
        "outliers_within_batch"
      ),
      "between_group_variability"
    ),
    "pooled_data_normal"
  )

  expect_equal(res_b$basis$value[res_b$basis$group == "CTD"],
               90.89, tolerance = 0.01)
  expect_equal(res_b$basis$value[res_b$basis$group == "RTD"],
               85.37, tolerance = 0.01)
  expect_equal(res_b$basis$value[res_b$basis$group == "ETW"],
               56.79, tolerance = 0.01)
  expect_equal(res_b$basis$value[res_b$basis$group == "ETW2"],
               50.55, tolerance = 0.01)

  expect_equal(res_b$n, 83)
  expect_equal(res_b$r, 4)
  expect_output(print(res_b), "b-basis", ignore.case = TRUE)
  expect_output(print(res_b), "pooled CV", ignore.case = TRUE)
  expect_output(print(res_b), "CTD.*90\\.8", ignore.case = TRUE)
  expect_output(print(res_b), "RTD.*85\\.3", ignore.case = TRUE)
  expect_output(print(res_b), "ETW.*56\\.7", ignore.case = TRUE)
  expect_output(print(res_b), "ETW2.*50\\.5", ignore.case = TRUE)

  res_a <- basis_pooled_cv(poolable_data, strength, condition,
                           p = 0.99, conf = 0.95,
                           override = c("pooled_data_normal",
                                        "outliers_within_batch",
                                        "between_group_variability"))
  expect_equal(res_a$basis$value[res_a$basis$group == "CTD"],
               81.62, tolerance = 0.01)
  expect_equal(res_a$basis$value[res_a$basis$group == "RTD"],
               76.67, tolerance = 0.01)
  expect_equal(res_a$basis$value[res_a$basis$group == "ETW"],
               50.98, tolerance = 0.01)
  expect_equal(res_a$basis$value[res_a$basis$group == "ETW2"],
               45.40, tolerance = 0.01)
  expect_output(print(res_a), "a-basis", ignore.case = TRUE)
})

test_that("Pooled data matches CMH17-STATS with mod CV, SD pooling", {
  # pooled SD modified CV results
  # pooled data fails Levene's test after mod CV transform
  # based on `poolable_data` dataset with ETW2 removed

  data <- filter(poolable_data, condition != "ETW2")

  res_b <- basis_pooled_sd(data, strength, condition, modcv = TRUE,
                           override = c("pooled_variance_equal",
                                        "outliers_within_batch",
                                        "between_group_variability"))

  expect_equal(res_b$basis$value[res_b$basis$group == "CTD"],
               92.25, tolerance = 0.01)
  expect_equal(res_b$basis$value[res_b$basis$group == "RTD"],
               85.91, tolerance = 0.01)
  expect_equal(res_b$basis$value[res_b$basis$group == "ETW"],
               52.97, tolerance = 0.01)

  expect_output(print(res_b), "Modified CV")

  res_a <- basis_pooled_sd(data, strength, condition,
                           p = 0.99, conf = 0.95, modcv = TRUE,
                           override = c("pooled_variance_equal",
                                        "outliers_within_batch",
                                        "between_group_variability"))
  expect_equal(res_a$basis$value[res_a$basis$group == "CTD"],
               83.81, tolerance = 0.01)
  expect_equal(res_a$basis$value[res_a$basis$group == "RTD"],
               77.48, tolerance = 0.01)
  expect_equal(res_a$basis$value[res_a$basis$group == "ETW"],
               44.47, tolerance = 0.01)
})

test_that("Pooled data matches CMH17-STATS with mod CV, CV pooling", {
  # pooled CV modified CV results
  # pooled data passes Levene's test after mod CV transform
  # based on `poolable_data` dataset with ETW2 removed

  data <- filter(poolable_data, condition != "ETW2")

  res_b <- basis_pooled_cv(data, strength, condition, modcv = TRUE,
                           override = c("outliers_within_batch",
                                        "between_group_variability"))

  expect_equal(res_b$basis$value[res_b$basis$group == "CTD"],
               90.31, tolerance = 0.01)
  expect_equal(res_b$basis$value[res_b$basis$group == "RTD"],
               84.83, tolerance = 0.01)
  expect_equal(res_b$basis$value[res_b$basis$group == "ETW"],
               56.43, tolerance = 0.01)

  expect_output(print(res_b), "Modified CV")

  res_a <- basis_pooled_cv(data, strength, condition,
                           p = 0.99, conf = 0.95, modcv = TRUE,
                           override = c("outliers_within_batch",
                                        "between_group_variability"))
  expect_equal(res_a$basis$value[res_a$basis$group == "CTD"],
               80.57, tolerance = 0.01)
  expect_equal(res_a$basis$value[res_a$basis$group == "RTD"],
               75.69, tolerance = 0.01)
  expect_equal(res_a$basis$value[res_a$basis$group == "ETW"],
               50.33, tolerance = 0.01)
})

vangel1994 <- tribble(
  ~n, ~z, ~p, ~conf,
  3,  13.976451, 0.90, 0.90,
  5,  4.1011886, 0.90, 0.90,
  7,  2.5440993, 0.90, 0.90,
  9,  1.9368858, 0.90, 0.90,
  11, 1.6127559, 0.90, 0.90,
  13, 1.4096961, 0.90, 0.90,
  15, 1.2695586, 0.90, 0.90,
  17, 1.1663923, 0.90, 0.90,
  19, 1.0868640, 0.90, 0.90,
  21, 1.0234110, 0.90, 0.90,
  3,  28.820048, 0.90, 0.95,
  5,  6.1981307, 0.90, 0.95,
  7,  3.4780112, 0.90, 0.95,
  9,  2.5168762, 0.90, 0.95,
  11, 2.0312134, 0.90, 0.95,
  13, 1.7377374, 0.90, 0.95,
  15, 1.5403989, 0.90, 0.95,
  17, 1.3979806, 0.90, 0.95,
  19, 1.2899172, 0.90, 0.95,
  21, 1.2048089, 0.90, 0.95,
  23, 1.1358259, 0.90, 0.95,
  25, 1.0786237, 0.90, 0.95,
  27, 1.0303046, 0.90, 0.95,
  3,  147.51275, 0.90, 0.99,
  5,  14.993461, 0.90, 0.99,
  7,  6.6442464, 0.90, 0.99,
  9,  4.2798170, 0.90, 0.99,
  11, 3.2197376, 0.90, 0.99,
  13, 2.6267547, 0.90, 0.99,
  15, 2.2493289, 0.90, 0.99,
  17, 1.9880239, 0.90, 0.99,
  19, 1.7961467, 0.90, 0.99,
  21, 1.6490109, 0.90, 0.99,
  23, 1.5323809, 0.90, 0.99,
  25, 1.4374854, 0.90, 0.99,
  27, 1.3586292, 0.90, 0.99,
  29, 1.2919549, 0.90, 0.99,
  31, 1.2347570, 0.90, 0.99,
  33, 1.1850813, 0.90, 0.99,
  35, 1.1414809, 0.90, 0.99,
  37, 1.1028613, 0.90, 0.99,
  39, 1.0683787, 0.90, 0.99,
  41, 1.0373720, 0.90, 0.99,
  43, 1.0093159, 0.90, 0.99,
  3,  20.478521, 0.95, 0.90,
  5,  5.8872014, 0.95, 0.90,
  7,  3.6322326, 0.95, 0.90,
  9,  2.7593956, 0.95, 0.90,
  11, 2.2953853, 0.95, 0.90,
  13, 2.0054547, 0.95, 0.90,
  15, 1.8057261, 0.95, 0.90,
  17, 1.6588820, 0.95, 0.90,
  19, 1.5457939, 0.95, 0.90,
  21, 1.4556317, 0.95, 0.90,
  23, 1.3817937, 0.95, 0.90,
  25, 1.3200198, 0.95, 0.90,
  27, 1.2674334, 0.95, 0.90,
  29, 1.2220187, 0.95, 0.90,
  31, 1.1823195, 0.95, 0.90,
  33, 1.1472560, 0.95, 0.90,
  35, 1.1160097, 0.95, 0.90,
  37, 1.0879479, 0.95, 0.90,
  39, 1.0625739, 0.95, 0.90,
  41, 1.0394913, 0.95, 0.90,
  43, 1.0183802, 0.95, 0.90,
  3,  42.149579, 0.95, 0.95,
  5,  8.8719351, 0.95, 0.95,
  7,  4.9501721, 0.95, 0.95,
  9,  3.5743714, 0.95, 0.95,
  11, 2.8819079, 0.95, 0.95,
  13, 2.4645176, 0.95, 0.95,
  15, 2.1843450, 0.95, 0.95,
  17, 1.9824011, 0.95, 0.95,
  19, 1.8293163, 0.95, 0.95,
  21, 1.7088376, 0.95, 0.95,
  23, 1.6112408, 0.95, 0.95,
  25, 1.5303474, 0.95, 0.95,
  27, 1.4620403, 0.95, 0.95,
  29, 1.4034674, 0.95, 0.95,
  31, 1.3525889, 0.95, 0.95,
  33, 1.3079057, 0.95, 0.95,
  35, 1.2682903, 0.95, 0.95,
  37, 1.2328780, 0.95, 0.95,
  39, 1.2009936, 0.95, 0.95,
  41, 1.1721022, 0.95, 0.95,
  43, 1.1457739, 0.95, 0.95,
  45, 1.1216591, 0.95, 0.95,
  47, 1.0994706, 0.95, 0.95,
  49, 1.0789699, 0.95, 0.95,
  51, 1.0599573, 0.95, 0.95,
  53, 1.0422642, 0.95, 0.95,
  55, 1.0257472, 0.95, 0.95,
  57, 1.0102836, 0.95, 0.95
)

test_that("Extended Hanson-Koopman matches median results from Vangel 1994", {
  # Vangel (1994) provides extensive tables of z for the case where i=1 and
  # j is the median observation. This test checks the results of this
  # package's function against those tables. Only the odd values of n
  # are checked so that the median is a single observation.

  vangel1994 %>%
    rowwise() %>%
    mutate(
      z_calc = hk_ext_z(n, 1, ceiling(n / 2), p, conf)
    ) %>%
    mutate(expect_equal(z, z_calc, tolerance = 0.00005,
                        label = paste0("Mismatch in `z` for n=", n,
                                       ", p=", p,
                                       " conf=", conf, ".\n",
                                       "z_vangel=", z,
                                       ", z_calc=", z_calc, "\n")))
})

cmh_17_1g_8_5_14 <- tribble(
  ~n, ~r, ~k,
  2, 2, 35.177,
  3, 3, 7.859,
  4, 4, 4.505,
  5, 4, 4.101,
  6, 5, 3.064,
  7, 5, 2.858,
  8, 6, 2.382,
  9, 6, 2.253,
  10, 6, 2.137,
  11, 7, 1.897,
  12, 7, 1.814,
  13, 7, 1.738,
  14, 8, 1.599,
  15, 8, 1.540,
  16, 8, 1.485,
  17, 8, 1.434,
  18, 9, 1.354,
  19, 9, 1.311,
  20, 10, 1.253,
  21, 10, 1.218,
  22, 10, 1.184,
  23, 11, 1.143,
  24, 11, 1.114,
  25, 11, 1.087,
  26, 11, 1.060,
  27, 11, 1.035,
  28, 12, 1.010
)

test_that("Extended HK matches CMH-17-1G Table 8.5.14", {
  # CMH-17-1G uses the optimal order statistic approach suggested by
  # Vangel (1994) for computing B-Basis values. There are a few values
  # of n where this package's implementation finds a different optimum order
  # statistic than CMH-17-1G uses. In these cases, the order statistic that
  # this package and CMH-17-1G are both very nearly optimal. These differences
  # are ignored in this test.

  cmh_17_1g_8_5_14 %>%
    rowwise() %>%
    mutate(z = hk_ext_z_j_opt(n, 0.90, 0.95)$z) %>%
    mutate(j = hk_ext_z_j_opt(n, 0.90, 0.95)$j) %>%
    filter(
      n != 17 & n != 20 & n != 23 & n != 24 & n != 28
      ) %>%
    mutate(expect_equal(j, r,
                        label = paste0("Mismatch in `j`/`r` for n=", n, ", ",
                                       "r_B_cmh=", r,
                                       ", j=", j, "\n"))) %>%
    mutate(expect_equal(z, k, tolerance = 0.005,
                        label = paste0("Mismatch in `k`/`z` for n=", n, ", ",
                                       "k_B_cmh=", k,
                                       ", z_calc=", z, "\n")))
})

test_that("Hanson-Koopman results match STAT17 for several values of n", {
  data <- c(
    139.6734,
    143.0032,
    130.4757,
    144.8327,
    138.7818,
    136.7693,
    148.636,
    131.0095,
    131.4933,
    142.8856,
    158.0198,
    145.2271,
    137.5991,
    139.8298,
    140.8557,
    137.6148,
    131.3614,
    152.7795,
    145.8792,
    152.9207,
    160.0989,
    145.192,
    128.6383,
    141.5992,
    122.5297,
    159.8209,
    151.672,
    159.0156
  )

  res <- basis_hk_ext(x = head(data, 28), p = 0.9, conf = 0.95,
                      method = "optimum-order", override = "all")
  expect_equal(res$basis, 122.36798, tolerance = 0.001)

  res <- basis_hk_ext(x = head(data, 27), p = 0.9, conf = 0.95,
                      method = "optimum-order", override = "all")
  expect_equal(res$basis, 121.96939, tolerance = 0.001)

  res <- basis_hk_ext(x = head(data, 26), p = 0.9, conf = 0.95,
                      method = "optimum-order", override = "all")
  expect_equal(res$basis, 121.57073, tolerance = 0.001)

  res <- basis_hk_ext(x = head(data, 23), p = 0.9, conf = 0.95,
                      method = "optimum-order", override = "all")
  expect_equal(res$basis, 127.11286, tolerance = 0.001)

  res <- basis_hk_ext(x = head(data, 22), p = 0.9, conf = 0.95,
                      method = "optimum-order", override = "all")
  expect_equal(res$basis, 128.82397, tolerance = 0.001)

  res <- basis_hk_ext(x = head(data, 21), p = 0.9, conf = 0.95,
                      method = "optimum-order", override = "all")
  expect_equal(res$basis, 128.52107, tolerance = 0.001)

  res <- basis_hk_ext(x = head(data, 20), p = 0.9, conf = 0.95,
                      method = "optimum-order", override = "all")
  expect_equal(res$basis, 128.20999, tolerance = 0.001)

  res <- basis_hk_ext(x = head(data, 19), p = 0.9, conf = 0.95,
                      method = "optimum-order", override = "all")
  expect_equal(res$basis, 127.74060, tolerance = 0.002)

  res <- basis_hk_ext(x = head(data, 18), p = 0.9, conf = 0.95,
                      method = "optimum-order", override = "all")
  expect_equal(res$basis, 127.36697, tolerance = 0.001)

  res <- basis_hk_ext(x = head(data, 17), p = 0.9, conf = 0.95,
                      method = "optimum-order", override = "all")
  expect_equal(res$basis, 127.02732, tolerance = 0.001)

  res <- basis_hk_ext(x = head(data, 16), p = 0.9, conf = 0.95,
                      method = "optimum-order", override = "all")
  expect_equal(res$basis, 126.23545, tolerance = 0.002)

  res <- basis_hk_ext(x = head(data, 15), p = 0.9, conf = 0.95,
                      method = "optimum-order", override = "all")
  expect_equal(res$basis, 125.68740, tolerance = 0.001)

  res <- basis_hk_ext(x = head(data, 14), p = 0.9, conf = 0.95,
                      method = "optimum-order", override = "all")
  expect_equal(res$basis, 125.17500, tolerance = 0.001)

  res <- basis_hk_ext(x = head(data, 13), p = 0.9, conf = 0.95,
                      method = "optimum-order", override = "all")
  expect_equal(res$basis, 124.07851, tolerance = 0.002)

  res <- basis_hk_ext(x = head(data, 12), p = 0.9, conf = 0.95,
                      method = "optimum-order", override = "all")
  expect_equal(res$basis, 121.17418, tolerance = 0.001)

  res <- basis_hk_ext(x = head(data, 11), p = 0.9, conf = 0.95,
                      method = "optimum-order", override = "all")
  expect_equal(res$basis, 120.26382, tolerance = 0.001)

  res <- basis_hk_ext(x = head(data, 10), p = 0.9, conf = 0.95,
                      method = "optimum-order", override = "all")
  expect_equal(res$basis, 120.75149, tolerance = 0.002)

  res <- basis_hk_ext(x = head(data, 9), p = 0.9, conf = 0.95,
                      method = "optimum-order", override = "all")
  expect_equal(res$basis, 119.80108, tolerance = 0.001)
})

cmh_17_1g_8_5_15 <- tribble(
  ~n, ~k,
  2, 80.0038,
  3, 16.9122,
  4, 9.49579,
  5, 6.89049,
  6, 5.57681,
  7, 4.78352,
  8, 4.25011,
  9, 3.86502,
  10, 3.57267,
  11, 3.34227,
  12, 3.1554,
  13, 3.00033,
  14, 2.86924,
  15, 2.75672,
  16, 2.65889,
  17, 2.5729,
  18, 2.4966,
  19, 2.42833,
  20, 2.36683,
  21, 2.31106,
  22, 2.2602,
  23, 2.21359,
  24, 2.17067,
  25, 2.131,
  26, 2.09419,
  27, 2.05991,
  28, 2.0279,
  29, 1.99791,
  30, 1.96975,
  31, 1.94324,
  32, 1.91822,
  33, 1.89457,
  34, 1.87215,
  35, 1.85088,
  36, 1.83065,
  37, 1.81139,
  38, 1.79301,
  39, 1.77546,
  40, 1.75868,
  41, 1.7426,
  42, 1.72718,
  43, 1.71239,
  44, 1.69817,
  45, 1.68449,
  46, 1.67132,
  47, 1.65862,
  48, 1.64638,
  49, 1.63456,
  50, 1.62313,
  52, 1.60139,
  54, 1.58101,
  56, 1.56184,
  58, 1.54377,
  60, 1.5267,
  62, 1.51053,
  64, 1.4952,
  66, 1.48063,
  68, 1.46675,
  70, 1.45352,
  72, 1.44089,
  74, 1.42881,
  76, 1.41724,
  78, 1.40614,
  80, 1.39549,
  82, 1.38525,
  84, 1.37541,
  86, 1.36592,
  88, 1.35678,
  90, 1.34796,
  92, 1.33944,
  94, 1.3312,
  96, 1.32324,
  98, 1.31553,
  100, 1.30806,
  105, 1.29036,
  110, 1.27392,
  115, 1.25859,
  120, 1.24425,
  125, 1.2308,
  130, 1.21814,
  135, 1.2062,
  140, 1.19491,
  145, 1.18421,
  150, 1.17406,
  155, 1.1644,
  160, 1.15519,
  165, 1.1464,
  170, 1.13801,
  175, 1.12997,
  180, 1.12226,
  185, 1.11486,
  190, 1.10776,
  195, 1.10092,
  200, 1.09434,
  205, 1.08799,
  210, 1.08187,
  215, 1.07595,
  220, 1.07024,
  225, 1.06471,
  230, 1.05935,
  235, 1.05417,
  240, 1.04914,
  245, 1.04426,
  250, 1.03952,
  275, 1.01773
)

test_that("Extended Hanson-Koopman matches CMH-17-1G Table 8.5.15", {
  # for A-Basis, CMH-17-1G uses the order statistics for 1 and n
  # to compute the tolerance limits. This test verifies that the code
  # in this package computes the same values of z (k, as CMH-17 calls it)
  cmh_17_1g_8_5_15 %>%
    rowwise() %>%
    mutate(z = hk_ext_z(n, 1, n, 0.99, 0.95)) %>%
    mutate(expect_equal(z, k, tolerance = 0.00005,
                        label = paste0("Mismatch in `k`/`z` for n=", n, ", ",
                                       "k_A_cmh=", k,
                                       ", z_calc=", z, "\n")))
})

cmh_17_1g_8_5_13 <- tribble(
  ~n, ~ra,
  299, 1,
  473, 2,
  628, 3,
  773, 4,
  913, 5,
  1049, 6,
  1182, 7,
  1312, 8,
  1441, 9,
  1568, 10,
  1693, 11,
  1818, 12,
  1941, 13,
  2064, 14,
  2185, 15,
  2306, 16,
  2426, 17,
  2546, 18,
  2665, 19,
  2784, 20,
  2902, 21,
  3020, 22,
  3137, 23,
  3254, 24,
  3371, 25,
  3487, 26,
  3603, 27,
  3719, 28,
  3834, 29,
  3949, 30,
  4064, 31,
  4179, 32,
  4293, 33,
  4407, 34,
  4521, 35,
  4635, 36,
  4749, 37,
  4862, 38,
  4975, 39,
  5088, 40,
  5201, 41,
  5314, 42,
  5427, 43,
  5539, 44,
  5651, 45,
  5764, 46,
  5876, 47,
  5988, 48,
  6099, 49,
  6211, 50,
  6323, 51,
  6434, 52,
  6545, 53,
  6657, 54,
  6769, 55,
  6879, 56,
  6990, 57,
  7100, 58,
  7211, 59,
  7322, 60,
  7432, 61,
  7543, 62,
  7653, 63,
  7763, 64,
  7874, 65,
  7984, 66,
  8094, 67,
  8204, 68,
  8314, 69,
  8423, 70,
  8533, 71,
  8643, 72,
  8753, 73,
  8862, 74,
  8972, 75,
  9081, 76,
  9190, 77,
  9300, 78,
  9409, 79,
  9518, 80,
  9627, 81,
  9736, 82,
  9854, 83,
  9954, 84,
  10063, 85,
  10172, 86,
  10281, 87,
  10390, 88,
  10498, 89,
  10607, 90,
  10716, 91,
  10824, 92,
  10933, 93,
  11041, 94,
  11150, 95,
  11258, 96,
  11366, 97,
  11475, 98,
  11583, 99,
  11691, 100
)

test_that("Non-parametric ranks for A-Basis match CMH-17-1G Table 8.5.13", {
  skip_on_cran()  # this test is a long-running test

  cmh_17_1g_8_5_13 %>%
    mutate(ra_lag = lag(ra)) %>%
    rowwise() %>%
    mutate(r_calc = nonpara_binomial_rank(n, 0.99, 0.95)) %>%
    mutate(expect_equal(ra, r_calc,
                        label = paste0(
                          "Mismatch in r for n=", n,
                          ". rA=", ra,
                          ", r_calc=", r_calc
                          ))) %>%
    filter(n > 299 & n < 6500) %>%
    # the rank for one sample larger should be the same
    mutate(r_calc_plus = nonpara_binomial_rank(n + 1, 0.99, 0.95)) %>%
    mutate(expect_equal(ra, r_calc_plus,
                        label = paste0(
                          "Mismatch in r for n=", n + 1,
                          ". rA=", ra, ", ",
                          "r_calc=", r_calc_plus
                        ))) %>%
    # the rank for one sample smaller should be the previous one
    mutate(r_calc_minus = nonpara_binomial_rank(n - 1, 0.99, 0.95)) %>%
    mutate(expect_equal(ra_lag, r_calc_minus,
                        label = paste0(
                          "Mismatch in r for n=", n - 1,
                          ". rA=", ra_lag, ", ",
                          "r_calc=", r_calc_minus
                        )))
})

test_that("nonpara_binomial_rank raises and error when sample too small", {
  expect_error(nonpara_binomial_rank(298, 0.99, 0.95),
               "p.*0\\.99.*conf.*0\\.95")
})

test_that("nonpara_binomial_rank raises an error when it can't converge", {
  expect_error(nonpara_binomial_rank(4000, 0.00001, 0.01),
               "p.*1e-05.*conf.*0\\.01")
})

cmh_17_1g_8_5_12 <- tribble(
  ~n, ~rb,
  29, 1,
  46, 2,
  61, 3,
  76, 4,
  89, 5,
  103, 6,
  116, 7,
  129, 8,
  142, 9,
  154, 10,
  167, 11,
  179, 12,
  191, 13,
  203, 14,
  215, 15,
  227, 16,
  239, 17,
  251, 18,
  263, 19,
  275, 20,
  298, 22,
  321, 24,
  345, 26,
  368, 28,
  391, 30,
  413, 32,
  436, 34,
  459, 36,
  481, 38,
  504, 40,
  526, 42,
  549, 44,
  571, 46,
  593, 48,
  615, 50,
  638, 52,
  660, 54,
  682, 56,
  704, 58,
  726, 60,
  781, 65,
  836, 70,
  890, 75,
  945, 80,
  999, 85,
  1053, 90,
  1107, 95,
  1161, 100,
  1269, 110,
  1376, 120,
  1483, 130,
  1590, 140,
  1696, 150,
  1803, 160,
  1909, 170,
  2015, 180,
  2120, 190,
  2226, 200,
  2331, 210,
  2437, 220,
  2542, 230,
  2647, 240,
  2752, 250,
  2857, 260,
  2962, 270,
  3066, 280,
  3171, 290,
  3276, 300,
  3380, 310,
  3484, 320,
  3589, 330,
  3693, 340,
  3797, 350,
  3901, 360,
  4005, 370,
  4109, 380,
  4213, 390,
  4317, 400,
  4421, 410,
  4525, 420,
  4629, 430,
  4733, 440,
  4836, 450,
  4940, 460,
  5044, 470,
  5147, 480,
  5251, 490,
  5354, 500,
  5613, 525,
  5871, 550,
  6130, 575,
  6388, 600,
  6645, 625,
  6903, 650,
  7161, 675,
  7418, 700,
  7727, 730,
  8036, 760,
  8344, 790,
  8652, 820,
  8960, 850,
  9268, 880,
  9576, 910,
  9884, 940,
  10191, 970,
  10499, 1000
)

test_that("Non-parametric ranks for BA-Basis match CMH-17-1G Table 8.5.12", {
  skip_on_cran()  # this test is a long-running test

  cmh_17_1g_8_5_12 %>%
    mutate(rb_lag = lag(rb)) %>%
    rowwise() %>%
    mutate(r_calc = nonpara_binomial_rank(n, 0.9, 0.95)) %>%
    mutate(expect_equal(rb, r_calc,
                        label = paste0(
                          "Mismatch in r for n=", n,
                          ". rB=", rb,
                          ", r_calc=", r_calc
                        ))) %>%
    # the rank for one sample larger should be the same
    mutate(r_calc_plus = nonpara_binomial_rank(n + 1, 0.9, 0.95)) %>%
    mutate(expect_equal(rb, r_calc_plus,
                        label = paste0(
                          "Mismatch in r for n=", n + 1,
                          ". rB=", rb, ", ",
                          "r_calc=", r_calc_plus
                        ))) %>%
    filter(n > 29 & n <= 275) %>%
    # the rank for one sample smaller should be the previous one
    # above n=275, Table 8.5.12 does not have consecutive ranks, so we can't
    # use the lag trick below to check sample sizes of n-1
    mutate(r_calc_minus = nonpara_binomial_rank(n - 1, 0.9, 0.95)) %>%
    mutate(expect_equal(rb_lag, r_calc_minus,
                        label = paste0(
                          "Mismatch in r for n=", n - 1,
                          ". rB=", rb_lag, ", ",
                          "r_calc=", r_calc_minus
                        )))
})

cmh_17_1g_8_3_11_1_1_etw2 <- tribble(
  ~batch, ~strength,
  1, 99.0239966,
  1, 103.341238,
  1, 100.30213,
  1, 98.4634133,
  1, 92.264728,
  1, 103.487693,
  1, 113.734763,
  2, 108.172659,
  2, 108.426732,
  2, 116.260375,
  2, 121.04961,
  2, 111.223082,
  2, 104.574843,
  2, 103.222552,
  3, 99.3918538,
  3, 87.3421658,
  3, 102.730741,
  3, 96.3694916,
  3, 99.5946088,
  3, 97.0712407
)

test_that("ANOVA results match STAT17 for sample data", {
  # Sample data from CMH-17-1G Section 8.3.11.2.2

  res <- cmh_17_1g_8_3_11_1_1_etw2 %>%
    basis_anova(strength, batch,
                override = c("equality_of_variance",
                             "number_of_groups"))

  expect_equal(res$basis, 63.2, tolerance = 0.05)
  expect_output(print(res), "b-basis.*63\\.2", ignore.case = TRUE)
  expect_output(print(res), "ANOVA", ignore.case = TRUE)
  expect_match(res$distribution, "ANOVA", ignore.case = TRUE)

  res <- cmh_17_1g_8_3_11_1_1_etw2 %>%
    basis_anova(strength, batch, p = 0.99, conf = 0.95,
                override = c("equality_of_variance",
                             "number_of_groups"))

  expect_equal(res$basis, 34.6, tolerance = 0.05)
  expect_output(print(res), "a-basis.*34\\.", ignore.case = TRUE)
  expect_output(print(res), "ANOVA", ignore.case = TRUE)
  expect_match(res$distribution, "ANOVA", ignore.case = TRUE)
})

test_that("ANOVA produces an error when there is only one group", {
  strength <- rep(10, 1)
  batch <- rep(10, 1)

  expect_error(
    basis_anova(x = strength, group = batch),
    "fewer than 2"
  )
})

test_that("anova basis values produce expected diagnostic failures", {
  set.seed(100)
  x <- c(rnorm(30, 100, 1), rnorm(30, 100, 10), 80)
  batch <- c(rep("A", 30), rep("B", 30), "A")

  expect_warning(
    expect_warning(
      expect_warning(
        res <- basis_anova(x = x, group = batch),
        "outliers_within_group"
      ),
      "equality_of_variance"
    ),
    "number_of_groups"
  )

  # Check that res$... contains the correct value
  expect_equal(res$group, batch)
  expect_equal(res$diagnostic_failures,
               c("outliers_within_group",
                 "equality_of_variance",
                 "number_of_groups"))
  expect_length(res$override, 0)

  # overriding the diagnostics should eliminate the warnings
  res <- basis_anova(x = x, group = batch,
                     override = c("outliers_within_group",
                                  "equality_of_variance",
                                  "number_of_groups"))

  expect_equal(res$override,
               c("outliers_within_group",
                 "equality_of_variance",
                 "number_of_groups"))
  expect_length(res$diagnostic_failures, 0)

  # overriding the diagnostics with "all" should do the same thing
  res <- basis_anova(x = x, group = batch,
                     override = "all")
  expect_equal(res$override,
               c("outliers_within_group",
                 "equality_of_variance",
                 "number_of_groups"))
  expect_length(res$diagnostic_failures, 0)
})

test_that("ANOVA method matches STAT17 when between-batch var. is small", {
  data <- tribble(
    ~x, ~batch,
    105.04953017290813, 1,
    105.74515635546253, 1,
    99.7549396676824, 1,
    107.44219439303261, 1,
    100.17657481474124, 1,
    106.601810738431, 1,
    101.15202811896768, 2,
    90.63466521331704, 2,
    106.93692070778634, 2,
    116.14555531325212, 2,
    100.20555336225114, 2,
    103.89002397699194, 2,
    110.50367678215923, 3,
    95.34690617376182, 3,
    105.03624331633935, 3,
    105.83852344481843, 3,
    105.8785931848096, 3,
    103.97623814685818, 3,
    94.92344509669459, 4,
    89.35739844589054, 4,
    110.45073142288507, 4,
    108.32807015574465, 4,
    104.35498641239826, 4,
    109.39785860273314, 4,
    102.88966425996772, 5,
    105.08208381529616, 5,
    109.82310733067601, 5,
    108.64289487358796, 5,
    99.87084985403291, 5,
    96.7651412720645, 5
  )
  res <- basis_anova(data, x, batch, override = "all")
  expect_equal(res$basis, 93.2, tolerance = 0.05)
})

test_that("glance.basis produces expected value", {
  # Sample data from CMH-17-1G Section 8.3.11.2.2

  res <- cmh_17_1g_8_3_11_1_1_etw2 %>%
    basis_anova(strength, batch,
                override = c("number_of_groups"))

  glance_res <- glance(res)

  expect_equal(glance_res[["p"]][1], 0.9)
  expect_equal(glance_res[["conf"]][1], 0.95)
  expect_equal(glance_res[["distribution"]][1], "ANOVA")
  expect_equal(glance_res[["n"]][1], nrow(cmh_17_1g_8_3_11_1_1_etw2))
  expect_equal(glance_res[["r"]][1], 3)
  expect_equal(glance_res[["basis"]][1], 63.2, tolerance = 0.05)

  glance_res_2 <- glance(res, TRUE)

  for (gn in names(glance_res)) {
    expect_equal(glance_res[[gn]], glance_res_2[[gn]])
  }

  expect_equal(glance_res_2[["outliers_within_group"]], "P")
  expect_equal(glance_res_2[["equality_of_variance"]], "P")
  expect_equal(glance_res_2[["number_of_groups"]], "O")

  expect_warning({
    glance_res_3 <- cmh_17_1g_8_3_11_1_1_etw2 %>%
      basis_anova(strength, batch) %>%
      glance(TRUE)
  })

  expect_equal(glance_res_3[["outliers_within_group"]], "P")
  expect_equal(glance_res_3[["equality_of_variance"]], "P")
  expect_equal(glance_res_3[["number_of_groups"]], "F")
})

test_that("glance for pooled methods works", {
  res <- carbon.fabric %>%
    filter(test == "WT") %>%
    basis_pooled_sd(strength, condition, batch,
                    override = c("outliers_within_batch")) %>%
    glance(TRUE)

  # 3 conditions should produce 3 basis values and hence 3 rows
  expect_equal(nrow(res), 3)
})

test_that("pooled methods process override='all'", {

  res <- basis_pooled_sd(poolable_data, strength, condition, modcv = TRUE,
                         override = "all")
  expect_equal(res$override,
               c("outliers_within_batch",
                 "between_group_variability",
                 "outliers_within_group",
                 "pooled_data_normal",
                 "pooled_variance_equal"))
  expect_length(res$diagnostic_failures, 0)

  res <- basis_pooled_cv(poolable_data, strength, condition, modcv = TRUE,
                         override = "all")
  expect_equal(res$override,
               c("outliers_within_batch",
                 "between_group_variability",
                 "outliers_within_group",
                 "pooled_data_normal",
                 "normalized_variance_equal"))
  expect_length(res$diagnostic_failures, 0)
})