Running ‘bugfixes.R’
  Running ‘comparisons.R’
  Running ‘regtest_2sample.R’
  Running ‘regtest_Ksample.R’
  Comparing ‘regtest_Ksample.Rout’ to ‘regtest_Ksample.Rout.save’ ... OK
  Running ‘regtest_MTP.R’
  Comparing ‘regtest_MTP.Rout’ to ‘regtest_MTP.Rout.save’ ... OK
  Running ‘regtest_contingency.R’
  Running ‘regtest_correlation.R’
  Comparing ‘regtest_correlation.Rout’ to ‘regtest_correlation.Rout.save’ ... OK
  Running ‘regtest_distribution.R’
  Running ‘regtest_helpers.R’
  Comparing ‘regtest_helpers.Rout’ to ‘regtest_helpers.Rout.save’ ... OK
  Running ‘regtest_midpvalue.R’
  Running ‘regtest_names.R’
  Comparing ‘regtest_names.Rout’ to ‘regtest_names.Rout.save’ ... OK
  Running ‘regtest_size.R’
  Comparing ‘regtest_size.Rout’ to ‘regtest_size.Rout.save’ ... OK
  Running ‘regtest_trafo.R’
  Comparing ‘regtest_trafo.Rout’ to ‘regtest_trafo.Rout.save’ ... OK
Running the tests in ‘tests/bugfixes.R’ failed.
Complete output:
  > ### Regression tests for fixed bugs
  > 
  > suppressWarnings(RNGversion("3.5.2"))
  > set.seed(290875)
  > library("coin")
  Loading required package: survival
  > isequal <- coin:::isequal
  > `%GE%` <- coin:::`%GE%`
  > options(useFancyQuotes = FALSE)
  > 
  > ### I() returns objects of class "AsIs" causing errors in 'trafo'
  > df <- data.frame(x1 = rnorm(100), x2 = rnorm(100), x3 = gl(2, 50),
  +                  x4 = I(letters[1:20]))
  > df <- within(df, {
  +     x5 <- x1 / x2
  +     x6 <- x1 < 0
  + })
  > it1 <- independence_test(I(x1 / x2) ~ x3, data = df)
  > it2 <- independence_test(x5 ~ x3, data = df)
  > stopifnot(identical(statistic(it1), statistic(it2))) # was OK
  > it3 <- independence_test(I(x1 < 0) ~ x3, data = df)
  > it4 <- independence_test(x6 ~ x3, data = df)
  > stopifnot(identical(statistic(it3), statistic(it4))) # wrong sign
  > try(independence_test(x4 ~ x3, data = df))           # coercion to numeric
  Error in ytrafo(object@y) : data class "character" is not supported
  > 
  > ### expectation was wrong when varonly = TRUE in case both
  > ### xtrafo and ytrafo were multivariate
  > if (require("multcomp")) {
  +     df <- data.frame(x = runif(30), y = runif(30), z = gl(3, 10))
  +     a <- independence_test(x + y ~ z, data = df,
  +          distribution = approximate(nresample = 19999),
  +          xtrafo = function(data) trafo(data, factor_trafo = function(x)
  +              model.matrix(~x - 1) %*% t(contrMat(table(x), "Tukey"))))
  +     b <- independence_test(x + y ~ z, data = df,
  +          xtrafo = function(data) trafo(data, factor_trafo = function(x)
  +              model.matrix(~x - 1) %*% t(contrMat(table(x), "Tukey"))))
  +     isequal(expectation(a), expectation(b))
  + }
  Loading required package: multcomp
  Loading required package: mvtnorm
  Loading required package: TH.data
  Loading required package: MASS
  
  Attaching package: 'TH.data'
  
  The following object is masked from 'package:MASS':
  
      geyser
  
  [1] TRUE
  > 
  > 
  > ### 'statistic' for linear and standardized statistics was wrong in case of
  > ### scores
  > data("jobsatisfaction")
  > stopifnot(unique(dim(statistic(lbl_test(jobsatisfaction), "linear"))) == 1)
  > stopifnot(unique(dim(statistic(lbl_test(jobsatisfaction), "standardized"))) == 1)
  > 
  > 
  > ### support() failed in most cases
  > df <- data.frame(x = runif(20), y = runif(20), z = gl(2, 10))
  > stopifnot(is.na(support(independence_test(x ~ z, data = df))))
  > stopifnot(is.na(support(independence_test(x ~ z, data = df, teststat = "quad"))))
  > ite <- independence_test(I(round(x, 1)) ~ z, data = df, dist = exact())
  > ae <- support(ite)
  > de <- dperm(ite, ae)
  > stopifnot(isequal(sum(de), 1))
  > pe <- pperm(ite, ae)
  > stopifnot(isequal(cumsum(de), pe))
  > ita <- independence_test(I(round(x, 1)) ~ z, data = df,
  +                          dist = approximate(nresample = 100000))
  > aa <- support(ita)
  > da <- dperm(ita, aa)
  > stopifnot(isequal(sum(da), 1))
  > pa <- pperm(ita, aa)
  > stopifnot(isequal(cumsum(da), pa))
  > mean(round(ae, 10) %in% round(aa, 10))
  [1] 1
  > 
  > plot(aa, da, type = "s", lty = 1)
  > lines(ae, de, type = "s", lty = 2)
  > itas <- independence_test(I(round(x, 1)) ~ z, data = df)
  > lines(ae[-1], diff(sapply(ae, function(x) pperm(itas, x))), lty = 3)
  > legend("topleft", lty = 1:3, legend = c("approx", "exact", "asympt"), bty = "n")
  > 
  > ### check correct handling of multiple censoring indicators (in modeltools)
  > ### was never wrong, just in case...
  > data("photocar", package = "coin")
  > i1 <- independence_test(Surv(time, event) + Surv(dmin, tumor) + ntumor ~ group,
  +                   data = photocar)
  > i2 <- independence_test(Surv(time, event) ~ group, data = photocar)
  > i3 <- independence_test(Surv(dmin, tumor) ~ group, data = photocar)
  > 
  > stopifnot(max(abs(statistic(i1, "standardized")[,1] -
  +                   statistic(i2, "stand"))) < sqrt(.Machine$double.eps))
  > stopifnot(max(abs(statistic(i1, "standardized")[,2] -
  +                   statistic(i3, "stand"))) < sqrt(.Machine$double.eps))
  > 
  > ### check new var_trafo argument
  > x <- rnorm(20)
  > y <- gl(2, 10)
  > a <- trafo(data.frame(x = x, y = y), numeric_trafo = normal_trafo)
  > b <- trafo(data.frame(x = x, y = y), var_trafo = list(x = normal_trafo))
  > stopifnot(isequal(a, b))
  > 
  > ### check for multiple ordered factors
  > mydf <- data.frame(x = ordered(gl(4, 5)), y = ordered(gl(5, 4)),
  +                    z = rnorm(20))
  > it1 <- independence_test(x + y ~ z , data = mydf)
  > stopifnot(isequal(drop(statistic(it1, "linear")),
  +           c(statistic(independence_test(x ~ z , data = mydf), "linear"),
  +             statistic(independence_test(y ~ z , data = mydf), "linear"))))
  > it1 <- independence_test(x + z ~ y , data = mydf)
  > stopifnot(isequal(drop(statistic(it1, "linear")),
  +           c(statistic(independence_test(x ~ y , data = mydf), "linear"),
  +             statistic(independence_test(z ~ y , data = mydf), "linear"))))
  > it1 <- independence_test(z ~ x + y , data = mydf)
  > stopifnot(isequal(drop(statistic(it1, "linear")),
  +           c(statistic(independence_test(z ~ x , data = mydf), "linear"),
  +             statistic(independence_test(z ~ y , data = mydf), "linear"))))
  > 
  > ### NA's and weights
  > mydf <- data.frame(x = 1:10, y = gl(2, 5), w = rep(2, 10))
  > s <- statistic(independence_test(x ~ y, data = mydf, weights = ~ w), "linear")
  > stopifnot(s == 30)
  > mydf$x[1] <- NA
  > s <- statistic(independence_test(x ~ y, data = mydf, weights = ~ w), "linear")
  > stopifnot(s == 28)
  > 
  > ### two observations only
  > mydf <- data.frame(x = 1:10, y = factor(rep(c(1, 2), 5)))
  > independence_test(y ~ x, data = mydf, subset = c(1, 6))
  
  	Asymptotic General Independence Test
  
  data:  y by x
  Z = -1, p-value = 0.3173
  alternative hypothesis: two.sided
  
  > independence_test(y ~ x, data = mydf, subset = c(1, 2))
  
  	Asymptotic General Independence Test
  
  data:  y by x
  Z = -1, p-value = 0.3173
  alternative hypothesis: two.sided
  
  > try(independence_test(y ~ x, data = mydf, subset = 1))
  Error in validObject(.Object) : 
    invalid class "IndependenceProblem" object: FALSE
  > 
  > ### names of expectation and covariance
  > YOY <- data.frame(length = c(46, 28, 46, 37, 32, 41, 42, 45, 38, 44,
  +                              42, 60, 32, 42, 45, 58, 27, 51, 42, 52,
  +                              38, 33, 26, 25, 28, 28, 26, 27, 27, 27,
  +                              31, 30, 27, 29, 30, 25, 25, 24, 27, 30),
  +                   site = factor(c(rep("I", 10), rep("II", 10),
  +                                   rep("III", 10), rep("IV", 10))))
  > 
  > it <- independence_test(length ~ site, data = YOY,
  +     ytrafo = function(data) trafo(data, numeric_trafo = rank_trafo),
  +     teststat = "quad")
  > expectation(it)
    I  II III  IV 
  205 205 205 205 
  > covariance(it)
              I        II       III        IV
  I   1019.0385 -339.6795 -339.6795 -339.6795
  II  -339.6795 1019.0385 -339.6795 -339.6795
  III -339.6795 -339.6795 1019.0385 -339.6795
  IV  -339.6795 -339.6795 -339.6795 1019.0385
  > 
  > mydf <- data.frame(x = rnorm(10), y = rnorm(10), z = gl(2, 5))
  > it <- independence_test(x + y ~ z, data = mydf)
  > statistic(it, "linear")
            x         y
  1 -2.066068 0.3292961
  > expectation(it)
         1:x        1:y 
  -1.6340021 -0.9750075 
  > covariance(it)
             1:x        1:y
  1:x  0.4417981 -0.8139286
  1:y -0.8139286  2.5231787
  > 
  > ### maxstat_trafo
  > n <- seq(from = 5, to = 100, by = 1)
  > for (i in n) {
  +    x <- round(rnorm(i) * 10, 1)
  +    xm <- maxstat_trafo(x)
  +    stopifnot(min(c(mean(xm[,1]), 1 - mean(xm[,ncol(xm)])) - 0.1) >
  +              -sqrt(.Machine$double.eps))
  + }
  > 
  > ### formula evaluation in 'parent.frame()', spotted by Z
  > foo <- function(x, y) independence_test(y ~ x)
  > a <- 1:10
  > b <- 1:10
  > foo(a, b)
  
  	Asymptotic General Independence Test
  
  data:  y by x
  Z = 3, p-value = 0.0027
  alternative hypothesis: two.sided
  
  > x <- 1
  > y <- 1
  > foo(a, b)
  
  	Asymptotic General Independence Test
  
  data:  y by x
  Z = 3, p-value = 0.0027
  alternative hypothesis: two.sided
  
  > 
  > ### factors with only one level
  > dat <- data.frame(y = rnorm(100), x1 = runif(100), x2 = factor(rep(0, 100)))
  > try(independence_test(y ~ x1  + x2, data = dat))
  Error in factor_trafo(x) : 
    can't deal with factors containing only one level
  > 
  > ### user specified g: names, MC
  > me <- as.table(matrix(c( 6,  8, 10,
  +                32, 47, 20), byrow = TRUE, nrow = 2,
  +     dimnames = list(group = c("In situ", "Control"),
  +                     genotype = c("AA", "AG", "GG"))))
  > medf <- as.data.frame(me)
  > 
  > add <- c(0, 1, 2)
  > dom <- c(0, 1, 1)
  > rez <- c(0, 0, 1)
  > g <- function(x) {
  +     x <- unlist(x)
  +     cbind(add[x], dom[x], rez[x])
  + }
  > it <- independence_test(group ~ genotype,
  +     data = medf, weights = ~ Freq, xtrafo = g)
  > statistic(it, "linear")
     In situ
  X1      28
  X2      18
  X3      10
  > 
  > it <- independence_test(group ~ genotype,
  +     data = medf, weights = ~ Freq, xtrafo = g,
  +     distribution = approximate(nresample = 49999))
  > pvalue(it)
  [1] 0.06848137
  99 percent confidence interval:
   0.06560306 0.07144217 
  
  > 
  > stopifnot(all.equal(statistic(independence_test(t(me), xtrafo = g), "linear"),
  +                     statistic(it, "linear")))
  > 
  > ### alternative trafo for ordered variables didn't work
  > ### spotted by Ludwig Hothorn <hothorn@biostat.uni-hannover.de>
  > tmp <- data.frame(x = ordered(sample(1:3, 20, replace = TRUE)), y = rnorm(20))
  > it1 <- independence_test(y ~ x, data = tmp, scores = list(x = c(1, 1, 2)))
  > g <- function(x) c(1, 1, 2)[unlist(x)]
  > it2 <- independence_test(y ~ x, data = tmp, xtrafo = g)
  > it3 <- independence_test(y ~ x, data = tmp,
  +     xtrafo = function(data) trafo(data, ordered_trafo = g))
  > stopifnot(all.equal(statistic(it1), statistic(it2)))
  > stopifnot(all.equal(statistic(it1), statistic(it3)))
  > 
  > ### precision problems in SR algorithm, >= <= issue
  > ### spotted by "Fay, Michael (NIH/NIAID) [E]" <mfay@niaid.nih.gov>
  > x <- c(1,2,3.1,4,5,6)
  > g <- factor(c(0,0,0,1,1,1))
  > it <- independence_test(x ~ g, distribution = exact())
  > stopifnot(pvalue(it) == 0.1)
  > itMC <- independence_test(x ~ g, distribution = approximate(99999))
  > ci <- attr(pvalue(itMC), "conf.int")
  > stopifnot(ci[1] < pvalue(it) && ci[2] > pvalue(it))
  > 
  > ### any() not applied to logicals
  > ### spotted by R-2.7.0 and Kaspar Daniel Hansen
  > x <- c(1.83,  0.50,  1.62,  2.48, 1.68, 1.88, 1.55, 3.06, 1.30)
  > y <- c(0.878, 0.647, 0.598, 2.05, 1.06, 1.29, 1.06, 3.14, 1.29)
  > ### must not give a warning
  > wilcoxsign_test(x ~ y, alternative = "greater",
  +                 distribution = exact())
  Error in NROW(x) : argument "x" is missing, with no default
  Calls: wilcoxsign_test ... asMethod -> new -> initialize -> initialize -> .local -> NROW
  Execution halted
