tests/testthat/test-htest.R

In TOSTER: Two One-Sided Tests (TOST) Equivalence Testing

# All htest extensions tested here

test_that("simple_htest: t-test & wilcox", {

  set.seed(3164964)

  expect_equal(t.test(1:10, y = c(7:20))$p.value,
               simple_htest(1:10, y = c(7:20))$p.value)
  expect_equal(t.test(1:10, y = c(7:20, 200))$p.value,
               simple_htest(1:10, y = c(7:20, 200))$p.value)

  testy1 = describe_htest(simple_htest(1:10, y = c(7:20, 200)))
  testy2 = describe_htest(t.test(1:10, y = c(7:20, 200)))
  expect_equal(testy1,testy2)

  testy1 = df_htest(simple_htest(1:10, y = c(7:20, 200)))
  testy2 = df_htest(t.test(1:10, y = c(7:20, 200)))
  expect_equal(testy1,testy2)

  testy1 = describe_htest(as_htest(t_TOST(1:10, y = c(7:20, 200), eqb = 3)))
  testy1 = describe_htest(as_htest(t_TOST(1:10, y = c(7:20, 200), eqb = 3,
                                          hypothesis = "MET")))

  # wilcox -- same data
  expect_equal(wilcox.test(1:10, y = c(7:20))$p.value,
               simple_htest(1:10, y = c(7:20), test = "w")$p.value)
  expect_equal(wilcox.test(1:10, y = c(7:20, 200))$p.value,
               simple_htest(1:10, y = c(7:20, 200), test = "w")$p.value)

  testy1 = describe_htest(simple_htest(1:10, y = c(7:20, 200), test = "w"))
  testy2 = describe_htest(wilcox.test(1:10, y = c(7:20, 200),
                                      conf.int = TRUE))
  expect_equal(testy1,testy2)

  testy1 = df_htest(simple_htest(1:10, y = c(7:20, 200), test = "w"))
  testy2 = df_htest(wilcox.test(1:10, y = c(7:20, 200),
                                      conf.int = TRUE))
  expect_equal(testy1,testy2)

  testy1 = describe_htest(as_htest(wilcox_TOST(1:10, y = c(7:20, 200), eqb = 3)))
  testy1 = describe_htest(as_htest(wilcox_TOST(1:10, y = c(7:20, 200), eqb = 3,
                                               hypothesis = "MET")))

  expect_error(simple_htest(1:10, y = c(7:20, 200), test = "w",
                            alternative = "e"))

  # Test equivalence -----
  testy1 = simple_htest(data = mtcars,
               mpg ~ am, test = "w",
               mu = 3,
               alternative = "e")
  testy2 = wilcox.test(data = mtcars,
                        mpg ~ am,
                        mu = -3,
                        alternative = "g")
  expect_equal(testy1$p.value,testy2$p.value)
  testy1 = simple_htest(data = mtcars,
                        mpg ~ am, test = "t",
                        mu = 3,
                        alternative = "e")
  testy2 = t.test(data = mtcars,
                       mpg ~ am,
                       mu = -3,
                       alternative = "g")
  expect_equal(testy1$p.value,testy2$p.value)


  # Test MET -----
  testy1 = simple_htest(data = mtcars,
                        mpg ~ am, test = "w",
                        mu = 3,
                        alternative = "m")
  testy2 = wilcox.test(data = mtcars,
                       mpg ~ am,
                       mu = -3,
                       alternative = "l")
  testydf = df_htest(testy1)
  expect_equal(testy1$p.value,testy2$p.value)
  expect_equal(testy1$p.value, testydf$p.value)
  testy1 = simple_htest(data = mtcars,
                        mpg ~ am, test = "t",
                        mu = 3,
                        alternative = "m")
  testy2 = t.test(data = mtcars,
                  mpg ~ am,
                  mu = -3,
                  alternative = "l")
  expect_equal(testy1$p.value,testy2$p.value)

  # Errors & Messages

  expect_error(df_htest(1),"htest must be of the class htest")
  expect_error(describe_htest(1),"htest must be of the class htest")
  testy2 = t.test(data = mtcars,
                  mpg ~ am,
                  mu = -3,
                  alternative = "l")
  testy2$p.value = NULL
  expect_error(describe_htest(testy2))
  ow = stats::oneway.test(extra ~ group, data = sleep)
  expect_message(describe_htest(ow))
  wil_noci =  wilcox.test(data = mtcars,
                                    mpg ~ am,
                                    mu = -3,
                                    alternative = "l")
  expect_message(describe_htest(wil_noci))

  expect_equal(df_htest(wil_noci, extract_names = FALSE)$statistic,
              73)

  expect_equal(df_htest(testy2, show_ci = FALSE)$conf.level,
               NULL)
  expect_equal(df_htest(testy2, test_statistics = FALSE)$t,
               NULL)
})

