R/brunnermunzel.permutation.test.R

In brunnermunzel: (Permuted) Brunner-Munzel Test

#' @useDynLib brunnermunzel
NULL

#' permuted Brunner-Munzel test
#'
#' This function performs the permuted Brunner-Munzel test.
#'
#' @return A list containing the following components:
#'  \item{method}{the characters ``permuted Brunner-Munzel Test''}
#'  \item{data.name}{a character string giving the name of the data.}
#'  \item{p.value}{the \eqn{p}-value of the test.}
#'  \item{estimate}{an estimate of the effect size}
#'
#' @references Karin Neubert and Edgar Brunner,
#' ``A studentized permutation test for the non-parametric
#' Behrens-Fisher problem'',
#' Computational Statistics and Data Analysis, Vol. 51, pp. 5192-5204 (2007).
#'
#' @seealso This function is made in reference to following cite (in Japanese):
#' Prof. Haruhiko Okumura
#' (\url{https://oku.edu.mie-u.ac.jp/~okumura/stat/brunner-munzel.html}).
#'
#' @note FORTRAN subroutine `combination` in combination.f is derived from
#' the program by shikino (\url{http://slpr.sakura.ne.jp/qp/combination})
#' (CC-BY-4.0).
#' Thanks to shikono for your useful subroutine.
#'
#' @examples
#' ## Hollander & Wolfe (1973), 29f.
#' ## Hamilton depression scale factor measurements in 9 patients with
#' ##  mixed anxiety and depression, taken at the first (x) and second
#' ##  (y) visit after initiation of a therapy (administration of a
#' ##  tranquilizer).
#' 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)
#'
#' brunnermunzel.permutation.test(x, y)
#' ##
#' ##       permuted Brunner-Munzel Test
#' ##
#' ## data:  x and y
#' ## p-value = 0.158
#' ## sample estimates:
#' ## P(X<Y)+.5*P(X=Y)
#' ##        0.2839506
#'
#' ## 'est' option
#' ## if 'est = "difference"' return P(X<Y) - P(X>Y)
#' brunnermunzel.permutation.test(x, y, est = "difference")
#' ##
#' ##       permuted Brunner-Munzel Test
#' ##
#' ## data:  x and y
#' ## p-value = 0.158
#' ## sample estimates:
#' ## P(X<Y)-P(X>Y)
#' ##    -0.4320988
#'
#'
#' ## Formula interface.
#' dat <- data.frame(
#'     value = c(x, y),
#'     group = factor(rep(c("x", "y"), c(length(x), length(y))),
#'                    levels = c("x", "y"))
#' )
#'
#' brunnermunzel.permutation.test(value ~ group, data = dat)
#' ##
#' ##       permuted Brunner-Munzel Test
#' ##
#' ## data:  value by group
#' ## p-value = 0.158
#' ## sample estimates:
#' ## P(X<Y)+.5*P(X=Y)
#' ##        0.2839506
#'
#'
#' ## Pain score on the third day after surgery for 14 patients under
#' ## the treatment Y and 11 patients under the treatment N
#' ## (see Brunner and Munzel, 2000; Neubert and Brunner, 2007).
#'
#' Y <- c(1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 2, 4, 1, 1)
#' N <- c(3, 3, 4, 3, 1, 2, 3, 1, 1, 5, 4)
#'
#' \donttest{
#' brunnermunzel.permutation.test(Y, N)
#' }
#' ##
#' ##       permuted Brunner-Munzel Test
#' ##
#' ## data:  Y and N
#' ## p-value = 0.008038
#' ## sample estimates:
#' ## P(X<Y)+.5*P(X=Y)
#' ##         0.788961
#'
#' ## Formula interface.
#' dat <- data.frame(
#'     value = c(Y, N),
#'     group = factor(rep(c("Y", "N"), c(length(Y), length(N))),
#'                    levels = c("Y", "N"))
#' )
#'
#' \donttest{
#' brunnermunzel.permutation.test(value ~ group, data = dat)
#' }
#' ##
#' ##       permuted Brunner-Munzel Test
#' ##
#' ## data:  value by group
#' ## p-value = 0.008038
#' ## sample estimates:
#' ## P(X<Y)+.5*P(X=Y)
#' ##         0.788961
#'
#'
#' ## Matrix or Table interface.
#' ##
#' dat1 <- matrix(c(4, 4, 2, 1, 5, 4), nr = 2, byrow = TRUE)
#' dat2 <- as.table(dat1)
#'
#' brunnermunzel.permutation.test(dat1)  # matrix
#' ##
#' ##       permuted Brunner-Munzel Test
#' ##
#' ## data:  Group1 and Group2
#' ## p-value = 0.1593
#' ## sample estimates:
#' ## P(X<Y)+.5*P(X=Y)
#' ##             0.68
#'
#' brunnermunzel.permutation.test(dat2)  # table
#' ##
#' ##       Brunner-Munzel Test
#' ##
#' ##       permuted Brunner-Munzel Test
#' ##
#' ## data:  A and B
#' ## p-value = 0.1593
#' ## sample estimates:
#' ## P(X<Y)+.5*P(X=Y)
#' ##             0.68
#'
#' @export
#'
brunnermunzel.permutation.test <-  function(x, ...)
    UseMethod("brunnermunzel.permutation.test")


