tests/testthat/test-basis-factor-hk.R

suppressMessages(library(dplyr))

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")))
})

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")))
})

Try the cmstatr package in your browser

Any scripts or data that you put into this service are public.

cmstatr documentation built on April 4, 2025, 1:46 a.m.