Running the tests in ‘tests/comparisons.R’ failed.
Complete output:
  > ### compare output for examples from StatXact 6 manual with
  > ### StatXact 6 output as given in the manual
  > 
  > library("coin")
  Loading required package: survival
  > set.seed(290875)
  > isequal <- coin:::isequal
  > options(useFancyQuotes = FALSE)
  > 
  > ### marginal homogenity
  > 
  > ### StatXact 6 manual, page 315
  > presidents <- as.table(matrix(c(28, 7, 13, 27), nrow = 2,
  +     dimnames = list(BeforeTVDebate = c("Carter", "Reagan"),
  +                     AfterTVDebate = c("Carter", "Reagan"))))
  > 
  > bta <- mh_test(presidents)
  Error in NROW(x) : argument "x" is missing, with no default
  Calls: mh_test ... asMethod -> new -> initialize -> initialize -> .local -> NROW
  Execution halted
Running the tests in ‘tests/regtest_2sample.R’ failed.
Complete output:
  > ### Regression tests for the 2 sample problem, i.e.,
  > ### testing the independence of a numeric variable
  > ### 'y' and a binary factor 'x' (possibly blocked)
  > 
  > suppressWarnings(RNGversion("3.5.2"))
  > set.seed(290875)
  > library("coin")
  Loading required package: survival
  > isequal <- coin:::isequal
  > options(useFancyQuotes = FALSE)
  > 
  > ### generate data
  > dat <- data.frame(x = gl(2, 50), y = rnorm(100), block = gl(5, 20))[sample(1:100, 75), ]
  > 
  > 
  > ### Wilcoxon Mann-Whitney Rank Sum Test
  > 
  > ### asymptotic distribution
  > ptwo <- wilcox.test(y ~ x, data = dat, correct = FALSE, exact = FALSE)$p.value
  > pless <- wilcox.test(y ~ x, data = dat, alternative = "less",
  +                      correct = FALSE, exact = FALSE)$p.value
  > pgreater <- wilcox.test(y ~ x, data = dat, alternative = "greater",
  +                         correct = FALSE, exact = FALSE)$p.value
  > 
  > stopifnot(isequal(pvalue(wilcox_test(y ~ x, data = dat)), ptwo))
  > stopifnot(isequal(pvalue(wilcox_test(y ~ x, data = dat, alternative = "less")),
  +                   pless))
  > stopifnot(isequal(pvalue(wilcox_test(y ~ x, data = dat, alternative = "greater")),
  +                   pgreater))
  > 
  > stopifnot(isequal(pvalue(oneway_test(y ~ x, data = dat, distribution = "asympt",
  +     ytrafo = function(data) trafo(data, numeric_trafo = rank_trafo))), ptwo))
  > ### check direct supply of a function via ytrafo
  > stopifnot(isequal(pvalue(oneway_test(y ~ x, data = dat, distribution = "asympt",
  +     ytrafo = rank_trafo)), ptwo))
  > stopifnot(isequal(pvalue(oneway_test(y ~ x, data = dat, distribution = "asympt",
  +     ytrafo = function(data) trafo(data, numeric_trafo = rank_trafo),
  +     alternative = "less")), pless))
  > stopifnot(isequal(pvalue(oneway_test(y ~ x, data = dat, distribution = "asympt",
  +     ytrafo = function(data) trafo(data, numeric_trafo = rank_trafo),
  +     alternative = "greater")), pgreater))
  > 
  > ### exact distribution
  > ptwo <- wilcox.test(y ~ x, data = dat, exact = TRUE)$p.value
  > pless <- wilcox.test(y ~ x, data = dat, alternative = "less", exact = TRUE)$p.value
  > pgreater <- wilcox.test(y ~ x, data = dat, alternative = "greater",
  +                         exact = TRUE)$p.value
  > 
  > stopifnot(isequal(pvalue(wilcox_test(y ~ x, data = dat, distribution = "exact")),
  +                   ptwo))
  > stopifnot(isequal(pvalue(wilcox_test(y ~ x, data = dat, alternative = "less",
  +                              distribution = "exact")), pless))
  > stopifnot(isequal(pvalue(wilcox_test(y ~ x, data = dat, alternative = "greater",
  +                              distribution = "exact")), pgreater))
  > 
  > stopifnot(isequal(pvalue(oneway_test(y ~ x, data = dat, distribution = "exact",
  +     ytrafo = function(data) trafo(data, numeric_trafo = rank_trafo))), ptwo))
  > stopifnot(isequal(pvalue(oneway_test(y ~ x, data = dat, distribution = "exact",
  +     ytrafo = function(data) trafo(data, numeric_trafo = rank_trafo),
  +     alternative = "less")), pless))
  > stopifnot(isequal(pvalue(oneway_test(y ~ x, data = dat, distribution = "exact",
  +     ytrafo = function(data) trafo(data, numeric_trafo = rank_trafo),
  +     alternative = "greater")), pgreater))
  > 
  > ### approximated distribution
  > rtwo <- pvalue(wilcox_test(y ~ x, data = dat, distribution = "approx")) / ptwo
  > rless <- pvalue(wilcox_test(y ~ x, data = dat, alternative = "less",
  +                    distribution = "approx")) / pless
  > rgreater <- pvalue(wilcox_test(y ~ x, data = dat, alternative = "greater",
  +                    distribution = "approx")) / pgreater
  > stopifnot(all(c(rtwo, rless, rgreater) > 0.9 &
  +               c(rtwo, rless, rgreater) < 1.1))
  > 
  > ### <FIXME> add block examples </FIXME>
  > 
  > pvalue(wilcox_test(y ~ x | block, data = dat, distribution = "asympt"))
  [1] 0.1794892
  > pvalue(wilcox_test(y ~ x | block, data = dat, distribution = "approx"))
  [1] 0.1956
  99 percent confidence interval:
   0.1854798 0.2060117 
  
  > pvalue(wilcox_test(y ~ x | block, data = dat, distribution = "exact"))
  [1] 0.1847028
  > 
  > pvalue(wilcox_test(y ~ x | block, data = dat, distribution = "asympt",
  +                    alternative = "less"))
  [1] 0.08974461
  > pvalue(wilcox_test(y ~ x | block, data = dat, distribution = "approx",
  +                    alternative = "less"))
  [1] 0.0885
  99 percent confidence interval:
   0.08133160 0.09606348 
  
  > pvalue(wilcox_test(y ~ x | block, data = dat, distribution = "exact",
  +                    alternative = "less"))
  [1] 0.09335664
  > 
  > pvalue(wilcox_test(y ~ x | block, data = dat, distribution = "asympt",
  +                    alternative = "greater"))
  [1] 0.9102554
  > pvalue(wilcox_test(y ~ x | block, data = dat, distribution = "approx",
  +                    alternative = "greater"))
  [1] 0.9158
  99 percent confidence interval:
   0.9083975 0.9228033 
  
  > pvalue(wilcox_test(y ~ x | block, data = dat, distribution = "exact",
  +                    alternative = "greater"))
  [1] 0.9098776
  > 
  > ### sanity checks
  > try(wilcox_test(x ~ y, data = dat))
  Error in check(object) : 
    'object' does not represent a two-sample problem (maybe the grouping variable is not a factor?)
  > try(wilcox_test(x ~ y | y, data = dat))
  Error in .local(.Object, ...) : 'block' is not a factor
  > 
  > 
  > ### Ansari-Bradley Test
  > 
  > ### asymptotic distribution
  > ptwo <- ansari.test(y ~ x, data = dat, correct = FALSE, exact = FALSE)$p.value
  > pless <- ansari.test(y ~ x, data = dat, alternative = "less",
  +                      correct = FALSE, exact = FALSE)$p.value
  > pgreater <- ansari.test(y ~ x, data = dat, alternative = "greater",
  +                         correct = FALSE, exact = FALSE)$p.value
  > 
  > stopifnot(isequal(pvalue(ansari_test(y ~ x, data = dat)), ptwo))
  > stopifnot(isequal(pvalue(ansari_test(y ~ x, data = dat, alternative = "less")),
  +                   pless))
  > stopifnot(isequal(pvalue(ansari_test(y ~ x, data = dat, alternative = "greater")),
  +                   pgreater))
  > 
  > stopifnot(isequal(pvalue(oneway_test(y ~ x, data = dat, distribution = "asympt",
  +     ytrafo = function(data) trafo(data, numeric_trafo = ansari_trafo))), ptwo))
  > stopifnot(isequal(pvalue(oneway_test(y ~ x, data = dat, distribution = "asympt",
  +     ytrafo = function(data) trafo(data, numeric_trafo = ansari_trafo),
  +     alternative = "greater")), pless))
  > stopifnot(isequal(pvalue(oneway_test(y ~ x, data = dat, distribution = "asympt",
  +     ytrafo = function(data) trafo(data, numeric_trafo = ansari_trafo),
  +     alternative = "less")), pgreater))
  > 
  > ### exact distribution
  > ptwo <- ansari.test(y ~ x, data = dat, exact = TRUE)$p.value
  > pless <- ansari.test(y ~ x, data = dat, alternative = "less", exact = TRUE)$p.value
  > pgreater <- ansari.test(y ~ x, data = dat, alternative = "greater",
  +                         exact = TRUE)$p.value
  > 
  > ### <FIXME>: Definition of two-sided P-values! </FIXME>
  > (isequal(pvalue(ansari_test(y ~ x, data = dat, distribution = "exact")),
  +                   ptwo))
  [1] 0.2125862712
  [1] 0.213574129
  [1] FALSE
  > stopifnot(isequal(pvalue(ansari_test(y ~ x, data = dat, alternative = "less",
  +                              distribution = "exact")), pless))
  > stopifnot(isequal(pvalue(ansari_test(y ~ x, data = dat, alternative = "greater",
  +                              distribution = "exact")), pgreater))
  > 
  > ### <FIXME>: Definition of two-sided P-values! </FIXME>
  > (isequal(pvalue(oneway_test(y ~ x, data = dat, distribution = "exact",
  +     ytrafo = function(data) trafo(data, numeric_trafo = ansari_trafo))), ptwo))
  [1] 0.2125862712
  [1] 0.213574129
  [1] FALSE
  > stopifnot(isequal(pvalue(oneway_test(y ~ x, data = dat, distribution = "exact",
  +     ytrafo = function(data) trafo(data, numeric_trafo = ansari_trafo),
  +     alternative = "greater")), pless))
  > stopifnot(isequal(pvalue(oneway_test(y ~ x, data = dat, distribution = "exact",
  +     ytrafo = function(data) trafo(data, numeric_trafo = ansari_trafo),
  +     alternative = "less")), pgreater))
  > 
  > ### approximated distribution
  > rtwo <- pvalue(ansari_test(y ~ x, data = dat, distribution = "approx")) / ptwo
  > rless <- pvalue(ansari_test(y ~ x, data = dat, alternative = "less",
  +                    distribution = "approx")) / pless
  > rgreater <- pvalue(ansari_test(y ~ x, data = dat, alternative = "greater",
  +                    distribution = "approx")) / pgreater
  > ### <FIXME> ??? </FIXME>
  > (all(c(rtwo, rless, rgreater) > 0.9 &
  +               c(rtwo, rless, rgreater) < 1.1))
  [1] TRUE
  > 
  > ### <FIXME> add block examples </FIXME>
  > 
  > ### sanity checks
  > try(ansari_test(x ~ y, data = dat))
  Error in check(object) : 
    'object' does not represent a K-sample problem (maybe the grouping variable is not a factor?)
  > try(ansari_test(x ~ y | y, data = dat))
  Error in .local(.Object, ...) : 'block' is not a factor
  > 
  > 
  > ### One-way Test
  > oneway_test(y ~ x, dat = dat)
  
  	Asymptotic Two-Sample Fisher-Pitman Permutation Test
  
  data:  y by x (1, 2)
  Z = -0.77735, p-value = 0.4369
  alternative hypothesis: true mu is not equal to 0
  
  > oneway_test(y ~ x, dat = dat, alternative = "less")
  
  	Asymptotic Two-Sample Fisher-Pitman Permutation Test
  
  data:  y by x (1, 2)
  Z = -0.77735, p-value = 0.2185
  alternative hypothesis: true mu is less than 0
  
  > oneway_test(y ~ x, dat = dat, alternative = "greater")
  
  	Asymptotic Two-Sample Fisher-Pitman Permutation Test
  
  data:  y by x (1, 2)
  Z = -0.77735, p-value = 0.7815
  alternative hypothesis: true mu is greater than 0
  
  > 
  > 
  > ### Normal Scores Test
  > normal_test(y ~ x, dat = dat)
  
  	Asymptotic Two-Sample van der Waerden (Normal Quantile) Test
  
  data:  y by x (1, 2)
  Z = -0.85965, p-value = 0.39
  alternative hypothesis: true mu is not equal to 0
  
  > normal_test(y ~ x, dat = dat, alternative = "less")
  
  	Asymptotic Two-Sample van der Waerden (Normal Quantile) Test
  
  data:  y by x (1, 2)
  Z = -0.85965, p-value = 0.195
  alternative hypothesis: true mu is less than 0
  
  > normal_test(y ~ x, dat = dat, alternative = "greater")
  
  	Asymptotic Two-Sample van der Waerden (Normal Quantile) Test
  
  data:  y by x (1, 2)
  Z = -0.85965, p-value = 0.805
  alternative hypothesis: true mu is greater than 0
  
  > 
  > 
  > ### Median Test
  > median_test(y ~ x, dat = dat)
  
  	Asymptotic Two-Sample Brown-Mood Median Test
  
  data:  y by x (1, 2)
  Z = -0.8015, p-value = 0.4228
  alternative hypothesis: true mu is not equal to 0
  
  > median_test(y ~ x, dat = dat, alternative = "less")
  
  	Asymptotic Two-Sample Brown-Mood Median Test
  
  data:  y by x (1, 2)
  Z = -0.8015, p-value = 0.2114
  alternative hypothesis: true mu is less than 0
  
  > median_test(y ~ x, dat = dat, alternative = "greater")
  
  	Asymptotic Two-Sample Brown-Mood Median Test
  
  data:  y by x (1, 2)
  Z = -0.8015, p-value = 0.7886
  alternative hypothesis: true mu is greater than 0
  
  > 
  > 
  > ### Savage Test
  > savage_test(y ~ x, dat = dat)
  
  	Asymptotic Two-Sample Savage Test
  
  data:  y by x (1, 2)
  Z = -1.0738, p-value = 0.2829
  alternative hypothesis: true mu is not equal to 0
  
  > savage_test(y ~ x, dat = dat, alternative = "less")
  
  	Asymptotic Two-Sample Savage Test
  
  data:  y by x (1, 2)
  Z = -1.0738, p-value = 0.1414
  alternative hypothesis: true mu is less than 0
  
  > savage_test(y ~ x, dat = dat, alternative = "greater")
  
  	Asymptotic Two-Sample Savage Test
  
  data:  y by x (1, 2)
  Z = -1.0738, p-value = 0.8586
  alternative hypothesis: true mu is greater than 0
  
  > 
  > 
  > ### Taha Test
  > taha_test(y ~ x, dat = dat)
  
  	Asymptotic Two-Sample Taha Test
  
  data:  y by x (1, 2)
  Z = -1.0835, p-value = 0.2786
  alternative hypothesis: true ratio of scales is not equal to 1
  
  > taha_test(y ~ x, dat = dat, alternative = "less")
  
  	Asymptotic Two-Sample Taha Test
  
  data:  y by x (1, 2)
  Z = -1.0835, p-value = 0.1393
  alternative hypothesis: true ratio of scales is less than 1
  
  > taha_test(y ~ x, dat = dat, alternative = "greater")
  
  	Asymptotic Two-Sample Taha Test
  
  data:  y by x (1, 2)
  Z = -1.0835, p-value = 0.8607
  alternative hypothesis: true ratio of scales is greater than 1
  
  > 
  > 
  > ### Klotz Test
  > klotz_test(y ~ x, dat = dat)
  
  	Asymptotic Two-Sample Klotz Test
  
  data:  y by x (1, 2)
  Z = -0.83541, p-value = 0.4035
  alternative hypothesis: true ratio of scales is not equal to 1
  
  > klotz_test(y ~ x, dat = dat, alternative = "less")
  
  	Asymptotic Two-Sample Klotz Test
  
  data:  y by x (1, 2)
  Z = -0.83541, p-value = 0.2017
  alternative hypothesis: true ratio of scales is less than 1
  
  > klotz_test(y ~ x, dat = dat, alternative = "greater")
  
  	Asymptotic Two-Sample Klotz Test
  
  data:  y by x (1, 2)
  Z = -0.83541, p-value = 0.7983
  alternative hypothesis: true ratio of scales is greater than 1
  
  > 
  > 
  > ### Mood Test
  > mood_test(y ~ x, dat = dat)
  
  	Asymptotic Two-Sample Mood Test
  
  data:  y by x (1, 2)
  Z = -1.1859, p-value = 0.2357
  alternative hypothesis: true ratio of scales is not equal to 1
  
  > mood_test(y ~ x, dat = dat, alternative = "less")
  
  	Asymptotic Two-Sample Mood Test
  
  data:  y by x (1, 2)
  Z = -1.1859, p-value = 0.1178
  alternative hypothesis: true ratio of scales is less than 1
  
  > mood_test(y ~ x, dat = dat, alternative = "greater")
  
  	Asymptotic Two-Sample Mood Test
  
  data:  y by x (1, 2)
  Z = -1.1859, p-value = 0.8822
  alternative hypothesis: true ratio of scales is greater than 1
  
  > 
  > 
  > ### Fligner-Killeen Test
  > fligner_test(y ~ x, dat = dat)
  
  	Asymptotic Two-Sample Fligner-Killeen Test
  
  data:  y by x (1, 2)
  Z = -0.94847, p-value = 0.3429
  alternative hypothesis: true ratio of scales is not equal to 1
  
  > fligner_test(y ~ x, dat = dat, alternative = "less")
  
  	Asymptotic Two-Sample Fligner-Killeen Test
  
  data:  y by x (1, 2)
  Z = -0.94847, p-value = 0.1714
  alternative hypothesis: true ratio of scales is less than 1
  
  > fligner_test(y ~ x, dat = dat, alternative = "greater")
  
  	Asymptotic Two-Sample Fligner-Killeen Test
  
  data:  y by x (1, 2)
  Z = -0.94847, p-value = 0.8286
  alternative hypothesis: true ratio of scales is greater than 1
  
  > 
  > 
  > ### Conover-Iman Test
  > conover_test(y ~ x, dat = dat)
  
  	Asymptotic Two-Sample Conover-Iman Test
  
  data:  y by x (1, 2)
  Z = -0.80308, p-value = 0.4219
  alternative hypothesis: true ratio of scales is not equal to 1
  
  > conover_test(y ~ x, dat = dat, alternative = "less")
  
  	Asymptotic Two-Sample Conover-Iman Test
  
  data:  y by x (1, 2)
  Z = -0.80308, p-value = 0.211
  alternative hypothesis: true ratio of scales is less than 1
  
  > conover_test(y ~ x, dat = dat, alternative = "greater")
  
  	Asymptotic Two-Sample Conover-Iman Test
  
  data:  y by x (1, 2)
  Z = -0.80308, p-value = 0.789
  alternative hypothesis: true ratio of scales is greater than 1
  
  > 
  > 
  > ### Logrank Test
  > logrank_test(Surv(y) ~ x, dat = dat)
  
  	Asymptotic Two-Sample Logrank Test
  
  data:  Surv(y) by x (1, 2)
  Z = -1.0738, p-value = 0.2829
  alternative hypothesis: true theta is not equal to 1
  
  > logrank_test(Surv(y) ~ x, dat = dat, alternative = "less")
  
  	Asymptotic Two-Sample Logrank Test
  
  data:  Surv(y) by x (1, 2)
  Z = -1.0738, p-value = 0.1414
  alternative hypothesis: true theta is less than 1
  
  > logrank_test(Surv(y) ~ x, dat = dat, alternative = "greater")
  
  	Asymptotic Two-Sample Logrank Test
  
  data:  Surv(y) by x (1, 2)
  Z = -1.0738, p-value = 0.8586
  alternative hypothesis: true theta is greater than 1
  
  > 
  > 
  > ### confidence intervals, cf Bauer 1972
  > 
  > ### Location Tests
  > location <- data.frame(y = c(6, 20, 27, 38, 46, 51, 54, 57,
  +                              10, 12, 15, 21, 32, 40, 41, 45),
  +                        x = gl(2, 8))
  > 
  > ### Wilcoxon Rank-Sum Test
  > wt <- wilcox_test(y ~ x, data = location, conf.int = TRUE,
  +                   distribution = "exact")
  > wt
  
  	Exact Wilcoxon-Mann-Whitney Test
  
  data:  y by x (1, 2)
  Z = 1.2603, p-value = 0.2345
  alternative hypothesis: true mu is not equal to 0
  95 percent confidence interval:
   -7 31
  sample estimates:
  difference in location 
                    11.5 
  
  > ci <- confint(wt)
  > wt0 <- wilcox.test(y ~ x, data = location, conf.int = TRUE)
  > stopifnot(isequal(wt0$confint, ci$confint))
  > stopifnot(isequal(wt0$estimate, ci$estimate))
  > 
  > wtx <- wilcox_test(y ~ x, data = location, conf.int = TRUE,
  +                    distribution = "approximate")
  > confint(wtx)
  95 percent confidence interval:
   -9 33 
  sample estimates:
  difference in location 
                    11.5 
  
  > 
  > wta <- wilcox_test(y ~ x, data = location, conf.int = TRUE)
  > confint(wta)
  95 percent confidence interval:
   -7 31 
  sample estimates:
  difference in location 
                11.20138 
  
  > 
  > ### Normal Scores Test
  > nt <- normal_test(y ~ x, data = location, conf.int = TRUE,
  +                   distribution = "exact")
  > nt
  
  	Exact Two-Sample van der Waerden (Normal Quantile) Test
  
  data:  y by x (1, 2)
  Z = 1.2278, p-value = 0.2275
  alternative hypothesis: true mu is not equal to 0
  95 percent confidence interval:
   -6 30
  sample estimates:
  difference in location 
                      11 
  
  > ci <- confint(nt)
  > stopifnot(isequal(ci$conf.int, c(-6, 30)))
  > stopifnot(isequal(ci$estimate, 11))
  > 
  > ntx <- normal_test(y ~ x, data = location, conf.int = TRUE,
  +                    distribution = "approximate")
  > confint(ntx)
  95 percent confidence interval:
   -7 30 
  sample estimates:
  difference in location 
                      11 
  
  > 
  > nta <- normal_test(y ~ x, data = location, conf.int = TRUE)
  > confint(nta)
  95 percent confidence interval:
   -7 30 
  sample estimates:
  difference in location 
                      11 
  
  > 
  > ### Median Test
  > mt <- median_test(y ~ x, data = location, conf.int = TRUE,
  +                   distribution = "exact")
  > mt
  
  	Exact Two-Sample Brown-Mood Median Test
  
  data:  y by x (1, 2)
  Z = 0.96825, p-value = 0.6193
  alternative hypothesis: true mu is not equal to 0
  95 percent confidence interval:
   -21  42
  sample estimates:
  difference in location 
                    15.5 
  
  > confint(mt)
  95 percent confidence interval:
   -21  42 
  sample estimates:
  difference in location 
                    15.5 
  
  > 
  > mtx <- median_test(y ~ x, data = location, conf.int = TRUE,
  +                    distribution = "approximate")
  > confint(mtx)
  95 percent confidence interval:
   -21  42 
  sample estimates:
  difference in location 
                    15.5 
  
  > 
  > mta <- median_test(y ~ x, data = location, conf.int = TRUE)
  > confint(mta)
  95 percent confidence interval:
   -21  42 
  sample estimates:
  difference in location 
                    12.6 
  
  > 
  > ### Savage Test
  > st <- savage_test(y ~ x, data = location, conf.int = TRUE,
  +                   distribution = "exact")
  > st
  
  	Exact Two-Sample Savage Test
  
  data:  y by x (1, 2)
  Z = 1.6095, p-value = 0.1091
  alternative hypothesis: true mu is not equal to 0
  95 percent confidence interval:
   -5 30
  sample estimates:
  difference in location 
                      12 
  
  > confint(st)
  95 percent confidence interval:
   -5 30 
  sample estimates:
  difference in location 
                      12 
  
  > 
  > stx <- savage_test(y ~ x, data = location, conf.int = TRUE,
  +                    distribution = "approximate")
  > confint(stx)
  95 percent confidence interval:
   -5 30 
  sample estimates:
  difference in location 
                      12 
  
  > 
  > sta <- savage_test(y ~ x, data = location, conf.int = TRUE)
  > confint(sta)
  95 percent confidence interval:
   -6 31 
  sample estimates:
  difference in location 
                      12 
  
  > 
  > 
  > ### Scale Tests
  > scale <- data.frame(y = c(-101, -35, -13, 10, 130, 236, 370, 556,
  +                           -145, -140, -40, -30, 2, 27, 68, 290),
  +                     x = gl(2, 8))
  > 
  > ### Ansari-Bradley Test
  > at <- ansari_test(y ~ x, data = scale, conf.int = TRUE, conf.level = 0.988,
  +                   distribution = "exact")
  > at
  
  	Exact Two-Sample Ansari-Bradley Test
  
  data:  y by x (1, 2)
  Z = -0.21129, p-value = 0.9173
  alternative hypothesis: true ratio of scales is not equal to 1
  98.8 percent confidence interval:
    0.1470588 20.5925926
  sample estimates:
  ratio of scales 
         1.221264 
  
  > ci <- confint(at)
  > stopifnot(isequal(ci$conf.int, c(10, 556) / c(68, 27)))
  > stopifnot(isequal(ci$estimate, mean(c(35 / 30, 370 / 290))))
  > 
  > atx <- ansari_test(y ~ x, data = scale, conf.int = TRUE, conf.level = 0.988,
  +                    distribution = "approximate")
  > confint(atx)
  98.8 percent confidence interval:
    0.09285714 20.59259259 
  sample estimates:
  ratio of scales 
         1.221264 
  
  > 
  > ata <- ansari_test(y ~ x, data = scale, conf.int = TRUE, conf.level = 0.988)
  > confint(ata)                                                  # wrong in < 1.3-0
  98.8 percent confidence interval:
    0.1470588 20.5925926 
  sample estimates:
  ratio of scales 
         1.236009 
  
  > 
  > ### Taha Test
  > tt <- taha_test(y ~ x, data = scale, conf.int = TRUE, conf.level = 0.51,
  +                 distribution = "exact")
  > tt
  
  	Exact Two-Sample Taha Test
  
  data:  y by x (1, 2)
  Z = 1.2612, p-value = 0.2157
  alternative hypothesis: true ratio of scales is not equal to 1
  51 percent confidence interval:
     4.814815 278.000000
  sample estimates:
  ratio of scales 
          13.7037 
  
  > confint(tt)
  51 percent confidence interval:
     4.814815 278.000000 
  sample estimates:
  ratio of scales 
          13.7037 
  
  > 
  > ttx <- taha_test(y ~ x, data = scale, conf.int = TRUE, conf.level = 0.51,
  +                  distribution = "approximate")
  > confint(ttx)
  51 percent confidence interval:
     4.814815 278.000000 
  sample estimates:
  ratio of scales 
          13.7037 
  
  > 
  > tta <- taha_test(y ~ x, data = scale, conf.int = TRUE, conf.level = 0.51)
  > confint(tta)                                                  # wrong in < 1.3-0
  51 percent confidence interval:
     4.814815 278.000000 
  sample estimates:
  ratio of scales 
          13.7037 
  
  > 
  > ### Klotz Test
  > kt <- klotz_test(y ~ x, data = scale, conf.int = TRUE, conf.level = 0.988,
  +                  distribution = "exact")
  > kt
  
  	Exact Two-Sample Klotz Test
  
  data:  y by x (1, 2)
  Z = 0.092301, p-value = 0.9243
  alternative hypothesis: true ratio of scales is not equal to 1
  98.8 percent confidence interval:
    0.1470588 13.7037037
  sample estimates:
  ratio of scales 
         1.221264 
  
  > confint(kt)
  98.8 percent confidence interval:
    0.1470588 13.7037037 
  sample estimates:
  ratio of scales 
         1.221264 
  
  > 
  > ktx <- klotz_test(y ~ x, data = scale, conf.int = TRUE, conf.level = 0.988,
  +                   distribution = "approximate")
  > confint(ktx)
  98.8 percent confidence interval:
    0.1470588 13.7037037 
  sample estimates:
  ratio of scales 
         1.221264 
  
  > 
  > kta <- klotz_test(y ~ x, data = scale, conf.int = TRUE, conf.level = 0.988)
  > confint(kta)                                                  # wrong in < 1.3-0
  98.8 percent confidence interval:
    0.08965518 20.59259260 
  sample estimates:
  ratio of scales 
         1.275862 
  
  > 
  > ### Mood Test
  > mt <- mood_test(y ~ x, data = scale, conf.int = TRUE, conf.level = 0.988,
  +                 distribution = "exact")
  > mt
  
  	Exact Two-Sample Mood Test
  
  data:  y by x (1, 2)
  Z = 0.15374, p-value = 0.9021
  alternative hypothesis: true ratio of scales is not equal to 1
  98.8 percent confidence interval:
    0.1470588 13.7037037
  sample estimates:
  ratio of scales 
         1.221264 
  
  > confint(mt)
  98.8 percent confidence interval:
    0.1470588 13.7037037 
  sample estimates:
  ratio of scales 
         1.221264 
  
  > 
  > mtx <- mood_test(y ~ x, data = scale, conf.int = TRUE, conf.level = 0.988,
  +                  distribution = "approximate")
  > confint(mtx)
  98.8 percent confidence interval:
    0.09285714 13.70370370 
  sample estimates:
  ratio of scales 
         1.221264 
  
  > 
  > mta <- mood_test(y ~ x, data = scale, conf.int = TRUE, conf.level = 0.988)
  > confint(mta)                                                  # wrong in < 1.3-0
  98.8 percent confidence interval:
    0.09285714 20.59259259 
  sample estimates:
  ratio of scales 
         1.267119 
  
  > 
  > ### Fligner-Killeen Test
  > ft <- fligner_test(y ~ x, data = scale, conf.int = TRUE, conf.level = 0.988,
  +                    distribution = "exact")
  > ft
  
  	Exact Two-Sample Fligner-Killeen Test
  
  data:  y by x (1, 2)
  Z = 1.4442, p-value = 0.1545
  alternative hypothesis: true ratio of scales is not equal to 1
  98.8 percent confidence interval:
    0.5460526 10.6875000
  sample estimates:
  ratio of scales 
         2.307692 
  
  > confint(ft)
  98.8 percent confidence interval:
    0.5460526 10.6875000 
  sample estimates:
  ratio of scales 
         2.307692 
  
  > 
  > ftx <- fligner_test(y ~ x, data = scale, conf.int = TRUE, conf.level = 0.988,
  +                     distribution = "approximate")
  > confint(ftx)
  98.8 percent confidence interval:
    0.5460526 10.3750000 
  sample estimates:
  ratio of scales 
         2.307692 
  
  > 
  > fta <- fligner_test(y ~ x, data = scale, conf.int = TRUE, conf.level = 0.988)
  > confint(fta)                                                  # wrong in < 1.3-0
  98.8 percent confidence interval:
    0.4761905 10.6875000 
  sample estimates:
  ratio of scales 
         2.290076 
  
  > 
  > ### Conover-Iman Test
  > ct <- conover_test(y ~ x, data = scale, conf.int = TRUE, conf.level = 0.988,
  +                    distribution = "exact")
  > ct
  
  	Exact Two-Sample Conover-Iman Test
  
  data:  y by x (1, 2)
  Z = 1.5314, p-value = 0.1304
  alternative hypothesis: true ratio of scales is not equal to 1
  98.8 percent confidence interval:
   0.7897727 9.8206522
  sample estimates:
  ratio of scales 
         3.048295 
  
  > confint(ct)
  98.8 percent confidence interval:
   0.7897727 9.8206522 
  sample estimates:
  ratio of scales 
         3.048295 
  
  > 
  > ctx <- conover_test(y ~ x, data = scale, conf.int = TRUE, conf.level = 0.988,
  +                     distribution = "approximate")
  > confint(ctx)
  98.8 percent confidence interval:
   0.7897727 9.8206522 
  sample estimates:
  ratio of scales 
         3.048295 
  
  > 
  > cta <- conover_test(y ~ x, data = scale, conf.int = TRUE, conf.level = 0.988)
  > confint(cta)                                                  # wrong in < 1.3-0
  98.8 percent confidence interval:
    0.5493881 10.6576087 
  sample estimates:
  ratio of scales 
         2.455078 
  
  > 
  > 
  > ### ties handling
  > y1 <- c(14, 18, 2, 4, -5, 14, -3, -1, 1, 6, 3, 3)
  > x1 <- c(8, 26, -7, -1, 2, 9, 0, -4, 13, 3, 3, 4)
  > pvalue(wilcoxsign_test(y1 ~ x1, alternative = "greater",
  +                        distribution = "exact", zero.method = "Wilcoxon"))
  Error in NROW(x) : argument "x" is missing, with no default
  Calls: pvalue ... asMethod -> new -> initialize -> initialize -> .local -> NROW
  Execution halted