#' @rdname brunnermunzel.permutation.test
#' @method brunnermunzel.permutation.test default
#'
#' @importFrom stats setNames terms
#'
#' @param x the numeric vector of data values from the sample 1,
#'  or 2 x n matrix of table
#' (number of row must be 2 and column is ordinal variables).
#' @param y the numeric vector of data values from the sample 2.
#'  If x is matrix or table, y must be missing.
#' @param alternative a character string specifying the alternative
#' hypothesis, must be one of \code{two.sided} (default), \code{greater} or
#' \code{less}. User can specify just the initial letter.
#' @param force
#'   \describe{
#'    \item{FALSE}{(default): If sample size is too large
#'      [number of combinations > 40116600 = choose(28, 14)],
#'      use \code{brunnermunzel.test}.}
#'    \item{TRUE}{: perform permuted Brunner-Munzel test
#'      regardless sample size.}
#'   }
#' @param est a method to calculate estimate and confidence interval,
#' must be either \code{original} (default) or \code{difference}.
#'    \describe{
#'      \item{original}{(default): return \eqn{p = P(X < Y) + 0.5 * P(X = Y)}}
#'      \item{difference}{: return mean difference.
#'                        i.e. \eqn{P(X < Y) - P(X > Y)} = 2 * p - 1}
#'    }
#' This change is proposed by Dr. Julian D. Karch.
#'
#' @export
#'
brunnermunzel.permutation.test.default <-
    function(x, y,
             alternative = c("two.sided", "greater", "less"),
             force = FALSE,
             est = c("original", "difference"),
             ...) {
        alternative <- match.arg(alternative)
        est <- match.arg(est)
        DNAME <-  paste(deparse(substitute(x)), "and",
                        deparse(substitute(y)))

        nx <- length(x); ny <- length(y)
        if (nx == 1L || ny == 1L) stop("not enough observations")
        n_nCr <- choose(nx + ny, nx)

        # switch brunnermunzel.test
        #  when sample number is too large
        if (!force && (n_nCr > 40116600L)) { # choose(28, 14)
            warning(c("Sample number is too large. ",
                      "Using 'brunnermunzel.test'\n"))
            res <- brunnermunzel.test(x, y,
                                      alternative = alternative,
                                      est = est,
                                      ...)
            return(res)
        }

        alter <- switch(alternative,
                        "two.sided" = 1L,
                        "greater" = 2L,
                        "less" = 3L)

        res <- .Fortran("bm_permutation_test",
                        n = as.integer(nx + ny),
                        r = as.integer(nx),
                        n_nCr = as.integer(n_nCr),
                        dat = as.double(c(x, y)),
                        alter = as.integer(alter),
                        pval = numeric(1),
                        pst = numeric(1),
                        PACKAGE = "brunnermunzel")

        if (est == "original") {
            ESTIMATE <- res$pst
            names(ESTIMATE) <- "P(X<Y)+.5*P(X=Y)"
        } else {
            ESTIMATE <- res$pst * 2 - 1
            names(ESTIMATE) <- "P(X<Y)-P(X>Y)"
        }

        structure(
            list(
                method = "permuted Brunner-Munzel Test",
                data.name = DNAME,
                p.value = res$pval,
                estimate = ESTIMATE),
            class = "htest")
    }


