tern coverage - 95.32%

Files
Source

#' Confidence Intervals for a Difference of Binomials
#'
#' @description `r lifecycle::badge("experimental")`
#'
#' Several confidence intervals for the difference between proportions.
#'
#' @param grp (`factor`)\cr vector assigning observations to one out of two groups
#'   (e.g. reference and treatment group).
#'
#' @name desctools_binom
NULL

#' Recycle List of Parameters
#'
#' This function recycles all supplied elements to the maximal dimension.
#'
#' @param ... (`any`)\cr Elements to recycle.
#'
#' @return A `list`.
#'
#' @keywords internal
#' @noRd
h_recycle <- function(...) {
  lst <- list(...)
  maxdim <- max(lengths(lst))
  res <- lapply(lst, rep, length.out = maxdim)
  attr(res, "maxdim") <- maxdim
  return(res)
}

#' @describeIn desctools_binom Several confidence intervals for the difference between proportions.
#'
#' @return A `matrix` of 3 values:
#'   * `est`: estimate of proportion difference.
#'   * `lwr.ci`: estimate of lower end of the confidence interval.
#'   * `upr.ci`: estimate of upper end of the confidence interval.
#'
#' @examples
#' # Internal function - desctools_binom
#' \dontrun{
#' set.seed(2)
#' rsp <- sample(c(TRUE, FALSE), replace = TRUE, size = 20)
#' grp <- factor(c(rep("A", 10), rep("B", 10)))
#' tbl <- table(grp, factor(rsp, levels = c(TRUE, FALSE)))
#' desctools_binom(
#'   tbl[1], sum(tbl[1], tbl[3]), tbl[2], sum(tbl[2], tbl[4]),
#'   conf.level = 0.90, method = "waldcc"
#' )
#' }
#'
#' @keywords internal
desctools_binom <- function(x1, n1, x2, n2, conf.level = 0.95, sides = c( # nolint
                              "two.sided",
                              "left", "right"
                            ), method = c(
                              "ac", "wald", "waldcc", "score",
                              "scorecc", "mn", "mee", "blj", "ha", "hal", "jp"
                            )) {
  if (missing(sides)) {
    sides <- match.arg(sides)
  }
  if (missing(method)) {
    method <- match.arg(method)
  }
  iBinomDiffCI <- function(x1, n1, x2, n2, conf.level, sides, # nolint
                           method) {
    if (sides != "two.sided") {
      conf.level <- 1 - 2 * (1 - conf.level) # nolint
    }
    alpha <- 1 - conf.level
    kappa <- stats::qnorm(1 - alpha / 2)
    p1_hat <- x1 / n1
    p2_hat <- x2 / n2
    est <- p1_hat - p2_hat
    switch(method,
      wald = {
        vd <- p1_hat * (1 - p1_hat) / n1 + p2_hat * (1 - p2_hat) / n2
        term2 <- kappa * sqrt(vd)
        ci_lwr <- max(-1, est - term2)
        ci_upr <- min(1, est + term2)
      },
      waldcc = {
        vd <- p1_hat * (1 - p1_hat) / n1 + p2_hat * (1 - p2_hat) / n2
        term2 <- kappa * sqrt(vd)
        term2 <- term2 + 0.5 * (1 / n1 + 1 / n2)
        ci_lwr <- max(-1, est - term2)
        ci_upr <- min(1, est + term2)
      },
      ac = {
        n1 <- n1 + 2
        n2 <- n2 + 2
        x1 <- x1 + 1
        x2 <- x2 + 1
        p1_hat <- x1 / n1
        p2_hat <- x2 / n2
        est1 <- p1_hat - p2_hat
        vd <- p1_hat * (1 - p1_hat) / n1 + p2_hat * (1 - p2_hat) / n2
        term2 <- kappa * sqrt(vd)
        ci_lwr <- max(-1, est1 - term2)
        ci_upr <- min(1, est1 + term2)
      },
      exact = {
        ci_lwr <- NA
        ci_upr <- NA
      },
      score = {
        w1 <- desctools_binomci(
          x = x1, n = n1, conf.level = conf.level,
          method = "wilson"
        )
        w2 <- desctools_binomci(
          x = x2, n = n2, conf.level = conf.level,
          method = "wilson"
        )
        l1 <- w1[2]
        u1 <- w1[3]
        l2 <- w2[2]
        u2 <- w2[3]
        ci_lwr <- est - kappa * sqrt(l1 * (1 - l1) / n1 +
          u2 * (1 - u2) / n2)
        ci_upr <- est + kappa * sqrt(u1 * (1 - u1) / n1 +
          l2 * (1 - l2) / n2)
      },
      scorecc = {
        w1 <- desctools_binomci(
          x = x1, n = n1, conf.level = conf.level,
          method = "wilsoncc"
        )
        w2 <- desctools_binomci(
          x = x2, n = n2, conf.level = conf.level,
          method = "wilsoncc"
        )
        l1 <- w1[2]
        u1 <- w1[3]
        l2 <- w2[2]
        u2 <- w2[3]
        ci_lwr <- max(-1, est - sqrt((p1_hat - l1)^2 +
          (u2 - p2_hat)^2))
        ci_upr <- min(1, est + sqrt((u1 - p1_hat)^2 + (p2_hat -
          l2)^2))
      },
      mee = {
        .score <- function(p1, n1, p2, n2, dif) {
          if (dif > 1) dif <- 1
          if (dif < -1) dif <- -1
          diff <- p1 - p2 - dif
          if (abs(diff) == 0) {
            res <- 0
          } else {
            t <- n2 / n1
            a <- 1 + t
            b <- -(1 + t + p1 + t * p2 + dif * (t + 2))
            c <- dif * dif + dif * (2 * p1 + t + 1) + p1 +
              t * p2
            d <- -p1 * dif * (1 + dif)
            v <- (b / a / 3)^3 - b * c / (6 * a * a) + d / a / 2
            if (abs(v) < .Machine$double.eps) v <- 0
            s <- sqrt((b / a / 3)^2 - c / a / 3)
            u <- ifelse(v > 0, 1, -1) * s
            w <- (3.141592654 + acos(v / u^3)) / 3
            p1d <- 2 * u * cos(w) - b / a / 3
            p2d <- p1d - dif
            n <- n1 + n2
            res <- (p1d * (1 - p1d) / n1 + p2d * (1 - p2d) / n2)
          }
          return(sqrt(res))
        }
        pval <- function(delta) {
          z <- (est - delta) / .score(
            p1_hat, n1, p2_hat,
            n2, delta
          )
          2 * min(stats::pnorm(z), 1 - stats::pnorm(z))
        }
        ci_lwr <- max(-1, stats::uniroot(function(delta) {
          pval(delta) -
            alpha
        }, interval = c(-1 + 1e-06, est - 1e-06))$root)
        ci_upr <- min(1, stats::uniroot(function(delta) {
          pval(delta) -
            alpha
        }, interval = c(est + 1e-06, 1 - 1e-06))$root)
      },
      blj = {
        p1_dash <- (x1 + 0.5) / (n1 + 1)
        p2_dash <- (x2 + 0.5) / (n2 + 1)
        vd <- p1_dash * (1 - p1_dash) / n1 + p2_dash * (1 -
          p2_dash) / n2
        term2 <- kappa * sqrt(vd)
        est_dash <- p1_dash - p2_dash
        ci_lwr <- max(-1, est_dash - term2)
        ci_upr <- min(1, est_dash + term2)
      },
      ha = {
        term2 <- 1 / (2 * min(n1, n2)) + kappa * sqrt(p1_hat *
          (1 - p1_hat) / (n1 - 1) + p2_hat * (1 - p2_hat) / (n2 -
            1))
        ci_lwr <- max(-1, est - term2)
        ci_upr <- min(1, est + term2)
      },
      mn = {
        .conf <- function(x1, n1, x2, n2, z, lower = FALSE) {
          p1 <- x1 / n1
          p2 <- x2 / n2
          p_hat <- p1 - p2
          dp <- 1 + ifelse(lower, 1, -1) * p_hat
          i <- 1
          while (i <= 50) {
            dp <- 0.5 * dp
            y <- p_hat + ifelse(lower, -1, 1) * dp
            score <- .score(p1, n1, p2, n2, y)
            if (score < z) {
              p_hat <- y
            }
            if ((dp < 1e-07) || (abs(z - score) < 1e-06)) {
              (break)()
            } else {
              i <- i +
                1
            }
          }
          return(y)
        }
        .score <- function(p1, n1, p2, n2, dif) {
          diff <- p1 - p2 - dif
          if (abs(diff) == 0) {
            res <- 0
          } else {
            t <- n2 / n1
            a <- 1 + t
            b <- -(1 + t + p1 + t * p2 + dif * (t + 2))
            c <- dif * dif + dif * (2 * p1 + t + 1) + p1 +
              t * p2
            d <- -p1 * dif * (1 + dif)
            v <- (b / a / 3)^3 - b * c / (6 * a * a) + d / a / 2
            s <- sqrt((b / a / 3)^2 - c / a / 3)
            u <- ifelse(v > 0, 1, -1) * s
            w <- (3.141592654 + acos(v / u^3)) / 3
            p1d <- 2 * u * cos(w) - b / a / 3
            p2d <- p1d - dif
            n <- n1 + n2
            var <- (p1d * (1 - p1d) / n1 + p2d * (1 - p2d) / n2) *
              n / (n - 1)
            res <- diff^2 / var
          }
          return(res)
        }
        z <- stats::qchisq(conf.level, 1)
        ci_lwr <- max(-1, .conf(x1, n1, x2, n2, z, TRUE))
        ci_upr <- min(1, .conf(x1, n1, x2, n2, z, FALSE))
      },
      beal = {
        a <- p1_hat + p2_hat
        b <- p1_hat - p2_hat
        u <- ((1 / n1) + (1 / n2)) / 4
        v <- ((1 / n1) - (1 / n2)) / 4
        V <- u * ((2 - a) * a - b^2) + 2 * v * (1 - a) * b # nolint
        z <- stats::qchisq(p = 1 - alpha / 2, df = 1)
        A <- sqrt(z * (V + z * u^2 * (2 - a) * a + z * v^2 * (1 - a)^2)) # nolint
        B <- (b + z * v * (1 - a)) / (1 + z * u) # nolint
        ci_lwr <- max(-1, B - A / (1 + z * u))
        ci_upr <- min(1, B + A / (1 + z * u))
      },
      hal = {
        psi <- (p1_hat + p2_hat) / 2
        u <- (1 / n1 + 1 / n2) / 4
        v <- (1 / n1 - 1 / n2) / 4
        z <- kappa
        theta <- ((p1_hat - p2_hat) + z^2 * v * (1 - 2 *
          psi)) / (1 + z^2 * u)
        w <- z / (1 + z^2 * u) * sqrt(u * (4 * psi * (1 - psi) -
          (p1_hat - p2_hat)^2) + 2 * v * (1 - 2 * psi) *
          (p1_hat - p2_hat) + 4 * z^2 * u^2 * (1 - psi) *
          psi + z^2 * v^2 * (1 - 2 * psi)^2)
        c(theta + w, theta - w)
        ci_lwr <- max(-1, theta - w)
        ci_upr <- min(1, theta + w)
      },
      jp = {
        psi <- 0.5 * ((x1 + 0.5) / (n1 + 1) + (x2 + 0.5) / (n2 +
          1))
        u <- (1 / n1 + 1 / n2) / 4
        v <- (1 / n1 - 1 / n2) / 4
        z <- kappa
        theta <- ((p1_hat - p2_hat) + z^2 * v * (1 - 2 *
          psi)) / (1 + z^2 * u)
        w <- z / (1 + z^2 * u) * sqrt(u * (4 * psi * (1 - psi) -
          (p1_hat - p2_hat)^2) + 2 * v * (1 - 2 * psi) *
          (p1_hat - p2_hat) + 4 * z^2 * u^2 * (1 - psi) *
          psi + z^2 * v^2 * (1 - 2 * psi)^2)
        c(theta + w, theta - w)
        ci_lwr <- max(-1, theta - w)
        ci_upr <- min(1, theta + w)
      },
    )
    ci <- c(
      est = est, lwr.ci = min(ci_lwr, ci_upr),
      upr.ci = max(ci_lwr, ci_upr)
    )
    if (sides == "left") {
      ci[3] <- 1
    } else if (sides == "right") {
      ci[2] <- -1
    }
    return(ci)
  }
  method <- match.arg(arg = method, several.ok = TRUE)
  sides <- match.arg(arg = sides, several.ok = TRUE)
  lst <- h_recycle(
    x1 = x1, n1 = n1, x2 = x2, n2 = n2, conf.level = conf.level,
    sides = sides, method = method
  )
  res <- t(sapply(1:attr(lst, "maxdim"), function(i) {
    iBinomDiffCI(
      x1 = lst$x1[i],
      n1 = lst$n1[i], x2 = lst$x2[i], n2 = lst$n2[i], conf.level = lst$conf.level[i],
      sides = lst$sides[i], method = lst$method[i]
    )
  }))
  lgn <- h_recycle(x1 = if (is.null(names(x1))) {
    paste("x1", seq_along(x1), sep = ".")
  } else {
    names(x1)
  }, n1 = if (is.null(names(n1))) {
    paste("n1", seq_along(n1), sep = ".")
  } else {
    names(n1)
  }, x2 = if (is.null(names(x2))) {
    paste("x2", seq_along(x2), sep = ".")
  } else {
    names(x2)
  }, n2 = if (is.null(names(n2))) {
    paste("n2", seq_along(n2), sep = ".")
  } else {
    names(n2)
  }, conf.level = conf.level, sides = sides, method = method)
  xn <- apply(as.data.frame(lgn[sapply(lgn, function(x) {
    length(unique(x)) !=
      1
  })]), 1, paste, collapse = ":")
  rownames(res) <- xn
  return(res)
}

#' @describeIn desctools_binom Compute confidence intervals for binomial proportions.
#'
#' @param x (`count`)\cr number of successes
#' @param n (`count`)\cr number of trials
#' @param conf.level (`proportion`)\cr confidence level, defaults to 0.95.
#' @param sides (`character`)\cr side of the confidence interval to compute. Must be one of "two-sided" (default),
#'   "left", or "right".
#' @param method (`character`)\cr method to use. Can be one out of: "wald", "wilson", "wilsoncc", "agresti-coull",
#'   "jeffreys", "modified wilson", "modified jeffreys", "clopper-pearson", "arcsine", "logit", "witting", "pratt",
#'   "midp", "lik", and "blaker".
#'
#' @return A `matrix` with 3 columns containing:
#'   * `est`: estimate of proportion difference.
#'   * `lwr.ci`: lower end of the confidence interval.
#'   * `upr.ci`: upper end of the confidence interval.
#'
#' @keywords internal
desctools_binomci <- function(x,
                              n,
                              conf.level = 0.95, # nolint
                              sides = c("two.sided", "left", "right"),
                              method = c(
                                "wilson", "wald", "waldcc", "agresti-coull",
                                "jeffreys", "modified wilson", "wilsoncc", "modified jeffreys",
                                "clopper-pearson", "arcsine", "logit", "witting", "pratt",
                                "midp", "lik", "blaker"
                              ),
                              rand = 123,
                              tol = 1e-05) {
  if (missing(method)) {
    method <- "wilson"
  }
  if (missing(sides)) {
    sides <- "two.sided"
  }
  iBinomCI <- function(x, n, conf.level = 0.95, sides = c( # nolint
                         "two.sided",
                         "left", "right"
                       ), method = c(
                         "wilson", "wilsoncc", "wald",
                         "waldcc", "agresti-coull", "jeffreys", "modified wilson",
                         "modified jeffreys", "clopper-pearson", "arcsine", "logit",
                         "witting", "pratt", "midp", "lik", "blaker"
                       ), rand = 123,
                       tol = 1e-05) {
    if (length(x) != 1) {
      stop("'x' has to be of length 1 (number of successes)")
    }
    if (length(n) != 1) {
      stop("'n' has to be of length 1 (number of trials)")
    }
    if (length(conf.level) != 1) {
      stop("'conf.level' has to be of length 1 (confidence level)")
    }
    if (conf.level < 0.5 || conf.level > 1) {
      stop("'conf.level' has to be in [0.5, 1]")
    }
    sides <- match.arg(sides, choices = c(
      "two.sided", "left",
      "right"
    ), several.ok = FALSE)
    if (sides != "two.sided") {
      conf.level <- 1 - 2 * (1 - conf.level) # nolint
    }
    alpha <- 1 - conf.level
    kappa <- stats::qnorm(1 - alpha / 2)
    p_hat <- x / n
    q_hat <- 1 - p_hat
    est <- p_hat
    switch(match.arg(arg = method, choices = c(
      "wilson",
      "wald", "waldcc", "wilsoncc", "agresti-coull", "jeffreys",
      "modified wilson", "modified jeffreys", "clopper-pearson",
      "arcsine", "logit", "witting", "pratt", "midp", "lik",
      "blaker"
    )),
    wald = {
      term2 <- kappa * sqrt(p_hat * q_hat) / sqrt(n)
      ci_lwr <- max(0, p_hat - term2)
      ci_upr <- min(1, p_hat + term2)
    },
    waldcc = {
      term2 <- kappa * sqrt(p_hat * q_hat) / sqrt(n)
      term2 <- term2 + 1 / (2 * n)
      ci_lwr <- max(0, p_hat - term2)
      ci_upr <- min(1, p_hat + term2)
    },
    wilson = {
      term1 <- (x + kappa^2 / 2) / (n + kappa^2)
      term2 <- kappa * sqrt(n) / (n + kappa^2) * sqrt(p_hat *
        q_hat + kappa^2 / (4 * n))
      ci_lwr <- max(0, term1 - term2)
      ci_upr <- min(1, term1 + term2)
    },
    wilsoncc = {
      lci <- (2 * x + kappa^2 - 1 - kappa * sqrt(kappa^2 -
        2 - 1 / n + 4 * p_hat * (n * q_hat + 1))) / (2 *
        (n + kappa^2))
      uci <- (2 * x + kappa^2 + 1 + kappa * sqrt(kappa^2 +
        2 - 1 / n + 4 * p_hat * (n * q_hat - 1))) / (2 *
        (n + kappa^2))
      ci_lwr <- max(0, ifelse(p_hat == 0, 0, lci))
      ci_upr <- min(1, ifelse(p_hat == 1, 1, uci))
    },
    `agresti-coull` = {
      x_tilde <- x + kappa^2 / 2
      n_tilde <- n + kappa^2
      p_tilde <- x_tilde / n_tilde
      q_tilde <- 1 - p_tilde
      est <- p_tilde
      term2 <- kappa * sqrt(p_tilde * q_tilde) / sqrt(n_tilde)
      ci_lwr <- max(0, p_tilde - term2)
      ci_upr <- min(1, p_tilde + term2)
    },
    jeffreys = {
      if (x == 0) {
        ci_lwr <- 0
      } else {
        ci_lwr <- stats::qbeta(
          alpha / 2,
          x + 0.5, n - x + 0.5
        )
      }
      if (x == n) {
        ci_upr <- 1
      } else {
        ci_upr <- stats::qbeta(1 -
          alpha / 2, x + 0.5, n - x + 0.5)
      }
    },
    `modified wilson` = {
      term1 <- (x + kappa^2 / 2) / (n + kappa^2)
      term2 <- kappa * sqrt(n) / (n + kappa^2) * sqrt(p_hat *
        q_hat + kappa^2 / (4 * n))
      if ((n <= 50 & x %in% c(1, 2)) | (n >= 51 & x %in%
        c(1:3))) {
        ci_lwr <- 0.5 * stats::qchisq(alpha, 2 *
          x) / n
      } else {
        ci_lwr <- max(0, term1 - term2)
      }
      if ((n <= 50 & x %in% c(n - 1, n - 2)) | (n >= 51 &
        x %in% c(n - (1:3)))) {
        ci_upr <- 1 - 0.5 * stats::qchisq(
          alpha,
          2 * (n - x)
        ) / n
      } else {
        ci_upr <- min(1, term1 +
          term2)
      }
    },
    `modified jeffreys` = {
      if (x == n) {
        ci_lwr <- (alpha / 2)^(1 / n)
      } else {
        if (x <= 1) {
          ci_lwr <- 0
        } else {
          ci_lwr <- stats::qbeta(
            alpha / 2,
            x + 0.5, n - x + 0.5
          )
        }
      }
      if (x == 0) {
        ci_upr <- 1 - (alpha / 2)^(1 / n)
      } else {
        if (x >= n - 1) {
          ci_upr <- 1
        } else {
          ci_upr <- stats::qbeta(1 -
            alpha / 2, x + 0.5, n - x + 0.5)
        }
      }
    },
    `clopper-pearson` = {
      ci_lwr <- stats::qbeta(alpha / 2, x, n - x + 1)
      ci_upr <- stats::qbeta(1 - alpha / 2, x + 1, n - x)
    },
    arcsine = {
      p_tilde <- (x + 0.375) / (n + 0.75)
      est <- p_tilde
      ci_lwr <- sin(asin(sqrt(p_tilde)) - 0.5 * kappa / sqrt(n))^2
      ci_upr <- sin(asin(sqrt(p_tilde)) + 0.5 * kappa / sqrt(n))^2
    },
    logit = {
      lambda_hat <- log(x / (n - x))
      V_hat <- n / (x * (n - x)) # nolint
      lambda_lower <- lambda_hat - kappa * sqrt(V_hat)
      lambda_upper <- lambda_hat + kappa * sqrt(V_hat)
      ci_lwr <- exp(lambda_lower) / (1 + exp(lambda_lower))
      ci_upr <- exp(lambda_upper) / (1 + exp(lambda_upper))
    },
    witting = {
      set.seed(rand)
      x_tilde <- x + stats::runif(1, min = 0, max = 1)
      pbinom_abscont <- function(q, size, prob) {
        v <- trunc(q)
        term1 <- stats::pbinom(v - 1, size = size, prob = prob)
        term2 <- (q - v) * stats::dbinom(v, size = size, prob = prob)
        return(term1 + term2)
      }
      qbinom_abscont <- function(p, size, x) {
        fun <- function(prob, size, x, p) {
          pbinom_abscont(x, size, prob) - p
        }
        stats::uniroot(fun,
          interval = c(0, 1), size = size,
          x = x, p = p
        )$root
      }
      ci_lwr <- qbinom_abscont(1 - alpha, size = n, x = x_tilde)
      ci_upr <- qbinom_abscont(alpha, size = n, x = x_tilde)
    },
    pratt = {
      if (x == 0) {
        ci_lwr <- 0
        ci_upr <- 1 - alpha^(1 / n)
      } else if (x == 1) {
        ci_lwr <- 1 - (1 - alpha / 2)^(1 / n)
        ci_upr <- 1 - (alpha / 2)^(1 / n)
      } else if (x == (n - 1)) {
        ci_lwr <- (alpha / 2)^(1 / n)
        ci_upr <- (1 - alpha / 2)^(1 / n)
      } else if (x == n) {
        ci_lwr <- alpha^(1 / n)
        ci_upr <- 1
      } else {
        z <- stats::qnorm(1 - alpha / 2)
        A <- ((x + 1) / (n - x))^2 # nolint
        B <- 81 * (x + 1) * (n - x) - 9 * n - 8 # nolint
        C <- (0 - 3) * z * sqrt(9 * (x + 1) * (n - x) * (9 * n + 5 - z^2) + n + 1) # nolint
        D <- 81 * (x + 1)^2 - 9 * (x + 1) * (2 + z^2) + 1 # nolint
        E <- 1 + A * ((B + C) / D)^3 # nolint
        ci_upr <- 1 / E
        A <- (x / (n - x - 1))^2 # nolint
        B <- 81 * x * (n - x - 1) - 9 * n - 8 # nolint
        C <- 3 * z * sqrt(9 * x * (n - x - 1) * (9 * n + 5 - z^2) + n + 1) # nolint
        D <- 81 * x^2 - 9 * x * (2 + z^2) + 1 # nolint
        E <- 1 + A * ((B + C) / D)^3 # nolint
        ci_lwr <- 1 / E
      }
    },
    midp = {
      f_low <- function(pi, x, n) {
        1 / 2 * stats::dbinom(x, size = n, prob = pi) + stats::pbinom(x,
          size = n, prob = pi, lower.tail = FALSE
        ) -
          (1 - conf.level) / 2
      }
      f_up <- function(pi, x, n) {
        1 / 2 * stats::dbinom(x, size = n, prob = pi) + stats::pbinom(x -
          1, size = n, prob = pi) - (1 - conf.level) / 2
      }
      ci_lwr <- 0
      ci_upr <- 1
      if (x != 0) {
        ci_lwr <- stats::uniroot(f_low,
          interval = c(0, p_hat),
          x = x, n = n
        )$root
      }
      if (x != n) {
        ci_upr <- stats::uniroot(f_up, interval = c(
          p_hat,
          1
        ), x = x, n = n)$root
      }
    },
    lik = {
      ci_lwr <- 0
      ci_upr <- 1
      z <- stats::qnorm(1 - alpha * 0.5)
      tol <- .Machine$double.eps^0.5
      BinDev <- function(y, x, mu, wt, bound = 0, tol = .Machine$double.eps^0.5, # nolint
                         ...) {
        ll_y <- ifelse(y %in% c(0, 1), 0, stats::dbinom(x, wt,
          y,
          log = TRUE
        ))
        ll_mu <- ifelse(mu %in% c(0, 1), 0, stats::dbinom(x,
          wt, mu,
          log = TRUE
        ))
        res <- ifelse(abs(y - mu) < tol, 0, sign(y -
          mu) * sqrt(-2 * (ll_y - ll_mu)))
        return(res - bound)
      }
      if (x != 0 && tol < p_hat) {
        ci_lwr <- if (BinDev(
          tol, x, p_hat, n, -z,
          tol
        ) <= 0) {
          stats::uniroot(
            f = BinDev, interval = c(tol, if (p_hat <
              tol || p_hat == 1) {
              1 - tol
            } else {
              p_hat
            }), bound = -z,
            x = x, mu = p_hat, wt = n
          )$root
        }
      }
      if (x != n && p_hat < (1 - tol)) {
        ci_upr <- if (BinDev(y = 1 - tol, x = x, mu = ifelse(p_hat >
          1 - tol, tol, p_hat), wt = n, bound = z, tol = tol) <
          0) {
          ci_lwr <- if (BinDev(
            tol, x, if (p_hat <
              tol || p_hat == 1) {
              1 - tol
            } else {
              p_hat
            }, n,
            -z, tol
          ) <= 0) {
            stats::uniroot(
              f = BinDev, interval = c(tol, p_hat),
              bound = -z, x = x, mu = p_hat, wt = n
            )$root
          }
        } else {
          stats::uniroot(
            f = BinDev, interval = c(if (p_hat >
              1 - tol) {
              tol
            } else {
              p_hat
            }, 1 - tol), bound = z,
            x = x, mu = p_hat, wt = n
          )$root
        }
      }
    },
    blaker = {
      acceptbin <- function(x, n, p) {
        p1 <- 1 - stats::pbinom(x - 1, n, p)
        p2 <- stats::pbinom(x, n, p)
        a1 <- p1 + stats::pbinom(stats::qbinom(p1, n, p) - 1, n, p)
        a2 <- p2 + 1 - stats::pbinom(
          stats::qbinom(1 - p2, n, p), n,
          p
        )
        return(min(a1, a2))
      }
      ci_lwr <- 0
      ci_upr <- 1
      if (x != 0) {
        ci_lwr <- stats::qbeta((1 - conf.level) / 2, x, n -
          x + 1)
        while (acceptbin(x, n, ci_lwr + tol) < (1 -
          conf.level)) {
          ci_lwr <- ci_lwr + tol
        }
      }
      if (x != n) {
        ci_upr <- stats::qbeta(1 - (1 - conf.level) / 2, x +
          1, n - x)
        while (acceptbin(x, n, ci_upr - tol) < (1 -
          conf.level)) {
          ci_upr <- ci_upr - tol
        }
      }
    }
    )
    ci <- c(est = est, lwr.ci = max(0, ci_lwr), upr.ci = min(
      1,
      ci_upr
    ))
    if (sides == "left") {
      ci[3] <- 1
    } else if (sides == "right") {
      ci[2] <- 0
    }
    return(ci)
  }
  lst <- list(
    x = x, n = n, conf.level = conf.level, sides = sides,
    method = method, rand = rand
  )
  maxdim <- max(unlist(lapply(lst, length)))
  lgp <- lapply(lst, rep, length.out = maxdim)
  lgn <- h_recycle(x = if (is.null(names(x))) {
    paste("x", seq_along(x), sep = ".")
  } else {
    names(x)
  }, n = if (is.null(names(n))) {
    paste("n", seq_along(n), sep = ".")
  } else {
    names(n)
  }, conf.level = conf.level, sides = sides, method = method)
  xn <- apply(as.data.frame(lgn[sapply(lgn, function(x) {
    length(unique(x)) !=
      1
  })]), 1, paste, collapse = ":")
  res <- t(sapply(1:maxdim, function(i) {
    iBinomCI(
      x = lgp$x[i],
      n = lgp$n[i], conf.level = lgp$conf.level[i], sides = lgp$sides[i],
      method = lgp$method[i], rand = lgp$rand[i]
    )
  }))
  colnames(res)[1] <- c("est")
  rownames(res) <- xn
  return(res)
}