Running the tests in ‘tests/regtest_contingency.R’ failed.
Complete output:
  > ### Regression tests for the r x c x K problem, i.e.,
  > ### testing the independence of a factor
  > ### 'y' and a factor factor 'x' (possibly blocked)
  > 
  > suppressWarnings(RNGversion("3.5.2"))
  > set.seed(290875)
  > library("coin")
  Loading required package: survival
  > isequal <- coin:::isequal
  > options(useFancyQuotes = FALSE)
  > 
  > 
  > ### generate data: 2 x 2 x K
  > dat <- data.frame(x = gl(2, 50), y = gl(2, 50)[sample(1:100)],
  +                   block = gl(10, 10)[sample(1:100)])[sample(1:100, 75),]
  > 
  > ### Pearsons Chisq Test, asymptotic distribution
  > ptwo <- chisq.test(table(dat$x, dat$y), correct = FALSE)$p.value
  > 
  > stopifnot(isequal(pvalue(chisq_test(y ~ x, data = dat)), ptwo))
  > stopifnot(isequal(pvalue(chisq_test(table(dat$y, dat$x))), ptwo))
  > 
  > ### Cochran-Mantel-Haenzel Test, asymptotic distribution
  > ptwo <- drop(mantelhaen.test(table(dat$x, dat$y, dat$block),
  +                              correct = FALSE)$p.value)
  > 
  > stopifnot(isequal(pvalue(cmh_test(y ~ x | block, data = dat)), ptwo))
  > stopifnot(isequal(pvalue(cmh_test(table(dat$y, dat$x, dat$block))), ptwo))
  > 
  > 
  > ### generate data: r x c x K
  > dat <- data.frame(x = gl(4, 25), y = gl(4, 25)[sample(1:100)],
  +                   block = gl(2, 50)[sample(1:100)])
  > 
  > ### Cochran-Mantel-Haenzel Test, asymptotic distribution
  > ### (was wrong in R < 2.1.0)
  > ptwo <- drop(mantelhaen.test(table(dat$y, dat$x, dat$block),
  +                              correct = FALSE)$p.value)
  > 
  > stopifnot(isequal(pvalue(cmh_test(y ~ x | block, data = dat)), ptwo))
  > stopifnot(isequal(pvalue(cmh_test(table(dat$y, dat$x, dat$block))), ptwo))
  > 
  > ### Linear-by-linear Test, asymptotic distribution
  > lbl_test(y ~ x | block, data = dat)
  
  	Asymptotic Linear-by-Linear Association Test
  
  data:  y (ordered) by x (1 < 2 < 3 < 4) 
  	 stratified by block
  Z = -0.35009, p-value = 0.7263
  alternative hypothesis: two.sided
  
  > lbl_test(table(dat$y, dat$x, dat$block))
  
  	Asymptotic Linear-by-Linear Association Test
  
  data:  Var2 (ordered) by
  	 Var1 (1 < 2 < 3 < 4) 
  	 stratified by Var3
  Z = -0.35009, p-value = 0.7263
  alternative hypothesis: two.sided
  
  > 
  > 
  > ### generate data: r x c x K
  > dat <- data.frame(x = gl(4, 25), y = gl(5, 20)[sample(1:100)],
  +                   block = gl(2, 50)[sample(1:100)])
  > 
  > ### Cochran-Mantel-Haenzel Test, asymptotic distribution
  > ### (was wrong in R < 2.1.0)
  > ptwo <- drop(mantelhaen.test(table(dat$y, dat$x, dat$block),
  +                              correct = FALSE)$p.value)
  > 
  > stopifnot(isequal(pvalue(cmh_test(y ~ x | block, data = dat)), ptwo))
  > stopifnot(isequal(pvalue(cmh_test(table(dat$y, dat$x, dat$block))), ptwo))
  > 
  > ### Linear-by-linear Test, asymptotic distribution
  > lbl_test(y ~ x | block, data = dat)
  
  	Asymptotic Linear-by-Linear Association Test
  
  data:  y (ordered) by x (1 < 2 < 3 < 4) 
  	 stratified by block
  Z = 1.1436, p-value = 0.2528
  alternative hypothesis: two.sided
  
  > lbl_test(table(dat$y, dat$x, dat$block))
  
  	Asymptotic Linear-by-Linear Association Test
  
  data:  Var2 (ordered) by
  	 Var1 (1 < 2 < 3 < 4 < 5) 
  	 stratified by Var3
  Z = 1.1436, p-value = 0.2528
  alternative hypothesis: two.sided
  
  > 
  > 
  > ### 2x2 table and maxstat
  > x <- c(rep(1,51), rep(2,49))
  > y <- factor(c(rep(0,49), rep(1,51)))[sample(1:100)]
  > sit <- as.vector(statistic(independence_test(table(x, y))))
  > stopifnot(isequal(as.vector(statistic(maxstat_test(y ~ x))), sit))
  > stopifnot(isequal(as.vector(statistic(maxstat_test(table(x, y)))), sit))
  > 
  > 
  > ### maxstat for multiple, ordered and unordered covariates
  > dat <- data.frame(w = rnorm(100), x = runif(100), y = gl(4, 25)[sample(1:100)],
  +                   z = ordered(gl(4, 25)[sample(1:100)]))
  > 
  > mt <- maxstat_test(w ~ x, data = dat)
  > mt
  
  	Asymptotic Generalized Maximally Selected Statistics
  
  data:  w by x
  maxT = 2.3738, p-value = 0.184
  alternative hypothesis: two.sided
  sample estimates:
    "best" cutpoint: <= 0.9060508
  
  > est <- mt@estimates$estimate$cutpoint
  > stopifnot(isequal(statistic(mt),
  +                   abs(statistic(independence_test(w ~ (x <= est), data = dat)))))
  > 
  > mt <- maxstat_test(w ~ y, data = dat)
  > mt
  
  	Asymptotic Generalized Maximally Selected Statistics
  
  data:  w by y (1, 2, 3, 4)
  maxT = 0.7294, p-value = 0.8905
  alternative hypothesis: two.sided
  sample estimates:
    "best" cutpoint: {1, 2, 3} vs. {4}
  
  > est <- mt@estimates$estimate$cutpoint
  > xx <- dat$y %in% est
  > stopifnot(isequal(statistic(mt),
  +                   abs(statistic(independence_test(w ~ xx, data = dat)))))
  > 
  > mt <- maxstat_test(w ~ z, data = dat)
  > mt
  
  	Asymptotic Generalized Maximally Selected Statistics
  
  data:  w by z (1 < 2 < 3 < 4)
  maxT = 0.52218, p-value = 0.9125
  alternative hypothesis: two.sided
  sample estimates:
    "best" cutpoint: {1} vs. {2, 3, 4}
  
  > est <- mt@estimates$estimate$cutpoint
  > xx <- dat$z <= est
  > stopifnot(isequal(statistic(mt),
  +                   abs(statistic(independence_test(w ~ xx, data = dat)))))
  > 
  > mt <- maxstat_test(w ~ x + y + z, data = dat)
  > mt
  
  	Asymptotic Generalized Maximally Selected Statistics
  
  data:  w by x, y, z(ordered)
  maxT = 2.3738, p-value = 0.286
  alternative hypothesis: two.sided
  sample estimates:
    "best" cutpoint: <= 0.9060508
         covariable: x
  
  > est <- mt@estimates$estimate
  > xsel <- dat[[est[[1]]]]
  > if (is.factor(xsel) && !is.ordered(xsel)) {
  +     xx <- xsel %in% est[[2]]
  + } else {
  +     xx <- xsel <= est[[2]]
  + }
  > stopifnot(isequal(statistic(mt),
  +                   abs(statistic(independence_test(w ~ xx, data = dat)))))
  > 
  > 
  > ### marginal homogeneity
  > rating <- c("low", "moderate", "high")
  > x <- as.table(matrix(c(20, 10,  5,
  +                        3, 30, 15,
  +                        0,  5, 40),
  +                      ncol = 3, byrow = TRUE,
  +                      dimnames = list(Rater1 = rating, Rater2 = rating)))
  > ### test statistic W_0 = 13.76
  > ### see http://ourworld.compuserve.com/homepages/jsuebersax/mcnemar.htm
  > stopifnot(all.equal(round(statistic(mh_test(x)), 2), 13.76))
  Error in NROW(x) : argument "x" is missing, with no default
  Calls: stopifnot ... asMethod -> new -> initialize -> initialize -> .local -> NROW
  Execution halted