test_that("brunner_munzel",{

  set.seed(2808)

  expect_equal(brunner_munzel(1:10, y = c(7:20))$p.value,
               simple_htest(1:10, y = c(7:20),
                            test = "b",
                            mu = .5)$p.value)
  expect_equal(brunner_munzel(1:10, y = c(7:20, 200))$p.value,
               simple_htest(1:10, y = c(7:20, 200),
                            test = "b",
                            mu = .5)$p.value)

  testy1 = describe_htest(simple_htest(1:10, y = c(7:20, 200),
                                       test = "b",
                                       mu = .5))
  testy2 = describe_htest(brunner_munzel(1:10, y = c(7:20, 200)))
  expect_equal(testy1,testy2)

  testy1 = df_htest(simple_htest(1:10, y = c(7:20, 200),
                                 test = "b",
                                 mu = .5))
  testy2 = df_htest(brunner_munzel(1:10, y = c(7:20, 200)))
  expect_equal(testy1,testy2)

  testy1 = simple_htest(data = mtcars,
                        mpg ~ am,
                        test = "b",
                        mu = .75,
                        alternative = "e")
  testy2 = brunner_munzel(data = mtcars,
                       mpg ~ am,
                       mu = .25,
                       alternative = "g")
  expect_equal(testy1$p.value,testy2$p.value)
  set.seed(1944)
  testy1 = simple_htest(data = mtcars,
                        mpg ~ am,
                        test = "b",
                        mu = .25,
                        alternative = "g",
                        perm = TRUE)
  set.seed(1944)
  testy2 = brunner_munzel(data = mtcars,
                          mpg ~ am,
                          mu = .25,
                          alternative = "g",
                          perm = TRUE)
  expect_equal(testy1$p.value,testy2$p.value)

  set.seed(1945)
  testy1 = simple_htest(data = mtcars,
                        mpg ~ am,
                        test = "b",
                        mu = .25,
                        alternative = "l",
                        perm = TRUE)
  set.seed(1945)
  testy2 = brunner_munzel(data = mtcars,
                          mpg ~ am,
                          mu = .25,
                          alternative = "l",
                          perm = TRUE)
  expect_equal(testy1$p.value,testy2$p.value)

  set.seed(1946)
  testy1 = simple_htest(data = mtcars,
                        mpg ~ am,
                        test = "b",
                        mu = .5,
                        alternative = "t",
                        perm = TRUE)
  set.seed(1946)
  testy2 = brunner_munzel(data = mtcars,
                          mpg ~ am,
                          mu = .5,
                          alternative = "t",
                          perm = TRUE)
  expect_equal(testy1$p.value,testy2$p.value)

  # Test equivalence -----
  testy1 = simple_htest(data = mtcars,
                        mpg ~ am, test = "b",
                        mu = .4,
                        alternative = "e")
  testy2 = brunner_munzel(data = mtcars,
                       mpg ~ am,
                       mu = .4,
                       alternative = "g")
  expect_equal(testy1$p.value,testy2$p.value)
  testy1 = simple_htest(data = mtcars,
                        mpg ~ am, test = "b",
                        mu = .3,
                        alternative = "e")
  testy2 = brunner_munzel(data = mtcars,
                  mpg ~ am,
                  mu = .3,
                  alternative = "g")
  expect_equal(testy1$p.value,testy2$p.value)


  # Test MET -----
  testy1 = simple_htest(data = mtcars,
                        mpg ~ am, test = "b",
                        mu = .3,
                        alternative = "m")
  testy2 = brunner_munzel(data = mtcars,
                       mpg ~ am,
                       mu = .3,
                       alternative = "l")
  testydf = df_htest(testy1)
  expect_equal(testy1$p.value,testy2$p.value)
  expect_equal(testy1$p.value, testydf$p.value)

  # Paired ------
  data(sleep)

  set.seed(1944)
  testy1 = simple_htest(data = sleep,
                        extra ~ group,
                        test = "b",
                        mu = .25,
                        alternative = "g",
                        perm = TRUE,
                        paired = TRUE)
  set.seed(1944)
  testy2 = brunner_munzel(data = sleep,
                          extra ~ group,
                          mu = .25,
                          alternative = "g",
                          perm = TRUE,
                          paired = TRUE)
  expect_equal(testy1$p.value,testy2$p.value)

  test_big = brunner_munzel(x = rnorm(100),
                            y = rnorm(100),
                          mu = .25,
                          alternative = "t",
                          perm = TRUE,
                          paired = TRUE)

  set.seed(1945)
  testy1 = simple_htest(data = sleep,
                        extra ~ group,
                        test = "b",
                        mu = .25,
                        alternative = "l",
                        perm = TRUE,
                        paired = TRUE)
  set.seed(1945)
  testy2 = brunner_munzel(data = sleep,
                          extra ~ group,
                          mu = .25,
                          alternative = "l",
                          perm = TRUE,
                          paired = TRUE)
  expect_equal(testy1$p.value,testy2$p.value)

  set.seed(1946)
  testy1 = simple_htest(data = sleep,
                        extra ~ group,
                        test = "b",
                        mu = .5,
                        alternative = "t",
                        perm = TRUE,
                        paired = TRUE)
  set.seed(1946)
  testy2 = brunner_munzel(data = sleep,
                          extra ~ group,
                          mu = .5,
                          alternative = "t",
                          perm = TRUE,
                          paired = TRUE)
  expect_equal(testy1$p.value,testy2$p.value)

  ## Equ ---
  testy1 = simple_htest(data = sleep,
                        extra ~ group, paired = TRUE,
                        test = "b",
                        mu = .3,
                        alternative = "e")
  testy2 = brunner_munzel(data = sleep,
                          extra ~ group,
                          paired = TRUE,
                          mu = .3,
                          alternative = "g")
  expect_equal(testy1$p.value,testy2$p.value)


  ## MET -----
  testy1 = simple_htest(data = sleep, paired = TRUE,
                        extra ~ group, test = "b",
                        mu = .3,
                        alternative = "m")
  testy2 = brunner_munzel(data = sleep,
                          extra ~ group,paired = TRUE,
                          mu = .3,
                          alternative = "l")
  testydf = df_htest(testy1)
  expect_equal(testy1$p.value,testy2$p.value)
  expect_equal(testy1$p.value, testydf$p.value)

  # Errors ----

  expect_error(brunner_munzel(x = rnorm(100),
                              #y = rnorm(100),
                              mu = .25,
                              alternative = "t",
                              perm = TRUE,
                              paired = TRUE))

  expect_error(brunner_munzel(x = "error",
                              #y = rnorm(100),
                              mu = .25,
                              alternative = "t",
                              perm = TRUE,
                              paired = TRUE))

  expect_error(brunner_munzel(x = rnorm(100),
                              y = rnorm(99),
                              mu = .25,
                              alternative = "t",
                              perm = TRUE,
                              paired = TRUE))

  expect_message(brunner_munzel(x = rnorm(300),
                              y = rnorm(300),
                              mu = .25,
                              alternative = "t",
                              perm = TRUE,
                              #paired = TRUE
                              ))



})