#' Missing Data
#'
#' @description `r lifecycle::badge("stable")`
#'
#' Substitute missing data with a string or factor level.
#'
#' @param x (`factor` or `character` vector)\cr values for which any missing values should be substituted.
#' @param label (`character`)\cr string that missing data should be replaced with.
#'
#' @return `x` with any `NA` values substituted by `label`.
#'
#' @examples
#' explicit_na(c(NA, "a", "b"))
#' is.na(explicit_na(c(NA, "a", "b")))
#'
#' explicit_na(factor(c(NA, "a", "b")))
#' is.na(explicit_na(factor(c(NA, "a", "b"))))
#'
#' explicit_na(sas_na(c("a", "")))
#'
#' @export
explicit_na <- function(x, label = "<Missing>") {
  checkmate::assert_string(label)

  if (is.factor(x)) {
    x <- forcats::fct_na_value_to_level(x, label)
    forcats::fct_drop(x, only = label)
  } else if (is.character(x)) {
    x[is.na(x)] <- label
    x
  } else {
    stop("only factors and character vectors allowed")
  }
}

#' Convert Strings to `NA`
#'
#' @description `r lifecycle::badge("stable")`
#'
#' SAS imports missing data as empty strings or strings with whitespaces only. This helper function can be used to
#' convert these values to `NA`s.
#'
#' @inheritParams explicit_na
#' @param empty (`logical`)\cr if `TRUE` empty strings get replaced by `NA`.
#' @param whitespaces (`logical`)\cr if `TRUE` then strings made from whitespaces only get replaced with `NA`.
#'
#' @return `x` with `""` and/or whitespace-only values substituted by `NA`, depending on the values of
#'   `empty` and `whitespaces`.
#'
#' @examples
#' sas_na(c("1", "", " ", "   ", "b"))
#' sas_na(factor(c("", " ", "b")))
#'
#' is.na(sas_na(c("1", "", " ", "   ", "b")))
#'
#' @export
sas_na <- function(x, empty = TRUE, whitespaces = TRUE) {
  checkmate::assert_flag(empty)
  checkmate::assert_flag(whitespaces)

  if (is.factor(x)) {
    empty_levels <- levels(x) == ""
    if (empty && any(empty_levels)) levels(x)[empty_levels] <- NA

    ws_levels <- grepl("^\\s+$", levels(x))
    if (whitespaces && any(ws_levels)) levels(x)[ws_levels] <- NA

    x
  } else if (is.character(x)) {
    if (empty) x[x == ""] <- NA_character_

    if (whitespaces) x[grepl("^\\s+$", x)] <- NA_character_

    x
  } else {
    stop("only factors and character vectors allowed")
  }
}

#' Compare Variables Between Groups
#'
#' @description `r lifecycle::badge("stable")`
#'
#' Comparison with a reference group for different `x` objects.
#'
#' @inheritParams argument_convention
#'
#' @note
#' * For factor variables, `denom` for factor proportions can only be `n` since the purpose is to compare proportions
#'   between columns, therefore a row-based proportion would not make sense. Proportion based on `N_col` would
#'   be difficult since we use counts for the chi-squared test statistic, therefore missing values should be accounted
#'   for as explicit factor levels.
#' * If factor variables contain `NA`, these `NA` values are excluded by default. To include `NA` values
#'   set `na.rm = FALSE` and missing values will be displayed as an `NA` level. Alternatively, an explicit
#'   factor level can be defined for `NA` values during pre-processing via [df_explicit_na()] - the
#'   default `na_level` (`"<Missing>"`) will also be excluded when `na.rm` is set to `TRUE`.
#' * For character variables, automatic conversion to factor does not guarantee that the table
#'   will be generated correctly. In particular for sparse tables this very likely can fail.
#'   Therefore it is always better to manually convert character variables to factors during pre-processing.
#' * For `compare_vars()`, the column split must define a reference group via `ref_group` so that the comparison
#'   is well defined.
#'
#' @seealso Relevant constructor function [create_afun_compare()], and [s_summary()] which is used internally
#'   to compute a summary within `s_compare()`.
#'
#' @name compare_variables
#' @include summarize_variables.R
NULL

#' @describeIn compare_variables S3 generic function to produce a comparison summary.
#'
#' @return
#' * `s_compare()` returns output of [s_summary()] and comparisons versus the reference group in the form of p-values.
#'
#' @export
s_compare <- function(x,
                      .ref_group,
                      .in_ref_col,
                      ...) {
  UseMethod("s_compare", x)
}

#' @describeIn compare_variables Method for `numeric` class. This uses the standard t-test
#'   to calculate the p-value.
#'
#' @method s_compare numeric
#'
#' @examples
#' # `s_compare.numeric`
#'
#' ## Usual case where both this and the reference group vector have more than 1 value.
#' s_compare(rnorm(10, 5, 1), .ref_group = rnorm(5, -5, 1), .in_ref_col = FALSE)
#'
#' ## If one group has not more than 1 value, then p-value is not calculated.
#' s_compare(rnorm(10, 5, 1), .ref_group = 1, .in_ref_col = FALSE)
#'
#' ## Empty numeric does not fail, it returns NA-filled items and no p-value.
#' s_compare(numeric(), .ref_group = numeric(), .in_ref_col = FALSE)
#'
#' @export
s_compare.numeric <- function(x,
                              .ref_group,
                              .in_ref_col,
                              ...) {
  checkmate::assert_numeric(x)
  checkmate::assert_numeric(.ref_group)
  checkmate::assert_flag(.in_ref_col)

  y <- s_summary.numeric(x = x, ...)

  y$pval <- if (!.in_ref_col && n_available(x) > 1 && n_available(.ref_group) > 1) {
    stats::t.test(x, .ref_group)$p.value
  } else {
    character()
  }

  y
}

#' @describeIn compare_variables Method for `factor` class. This uses the chi-squared test
#'   to calculate the p-value.
#'
#' @param denom (`string`)\cr choice of denominator for factor proportions,
#'   can only be `n` (number of values in this row and column intersection).
#'
#' @method s_compare factor
#'
#' @examples
#' # `s_compare.factor`
#'
#' ## Basic usage:
#' x <- factor(c("a", "a", "b", "c", "a"))
#' y <- factor(c("a", "b", "c"))
#' s_compare(x = x, .ref_group = y, .in_ref_col = FALSE)
#'
#' ## Management of NA values.
#' x <- explicit_na(factor(c("a", "a", "b", "c", "a", NA, NA)))
#' y <- explicit_na(factor(c("a", "b", "c", NA)))
#' s_compare(x = x, .ref_group = y, .in_ref_col = FALSE, na.rm = TRUE)
#' s_compare(x = x, .ref_group = y, .in_ref_col = FALSE, na.rm = FALSE)
#'
#' @export
s_compare.factor <- function(x,
                             .ref_group,
                             .in_ref_col,
                             denom = "n",
                             na.rm = TRUE, # nolint
                             ...) {
  checkmate::assert_flag(.in_ref_col)
  assert_valid_factor(x)
  assert_valid_factor(.ref_group)
  denom <- match.arg(denom)

  y <- s_summary.factor(
    x = x,
    denom = denom,
    na.rm = na.rm,
    ...
  )

  if (na.rm) {
    x <- x[!is.na(x)] %>% fct_discard("<Missing>")
    .ref_group <- .ref_group[!is.na(.ref_group)] %>% fct_discard("<Missing>")
  } else {
    x <- x %>% explicit_na(label = "NA")
    .ref_group <- .ref_group %>% explicit_na(label = "NA")
  }

  checkmate::assert_factor(x, levels = levels(.ref_group), min.levels = 2)

  y$pval <- if (!.in_ref_col && length(x) > 0 && length(.ref_group) > 0) {
    tab <- rbind(table(x), table(.ref_group))
    res <- suppressWarnings(stats::chisq.test(tab))
    res$p.value
  } else {
    character()
  }

  y
}

#' @describeIn compare_variables Method for `character` class. This makes an automatic
#'   conversion to `factor` (with a warning) and then forwards to the method for factors.
#'
#' @param verbose (`logical`)\cr Whether warnings and messages should be printed. Mainly used
#'   to print out information about factor casting. Defaults to `TRUE`.
#'
#' @method s_compare character
#'
#' @examples
#' # `s_compare.character`
#'
#' ## Basic usage:
#' x <- c("a", "a", "b", "c", "a")
#' y <- c("a", "b", "c")
#' s_compare(x, .ref_group = y, .in_ref_col = FALSE, .var = "x", verbose = FALSE)
#'
#' ## Note that missing values handling can make a large difference:
#' x <- c("a", "a", "b", "c", "a", NA)
#' y <- c("a", "b", "c", rep(NA, 20))
#' s_compare(x,
#'   .ref_group = y, .in_ref_col = FALSE,
#'   .var = "x", verbose = FALSE
#' )
#' s_compare(x,
#'   .ref_group = y, .in_ref_col = FALSE, .var = "x",
#'   na.rm = FALSE, verbose = FALSE
#' )
#'
#' @export
s_compare.character <- function(x,
                                .ref_group,
                                .in_ref_col,
                                denom = "n",
                                na.rm = TRUE, # nolint
                                .var,
                                verbose = TRUE,
                                ...) {
  x <- as_factor_keep_attributes(x, x_name = .var, verbose = verbose)
  .ref_group <- as_factor_keep_attributes(.ref_group, x_name = .var, verbose = verbose)
  s_compare(
    x = x,
    .ref_group = .ref_group,
    .in_ref_col = .in_ref_col,
    denom = denom,
    na.rm = na.rm,
    ...
  )
}

#' @describeIn compare_variables Method for `logical` class. A chi-squared test
#'   is used. If missing values are not removed, then they are counted as `FALSE`.
#'
#' @method s_compare logical
#'
#' @examples
#' # `s_compare.logical`
#'
#' ## Basic usage:
#' x <- c(TRUE, FALSE, TRUE, TRUE)
#' y <- c(FALSE, FALSE, TRUE)
#' s_compare(x, .ref_group = y, .in_ref_col = FALSE)
#'
#' ## Management of NA values.
#' x <- c(NA, TRUE, FALSE)
#' y <- c(NA, NA, NA, NA, FALSE)
#' s_compare(x, .ref_group = y, .in_ref_col = FALSE, na.rm = TRUE)
#' s_compare(x, .ref_group = y, .in_ref_col = FALSE, na.rm = FALSE)
#'
#' @export
s_compare.logical <- function(x,
                              .ref_group,
                              .in_ref_col,
                              na.rm = TRUE, # nolint
                              denom = "n",
                              ...) {
  denom <- match.arg(denom)

  y <- s_summary.logical(
    x = x,
    na.rm = na.rm,
    denom = denom,
    ...
  )

  if (na.rm) {
    x <- stats::na.omit(x)
    .ref_group <- stats::na.omit(.ref_group)
  } else {
    x[is.na(x)] <- FALSE
    .ref_group[is.na(.ref_group)] <- FALSE
  }

  y$pval <- if (!.in_ref_col && length(x) > 0 && length(.ref_group) > 0) {
    x <- factor(x, levels = c(TRUE, FALSE))
    .ref_group <- factor(.ref_group, levels = c(TRUE, FALSE))
    tbl <- rbind(table(x), table(.ref_group))
    suppressWarnings(prop_chisq(tbl))
  } else {
    character()
  }

  y
}

#' @describeIn compare_variables Formatted analysis function which is used as `afun`
#'   in `compare_vars()`.
#'
#' @return
#' * `a_compare()` returns the corresponding list with formatted [rtables::CellValue()].
#'
#' @export
a_compare <- function(x,
                      .ref_group,
                      .in_ref_col,
                      ...,
                      .var) {
  UseMethod("a_compare", x)
}

#' @describeIn compare_variables Formatted analysis function method for `numeric` class.
#'
#' @examples
#' # `a_compare.numeric`
#' a_compare(
#'   rnorm(10, 5, 1),
#'   .ref_group = rnorm(20, -5, 1),
#'   .in_ref_col = FALSE,
#'   .var = "bla"
#' )
#'
#' @export
a_compare.numeric <- make_afun(
  s_compare.numeric,
  .formats = c(
    .a_summary_numeric_formats,
    pval = "x.xxxx | (<0.0001)"
  ),
  .labels = c(
    .a_summary_numeric_labels,
    pval = "p-value (t-test)"
  ),
  .null_ref_cells = FALSE
)

.a_compare_counts_formats <- c(
  .a_summary_counts_formats,
  pval = "x.xxxx | (<0.0001)"
)

.a_compare_counts_labels <- c(
  pval = "p-value (chi-squared test)"
)

#' @describeIn compare_variables Formatted analysis function method for `factor` class.
#'
#' @examples
#' # `a_compare.factor`
#' # We need to ungroup `count` and `count_fraction` first so that the `rtables` formatting
#' # functions can be applied correctly.
#' afun <- make_afun(
#'   getS3method("a_compare", "factor"),
#'   .ungroup_stats = c("count", "count_fraction")
#' )
#' x <- factor(c("a", "a", "b", "c", "a"))
#' y <- factor(c("a", "a", "b", "c"))
#' afun(x, .ref_group = y, .in_ref_col = FALSE)
#'
#' @export
a_compare.factor <- make_afun(
  s_compare.factor,
  .formats = .a_compare_counts_formats,
  .labels = .a_compare_counts_labels,
  .null_ref_cells = FALSE
)

#' @describeIn compare_variables Formatted analysis function method for `character` class.
#'
#' @examples
#' # `a_compare.character`
#' afun <- make_afun(
#'   getS3method("a_compare", "character"),
#'   .ungroup_stats = c("count", "count_fraction")
#' )
#' x <- c("A", "B", "A", "C")
#' y <- c("B", "A", "C")
#' afun(x, .ref_group = y, .in_ref_col = FALSE, .var = "x", verbose = FALSE)
#'
#' @export
a_compare.character <- make_afun(
  s_compare.character,
  .formats = .a_compare_counts_formats,
  .labels = .a_compare_counts_labels,
  .null_ref_cells = FALSE
)

#' @describeIn compare_variables Formatted analysis function method for `logical` class.
#'
#' @examples
#' # `a_compare.logical`
#' afun <- make_afun(
#'   getS3method("a_compare", "logical")
#' )
#' x <- c(TRUE, FALSE, FALSE, TRUE, TRUE)
#' y <- c(TRUE, FALSE)
#' afun(x, .ref_group = y, .in_ref_col = FALSE)
#'
#' @export
a_compare.logical <- make_afun(
  s_compare.logical,
  .formats = .a_compare_counts_formats,
  .labels = .a_compare_counts_labels,
  .null_ref_cells = FALSE
)

#' Constructor Function for [compare_vars()]
#'
#' @description `r lifecycle::badge("stable")`
#'
#' Constructor function which creates a combined formatted analysis function.
#'
#' @inheritParams argument_convention
#' @param .indent_mods (named `vector` of `integer`)\cr indent modifiers for the labels. Each element of the vector
#'   should be a name-value pair with name corresponding to a statistic specified in `.stats` and value the indentation
#'   for that statistic's row label.
#'
#' @return Combined formatted analysis function for use in [compare_vars()].
#'
#' @note Since [a_compare()] is generic and we want customization of the formatting arguments
#'   via [rtables::make_afun()], we need to create another temporary generic function, with
#'   corresponding customized methods. Then in order for the methods to be found,
#'   we need to wrap them in a combined `afun`. Since this is required by two layout creating
#'   functions (and possibly others in the future), we provide a constructor that does this:
#'   [create_afun_compare()].
#'
#' @seealso [compare_vars()]
#'
#' @examples
#' # `create_afun_compare()` to create combined `afun`
#'
#' afun <- create_afun_compare(
#'   .stats = c("n", "count_fraction", "mean_sd", "pval"),
#'   .indent_mods = c(pval = 1L)
#' )
#'
#' lyt <- basic_table() %>%
#'   split_cols_by("ARMCD", ref_group = "ARM A") %>%
#'   analyze(
#'     "AGE",
#'     afun = afun,
#'     show_labels = "visible"
#'   )
#' build_table(lyt, df = tern_ex_adsl)
#'
#' lyt <- basic_table() %>%
#'   split_cols_by("ARMCD", ref_group = "ARM A") %>%
#'   analyze(
#'     "SEX",
#'     afun = afun,
#'     show_labels = "visible"
#'   )
#' build_table(lyt, df = tern_ex_adsl)
#'
#' @export
create_afun_compare <- function(.stats = NULL,
                                .formats = NULL,
                                .labels = NULL,
                                .indent_mods = NULL) {
  function(x,
           .ref_group,
           .in_ref_col,
           ...,
           .var) {
    afun <- function(x, ...) {
      UseMethod("afun", x)
    }

    numeric_stats <- afun_selected_stats(
      .stats,
      all_stats = c(names(.a_summary_numeric_formats), "pval")
    )
    afun.numeric <- make_afun( # nolint
      a_compare.numeric,
      .stats = numeric_stats,
      .formats = extract_by_name(.formats, numeric_stats),
      .labels = extract_by_name(.labels, numeric_stats),
      .indent_mods = extract_by_name(.indent_mods, numeric_stats),
      .null_ref_cells = FALSE
    )

    factor_stats <- afun_selected_stats(
      .stats,
      all_stats = names(.a_compare_counts_formats)
    )
    ungroup_stats <- afun_selected_stats(.stats, c("count", "count_fraction"))
    afun.factor <- make_afun( # nolint
      a_compare.factor,
      .stats = factor_stats,
      .formats = extract_by_name(.formats, factor_stats),
      .labels = extract_by_name(.labels, factor_stats),
      .indent_mods = extract_by_name(.indent_mods, factor_stats),
      .ungroup_stats = ungroup_stats,
      .null_ref_cells = FALSE
    )

    afun.character <- make_afun( # nolint
      a_compare.character,
      .stats = factor_stats,
      .formats = extract_by_name(.formats, factor_stats),
      .labels = extract_by_name(.labels, factor_stats),
      .indent_mods = extract_by_name(.indent_mods, factor_stats),
      .ungroup_stats = ungroup_stats,
      .null_ref_cells = FALSE
    )

    afun.logical <- make_afun( # nolint
      a_compare.logical,
      .stats = factor_stats,
      .formats = extract_by_name(.formats, factor_stats),
      .labels = extract_by_name(.labels, factor_stats),
      .indent_mods = extract_by_name(.indent_mods, factor_stats),
      .null_ref_cells = FALSE
    )

    afun(
      x = x,
      .ref_group = .ref_group,
      .in_ref_col = .in_ref_col,
      ...,
      .var = .var
    )
  }
}

#' @describeIn compare_variables Layout-creating function which can take statistics function arguments
#'   and additional format arguments. This function is a wrapper for [rtables::analyze()].
#'
#' @param ... arguments passed to `s_compare()`.
#' @param .indent_mods (named `vector` of `integer`)\cr indent modifiers for the labels. Each element of the vector
#'   should be a name-value pair with name corresponding to a statistic specified in `.stats` and value the indentation
#'   for that statistic's row label.
#'
#' @return
#' * `compare_vars()` returns a layout object suitable for passing to further layouting functions,
#'   or to [rtables::build_table()]. Adding this function to an `rtable` layout will add formatted rows containing
#'   the statistics from `s_compare()` to the table layout.
#'
#' @examples
#' # `compare_vars()` in `rtables` pipelines
#'
#' ## Default output within a `rtables` pipeline.
#' lyt <- basic_table() %>%
#'   split_cols_by("ARMCD", ref_group = "ARM B") %>%
#'   compare_vars(c("AGE", "SEX"))
#' build_table(lyt, tern_ex_adsl)
#'
#' ## Select and format statistics output.
#' lyt <- basic_table() %>%
#'   split_cols_by("ARMCD", ref_group = "ARM C") %>%
#'   compare_vars(
#'     vars = "AGE",
#'     .stats = c("mean_sd", "pval"),
#'     .formats = c(mean_sd = "xx.x, xx.x"),
#'     .labels = c(mean_sd = "Mean, SD")
#'   )
#' build_table(lyt, df = tern_ex_adsl)
#'
#' @export
compare_vars <- function(lyt,
                         vars,
                         var_labels = vars,
                         nested = TRUE,
                         ...,
                         na_level = NA_character_,
                         show_labels = "default",
                         table_names = vars,
                         .stats = c("n", "mean_sd", "count_fraction", "pval"),
                         .formats = NULL,
                         .labels = NULL,
                         .indent_mods = NULL) {
  afun <- create_afun_compare(.stats, .formats, .labels, .indent_mods)

  analyze(
    lyt = lyt,
    vars = vars,
    var_labels = var_labels,
    afun = afun,
    nested = nested,
    extra_args = list(...),
    na_str = na_level,
    inclNAs = TRUE,
    show_labels = show_labels,
    table_names = table_names
  )
}

#' Summary numeric variables in columns
#'
#' @description `r lifecycle::badge("experimental")`
#'
#' Layout-creating function which can be used for creating column-wise summary tables.
#' This function sets the analysis methods as column labels and is a wrapper for
#' [rtables::analyze_colvars()]. It was designed principally for PK tables.
#'
#' @inheritParams argument_convention
#' @inheritParams rtables::analyze_colvars
#' @param row_labels (`character`)\cr as this function works in columns space, usual `.labels`
#'   character vector applies on the column space. You can change the row labels by defining this
#'   parameter to a named character vector with names corresponding to the split values. It defaults
#'   to `NULL` and if it contains only one `string`, it will duplicate that as a row label.
#' @param do_summarize_row_groups (`flag`)\cr defaults to `FALSE` and applies the analysis to the current
#'   label rows. This is a wrapper of [rtables::summarize_row_groups()] and it can accept `labelstr`
#'   to define row labels. This behavior is not supported as we never need to overload row labels.
#' @param split_col_vars (`flag`)\cr defaults to `TRUE` and puts the analysis results onto the columns.
#'   This option allows you to add multiple instances of this functions, also in a nested fashion,
#'   without adding more splits. This split must happen only one time on a single layout.
#'
#' @return
#' A layout object suitable for passing to further layouting functions, or to [rtables::build_table()].
#' Adding this function to an `rtable` layout will summarize the given variables, arrange the output
#' in columns, and add it to the table layout.
#'
#' @note This is an experimental implementation of [rtables::summarize_row_groups()] and
#'   [rtables::analyze_colvars()] that may be subjected to changes as `rtables` extends its
#'   support to more complex analysis pipelines on the column space. For the same reasons,
#'   we encourage to read the examples carefully and file issues for cases that differ from
#'   them.
#'
#'   Here `labelstr` behaves differently than usual. If it is not defined (default as `NULL`),
#'   row labels are assigned automatically to the split values in case of `rtables::analyze_colvars`
#'   (`do_summarize_row_groups = FALSE`, the default), and to the group label for
#'   `do_summarize_row_groups = TRUE`.
#'
#' @seealso [summarize_vars()], [rtables::analyze_colvars()].
#'
#' @examples
#' library(dplyr)
#'
#' # Data preparation
#' adpp <- tern_ex_adpp %>% h_pkparam_sort()
#'
#' lyt <- basic_table() %>%
#'   split_rows_by(var = "STRATA1", label_pos = "topleft") %>%
#'   split_rows_by(
#'     var = "SEX",
#'     label_pos = "topleft",
#'     child_label = "hidden"
#'   ) %>% # Removes duplicated labels
#'   analyze_vars_in_cols(vars = "AGE")
#' result <- build_table(lyt = lyt, df = adpp)
#' result
#'
#' # By selecting just some statistics and ad-hoc labels
#' lyt <- basic_table() %>%
#'   split_rows_by(var = "ARM", label_pos = "topleft") %>%
#'   split_rows_by(
#'     var = "SEX",
#'     label_pos = "topleft",
#'     child_labels = "hidden",
#'     split_fun = drop_split_levels
#'   ) %>%
#'   analyze_vars_in_cols(
#'     vars = "AGE",
#'     .stats = c("n", "cv", "geom_mean"),
#'     .labels = c(
#'       n = "aN",
#'       cv = "aCV",
#'       geom_mean = "aGeomMean"
#'     )
#'   )
#' result <- build_table(lyt = lyt, df = adpp)
#' result
#'
#' # Changing row labels
#' lyt <- basic_table() %>%
#'   analyze_vars_in_cols(
#'     vars = "AGE",
#'     row_labels = "some custom label"
#'   )
#' result <- build_table(lyt, df = adpp)
#' result
#'
#' # Pharmacokinetic parameters
#' lyt <- basic_table() %>%
#'   split_rows_by(
#'     var = "TLG_DISPLAY",
#'     split_label = "PK Parameter",
#'     label_pos = "topleft",
#'     child_label = "hidden"
#'   ) %>%
#'   analyze_vars_in_cols(
#'     vars = "AVAL"
#'   )
#' result <- build_table(lyt, df = adpp)
#' result
#'
#' # Multiple calls (summarize label and analyze underneath)
#' lyt <- basic_table() %>%
#'   split_rows_by(
#'     var = "TLG_DISPLAY",
#'     split_label = "PK Parameter",
#'     label_pos = "topleft"
#'   ) %>%
#'   analyze_vars_in_cols(
#'     vars = "AVAL",
#'     do_summarize_row_groups = TRUE # does a summarize level
#'   ) %>%
#'   split_rows_by("SEX",
#'     child_label = "hidden",
#'     label_pos = "topleft"
#'   ) %>%
#'   analyze_vars_in_cols(
#'     vars = "AVAL",
#'     split_col_vars = FALSE # avoids re-splitting the columns
#'   )
#' result <- build_table(lyt, df = adpp)
#' result
#'
#' @export
analyze_vars_in_cols <- function(lyt,
                                 vars,
                                 ...,
                                 .stats = c(
                                   "n",
                                   "mean",
                                   "sd",
                                   "se",
                                   "cv",
                                   "geom_cv"
                                 ),
                                 .labels = c(
                                   n = "n",
                                   mean = "Mean",
                                   sd = "SD",
                                   se = "SE",
                                   cv = "CV (%)",
                                   geom_cv = "CV % Geometric Mean"
                                 ),
                                 row_labels = NULL,
                                 do_summarize_row_groups = FALSE,
                                 split_col_vars = TRUE,
                                 .indent_mods = NULL,
                                 nested = TRUE,
                                 na_level = NULL,
                                 .formats = NULL) {
  checkmate::assert_string(na_level, null.ok = TRUE)
  checkmate::assert_character(row_labels, null.ok = TRUE)
  checkmate::assert_int(.indent_mods, null.ok = TRUE)
  checkmate::assert_flag(nested)
  checkmate::assert_flag(split_col_vars)
  checkmate::assert_flag(do_summarize_row_groups)

  # Automatic assignment of formats
  if (is.null(.formats)) {
    # General values
    sf_numeric <- summary_formats("numeric")
    sf_counts <- summary_formats("counts")[-1]
    formats_v <- c(sf_numeric, sf_counts)
  } else {
    formats_v <- .formats
  }

  # Check for vars in the case that one or more are used
  if (length(vars) == 1) {
    vars <- rep(vars, length(.stats))
  } else if (length(vars) != length(.stats)) {
    stop(
      "Analyzed variables (vars) does not have the same ",
      "number of elements of specified statistics (.stats)."
    )
  }

  if (split_col_vars) {
    # Checking there is not a previous identical column split
    clyt <- tail(clayout(lyt), 1)[[1]]

    dummy_lyt <- split_cols_by_multivar(
      lyt = basic_table(),
      vars = vars,
      varlabels = .labels[.stats]
    )

    if (any(sapply(clyt, identical, y = get_last_col_split(dummy_lyt)))) {
      stop(
        "Column split called again with the same values. ",
        "This can create many unwanted columns. Please consider adding ",
        "split_col_vars = FALSE to the last call of ",
        deparse(sys.calls()[[sys.nframe() - 1]]), "."
      )
    }

    # Main col split
    lyt <- split_cols_by_multivar(
      lyt = lyt,
      vars = vars,
      varlabels = .labels[.stats]
    )
  }

  if (do_summarize_row_groups) {
    if (length(unique(vars)) > 1) {
      stop("When using do_summarize_row_groups only one label level var should be inserted.")
    }

    # Function list for do_summarize_row_groups. Slightly different handling of labels
    cfun_list <- Map(
      function(stat) {
        function(u, .spl_context, labelstr, ...) {
          # Statistic
          res <- s_summary(u, ...)[[stat]]

          # Label check and replacement
          if (length(row_labels) > 1) {
            if (!(labelstr %in% names(row_labels))) {
              stop(
                "Replacing the labels in do_summarize_row_groups needs a named vector",
                "that contains the split values. In the current split variable ",
                .spl_context$split[nrow(.spl_context)],
                " the labelstr value (split value by default) ", labelstr, " is not in",
                " row_labels names: ", names(row_labels)
              )
            }
            lbl <- unlist(row_labels[labelstr])
          } else {
            lbl <- labelstr
          }

          # Cell creation
          rcell(res,
            label = lbl,
            format = formats_v[names(formats_v) == stat][[1]],
            format_na_str = na_level,
            indent_mod = ifelse(is.null(.indent_mods), 0L, .indent_mods)
          )
        }
      },
      stat = .stats
    )

    # Main call to rtables
    summarize_row_groups(
      lyt = lyt,
      var = unique(vars),
      cfun = cfun_list,
      extra_args = list(...)
    )
  } else {
    # Function list for analyze_colvars
    afun_list <- Map(
      function(stat) {
        function(u, .spl_context, ...) {
          # Main statistics
          res <- s_summary(u, ...)[[stat]]

          # Label from context
          label_from_context <- .spl_context$value[nrow(.spl_context)]

          # Label switcher
          if (is.null(row_labels)) {
            lbl <- label_from_context
          } else {
            if (length(row_labels) > 1) {
              if (!(label_from_context %in% names(row_labels))) {
                stop(
                  "Replacing the labels in do_summarize_row_groups needs a named vector",
                  "that contains the split values. In the current split variable ",
                  .spl_context$split[nrow(.spl_context)],
                  " the split value ", label_from_context, " is not in",
                  " row_labels names: ", names(row_labels)
                )
              }
              lbl <- unlist(row_labels[label_from_context])
            } else {
              lbl <- row_labels
            }
          }

          # Cell creation
          rcell(res,
            label = lbl,
            format = formats_v[names(formats_v) == stat][[1]],
            format_na_str = na_level,
            indent_mod = ifelse(is.null(.indent_mods), 0L, .indent_mods)
          )
        }
      },
      stat = .stats
    )

    # Main call to rtables
    analyze_colvars(lyt,
      afun = afun_list,
      nested = nested,
      extra_args = list(...)
    )
  }
}