Running the tests in ‘tests/regtest_distribution.R’ failed.
Complete output:
  > ### Regression tests for the distribution functions
  > 
  > suppressWarnings(RNGversion("3.5.2"))
  > set.seed(290875)
  > library("coin")
  Loading required package: survival
  > isequal <- coin:::isequal
  > options(useFancyQuotes = FALSE)
  > 
  > 
  > ### generate independent two-sample data
  > dta <- data.frame(y = rnorm(20), x = gl(2, 10), b = factor(rep(1:4, 5)),
  +                   w = rep(1:3, length = 20))
  > dta$y5 <- round(dta$y, 5)
  > dta$y3 <- round(dta$y, 3)
  > 
  > ### generate paired two-sample data
  > dta2 <- data.frame(y = as.vector(rbind(abs(dta$y) * (dta$y >= 0),
  +                                        abs(dta$y) * (dta$y < 0))),
  +                    x = factor(rep(0:1, length(dta$y)),
  +                               labels = c("pos", "neg")),
  +                    b = gl(length(dta$y), 2))
  > dta2$y5 <- round(dta2$y, 5)
  > dta2$y3 <- round(dta2$y, 3)
  > 
  > 
  > ### check 'algorithm = "auto"'
  > 
  > ### scores that can't be mapped into integers
  > 
  > ### two-sample with block
  > try(independence_test(y ~ x | b, data = dta,
  +                       distribution = exact(algo = "auto")))
  Error in vdW_split_up_2sample(object) : 
    cannot compute exact p-values with blocks
  > 
  > try(independence_test(y ~ x | b, data = dta,
  +                       distribution = exact(algo = "shift")))
  Error in .local(object, ...) : 
    cannot compute exact distribution with real-valued scores
  > 
  > try(independence_test(y ~ x | b, data = dta,
  +                       distribution = exact(algo = "split-up")))
  Error in vdW_split_up_2sample(object) : 
    cannot compute exact p-values with blocks
  > 
  > ### two-sample without block
  > it <- independence_test(y ~ x, data = dta,
  +                         distribution = exact(algo = "auto"))
  > it@distribution@name
  [1] "Exact Distribution for Independent Two-Sample Tests (van de Wiel Split-Up Algorithm)"
  > 
  > try(independence_test(y ~ x, data = dta,
  +                       distribution = exact(algo = "shift")))
  Error in .local(object, ...) : 
    cannot compute exact distribution with real-valued scores
  > 
  > it <- independence_test(y ~ x, data = dta,
  +                         distribution = exact(algo = "split-up"))
  > it@distribution@name
  [1] "Exact Distribution for Independent Two-Sample Tests (van de Wiel Split-Up Algorithm)"
  > pvalue(it)
  [1] 0.9417069
  > midpvalue(it)
  [1] NA
  > pvalue_interval(it)
  [1] NA
  > 
  > ### paired two-sample
  > try(symmetry_test(y ~ x | b, data = dta2, paired = TRUE,
  +                   distribution = exact(algo = "auto")))
  Error in NROW(x) : argument "x" is missing, with no default
  > 
  > try(symmetry_test(y ~ x | b, data = dta2, paired = TRUE,
  +                   distribution = exact(algo = "shift")))
  Error in NROW(x) : argument "x" is missing, with no default
  > 
  > try(symmetry_test(y ~ x | b, data = dta2, paired = TRUE,
  +                   distribution = exact(algo = "split-up")))
  Error in NROW(x) : argument "x" is missing, with no default
  > 
  > ### mapped into integers using 'fact'
  > 
  > ### two-sample with block
  > it <- independence_test(y5 ~ x | b, data = dta,
  +                         distribution = exact(algo = "auto", fact = 1e5))
  > it@distribution@name
  [1] "Exact Distribution for Independent Two-Sample Tests (Streitberg-Roehmel Shift Algorithm)"
  > 
  > it <- independence_test(y5 ~ x | b, data = dta,
  +                         distribution = exact(algo = "shift", fact = 1e5))
  > it@distribution@name
  [1] "Exact Distribution for Independent Two-Sample Tests (Streitberg-Roehmel Shift Algorithm)"
  > pvalue(it)
  [1] 0.9247
  > midpvalue(it)
  [1] 0.92465
  > pvalue_interval(it)
     p_0    p_1 
  0.9246 0.9247 
  > s <- sample(support(it), 5)
  > pvalue(it@distribution, s)
  [1] 0.9608 0.1270 0.0178 0.3248 0.4035
  > midpvalue(it@distribution, s)
  [1] 0.96075 0.12695 0.01775 0.32475 0.40345
  > pvalue_interval(it@distribution, s)
        [,1]   [,2]   [,3]   [,4]   [,5]
  p_0 0.9607 0.1269 0.0177 0.3247 0.4034
  p_1 0.9608 0.1270 0.0178 0.3248 0.4035
  > 
  > try(independence_test(y5 ~ x | b, data = dta,
  +                         distribution = exact(algo = "split-up")))
  Error in vdW_split_up_2sample(object) : 
    cannot compute exact p-values with blocks
  > 
  > ### two-sample without block
  > it <- independence_test(y5 ~ x, data = dta,
  +                         distribution = exact(algo = "auto", fact = 1e5))
  > it@distribution@name
  [1] "Exact Distribution for Independent Two-Sample Tests (Streitberg-Roehmel Shift Algorithm)"
  > 
  > it <- independence_test(y5 ~ x, data = dta,
  +                         distribution = exact(algo = "shift", fact = 1e5))
  > it@distribution@name
  [1] "Exact Distribution for Independent Two-Sample Tests (Streitberg-Roehmel Shift Algorithm)"
  > pvalue(it)
  [1] 0.9417177
  > midpvalue(it)
  [1] 0.9417069
  > pvalue_interval(it)
        p_0       p_1 
  0.9416961 0.9417177 
  > s <- sample(support(it), 5)
  > pvalue(it@distribution, s)
  [1] 0.29646669 0.07164043 0.29003659 0.11804759 0.92572907
  > midpvalue(it@distribution, s)
  [1] 0.29646128 0.07163502 0.29003118 0.11804217 0.92572366
  > pvalue_interval(it@distribution, s)
           [,1]       [,2]      [,3]      [,4]      [,5]
  p_0 0.2964559 0.07162961 0.2900258 0.1180368 0.9257182
  p_1 0.2964667 0.07164043 0.2900366 0.1180476 0.9257291
  > 
  > it <- independence_test(y5 ~ x, data = dta,
  +                         distribution = exact(algo = "split-up"))
  > it@distribution@name
  [1] "Exact Distribution for Independent Two-Sample Tests (van de Wiel Split-Up Algorithm)"
  > pvalue(it)
  [1] 0.9417177
  > midpvalue(it)
  [1] NA
  > pvalue_interval(it)
  [1] NA
  > pvalue(it@distribution, s)
  [1] 0.29646669 0.07164043 0.29003659 0.11804759 0.92572907
  > midpvalue(it@distribution, s)
  [1] NA
  > pvalue_interval(it@distribution, s)
  [1] NA
  > 
  > ### paired two-sample
  > st <- symmetry_test(y5 ~ x | b, data = dta2, paired = TRUE,
  +                     distribution = exact(algo = "auto", fact = 1e5))
  Error in NROW(x) : argument "x" is missing, with no default
  Calls: symmetry_test ... asMethod -> new -> initialize -> initialize -> .local -> NROW
  Execution halted
