R/ztest.paired.R

In powertools: Power and Sample Size Tools

#' Power calculation for paired z test
#'
#' @description
#' This function performs power and sample size calculations for a paired z test which is
#' analogous to a paired t test with variance assumed to be known.
#' This function is provided largely for pedagogical purposes; in general,
#' for real studies, the paired t test procedure should be used.
#'
#'
#' @param N The sample size; the number of pairs.
#' @param delta DeltaA (the true mean difference) - Delta0 (the mean difference under the null).
#' @param sd1 The pre standard deviation; defaults to 1.
#' @param sd2 The post standard deviation; defaults to 1.
#' @param rho The correlation between pre and post measurements on the same individual; defaults to 0.
#' @param alpha The significance level (type 1 error rate); defaults to 0.05.
#' @param power The specified level of power.
#' @param sides Either 1 or 2 (default) to specify a one- or two- sided hypothesis test.
#' @param v Either TRUE for verbose output or FALSE (default) to output computed argument only.
#'
#' @return A list of the arguments (including the computed one).
#' @export
#'
#' @examples
#' ztest.paired(N = NULL, delta = 4, sd1 = 10, sd2 = 10, rho = 0.4, power = 0.8, sides = 2)

ztest.paired <- function (N = NULL, delta = NULL,
                          sd1 = 1, sd2 = 1, rho = NULL,
                          alpha = 0.05, power = NULL, sides = 2,
                          v = FALSE) {

  # Check if the arguments are specified correctly
  check.many(list(N, delta, alpha, power), "oneof")
  check.param(N, "pos")
  check.param(sd1, "req"); check.param(sd1, "pos")
  check.param(sd2, "req"); check.param(sd2, "pos")
  check.param(rho, "req"); check.param(rho, "uniti")
  check.param(delta, "num")
  check.param(alpha, "unit")
  check.param(power, "unit")
  check.param(sides, "req"); check.param(sides, "vals", valslist = c(1, 2))
  check.param(v, "req"); check.param(v, "bool")

  # Calculate the standard deviation of differences within pairs
  sigmad <- sqrt(sd1^2 + sd2^2 - 2 * rho * sd1 * sd2)

  # Calculate power
  if (sides == 1)
    p.body <- quote({
      d <- abs(delta) / sigmad
      stats::pnorm(stats::qnorm(alpha) + sqrt(N) * d)
    })

  else if (sides == 2)
    p.body <- quote({
      d <- delta / sigmad
      stats::pnorm(stats::qnorm(alpha / 2) + sqrt(N) * d) +
        stats::pnorm(stats::qnorm(alpha / 2) - sqrt(N) * d)
    })

  NOTE <- "N is the number of pairs"
  if (!v) print(paste("NOTE:", NOTE))

  # Use stats::uniroot function to calculate missing argument
  if (is.null(power)) {
    power <- eval(p.body)
    if (!v) return(power)
  }
  else if (is.null(N)) {
    N <- stats::uniroot(function(N) eval(p.body) - power, c(2, 1e+07))$root
    if (!v) return(N)
  }
  else if (is.null(delta)) {
    delta <- stats::uniroot(function(delta) eval(p.body) - power,  c(1e-07, 1e+07))$root
    if (!v) return(delta)
  }
  else if (is.null(alpha)) {
    alpha <- stats::uniroot(function(alpha) eval(p.body) - power, c(1e-10, 1 - 1e-10))$root
    if (!v) return(alpha)
  }
  else stop("internal error")

  # Generate output text
  METHOD <- "Paired z test power calculation"
  sd <- c(sd1, sd2)

  # Print output as a power.htest object
  structure(list(N = N, delta = delta, `sd1, sd2` = sd, rho = rho,
                 alpha = alpha, power = power, sides = sides,
                 method = METHOD, note = NOTE), class = "power.htest")
}