# Help function
get_last_col_split <- function(lyt) {
  tail(tail(clayout(lyt), 1)[[1]], 1)[[1]]
}

#' Controls for Cox Regression
#'
#' @description `r lifecycle::badge("stable")`
#'
#' Sets a list of parameters for Cox regression fit. Used internally.
#'
#' @inheritParams argument_convention
#' @param pval_method (`string`)\cr the method used for estimation of p.values; `wald` (default) or `likelihood`.
#' @param interaction (`flag`)\cr if `TRUE`, the model includes the interaction between the studied
#'   treatment and candidate covariate. Note that for univariate models without treatment arm, and
#'   multivariate models, no interaction can be used so that this needs to be `FALSE`.
#' @param ties (`string`)\cr among `exact` (equivalent to `DISCRETE` in SAS), `efron` and `breslow`,
#'   see [survival::coxph()]. Note: there is no equivalent of SAS `EXACT` method in R.
#'
#' @return A `list` of items with names corresponding to the arguments.
#'
#' @seealso [fit_coxreg_univar()] and [fit_coxreg_multivar()].
#'
#' @examples
#' control_coxreg()
#'
#' @export
control_coxreg <- function(pval_method = c("wald", "likelihood"),
                           ties = c("exact", "efron", "breslow"),
                           conf_level = 0.95,
                           interaction = FALSE) {
  pval_method <- match.arg(pval_method)
  ties <- match.arg(ties)
  checkmate::assert_flag(interaction)
  assert_proportion_value(conf_level)
  list(
    pval_method = pval_method,
    ties = ties,
    conf_level = conf_level,
    interaction = interaction
  )
}

#' Custom Tidy Methods for Cox Regression
#'
#' @description `r lifecycle::badge("stable")`
#'
#' @inheritParams argument_convention
#' @param x (`list`)\cr Result of the Cox regression model fitted by [fit_coxreg_univar()] (for univariate models)
#'   or [fit_coxreg_multivar()] (for multivariate models).
#'
#' @return [tidy()] returns:
#' * For `summary.coxph` objects,  a `data.frame` with columns: `Pr(>|z|)`, `exp(coef)`, `exp(-coef)`, `lower .95`,
#'   `upper .95`, `level`, and `n`.
#' * For `coxreg.univar` objects, a `data.frame` with columns: `effect`, `term`, `term_label`, `level`, `n`, `hr`,
#'   `lcl`, `ucl`, `pval`, and `ci`.
#' * For `coxreg.multivar` objects, a `data.frame` with columns: `term`, `pval`, `term_label`, `hr`, `lcl`, `ucl`,
#'   `level`, and `ci`.
#'
#' @seealso [cox_regression]
#'
#' @name tidy_coxreg
NULL

#' @describeIn tidy_coxreg Custom tidy method for [survival::coxph()] summary results.
#'
#' Tidy the [survival::coxph()] results into a `data.frame` to extract model results.
#'
#' @method tidy summary.coxph
#'
#' @examples
#' library(survival)
#' library(broom)
#'
#' set.seed(1, kind = "Mersenne-Twister")
#'
#' dta_bladder <- with(
#'   data = bladder[bladder$enum < 5, ],
#'   data.frame(
#'     time = stop,
#'     status = event,
#'     armcd = as.factor(rx),
#'     covar1 = as.factor(enum),
#'     covar2 = factor(
#'       sample(as.factor(enum)),
#'       levels = 1:4, labels = c("F", "F", "M", "M")
#'     )
#'   )
#' )
#' labels <- c("armcd" = "ARM", "covar1" = "A Covariate Label", "covar2" = "Sex (F/M)")
#' formatters::var_labels(dta_bladder)[names(labels)] <- labels
#' dta_bladder$age <- sample(20:60, size = nrow(dta_bladder), replace = TRUE)
#'
#' formula <- "survival::Surv(time, status) ~ armcd + covar1"
#' msum <- summary(coxph(stats::as.formula(formula), data = dta_bladder))
#' tidy(msum)
#'
#' @export
tidy.summary.coxph <- function(x, # nolint
                               ...) {
  checkmate::assert_class(x, "summary.coxph")
  pval <- x$coefficients
  confint <- x$conf.int
  levels <- rownames(pval)

  pval <- tibble::as_tibble(pval)
  confint <- tibble::as_tibble(confint)

  ret <- cbind(pval[, grepl("Pr", names(pval))], confint)
  ret$level <- levels
  ret$n <- x[["n"]]
  ret
}

#' @describeIn tidy_coxreg Custom tidy method for a univariate Cox regression.
#'
#' Tidy up the result of a Cox regression model fitted by [fit_coxreg_univar()].
#'
#' @method tidy coxreg.univar
#'
#' @examples
#' ## Cox regression: arm + 1 covariate.
#' mod1 <- fit_coxreg_univar(
#'   variables = list(
#'     time = "time", event = "status", arm = "armcd",
#'     covariates = "covar1"
#'   ),
#'   data = dta_bladder,
#'   control = control_coxreg(conf_level = 0.91)
#' )
#'
#' ## Cox regression: arm + 1 covariate + interaction, 2 candidate covariates.
#' mod2 <- fit_coxreg_univar(
#'   variables = list(
#'     time = "time", event = "status", arm = "armcd",
#'     covariates = c("covar1", "covar2")
#'   ),
#'   data = dta_bladder,
#'   control = control_coxreg(conf_level = 0.91, interaction = TRUE)
#' )
#'
#' tidy(mod1)
#' tidy(mod2)
#'
#' @export
tidy.coxreg.univar <- function(x, # nolint
                               ...) {
  checkmate::assert_class(x, "coxreg.univar")
  mod <- x$mod
  vars <- c(x$vars$arm, x$vars$covariates)
  has_arm <- "arm" %in% names(x$vars)

  result <- if (!has_arm) {
    Map(
      mod = mod, vars = vars,
      f = function(mod, vars) {
        h_coxreg_multivar_extract(
          var = vars,
          data = x$data,
          mod = mod,
          control = x$control
        )
      }
    )
  } else if (x$control$interaction) {
    Map(
      mod = mod, covar = vars,
      f = function(mod, covar) {
        h_coxreg_extract_interaction(
          effect = x$vars$arm, covar = covar, mod = mod, data = x$data,
          at = x$at, control = x$control
        )
      }
    )
  } else {
    Map(
      mod = mod, vars = vars,
      f = function(mod, vars) {
        h_coxreg_univar_extract(
          effect = x$vars$arm, covar = vars, data = x$data, mod = mod,
          control = x$control
        )
      }
    )
  }
  result <- do.call(rbind, result)

  result$ci <- Map(lcl = result$lcl, ucl = result$ucl, f = function(lcl, ucl) c(lcl, ucl))
  result$n <- lapply(result$n, empty_vector_if_na)
  result$ci <- lapply(result$ci, empty_vector_if_na)
  result$hr <- lapply(result$hr, empty_vector_if_na)
  if (x$control$interaction) {
    result$pval_inter <- lapply(result$pval_inter, empty_vector_if_na)
    # Remove interaction p-values due to change in specifications.
    result$pval[result$effect != "Treatment:"] <- NA
  }
  result$pval <- lapply(result$pval, empty_vector_if_na)
  attr(result, "conf_level") <- x$control$conf_level
  result
}

#' @describeIn tidy_coxreg Custom tidy method for a multivariate Cox regression.
#'
#' Tidy up the result of a Cox regression model fitted by [fit_coxreg_multivar()].
#'
#' @method tidy coxreg.multivar
#'
#' @examples
#' multivar_model <- fit_coxreg_multivar(
#'   variables = list(
#'     time = "time", event = "status", arm = "armcd",
#'     covariates = c("covar1", "covar2")
#'   ),
#'   data = dta_bladder
#' )
#' broom::tidy(multivar_model)
#'
#' @export
tidy.coxreg.multivar <- function(x, # nolint
                                 ...) {
  checkmate::assert_class(x, "coxreg.multivar")
  vars <- c(x$vars$arm, x$vars$covariates)

  # Convert the model summaries to data.
  result <- Map(
    vars = vars,
    f = function(vars) {
      h_coxreg_multivar_extract(
        var = vars, data = x$data,
        mod = x$mod, control = x$control
      )
    }
  )
  result <- do.call(rbind, result)

  result$ci <- Map(lcl = result$lcl, ucl = result$ucl, f = function(lcl, ucl) c(lcl, ucl))
  result$ci <- lapply(result$ci, empty_vector_if_na)
  result$hr <- lapply(result$hr, empty_vector_if_na)
  result$pval <- lapply(result$pval, empty_vector_if_na)
  result <- result[, names(result) != "n"]
  attr(result, "conf_level") <- x$control$conf_level

  result
}

#' Fits for Cox Proportional Hazards Regression
#'
#' @description `r lifecycle::badge("stable")`
#'
#' Fitting functions for univariate and multivariate Cox regression models.
#'
#' @param variables (`list`)\cr a named list corresponds to the names of variables found in `data`, passed as a named
#'   list and corresponding to `time`, `event`, `arm`, `strata`, and `covariates` terms. If `arm` is missing from
#'   `variables`, then only Cox model(s) including the `covariates` will be fitted and the corresponding effect
#'   estimates will be tabulated later.
#' @param data (`data.frame`)\cr the dataset containing the variables to fit the models.
#' @param at (`list` of `numeric`)\cr when the candidate covariate is a `numeric`, use `at` to specify
#'   the value of the covariate at which the effect should be estimated.
#' @param control (`list`)\cr a list of parameters as returned by the helper function [control_coxreg()].
#'
#' @seealso [h_cox_regression] for relevant helper functions, [cox_regression].
#'
#' @examples
#' library(survival)
#'
#' set.seed(1, kind = "Mersenne-Twister")
#'
#' # Testing dataset [survival::bladder].
#' dta_bladder <- with(
#'   data = bladder[bladder$enum < 5, ],
#'   data.frame(
#'     time = stop,
#'     status = event,
#'     armcd = as.factor(rx),
#'     covar1 = as.factor(enum),
#'     covar2 = factor(
#'       sample(as.factor(enum)),
#'       levels = 1:4, labels = c("F", "F", "M", "M")
#'     )
#'   )
#' )
#' labels <- c("armcd" = "ARM", "covar1" = "A Covariate Label", "covar2" = "Sex (F/M)")
#' formatters::var_labels(dta_bladder)[names(labels)] <- labels
#' dta_bladder$age <- sample(20:60, size = nrow(dta_bladder), replace = TRUE)
#'
#' plot(
#'   survfit(Surv(time, status) ~ armcd + covar1, data = dta_bladder),
#'   lty = 2:4,
#'   xlab = "Months",
#'   col = c("blue1", "blue2", "blue3", "blue4", "red1", "red2", "red3", "red4")
#' )
#'
#' @name fit_coxreg
NULL

#' @describeIn fit_coxreg Fit a series of univariate Cox regression models given the inputs.
#'
#' @return
#' * `fit_coxreg_univar()` returns a `coxreg.univar` class object which is a named `list`
#'   with 5 elements:
#'   * `mod`: Cox regression models fitted by [survival::coxph()].
#'   * `data`: The original data frame input.
#'   * `control`: The original control input.
#'   * `vars`: The variables used in the model.
#'   * `at`: Value of the covariate at which the effect should be estimated.
#'
#' @note When using `fit_coxreg_univar` there should be two study arms.
#'
#' @examples
#' # fit_coxreg_univar
#'
#' ## Cox regression: arm + 1 covariate.
#' mod1 <- fit_coxreg_univar(
#'   variables = list(
#'     time = "time", event = "status", arm = "armcd",
#'     covariates = "covar1"
#'   ),
#'   data = dta_bladder,
#'   control = control_coxreg(conf_level = 0.91)
#' )
#'
#' ## Cox regression: arm + 1 covariate + interaction, 2 candidate covariates.
#' mod2 <- fit_coxreg_univar(
#'   variables = list(
#'     time = "time", event = "status", arm = "armcd",
#'     covariates = c("covar1", "covar2")
#'   ),
#'   data = dta_bladder,
#'   control = control_coxreg(conf_level = 0.91, interaction = TRUE)
#' )
#'
#' ## Cox regression: arm + 1 covariate, stratified analysis.
#' mod3 <- fit_coxreg_univar(
#'   variables = list(
#'     time = "time", event = "status", arm = "armcd", strata = "covar2",
#'     covariates = c("covar1")
#'   ),
#'   data = dta_bladder,
#'   control = control_coxreg(conf_level = 0.91)
#' )
#'
#' ## Cox regression: no arm, only covariates.
#' mod4 <- fit_coxreg_univar(
#'   variables = list(
#'     time = "time", event = "status",
#'     covariates = c("covar1", "covar2")
#'   ),
#'   data = dta_bladder
#' )
#'
#' @export
fit_coxreg_univar <- function(variables,
                              data,
                              at = list(),
                              control = control_coxreg()) {
  checkmate::assert_list(variables, names = "named")
  has_arm <- "arm" %in% names(variables)
  arm_name <- if (has_arm) "arm" else NULL

  checkmate::assert_character(variables$covariates, null.ok = TRUE)

  assert_df_with_variables(data, variables)
  assert_list_of_variables(variables[c(arm_name, "event", "time")])

  if (!is.null(variables$strata)) {
    checkmate::assert_disjunct(control$pval_method, "likelihood")
  }
  if (has_arm) {
    assert_df_with_factors(data, list(val = variables$arm), min.levels = 2, max.levels = 2)
  }
  vars <- unlist(variables[c(arm_name, "covariates", "strata")], use.names = FALSE)
  for (i in vars) {
    if (is.factor(data[[i]])) {
      attr(data[[i]], "levels") <- levels(droplevels(data[[i]]))
    }
  }
  forms <- h_coxreg_univar_formulas(variables, interaction = control$interaction)
  mod <- lapply(
    forms, function(x) {
      survival::coxph(formula = stats::as.formula(x), data = data, ties = control$ties)
    }
  )
  structure(
    list(
      mod = mod,
      data = data,
      control = control,
      vars = variables,
      at = at
    ),
    class = "coxreg.univar"
  )
}

#' @describeIn fit_coxreg Fit a multivariate Cox regression model.
#'
#' @return
#' * `fit_coxreg_multivar()` returns a `coxreg.multivar` class object which is a named list
#'   with 4 elements:
#'   * `mod`: Cox regression model fitted by [survival::coxph()].
#'   * `data`: The original data frame input.
#'   * `control`: The original control input.
#'   * `vars`: The variables used in the model.
#'
#' @examples
#' # fit_coxreg_multivar
#'
#' ## Cox regression: multivariate Cox regression.
#' multivar_model <- fit_coxreg_multivar(
#'   variables = list(
#'     time = "time", event = "status", arm = "armcd",
#'     covariates = c("covar1", "covar2")
#'   ),
#'   data = dta_bladder
#' )
#'
#' # Example without treatment arm.
#' multivar_covs_model <- fit_coxreg_multivar(
#'   variables = list(
#'     time = "time", event = "status",
#'     covariates = c("covar1", "covar2")
#'   ),
#'   data = dta_bladder
#' )
#'
#' @export
fit_coxreg_multivar <- function(variables,
                                data,
                                control = control_coxreg()) {
  checkmate::assert_list(variables, names = "named")
  has_arm <- "arm" %in% names(variables)
  arm_name <- if (has_arm) "arm" else NULL

  if (!is.null(variables$covariates)) {
    checkmate::assert_character(variables$covariates)
  }

  checkmate::assert_false(control$interaction)
  assert_df_with_variables(data, variables)
  assert_list_of_variables(variables[c(arm_name, "event", "time")])

  if (!is.null(variables$strata)) {
    checkmate::assert_disjunct(control$pval_method, "likelihood")
  }

  form <- h_coxreg_multivar_formula(variables)
  mod <- survival::coxph(
    formula = stats::as.formula(form),
    data = data,
    ties = control$ties
  )
  structure(
    list(
      mod = mod,
      data = data,
      control = control,
      vars = variables
    ),
    class = "coxreg.multivar"
  )
}

#' Muffled `car::Anova`
#'
#' Applied on survival models, [car::Anova()] signal that the `strata` terms is dropped from the model formula when
#' present, this function deliberately muffles this message.
#'
#' @param mod (`coxph`)\cr Cox regression model fitted by [survival::coxph()].
#' @param test_statistic (`string`)\cr the method used for estimation of p.values; `wald` (default) or `likelihood`.
#'
#' @return Returns the output of [car::Anova()], with convergence message muffled.
#'
#' @keywords internal
muffled_car_anova <- function(mod, test_statistic) {
  tryCatch(
    withCallingHandlers(
      expr = {
        car::Anova(
          mod,
          test.statistic = test_statistic,
          type = "III"
        )
      },
      message = function(m) invokeRestart("muffleMessage"),
      error = function(e) {
        stop(paste(
          "the model seems to have convergence problems, please try to change",
          "the configuration of covariates or strata variables, e.g.",
          "- original error:", e
        ))
      }
    )
  )
}

#' Cox Regression Helper: Interactions
#'
#' @description `r lifecycle::badge("stable")`
#'
#' Test and estimate the effect of a treatment in interaction with a covariate.
#' The effect is estimated as the HR of the tested treatment for a given level
#' of the covariate, in comparison to the treatment control.
#'
#' @inheritParams argument_convention
#' @param x (`numeric` or `factor`)\cr the values of the effect to be tested.
#' @param effect (`string`)\cr the name of the effect to be tested and estimated.
#' @param covar (`string`)\cr the name of the covariate in the model.
#' @param mod (`coxph`)\cr the Cox regression model.
#' @param label (`string`)\cr the label to be returned as `term_label`.
#' @param control (`list`)\cr a list of controls as returned by [control_coxreg()].
#' @param ... see methods.
#'
#' @examples
#' library(survival)
#'
#' set.seed(1, kind = "Mersenne-Twister")
#'
#' # Testing dataset [survival::bladder].
#' dta_bladder <- with(
#'   data = bladder[bladder$enum < 5, ],
#'   data.frame(
#'     time = stop,
#'     status = event,
#'     armcd = as.factor(rx),
#'     covar1 = as.factor(enum),
#'     covar2 = factor(
#'       sample(as.factor(enum)),
#'       levels = 1:4,
#'       labels = c("F", "F", "M", "M")
#'     )
#'   )
#' )
#' labels <- c("armcd" = "ARM", "covar1" = "A Covariate Label", "covar2" = "Sex (F/M)")
#' formatters::var_labels(dta_bladder)[names(labels)] <- labels
#' dta_bladder$age <- sample(20:60, size = nrow(dta_bladder), replace = TRUE)
#'
#' plot(
#'   survfit(Surv(time, status) ~ armcd + covar1, data = dta_bladder),
#'   lty = 2:4,
#'   xlab = "Months",
#'   col = c("blue1", "blue2", "blue3", "blue4", "red1", "red2", "red3", "red4")
#' )
#'
#' @name cox_regression_inter
NULL

#' @describeIn cox_regression_inter S3 generic helper function to determine interaction effect.
#'
#' @return
#' * `h_coxreg_inter_effect()` returns a `data.frame` of covariate interaction effects consisting of the following
#'   variables: `effect`, `term`, `term_label`, `level`, `n`, `hr`, `lcl`, `ucl`, `pval`, and `pval_inter`.
#'
#' @export
h_coxreg_inter_effect <- function(x,
                                  effect,
                                  covar,
                                  mod,
                                  label,
                                  control,
                                  ...) {
  UseMethod("h_coxreg_inter_effect", x)
}

#' @describeIn cox_regression_inter Estimate the interaction with a `numeric` covariate.
#'
#' @param at (`list`)\cr a list with items named after the covariate, every
#'   item is a vector of levels at which the interaction should be estimated.
#'
#' @export
h_coxreg_inter_effect.numeric <- function(x,
                                          effect,
                                          covar,
                                          mod,
                                          label,
                                          control,
                                          at,
                                          ...) {
  betas <- stats::coef(mod)
  attrs <- attr(stats::terms(mod), "term.labels")
  term_indices <- grep(
    pattern = effect,
    x = attrs[!grepl("strata\\(", attrs)]
  )
  checkmate::assert_vector(term_indices, len = 2)
  betas <- betas[term_indices]
  betas_var <- diag(stats::vcov(mod))[term_indices]
  betas_cov <- stats::vcov(mod)[term_indices[1], term_indices[2]]
  xval <- if (is.null(at[[covar]])) {
    stats::median(x)
  } else {
    at[[covar]]
  }
  effect_index <- !grepl(covar, names(betas))
  coef_hat <- betas[effect_index] + xval * betas[!effect_index]
  coef_se <- sqrt(
    betas_var[effect_index] +
      xval ^ 2 * betas_var[!effect_index] + # styler: off
      2 * xval * betas_cov
  )
  q_norm <- stats::qnorm((1 + control$conf_level) / 2)
  data.frame(
    effect = "Covariate:",
    term = rep(covar, length(xval)),
    term_label = paste0("  ", xval),
    level = as.character(xval),
    n = NA,
    hr = exp(coef_hat),
    lcl = exp(coef_hat - q_norm * coef_se),
    ucl = exp(coef_hat + q_norm * coef_se),
    pval = NA,
    pval_inter = NA,
    stringsAsFactors = FALSE
  )
}

#' @describeIn cox_regression_inter Estimate the interaction with a `factor` covariate.
#'
#' @param data (`data.frame`)\cr the data frame on which the model was fit.
#'
#' @export
h_coxreg_inter_effect.factor <- function(x,
                                         effect,
                                         covar,
                                         mod,
                                         label,
                                         control,
                                         data,
                                         ...) {
  y <- h_coxreg_inter_estimations(
    variable = effect, given = covar,
    lvl_var = levels(data[[effect]]),
    lvl_given = levels(data[[covar]]),
    mod = mod,
    conf_level = 0.95
  )[[1]]

  data.frame(
    effect = "Covariate:",
    term = rep(covar, nrow(y)),
    term_label = as.character(paste0("  ", levels(data[[covar]]))),
    level = as.character(levels(data[[covar]])),
    n = NA,
    hr = y[, "hr"],
    lcl = y[, "lcl"],
    ucl = y[, "ucl"],
    pval = NA,
    pval_inter = NA,
    stringsAsFactors = FALSE
  )
}

#' @describeIn cox_regression_inter A higher level function to get
#'   the results of the interaction test and the estimated values.
#'
#' @return
#' * `h_coxreg_extract_interaction()` returns the result of an interaction test and the estimated values. If
#'   no interaction, [h_coxreg_univar_extract()] is applied instead.
#'
#' @examples
#' mod <- coxph(Surv(time, status) ~ armcd * covar1, data = dta_bladder)
#' h_coxreg_extract_interaction(
#'   mod = mod, effect = "armcd", covar = "covar1", data = dta_bladder,
#'   control = control_coxreg()
#' )
#'
#' @export
h_coxreg_extract_interaction <- function(effect,
                                         covar,
                                         mod,
                                         data,
                                         at,
                                         control) {
  if (!any(attr(stats::terms(mod), "order") == 2)) {
    y <- h_coxreg_univar_extract(
      effect = effect, covar = covar, mod = mod, data = data, control = control
    )
    y$pval_inter <- NA
    y
  } else {
    test_statistic <- c(wald = "Wald", likelihood = "LR")[control$pval_method]

    # Test the main treatment effect.
    mod_aov <- muffled_car_anova(mod, test_statistic)
    sum_anova <- broom::tidy(mod_aov)
    pval <- sum_anova[sum_anova$term == effect, ][["p.value"]]

    # Test the interaction effect.
    pval_inter <- sum_anova[grep(":", sum_anova$term), ][["p.value"]]
    covar_test <- data.frame(
      effect = "Covariate:",
      term = covar,
      term_label = unname(labels_or_names(data[covar])),
      level = "",
      n = mod$n, hr = NA, lcl = NA, ucl = NA, pval = pval,
      pval_inter = pval_inter,
      stringsAsFactors = FALSE
    )
    # Estimate the interaction.
    y <- h_coxreg_inter_effect(
      data[[covar]],
      covar = covar,
      effect = effect,
      mod = mod,
      label = unname(labels_or_names(data[covar])),
      at = at,
      control = control,
      data = data
    )
    rbind(covar_test, y)
  }
}