test_that("All other htests",{
  # Quade ---

  ## Conover (1999, p. 375f):
  ## Numbers of five brands of a new hand lotion sold in seven stores
  ## during one week.
  y <- matrix(c( 5,  4,  7, 10, 12,
                 1,  3,  1,  0,  2,
                 16, 12, 22, 22, 35,
                 5,  4,  3,  5,  4,
                 10,  9,  7, 13, 10,
                 19, 18, 28, 37, 58,
                 10,  7,  6,  8,  7),
              nrow = 7, byrow = TRUE,
              dimnames =
                list(Store = as.character(1:7),
                     Brand = LETTERS[1:5]))

  htest <- stats::quade.test(y)
  df1 = df_htest(htest = htest)
  pr = describe_htest(htest, alpha = .05)
  expect_equal(htest$p.value,df1$p.value)

  # McNemar ----

  ## Agresti (1990), p. 350.
  ## Presidential Approval Ratings.
  ##  Approval of the President's performance in office in two surveys,
  ##  one month apart, for a random sample of 1600 voting-age Americans.
  Performance <-
    matrix(c(794, 86, 150, 570),
           nrow = 2,
           dimnames = list("1st Survey" = c("Approve", "Disapprove"),
                           "2nd Survey" = c("Approve", "Disapprove")))

  htest = stats::mcnemar.test(Performance)
  df1 = df_htest(htest)
  pr = describe_htest(htest, alpha = .05)
  expect_equal(htest$p.value,df1$p.value)

  # Friedman ---
  ## Hollander & Wolfe (1973), p. 140ff.
  ## Comparison of three methods ("round out", "narrow angle", and
  ##  "wide angle") for rounding first base.  For each of 18 players
  ##  and the three method, the average time of two runs from a point on
  ##  the first base line 35ft from home plate to a point 15ft short of
  ##  second base is recorded.
  RoundingTimes <-
    matrix(c(5.40, 5.50, 5.55,
             5.85, 5.70, 5.75,
             5.20, 5.60, 5.50,
             5.55, 5.50, 5.40,
             5.90, 5.85, 5.70,
             5.45, 5.55, 5.60,
             5.40, 5.40, 5.35,
             5.45, 5.50, 5.35,
             5.25, 5.15, 5.00,
             5.85, 5.80, 5.70,
             5.25, 5.20, 5.10,
             5.65, 5.55, 5.45,
             5.60, 5.35, 5.45,
             5.05, 5.00, 4.95,
             5.50, 5.50, 5.40,
             5.45, 5.55, 5.50,
             5.55, 5.55, 5.35,
             5.45, 5.50, 5.55,
             5.50, 5.45, 5.25,
             5.65, 5.60, 5.40,
             5.70, 5.65, 5.55,
             6.30, 6.30, 6.25),
           nrow = 22,
           byrow = TRUE,
           dimnames = list(1 : 22,
                           c("Round Out", "Narrow Angle", "Wide Angle")))
  htest = friedman.test(RoundingTimes)
  df1 = df_htest(htest)
  pr = describe_htest(htest, alpha = .05)
  expect_equal(htest$p.value,df1$p.value)

  # Kruskal-Wallis -----

  ## Hollander & Wolfe (1973), 116.
  ## Mucociliary efficiency from the rate of removal of dust in normal
  ##  subjects, subjects with obstructive airway disease, and subjects
  ##  with asbestosis.
  x <- c(2.9, 3.0, 2.5, 2.6, 3.2) # normal subjects
  y <- c(3.8, 2.7, 4.0, 2.4)      # with obstructive airway disease
  z <- c(2.8, 3.4, 3.7, 2.2, 2.0) # with asbestosis
  htest = kruskal.test(list(x, y, z))
  df1 = df_htest(htest)
  pr = describe_htest(htest, alpha = .05)
  expect_equal(htest$p.value,df1$p.value)

  # Chi-squared -----

  ## Effect of simulating p-values
  x <- matrix(c(12, 5, 7, 7), ncol = 2)
  #chisq.test(x)$p.value           # 0.4233
  htest = chisq.test(x, simulate.p.value = TRUE, B = 10000)
  df1 = df_htest(htest)
  pr = describe_htest(htest, alpha = .05)
  expect_equal(htest$p.value,df1$p.value)

  # Fisher exact ----
  ## Agresti (1990, p. 61f; 2002, p. 91) Fisher's Tea Drinker
  ## A British woman claimed to be able to distinguish whether milk or
  ##  tea was added to the cup first.  To test, she was given 8 cups of
  ##  tea, in four of which milk was added first.  The null hypothesis
  ##  is that there is no association between the true order of pouring
  ##  and the woman's guess, the alternative that there is a positive
  ##  association (that the odds ratio is greater than 1).
  TeaTasting <-
    matrix(c(3, 1, 1, 3),
           nrow = 2,
           dimnames = list(Guess = c("Milk", "Tea"),
                           Truth = c("Milk", "Tea")))
  htest = fisher.test(TeaTasting, alternative = "greater")
  df1 = df_htest(htest)
  pr = describe_htest(htest, alpha = .05)
  expect_equal(htest$p.value,df1$p.value)

  # proportions test -----

  ## Data from Fleiss (1981), p. 139.
  ## H0: The null hypothesis is that the four populations from which
  ##     the patients were drawn have the same true proportion of smokers.
  ## A:  The alternative is that this proportion is different in at
  ##     least one of the populations.

  smokers  <- c( 83, 90, 129, 70 )
  patients <- c( 86, 93, 136, 82 )
  htest = prop.test(smokers, patients)
  df1 = df_htest(htest)
  pr = describe_htest(htest, alpha = .05)
  expect_equal(htest$p.value,df1$p.value)

  # binomial test -----

  ## Conover (1971), p. 97f.
  ## Under (the assumption of) simple Mendelian inheritance, a cross
  ##  between plants of two particular genotypes produces progeny 1/4 of
  ##  which are "dwarf" and 3/4 of which are "giant", respectively.
  ##  In an experiment to determine if this assumption is reasonable, a
  ##  cross results in progeny having 243 dwarf and 682 giant plants.
  ##  If "giant" is taken as success, the null hypothesis is that p =
  ##  3/4 and the alternative that p != 3/4.
  htest = binom.test(c(682, 243), p = 3/4)
  df1 = df_htest(htest)
  pr = describe_htest(htest, alpha = .05)
  expect_equal(htest$p.value,df1$p.value)

})