Running the tests in ‘tests/regtest_midpvalue.R’ failed.
Complete output:
  > ### Regression tests for mid-p confidence intervals and mid-p-values
  > 
  > set.seed(290875)
  > library("coin")
  Loading required package: survival
  > isequal <- coin:::isequal
  > options(useFancyQuotes = FALSE)
  > 
  > 
  > ###
  > ### Mid-p confidence intervals
  > ###
  > 
  > 
  > ### Berry and Armitage (1995, p. 420)
  > mpci <- coin:::confint_binom(5, 20, level = 0.95, method = "mid-p")
  > stopifnot(isequal(round(mpci, 3), c(0.098, 0.470)))
  > 
  > 
  > ### Agresti and Gottard (2001, p. 369)
  > mpci <- coin:::confint_binom(4, 10, level = 0.95, method = "mid-p")
  > stopifnot(isequal(round(mpci, 3), c(0.142, 0.709)))
  > 
  > mpci <- coin:::confint_binom(0, 10, level = 0.95, method = "mid-p")
  > stopifnot(isequal(round(mpci[1], 3), 0))
  > 
  > mpci <- coin:::confint_binom(10, 10, level = 0.95, method = "mid-p")
  > stopifnot(isequal(round(mpci[2], 3), 1))
  > 
  > 
  > ### Newcombe (1998, p. 861, Tab. I)
  > mpci <- coin:::confint_binom(81, 263, level = 0.95, method = "mid-p")
  > stopifnot(isequal(round(mpci, 4), c(0.2544, 0.3658)))
  > 
  > mpci <- coin:::confint_binom(15, 148, level = 0.95, method = "mid-p")
  > stopifnot(isequal(round(mpci, 4), c(0.0601, 0.1581)))
  > 
  > mpci <- coin:::confint_binom(0, 20, level = 0.95, method = "mid-p")
  > stopifnot(isequal(round(mpci, 4), c(0.0000, 0.1391)))
  > 
  > mpci <- coin:::confint_binom(1, 29, level = 0.95, method = "mid-p")
  > stopifnot(isequal(round(mpci, 4), c(0.0017, 0.1585)))
  > 
  > 
  > ###
  > ### Mid-p-pvalues
  > ###
  > 
  > 
  > ### Data from Hwang and Yang (2001, p. 810)
  > tea <- matrix(c(3, 1,
  +                 1, 3),
  +               nrow = 2, byrow = TRUE)
  > 
  > ### Results from Hwang and Yang (2001, p. 810)
  > ct_e <- chisq_test(as.table(tea),
  +                    distribution = "exact")
  > stopifnot(isequal(round(pvalue(ct_e), 3), 0.486))
  > stopifnot(isequal(round(midpvalue(ct_e), 3), 0.257))
  > 
  > it_e_s <- independence_test(as.table(tea),
  +                             distribution = "exact",
  +                             teststat = "scalar")
  > stopifnot(isequal(pvalue(it_e_s), pvalue(ct_e)))
  > stopifnot(isequal(midpvalue(it_e_s), midpvalue(ct_e)))
  > 
  > it_e_q <- independence_test(as.table(tea),
  +                             distribution = "exact",
  +                             teststat = "quad")
  > stopifnot(isequal(pvalue(it_e_q), pvalue(ct_e)))
  > stopifnot(isequal(midpvalue(it_e_q), midpvalue(ct_e)))
  > 
  > it_e_s_gr <- independence_test(as.table(tea),
  +                                distribution = "exact",
  +                                teststat = "scalar",
  +                                alternative = "greater")
  > stopifnot(isequal(round(pvalue(it_e_s_gr), 3), 0.243))
  > stopifnot(isequal(round(midpvalue(it_e_s_gr), 4), 0.1286))
  >                                  # p = 0.1285 according to Hwang and Yang (2001)
  > 
  > ### Additional results: Monte Carlo
  > set.seed(290875)
  > ct_m <- chisq_test(as.table(tea),
  +                    distribution = "approximate")
  > (p <- pvalue(ct_m))
  [1] 0.4835
  99 percent confidence interval:
   0.4705882 0.4964275 
  
  > pci <- attr(p, "conf.int")
  > stopifnot(pci[1] < 0.486 & pci[2] > 0.486)
  > (mp <- midpvalue(ct_m))
  [1] 0.25755
  99 percent confidence interval:
   0.2463534 0.2688772 
  
  > mpci <- attr(mp, "conf.int")
  > stopifnot(mpci[1] < 0.257 & mpci[2] > 0.257)
  > 
  > set.seed(290875)
  > it_m_s <- independence_test(as.table(tea),
  +                             distribution = approximate(nresample = 10000),
  +                             teststat = "scalar")
  > stopifnot(isequal(pvalue(it_m_s), pvalue(ct_m)))
  > stopifnot(isequal(midpvalue(it_m_s), midpvalue(ct_m)))
  > 
  > set.seed(290875)
  > it_m_q <- independence_test(as.table(tea),
  +                             distribution = approximate(nresample = 10000),
  +                             teststat = "quad")
  > stopifnot(isequal(pvalue(it_m_q), pvalue(ct_m)))
  > stopifnot(isequal(midpvalue(it_m_q), midpvalue(ct_m)))
  > 
  > set.seed(290875)
  > it_m_s_gr <- independence_test(as.table(tea),
  +                                distribution = approximate(nresample = 10000),
  +                                teststat = "scalar",
  +                                alternative = "greater")
  > (p <- pvalue(it_m_s_gr))
  [1] 0.2419
  99 percent confidence interval:
   0.2309447 0.2531023 
  
  > pci <- attr(p, "conf.int")
  > stopifnot(pci[1] < 0.243 & pci[2] > 0.243)
  > (mp <- midpvalue(it_m_s_gr))
  [1] 0.1288
  99 percent confidence interval:
   0.1203477 0.1376053 
  
  > mpci <- attr(mp, "conf.int")
  > stopifnot(mpci[1] < 0.1286 & mpci[2] > 0.1286)
  > 
  > 
  > ### Data from Lydersen and Laake (2003, p. 3862)
  > davis <- matrix(c(3,  6,
  +                   2, 19),
  +                 nrow = 2, byrow = TRUE)
  > 
  > ### Results from Lydersen and Laake (2003, p. 3863, Tab. II)
  > ct_e <- chisq_test(as.table(davis),
  +                    distribution = "exact")
  > stopifnot(isequal(round(pvalue(ct_e), 4), 0.2860))
  > stopifnot(isequal(round(midpvalue(ct_e), 4), 0.1527))
  > 
  > it_e_s <- independence_test(as.table(davis),
  +                             distribution = "exact",
  +                             teststat = "scalar")
  > stopifnot(isequal(pvalue(it_e_s), pvalue(ct_e)))
  > stopifnot(isequal(midpvalue(it_e_s), midpvalue(ct_e)))
  > 
  > it_e_q <- independence_test(as.table(davis),
  +                             distribution = "exact",
  +                             teststat = "quad")
  > stopifnot(isequal(pvalue(it_e_q), pvalue(ct_e)))
  > stopifnot(isequal(midpvalue(it_e_q), midpvalue(ct_e)))
  > 
  > ### Additional results: Monte Carlo
  > set.seed(290875)
  > ct_m <- chisq_test(as.table(davis),
  +                    distribution = "approximate")
  > (p <- pvalue(ct_m))
  [1] 0.2847
  99 percent confidence interval:
   0.273131 0.296475 
  
  > pci <- attr(p, "conf.int")
  > stopifnot(pci[1] < 0.2860 & pci[2] > 0.2860)
  > (mp <- midpvalue(ct_m))
  [1] 0.15175
  99 percent confidence interval:
   0.1427229 0.1612083 
  
  > mpci <- attr(mp, "conf.int")
  > stopifnot(mpci[1] < 0.1527 & mpci[2] > 0.1527)
  > 
  > set.seed(290875)
  > it_m_s <- independence_test(as.table(davis),
  +                             distribution = approximate(nresample = 10000),
  +                             teststat = "scalar")
  > stopifnot(isequal(pvalue(it_m_s), pvalue(ct_m)))
  > stopifnot(isequal(midpvalue(it_m_s), midpvalue(ct_m)))
  > 
  > set.seed(290875)
  > it_m_q <- independence_test(as.table(davis),
  +                             distribution = approximate(nresample = 10000),
  +                             teststat = "quad")
  > stopifnot(isequal(pvalue(it_m_q), pvalue(ct_m)))
  > stopifnot(isequal(midpvalue(it_m_q), midpvalue(ct_m)))
  > 
  > 
  > ### Data from Lydersen, Fagerland and Laake (2009, p. 1160, Tab. I)
  > cardiac <- matrix(c(1, 33,
  +                     7, 27),
  +                   nrow = 2, byrow = TRUE)
  > 
  > ### Results from Lydersen, Fagerland and Laake (2009, pp. 1171--1172)
  > ct_e <- chisq_test(as.table(cardiac),
  +                    distribution = "exact")
  > stopifnot(isequal(round(pvalue(ct_e), 4), 0.0544))
  > stopifnot(isequal(round(midpvalue(ct_e), 4), 0.0297))
  > 
  > it_e_s <- independence_test(as.table(cardiac),
  +                             distribution = "exact",
  +                             teststat = "scalar")
  > stopifnot(isequal(pvalue(it_e_s), pvalue(ct_e)))
  > stopifnot(isequal(midpvalue(it_e_s), midpvalue(ct_e)))
  > 
  > it_e_q <- independence_test(as.table(cardiac),
  +                             distribution = "exact",
  +                             teststat = "quad")
  > stopifnot(isequal(pvalue(it_e_q), pvalue(ct_e)))
  > stopifnot(isequal(midpvalue(it_e_q), midpvalue(ct_e)))
  > 
  > ### Additional results: Monte Carlo
  > set.seed(290875)
  > ct_m <- chisq_test(as.table(cardiac),
  +                    distribution = approximate(nresample = 10000))
  > (p <- pvalue(ct_m))
  [1] 0.0557
  99 percent confidence interval:
   0.04995533 0.06187217 
  
  > pci <- attr(p, "conf.int")
  > stopifnot(pci[1] < 0.0544 & pci[2] > 0.0544)
  > (mp <- midpvalue(ct_m))
  [1] 0.0302
  99 percent confidence interval:
   0.02601100 0.03483585 
  
  > mpci <- attr(mp, "conf.int")
  > stopifnot(mpci[1] < 0.0297 & mpci[2] > 0.0297)
  > 
  > set.seed(290875)
  > it_m_s <- independence_test(as.table(cardiac),
  +                             distribution = approximate(nresample = 10000),
  +                             teststat = "scalar")
  > stopifnot(isequal(pvalue(it_m_s), pvalue(ct_m)))
  > stopifnot(isequal(midpvalue(it_m_s), midpvalue(ct_m)))
  > 
  > set.seed(290875)
  > it_m_q <- independence_test(as.table(cardiac),
  +                             distribution = approximate(nresample = 10000),
  +                             teststat = "quad")
  > stopifnot(isequal(pvalue(it_m_q), pvalue(ct_m)))
  > stopifnot(isequal(midpvalue(it_m_q), midpvalue(ct_m)))
  > 
  > 
  > ### Data from Lydersen, Fagerland and Laake (2009, p. 1160, Tab. II)
  > exfoliative <- matrix(c( 0, 16,
  +                         15, 57),
  +                       nrow = 2, byrow = TRUE)
  > 
  > ### Results from Lydersen, Fagerland and Laake (2009, p. 1173)
  > ct_e <- chisq_test(as.table(exfoliative),
  +                    distribution = "exact")
  > stopifnot(isequal(round(pvalue(ct_e), 4), 0.0629))
  > stopifnot(isequal(round(midpvalue(ct_e), 4), 0.0447))
  > 
  > it_e_s <- independence_test(as.table(exfoliative),
  +                             distribution = "exact",
  +                             teststat = "scalar")
  > stopifnot(isequal(pvalue(it_e_s), pvalue(ct_e)))
  > stopifnot(isequal(midpvalue(it_e_s), midpvalue(ct_e)))
  > 
  > it_e_q <- independence_test(as.table(exfoliative),
  +                             distribution = "exact",
  +                             teststat = "quad")
  > stopifnot(isequal(pvalue(it_e_q), pvalue(ct_e)))
  > stopifnot(isequal(midpvalue(it_e_q), midpvalue(ct_e)))
  > 
  > ### Additional results: Monte Carlo
  > set.seed(290875)
  > ct_m <- chisq_test(as.table(exfoliative),
  +                    distribution = approximate(nresample = 10000))
  > (p <- pvalue(ct_m))
  [1] 0.0643
  99 percent confidence interval:
   0.05814061 0.07087830 
  
  > pci <- attr(p, "conf.int")
  > stopifnot(pci[1] < 0.0629 & pci[2] > 0.0629)
  > (mp <- midpvalue(ct_m))
  [1] 0.0444
  99 percent confidence interval:
   0.03930789 0.04992546 
  
  > mpci <- attr(mp, "conf.int")
  > stopifnot(mpci[1] < 0.0447 & mpci[2] > 0.0447)
  > 
  > set.seed(290875)
  > it_m_s <- independence_test(as.table(exfoliative),
  +                             distribution = approximate(nresample = 10000),
  +                             teststat = "scalar")
  > stopifnot(isequal(pvalue(it_m_s), pvalue(ct_m)))
  > stopifnot(isequal(midpvalue(it_m_s), midpvalue(ct_m)))
  > 
  > set.seed(290875)
  > it_m_q <- independence_test(as.table(exfoliative),
  +                             distribution = approximate(nresample = 10000),
  +                             teststat = "quad")
  > stopifnot(isequal(pvalue(it_m_q), pvalue(ct_m)))
  > stopifnot(isequal(midpvalue(it_m_q), midpvalue(ct_m)))
  > 
  > 
  > ### Data from Fagerland, Lydersen and Laake (2013, p. 2, Tab. 1)
  > ahr <- matrix(c(1,  1,
  +                 7, 12),
  +               nrow = 2, byrow = TRUE)
  > 
  > ### Results from Fagerland, Lydersen and Laake (2013, p. 7, Tab. 6)
  > mt_e <- mh_test(as.table(ahr),
  +                 distribution = "exact")
  Error in NROW(x) : argument "x" is missing, with no default
  Calls: mh_test ... asMethod -> new -> initialize -> initialize -> .local -> NROW
  Execution halted