#' @describeIn cox_regression_inter Hazard ratio estimation in interactions.
#'
#' @param variable,given (`string`)\cr the name of variables in interaction. We seek the estimation
#'   of the levels of `variable` given the levels of `given`.
#' @param lvl_var,lvl_given (`character`)\cr corresponding levels has given by [levels()].
#' @param mod (`coxph`)\cr a fitted Cox regression model (see [survival::coxph()]).
#'
#' @details Given the cox regression investigating the effect of Arm (A, B, C; reference A)
#'   and Sex (F, M; reference Female) and the model being abbreviated: y ~ Arm + Sex + Arm:Sex.
#'   The cox regression estimates the coefficients along with a variance-covariance matrix for:
#'
#'   - b1 (arm b), b2 (arm c)
#'   - b3 (sex m)
#'   - b4 (arm b: sex m), b5 (arm c: sex m)
#'
#'   The estimation of the Hazard Ratio for arm C/sex M is given in reference
#'   to arm A/Sex M by exp(b2 + b3 + b5)/ exp(b3) = exp(b2 + b5).
#'   The interaction coefficient is deduced by b2 + b5 while the standard error
#'   is obtained as $sqrt(Var b2 + Var b5 + 2 * covariance (b2,b5))$.
#'
#' @return
#' * `h_coxreg_inter_estimations()` returns a list of matrices (one per level of variable) with rows corresponding
#'   to the combinations of `variable` and `given`, with columns:
#'   * `coef_hat`: Estimation of the coefficient.
#'   * `coef_se`: Standard error of the estimation.
#'   * `hr`: Hazard ratio.
#'   * `lcl, ucl`: Lower/upper confidence limit of the hazard ratio.
#'
#' @examples
#' mod <- coxph(Surv(time, status) ~ armcd * covar1, data = dta_bladder)
#' result <- h_coxreg_inter_estimations(
#'   variable = "armcd", given = "covar1",
#'   lvl_var = levels(dta_bladder$armcd),
#'   lvl_given = levels(dta_bladder$covar1),
#'   mod = mod, conf_level = .95
#' )
#' result
#'
#' @export
h_coxreg_inter_estimations <- function(variable,
                                       given,
                                       lvl_var,
                                       lvl_given,
                                       mod,
                                       conf_level = 0.95) {
  var_lvl <- paste0(variable, lvl_var[-1]) # [-1]: reference level
  giv_lvl <- paste0(given, lvl_given)
  design_mat <- expand.grid(variable = var_lvl, given = giv_lvl)
  design_mat <- design_mat[order(design_mat$variable, design_mat$given), ]
  design_mat <- within(
    data = design_mat,
    expr = {
      inter <- paste0(variable, ":", given)
      rev_inter <- paste0(given, ":", variable)
    }
  )
  split_by_variable <- design_mat$variable
  interaction_names <- paste(design_mat$variable, design_mat$given, sep = "/")

  mmat <- stats::model.matrix(mod)[1, ]
  mmat[!mmat == 0] <- 0

  design_mat <- apply(
    X = design_mat, MARGIN = 1, FUN = function(x) {
      mmat[names(mmat) %in% x[-which(names(x) == "given")]] <- 1
      mmat
    }
  )
  colnames(design_mat) <- interaction_names

  coef <- stats::coef(mod)
  vcov <- stats::vcov(mod)
  betas <- as.matrix(coef)
  coef_hat <- t(design_mat) %*% betas
  dimnames(coef_hat)[2] <- "coef"
  coef_se <- apply(
    design_mat, 2,
    function(x) {
      vcov_el <- as.logical(x)
      y <- vcov[vcov_el, vcov_el]
      y <- sum(y)
      y <- sqrt(y)
      return(y)
    }
  )
  q_norm <- stats::qnorm((1 + conf_level) / 2)
  y <- cbind(coef_hat, `se(coef)` = coef_se)
  y <- apply(y, 1, function(x) {
    x["hr"] <- exp(x["coef"])
    x["lcl"] <- exp(x["coef"] - q_norm * x["se(coef)"])
    x["ucl"] <- exp(x["coef"] + q_norm * x["se(coef)"])
    x
  })
  y <- t(y)
  y <- by(y, split_by_variable, identity)
  y <- lapply(y, as.matrix)
  attr(y, "details") <- paste0(
    "Estimations of ", variable,
    " hazard ratio given the level of ", given, " compared to ",
    variable, " level ", lvl_var[1], "."
  )
  y
}

#' Encode Categorical Missing Values in a Data Frame
#'
#' @description `r lifecycle::badge("stable")`
#'
#' This is a helper function to encode missing entries across groups of categorical
#' variables in a data frame.
#'
#' @details Missing entries are those with `NA` or empty strings and will
#'   be replaced with a specified value. If factor variables include missing
#'   values, the missing value will be inserted as the last level.
#'   Similarly, in case character or logical variables should be converted to factors
#'   with the `char_as_factor` or `logical_as_factor` options, the missing values will
#'   be set as the last level.
#'
#' @param data (`data.frame`)\cr data set.
#' @param omit_columns (`character`)\cr names of variables from `data` that should
#'   not be modified by this function.
#' @param char_as_factor (`flag`)\cr whether to convert character variables
#'   in `data` to factors.
#' @param logical_as_factor (`flag`)\cr whether to convert logical variables
#'   in `data` to factors.
#' @param na_level (`string`)\cr used to replace all `NA` or empty
#'   values inside non-`omit_columns` columns.
#'
#' @return A `data.frame` with the chosen modifications applied.
#'
#' @seealso [sas_na()] and [explicit_na()] for other missing data helper functions.
#'
#' @examples
#' my_data <- data.frame(
#'   u = c(TRUE, FALSE, NA, TRUE),
#'   v = factor(c("A", NA, NA, NA), levels = c("Z", "A")),
#'   w = c("A", "B", NA, "C"),
#'   x = c("D", "E", "F", NA),
#'   y = c("G", "H", "I", ""),
#'   z = c(1, 2, 3, 4),
#'   stringsAsFactors = FALSE
#' )
#'
#' # Example 1
#' # Encode missing values in all character or factor columns.
#' df_explicit_na(my_data)
#' # Also convert logical columns to factor columns.
#' df_explicit_na(my_data, logical_as_factor = TRUE)
#' # Encode missing values in a subset of columns.
#' df_explicit_na(my_data, omit_columns = c("x", "y"))
#'
#' # Example 2
#' # Here we purposefully convert all `M` values to `NA` in the `SEX` variable.
#' # After running `df_explicit_na` the `NA` values are encoded as `<Missing>` but they are not
#' # included when generating `rtables`.
#' adsl <- tern_ex_adsl
#' adsl$SEX[adsl$SEX == "M"] <- NA
#' adsl <- df_explicit_na(adsl)
#'
#' # If you want the `Na` values to be displayed in the table use the `na_level` argument.
#' adsl <- tern_ex_adsl
#' adsl$SEX[adsl$SEX == "M"] <- NA
#' adsl <- df_explicit_na(adsl, na_level = "Missing Values")
#'
#' # Example 3
#' # Numeric variables that have missing values are not altered. This means that any `NA` value in
#' # a numeric variable will not be included in the summary statistics, nor will they be included
#' # in the denominator value for calculating the percent values.
#' adsl <- tern_ex_adsl
#' adsl$AGE[adsl$AGE < 30] <- NA
#' adsl <- df_explicit_na(adsl)
#'
#' @export
df_explicit_na <- function(data,
                           omit_columns = NULL,
                           char_as_factor = TRUE,
                           logical_as_factor = FALSE,
                           na_level = "<Missing>") {
  checkmate::assert_character(omit_columns, null.ok = TRUE, min.len = 1, any.missing = FALSE)
  checkmate::assert_data_frame(data)
  checkmate::assert_flag(char_as_factor)
  checkmate::assert_flag(logical_as_factor)
  checkmate::assert_string(na_level)

  target_vars <- if (is.null(omit_columns)) {
    names(data)
  } else {
    setdiff(names(data), omit_columns) # May have duplicates.
  }
  if (length(target_vars) == 0) {
    return(data)
  }

  l_target_vars <- split(target_vars, target_vars)

  # Makes sure target_vars exist in data and names are not duplicated.
  assert_df_with_variables(data, l_target_vars)

  for (x in target_vars) {
    xi <- data[[x]]
    xi_label <- obj_label(xi)

    # Determine whether to convert character or logical input.
    do_char_conversion <- is.character(xi) && char_as_factor
    do_logical_conversion <- is.logical(xi) && logical_as_factor

    # Pre-convert logical to character to deal correctly with replacing NA
    # values below.
    if (do_logical_conversion) {
      xi <- as.character(xi)
    }

    if (is.factor(xi) || is.character(xi)) {
      # Handle empty strings and NA values.
      xi <- explicit_na(sas_na(xi), label = na_level)

      # Convert to factors if requested for the original type,
      # set na_level as the last value.
      if (do_char_conversion || do_logical_conversion) {
        levels_xi <- setdiff(sort(unique(xi)), na_level)
        if (na_level %in% unique(xi)) {
          levels_xi <- c(levels_xi, na_level)
        }

        xi <- factor(xi, levels = levels_xi)
      }

      data[, x] <- formatters::with_label(xi, label = xi_label)
    }
  }
  return(data)
}

#' `rtables` Access Helper Functions
#'
#' @description `r lifecycle::badge("stable")`
#'
#' These are a couple of functions that help with accessing the data in `rtables` objects.
#' Currently these work for occurrence tables, which are defined as having a count as the first
#' element and a fraction as the second element in each cell.
#'
#' @seealso [prune_occurrences] for usage of these functions.
#'
#' @name rtables_access
NULL

#' @describeIn rtables_access Helper function to extract the first values from each content
#'   cell and from specified columns in a `TableRow`. Defaults to all columns.
#'
#' @param table_row (`TableRow`)\cr an analysis row in a occurrence table.
#' @param col_names (`character`)\cr the names of the columns to extract from.
#' @param col_indices (`integer`)\cr the indices of the columns to extract from. If `col_names` are provided,
#'   then these are inferred from the names of `table_row`. Note that this currently only works well with a single
#'   column split.
#'
#' @return
#' * `h_row_first_values()` returns a `vector` of numeric values.
#'
#' @examples
#' tbl <- basic_table() %>%
#'   split_cols_by("ARM") %>%
#'   split_rows_by("RACE") %>%
#'   analyze("AGE", function(x) {
#'     list(
#'       "mean (sd)" = rcell(c(mean(x), sd(x)), format = "xx.x (xx.x)"),
#'       "n" = length(x),
#'       "frac" = rcell(c(0.1, 0.1), format = "xx (xx)")
#'     )
#'   }) %>%
#'   build_table(tern_ex_adsl) %>%
#'   prune_table()
#' tree_row_elem <- collect_leaves(tbl[2, ])[[1]]
#' result <- max(h_row_first_values(tree_row_elem))
#' result
#'
#' @export
h_row_first_values <- function(table_row,
                               col_names = NULL,
                               col_indices = NULL) {
  col_indices <- check_names_indices(table_row, col_names, col_indices)
  checkmate::assert_integerish(col_indices)
  checkmate::assert_subset(col_indices, seq_len(ncol(table_row)))

  # Main values are extracted
  row_vals <- row_values(table_row)[col_indices]

  # Main return
  vapply(row_vals, function(rv) {
    if (is.null(rv)) {
      NA_real_
    } else {
      rv[1L]
    }
  }, FUN.VALUE = numeric(1))
}

#' @describeIn rtables_access Helper function that extracts row values and checks if they are
#'   convertible to integers (`integerish` values).
#'
#' @return
#' * `h_row_counts()` returns a `vector` of numeric values.
#'
#' @examples
#' # Row counts (integer values)
#' \dontrun{
#' h_row_counts(tree_row_elem) # Fails because there are no integers
#' }
#' # Using values with integers
#' tree_row_elem <- collect_leaves(tbl[3, ])[[1]]
#' result <- h_row_counts(tree_row_elem)
#' # result
#'
#' @export
h_row_counts <- function(table_row,
                         col_names = NULL,
                         col_indices = NULL) {
  counts <- h_row_first_values(table_row, col_names, col_indices)
  checkmate::assert_integerish(counts)
  counts
}

#' @describeIn rtables_access helper function to extract fractions from specified columns in a `TableRow`.
#'   More specifically it extracts the second values from each content cell and checks it is a fraction.
#'
#' @return
#' * `h_row_fractions()` returns a `vector` of proportions.
#'
#' @examples
#' # Row fractions
#' tree_row_elem <- collect_leaves(tbl[4, ])[[1]]
#' h_row_fractions(tree_row_elem)
#'
#' @export
h_row_fractions <- function(table_row,
                            col_names = NULL,
                            col_indices = NULL) {
  col_indices <- check_names_indices(table_row, col_names, col_indices)
  row_vals <- row_values(table_row)[col_indices]
  fractions <- sapply(row_vals, "[", 2L)
  checkmate::assert_numeric(fractions, lower = 0, upper = 1)
  fractions
}

#' @describeIn rtables_access Helper function to extract column counts from specified columns in a table.
#'
#' @param table (`VTableNodeInfo`)\cr an occurrence table or row.
#'
#' @return
#' * `h_col_counts()` returns a `vector` of column counts.
#'
#' @export
h_col_counts <- function(table,
                         col_names = NULL,
                         col_indices = NULL) {
  col_indices <- check_names_indices(table, col_names, col_indices)
  counts <- col_counts(table)[col_indices]
  stats::setNames(counts, col_names)
}

#' @describeIn rtables_access Helper function to get first row of content table of current table.
#'
#' @return
#' * `h_content_first_row()` returns a row from an `rtables` table.
#'
#' @export
h_content_first_row <- function(table) {
  ct <- content_table(table)
  tree_children(ct)[[1]]
}

#' @describeIn rtables_access Helper function which says whether current table is a leaf in the tree.
#'
#' @return
#' * `is_leaf_table()` returns a `logical` value indicating whether current table is a leaf.
#'
#' @keywords internal
is_leaf_table <- function(table) {
  children <- tree_children(table)
  child_classes <- unique(sapply(children, class))
  identical(child_classes, "ElementaryTable")
}

#' @describeIn rtables_access Internal helper function that tests standard inputs for column indices.
#'
#' @return
#' * `check_names_indices` returns column indices.
#'
#' @keywords internal
check_names_indices <- function(table_row,
                                col_names = NULL,
                                col_indices = NULL) {
  if (!is.null(col_names)) {
    if (!is.null(col_indices)) {
      stop(
        "Inserted both col_names and col_indices when selecting row values. ",
        "Please choose one."
      )
    }
    col_indices <- h_col_indices(table_row, col_names)
  }
  if (is.null(col_indices)) {
    ll <- ifelse(is.null(ncol(table_row)), length(table_row), ncol(table_row))
    col_indices <- seq_len(ll)
  }

  return(col_indices)
}

#' Helper Functions for Multivariate Logistic Regression
#'
#' @description `r lifecycle::badge("stable")`
#'
#' Helper functions used in calculations for logistic regression.
#'
#' @inheritParams argument_convention
#' @param fit_glm (`glm`)\cr logistic regression model fitted by [stats::glm()] with "binomial" family.
#'   Limited functionality is also available for conditional logistic regression models fitted by
#'   [survival::clogit()], currently this is used only by [extract_rsp_biomarkers()].
#' @param x (`string` or `character`)\cr a variable or interaction term in `fit_glm` (depending on the
#'   helper function).
#'
#' @examples
#' library(dplyr)
#' library(broom)
#'
#' adrs_f <- tern_ex_adrs %>%
#'   filter(PARAMCD == "BESRSPI") %>%
#'   filter(RACE %in% c("ASIAN", "WHITE", "BLACK OR AFRICAN AMERICAN")) %>%
#'   mutate(
#'     Response = case_when(AVALC %in% c("PR", "CR") ~ 1, TRUE ~ 0),
#'     RACE = factor(RACE),
#'     SEX = factor(SEX)
#'   )
#' formatters::var_labels(adrs_f) <- c(formatters::var_labels(tern_ex_adrs), Response = "Response")
#' mod1 <- fit_logistic(
#'   data = adrs_f,
#'   variables = list(
#'     response = "Response",
#'     arm = "ARMCD",
#'     covariates = c("AGE", "RACE")
#'   )
#' )
#' mod2 <- fit_logistic(
#'   data = adrs_f,
#'   variables = list(
#'     response = "Response",
#'     arm = "ARMCD",
#'     covariates = c("AGE", "RACE"),
#'     interaction = "AGE"
#'   )
#' )
#'
#' @name h_logistic_regression
NULL

#' @describeIn h_logistic_regression Helper function to extract interaction variable names from a fitted
#'   model assuming only one interaction term.
#'
#' @return Vector of names of interaction variables.
#'
#' @export
h_get_interaction_vars <- function(fit_glm) {
  checkmate::assert_class(fit_glm, "glm")
  terms_name <- attr(stats::terms(fit_glm), "term.labels")
  terms_order <- attr(stats::terms(fit_glm), "order")
  interaction_term <- terms_name[terms_order == 2]
  checkmate::assert_string(interaction_term)
  strsplit(interaction_term, split = ":")[[1]]
}

#' @describeIn h_logistic_regression Helper function to get the right coefficient name from the
#'   interaction variable names and the given levels. The main value here is that the order
#'   of first and second variable is checked in the `interaction_vars` input.
#'
#' @param interaction_vars (`character` of length 2)\cr interaction variable names.
#' @param first_var_with_level (`character` of length 2)\cr the first variable name with
#'   the interaction level.
#' @param second_var_with_level (`character` of length 2)\cr the second variable name with
#'   the interaction level.
#'
#' @return Name of coefficient.
#'
#' @export
h_interaction_coef_name <- function(interaction_vars,
                                    first_var_with_level,
                                    second_var_with_level) {
  checkmate::assert_character(interaction_vars, len = 2, any.missing = FALSE)
  checkmate::assert_character(first_var_with_level, len = 2, any.missing = FALSE)
  checkmate::assert_character(second_var_with_level, len = 2, any.missing = FALSE)
  checkmate::assert_subset(c(first_var_with_level[1], second_var_with_level[1]), interaction_vars)

  first_name <- paste(first_var_with_level, collapse = "")
  second_name <- paste(second_var_with_level, collapse = "")
  if (first_var_with_level[1] == interaction_vars[1]) {
    paste(first_name, second_name, sep = ":")
  } else if (second_var_with_level[1] == interaction_vars[1]) {
    paste(second_name, first_name, sep = ":")
  }
}

#' @describeIn h_logistic_regression Helper function to calculate the odds ratio estimates
#'   for the case when both the odds ratio and the interaction variable are categorical.
#'
#' @param odds_ratio_var (`string`)\cr the odds ratio variable.
#' @param interaction_var (`string`)\cr the interaction variable.
#'
#' @return Odds ratio.
#'
#' @export
h_or_cat_interaction <- function(odds_ratio_var,
                                 interaction_var,
                                 fit_glm,
                                 conf_level = 0.95) {
  interaction_vars <- h_get_interaction_vars(fit_glm)
  checkmate::assert_string(odds_ratio_var)
  checkmate::assert_string(interaction_var)
  checkmate::assert_subset(c(odds_ratio_var, interaction_var), interaction_vars)
  checkmate::assert_vector(interaction_vars, len = 2)

  xs_level <- fit_glm$xlevels
  xs_coef <- stats::coef(fit_glm)
  xs_vcov <- stats::vcov(fit_glm)
  y <- list()
  for (var_level in xs_level[[odds_ratio_var]][-1]) {
    x <- list()
    for (ref_level in xs_level[[interaction_var]]) {
      coef_names <- paste0(odds_ratio_var, var_level)
      if (ref_level != xs_level[[interaction_var]][1]) {
        interaction_coef_name <- h_interaction_coef_name(
          interaction_vars,
          c(odds_ratio_var, var_level),
          c(interaction_var, ref_level)
        )
        coef_names <- c(
          coef_names,
          interaction_coef_name
        )
      }
      if (length(coef_names) > 1) {
        ones <- t(c(1, 1))
        est <- as.numeric(ones %*% xs_coef[coef_names])
        se <- sqrt(as.numeric(ones %*% xs_vcov[coef_names, coef_names] %*% t(ones)))
      } else {
        est <- xs_coef[coef_names]
        se <- sqrt(as.numeric(xs_vcov[coef_names, coef_names]))
      }
      or <- exp(est)
      ci <- exp(est + c(lcl = -1, ucl = 1) * stats::qnorm((1 + conf_level) / 2) * se)
      x[[ref_level]] <- list(or = or, ci = ci)
    }
    y[[var_level]] <- x
  }
  y
}

#' @describeIn h_logistic_regression Helper function to calculate the odds ratio estimates
#'   for the case when either the odds ratio or the interaction variable is continuous.
#'
#' @param at (`NULL` or `numeric`)\cr optional values for the interaction variable. Otherwise
#'   the median is used.
#'
#' @return Odds ratio.
#'
#' @note We don't provide a function for the case when both variables are continuous because
#'   this does not arise in this table, as the treatment arm variable will always be involved
#'   and categorical.
#'
#' @export
h_or_cont_interaction <- function(odds_ratio_var,
                                  interaction_var,
                                  fit_glm,
                                  at = NULL,
                                  conf_level = 0.95) {
  interaction_vars <- h_get_interaction_vars(fit_glm)
  checkmate::assert_string(odds_ratio_var)
  checkmate::assert_string(interaction_var)
  checkmate::assert_subset(c(odds_ratio_var, interaction_var), interaction_vars)
  checkmate::assert_vector(interaction_vars, len = 2)
  checkmate::assert_numeric(at, min.len = 1, null.ok = TRUE, any.missing = FALSE)
  xs_level <- fit_glm$xlevels
  xs_coef <- stats::coef(fit_glm)
  xs_vcov <- stats::vcov(fit_glm)
  xs_class <- attr(fit_glm$terms, "dataClasses")
  model_data <- fit_glm$model
  if (!is.null(at)) {
    checkmate::assert_set_equal(xs_class[interaction_var], "numeric")
  }
  y <- list()
  if (xs_class[interaction_var] == "numeric") {
    if (is.null(at)) {
      at <- ceiling(stats::median(model_data[[interaction_var]]))
    }

    for (var_level in xs_level[[odds_ratio_var]][-1]) {
      x <- list()
      for (increment in at) {
        coef_names <- paste0(odds_ratio_var, var_level)
        if (increment != 0) {
          interaction_coef_name <- h_interaction_coef_name(
            interaction_vars,
            c(odds_ratio_var, var_level),
            c(interaction_var, "")
          )
          coef_names <- c(
            coef_names,
            interaction_coef_name
          )
        }
        if (length(coef_names) > 1) {
          xvec <- t(c(1, increment))
          est <- as.numeric(xvec %*% xs_coef[coef_names])
          se <- sqrt(as.numeric(xvec %*% xs_vcov[coef_names, coef_names] %*% t(xvec)))
        } else {
          est <- xs_coef[coef_names]
          se <- sqrt(as.numeric(xs_vcov[coef_names, coef_names]))
        }
        or <- exp(est)
        ci <- exp(est + c(lcl = -1, ucl = 1) * stats::qnorm((1 + conf_level) / 2) * se)
        x[[as.character(increment)]] <- list(or = or, ci = ci)
      }
      y[[var_level]] <- x
    }
  } else {
    checkmate::assert_set_equal(xs_class[odds_ratio_var], "numeric")
    checkmate::assert_set_equal(xs_class[interaction_var], "factor")
    for (var_level in xs_level[[interaction_var]]) {
      coef_names <- odds_ratio_var
      if (var_level != xs_level[[interaction_var]][1]) {
        interaction_coef_name <- h_interaction_coef_name(
          interaction_vars,
          c(odds_ratio_var, ""),
          c(interaction_var, var_level)
        )
        coef_names <- c(
          coef_names,
          interaction_coef_name
        )
      }
      if (length(coef_names) > 1) {
        xvec <- t(c(1, 1))
        est <- as.numeric(xvec %*% xs_coef[coef_names])
        se <- sqrt(as.numeric(xvec %*% xs_vcov[coef_names, coef_names] %*% t(xvec)))
      } else {
        est <- xs_coef[coef_names]
        se <- sqrt(as.numeric(xs_vcov[coef_names, coef_names]))
      }
      or <- exp(est)
      ci <- exp(est + c(lcl = -1, ucl = 1) * stats::qnorm((1 + conf_level) / 2) * se)
      y[[var_level]] <- list(or = or, ci = ci)
    }
  }
  y
}

#' @describeIn h_logistic_regression Helper function to calculate the odds ratio estimates
#'   in case of an interaction. This is a wrapper for [h_or_cont_interaction()] and
#'   [h_or_cat_interaction()].
#'
#' @return Odds ratio.
#'
#' @export
h_or_interaction <- function(odds_ratio_var,
                             interaction_var,
                             fit_glm,
                             at = NULL,
                             conf_level = 0.95) {
  xs_class <- attr(fit_glm$terms, "dataClasses")
  if (any(xs_class[c(odds_ratio_var, interaction_var)] == "numeric")) {
    h_or_cont_interaction(
      odds_ratio_var,
      interaction_var,
      fit_glm,
      at = at,
      conf_level = conf_level
    )
  } else if (all(xs_class[c(odds_ratio_var, interaction_var)] == "factor")) {
    h_or_cat_interaction(
      odds_ratio_var,
      interaction_var,
      fit_glm,
      conf_level = conf_level
    )
  } else {
    stop("wrong interaction variable class, the interaction variable is not a numeric nor a factor")
  }
}

#' @describeIn h_logistic_regression Helper function to construct term labels from simple terms and the table
#'   of numbers of patients.
#'
#' @param terms (`character`)\cr simple terms.
#' @param table (`table`)\cr table containing numbers for terms.
#'
#' @return Term labels containing numbers of patients.
#'
#' @export
h_simple_term_labels <- function(terms,
                                 table) {
  checkmate::assert_true(is.table(table))
  checkmate::assert_multi_class(terms, classes = c("factor", "character"))
  terms <- as.character(terms)
  term_n <- table[terms]
  paste0(terms, ", n = ", term_n)
}

#' @describeIn h_logistic_regression Helper function to construct term labels from interaction terms and the table
#'   of numbers of patients.
#'
#' @param terms1 (`character`)\cr terms for first dimension (rows).
#' @param terms2 (`character`)\cr terms for second dimension (rows).
#' @param any (`flag`)\cr whether any of `term1` and `term2` can be fulfilled to count the
#'   number of patients. In that case they can only be scalar (strings).
#'
#' @return Term labels containing numbers of patients.
#'
#' @export
h_interaction_term_labels <- function(terms1,
                                      terms2,
                                      table,
                                      any = FALSE) {
  checkmate::assert_true(is.table(table))
  checkmate::assert_flag(any)
  checkmate::assert_multi_class(terms1, classes = c("factor", "character"))
  checkmate::assert_multi_class(terms2, classes = c("factor", "character"))
  terms1 <- as.character(terms1)
  terms2 <- as.character(terms2)
  if (any) {
    checkmate::assert_scalar(terms1)
    checkmate::assert_scalar(terms2)
    paste0(
      terms1, " or ", terms2, ", n = ",
      # Note that we double count in the initial sum the cell [terms1, terms2], therefore subtract.
      sum(c(table[terms1, ], table[, terms2])) - table[terms1, terms2]
    )
  } else {
    term_n <- table[cbind(terms1, terms2)]
    paste0(terms1, " * ", terms2, ", n = ", term_n)
  }
}

