R/rateratio_test.R

In lucavd/ratetest: Test if the rate of a event in two groups is different

# this function is from the package 'rateratio.test' v. 1.0-2 by Michael Fay
# https://cran.rstudio.com/web/packages/rateratio.test/index.html
# included as is to avoid dependencies and enhance reproducibility over time

#' An Exact Rate Ratio Test Assuming Poisson Counts
#' Performs the uniformy most powerful unbiased test on the ratio of rates of two Poisson counts with given time (e.g., perons-years) at risk for each count.
#' @param x
#' @param n
#' @param RR
#' @param alternative
#' @param conf.level
#'
#' @return An object of class ‘htest’
#' @export
#'
#' @examples
#' rateratio.test(c(2,9),c(17877,16660))
#'
rateratio.test <- function (x, n, RR = 1, alternative = c("two.sided", "less",
                                        "greater"), conf.level = 0.95)
{
  DNAME <- deparse(substitute(x))
  if (is.matrix(x)) {
    stop("'x' must be a vector")
  }
  else {
    DNAME <- paste(DNAME, "with time of", deparse(substitute(n)))
    if ((l <- length(x)) != length(n))
      stop("'x' and 'n' must have the same length")
  }
  OK <- complete.cases(x, n)
  x <- x[OK]
  n <- n[OK]
  if ((k <- length(x)) != 2)
    stop("x must have a length 2")
  if (any(n <= 0))
    stop("elements of 'n' must be positive")
  if (any(x < 0))
    stop("elements of 'x' must be nonnegative")
  if (!is.null(RR)) {
    DNAME <- paste(DNAME, ", null rate ratio ", deparse(substitute(RR)),
                   sep = "")
    if (RR <= 0)
      stop("RR must be greater than 0")
  }
  alternative <- match.arg(alternative)
  if ((length(conf.level) != 1) || is.na(conf.level) || (conf.level <=
                                                         0) || (conf.level >= 1))
    stop("'conf.level' must be a single number between 0 and 1")
  Y <- x[1]
  N <- n[1]
  X <- x[2]
  M <- n[2]
  ESTIMATE <- c((Y/N)/(X/M), Y/N, X/M)
  names(ESTIMATE) <- c("Rate Ratio", "Rate 1",
                       "Rate 2")
  pRR <- (N * RR)/(N * RR + M)
  pval.less <- pbinom(Y, X + Y, pRR)
  pval.greater <- pbinom(Y - 1, Y + X, pRR, lower.tail = FALSE)
  PVAL <- switch(alternative, less = pval.less, greater = pval.greater,
                 two.sided = min(1, 2 * min(pval.less, pval.greater)))
  p.L <- function(x, n, alpha) {
    if (x == 0)
      0
    else qbeta(alpha, x, n - x + 1)
  }
  p.U <- function(x, n, alpha) {
    if (x == n)
      1
    else qbeta(1 - alpha, x + 1, n - x)
  }
  CINT <- switch(alternative, less = c(0, (p.U(Y, X + Y, 1 -
                                                 conf.level) * M)/(N * (1 - p.U(Y, X + Y, 1 - conf.level)))),
                 greater = c((p.L(Y, X + Y, 1 - conf.level) * M)/(N *
                                                                    (1 - p.L(Y, X + Y, 1 - conf.level))), Inf), two.sided = c((p.L(Y,
                                                                                                                                   X + Y, (1 - conf.level)/2) * M)/(N * (1 - p.L(Y,
                                                                                                                                                                                 X + Y, (1 - conf.level)/2))), (p.U(Y, X + Y, (1 -
                                                                                                                                                                                                                                 conf.level)/2) * M)/(N * (1 - p.U(Y, X + Y, (1 -
                                                                                                                                                                                                                                                                                conf.level)/2)))))
  attr(CINT, "conf.level") <- conf.level
  names(RR) <- "rate ratio"
  RVAL <- list(p.value = as.numeric(PVAL), estimate = ESTIMATE,
               null.value = RR, conf.int = CINT, alternative = alternative,
               method = "Exact Rate Ratio Test, assuming Poisson counts",
               data.name = DNAME)
  class(RVAL) <- "htest"
  return(RVAL)
}