#' @rdname brunnermunzel.permutation.test
#' @method brunnermunzel.permutation.test formula
#'
#' @param formula a formula of the form \code{lhs ~ rhs} where \code{lhs}
#' is a numeric variable giving the data values and \code{rhs} a factor
#' with two levels giving the corresponding groups.
#' @param data an optional matrix or data frame (or similar: see
#' \code{\link{model.frame}}) containing the variables in the
#' formula \code{formula}.  By default the variables are taken from
#' \code{environment(formula)}.
#' @param subset an optional vector specifying a subset of observations
#' to be used.
#' @param na.action a function which indicates what should happen when
#' the data contain \code{NA}s.  Defaults to
#' \code{getOption("na.action")}.
#' @param \dots further arguments to be passed to or from methods
#' (This argument is for only formula).
#'
#' @export
#'
brunnermunzel.permutation.test.formula <-
    function(formula, data, subset = NULL, na.action, ...)
{
    if (missing(formula)
       || (length(formula) != 3L)
       || (length(attr(terms(formula[-2L]), "term.labels")) != 1L))
        stop("'formula' missing or incorrect")

    m <- match.call(expand.dots = FALSE)
    if (is.matrix(eval(m$data, parent.frame())))
        m$data <- as.data.frame(data)
    ## need stats:: for non-standard evaluation
    m[[1L]] <- quote(stats::model.frame)
    m$... <- NULL

    mf <- eval(m, parent.frame())
    DNAME <- paste(names(mf), collapse = " by ")
    names(mf) <- NULL
    response <- attr(attr(mf, "terms"), "response")

    g <- factor(mf[[-response]])
    if (nlevels(g) != 2L)
        stop("grouping factor must have exactly 2 levels")

    DATA <- setNames(split(mf[[response]], g), c("x", "y"))
    y <- do.call("brunnermunzel.permutation.test", c(DATA, list(...)))
    y$data.name <- DNAME
    y
}


#' @rdname brunnermunzel.permutation.test
#' @method brunnermunzel.permutation.test matrix
#'
#' @export
#'
brunnermunzel.permutation.test.matrix <- function(x, ...) {
    if (!is.matrix(x)) {
        stop("Class of data must be 'matrix' or 'table'")
    }

    if (nrow(x) != 2L) {
        stop("Number of rows must be 2")
    }

    rname <- rownames(x)
    if (is.null(rname)) {
        rname <- c("Group1", "Group2")
    }
    DNAME <- paste(rname, collapse = " and ")

    lv <- seq_len(ncol(x))
    g1 <- rep(lv, x[1L,])
    g2 <- rep(lv, x[2L,])

    z <- do.call("brunnermunzel.permutation.test",
                 c(list(g1, g2, ...)))
    z$data.name <- DNAME
    z
}


#' @rdname brunnermunzel.permutation.test
#' @method brunnermunzel.permutation.test table
#'
#' @export
#'
brunnermunzel.permutation.test.table <- function(x, ...) {
    brunnermunzel.permutation.test.matrix(x, ...)
}