#' @describeIn h_logistic_regression Helper function to tabulate the main effect
#'   results of a (conditional) logistic regression model.
#'
#' @return Tabulated main effect results from a logistic regression model.
#'
#' @examples
#' h_glm_simple_term_extract("AGE", mod1)
#' h_glm_simple_term_extract("ARMCD", mod1)
#'
#' @export
h_glm_simple_term_extract <- function(x, fit_glm) {
  checkmate::assert_multi_class(fit_glm, c("glm", "clogit"))
  checkmate::assert_string(x)

  xs_class <- attr(fit_glm$terms, "dataClasses")
  xs_level <- fit_glm$xlevels
  xs_coef <- summary(fit_glm)$coefficients
  stats <- if (inherits(fit_glm, "glm")) {
    c("estimate" = "Estimate", "std_error" = "Std. Error", "pvalue" = "Pr(>|z|)")
  } else {
    c("estimate" = "coef", "std_error" = "se(coef)", "pvalue" = "Pr(>|z|)")
  }
  # Make sure x is not an interaction term.
  checkmate::assert_subset(x, names(xs_class))
  x_sel <- if (xs_class[x] == "numeric") x else paste0(x, xs_level[[x]][-1])
  x_stats <- as.data.frame(xs_coef[x_sel, stats, drop = FALSE], stringsAsFactors = FALSE)
  colnames(x_stats) <- names(stats)
  x_stats$estimate <- as.list(x_stats$estimate)
  x_stats$std_error <- as.list(x_stats$std_error)
  x_stats$pvalue <- as.list(x_stats$pvalue)
  x_stats$df <- as.list(1)
  if (xs_class[x] == "numeric") {
    x_stats$term <- x
    x_stats$term_label <- if (inherits(fit_glm, "glm")) {
      formatters::var_labels(fit_glm$data[x], fill = TRUE)
    } else {
      # We just fill in here with the `term` itself as we don't have the data available.
      x
    }
    x_stats$is_variable_summary <- FALSE
    x_stats$is_term_summary <- TRUE
  } else {
    checkmate::assert_class(fit_glm, "glm")
    # The reason is that we don't have the original data set in the `clogit` object
    # and therefore cannot determine the `x_numbers` here.
    x_numbers <- table(fit_glm$data[[x]])
    x_stats$term <- xs_level[[x]][-1]
    x_stats$term_label <- h_simple_term_labels(x_stats$term, x_numbers)
    x_stats$is_variable_summary <- FALSE
    x_stats$is_term_summary <- TRUE
    main_effects <- car::Anova(fit_glm, type = 3, test.statistic = "Wald")
    x_main <- data.frame(
      pvalue = main_effects[x, "Pr(>Chisq)", drop = TRUE],
      term = xs_level[[x]][1],
      term_label = paste("Reference", h_simple_term_labels(xs_level[[x]][1], x_numbers)),
      df = main_effects[x, "Df", drop = TRUE],
      stringsAsFactors = FALSE
    )
    x_main$pvalue <- as.list(x_main$pvalue)
    x_main$df <- as.list(x_main$df)
    x_main$estimate <- list(numeric(0))
    x_main$std_error <- list(numeric(0))
    if (length(xs_level[[x]][-1]) == 1) {
      x_main$pvalue <- list(numeric(0))
      x_main$df <- list(numeric(0))
    }
    x_main$is_variable_summary <- TRUE
    x_main$is_term_summary <- FALSE
    x_stats <- rbind(x_main, x_stats)
  }
  x_stats$variable <- x
  x_stats$variable_label <- if (inherits(fit_glm, "glm")) {
    formatters::var_labels(fit_glm$data[x], fill = TRUE)
  } else {
    x
  }
  x_stats$interaction <- ""
  x_stats$interaction_label <- ""
  x_stats$reference <- ""
  x_stats$reference_label <- ""
  rownames(x_stats) <- NULL
  x_stats[c(
    "variable",
    "variable_label",
    "term",
    "term_label",
    "interaction",
    "interaction_label",
    "reference",
    "reference_label",
    "estimate",
    "std_error",
    "df",
    "pvalue",
    "is_variable_summary",
    "is_term_summary"
  )]
}

#' @describeIn h_logistic_regression Helper function to tabulate the interaction term
#'   results of a logistic regression model.
#'
#' @return Tabulated interaction term results from a logistic regression model.
#'
#' @examples
#' h_glm_interaction_extract("ARMCD:AGE", mod2)
#'
#' @export
h_glm_interaction_extract <- function(x, fit_glm) {
  vars <- h_get_interaction_vars(fit_glm)
  xs_class <- attr(fit_glm$terms, "dataClasses")

  checkmate::assert_string(x)

  # Only take two-way interaction
  checkmate::assert_vector(vars, len = 2)

  # Only consider simple case: first variable in interaction is arm, a categorical variable
  checkmate::assert_disjunct(xs_class[vars[1]], "numeric")

  xs_level <- fit_glm$xlevels
  xs_coef <- summary(fit_glm)$coefficients
  main_effects <- car::Anova(fit_glm, type = 3, test.statistic = "Wald")
  stats <- c("estimate" = "Estimate", "std_error" = "Std. Error", "pvalue" = "Pr(>|z|)")
  v1_comp <- xs_level[[vars[1]]][-1]
  if (xs_class[vars[2]] == "numeric") {
    x_stats <- as.data.frame(
      xs_coef[paste0(vars[1], v1_comp, ":", vars[2]), stats, drop = FALSE],
      stringsAsFactors = FALSE
    )
    colnames(x_stats) <- names(stats)
    x_stats$term <- v1_comp
    x_numbers <- table(fit_glm$data[[vars[1]]])
    x_stats$term_label <- h_simple_term_labels(v1_comp, x_numbers)
    v1_ref <- xs_level[[vars[1]]][1]
    term_main <- v1_ref
    ref_label <- h_simple_term_labels(v1_ref, x_numbers)
  } else if (xs_class[vars[2]] != "numeric") {
    v2_comp <- xs_level[[vars[2]]][-1]
    v1_v2_grid <- expand.grid(v1 = v1_comp, v2 = v2_comp)
    x_sel <- paste(
      paste0(vars[1], v1_v2_grid$v1),
      paste0(vars[2], v1_v2_grid$v2),
      sep = ":"
    )
    x_stats <- as.data.frame(xs_coef[x_sel, stats, drop = FALSE], stringsAsFactors = FALSE)
    colnames(x_stats) <- names(stats)
    x_stats$term <- paste(v1_v2_grid$v1, "*", v1_v2_grid$v2)
    x_numbers <- table(fit_glm$data[[vars[1]]], fit_glm$data[[vars[2]]])
    x_stats$term_label <- h_interaction_term_labels(v1_v2_grid$v1, v1_v2_grid$v2, x_numbers)
    v1_ref <- xs_level[[vars[1]]][1]
    v2_ref <- xs_level[[vars[2]]][1]
    term_main <- paste(vars[1], vars[2], sep = " * ")
    ref_label <- h_interaction_term_labels(v1_ref, v2_ref, x_numbers, any = TRUE)
  }
  x_stats$df <- as.list(1)
  x_stats$pvalue <- as.list(x_stats$pvalue)
  x_stats$is_variable_summary <- FALSE
  x_stats$is_term_summary <- TRUE
  x_main <- data.frame(
    pvalue = main_effects[x, "Pr(>Chisq)", drop = TRUE],
    term = term_main,
    term_label = paste("Reference", ref_label),
    df = main_effects[x, "Df", drop = TRUE],
    stringsAsFactors = FALSE
  )
  x_main$pvalue <- as.list(x_main$pvalue)
  x_main$df <- as.list(x_main$df)
  x_main$estimate <- list(numeric(0))
  x_main$std_error <- list(numeric(0))
  x_main$is_variable_summary <- TRUE
  x_main$is_term_summary <- FALSE

  x_stats <- rbind(x_main, x_stats)
  x_stats$variable <- x
  x_stats$variable_label <- paste(
    "Interaction of",
    formatters::var_labels(fit_glm$data[vars[1]], fill = TRUE),
    "*",
    formatters::var_labels(fit_glm$data[vars[2]], fill = TRUE)
  )
  x_stats$interaction <- ""
  x_stats$interaction_label <- ""
  x_stats$reference <- ""
  x_stats$reference_label <- ""
  rownames(x_stats) <- NULL
  x_stats[c(
    "variable",
    "variable_label",
    "term",
    "term_label",
    "interaction",
    "interaction_label",
    "reference",


1		#' Missing Data
2		#'
3		#' @description `r lifecycle::badge("stable")`
4		#'
5		#' Substitute missing data with a string or factor level.
6		#'
7		#' @param x (`factor` or `character` vector)\cr values for which any missing values should be substituted.
8		#' @param label (`character`)\cr string that missing data should be replaced with.
9		#'
10		#' @return `x` with any `NA` values substituted by `label`.
11		#'
12		#' @examples
13		#' explicit_na(c(NA, "a", "b"))
14		#' is.na(explicit_na(c(NA, "a", "b")))
15		#'
16		#' explicit_na(factor(c(NA, "a", "b")))
17		#' is.na(explicit_na(factor(c(NA, "a", "b"))))
18		#'
19		#' explicit_na(sas_na(c("a", "")))
20		#'
21		#' @export
22		explicit_na <- function(x, label = "<Missing>") {
23	409x	checkmate::assert_string(label)
24
25	409x	if (is.factor(x)) {
26	307x	x <- forcats::fct_na_value_to_level(x, label)
27	307x	forcats::fct_drop(x, only = label)
28	102x	} else if (is.character(x)) {
29	102x	x[is.na(x)] <- label
30	102x	x
31		} else {
32	!	stop("only factors and character vectors allowed")
33		}
34		}
35
36		#' Convert Strings to `NA`
37		#'
38		#' @description `r lifecycle::badge("stable")`
39		#'
40		#' SAS imports missing data as empty strings or strings with whitespaces only. This helper function can be used to
41		#' convert these values to `NA`s.
42		#'
43		#' @inheritParams explicit_na
44		#' @param empty (`logical`)\cr if `TRUE` empty strings get replaced by `NA`.
45		#' @param whitespaces (`logical`)\cr if `TRUE` then strings made from whitespaces only get replaced with `NA`.
46		#'
47		#' @return `x` with `""` and/or whitespace-only values substituted by `NA`, depending on the values of
48		#' `empty` and `whitespaces`.
49		#'
50		#' @examples
51		#' sas_na(c("1", "", " ", " ", "b"))
52		#' sas_na(factor(c("", " ", "b")))
53		#'
54		#' is.na(sas_na(c("1", "", " ", " ", "b")))
55		#'
56		#' @export
57		sas_na <- function(x, empty = TRUE, whitespaces = TRUE) {
58	406x	checkmate::assert_flag(empty)
59	406x	checkmate::assert_flag(whitespaces)
60
61	406x	if (is.factor(x)) {
62	300x	empty_levels <- levels(x) == ""
63	11x	if (empty && any(empty_levels)) levels(x)[empty_levels] <- NA
64
65	300x	ws_levels <- grepl("^\\s+$", levels(x))
66	!	if (whitespaces && any(ws_levels)) levels(x)[ws_levels] <- NA
67
68	300x	x
69	106x	} else if (is.character(x)) {
70	106x	if (empty) x[x == ""] <- NA_character_
71
72	106x	if (whitespaces) x[grepl("^\\s+$", x)] <- NA_character_
73
74	106x	x
75		} else {
76	!	stop("only factors and character vectors allowed")
77		}
78		}

1		#' Compare Variables Between Groups
2		#'
3		#' @description `r lifecycle::badge("stable")`
4		#'
5		#' Comparison with a reference group for different `x` objects.
6		#'
7		#' @inheritParams argument_convention
8		#'
9		#' @note
10		#' * For factor variables, `denom` for factor proportions can only be `n` since the purpose is to compare proportions
11		#' between columns, therefore a row-based proportion would not make sense. Proportion based on `N_col` would
12		#' be difficult since we use counts for the chi-squared test statistic, therefore missing values should be accounted
13		#' for as explicit factor levels.
14		#' * If factor variables contain `NA`, these `NA` values are excluded by default. To include `NA` values
15		#' set `na.rm = FALSE` and missing values will be displayed as an `NA` level. Alternatively, an explicit
16		#' factor level can be defined for `NA` values during pre-processing via [df_explicit_na()] - the
17		#' default `na_level` (`"<Missing>"`) will also be excluded when `na.rm` is set to `TRUE`.
18		#' * For character variables, automatic conversion to factor does not guarantee that the table
19		#' will be generated correctly. In particular for sparse tables this very likely can fail.
20		#' Therefore it is always better to manually convert character variables to factors during pre-processing.
21		#' * For `compare_vars()`, the column split must define a reference group via `ref_group` so that the comparison
22		#' is well defined.
23		#'
24		#' @seealso Relevant constructor function [create_afun_compare()], and [s_summary()] which is used internally
25		#' to compute a summary within `s_compare()`.
26		#'
27		#' @name compare_variables
28		#' @include summarize_variables.R
29		NULL
30
31		#' @describeIn compare_variables S3 generic function to produce a comparison summary.
32		#'
33		#' @return
34		#' * `s_compare()` returns output of [s_summary()] and comparisons versus the reference group in the form of p-values.
35		#'
36		#' @export
37		s_compare <- function(x,
38		.ref_group,
39		.in_ref_col,
40		...) {
41	9x	UseMethod("s_compare", x)
42		}
43
44		#' @describeIn compare_variables Method for `numeric` class. This uses the standard t-test
45		#' to calculate the p-value.
46		#'
47		#' @method s_compare numeric
48		#'
49		#' @examples
50		#' # `s_compare.numeric`
51		#'
52		#' ## Usual case where both this and the reference group vector have more than 1 value.
53		#' s_compare(rnorm(10, 5, 1), .ref_group = rnorm(5, -5, 1), .in_ref_col = FALSE)
54		#'
55		#' ## If one group has not more than 1 value, then p-value is not calculated.
56		#' s_compare(rnorm(10, 5, 1), .ref_group = 1, .in_ref_col = FALSE)
57		#'
58		#' ## Empty numeric does not fail, it returns NA-filled items and no p-value.
59		#' s_compare(numeric(), .ref_group = numeric(), .in_ref_col = FALSE)
60		#'
61		#' @export
62		s_compare.numeric <- function(x,
63		.ref_group,
64		.in_ref_col,
65		...) {
66	2x	checkmate::assert_numeric(x)
67	2x	checkmate::assert_numeric(.ref_group)
68	2x	checkmate::assert_flag(.in_ref_col)
69
70	2x	y <- s_summary.numeric(x = x, ...)
71
72	2x	y$pval <- if (!.in_ref_col && n_available(x) > 1 && n_available(.ref_group) > 1) {
73	1x	stats::t.test(x, .ref_group)$p.value
74		} else {
75	1x	character()
76		}
77
78	2x	y
79		}
80
81		#' @describeIn compare_variables Method for `factor` class. This uses the chi-squared test
82		#' to calculate the p-value.
83		#'
84		#' @param denom (`string`)\cr choice of denominator for factor proportions,
85		#' can only be `n` (number of values in this row and column intersection).
86		#'
87		#' @method s_compare factor
88		#'
89		#' @examples
90		#' # `s_compare.factor`
91		#'
92		#' ## Basic usage:
93		#' x <- factor(c("a", "a", "b", "c", "a"))
94		#' y <- factor(c("a", "b", "c"))
95		#' s_compare(x = x, .ref_group = y, .in_ref_col = FALSE)
96		#'
97		#' ## Management of NA values.
98		#' x <- explicit_na(factor(c("a", "a", "b", "c", "a", NA, NA)))
99		#' y <- explicit_na(factor(c("a", "b", "c", NA)))
100		#' s_compare(x = x, .ref_group = y, .in_ref_col = FALSE, na.rm = TRUE)
101		#' s_compare(x = x, .ref_group = y, .in_ref_col = FALSE, na.rm = FALSE)
102		#'
103		#' @export
104		s_compare.factor <- function(x,
105		.ref_group,
106		.in_ref_col,
107		denom = "n",
108		na.rm = TRUE, # nolint
109		...) {
110	3x	checkmate::assert_flag(.in_ref_col)
111	3x	assert_valid_factor(x)
112	3x	assert_valid_factor(.ref_group)
113	3x	denom <- match.arg(denom)
114
115	3x	y <- s_summary.factor(
116	3x	x = x,
117	3x	denom = denom,
118	3x	na.rm = na.rm,
119		...
120		)
121
122	3x	if (na.rm) {
123	3x	x <- x[!is.na(x)] %>% fct_discard("<Missing>")
124	3x	.ref_group <- .ref_group[!is.na(.ref_group)] %>% fct_discard("<Missing>")
125		} else {
126	!	x <- x %>% explicit_na(label = "NA")
127	!	.ref_group <- .ref_group %>% explicit_na(label = "NA")
128		}
129
130	3x	checkmate::assert_factor(x, levels = levels(.ref_group), min.levels = 2)
131
132	3x	y$pval <- if (!.in_ref_col && length(x) > 0 && length(.ref_group) > 0) {
133	3x	tab <- rbind(table(x), table(.ref_group))
134	3x	res <- suppressWarnings(stats::chisq.test(tab))
135	3x	res$p.value
136		} else {
137	!	character()
138		}
139
140	3x	y
141		}
142
143		#' @describeIn compare_variables Method for `character` class. This makes an automatic
144		#' conversion to `factor` (with a warning) and then forwards to the method for factors.
145		#'
146		#' @param verbose (`logical`)\cr Whether warnings and messages should be printed. Mainly used
147		#' to print out information about factor casting. Defaults to `TRUE`.
148		#'
149		#' @method s_compare character
150		#'
151		#' @examples
152		#' # `s_compare.character`
153		#'
154		#' ## Basic usage:
155		#' x <- c("a", "a", "b", "c", "a")
156		#' y <- c("a", "b", "c")
157		#' s_compare(x, .ref_group = y, .in_ref_col = FALSE, .var = "x", verbose = FALSE)
158		#'
159		#' ## Note that missing values handling can make a large difference:
160		#' x <- c("a", "a", "b", "c", "a", NA)
161		#' y <- c("a", "b", "c", rep(NA, 20))
162		#' s_compare(x,
163		#' .ref_group = y, .in_ref_col = FALSE,
164		#' .var = "x", verbose = FALSE
165		#' )
166		#' s_compare(x,
167		#' .ref_group = y, .in_ref_col = FALSE, .var = "x",
168		#' na.rm = FALSE, verbose = FALSE
169		#' )
170		#'
171		#' @export
172		s_compare.character <- function(x,
173		.ref_group,
174		.in_ref_col,
175		denom = "n",
176		na.rm = TRUE, # nolint
177		.var,
178		verbose = TRUE,
179		...) {
180	1x	x <- as_factor_keep_attributes(x, x_name = .var, verbose = verbose)
181	1x	.ref_group <- as_factor_keep_attributes(.ref_group, x_name = .var, verbose = verbose)
182	1x	s_compare(
183	1x	x = x,
184	1x	.ref_group = .ref_group,
185	1x	.in_ref_col = .in_ref_col,
186	1x	denom = denom,
187	1x	na.rm = na.rm,
188		...
189		)
190		}
191
192		#' @describeIn compare_variables Method for `logical` class. A chi-squared test
193		#' is used. If missing values are not removed, then they are counted as `FALSE`.
194		#'
195		#' @method s_compare logical
196		#'
197		#' @examples
198		#' # `s_compare.logical`
199		#'
200		#' ## Basic usage:
201		#' x <- c(TRUE, FALSE, TRUE, TRUE)
202		#' y <- c(FALSE, FALSE, TRUE)
203		#' s_compare(x, .ref_group = y, .in_ref_col = FALSE)
204		#'
205		#' ## Management of NA values.
206		#' x <- c(NA, TRUE, FALSE)
207		#' y <- c(NA, NA, NA, NA, FALSE)
208		#' s_compare(x, .ref_group = y, .in_ref_col = FALSE, na.rm = TRUE)
209		#' s_compare(x, .ref_group = y, .in_ref_col = FALSE, na.rm = FALSE)
210		#'
211		#' @export
212		s_compare.logical <- function(x,
213		.ref_group,
214		.in_ref_col,
215		na.rm = TRUE, # nolint
216		denom = "n",
217		...) {
218	3x	denom <- match.arg(denom)
219
220	3x	y <- s_summary.logical(
221	3x	x = x,
222	3x	na.rm = na.rm,
223	3x	denom = denom,
224		...
225		)
226
227	3x	if (na.rm) {
228	2x	x <- stats::na.omit(x)
229	2x	.ref_group <- stats::na.omit(.ref_group)
230		} else {
231	1x	x[is.na(x)] <- FALSE
232	1x	.ref_group[is.na(.ref_group)] <- FALSE
233		}
234
235	3x	y$pval <- if (!.in_ref_col && length(x) > 0 && length(.ref_group) > 0) {
236	3x	x <- factor(x, levels = c(TRUE, FALSE))
237	3x	.ref_group <- factor(.ref_group, levels = c(TRUE, FALSE))
238	3x	tbl <- rbind(table(x), table(.ref_group))
239	3x	suppressWarnings(prop_chisq(tbl))
240		} else {
241	!	character()
242		}
243
244	3x	y
245		}
246
247		#' @describeIn compare_variables Formatted analysis function which is used as `afun`
248		#' in `compare_vars()`.
249		#'
250		#' @return
251		#' * `a_compare()` returns the corresponding list with formatted [rtables::CellValue()].
252		#'
253		#' @export
254		a_compare <- function(x,
255		.ref_group,
256		.in_ref_col,
257		...,
258		.var) {
259	!	UseMethod("a_compare", x)
260		}
261
262		#' @describeIn compare_variables Formatted analysis function method for `numeric` class.
263		#'
264		#' @examples
265		#' # `a_compare.numeric`
266		#' a_compare(
267		#' rnorm(10, 5, 1),
268		#' .ref_group = rnorm(20, -5, 1),
269		#' .in_ref_col = FALSE,
270		#' .var = "bla"
271		#' )
272		#'
273		#' @export
274		a_compare.numeric <- make_afun(
275		s_compare.numeric,
276		.formats = c(
277		.a_summary_numeric_formats,
278		pval = "x.xxxx \| (<0.0001)"
279		),
280		.labels = c(
281		.a_summary_numeric_labels,
282		pval = "p-value (t-test)"
283		),
284		.null_ref_cells = FALSE
285		)
286
287		.a_compare_counts_formats <- c(
288		.a_summary_counts_formats,
289		pval = "x.xxxx \| (<0.0001)"
290		)
291
292		.a_compare_counts_labels <- c(
293		pval = "p-value (chi-squared test)"
294		)
295
296		#' @describeIn compare_variables Formatted analysis function method for `factor` class.
297		#'
298		#' @examples
299		#' # `a_compare.factor`
300		#' # We need to ungroup `count` and `count_fraction` first so that the `rtables` formatting
301		#' # functions can be applied correctly.
302		#' afun <- make_afun(
303		#' getS3method("a_compare", "factor"),
304		#' .ungroup_stats = c("count", "count_fraction")
305		#' )
306		#' x <- factor(c("a", "a", "b", "c", "a"))
307		#' y <- factor(c("a", "a", "b", "c"))
308		#' afun(x, .ref_group = y, .in_ref_col = FALSE)
309		#'
310		#' @export
311		a_compare.factor <- make_afun(
312		s_compare.factor,
313		.formats = .a_compare_counts_formats,
314		.labels = .a_compare_counts_labels,
315		.null_ref_cells = FALSE
316		)
317
318		#' @describeIn compare_variables Formatted analysis function method for `character` class.
319		#'
320		#' @examples
321		#' # `a_compare.character`
322		#' afun <- make_afun(
323		#' getS3method("a_compare", "character"),
324		#' .ungroup_stats = c("count", "count_fraction")
325		#' )
326		#' x <- c("A", "B", "A", "C")
327		#' y <- c("B", "A", "C")
328		#' afun(x, .ref_group = y, .in_ref_col = FALSE, .var = "x", verbose = FALSE)
329		#'
330		#' @export
331		a_compare.character <- make_afun(
332		s_compare.character,
333		.formats = .a_compare_counts_formats,
334		.labels = .a_compare_counts_labels,
335		.null_ref_cells = FALSE
336		)
337
338		#' @describeIn compare_variables Formatted analysis function method for `logical` class.
339		#'
340		#' @examples
341		#' # `a_compare.logical`
342		#' afun <- make_afun(
343		#' getS3method("a_compare", "logical")
344		#' )
345		#' x <- c(TRUE, FALSE, FALSE, TRUE, TRUE)
346		#' y <- c(TRUE, FALSE)
347		#' afun(x, .ref_group = y, .in_ref_col = FALSE)
348		#'
349		#' @export
350		a_compare.logical <- make_afun(
351		s_compare.logical,
352		.formats = .a_compare_counts_formats,
353		.labels = .a_compare_counts_labels,
354		.null_ref_cells = FALSE
355		)
356
357		#' Constructor Function for [compare_vars()]
358		#'
359		#' @description `r lifecycle::badge("stable")`
360		#'
361		#' Constructor function which creates a combined formatted analysis function.
362		#'
363		#' @inheritParams argument_convention
364		#' @param .indent_mods (named `vector` of `integer`)\cr indent modifiers for the labels. Each element of the vector
365		#' should be a name-value pair with name corresponding to a statistic specified in `.stats` and value the indentation
366		#' for that statistic's row label.
367		#'
368		#' @return Combined formatted analysis function for use in [compare_vars()].
369		#'
370		#' @note Since [a_compare()] is generic and we want customization of the formatting arguments
371		#' via [rtables::make_afun()], we need to create another temporary generic function, with
372		#' corresponding customized methods. Then in order for the methods to be found,
373		#' we need to wrap them in a combined `afun`. Since this is required by two layout creating
374		#' functions (and possibly others in the future), we provide a constructor that does this:
375		#' [create_afun_compare()].
376		#'
377		#' @seealso [compare_vars()]
378		#'
379		#' @examples
380		#' # `create_afun_compare()` to create combined `afun`
381		#'
382		#' afun <- create_afun_compare(
383		#' .stats = c("n", "count_fraction", "mean_sd", "pval"),
384		#' .indent_mods = c(pval = 1L)
385		#' )
386		#'
387		#' lyt <- basic_table() %>%
388		#' split_cols_by("ARMCD", ref_group = "ARM A") %>%
389		#' analyze(
390		#' "AGE",
391		#' afun = afun,
392		#' show_labels = "visible"
393		#' )
394		#' build_table(lyt, df = tern_ex_adsl)
395		#'
396		#' lyt <- basic_table() %>%
397		#' split_cols_by("ARMCD", ref_group = "ARM A") %>%
398		#' analyze(
399		#' "SEX",
400		#' afun = afun,
401		#' show_labels = "visible"
402		#' )
403		#' build_table(lyt, df = tern_ex_adsl)
404		#'
405		#' @export
406		create_afun_compare <- function(.stats = NULL,
407		.formats = NULL,
408		.labels = NULL,
409		.indent_mods = NULL) {
410	3x	function(x,
411	3x	.ref_group,
412	3x	.in_ref_col,
413		...,
414	3x	.var) {
415	15x	afun <- function(x, ...) {
416	15x	UseMethod("afun", x)
417		}
418
419	15x	numeric_stats <- afun_selected_stats(
420	15x	.stats,
421	15x	all_stats = c(names(.a_summary_numeric_formats), "pval")
422		)
423	15x	afun.numeric <- make_afun( # nolint
424	15x	a_compare.numeric,
425	15x	.stats = numeric_stats,
426	15x	.formats = extract_by_name(.formats, numeric_stats),
427	15x	.labels = extract_by_name(.labels, numeric_stats),
428	15x	.indent_mods = extract_by_name(.indent_mods, numeric_stats),
429	15x	.null_ref_cells = FALSE
430		)
431
432	15x	factor_stats <- afun_selected_stats(
433	15x	.stats,
434	15x	all_stats = names(.a_compare_counts_formats)
435		)
436	15x	ungroup_stats <- afun_selected_stats(.stats, c("count", "count_fraction"))
437	15x	afun.factor <- make_afun( # nolint
438	15x	a_compare.factor,
439	15x	.stats = factor_stats,
440	15x	.formats = extract_by_name(.formats, factor_stats),
441	15x	.labels = extract_by_name(.labels, factor_stats),
442	15x	.indent_mods = extract_by_name(.indent_mods, factor_stats),
443	15x	.ungroup_stats = ungroup_stats,
444	15x	.null_ref_cells = FALSE
445		)
446
447	15x	afun.character <- make_afun( # nolint
448	15x	a_compare.character,
449	15x	.stats = factor_stats,
450	15x	.formats = extract_by_name(.formats, factor_stats),
451	15x	.labels = extract_by_name(.labels, factor_stats),
452	15x	.indent_mods = extract_by_name(.indent_mods, factor_stats),
453	15x	.ungroup_stats = ungroup_stats,
454	15x	.null_ref_cells = FALSE
455		)
456
457	15x	afun.logical <- make_afun( # nolint
458	15x	a_compare.logical,
459	15x	.stats = factor_stats,
460	15x	.formats = extract_by_name(.formats, factor_stats),
461	15x	.labels = extract_by_name(.labels, factor_stats),
462	15x	.indent_mods = extract_by_name(.indent_mods, factor_stats),
463	15x	.null_ref_cells = FALSE
464		)
465
466	15x	afun(
467	15x	x = x,
468	15x	.ref_group = .ref_group,
469	15x	.in_ref_col = .in_ref_col,
470		...,
471	15x	.var = .var
472		)
473		}
474		}
475
476		#' @describeIn compare_variables Layout-creating function which can take statistics function arguments
477		#' and additional format arguments. This function is a wrapper for [rtables::analyze()].
478		#'
479		#' @param ... arguments passed to `s_compare()`.
480		#' @param .indent_mods (named `vector` of `integer`)\cr indent modifiers for the labels. Each element of the vector
481		#' should be a name-value pair with name corresponding to a statistic specified in `.stats` and value the indentation
482		#' for that statistic's row label.
483		#'
484		#' @return
485		#' * `compare_vars()` returns a layout object suitable for passing to further layouting functions,
486		#' or to [rtables::build_table()]. Adding this function to an `rtable` layout will add formatted rows containing
487		#' the statistics from `s_compare()` to the table layout.
488		#'
489		#' @examples
490		#' # `compare_vars()` in `rtables` pipelines
491		#'
492		#' ## Default output within a `rtables` pipeline.
493		#' lyt <- basic_table() %>%
494		#' split_cols_by("ARMCD", ref_group = "ARM B") %>%
495		#' compare_vars(c("AGE", "SEX"))
496		#' build_table(lyt, tern_ex_adsl)
497		#'
498		#' ## Select and format statistics output.
499		#' lyt <- basic_table() %>%
500		#' split_cols_by("ARMCD", ref_group = "ARM C") %>%
501		#' compare_vars(
502		#' vars = "AGE",
503		#' .stats = c("mean_sd", "pval"),
504		#' .formats = c(mean_sd = "xx.x, xx.x"),
505		#' .labels = c(mean_sd = "Mean, SD")
506		#' )
507		#' build_table(lyt, df = tern_ex_adsl)
508		#'
509		#' @export
510		compare_vars <- function(lyt,
511		vars,
512		var_labels = vars,
513		nested = TRUE,
514		...,
515		na_level = NA_character_,
516		show_labels = "default",
517		table_names = vars,
518		.stats = c("n", "mean_sd", "count_fraction", "pval"),
519		.formats = NULL,
520		.labels = NULL,
521		.indent_mods = NULL) {
522	3x	afun <- create_afun_compare(.stats, .formats, .labels, .indent_mods)
523
524	3x	analyze(
525	3x	lyt = lyt,
526	3x	vars = vars,
527	3x	var_labels = var_labels,
528	3x	afun = afun,
529	3x	nested = nested,
530	3x	extra_args = list(...),
531	3x	na_str = na_level,
532	3x	inclNAs = TRUE,
533	3x	show_labels = show_labels,
534	3x	table_names = table_names
535		)
536		}

1		#' Summary numeric variables in columns
2		#'
3		#' @description `r lifecycle::badge("experimental")`
4		#'
5		#' Layout-creating function which can be used for creating column-wise summary tables.
6		#' This function sets the analysis methods as column labels and is a wrapper for
7		#' [rtables::analyze_colvars()]. It was designed principally for PK tables.
8		#'
9		#' @inheritParams argument_convention
10		#' @inheritParams rtables::analyze_colvars
11		#' @param row_labels (`character`)\cr as this function works in columns space, usual `.labels`
12		#' character vector applies on the column space. You can change the row labels by defining this
13		#' parameter to a named character vector with names corresponding to the split values. It defaults
14		#' to `NULL` and if it contains only one `string`, it will duplicate that as a row label.
15		#' @param do_summarize_row_groups (`flag`)\cr defaults to `FALSE` and applies the analysis to the current
16		#' label rows. This is a wrapper of [rtables::summarize_row_groups()] and it can accept `labelstr`
17		#' to define row labels. This behavior is not supported as we never need to overload row labels.
18		#' @param split_col_vars (`flag`)\cr defaults to `TRUE` and puts the analysis results onto the columns.
19		#' This option allows you to add multiple instances of this functions, also in a nested fashion,
20		#' without adding more splits. This split must happen only one time on a single layout.
21		#'
22		#' @return
23		#' A layout object suitable for passing to further layouting functions, or to [rtables::build_table()].
24		#' Adding this function to an `rtable` layout will summarize the given variables, arrange the output
25		#' in columns, and add it to the table layout.
26		#'
27		#' @note This is an experimental implementation of [rtables::summarize_row_groups()] and
28		#' [rtables::analyze_colvars()] that may be subjected to changes as `rtables` extends its
29		#' support to more complex analysis pipelines on the column space. For the same reasons,
30		#' we encourage to read the examples carefully and file issues for cases that differ from
31		#' them.
32		#'
33		#' Here `labelstr` behaves differently than usual. If it is not defined (default as `NULL`),
34		#' row labels are assigned automatically to the split values in case of `rtables::analyze_colvars`
35		#' (`do_summarize_row_groups = FALSE`, the default), and to the group label for
36		#' `do_summarize_row_groups = TRUE`.
37		#'
38		#' @seealso [summarize_vars()], [rtables::analyze_colvars()].
39		#'
40		#' @examples
41		#' library(dplyr)
42		#'
43		#' # Data preparation
44		#' adpp <- tern_ex_adpp %>% h_pkparam_sort()
45		#'
46		#' lyt <- basic_table() %>%
47		#' split_rows_by(var = "STRATA1", label_pos = "topleft") %>%
48		#' split_rows_by(
49		#' var = "SEX",
50		#' label_pos = "topleft",
51		#' child_label = "hidden"
52		#' ) %>% # Removes duplicated labels
53		#' analyze_vars_in_cols(vars = "AGE")
54		#' result <- build_table(lyt = lyt, df = adpp)
55		#' result
56		#'
57		#' # By selecting just some statistics and ad-hoc labels
58		#' lyt <- basic_table() %>%
59		#' split_rows_by(var = "ARM", label_pos = "topleft") %>%
60		#' split_rows_by(
61		#' var = "SEX",
62		#' label_pos = "topleft",
63		#' child_labels = "hidden",
64		#' split_fun = drop_split_levels
65		#' ) %>%
66		#' analyze_vars_in_cols(
67		#' vars = "AGE",
68		#' .stats = c("n", "cv", "geom_mean"),
69		#' .labels = c(
70		#' n = "aN",
71		#' cv = "aCV",
72		#' geom_mean = "aGeomMean"
73		#' )
74		#' )
75		#' result <- build_table(lyt = lyt, df = adpp)
76		#' result
77		#'
78		#' # Changing row labels
79		#' lyt <- basic_table() %>%
80		#' analyze_vars_in_cols(
81		#' vars = "AGE",
82		#' row_labels = "some custom label"
83		#' )
84		#' result <- build_table(lyt, df = adpp)
85		#' result
86		#'
87		#' # Pharmacokinetic parameters
88		#' lyt <- basic_table() %>%
89		#' split_rows_by(
90		#' var = "TLG_DISPLAY",
91		#' split_label = "PK Parameter",
92		#' label_pos = "topleft",
93		#' child_label = "hidden"
94		#' ) %>%
95		#' analyze_vars_in_cols(
96		#' vars = "AVAL"
97		#' )
98		#' result <- build_table(lyt, df = adpp)
99		#' result
100		#'
101		#' # Multiple calls (summarize label and analyze underneath)
102		#' lyt <- basic_table() %>%
103		#' split_rows_by(
104		#' var = "TLG_DISPLAY",
105		#' split_label = "PK Parameter",
106		#' label_pos = "topleft"
107		#' ) %>%
108		#' analyze_vars_in_cols(
109		#' vars = "AVAL",
110		#' do_summarize_row_groups = TRUE # does a summarize level
111		#' ) %>%
112		#' split_rows_by("SEX",
113		#' child_label = "hidden",
114		#' label_pos = "topleft"
115		#' ) %>%
116		#' analyze_vars_in_cols(
117		#' vars = "AVAL",
118		#' split_col_vars = FALSE # avoids re-splitting the columns
119		#' )
120		#' result <- build_table(lyt, df = adpp)
121		#' result
122		#'
123		#' @export
124		analyze_vars_in_cols <- function(lyt,
125		vars,
126		...,
127		.stats = c(
128		"n",
129		"mean",
130		"sd",
131		"se",
132		"cv",
133		"geom_cv"
134		),
135		.labels = c(
136		n = "n",
137		mean = "Mean",
138		sd = "SD",
139		se = "SE",
140		cv = "CV (%)",
141		geom_cv = "CV % Geometric Mean"
142		),
143		row_labels = NULL,
144		do_summarize_row_groups = FALSE,
145		split_col_vars = TRUE,
146		.indent_mods = NULL,
147		nested = TRUE,
148		na_level = NULL,
149		.formats = NULL) {
150	6x	checkmate::assert_string(na_level, null.ok = TRUE)
151	6x	checkmate::assert_character(row_labels, null.ok = TRUE)
152	6x	checkmate::assert_int(.indent_mods, null.ok = TRUE)
153	6x	checkmate::assert_flag(nested)
154	6x	checkmate::assert_flag(split_col_vars)
155	6x	checkmate::assert_flag(do_summarize_row_groups)
156
157		# Automatic assignment of formats
158	6x	if (is.null(.formats)) {
159		# General values
160	6x	sf_numeric <- summary_formats("numeric")
161	6x	sf_counts <- summary_formats("counts")[-1]
162	6x	formats_v <- c(sf_numeric, sf_counts)
163		} else {
164	!	formats_v <- .formats
165		}
166
167		# Check for vars in the case that one or more are used
168	6x	if (length(vars) == 1) {
169	5x	vars <- rep(vars, length(.stats))
170	1x	} else if (length(vars) != length(.stats)) {
171	1x	stop(
172	1x	"Analyzed variables (vars) does not have the same ",
173	1x	"number of elements of specified statistics (.stats)."
174		)
175		}
176
177	5x	if (split_col_vars) {
178		# Checking there is not a previous identical column split
179	4x	clyt <- tail(clayout(lyt), 1)[[1]]
180
181	4x	dummy_lyt <- split_cols_by_multivar(
182	4x	lyt = basic_table(),
183	4x	vars = vars,
184	4x	varlabels = .labels[.stats]
185		)
186
187	4x	if (any(sapply(clyt, identical, y = get_last_col_split(dummy_lyt)))) {
188	!	stop(
189	!	"Column split called again with the same values. ",
190	!	"This can create many unwanted columns. Please consider adding ",
191	!	"split_col_vars = FALSE to the last call of ",
192	!	deparse(sys.calls()[[sys.nframe() - 1]]), "."
193		)
194		}
195
196		# Main col split
197	4x	lyt <- split_cols_by_multivar(
198	4x	lyt = lyt,
199	4x	vars = vars,
200	4x	varlabels = .labels[.stats]
201		)
202		}
203
204	5x	if (do_summarize_row_groups) {
205	2x	if (length(unique(vars)) > 1) {
206	!	stop("When using do_summarize_row_groups only one label level var should be inserted.")
207		}
208
209		# Function list for do_summarize_row_groups. Slightly different handling of labels
210	2x	cfun_list <- Map(
211	2x	function(stat) {
212	12x	function(u, .spl_context, labelstr, ...) {
213		# Statistic
214	24x	res <- s_summary(u, ...)[[stat]]
215
216		# Label check and replacement
217	24x	if (length(row_labels) > 1) {
218	12x	if (!(labelstr %in% names(row_labels))) {
219	!	stop(
220	!	"Replacing the labels in do_summarize_row_groups needs a named vector",
221	!	"that contains the split values. In the current split variable ",
222	!	.spl_context$split[nrow(.spl_context)],
223	!	" the labelstr value (split value by default) ", labelstr, " is not in",
224	!	" row_labels names: ", names(row_labels)
225		)
226		}
227	12x	lbl <- unlist(row_labels[labelstr])
228		} else {
229	12x	lbl <- labelstr
230		}
231
232		# Cell creation
233	24x	rcell(res,
234	24x	label = lbl,
235	24x	format = formats_v[names(formats_v) == stat][[1]],
236	24x	format_na_str = na_level,
237	24x	indent_mod = ifelse(is.null(.indent_mods), 0L, .indent_mods)
238		)
239		}
240		},
241	2x	stat = .stats
242		)
243
244		# Main call to rtables
245	2x	summarize_row_groups(
246	2x	lyt = lyt,
247	2x	var = unique(vars),
248	2x	cfun = cfun_list,
249	2x	extra_args = list(...)
250		)
251		} else {
252		# Function list for analyze_colvars
253	3x	afun_list <- Map(
254	3x	function(stat) {
255	15x	function(u, .spl_context, ...) {
256		# Main statistics
257	78x	res <- s_summary(u, ...)[[stat]]
258
259		# Label from context
260	78x	label_from_context <- .spl_context$value[nrow(.spl_context)]
261
262		# Label switcher
263	78x	if (is.null(row_labels)) {
264	18x	lbl <- label_from_context
265		} else {
266	60x	if (length(row_labels) > 1) {
267	48x	if (!(label_from_context %in% names(row_labels))) {
268	!	stop(
269	!	"Replacing the labels in do_summarize_row_groups needs a named vector",
270	!	"that contains the split values. In the current split variable ",
271	!	.spl_context$split[nrow(.spl_context)],
272	!	" the split value ", label_from_context, " is not in",
273	!	" row_labels names: ", names(row_labels)
274		)
275		}
276	48x	lbl <- unlist(row_labels[label_from_context])
277		} else {
278	12x	lbl <- row_labels
279		}
280		}
281
282		# Cell creation
283	78x	rcell(res,
284	78x	label = lbl,
285	78x	format = formats_v[names(formats_v) == stat][[1]],
286	78x	format_na_str = na_level,
287	78x	indent_mod = ifelse(is.null(.indent_mods), 0L, .indent_mods)
288		)
289		}
290		},
291	3x	stat = .stats
292		)
293
294		# Main call to rtables
295	3x	analyze_colvars(lyt,
296	3x	afun = afun_list,
297	3x	nested = nested,
298	3x	extra_args = list(...)
299		)
300		}
301		}
302
303		# Help function
304		get_last_col_split <- function(lyt) {
305	!	tail(tail(clayout(lyt), 1)[[1]], 1)[[1]]
306		}

1		#' Controls for Cox Regression
2		#'
3		#' @description `r lifecycle::badge("stable")`
4		#'
5		#' Sets a list of parameters for Cox regression fit. Used internally.
6		#'
7		#' @inheritParams argument_convention
8		#' @param pval_method (`string`)\cr the method used for estimation of p.values; `wald` (default) or `likelihood`.
9		#' @param interaction (`flag`)\cr if `TRUE`, the model includes the interaction between the studied
10		#' treatment and candidate covariate. Note that for univariate models without treatment arm, and
11		#' multivariate models, no interaction can be used so that this needs to be `FALSE`.
12		#' @param ties (`string`)\cr among `exact` (equivalent to `DISCRETE` in SAS), `efron` and `breslow`,
13		#' see [survival::coxph()]. Note: there is no equivalent of SAS `EXACT` method in R.
14		#'
15		#' @return A `list` of items with names corresponding to the arguments.
16		#'
17		#' @seealso [fit_coxreg_univar()] and [fit_coxreg_multivar()].
18		#'
19		#' @examples
20		#' control_coxreg()
21		#'
22		#' @export
23		control_coxreg <- function(pval_method = c("wald", "likelihood"),
24		ties = c("exact", "efron", "breslow"),
25		conf_level = 0.95,
26		interaction = FALSE) {
27	40x	pval_method <- match.arg(pval_method)
28	40x	ties <- match.arg(ties)
29	40x	checkmate::assert_flag(interaction)
30	40x	assert_proportion_value(conf_level)
31	40x	list(
32	40x	pval_method = pval_method,
33	40x	ties = ties,
34	40x	conf_level = conf_level,
35	40x	interaction = interaction
36		)
37		}
38
39		#' Custom Tidy Methods for Cox Regression
40		#'
41		#' @description `r lifecycle::badge("stable")`
42		#'
43		#' @inheritParams argument_convention
44		#' @param x (`list`)\cr Result of the Cox regression model fitted by [fit_coxreg_univar()] (for univariate models)
45		#' or [fit_coxreg_multivar()] (for multivariate models).
46		#'
47		#' @return [tidy()] returns:
48		#' * For `summary.coxph` objects, a `data.frame` with columns: `Pr(>\|z\|)`, `exp(coef)`, `exp(-coef)`, `lower .95`,
49		#' `upper .95`, `level`, and `n`.
50		#' * For `coxreg.univar` objects, a `data.frame` with columns: `effect`, `term`, `term_label`, `level`, `n`, `hr`,
51		#' `lcl`, `ucl`, `pval`, and `ci`.
52		#' * For `coxreg.multivar` objects, a `data.frame` with columns: `term`, `pval`, `term_label`, `hr`, `lcl`, `ucl`,
53		#' `level`, and `ci`.
54		#'
55		#' @seealso [cox_regression]
56		#'
57		#' @name tidy_coxreg
58		NULL
59
60		#' @describeIn tidy_coxreg Custom tidy method for [survival::coxph()] summary results.
61		#'
62		#' Tidy the [survival::coxph()] results into a `data.frame` to extract model results.
63		#'
64		#' @method tidy summary.coxph
65		#'
66		#' @examples
67		#' library(survival)
68		#' library(broom)
69		#'
70		#' set.seed(1, kind = "Mersenne-Twister")
71		#'
72		#' dta_bladder <- with(
73		#' data = bladder[bladder$enum < 5, ],
74		#' data.frame(
75		#' time = stop,
76		#' status = event,
77		#' armcd = as.factor(rx),
78		#' covar1 = as.factor(enum),
79		#' covar2 = factor(
80		#' sample(as.factor(enum)),
81		#' levels = 1:4, labels = c("F", "F", "M", "M")
82		#' )
83		#' )
84		#' )
85		#' labels <- c("armcd" = "ARM", "covar1" = "A Covariate Label", "covar2" = "Sex (F/M)")
86		#' formatters::var_labels(dta_bladder)[names(labels)] <- labels
87		#' dta_bladder$age <- sample(20:60, size = nrow(dta_bladder), replace = TRUE)
88		#'
89		#' formula <- "survival::Surv(time, status) ~ armcd + covar1"
90		#' msum <- summary(coxph(stats::as.formula(formula), data = dta_bladder))
91		#' tidy(msum)
92		#'
93		#' @export
94		tidy.summary.coxph <- function(x, # nolint
95		...) {
96	120x	checkmate::assert_class(x, "summary.coxph")
97	120x	pval <- x$coefficients
98	120x	confint <- x$conf.int
99	120x	levels <- rownames(pval)
100
101	120x	pval <- tibble::as_tibble(pval)
102	120x	confint <- tibble::as_tibble(confint)
103
104	120x	ret <- cbind(pval[, grepl("Pr", names(pval))], confint)
105	120x	ret$level <- levels
106	120x	ret$n <- x[["n"]]
107	120x	ret
108		}
109
110		#' @describeIn tidy_coxreg Custom tidy method for a univariate Cox regression.
111		#'
112		#' Tidy up the result of a Cox regression model fitted by [fit_coxreg_univar()].
113		#'
114		#' @method tidy coxreg.univar
115		#'
116		#' @examples
117		#' ## Cox regression: arm + 1 covariate.
118		#' mod1 <- fit_coxreg_univar(
119		#' variables = list(
120		#' time = "time", event = "status", arm = "armcd",
121		#' covariates = "covar1"
122		#' ),
123		#' data = dta_bladder,
124		#' control = control_coxreg(conf_level = 0.91)
125		#' )
126		#'
127		#' ## Cox regression: arm + 1 covariate + interaction, 2 candidate covariates.
128		#' mod2 <- fit_coxreg_univar(
129		#' variables = list(
130		#' time = "time", event = "status", arm = "armcd",
131		#' covariates = c("covar1", "covar2")
132		#' ),
133		#' data = dta_bladder,
134		#' control = control_coxreg(conf_level = 0.91, interaction = TRUE)
135		#' )
136		#'
137		#' tidy(mod1)
138		#' tidy(mod2)
139		#'
140		#' @export
141		tidy.coxreg.univar <- function(x, # nolint
142		...) {
143	26x	checkmate::assert_class(x, "coxreg.univar")
144	26x	mod <- x$mod
145	26x	vars <- c(x$vars$arm, x$vars$covariates)
146	26x	has_arm <- "arm" %in% names(x$vars)
147
148	26x	result <- if (!has_arm) {
149	5x	Map(
150	5x	mod = mod, vars = vars,
151	5x	f = function(mod, vars) {
152	6x	h_coxreg_multivar_extract(
153	6x	var = vars,
154	6x	data = x$data,
155	6x	mod = mod,
156	6x	control = x$control
157		)
158		}
159		)
160	26x	} else if (x$control$interaction) {
161	8x	Map(
162	8x	mod = mod, covar = vars,
163	8x	f = function(mod, covar) {
164	17x	h_coxreg_extract_interaction(
165	17x	effect = x$vars$arm, covar = covar, mod = mod, data = x$data,
166	17x	at = x$at, control = x$control
167		)
168		}
169		)
170		} else {
171	13x	Map(
172	13x	mod = mod, vars = vars,
173	13x	f = function(mod, vars) {
174	34x	h_coxreg_univar_extract(
175	34x	effect = x$vars$arm, covar = vars, data = x$data, mod = mod,
176	34x	control = x$control
177		)
178		}
179		)
180		}
181	26x	result <- do.call(rbind, result)
182
183	26x	result$ci <- Map(lcl = result$lcl, ucl = result$ucl, f = function(lcl, ucl) c(lcl, ucl))
184	26x	result$n <- lapply(result$n, empty_vector_if_na)
185	26x	result$ci <- lapply(result$ci, empty_vector_if_na)
186	26x	result$hr <- lapply(result$hr, empty_vector_if_na)
187	26x	if (x$control$interaction) {
188	8x	result$pval_inter <- lapply(result$pval_inter, empty_vector_if_na)
189		# Remove interaction p-values due to change in specifications.
190	8x	result$pval[result$effect != "Treatment:"] <- NA
191		}
192	26x	result$pval <- lapply(result$pval, empty_vector_if_na)
193	26x	attr(result, "conf_level") <- x$control$conf_level
194	26x	result
195		}
196
197		#' @describeIn tidy_coxreg Custom tidy method for a multivariate Cox regression.
198		#'
199		#' Tidy up the result of a Cox regression model fitted by [fit_coxreg_multivar()].
200		#'
201		#' @method tidy coxreg.multivar
202		#'
203		#' @examples
204		#' multivar_model <- fit_coxreg_multivar(
205		#' variables = list(
206		#' time = "time", event = "status", arm = "armcd",
207		#' covariates = c("covar1", "covar2")
208		#' ),
209		#' data = dta_bladder
210		#' )
211		#' broom::tidy(multivar_model)
212		#'
213		#' @export
214		tidy.coxreg.multivar <- function(x, # nolint
215		...) {
216	8x	checkmate::assert_class(x, "coxreg.multivar")
217	8x	vars <- c(x$vars$arm, x$vars$covariates)
218
219		# Convert the model summaries to data.
220	8x	result <- Map(
221	8x	vars = vars,
222	8x	f = function(vars) {
223	28x	h_coxreg_multivar_extract(
224	28x	var = vars, data = x$data,
225	28x	mod = x$mod, control = x$control
226		)
227		}
228		)
229	8x	result <- do.call(rbind, result)
230
231	8x	result$ci <- Map(lcl = result$lcl, ucl = result$ucl, f = function(lcl, ucl) c(lcl, ucl))
232	8x	result$ci <- lapply(result$ci, empty_vector_if_na)
233	8x	result$hr <- lapply(result$hr, empty_vector_if_na)
234	8x	result$pval <- lapply(result$pval, empty_vector_if_na)
235	8x	result <- result[, names(result) != "n"]
236	8x	attr(result, "conf_level") <- x$control$conf_level
237
238	8x	result
239		}
240
241		#' Fits for Cox Proportional Hazards Regression
242		#'
243		#' @description `r lifecycle::badge("stable")`
244		#'
245		#' Fitting functions for univariate and multivariate Cox regression models.
246		#'
247		#' @param variables (`list`)\cr a named list corresponds to the names of variables found in `data`, passed as a named
248		#' list and corresponding to `time`, `event`, `arm`, `strata`, and `covariates` terms. If `arm` is missing from
249		#' `variables`, then only Cox model(s) including the `covariates` will be fitted and the corresponding effect
250		#' estimates will be tabulated later.
251		#' @param data (`data.frame`)\cr the dataset containing the variables to fit the models.
252		#' @param at (`list` of `numeric`)\cr when the candidate covariate is a `numeric`, use `at` to specify
253		#' the value of the covariate at which the effect should be estimated.
254		#' @param control (`list`)\cr a list of parameters as returned by the helper function [control_coxreg()].
255		#'
256		#' @seealso [h_cox_regression] for relevant helper functions, [cox_regression].
257		#'
258		#' @examples
259		#' library(survival)
260		#'
261		#' set.seed(1, kind = "Mersenne-Twister")
262		#'
263		#' # Testing dataset [survival::bladder].
264		#' dta_bladder <- with(
265		#' data = bladder[bladder$enum < 5, ],
266		#' data.frame(
267		#' time = stop,
268		#' status = event,
269		#' armcd = as.factor(rx),
270		#' covar1 = as.factor(enum),
271		#' covar2 = factor(
272		#' sample(as.factor(enum)),
273		#' levels = 1:4, labels = c("F", "F", "M", "M")
274		#' )
275		#' )
276		#' )
277		#' labels <- c("armcd" = "ARM", "covar1" = "A Covariate Label", "covar2" = "Sex (F/M)")
278		#' formatters::var_labels(dta_bladder)[names(labels)] <- labels
279		#' dta_bladder$age <- sample(20:60, size = nrow(dta_bladder), replace = TRUE)
280		#'
281		#' plot(
282		#' survfit(Surv(time, status) ~ armcd + covar1, data = dta_bladder),
283		#' lty = 2:4,
284		#' xlab = "Months",
285		#' col = c("blue1", "blue2", "blue3", "blue4", "red1", "red2", "red3", "red4")
286		#' )
287		#'
288		#' @name fit_coxreg
289		NULL
290
291		#' @describeIn fit_coxreg Fit a series of univariate Cox regression models given the inputs.
292		#'
293		#' @return
294		#' * `fit_coxreg_univar()` returns a `coxreg.univar` class object which is a named `list`
295		#' with 5 elements:
296		#' * `mod`: Cox regression models fitted by [survival::coxph()].
297		#' * `data`: The original data frame input.
298		#' * `control`: The original control input.
299		#' * `vars`: The variables used in the model.
300		#' * `at`: Value of the covariate at which the effect should be estimated.
301		#'
302		#' @note When using `fit_coxreg_univar` there should be two study arms.
303		#'
304		#' @examples
305		#' # fit_coxreg_univar
306		#'
307		#' ## Cox regression: arm + 1 covariate.
308		#' mod1 <- fit_coxreg_univar(
309		#' variables = list(
310		#' time = "time", event = "status", arm = "armcd",
311		#' covariates = "covar1"
312		#' ),
313		#' data = dta_bladder,
314		#' control = control_coxreg(conf_level = 0.91)
315		#' )
316		#'
317		#' ## Cox regression: arm + 1 covariate + interaction, 2 candidate covariates.
318		#' mod2 <- fit_coxreg_univar(
319		#' variables = list(
320		#' time = "time", event = "status", arm = "armcd",
321		#' covariates = c("covar1", "covar2")
322		#' ),
323		#' data = dta_bladder,
324		#' control = control_coxreg(conf_level = 0.91, interaction = TRUE)
325		#' )
326		#'
327		#' ## Cox regression: arm + 1 covariate, stratified analysis.
328		#' mod3 <- fit_coxreg_univar(
329		#' variables = list(
330		#' time = "time", event = "status", arm = "armcd", strata = "covar2",
331		#' covariates = c("covar1")
332		#' ),
333		#' data = dta_bladder,
334		#' control = control_coxreg(conf_level = 0.91)
335		#' )
336		#'
337		#' ## Cox regression: no arm, only covariates.
338		#' mod4 <- fit_coxreg_univar(
339		#' variables = list(
340		#' time = "time", event = "status",
341		#' covariates = c("covar1", "covar2")
342		#' ),
343		#' data = dta_bladder
344		#' )
345		#'
346		#' @export
347		fit_coxreg_univar <- function(variables,
348		data,
349		at = list(),
350		control = control_coxreg()) {
351	31x	checkmate::assert_list(variables, names = "named")
352	31x	has_arm <- "arm" %in% names(variables)
353	31x	arm_name <- if (has_arm) "arm" else NULL
354
355	31x	checkmate::assert_character(variables$covariates, null.ok = TRUE)
356
357	31x	assert_df_with_variables(data, variables)
358	31x	assert_list_of_variables(variables[c(arm_name, "event", "time")])
359
360	31x	if (!is.null(variables$strata)) {
361	4x	checkmate::assert_disjunct(control$pval_method, "likelihood")
362		}
363	30x	if (has_arm) {
364	24x	assert_df_with_factors(data, list(val = variables$arm), min.levels = 2, max.levels = 2)
365		}
366	29x	vars <- unlist(variables[c(arm_name, "covariates", "strata")], use.names = FALSE)
367	29x	for (i in vars) {
368	66x	if (is.factor(data[[i]])) {
369	58x	attr(data[[i]], "levels") <- levels(droplevels(data[[i]]))
370		}
371		}
372	29x	forms <- h_coxreg_univar_formulas(variables, interaction = control$interaction)
373	29x	mod <- lapply(
374	29x	forms, function(x) {
375	62x	survival::coxph(formula = stats::as.formula(x), data = data, ties = control$ties)
376		}
377		)
378	29x	structure(
379	29x	list(
380	29x	mod = mod,
381	29x	data = data,
382	29x	control = control,
383	29x	vars = variables,
384	29x	at = at
385		),
386	29x	class = "coxreg.univar"
387		)
388		}
389
390		#' @describeIn fit_coxreg Fit a multivariate Cox regression model.
391		#'
392		#' @return
393		#' * `fit_coxreg_multivar()` returns a `coxreg.multivar` class object which is a named list
394		#' with 4 elements:
395		#' * `mod`: Cox regression model fitted by [survival::coxph()].
396		#' * `data`: The original data frame input.
397		#' * `control`: The original control input.
398		#' * `vars`: The variables used in the model.
399		#'
400		#' @examples
401		#' # fit_coxreg_multivar
402		#'
403		#' ## Cox regression: multivariate Cox regression.
404		#' multivar_model <- fit_coxreg_multivar(
405		#' variables = list(
406		#' time = "time", event = "status", arm = "armcd",
407		#' covariates = c("covar1", "covar2")
408		#' ),
409		#' data = dta_bladder
410		#' )
411		#'
412		#' # Example without treatment arm.
413		#' multivar_covs_model <- fit_coxreg_multivar(
414		#' variables = list(
415		#' time = "time", event = "status",
416		#' covariates = c("covar1", "covar2")
417		#' ),
418		#' data = dta_bladder
419		#' )
420		#'
421		#' @export
422		fit_coxreg_multivar <- function(variables,
423		data,
424		control = control_coxreg()) {
425	51x	checkmate::assert_list(variables, names = "named")
426	51x	has_arm <- "arm" %in% names(variables)
427	51x	arm_name <- if (has_arm) "arm" else NULL
428
429	51x	if (!is.null(variables$covariates)) {
430	13x	checkmate::assert_character(variables$covariates)
431		}
432
433	51x	checkmate::assert_false(control$interaction)
434	51x	assert_df_with_variables(data, variables)
435	51x	assert_list_of_variables(variables[c(arm_name, "event", "time")])
436
437	51x	if (!is.null(variables$strata)) {
438	3x	checkmate::assert_disjunct(control$pval_method, "likelihood")
439		}
440
441	50x	form <- h_coxreg_multivar_formula(variables)
442	50x	mod <- survival::coxph(
443	50x	formula = stats::as.formula(form),
444	50x	data = data,
445	50x	ties = control$ties
446		)
447	50x	structure(
448	50x	list(
449	50x	mod = mod,
450	50x	data = data,
451	50x	control = control,
452	50x	vars = variables
453		),
454	50x	class = "coxreg.multivar"
455		)
456		}
457
458		#' Muffled `car::Anova`
459		#'
460		#' Applied on survival models, [car::Anova()] signal that the `strata` terms is dropped from the model formula when
461		#' present, this function deliberately muffles this message.
462		#'
463		#' @param mod (`coxph`)\cr Cox regression model fitted by [survival::coxph()].
464		#' @param test_statistic (`string`)\cr the method used for estimation of p.values; `wald` (default) or `likelihood`.
465		#'
466		#' @return Returns the output of [car::Anova()], with convergence message muffled.
467		#'
468		#' @keywords internal
469		muffled_car_anova <- function(mod, test_statistic) {
470	134x	tryCatch(
471	134x	withCallingHandlers(
472	134x	expr = {
473	134x	car::Anova(
474	134x	mod,
475	134x	test.statistic = test_statistic,
476	134x	type = "III"
477		)
478		},
479	134x	message = function(m) invokeRestart("muffleMessage"),
480	134x	error = function(e) {
481	1x	stop(paste(
482	1x	"the model seems to have convergence problems, please try to change",
483	1x	"the configuration of covariates or strata variables, e.g.",
484	1x	"- original error:", e
485		))
486		}
487		)
488		)
489		}

1		#' Cox Regression Helper: Interactions
2		#'
3		#' @description `r lifecycle::badge("stable")`
4		#'
5		#' Test and estimate the effect of a treatment in interaction with a covariate.
6		#' The effect is estimated as the HR of the tested treatment for a given level
7		#' of the covariate, in comparison to the treatment control.
8		#'
9		#' @inheritParams argument_convention
10		#' @param x (`numeric` or `factor`)\cr the values of the effect to be tested.
11		#' @param effect (`string`)\cr the name of the effect to be tested and estimated.
12		#' @param covar (`string`)\cr the name of the covariate in the model.
13		#' @param mod (`coxph`)\cr the Cox regression model.
14		#' @param label (`string`)\cr the label to be returned as `term_label`.
15		#' @param control (`list`)\cr a list of controls as returned by [control_coxreg()].
16		#' @param ... see methods.
17		#'
18		#' @examples
19		#' library(survival)
20		#'
21		#' set.seed(1, kind = "Mersenne-Twister")
22		#'
23		#' # Testing dataset [survival::bladder].
24		#' dta_bladder <- with(
25		#' data = bladder[bladder$enum < 5, ],
26		#' data.frame(
27		#' time = stop,
28		#' status = event,
29		#' armcd = as.factor(rx),
30		#' covar1 = as.factor(enum),
31		#' covar2 = factor(
32		#' sample(as.factor(enum)),
33		#' levels = 1:4,
34		#' labels = c("F", "F", "M", "M")
35		#' )
36		#' )
37		#' )
38		#' labels <- c("armcd" = "ARM", "covar1" = "A Covariate Label", "covar2" = "Sex (F/M)")
39		#' formatters::var_labels(dta_bladder)[names(labels)] <- labels
40		#' dta_bladder$age <- sample(20:60, size = nrow(dta_bladder), replace = TRUE)
41		#'
42		#' plot(
43		#' survfit(Surv(time, status) ~ armcd + covar1, data = dta_bladder),
44		#' lty = 2:4,
45		#' xlab = "Months",
46		#' col = c("blue1", "blue2", "blue3", "blue4", "red1", "red2", "red3", "red4")
47		#' )
48		#'
49		#' @name cox_regression_inter
50		NULL
51
52		#' @describeIn cox_regression_inter S3 generic helper function to determine interaction effect.
53		#'
54		#' @return
55		#' * `h_coxreg_inter_effect()` returns a `data.frame` of covariate interaction effects consisting of the following
56		#' variables: `effect`, `term`, `term_label`, `level`, `n`, `hr`, `lcl`, `ucl`, `pval`, and `pval_inter`.
57		#'
58		#' @export
59		h_coxreg_inter_effect <- function(x,
60		effect,
61		covar,
62		mod,
63		label,
64		control,
65		...) {
66	16x	UseMethod("h_coxreg_inter_effect", x)
67		}
68
69		#' @describeIn cox_regression_inter Estimate the interaction with a `numeric` covariate.
70		#'
71		#' @param at (`list`)\cr a list with items named after the covariate, every
72		#' item is a vector of levels at which the interaction should be estimated.
73		#'
74		#' @export
75		h_coxreg_inter_effect.numeric <- function(x,
76		effect,
77		covar,
78		mod,
79		label,
80		control,
81		at,
82		...) {
83	7x	betas <- stats::coef(mod)
84	7x	attrs <- attr(stats::terms(mod), "term.labels")
85	7x	term_indices <- grep(
86	7x	pattern = effect,
87	7x	x = attrs[!grepl("strata\\(", attrs)]
88		)
89	7x	checkmate::assert_vector(term_indices, len = 2)
90	7x	betas <- betas[term_indices]
91	7x	betas_var <- diag(stats::vcov(mod))[term_indices]
92	7x	betas_cov <- stats::vcov(mod)[term_indices[1], term_indices[2]]
93	7x	xval <- if (is.null(at[[covar]])) {
94	6x	stats::median(x)
95		} else {
96	1x	at[[covar]]
97		}
98	7x	effect_index <- !grepl(covar, names(betas))
99	7x	coef_hat <- betas[effect_index] + xval * betas[!effect_index]
100	7x	coef_se <- sqrt(
101	7x	betas_var[effect_index] +
102	7x	xval ^ 2 * betas_var[!effect_index] + # styler: off
103	7x	2 * xval * betas_cov
104		)
105	7x	q_norm <- stats::qnorm((1 + control$conf_level) / 2)
106	7x	data.frame(
107	7x	effect = "Covariate:",
108	7x	term = rep(covar, length(xval)),
109	7x	term_label = paste0(" ", xval),
110	7x	level = as.character(xval),
111	7x	n = NA,
112	7x	hr = exp(coef_hat),
113	7x	lcl = exp(coef_hat - q_norm * coef_se),
114	7x	ucl = exp(coef_hat + q_norm * coef_se),
115	7x	pval = NA,
116	7x	pval_inter = NA,
117	7x	stringsAsFactors = FALSE
118		)
119		}
120
121		#' @describeIn cox_regression_inter Estimate the interaction with a `factor` covariate.
122		#'
123		#' @param data (`data.frame`)\cr the data frame on which the model was fit.
124		#'
125		#' @export
126		h_coxreg_inter_effect.factor <- function(x,
127		effect,
128		covar,
129		mod,
130		label,
131		control,
132		data,
133		...) {
134	9x	y <- h_coxreg_inter_estimations(
135	9x	variable = effect, given = covar,
136	9x	lvl_var = levels(data[[effect]]),
137	9x	lvl_given = levels(data[[covar]]),
138	9x	mod = mod,
139	9x	conf_level = 0.95
140	9x	)[[1]]
141
142	9x	data.frame(
143	9x	effect = "Covariate:",
144	9x	term = rep(covar, nrow(y)),
145	9x	term_label = as.character(paste0(" ", levels(data[[covar]]))),
146	9x	level = as.character(levels(data[[covar]])),
147	9x	n = NA,
148	9x	hr = y[, "hr"],
149	9x	lcl = y[, "lcl"],
150	9x	ucl = y[, "ucl"],
151	9x	pval = NA,
152	9x	pval_inter = NA,
153	9x	stringsAsFactors = FALSE
154		)
155		}
156
157		#' @describeIn cox_regression_inter A higher level function to get
158		#' the results of the interaction test and the estimated values.
159		#'
160		#' @return
161		#' * `h_coxreg_extract_interaction()` returns the result of an interaction test and the estimated values. If
162		#' no interaction, [h_coxreg_univar_extract()] is applied instead.
163		#'
164		#' @examples
165		#' mod <- coxph(Surv(time, status) ~ armcd * covar1, data = dta_bladder)
166		#' h_coxreg_extract_interaction(
167		#' mod = mod, effect = "armcd", covar = "covar1", data = dta_bladder,
168		#' control = control_coxreg()
169		#' )
170		#'
171		#' @export
172		h_coxreg_extract_interaction <- function(effect,
173		covar,
174		mod,
175		data,
176		at,
177		control) {
178	21x	if (!any(attr(stats::terms(mod), "order") == 2)) {
179	8x	y <- h_coxreg_univar_extract(
180	8x	effect = effect, covar = covar, mod = mod, data = data, control = control
181		)
182	8x	y$pval_inter <- NA
183	8x	y
184		} else {
185	13x	test_statistic <- c(wald = "Wald", likelihood = "LR")[control$pval_method]
186
187		# Test the main treatment effect.
188	13x	mod_aov <- muffled_car_anova(mod, test_statistic)
189	13x	sum_anova <- broom::tidy(mod_aov)
190	13x	pval <- sum_anova[sum_anova$term == effect, ][["p.value"]]
191
192		# Test the interaction effect.
193	13x	pval_inter <- sum_anova[grep(":", sum_anova$term), ][["p.value"]]
194	13x	covar_test <- data.frame(
195	13x	effect = "Covariate:",
196	13x	term = covar,
197	13x	term_label = unname(labels_or_names(data[covar])),
198	13x	level = "",
199	13x	n = mod$n, hr = NA, lcl = NA, ucl = NA, pval = pval,
200	13x	pval_inter = pval_inter,
201	13x	stringsAsFactors = FALSE
202		)
203		# Estimate the interaction.
204	13x	y <- h_coxreg_inter_effect(
205	13x	data[[covar]],
206	13x	covar = covar,
207	13x	effect = effect,
208	13x	mod = mod,
209	13x	label = unname(labels_or_names(data[covar])),
210	13x	at = at,
211	13x	control = control,
212	13x	data = data
213		)
214	13x	rbind(covar_test, y)
215		}
216		}
217
218		#' @describeIn cox_regression_inter Hazard ratio estimation in interactions.
219		#'
220		#' @param variable,given (`string`)\cr the name of variables in interaction. We seek the estimation
221		#' of the levels of `variable` given the levels of `given`.
222		#' @param lvl_var,lvl_given (`character`)\cr corresponding levels has given by [levels()].
223		#' @param mod (`coxph`)\cr a fitted Cox regression model (see [survival::coxph()]).
224		#'
225		#' @details Given the cox regression investigating the effect of Arm (A, B, C; reference A)
226		#' and Sex (F, M; reference Female) and the model being abbreviated: y ~ Arm + Sex + Arm:Sex.
227		#' The cox regression estimates the coefficients along with a variance-covariance matrix for:
228		#'
229		#' - b1 (arm b), b2 (arm c)
230		#' - b3 (sex m)
231		#' - b4 (arm b: sex m), b5 (arm c: sex m)
232		#'
233		#' The estimation of the Hazard Ratio for arm C/sex M is given in reference
234		#' to arm A/Sex M by exp(b2 + b3 + b5)/ exp(b3) = exp(b2 + b5).
235		#' The interaction coefficient is deduced by b2 + b5 while the standard error
236		#' is obtained as $sqrt(Var b2 + Var b5 + 2 * covariance (b2,b5))$.
237		#'
238		#' @return
239		#' * `h_coxreg_inter_estimations()` returns a list of matrices (one per level of variable) with rows corresponding
240		#' to the combinations of `variable` and `given`, with columns:
241		#' * `coef_hat`: Estimation of the coefficient.
242		#' * `coef_se`: Standard error of the estimation.
243		#' * `hr`: Hazard ratio.
244		#' * `lcl, ucl`: Lower/upper confidence limit of the hazard ratio.
245		#'
246		#' @examples
247		#' mod <- coxph(Surv(time, status) ~ armcd * covar1, data = dta_bladder)
248		#' result <- h_coxreg_inter_estimations(
249		#' variable = "armcd", given = "covar1",
250		#' lvl_var = levels(dta_bladder$armcd),
251		#' lvl_given = levels(dta_bladder$covar1),
252		#' mod = mod, conf_level = .95
253		#' )
254		#' result
255		#'
256		#' @export
257		h_coxreg_inter_estimations <- function(variable,
258		given,
259		lvl_var,
260		lvl_given,
261		mod,
262		conf_level = 0.95) {
263	10x	var_lvl <- paste0(variable, lvl_var[-1]) # [-1]: reference level
264	10x	giv_lvl <- paste0(given, lvl_given)
265	10x	design_mat <- expand.grid(variable = var_lvl, given = giv_lvl)
266	10x	design_mat <- design_mat[order(design_mat$variable, design_mat$given), ]
267	10x	design_mat <- within(
268	10x	data = design_mat,
269	10x	expr = {
270	10x	inter <- paste0(variable, ":", given)
271	10x	rev_inter <- paste0(given, ":", variable)
272		}
273		)
274	10x	split_by_variable <- design_mat$variable
275	10x	interaction_names <- paste(design_mat$variable, design_mat$given, sep = "/")
276
277	10x	mmat <- stats::model.matrix(mod)[1, ]
278	10x	mmat[!mmat == 0] <- 0
279
280	10x	design_mat <- apply(
281	10x	X = design_mat, MARGIN = 1, FUN = function(x) {
282	32x	mmat[names(mmat) %in% x[-which(names(x) == "given")]] <- 1
283	32x	mmat
284		}
285		)
286	10x	colnames(design_mat) <- interaction_names
287
288	10x	coef <- stats::coef(mod)
289	10x	vcov <- stats::vcov(mod)
290	10x	betas <- as.matrix(coef)
291	10x	coef_hat <- t(design_mat) %*% betas
292	10x	dimnames(coef_hat)[2] <- "coef"
293	10x	coef_se <- apply(
294	10x	design_mat, 2,
295	10x	function(x) {
296	32x	vcov_el <- as.logical(x)
297	32x	y <- vcov[vcov_el, vcov_el]
298	32x	y <- sum(y)
299	32x	y <- sqrt(y)
300	32x	return(y)
301		}
302		)
303	10x	q_norm <- stats::qnorm((1 + conf_level) / 2)
304	10x	y <- cbind(coef_hat, `se(coef)` = coef_se)
305	10x	y <- apply(y, 1, function(x) {
306	32x	x["hr"] <- exp(x["coef"])
307	32x	x["lcl"] <- exp(x["coef"] - q_norm * x["se(coef)"])
308	32x	x["ucl"] <- exp(x["coef"] + q_norm * x["se(coef)"])
309	32x	x
310		})
311	10x	y <- t(y)
312	10x	y <- by(y, split_by_variable, identity)
313	10x	y <- lapply(y, as.matrix)
314	10x	attr(y, "details") <- paste0(
315	10x	"Estimations of ", variable,
316	10x	" hazard ratio given the level of ", given, " compared to ",
317	10x	variable, " level ", lvl_var[1], "."
318		)
319	10x	y
320		}

1		#' Encode Categorical Missing Values in a Data Frame
2		#'
3		#' @description `r lifecycle::badge("stable")`
4		#'
5		#' This is a helper function to encode missing entries across groups of categorical
6		#' variables in a data frame.
7		#'
8		#' @details Missing entries are those with `NA` or empty strings and will
9		#' be replaced with a specified value. If factor variables include missing
10		#' values, the missing value will be inserted as the last level.
11		#' Similarly, in case character or logical variables should be converted to factors
12		#' with the `char_as_factor` or `logical_as_factor` options, the missing values will
13		#' be set as the last level.
14		#'
15		#' @param data (`data.frame`)\cr data set.
16		#' @param omit_columns (`character`)\cr names of variables from `data` that should
17		#' not be modified by this function.
18		#' @param char_as_factor (`flag`)\cr whether to convert character variables
19		#' in `data` to factors.
20		#' @param logical_as_factor (`flag`)\cr whether to convert logical variables
21		#' in `data` to factors.
22		#' @param na_level (`string`)\cr used to replace all `NA` or empty
23		#' values inside non-`omit_columns` columns.
24		#'
25		#' @return A `data.frame` with the chosen modifications applied.
26		#'
27		#' @seealso [sas_na()] and [explicit_na()] for other missing data helper functions.
28		#'
29		#' @examples
30		#' my_data <- data.frame(
31		#' u = c(TRUE, FALSE, NA, TRUE),
32		#' v = factor(c("A", NA, NA, NA), levels = c("Z", "A")),
33		#' w = c("A", "B", NA, "C"),
34		#' x = c("D", "E", "F", NA),
35		#' y = c("G", "H", "I", ""),
36		#' z = c(1, 2, 3, 4),
37		#' stringsAsFactors = FALSE
38		#' )
39		#'
40		#' # Example 1
41		#' # Encode missing values in all character or factor columns.
42		#' df_explicit_na(my_data)
43		#' # Also convert logical columns to factor columns.
44		#' df_explicit_na(my_data, logical_as_factor = TRUE)
45		#' # Encode missing values in a subset of columns.
46		#' df_explicit_na(my_data, omit_columns = c("x", "y"))
47		#'
48		#' # Example 2
49		#' # Here we purposefully convert all `M` values to `NA` in the `SEX` variable.
50		#' # After running `df_explicit_na` the `NA` values are encoded as `<Missing>` but they are not
51		#' # included when generating `rtables`.
52		#' adsl <- tern_ex_adsl
53		#' adsl$SEX[adsl$SEX == "M"] <- NA
54		#' adsl <- df_explicit_na(adsl)
55		#'
56		#' # If you want the `Na` values to be displayed in the table use the `na_level` argument.
57		#' adsl <- tern_ex_adsl
58		#' adsl$SEX[adsl$SEX == "M"] <- NA
59		#' adsl <- df_explicit_na(adsl, na_level = "Missing Values")
60		#'
61		#' # Example 3
62		#' # Numeric variables that have missing values are not altered. This means that any `NA` value in
63		#' # a numeric variable will not be included in the summary statistics, nor will they be included
64		#' # in the denominator value for calculating the percent values.
65		#' adsl <- tern_ex_adsl
66		#' adsl$AGE[adsl$AGE < 30] <- NA
67		#' adsl <- df_explicit_na(adsl)
68		#'
69		#' @export
70		df_explicit_na <- function(data,
71		omit_columns = NULL,
72		char_as_factor = TRUE,
73		logical_as_factor = FALSE,
74		na_level = "<Missing>") {
75	27x	checkmate::assert_character(omit_columns, null.ok = TRUE, min.len = 1, any.missing = FALSE)
76	26x	checkmate::assert_data_frame(data)
77	25x	checkmate::assert_flag(char_as_factor)
78	24x	checkmate::assert_flag(logical_as_factor)
79	24x	checkmate::assert_string(na_level)
80
81	22x	target_vars <- if (is.null(omit_columns)) {
82	20x	names(data)
83		} else {
84	2x	setdiff(names(data), omit_columns) # May have duplicates.
85		}
86	22x	if (length(target_vars) == 0) {
87	1x	return(data)
88		}
89
90	21x	l_target_vars <- split(target_vars, target_vars)
91
92		# Makes sure target_vars exist in data and names are not duplicated.
93	21x	assert_df_with_variables(data, l_target_vars)
94
95	21x	for (x in target_vars) {
96	514x	xi <- data[[x]]
97	514x	xi_label <- obj_label(xi)
98
99		# Determine whether to convert character or logical input.
100	514x	do_char_conversion <- is.character(xi) && char_as_factor
101	514x	do_logical_conversion <- is.logical(xi) && logical_as_factor
102
103		# Pre-convert logical to character to deal correctly with replacing NA
104		# values below.
105	514x	if (do_logical_conversion) {
106	2x	xi <- as.character(xi)
107		}
108
109	514x	if (is.factor(xi) \|\| is.character(xi)) {
110		# Handle empty strings and NA values.
111	387x	xi <- explicit_na(sas_na(xi), label = na_level)
112
113		# Convert to factors if requested for the original type,
114		# set na_level as the last value.
115	387x	if (do_char_conversion \|\| do_logical_conversion) {
116	81x	levels_xi <- setdiff(sort(unique(xi)), na_level)
117	81x	if (na_level %in% unique(xi)) {
118	21x	levels_xi <- c(levels_xi, na_level)
119		}
120
121	81x	xi <- factor(xi, levels = levels_xi)
122		}
123
124	387x	data[, x] <- formatters::with_label(xi, label = xi_label)
125		}
126		}
127	21x	return(data)
128		}

1		#' `rtables` Access Helper Functions
2		#'
3		#' @description `r lifecycle::badge("stable")`
4		#'
5		#' These are a couple of functions that help with accessing the data in `rtables` objects.
6		#' Currently these work for occurrence tables, which are defined as having a count as the first
7		#' element and a fraction as the second element in each cell.
8		#'
9		#' @seealso [prune_occurrences] for usage of these functions.
10		#'
11		#' @name rtables_access
12		NULL
13
14		#' @describeIn rtables_access Helper function to extract the first values from each content
15		#' cell and from specified columns in a `TableRow`. Defaults to all columns.
16		#'
17		#' @param table_row (`TableRow`)\cr an analysis row in a occurrence table.
18		#' @param col_names (`character`)\cr the names of the columns to extract from.
19		#' @param col_indices (`integer`)\cr the indices of the columns to extract from. If `col_names` are provided,
20		#' then these are inferred from the names of `table_row`. Note that this currently only works well with a single
21		#' column split.
22		#'
23		#' @return
24		#' * `h_row_first_values()` returns a `vector` of numeric values.
25		#'
26		#' @examples
27		#' tbl <- basic_table() %>%
28		#' split_cols_by("ARM") %>%
29		#' split_rows_by("RACE") %>%
30		#' analyze("AGE", function(x) {
31		#' list(
32		#' "mean (sd)" = rcell(c(mean(x), sd(x)), format = "xx.x (xx.x)"),
33		#' "n" = length(x),
34		#' "frac" = rcell(c(0.1, 0.1), format = "xx (xx)")
35		#' )
36		#' }) %>%
37		#' build_table(tern_ex_adsl) %>%
38		#' prune_table()
39		#' tree_row_elem <- collect_leaves(tbl[2, ])[[1]]
40		#' result <- max(h_row_first_values(tree_row_elem))
41		#' result
42		#'
43		#' @export
44		h_row_first_values <- function(table_row,
45		col_names = NULL,
46		col_indices = NULL) {
47	727x	col_indices <- check_names_indices(table_row, col_names, col_indices)
48	727x	checkmate::assert_integerish(col_indices)
49	727x	checkmate::assert_subset(col_indices, seq_len(ncol(table_row)))
50
51		# Main values are extracted
52	727x	row_vals <- row_values(table_row)[col_indices]
53
54		# Main return
55	727x	vapply(row_vals, function(rv) {
56	2066x	if (is.null(rv)) {
57	727x	NA_real_
58		} else {
59	2063x	rv[1L]
60		}
61	727x	}, FUN.VALUE = numeric(1))
62		}
63
64		#' @describeIn rtables_access Helper function that extracts row values and checks if they are
65		#' convertible to integers (`integerish` values).
66		#'
67		#' @return
68		#' * `h_row_counts()` returns a `vector` of numeric values.
69		#'
70		#' @examples
71		#' # Row counts (integer values)
72		#' \dontrun{
73		#' h_row_counts(tree_row_elem) # Fails because there are no integers
74		#' }
75		#' # Using values with integers
76		#' tree_row_elem <- collect_leaves(tbl[3, ])[[1]]
77		#' result <- h_row_counts(tree_row_elem)
78		#' # result
79		#'
80		#' @export
81		h_row_counts <- function(table_row,
82		col_names = NULL,
83		col_indices = NULL) {
84	727x	counts <- h_row_first_values(table_row, col_names, col_indices)
85	727x	checkmate::assert_integerish(counts)
86	727x	counts
87		}
88
89		#' @describeIn rtables_access helper function to extract fractions from specified columns in a `TableRow`.
90		#' More specifically it extracts the second values from each content cell and checks it is a fraction.
91		#'
92		#' @return
93		#' * `h_row_fractions()` returns a `vector` of proportions.
94		#'
95		#' @examples
96		#' # Row fractions
97		#' tree_row_elem <- collect_leaves(tbl[4, ])[[1]]
98		#' h_row_fractions(tree_row_elem)
99		#'
100		#' @export
101		h_row_fractions <- function(table_row,
102		col_names = NULL,
103		col_indices = NULL) {
104	243x	col_indices <- check_names_indices(table_row, col_names, col_indices)
105	243x	row_vals <- row_values(table_row)[col_indices]
106	243x	fractions <- sapply(row_vals, "[", 2L)
107	243x	checkmate::assert_numeric(fractions, lower = 0, upper = 1)
108	243x	fractions
109		}
110
111		#' @describeIn rtables_access Helper function to extract column counts from specified columns in a table.
112		#'
113		#' @param table (`VTableNodeInfo`)\cr an occurrence table or row.
114		#'
115		#' @return
116		#' * `h_col_counts()` returns a `vector` of column counts.
117		#'
118		#' @export
119		h_col_counts <- function(table,
120		col_names = NULL,
121		col_indices = NULL) {
122	304x	col_indices <- check_names_indices(table, col_names, col_indices)
123	304x	counts <- col_counts(table)[col_indices]
124	304x	stats::setNames(counts, col_names)
125		}
126
127		#' @describeIn rtables_access Helper function to get first row of content table of current table.
128		#'
129		#' @return
130		#' * `h_content_first_row()` returns a row from an `rtables` table.
131		#'
132		#' @export
133		h_content_first_row <- function(table) {
134	27x	ct <- content_table(table)
135	27x	tree_children(ct)[[1]]
136		}
137
138		#' @describeIn rtables_access Helper function which says whether current table is a leaf in the tree.
139		#'
140		#' @return
141		#' * `is_leaf_table()` returns a `logical` value indicating whether current table is a leaf.
142		#'
143		#' @keywords internal
144		is_leaf_table <- function(table) {
145	168x	children <- tree_children(table)
146	168x	child_classes <- unique(sapply(children, class))
147	168x	identical(child_classes, "ElementaryTable")
148		}
149
150		#' @describeIn rtables_access Internal helper function that tests standard inputs for column indices.
151		#'
152		#' @return
153		#' * `check_names_indices` returns column indices.
154		#'
155		#' @keywords internal
156		check_names_indices <- function(table_row,
157		col_names = NULL,
158		col_indices = NULL) {
159	1274x	if (!is.null(col_names)) {
160	1231x	if (!is.null(col_indices)) {
161	!	stop(
162	!	"Inserted both col_names and col_indices when selecting row values. ",
163	!	"Please choose one."
164		)
165		}
166	1231x	col_indices <- h_col_indices(table_row, col_names)
167		}
168	1274x	if (is.null(col_indices)) {
169	37x	ll <- ifelse(is.null(ncol(table_row)), length(table_row), ncol(table_row))
170	37x	col_indices <- seq_len(ll)
171		}
172
173	1274x	return(col_indices)
174		}

1		#' Helper Functions for Multivariate Logistic Regression
2		#'
3		#' @description `r lifecycle::badge("stable")`
4		#'
5		#' Helper functions used in calculations for logistic regression.
6		#'
7		#' @inheritParams argument_convention
8		#' @param fit_glm (`glm`)\cr logistic regression model fitted by [stats::glm()] with "binomial" family.
9		#' Limited functionality is also available for conditional logistic regression models fitted by
10		#' [survival::clogit()], currently this is used only by [extract_rsp_biomarkers()].
11		#' @param x (`string` or `character`)\cr a variable or interaction term in `fit_glm` (depending on the
12		#' helper function).
13		#'
14		#' @examples
15		#' library(dplyr)
16		#' library(broom)
17		#'
18		#' adrs_f <- tern_ex_adrs %>%
19		#' filter(PARAMCD == "BESRSPI") %>%
20		#' filter(RACE %in% c("ASIAN", "WHITE", "BLACK OR AFRICAN AMERICAN")) %>%
21		#' mutate(
22		#' Response = case_when(AVALC %in% c("PR", "CR") ~ 1, TRUE ~ 0),
23		#' RACE = factor(RACE),
24		#' SEX = factor(SEX)
25		#' )
26		#' formatters::var_labels(adrs_f) <- c(formatters::var_labels(tern_ex_adrs), Response = "Response")
27		#' mod1 <- fit_logistic(
28		#' data = adrs_f,
29		#' variables = list(
30		#' response = "Response",
31		#' arm = "ARMCD",
32		#' covariates = c("AGE", "RACE")
33		#' )
34		#' )
35		#' mod2 <- fit_logistic(
36		#' data = adrs_f,
37		#' variables = list(
38		#' response = "Response",
39		#' arm = "ARMCD",
40		#' covariates = c("AGE", "RACE"),
41		#' interaction = "AGE"
42		#' )
43		#' )
44		#'
45		#' @name h_logistic_regression
46		NULL
47
48		#' @describeIn h_logistic_regression Helper function to extract interaction variable names from a fitted
49		#' model assuming only one interaction term.
50		#'
51		#' @return Vector of names of interaction variables.
52		#'
53		#' @export
54		h_get_interaction_vars <- function(fit_glm) {
55	27x	checkmate::assert_class(fit_glm, "glm")
56	27x	terms_name <- attr(stats::terms(fit_glm), "term.labels")
57	27x	terms_order <- attr(stats::terms(fit_glm), "order")
58	27x	interaction_term <- terms_name[terms_order == 2]
59	27x	checkmate::assert_string(interaction_term)
60	27x	strsplit(interaction_term, split = ":")[[1]]
61		}
62
63		#' @describeIn h_logistic_regression Helper function to get the right coefficient name from the
64		#' interaction variable names and the given levels. The main value here is that the order
65		#' of first and second variable is checked in the `interaction_vars` input.
66		#'
67		#' @param interaction_vars (`character` of length 2)\cr interaction variable names.
68		#' @param first_var_with_level (`character` of length 2)\cr the first variable name with
69		#' the interaction level.
70		#' @param second_var_with_level (`character` of length 2)\cr the second variable name with
71		#' the interaction level.
72		#'
73		#' @return Name of coefficient.
74		#'
75		#' @export
76		h_interaction_coef_name <- function(interaction_vars,
77		first_var_with_level,
78		second_var_with_level) {
79	45x	checkmate::assert_character(interaction_vars, len = 2, any.missing = FALSE)
80	45x	checkmate::assert_character(first_var_with_level, len = 2, any.missing = FALSE)
81	45x	checkmate::assert_character(second_var_with_level, len = 2, any.missing = FALSE)
82	45x	checkmate::assert_subset(c(first_var_with_level[1], second_var_with_level[1]), interaction_vars)
83
84	45x	first_name <- paste(first_var_with_level, collapse = "")
85	45x	second_name <- paste(second_var_with_level, collapse = "")
86	45x	if (first_var_with_level[1] == interaction_vars[1]) {
87	34x	paste(first_name, second_name, sep = ":")
88	11x	} else if (second_var_with_level[1] == interaction_vars[1]) {
89	11x	paste(second_name, first_name, sep = ":")
90		}
91		}
92
93		#' @describeIn h_logistic_regression Helper function to calculate the odds ratio estimates
94		#' for the case when both the odds ratio and the interaction variable are categorical.
95		#'
96		#' @param odds_ratio_var (`string`)\cr the odds ratio variable.
97		#' @param interaction_var (`string`)\cr the interaction variable.
98		#'
99		#' @return Odds ratio.
100		#'
101		#' @export
102		h_or_cat_interaction <- function(odds_ratio_var,
103		interaction_var,
104		fit_glm,
105		conf_level = 0.95) {
106	7x	interaction_vars <- h_get_interaction_vars(fit_glm)
107	7x	checkmate::assert_string(odds_ratio_var)
108	7x	checkmate::assert_string(interaction_var)
109	7x	checkmate::assert_subset(c(odds_ratio_var, interaction_var), interaction_vars)
110	7x	checkmate::assert_vector(interaction_vars, len = 2)
111
112	7x	xs_level <- fit_glm$xlevels
113	7x	xs_coef <- stats::coef(fit_glm)
114	7x	xs_vcov <- stats::vcov(fit_glm)
115	7x	y <- list()
116	7x	for (var_level in xs_level[[odds_ratio_var]][-1]) {
117	12x	x <- list()
118	12x	for (ref_level in xs_level[[interaction_var]]) {
119	32x	coef_names <- paste0(odds_ratio_var, var_level)
120	32x	if (ref_level != xs_level[[interaction_var]][1]) {
121	20x	interaction_coef_name <- h_interaction_coef_name(
122	20x	interaction_vars,
123	20x	c(odds_ratio_var, var_level),
124	20x	c(interaction_var, ref_level)
125		)
126	20x	coef_names <- c(
127	20x	coef_names,
128	20x	interaction_coef_name
129		)
130		}
131	32x	if (length(coef_names) > 1) {
132	20x	ones <- t(c(1, 1))
133	20x	est <- as.numeric(ones %*% xs_coef[coef_names])
134	20x	se <- sqrt(as.numeric(ones %% xs_vcov[coef_names, coef_names] %% t(ones)))
135		} else {
136	12x	est <- xs_coef[coef_names]
137	12x	se <- sqrt(as.numeric(xs_vcov[coef_names, coef_names]))
138		}
139	32x	or <- exp(est)
140	32x	ci <- exp(est + c(lcl = -1, ucl = 1) * stats::qnorm((1 + conf_level) / 2) * se)
141	32x	x[[ref_level]] <- list(or = or, ci = ci)
142		}
143	12x	y[[var_level]] <- x
144		}
145	7x	y
146		}
147
148		#' @describeIn h_logistic_regression Helper function to calculate the odds ratio estimates
149		#' for the case when either the odds ratio or the interaction variable is continuous.
150		#'
151		#' @param at (`NULL` or `numeric`)\cr optional values for the interaction variable. Otherwise
152		#' the median is used.
153		#'
154		#' @return Odds ratio.
155		#'
156		#' @note We don't provide a function for the case when both variables are continuous because
157		#' this does not arise in this table, as the treatment arm variable will always be involved
158		#' and categorical.
159		#'
160		#' @export
161		h_or_cont_interaction <- function(odds_ratio_var,
162		interaction_var,
163		fit_glm,
164		at = NULL,
165		conf_level = 0.95) {
166	9x	interaction_vars <- h_get_interaction_vars(fit_glm)
167	9x	checkmate::assert_string(odds_ratio_var)
168	9x	checkmate::assert_string(interaction_var)
169	9x	checkmate::assert_subset(c(odds_ratio_var, interaction_var), interaction_vars)
170	9x	checkmate::assert_vector(interaction_vars, len = 2)
171	9x	checkmate::assert_numeric(at, min.len = 1, null.ok = TRUE, any.missing = FALSE)
172	9x	xs_level <- fit_glm$xlevels
173	9x	xs_coef <- stats::coef(fit_glm)
174	9x	xs_vcov <- stats::vcov(fit_glm)
175	9x	xs_class <- attr(fit_glm$terms, "dataClasses")
176	9x	model_data <- fit_glm$model
177	9x	if (!is.null(at)) {
178	2x	checkmate::assert_set_equal(xs_class[interaction_var], "numeric")
179		}
180	9x	y <- list()
181	9x	if (xs_class[interaction_var] == "numeric") {
182	6x	if (is.null(at)) {
183	4x	at <- ceiling(stats::median(model_data[[interaction_var]]))
184		}
185
186	6x	for (var_level in xs_level[[odds_ratio_var]][-1]) {
187	12x	x <- list()
188	12x	for (increment in at) {
189	18x	coef_names <- paste0(odds_ratio_var, var_level)
190	18x	if (increment != 0) {
191	18x	interaction_coef_name <- h_interaction_coef_name(
192	18x	interaction_vars,
193	18x	c(odds_ratio_var, var_level),
194	18x	c(interaction_var, "")
195		)
196	18x	coef_names <- c(
197	18x	coef_names,
198	18x	interaction_coef_name
199		)
200		}
201	18x	if (length(coef_names) > 1) {
202	18x	xvec <- t(c(1, increment))
203	18x	est <- as.numeric(xvec %*% xs_coef[coef_names])
204	18x	se <- sqrt(as.numeric(xvec %% xs_vcov[coef_names, coef_names] %% t(xvec)))
205		} else {
206	!	est <- xs_coef[coef_names]
207	!	se <- sqrt(as.numeric(xs_vcov[coef_names, coef_names]))
208		}
209	18x	or <- exp(est)
210	18x	ci <- exp(est + c(lcl = -1, ucl = 1) * stats::qnorm((1 + conf_level) / 2) * se)
211	18x	x[[as.character(increment)]] <- list(or = or, ci = ci)
212		}
213	12x	y[[var_level]] <- x
214		}
215		} else {
216	3x	checkmate::assert_set_equal(xs_class[odds_ratio_var], "numeric")
217	3x	checkmate::assert_set_equal(xs_class[interaction_var], "factor")
218	3x	for (var_level in xs_level[[interaction_var]]) {
219	9x	coef_names <- odds_ratio_var
220	9x	if (var_level != xs_level[[interaction_var]][1]) {
221	6x	interaction_coef_name <- h_interaction_coef_name(
222	6x	interaction_vars,
223	6x	c(odds_ratio_var, ""),
224	6x	c(interaction_var, var_level)
225		)
226	6x	coef_names <- c(
227	6x	coef_names,
228	6x	interaction_coef_name
229		)
230		}
231	9x	if (length(coef_names) > 1) {
232	6x	xvec <- t(c(1, 1))
233	6x	est <- as.numeric(xvec %*% xs_coef[coef_names])
234	6x	se <- sqrt(as.numeric(xvec %% xs_vcov[coef_names, coef_names] %% t(xvec)))
235		} else {
236	3x	est <- xs_coef[coef_names]
237	3x	se <- sqrt(as.numeric(xs_vcov[coef_names, coef_names]))
238		}
239	9x	or <- exp(est)
240	9x	ci <- exp(est + c(lcl = -1, ucl = 1) * stats::qnorm((1 + conf_level) / 2) * se)
241	9x	y[[var_level]] <- list(or = or, ci = ci)
242		}
243		}
244	9x	y
245		}
246
247		#' @describeIn h_logistic_regression Helper function to calculate the odds ratio estimates
248		#' in case of an interaction. This is a wrapper for [h_or_cont_interaction()] and
249		#' [h_or_cat_interaction()].
250		#'
251		#' @return Odds ratio.
252		#'
253		#' @export
254		h_or_interaction <- function(odds_ratio_var,
255		interaction_var,
256		fit_glm,
257		at = NULL,
258		conf_level = 0.95) {
259	13x	xs_class <- attr(fit_glm$terms, "dataClasses")
260	13x	if (any(xs_class[c(odds_ratio_var, interaction_var)] == "numeric")) {
261	7x	h_or_cont_interaction(
262	7x	odds_ratio_var,
263	7x	interaction_var,
264	7x	fit_glm,
265	7x	at = at,
266	7x	conf_level = conf_level
267		)
268	6x	} else if (all(xs_class[c(odds_ratio_var, interaction_var)] == "factor")) {
269	6x	h_or_cat_interaction(
270	6x	odds_ratio_var,
271	6x	interaction_var,
272	6x	fit_glm,
273	6x	conf_level = conf_level
274		)
275		} else {
276	!	stop("wrong interaction variable class, the interaction variable is not a numeric nor a factor")
277		}
278		}
279
280		#' @describeIn h_logistic_regression Helper function to construct term labels from simple terms and the table
281		#' of numbers of patients.
282		#'
283		#' @param terms (`character`)\cr simple terms.
284		#' @param table (`table`)\cr table containing numbers for terms.
285		#'
286		#' @return Term labels containing numbers of patients.
287		#'
288		#' @export
289		h_simple_term_labels <- function(terms,
290		table) {
291	45x	checkmate::assert_true(is.table(table))
292	45x	checkmate::assert_multi_class(terms, classes = c("factor", "character"))
293	45x	terms <- as.character(terms)
294	45x	term_n <- table[terms]
295	45x	paste0(terms, ", n = ", term_n)
296		}
297
298		#' @describeIn h_logistic_regression Helper function to construct term labels from interaction terms and the table
299		#' of numbers of patients.
300		#'
301		#' @param terms1 (`character`)\cr terms for first dimension (rows).
302		#' @param terms2 (`character`)\cr terms for second dimension (rows).
303		#' @param any (`flag`)\cr whether any of `term1` and `term2` can be fulfilled to count the
304		#' number of patients. In that case they can only be scalar (strings).
305		#'
306		#' @return Term labels containing numbers of patients.
307		#'
308		#' @export
309		h_interaction_term_labels <- function(terms1,
310		terms2,
311		table,
312		any = FALSE) {
313	8x	checkmate::assert_true(is.table(table))
314	8x	checkmate::assert_flag(any)
315	8x	checkmate::assert_multi_class(terms1, classes = c("factor", "character"))
316	8x	checkmate::assert_multi_class(terms2, classes = c("factor", "character"))
317	8x	terms1 <- as.character(terms1)
318	8x	terms2 <- as.character(terms2)
319	8x	if (any) {
320	4x	checkmate::assert_scalar(terms1)
321	4x	checkmate::assert_scalar(terms2)
322	4x	paste0(
323	4x	terms1, " or ", terms2, ", n = ",
324		# Note that we double count in the initial sum the cell [terms1, terms2], therefore subtract.
325	4x	sum(c(table[terms1, ], table[, terms2])) - table[terms1, terms2]
326		)
327		} else {
328	4x	term_n <- table[cbind(terms1, terms2)]
329	4x	paste0(terms1, " * ", terms2, ", n = ", term_n)
330		}
331		}
332
333		#' @describeIn h_logistic_regression Helper function to tabulate the main effect
334		#' results of a (conditional) logistic regression model.
335		#'
336		#' @return Tabulated main effect results from a logistic regression model.
337		#'
338		#' @examples
339		#' h_glm_simple_term_extract("AGE", mod1)
340		#' h_glm_simple_term_extract("ARMCD", mod1)
341		#'
342		#' @export
343		h_glm_simple_term_extract <- function(x, fit_glm) {
344	61x	checkmate::assert_multi_class(fit_glm, c("glm", "clogit"))
345	61x	checkmate::assert_string(x)
346
347	61x	xs_class <- attr(fit_glm$terms, "dataClasses")
348	61x	xs_level <- fit_glm$xlevels
349	61x	xs_coef <- summary(fit_glm)$coefficients
350	61x	stats <- if (inherits(fit_glm, "glm")) {
351	49x	c("estimate" = "Estimate", "std_error" = "Std. Error", "pvalue" = "Pr(>\|z\|)")
352		} else {
353	12x	c("estimate" = "coef", "std_error" = "se(coef)", "pvalue" = "Pr(>\|z\|)")
354		}
355		# Make sure x is not an interaction term.
356	61x	checkmate::assert_subset(x, names(xs_class))
357	61x	x_sel <- if (xs_class[x] == "numeric") x else paste0(x, xs_level[[x]][-1])
358	61x	x_stats <- as.data.frame(xs_coef[x_sel, stats, drop = FALSE], stringsAsFactors = FALSE)
359	61x	colnames(x_stats) <- names(stats)
360	61x	x_stats$estimate <- as.list(x_stats$estimate)
361	61x	x_stats$std_error <- as.list(x_stats$std_error)
362	61x	x_stats$pvalue <- as.list(x_stats$pvalue)
363	61x	x_stats$df <- as.list(1)
364	61x	if (xs_class[x] == "numeric") {
365	46x	x_stats$term <- x
366	46x	x_stats$term_label <- if (inherits(fit_glm, "glm")) {
367	34x	formatters::var_labels(fit_glm$data[x], fill = TRUE)
368		} else {
369		# We just fill in here with the `term` itself as we don't have the data available.
370	12x	x
371		}
372	46x	x_stats$is_variable_summary <- FALSE
373	46x	x_stats$is_term_summary <- TRUE
374		} else {
375	15x	checkmate::assert_class(fit_glm, "glm")
376		# The reason is that we don't have the original data set in the `clogit` object
377		# and therefore cannot determine the `x_numbers` here.
378	15x	x_numbers <- table(fit_glm$data[[x]])
379	15x	x_stats$term <- xs_level[[x]][-1]
380	15x	x_stats$term_label <- h_simple_term_labels(x_stats$term, x_numbers)
381	15x	x_stats$is_variable_summary <- FALSE
382	15x	x_stats$is_term_summary <- TRUE
383	15x	main_effects <- car::Anova(fit_glm, type = 3, test.statistic = "Wald")
384	15x	x_main <- data.frame(
385	15x	pvalue = main_effects[x, "Pr(>Chisq)", drop = TRUE],
386	15x	term = xs_level[[x]][1],
387	15x	term_label = paste("Reference", h_simple_term_labels(xs_level[[x]][1], x_numbers)),
388	15x	df = main_effects[x, "Df", drop = TRUE],
389	15x	stringsAsFactors = FALSE
390		)
391	15x	x_main$pvalue <- as.list(x_main$pvalue)
392	15x	x_main$df <- as.list(x_main$df)
393	15x	x_main$estimate <- list(numeric(0))
394	15x	x_main$std_error <- list(numeric(0))
395	15x	if (length(xs_level[[x]][-1]) == 1) {
396	6x	x_main$pvalue <- list(numeric(0))
397	6x	x_main$df <- list(numeric(0))
398		}
399	15x	x_main$is_variable_summary <- TRUE
400	15x	x_main$is_term_summary <- FALSE
401	15x	x_stats <- rbind(x_main, x_stats)
402		}
403	61x	x_stats$variable <- x
404	61x	x_stats$variable_label <- if (inherits(fit_glm, "glm")) {
405	49x	formatters::var_labels(fit_glm$data[x], fill = TRUE)
406		} else {
407	12x	x
408		}
409	61x	x_stats$interaction <- ""
410	61x	x_stats$interaction_label <- ""
411	61x	x_stats$reference <- ""
412	61x	x_stats$reference_label <- ""
413	61x	rownames(x_stats) <- NULL
414	61x	x_stats[c(
415	61x	"variable",
416	61x	"variable_label",
417	61x	"term",
418	61x	"term_label",
419	61x	"interaction",
420	61x	"interaction_label",
421	61x	"reference",
422	61x	"reference_label",
423	61x	"estimate",
424	61x	"std_error",
425	61x	"df",
426	61x	"pvalue",
427	61x	"is_variable_summary",
428	61x	"is_term_summary"
429		)]
430		}
431
432		#' @describeIn h_logistic_regression Helper function to tabulate the interaction term
433		#' results of a logistic regression model.
434		#'
435		#' @return Tabulated interaction term results from a logistic regression model.
436		#'
437		#' @examples
438		#' h_glm_interaction_extract("ARMCD:AGE", mod2)
439		#'
440		#' @export
441		h_glm_interaction_extract <- function(x, fit_glm) {
442	6x	vars <- h_get_interaction_vars(fit_glm)
443	6x	xs_class <- attr(fit_glm$terms, "dataClasses")
444
445	6x	checkmate::assert_string(x)
446
447		# Only take two-way interaction
448	6x	checkmate::assert_vector(vars, len = 2)
449
450		# Only consider simple case: first variable in interaction is arm, a categorical variable
451	6x	checkmate::assert_disjunct(xs_class[vars[1]], "numeric")
452
453	6x	xs_level <- fit_glm$xlevels
454	6x	xs_coef <- summary(fit_glm)$coefficients
455	6x	main_effects <- car::Anova(fit_glm, type = 3, test.statistic = "Wald")
456	6x	stats <- c("estimate" = "Estimate", "std_error" = "Std. Error", "pvalue" = "Pr(>\|z\|)")
457	6x	v1_comp <- xs_level[[vars[1]]][-1]
458	6x	if (xs_class[vars[2]] == "numeric") {
459	3x	x_stats <- as.data.frame(
460	3x	xs_coef[paste0(vars[1], v1_comp, ":", vars[2]), stats, drop = FALSE],
461	3x	stringsAsFactors = FALSE
462		)
463	3x	colnames(x_stats) <- names(stats)
464	3x	x_stats$term <- v1_comp
465	3x	x_numbers <- table(fit_glm$data[[vars[1]]])
466	3x	x_stats$term_label <- h_simple_term_labels(v1_comp, x_numbers)
467	3x	v1_ref <- xs_level[[vars[1]]][1]
468	3x	term_main <- v1_ref
469	3x	ref_label <- h_simple_term_labels(v1_ref, x_numbers)
470	3x	} else if (xs_class[vars[2]] != "numeric") {
471	3x	v2_comp <- xs_level[[vars[2]]][-1]
472	3x	v1_v2_grid <- expand.grid(v1 = v1_comp, v2 = v2_comp)
473	3x	x_sel <- paste(
474	3x	paste0(vars[1], v1_v2_grid$v1),
475	3x	paste0(vars[2], v1_v2_grid$v2),
476	3x	sep = ":"
477		)
478	3x	x_stats <- as.data.frame(xs_coef[x_sel, stats, drop = FALSE], stringsAsFactors = FALSE)
479	3x	colnames(x_stats) <- names(stats)
480	3x	x_stats$term <- paste(v1_v2_grid$v1, "*", v1_v2_grid$v2)
481	3x	x_numbers <- table(fit_glm$data[[vars[1]]], fit_glm$data[[vars[2]]])
482	3x	x_stats$term_label <- h_interaction_term_labels(v1_v2_grid$v1, v1_v2_grid$v2, x_numbers)
483	3x	v1_ref <- xs_level[[vars[1]]][1]
484	3x	v2_ref <- xs_level[[vars[2]]][1]
485	3x	term_main <- paste(vars[1], vars[2], sep = " * ")
486	3x	ref_label <- h_interaction_term_labels(v1_ref, v2_ref, x_numbers, any = TRUE)
487		}
488	6x	x_stats$df <- as.list(1)
489	6x	x_stats$pvalue <- as.list(x_stats$pvalue)
490	6x	x_stats$is_variable_summary <- FALSE
491	6x	x_stats$is_term_summary <- TRUE
492	6x	x_main <- data.frame(
493	6x	pvalue = main_effects[x, "Pr(>Chisq)", drop = TRUE],
494	6x	term = term_main,
495	6x	term_label = paste("Reference", ref_label),
496	6x	df = main_effects[x, "Df", drop = TRUE],
497	6x	stringsAsFactors = FALSE
498		)
499	6x	x_main$pvalue <- as.list(x_main$pvalue)
500	6x	x_main$df <- as.list(x_main$df)
501	6x	x_main$estimate <- list(numeric(0))
502	6x	x_main$std_error <- list(numeric(0))
503	6x	x_main$is_variable_summary <- TRUE
504	6x	x_main$is_term_summary <- FALSE
505
506	6x	x_stats <- rbind(x_main, x_stats)
507	6x	x_stats$variable <- x
508	6x	x_stats$variable_label <- paste(
509	6x	"Interaction of",
510	6x	formatters::var_labels(fit_glm$data[vars[1]], fill = TRUE),
511		"*",
512	6x	formatters::var_labels(fit_glm$data[vars[2]], fill = TRUE)
513		)
514	6x	x_stats$interaction <- ""
515	6x	x_stats$interaction_label <- ""
516	6x	x_stats$reference <- ""
517	6x	x_stats$reference_label <- ""
518	6x	rownames(x_stats) <- NULL
519	6x	x_stats[c(
520	6x	"variable",
521	6x	"variable_label",
522	6x	"term",
523	6x	"term_label",
524	6x	"interaction",
525	6x	"interaction_label",
526	6x	"reference",
527	6x