tern coverage - 95.08%

Files
Source

#' 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)
}

#' Odds Ratio Estimation
#'
#' @description `r lifecycle::badge("stable")`
#'
#' Compares bivariate responses between two groups in terms of odds ratios
#' along with a confidence interval.
#'
#' @inheritParams argument_convention
#'
#' @details This function uses either logistic regression for unstratified
#'   analyses, or conditional logistic regression for stratified analyses.
#'   The Wald confidence interval with the specified confidence level is
#'   calculated.
#'
#' @note For stratified analyses, there is currently no implementation for conditional
#'   likelihood confidence intervals, therefore the likelihood confidence interval is not
#'   yet available as an option. Besides, when `rsp` contains only responders or non-responders,
#'   then the result values will be `NA`, because no odds ratio estimation is possible.
#'
#' @seealso Relevant helper function [h_odds_ratio()].
#'
#' @name odds_ratio
NULL

#' @describeIn odds_ratio Statistics function which estimates the odds ratio
#'   between a treatment and a control. A `variables` list with `arm` and `strata`
#'   variable names must be passed if a stratified analysis is required.
#'
#' @inheritParams split_cols_by_groups
#'
#' @return
#' * `s_odds_ratio()` returns a named list with the statistics `or_ci`
#'   (containing `est`, `lcl`, and `ucl`) and `n_tot`.
#'
#' @examples
#' set.seed(12)
#' dta <- data.frame(
#'   rsp = sample(c(TRUE, FALSE), 100, TRUE),
#'   grp = factor(rep(c("A", "B"), each = 50), levels = c("B", "A")),
#'   strata = factor(sample(c("C", "D"), 100, TRUE))
#' )
#'
#' # Unstratified analysis.
#' s_odds_ratio(
#'   df = subset(dta, grp == "A"),
#'   .var = "rsp",
#'   .ref_group = subset(dta, grp == "B"),
#'   .in_ref_col = FALSE,
#'   .df_row = dta
#' )
#'
#' # Stratified analysis.
#' s_odds_ratio(
#'   df = subset(dta, grp == "A"),
#'   .var = "rsp",
#'   .ref_group = subset(dta, grp == "B"),
#'   .in_ref_col = FALSE,
#'   .df_row = dta,
#'   variables = list(arm = "grp", strata = "strata")
#' )
#'
#' @export
s_odds_ratio <- function(df,
                         .var,
                         .ref_group,
                         .in_ref_col,
                         .df_row,
                         variables = list(arm = NULL, strata = NULL),
                         conf_level = 0.95,
                         groups_list = NULL) {
  y <- list(or_ci = "", n_tot = "")

  if (!.in_ref_col) {
    assert_proportion_value(conf_level)
    assert_df_with_variables(df, list(rsp = .var))
    assert_df_with_variables(.ref_group, list(rsp = .var))

    if (is.null(variables$strata)) {
      data <- data.frame(
        rsp = c(.ref_group[[.var]], df[[.var]]),
        grp = factor(
          rep(c("ref", "Not-ref"), c(nrow(.ref_group), nrow(df))),
          levels = c("ref", "Not-ref")
        )
      )
      y <- or_glm(data, conf_level = conf_level)
    } else {
      assert_df_with_variables(.df_row, c(list(rsp = .var), variables))

      # The group variable prepared for clogit must be synchronised with combination groups definition.
      if (is.null(groups_list)) {
        ref_grp <- as.character(unique(.ref_group[[variables$arm]]))
        trt_grp <- as.character(unique(df[[variables$arm]]))
        grp <- stats::relevel(factor(.df_row[[variables$arm]]), ref = ref_grp)
      } else {
        # If more than one level in reference col.
        reference <- as.character(unique(.ref_group[[variables$arm]]))
        grp_ref_flag <- vapply(
          X = groups_list,
          FUN.VALUE = TRUE,
          FUN = function(x) all(reference %in% x)
        )
        ref_grp <- names(groups_list)[grp_ref_flag]

        # If more than one level in treatment col.
        treatment <- as.character(unique(df[[variables$arm]]))
        grp_trt_flag <- vapply(
          X = groups_list,
          FUN.VALUE = TRUE,
          FUN = function(x) all(treatment %in% x)
        )
        trt_grp <- names(groups_list)[grp_trt_flag]

        grp <- combine_levels(.df_row[[variables$arm]], levels = reference, new_level = ref_grp)
        grp <- combine_levels(grp, levels = treatment, new_level = trt_grp)
      }

      # The reference level in `grp` must be the same as in the `rtables` column split.
      data <- data.frame(
        rsp = .df_row[[.var]],
        grp = grp,
        strata = interaction(.df_row[variables$strata])
      )
      y_all <- or_clogit(data, conf_level = conf_level)
      checkmate::assert_string(trt_grp)
      checkmate::assert_subset(trt_grp, names(y_all$or_ci))
      y$or_ci <- y_all$or_ci[[trt_grp]]
      y$n_tot <- y_all$n_tot
    }
  }

  y$or_ci <- formatters::with_label(
    x = y$or_ci,
    label = paste0("Odds Ratio (", 100 * conf_level, "% CI)")
  )

  y$n_tot <- formatters::with_label(
    x = y$n_tot,
    label = "Total n"
  )

  y
}

#' @describeIn odds_ratio Formatted analysis function which is used as `afun` in `estimate_odds_ratio()`.
#'
#' @return
#' * `a_odds_ratio()` returns the corresponding list with formatted [rtables::CellValue()].
#'
#' @examples
#' a_odds_ratio(
#'   df = subset(dta, grp == "A"),
#'   .var = "rsp",
#'   .ref_group = subset(dta, grp == "B"),
#'   .in_ref_col = FALSE,
#'   .df_row = dta
#' )
#'
#' @export
a_odds_ratio <- make_afun(
  s_odds_ratio,
  .formats = c(or_ci = "xx.xx (xx.xx - xx.xx)"),
  .indent_mods = c(or_ci = 1L)
)

#' @describeIn odds_ratio 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_odds_ratio()`.
#'
#' @return
#' * `estimate_odds_ratio()` 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_odds_ratio()` to the table layout.
#'
#' @examples
#' dta <- data.frame(
#'   rsp = sample(c(TRUE, FALSE), 100, TRUE),
#'   grp = factor(rep(c("A", "B"), each = 50))
#' )
#'
#' l <- basic_table() %>%
#'   split_cols_by(var = "grp", ref_group = "B") %>%
#'   estimate_odds_ratio(vars = "rsp")
#'
#' build_table(l, df = dta)
#'
#' @export
estimate_odds_ratio <- function(lyt,
                                vars,
                                nested = TRUE,
                                ...,
                                show_labels = "hidden",
                                table_names = vars,
                                .stats = "or_ci",
                                .formats = NULL,
                                .labels = NULL,
                                .indent_mods = NULL) {
  afun <- make_afun(
    a_odds_ratio,
    .stats = .stats,
    .formats = .formats,
    .labels = .labels,
    .indent_mods = .indent_mods
  )

  analyze(
    lyt,
    vars,
    afun = afun,
    nested = nested,
    extra_args = list(...),
    show_labels = show_labels,
    table_names = table_names
  )
}

#' Helper Functions for Odds Ratio Estimation
#'
#' @description `r lifecycle::badge("stable")`
#'
#' Functions to calculate odds ratios in [estimate_odds_ratio()].
#'
#' @inheritParams argument_convention
#' @param data (`data.frame`)\cr data frame containing at least the variables `rsp` and `grp`, and optionally
#'   `strata` for [or_clogit()].
#'
#' @return A named `list` of elements `or_ci` and `n_tot`.
#'
#' @seealso [odds_ratio]
#'
#' @name h_odds_ratio
NULL

#' @describeIn h_odds_ratio Estimates the odds ratio based on [stats::glm()]. Note that there must be
#'   exactly 2 groups in `data` as specified by the `grp` variable.
#'
#' @examples
#' # Data with 2 groups.
#' data <- data.frame(
#'   rsp = as.logical(c(1, 1, 0, 1, 0, 0, 1, 1)),
#'   grp = letters[c(1, 1, 1, 2, 2, 2, 1, 2)],
#'   strata = letters[c(1, 2, 1, 2, 2, 2, 1, 2)],
#'   stringsAsFactors = TRUE
#' )
#'
#' # Odds ratio based on glm.
#' or_glm(data, conf_level = 0.95)
#'
#' @export
or_glm <- function(data, conf_level) {
  checkmate::assert_logical(data$rsp)
  assert_proportion_value(conf_level)
  assert_df_with_variables(data, list(rsp = "rsp", grp = "grp"))
  checkmate::assert_multi_class(data$grp, classes = c("factor", "character"))

  data$grp <- as_factor_keep_attributes(data$grp)
  assert_df_with_factors(data, list(val = "grp"), min.levels = 2, max.levels = 2)
  formula <- stats::as.formula("rsp ~ grp")
  model_fit <- stats::glm(
    formula = formula, data = data,
    family = stats::binomial(link = "logit")
  )

  # Note that here we need to discard the intercept.
  or <- exp(stats::coef(model_fit)[-1])
  or_ci <- exp(
    stats::confint.default(model_fit, level = conf_level)[-1, , drop = FALSE]
  )

  values <- stats::setNames(c(or, or_ci), c("est", "lcl", "ucl"))
  n_tot <- stats::setNames(nrow(model_fit$model), "n_tot")

  list(or_ci = values, n_tot = n_tot)
}

#' @describeIn h_odds_ratio estimates the odds ratio based on [survival::clogit()]. This is done for
#'   the whole data set including all groups, since the results are not the same as when doing
#'   pairwise comparisons between the groups.
#'
#' @examples
#' # Data with 3 groups.
#' data <- data.frame(
#'   rsp = as.logical(c(1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0)),
#'   grp = letters[c(1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3)],
#'   strata = LETTERS[c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2)],
#'   stringsAsFactors = TRUE
#' )
#'
#' # Odds ratio based on stratified estimation by conditional logistic regression.
#' or_clogit(data, conf_level = 0.95)
#'
#' @export
or_clogit <- function(data, conf_level) {
  checkmate::assert_logical(data$rsp)
  assert_proportion_value(conf_level)
  assert_df_with_variables(data, list(rsp = "rsp", grp = "grp", strata = "strata"))
  checkmate::assert_multi_class(data$grp, classes = c("factor", "character"))
  checkmate::assert_multi_class(data$strata, classes = c("factor", "character"))

  data$grp <- as_factor_keep_attributes(data$grp)
  data$strata <- as_factor_keep_attributes(data$strata)

  # Deviation from convention: `survival::strata` must be simply `strata`.
  formula <- stats::as.formula("rsp ~ grp + strata(strata)")
  model_fit <- clogit_with_tryCatch(formula = formula, data = data)

  # Create a list with one set of OR estimates and CI per coefficient, i.e.
  # comparison of one group vs. the reference group.
  coef_est <- stats::coef(model_fit)
  ci_est <- stats::confint(model_fit, level = conf_level)
  or_ci <- list()
  for (coef_name in names(coef_est)) {
    grp_name <- gsub("^grp", "", x = coef_name)
    or_ci[[grp_name]] <- stats::setNames(
      object = exp(c(coef_est[coef_name], ci_est[coef_name, , drop = TRUE])),
      nm = c("est", "lcl", "ucl")
    )
  }
  list(or_ci = or_ci, n_tot = c(n_tot = model_fit$n))
}

#' Incidence Rate
#'
#' @description `r lifecycle::badge("stable")`
#'
#' Estimate the event rate adjusted for person-years at risk, otherwise known
#' as incidence rate. Primary analysis variable is the person-years at risk.
#'
#' @inheritParams argument_convention
#' @param control (`list`)\cr parameters for estimation details, specified by using
#'   the helper function [control_incidence_rate()]. Possible parameter options are:
#'   * `conf_level` (`proportion`)\cr confidence level for the estimated incidence rate.
#'   * `conf_type` (`string`)\cr `normal` (default), `normal_log`, `exact`, or `byar`
#'     for confidence interval type.
#'   * `input_time_unit` (`string`)\cr `day`, `week`, `month`, or `year` (default)
#'     indicating time unit for data input.
#'   * `num_pt_year` (`numeric`)\cr time unit for desired output (in person-years).
#' @param n_events (`integer`)\cr number of events observed.
#'
#' @seealso [control_incidence_rate()] and helper functions [h_incidence_rate].
#'
#' @name incidence_rate
NULL

#' @describeIn incidence_rate Statistics function which estimates the incidence rate and the
#'   associated confidence interval.
#'
#' @return
#' * `s_incidence_rate()` returns the following statistics:
#'   - `person_years`: Total person-years at risk.
#'   - `n_events`: Total number of events observed.
#'   - `rate`: Estimated incidence rate.
#'   - `rate_ci`: Confidence interval for the incidence rate.
#'
#' @examples
#' library(dplyr)
#'
#' df <- data.frame(
#'   USUBJID = as.character(seq(6)),
#'   CNSR = c(0, 1, 1, 0, 0, 0),
#'   AVAL = c(10.1, 20.4, 15.3, 20.8, 18.7, 23.4),
#'   ARM = factor(c("A", "A", "A", "B", "B", "B"))
#' ) %>%
#'   mutate(is_event = CNSR == 0) %>%
#'   mutate(n_events = as.integer(is_event))
#'
#' @keywords internal
s_incidence_rate <- function(df,
                             .var,
                             n_events,
                             is_event,
                             control = control_incidence_rate()) {
  if (!missing(is_event)) {
    warning("argument is_event will be deprecated. Please use n_events.")

    if (missing(n_events)) {
      assert_df_with_variables(df, list(tte = .var, is_event = is_event))
      checkmate::assert_string(.var)
      checkmate::assert_logical(df[[is_event]], any.missing = FALSE)
      checkmate::assert_numeric(df[[.var]], any.missing = FALSE)
      n_events <- is_event
    }
  } else {
    assert_df_with_variables(df, list(tte = .var, n_events = n_events))
    checkmate::assert_string(.var)
    checkmate::assert_numeric(df[[.var]], any.missing = FALSE)
    checkmate::assert_integer(df[[n_events]], any.missing = FALSE)
  }

  input_time_unit <- control$input_time_unit
  num_pt_year <- control$num_pt_year
  conf_level <- control$conf_level
  person_years <- sum(df[[.var]], na.rm = TRUE) * (
    1 * (input_time_unit == "year") +
      1 / 12 * (input_time_unit == "month") +
      1 / 52.14 * (input_time_unit == "week") +
      1 / 365.24 * (input_time_unit == "day")
  )
  n_events <- sum(df[[n_events]], na.rm = TRUE)

  result <- h_incidence_rate(
    person_years,
    n_events,
    control
  )
  list(
    person_years = formatters::with_label(person_years, "Total patient-years at risk"),
    n_events = formatters::with_label(n_events, "Number of adverse events observed"),
    rate = formatters::with_label(result$rate, paste("AE rate per", num_pt_year, "patient-years")),
    rate_ci = formatters::with_label(result$rate_ci, f_conf_level(conf_level))
  )
}

#' @describeIn incidence_rate Formatted analysis function which is used as `afun`
#'   in `estimate_incidence_rate()`.
#'
#' @return
#' * `a_incidence_rate()` returns the corresponding list with formatted [rtables::CellValue()].
#'
#'
#' @keywords internal
a_incidence_rate <- make_afun(
  s_incidence_rate,
  .formats = c(
    "person_years" = "xx.x",
    "n_events" = "xx",
    "rate" = "xx.xx",
    "rate_ci" = "(xx.xx, xx.xx)"
  )
)

#' @describeIn incidence_rate Layout-creating function which can take statistics function arguments
#'   and additional format arguments. This function is a wrapper for [rtables::analyze()].
#'
#' @return
#' * `estimate_incidence_rate()` 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_incidence_rate()` to the table layout.
#'
#' @examples
#' basic_table() %>%
#'   split_cols_by("ARM") %>%
#'   add_colcounts() %>%
#'   estimate_incidence_rate(
#'     vars = "AVAL",
#'     n_events = "n_events",
#'     control = control_incidence_rate(
#'       input_time_unit = "month",
#'       num_pt_year = 100
#'     )
#'   ) %>%
#'   build_table(df)
#'
#' @export
estimate_incidence_rate <- function(lyt,
                                    vars,
                                    nested = TRUE,
                                    ...,
                                    show_labels = "hidden",
                                    table_names = vars,
                                    .stats = NULL,
                                    .formats = NULL,
                                    .labels = NULL,
                                    .indent_mods = NULL) {
  afun <- make_afun(
    a_incidence_rate,
    .stats = .stats,
    .formats = .formats,
    .labels = .labels,
    .indent_mods = .indent_mods
  )

  analyze(
    lyt,
    vars,
    show_labels = show_labels,
    table_names = table_names,
    afun = afun,
    nested = nested,
    extra_args = list(...)
  )
}

#' Helper Functions for Incidence Rate
#'
#' @description `r lifecycle::badge("stable")`
#'
#' @param control (`list`)\cr parameters for estimation details, specified by using
#'   the helper function [control_incidence_rate()]. Possible parameter options are:
#'   * `conf_level`: (`proportion`)\cr confidence level for the estimated incidence rate.
#'   * `conf_type`: (`string`)\cr `normal` (default), `normal_log`, `exact`, or `byar`
#'     for confidence interval type.
#'   * `input_time_unit`: (`string`)\cr `day`, `week`, `month`, or `year` (default)
#'     indicating time unit for data input.
#'   * `num_pt_year`: (`numeric`)\cr time unit for desired output (in person-years).
#' @param person_years (`numeric`)\cr total person-years at risk.
#' @param alpha (`numeric`)\cr two-sided alpha-level for confidence interval.
#' @param n_events (`integer`)\cr number of events observed.
#'
#' @return Estimated incidence rate `rate` and associated confidence interval `rate_ci`.
#'
#' @seealso [incidence_rate]
#'
#' @name h_incidence_rate
NULL

#' @describeIn h_incidence_rate Helper function to estimate the incidence rate and
#'   associated confidence interval based on the normal approximation for the
#'   incidence rate. Unit is one person-year.
#'
#' @examples
#' h_incidence_rate_normal(200, 2)
#'
#' @export
h_incidence_rate_normal <- function(person_years,
                                    n_events,
                                    alpha = 0.05) {
  checkmate::assert_number(person_years)
  checkmate::assert_number(n_events)
  assert_proportion_value(alpha)

  est <- n_events / person_years
  se <- sqrt(est / person_years)
  ci <- est + c(-1, 1) * stats::qnorm(1 - alpha / 2) * se

  list(rate = est, rate_ci = ci)
}

#' @describeIn h_incidence_rate Helper function to estimate the incidence rate and
#'   associated confidence interval based on the normal approximation for the
#'   logarithm of the incidence rate. Unit is one person-year.
#'
#' @examples
#' h_incidence_rate_normal_log(200, 2)
#'
#' @export
h_incidence_rate_normal_log <- function(person_years,
                                        n_events,
                                        alpha = 0.05) {
  checkmate::assert_number(person_years)
  checkmate::assert_number(n_events)
  assert_proportion_value(alpha)

  rate_est <- n_events / person_years
  rate_se <- sqrt(rate_est / person_years)
  lrate_est <- log(rate_est)
  lrate_se <- rate_se / rate_est
  ci <- exp(lrate_est + c(-1, 1) * stats::qnorm(1 - alpha / 2) * lrate_se)

  list(rate = rate_est, rate_ci = ci)
}

#' @describeIn h_incidence_rate Helper function to estimate the incidence rate and
#'   associated exact confidence interval. Unit is one person-year.
#'
#' @examples
#' h_incidence_rate_exact(200, 2)
#'
#' @export
h_incidence_rate_exact <- function(person_years,
                                   n_events,
                                   alpha = 0.05) {
  checkmate::assert_number(person_years)
  checkmate::assert_number(n_events)
  assert_proportion_value(alpha)

  est <- n_events / person_years
  lcl <- stats::qchisq(p = (alpha) / 2, df = 2 * n_events) / (2 * person_years)
  ucl <- stats::qchisq(p = 1 - (alpha) / 2, df = 2 * n_events + 2) / (2 * person_years)

  list(rate = est, rate_ci = c(lcl, ucl))
}

#' @describeIn h_incidence_rate Helper function to estimate the incidence rate and
#'   associated `Byar`'s confidence interval. Unit is one person-year.
#'
#' @examples
#' h_incidence_rate_byar(200, 2)
#'
#' @export
h_incidence_rate_byar <- function(person_years,
                                  n_events,
                                  alpha = 0.05) {
  checkmate::assert_number(person_years)
  checkmate::assert_number(n_events)
  assert_proportion_value(alpha)

  est <- n_events / person_years
  seg_1 <- n_events + 0.5
  seg_2 <- 1 - 1 / (9 * (n_events + 0.5))
  seg_3 <- stats::qnorm(1 - alpha / 2) * sqrt(1 / (n_events + 0.5)) / 3
  lcl <- seg_1 * ((seg_2 - seg_3)^3) / person_years
  ucl <- seg_1 * ((seg_2 + seg_3) ^ 3) / person_years # styler: off

  list(rate = est, rate_ci = c(lcl, ucl))
}

#' @describeIn h_incidence_rate Helper function to estimate the incidence rate and
#'   associated confidence interval.
#'
#'
#' @keywords internal
h_incidence_rate <- function(person_years,
                             n_events,
                             control = control_incidence_rate()) {
  alpha <- 1 - control$conf_level
  est <- switch(control$conf_type,
    normal = h_incidence_rate_normal(person_years, n_events, alpha),
    normal_log = h_incidence_rate_normal_log(person_years, n_events, alpha),
    exact = h_incidence_rate_exact(person_years, n_events, alpha),
    byar = h_incidence_rate_byar(person_years, n_events, alpha)
  )

  num_pt_year <- control$num_pt_year
  list(
    rate = est$rate * num_pt_year,
    rate_ci = est$rate_ci * num_pt_year
  )
}

#' Kaplan-Meier Plot
#'
#' @description `r lifecycle::badge("stable")`
#'
#' From a survival model, a graphic is rendered along with tabulated annotation
#' including the number of patient at risk at given time and the median survival
#' per group.
#'
#' @inheritParams grid::gTree
#' @inheritParams argument_convention
#' @param df (`data.frame`)\cr data set containing all analysis variables.
#' @param variables (named `list`)\cr variable names. Details are:
#'   * `tte` (`numeric`)\cr variable indicating time-to-event duration values.
#'   * `is_event` (`logical`)\cr event variable. `TRUE` if event, `FALSE` if time to event is censored.
#'   * `arm` (`factor`)\cr the treatment group variable.
#'   * `strat` (`character` or `NULL`)\cr variable names indicating stratification factors.
#' @param control_surv (`list`)\cr parameters for comparison details, specified by using
#'   the helper function [control_surv_timepoint()]. Some possible parameter options are:
#'   * `conf_level` (`proportion`)\cr confidence level of the interval for survival rate.
#'   * `conf_type` (`string`)\cr `"plain"` (default), `"log"`, `"log-log"` for confidence interval type,
#'     see more in [survival::survfit()]. Note that the option "none" is no longer supported.
#' @param xticks (`numeric`, `number`, or `NULL`)\cr numeric vector of ticks or single number with spacing
#'   between ticks on the x axis. If `NULL` (default), [labeling::extended()] is used to determine
#'   an optimal tick position on the x axis.
#' @param yval (`string`)\cr value of y-axis. Options are `Survival` (default) and `Failure` probability.
#' @param censor_show (`flag`)\cr whether to show censored.
#' @param xlab (`string`)\cr label of x-axis.
#' @param ylab (`string`)\cr label of y-axis.
#' @param ylim (`vector` of `numeric`)\cr vector of length 2 containing lower and upper limits for the y-axis.
#'   If `NULL` (default), the minimum and maximum y-values displayed are used as limits.
#' @param title (`string`)\cr title for plot.
#' @param footnotes (`string`)\cr footnotes for plot.
#' @param col (`character`)\cr lines colors. Length of a vector should be equal
#'   to number of strata from [survival::survfit()].
#' @param lty (`numeric`)\cr line type. Length of a vector should be equal
#'   to number of strata from [survival::survfit()].
#' @param lwd (`numeric`)\cr line width. Length of a vector should be equal
#'   to number of strata from [survival::survfit()].
#' @param pch (`numeric`, `string`)\cr value or character of points symbol to indicate censored cases.
#' @param size (`numeric`)\cr size of censored point, a class of `unit`.
#' @param max_time (`numeric`)\cr maximum value to show on X axis. Only data values less than or up to
#'   this threshold value will be plotted (defaults to `NULL`).
#' @param font_size (`number`)\cr font size to be used.
#' @param ci_ribbon (`flag`)\cr draw the confidence interval around the Kaplan-Meier curve.
#' @param ggtheme (`theme`)\cr a graphical theme as provided by `ggplot2` to control outlook of the Kaplan-Meier curve.
#' @param annot_at_risk (`flag`)\cr compute and add the annotation table reporting the number of patient at risk
#'   matching the main grid of the Kaplan-Meier curve.
#' @param annot_at_risk_title (`flag`)\cr whether the "Patients at Risk" title should be added above the `annot_at_risk`
#'   table. Has no effect if `annot_at_risk` is `FALSE`. Defaults to `TRUE`.
#' @param annot_surv_med (`flag`)\cr compute and add the annotation table on the Kaplan-Meier curve estimating the
#'   median survival time per group.
#' @param annot_coxph (`flag`)\cr add the annotation table from a [survival::coxph()] model.
#' @param annot_stats (`string`)\cr statistics annotations to add to the plot. Options are
#'   `median` (median survival follow-up time) and `min` (minimum survival follow-up time).
#' @param annot_stats_vlines (`flag`)\cr add vertical lines corresponding to each of the statistics
#'   specified by `annot_stats`. If `annot_stats` is `NULL` no lines will be added.
#' @param control_coxph_pw (`list`)\cr parameters for comparison details, specified by using
#'   the helper function [control_coxph()]. Some possible parameter options are:
#'   * `pval_method` (`string`)\cr p-value method for testing hazard ratio = 1.
#'     Default method is `"log-rank"`, can also be set to `"wald"` or `"likelihood"`.
#'   * `ties` (`string`)\cr method for tie handling. Default is `"efron"`,
#'     can also be set to `"breslow"` or `"exact"`. See more in [survival::coxph()]
#'   * `conf_level` (`proportion`)\cr confidence level of the interval for HR.
#' @param position_coxph (`numeric`)\cr x and y positions for plotting [survival::coxph()] model.
#' @param position_surv_med (`numeric`)\cr x and y positions for plotting annotation table estimating median survival
#'   time per group.
#' @param width_annots (named `list` of `unit`s)\cr a named list of widths for annotation tables with names `surv_med`
#'   (median survival time table) and `coxph` ([survival::coxph()] model table), where each value is the width
#'   (in units) to implement when printing the annotation table.
#'
#' @return A `grob` of class `gTree`.
#'
#' @examples
#' \donttest{
#' library(dplyr)
#' library(ggplot2)
#' library(survival)
#' library(grid)
#' library(nestcolor)
#'
#' df <- tern_ex_adtte %>%
#'   filter(PARAMCD == "OS") %>%
#'   mutate(is_event = CNSR == 0)
#' variables <- list(tte = "AVAL", is_event = "is_event", arm = "ARMCD")
#'
#' # 1. Example - basic option
#'
#' res <- g_km(df = df, variables = variables)
#' res <- g_km(df = df, variables = variables, yval = "Failure")
#' res <- g_km(
#'   df = df,
#'   variables = variables,
#'   control_surv = control_surv_timepoint(conf_level = 0.9),
#'   col = c("grey25", "grey50", "grey75"),
#'   annot_at_risk_title = FALSE
#' )
#' res <- g_km(df = df, variables = variables, ggtheme = theme_minimal())
#' res <- g_km(df = df, variables = variables, ggtheme = theme_minimal(), lty = 1:3)
#' res <- g_km(df = df, variables = variables, max = 2000)
#' res <- g_km(
#'   df = df,
#'   variables = variables,
#'   annot_stats = c("min", "median"),
#'   annot_stats_vlines = TRUE
#' )
#'
#' # 2. Example - Arrange several KM curve on a single graph device
#'
#' # 2.1 Use case: A general graph on the top, a zoom on the bottom.
#' grid.newpage()
#' lyt <- grid.layout(nrow = 2, ncol = 1) %>%
#'   viewport(layout = .) %>%
#'   pushViewport()
#'
#' res <- g_km(
#'   df = df, variables = variables, newpage = FALSE, annot_surv_med = FALSE,
#'   vp = viewport(layout.pos.row = 1, layout.pos.col = 1)
#' )
#' res <- g_km(
#'   df = df, variables = variables, max = 1000, newpage = FALSE, annot_surv_med = FALSE,
#'   ggtheme = theme_dark(),
#'   vp = viewport(layout.pos.row = 2, layout.pos.col = 1)
#' )
#'
#' # 2.1 Use case: No annotations on top, annotated graph on bottom
#' grid.newpage()
#' lyt <- grid.layout(nrow = 2, ncol = 1) %>%
#'   viewport(layout = .) %>%
#'   pushViewport()
#'
#' res <- g_km(
#'   df = df, variables = variables, newpage = FALSE,
#'   annot_surv_med = FALSE, annot_at_risk = FALSE,
#'   vp = viewport(layout.pos.row = 1, layout.pos.col = 1)
#' )
#' res <- g_km(
#'   df = df, variables = variables, max = 2000, newpage = FALSE, annot_surv_med = FALSE,
#'   annot_at_risk = TRUE,
#'   ggtheme = theme_dark(),
#'   vp = viewport(layout.pos.row = 2, layout.pos.col = 1)
#' )
#'
#' # Add annotation from a pairwise coxph analysis
#' g_km(
#'   df = df, variables = variables,
#'   annot_coxph = TRUE
#' )
#'
#' # Change widths/sizes of surv_med and coxph annotation tables.
#' g_km(
#'   df = df, variables = c(variables, list(strat = "SEX")),
#'   annot_coxph = TRUE,
#'   width_annots = list(surv_med = grid::unit(2, "in"), coxph = grid::unit(3, "in"))
#' )
#'
#' g_km(
#'   df = df, variables = c(variables, list(strat = "SEX")),
#'   font_size = 15,
#'   annot_coxph = TRUE,
#'   control_coxph = control_coxph(pval_method = "wald", ties = "exact", conf_level = 0.99),
#'   position_coxph = c(0.5, 0.5)
#' )
#'
#' # Change position of the treatment group annotation table.
#' g_km(
#'   df = df, variables = c(variables, list(strat = "SEX")),
#'   font_size = 15,
#'   annot_coxph = TRUE,
#'   control_coxph = control_coxph(pval_method = "wald", ties = "exact", conf_level = 0.99),
#'   position_surv_med = c(1, 0.7)
#' )
#' }
#'
#' @export
g_km <- function(df,
                 variables,
                 control_surv = control_surv_timepoint(),
                 col = NULL,
                 lty = NULL,
                 lwd = .5,
                 censor_show = TRUE,
                 pch = 3,
                 size = 2,
                 max_time = NULL,
                 xticks = NULL,
                 xlab = "Days",
                 yval = c("Survival", "Failure"),
                 ylab = paste(yval, "Probability"),
                 ylim = NULL,
                 title = NULL,
                 footnotes = NULL,
                 draw = TRUE,
                 newpage = TRUE,
                 gp = NULL,
                 vp = NULL,
                 name = NULL,
                 font_size = 12,
                 ci_ribbon = FALSE,
                 ggtheme = nestcolor::theme_nest(),
                 annot_at_risk = TRUE,
                 annot_at_risk_title = TRUE,
                 annot_surv_med = TRUE,
                 annot_coxph = FALSE,
                 annot_stats = NULL,
                 annot_stats_vlines = FALSE,
                 control_coxph_pw = control_coxph(),
                 position_coxph = c(-0.03, -0.02),
                 position_surv_med = c(0.95, 0.9),
                 width_annots = list(surv_med = grid::unit(0.3, "npc"), coxph = grid::unit(0.4, "npc"))) {
  checkmate::assert_list(variables)
  checkmate::assert_subset(c("tte", "arm", "is_event"), names(variables))
  checkmate::assert_string(title, null.ok = TRUE)
  checkmate::assert_string(footnotes, null.ok = TRUE)
  checkmate::assert_character(col, null.ok = TRUE)
  checkmate::assert_subset(annot_stats, c("median", "min"))
  checkmate::assert_logical(annot_stats_vlines)
  checkmate::assert_true(all(sapply(width_annots, grid::is.unit)))

  tte <- variables$tte
  is_event <- variables$is_event
  arm <- variables$arm

  assert_valid_factor(df[[arm]])
  assert_df_with_variables(df, list(tte = tte, is_event = is_event, arm = arm))
  checkmate::assert_logical(df[[is_event]], min.len = 1, any.missing = FALSE)
  checkmate::assert_numeric(df[[tte]], min.len = 1, any.missing = FALSE)

  armval <- as.character(unique(df[[arm]]))
  if (annot_coxph && length(armval) < 2) {
    stop(paste(
      "When `annot_coxph` = TRUE, `df` must contain at least 2 levels of `variables$arm`",
      "in order to calculate the hazard ratio."
    ))
  } else if (length(armval) > 1) {
    armval <- NULL
  }
  yval <- match.arg(yval)
  formula <- stats::as.formula(paste0("survival::Surv(", tte, ", ", is_event, ") ~ ", arm))
  fit_km <- survival::survfit(
    formula = formula,
    data = df,
    conf.int = control_surv$conf_level,
    conf.type = control_surv$conf_type
  )
  data_plot <- h_data_plot(
    fit_km = fit_km,
    armval = armval,
    max_time = max_time
  )

  xticks <- h_xticks(data = data_plot, xticks = xticks, max_time = max_time)
  gg <- h_ggkm(
    data = data_plot,
    censor_show = censor_show,
    pch = pch,
    size = size,
    xticks = xticks,
    xlab = xlab,
    yval = yval,
    ylab = ylab,
    ylim = ylim,
    title = title,
    footnotes = footnotes,
    max_time = max_time,
    lwd = lwd,
    lty = lty,
    col = col,
    ggtheme = ggtheme,
    ci_ribbon = ci_ribbon
  )

  if (!is.null(annot_stats)) {
    if ("median" %in% annot_stats) {
      fit_km_all <- survival::survfit(
        formula = stats::as.formula(paste0("survival::Surv(", tte, ", ", is_event, ") ~ ", 1)),
        data = df,
        conf.int = control_surv$conf_level,
        conf.type = control_surv$conf_type
      )
      gg <- gg +
        geom_text(
          size = 8 / ggplot2::.pt, col = 1,
          x = stats::median(fit_km_all) + 0.065 * max(data_plot$time),
          y = ifelse(yval == "Survival", 0.62, 0.38),
          label = paste("Median F/U:\n", round(stats::median(fit_km_all), 1), tolower(df$AVALU[1]))
        )
      if (annot_stats_vlines) {
        gg <- gg +
          geom_segment(aes(x = stats::median(fit_km_all), xend = stats::median(fit_km_all), y = -Inf, yend = Inf),
            linetype = 2, col = "darkgray"
          )
      }
    }
    if ("min" %in% annot_stats) {
      min_fu <- min(df[[tte]])
      gg <- gg +
        geom_text(
          size = 8 / ggplot2::.pt, col = 1,
          x = min_fu + max(data_plot$time) * ifelse(yval == "Survival", 0.05, 0.07),
          y = ifelse(yval == "Survival", 1.0, 0.05),
          label = paste("Min. F/U:\n", round(min_fu, 1), tolower(df$AVALU[1]))
        )
      if (annot_stats_vlines) {
        gg <- gg +
          geom_segment(aes(x = min_fu, xend = min_fu, y = Inf, yend = -Inf), linetype = 2, col = "darkgray")
      }
    }
    gg <- gg + ggplot2::guides(fill = ggplot2::guide_legend(override.aes = list(shape = NA, label = "")))
  }

  g_el <- h_decompose_gg(gg)

  if (annot_at_risk) {
    # This is the content of the table that will be below the graph.
    annot_tbl <- summary(fit_km, time = xticks)
    annot_tbl <- if (is.null(fit_km$strata)) {
      data.frame(
        n.risk = annot_tbl$n.risk,
        time = annot_tbl$time,
        strata = as.factor(armval)
      )
    } else {
      strata_lst <- strsplit(sub("=", "equals", levels(annot_tbl$strata)), "equals")
      levels(annot_tbl$strata) <- matrix(unlist(strata_lst), ncol = 2, byrow = TRUE)[, 2]
      data.frame(
        n.risk = annot_tbl$n.risk,
        time = annot_tbl$time,
        strata = annot_tbl$strata
      )
    }

    grobs_patient <- h_grob_tbl_at_risk(
      data = data_plot,
      annot_tbl = annot_tbl,
      xlim = max(max_time, data_plot$time, xticks),
      title = annot_at_risk_title
    )
  }

  if (annot_at_risk || annot_surv_med || annot_coxph) {
    lyt <- h_km_layout(
      data = data_plot, g_el = g_el, title = title, footnotes = footnotes,
      annot_at_risk = annot_at_risk, annot_at_risk_title = annot_at_risk_title
    )
    at_risk_ttl <- as.numeric(annot_at_risk_title)
    ttl_row <- as.numeric(!is.null(title))
    foot_row <- as.numeric(!is.null(footnotes))
    km_grob <- grid::gTree(
      vp = grid::viewport(layout = lyt, height = .95, width = .95),
      children = grid::gList(
        # Title.
        if (ttl_row == 1) {
          grid::gTree(
            vp = grid::viewport(layout.pos.row = 1, layout.pos.col = 2),
            children = grid::gList(grid::textGrob(label = title, x = grid::unit(0, "npc"), hjust = 0))
          )
        },

        # The Kaplan - Meier curve (top-right corner).
        grid::gTree(
          vp = grid::viewport(layout.pos.row = 1 + ttl_row, layout.pos.col = 2),
          children = grid::gList(g_el$panel)
        ),

        # Survfit summary table (top-right corner).
        if (annot_surv_med) {
          grid::gTree(
            vp = grid::viewport(layout.pos.row = 1 + ttl_row, layout.pos.col = 2),
            children = h_grob_median_surv(
              fit_km = fit_km,
              armval = armval,
              x = position_surv_med[1],
              y = position_surv_med[2],
              width = if (!is.null(width_annots[["surv_med"]])) width_annots[["surv_med"]] else grid::unit(0.3, "npc"),
              ttheme = gridExtra::ttheme_default(base_size = font_size)
            )
          )
        },
        if (annot_coxph) {
          grid::gTree(
            vp = grid::viewport(layout.pos.row = 1 + ttl_row, layout.pos.col = 2),
            children = h_grob_coxph(
              df = df,
              variables = variables,
              control_coxph_pw = control_coxph_pw,
              x = position_coxph[1],
              y = position_coxph[2],
              width = if (!is.null(width_annots[["coxph"]])) width_annots[["coxph"]] else grid::unit(0.4, "npc"),
              ttheme = gridExtra::ttheme_default(
                base_size = font_size,
                padding = grid::unit(c(1, .5), "lines"),
                core = list(bg_params = list(fill = c("grey95", "grey90"), alpha = .5))
              )
            )
          )
        },

        # Add the y-axis annotation (top-left corner).
        grid::gTree(
          vp = grid::viewport(layout.pos.row = 1 + ttl_row, layout.pos.col = 1),
          children = h_grob_y_annot(ylab = g_el$ylab, yaxis = g_el$yaxis)
        ),

        # Add the x-axis annotation (second row below the Kaplan Meier Curve).
        grid::gTree(
          vp = grid::viewport(layout.pos.row = 2 + ttl_row, layout.pos.col = 2),
          children = grid::gList(rbind(g_el$xaxis, g_el$xlab))
        ),

        # Add the legend.
        grid::gTree(
          vp = grid::viewport(layout.pos.row = 3 + ttl_row, layout.pos.col = 2),
          children = grid::gList(g_el$guide)
        ),

        # Add the table with patient-at-risk numbers.
        if (annot_at_risk && annot_at_risk_title) {
          grid::gTree(
            vp = grid::viewport(layout.pos.row = 4 + ttl_row, layout.pos.col = 1),
            children = grobs_patient$title
          )
        },
        if (annot_at_risk) {
          grid::gTree(
            vp = grid::viewport(layout.pos.row = 4 + at_risk_ttl + ttl_row, layout.pos.col = 2),
            children = grobs_patient$at_risk
          )
        },
        if (annot_at_risk) {
          grid::gTree(
            vp = grid::viewport(layout.pos.row = 4 + at_risk_ttl + ttl_row, layout.pos.col = 1),
            children = grobs_patient$label
          )
        },
        if (annot_at_risk) {
          # Add the x-axis for the table.
          grid::gTree(
            vp = grid::viewport(layout.pos.row = 5 + at_risk_ttl + ttl_row, layout.pos.col = 2),
            children = grid::gList(rbind(g_el$xaxis, g_el$xlab))
          )
        },

        # Footnotes.
        if (foot_row == 1) {
          grid::gTree(
            vp = grid::viewport(
              layout.pos.row = ifelse(annot_at_risk, 6 + at_risk_ttl + ttl_row, 4 + ttl_row),
              layout.pos.col = 2
            ),
            children = grid::gList(grid::textGrob(label = footnotes, x = grid::unit(0, "npc"), hjust = 0))
          )
        }
      )
    )

    result <- grid::gTree(
      vp = vp,
      gp = gp,
      name = name,
      children = grid::gList(km_grob)
    )
  } else {
    result <- grid::gTree(
      vp = vp,
      gp = gp,
      name = name,
      children = grid::gList(ggplot2::ggplotGrob(gg))
    )
  }

  if (newpage && draw) grid::grid.newpage()
  if (draw) grid::grid.draw(result)
  invisible(result)
}

#' Helper function: tidy survival fit
#'
#' @description `r lifecycle::badge("stable")`
#'
#' Convert the survival fit data into a data frame designed for plotting
#' within `g_km`.
#'
#' This starts from the [broom::tidy()] result, and then:
#'   * Post-processes the `strata` column into a factor.
#'   * Extends each stratum by an additional first row with time 0 and probability 1 so that
#'     downstream plot lines start at those coordinates.
#'   * Adds a `censor` column.
#'   * Filters the rows before `max_time`.
#'
#' @inheritParams g_km
#' @param fit_km (`survfit`)\cr result of [survival::survfit()].
#' @param armval (`string`)\cr used as strata name when treatment arm variable only has one level. Default is `"All"`.
#'
#' @return A `tibble` with columns `time`, `n.risk`, `n.event`, `n.censor`, `estimate`, `std.error`, `conf.high`,
#'   `conf.low`, `strata`, and `censor`.
#'
#' @examples
#' \donttest{
#' library(dplyr)
#' library(survival)
#'
#' # Test with multiple arms
#' tern_ex_adtte %>%
#'   filter(PARAMCD == "OS") %>%
#'   survfit(form = Surv(AVAL, 1 - CNSR) ~ ARMCD, data = .) %>%
#'   h_data_plot()
#'
#' # Test with single arm
#' tern_ex_adtte %>%
#'   filter(PARAMCD == "OS", ARMCD == "ARM B") %>%
#'   survfit(form = Surv(AVAL, 1 - CNSR) ~ ARMCD, data = .) %>%
#'   h_data_plot(armval = "ARM B")
#' }
#'
#' @export
h_data_plot <- function(fit_km,
                        armval = "All",
                        max_time = NULL) {
  y <- broom::tidy(fit_km)

  if (!is.null(fit_km$strata)) {
    fit_km_var_level <- strsplit(sub("=", "equals", names(fit_km$strata)), "equals")
    strata_levels <- vapply(fit_km_var_level, FUN = "[", FUN.VALUE = "a", i = 2)
    strata_var_level <- strsplit(sub("=", "equals", y$strata), "equals")
    y$strata <- factor(
      vapply(strata_var_level, FUN = "[", FUN.VALUE = "a", i = 2),
      levels = strata_levels
    )
  } else {
    y$strata <- armval
  }

  y_by_strata <- split(y, y$strata)
  y_by_strata_extended <- lapply(
    y_by_strata,
    FUN = function(tbl) {
      first_row <- tbl[1L, ]
      first_row$time <- 0
      first_row$n.risk <- sum(first_row[, c("n.risk", "n.event", "n.censor")])
      first_row$n.event <- first_row$n.censor <- 0
      first_row$estimate <- first_row$conf.high <- first_row$conf.low <- 1
      first_row$std.error <- 0
      rbind(
        first_row,
        tbl
      )
    }
  )
  y <- do.call(rbind, y_by_strata_extended)

  y$censor <- ifelse(y$n.censor > 0, y$estimate, NA)
  if (!is.null(max_time)) {
    y <- y[y$time <= max(max_time), ]
  }
  y
}

#' Helper function: x tick positions
#'
#' @description `r lifecycle::badge("stable")`
#'
#' Calculate the positions of ticks on the x-axis. However, if `xticks` already
#' exists it is kept as is. It is based on the same function `ggplot2` relies on,
#' and is required in the graphic and the patient-at-risk annotation table.
#'
#' @inheritParams g_km
#' @inheritParams h_ggkm
#'
#' @return A vector of positions to use for x-axis ticks on a `ggplot` object.
#'
#' @examples
#' \donttest{
#' library(dplyr)
#' library(survival)
#'
#' data <- tern_ex_adtte %>%
#'   filter(PARAMCD == "OS") %>%
#'   survfit(form = Surv(AVAL, 1 - CNSR) ~ ARMCD, data = .) %>%
#'   h_data_plot()
#'
#' h_xticks(data)
#' h_xticks(data, xticks = seq(0, 3000, 500))
#' h_xticks(data, xticks = 500)
#' h_xticks(data, xticks = 500, max_time = 6000)
#' h_xticks(data, xticks = c(0, 500), max_time = 300)
#' h_xticks(data, xticks = 500, max_time = 300)
#' }
#'
#' @export
h_xticks <- function(data, xticks = NULL, max_time = NULL) {
  if (is.null(xticks)) {
    if (is.null(max_time)) {
      labeling::extended(range(data$time)[1], range(data$time)[2], m = 5)
    } else {
      labeling::extended(range(data$time)[1], max(range(data$time)[2], max_time), m = 5)
    }
  } else if (checkmate::test_number(xticks)) {
    if (is.null(max_time)) {
      seq(0, max(data$time), xticks)
    } else {
      seq(0, max(data$time, max_time), xticks)
    }
  } else if (is.numeric(xticks)) {
    xticks
  } else {
    stop(
      paste(
        "xticks should be either `NULL`",
        "or a single number (interval between x ticks)",
        "or a numeric vector (position of ticks on the x axis)"
      )
    )
  }
}

#' Helper function: KM plot
#'
#' @description `r lifecycle::badge("stable")`
#'
#' Draw the Kaplan-Meier plot using `ggplot2`.
#'
#' @inheritParams g_km
#' @param data (`data.frame`)\cr survival data as pre-processed by `h_data_plot`.
#'
#' @return A `ggplot` object.
#'
#' @examples
#' \donttest{
#' library(dplyr)
#' library(survival)
#'
#' fit_km <- tern_ex_adtte %>%
#'   filter(PARAMCD == "OS") %>%
#'   survfit(form = Surv(AVAL, 1 - CNSR) ~ ARMCD, data = .)
#' data_plot <- h_data_plot(fit_km = fit_km)
#' xticks <- h_xticks(data = data_plot)
#' gg <- h_ggkm(
#'   data = data_plot,
#'   censor_show = TRUE,
#'   xticks = xticks,
#'   xlab = "Days",
#'   yval = "Survival",
#'   ylab = "Survival Probability",
#'   title = "Survival"
#' )
#' gg
#' }
#'
#' @export
h_ggkm <- function(data,
                   xticks = NULL,
                   yval = "Survival",
                   censor_show,
                   xlab,
                   ylab,
                   ylim = NULL,
                   title,
                   footnotes = NULL,
                   max_time = NULL,
                   lwd = 1,
                   lty = NULL,
                   pch = 3,
                   size = 2,
                   col = NULL,
                   ci_ribbon = FALSE,
                   ggtheme = nestcolor::theme_nest()) {
  checkmate::assert_numeric(lty, null.ok = TRUE)
  checkmate::assert_character(col, null.ok = TRUE)

  if (is.null(ylim)) {
    data_lims <- data
    if (yval == "Failure") data_lims[["estimate"]] <- 1 - data_lims[["estimate"]]
    if (!is.null(max_time)) {
      y_lwr <- min(data_lims[data_lims$time < max_time, ][["estimate"]])
      y_upr <- max(data_lims[data_lims$time < max_time, ][["estimate"]])
    } else {
      y_lwr <- min(data_lims[["estimate"]])
      y_upr <- max(data_lims[["estimate"]])
    }
    ylim <- c(y_lwr, y_upr)
  }
  checkmate::assert_numeric(ylim, finite = TRUE, any.missing = FALSE, len = 2, sorted = TRUE)

  # change estimates of survival to estimates of failure (1 - survival)
  if (yval == "Failure") {
    data$estimate <- 1 - data$estimate
    data[c("conf.high", "conf.low")] <- list(1 - data$conf.low, 1 - data$conf.high)
    data$censor <- 1 - data$censor
  }

  gg <- {
    ggplot2::ggplot(
      data = data,
      mapping = ggplot2::aes(
        x = .data[["time"]],
        y = .data[["estimate"]],
        ymin = .data[["conf.low"]],
        ymax = .data[["conf.high"]],
        color = .data[["strata"]],
        fill = .data[["strata"]]
      )
    ) +
      ggplot2::geom_hline(yintercept = 0)
  }

  if (ci_ribbon) {
    gg <- gg + ggplot2::geom_ribbon(alpha = .3, lty = 0)
  }

  gg <- if (is.null(lty)) {
    gg +
      ggplot2::geom_step(linewidth = lwd)
  } else if (checkmate::test_number(lty)) {
    gg +
      ggplot2::geom_step(linewidth = lwd, lty = lty)
  } else if (is.numeric(lty)) {
    gg +
      ggplot2::geom_step(mapping = ggplot2::aes(linetype = .data[["strata"]]), linewidth = lwd) +
      ggplot2::scale_linetype_manual(values = lty)
  }

  gg <- gg +
    ggplot2::coord_cartesian(ylim = ylim) +
    ggplot2::labs(x = xlab, y = ylab, title = title, caption = footnotes)

  if (!is.null(col)) {
    gg <- gg +
      ggplot2::scale_color_manual(values = col) +
      ggplot2::scale_fill_manual(values = col)
  }
  if (censor_show) {
    dt <- data[data$n.censor != 0, ]
    dt$censor_lbl <- factor("Censored")

    gg <- gg + ggplot2::geom_point(
      data = dt,
      ggplot2::aes(
        x = .data[["time"]],
        y = .data[["censor"]],
        shape = .data[["censor_lbl"]]
      ),
      size = size,
      show.legend = TRUE,
      inherit.aes = TRUE
    ) +
      ggplot2::scale_shape_manual(name = NULL, values = pch) +
      ggplot2::guides(
        shape = ggplot2::guide_legend(override.aes = list(linetype = NA)),
        fill = ggplot2::guide_legend(override.aes = list(shape = NA))
      )
  }

  if (!is.null(max_time) && !is.null(xticks)) {
    gg <- gg + ggplot2::scale_x_continuous(breaks = xticks, limits = c(min(0, xticks), max(c(xticks, max_time))))
  } else if (!is.null(xticks)) {
    if (max(data$time) <= max(xticks)) {
      gg <- gg + ggplot2::scale_x_continuous(breaks = xticks, limits = c(min(0, min(xticks)), max(xticks)))
    } else {
      gg <- gg + ggplot2::scale_x_continuous(breaks = xticks)
    }
  } else if (!is.null(max_time)) {
    gg <- gg + ggplot2::scale_x_continuous(limits = c(0, max_time))
  }

  if (!is.null(ggtheme)) {
    gg <- gg + ggtheme
  }

  gg + ggplot2::theme(
    legend.position = "bottom",
    legend.title = ggplot2::element_blank(),
    legend.key.height = unit(0.02, "npc"),
    panel.grid.major.x = ggplot2::element_line(linewidth = 2)
  )
}

#' `ggplot` Decomposition
#'
#' @description `r lifecycle::badge("stable")`
#'
#' The elements composing the `ggplot` are extracted and organized in a `list`.
#'
#' @param gg (`ggplot`)\cr a graphic to decompose.
#'
#' @return A named `list` with elements:
#'   * `panel`: The panel.
#'   * `yaxis`: The y-axis.
#'   * `xaxis`: The x-axis.
#'   * `xlab`: The x-axis label.
#'   * `ylab`: The y-axis label.
#'   * `guide`: The legend.
#'
#' @examples
#' \donttest{
#' library(dplyr)
#' library(survival)
#' library(grid)
#'
#' fit_km <- tern_ex_adtte %>%
#'   filter(PARAMCD == "OS") %>%
#'   survfit(form = Surv(AVAL, 1 - CNSR) ~ ARMCD, data = .)
#' data_plot <- h_data_plot(fit_km = fit_km)
#' xticks <- h_xticks(data = data_plot)
#' gg <- h_ggkm(
#'   data = data_plot,
#'   yval = "Survival",
#'   censor_show = TRUE,
#'   xticks = xticks, xlab = "Days", ylab = "Survival Probability",
#'   title = "tt",
#'   footnotes = "ff"
#' )
#'
#' g_el <- h_decompose_gg(gg)
#' grid::grid.newpage()
#' grid.rect(gp = grid::gpar(lty = 1, col = "red", fill = "gray85", lwd = 5))
#' grid::grid.draw(g_el$panel)
#'
#' grid::grid.newpage()
#' grid.rect(gp = grid::gpar(lty = 1, col = "royalblue", fill = "gray85", lwd = 5))
#' grid::grid.draw(with(g_el, cbind(ylab, yaxis)))
#' }
#'
#' @export
h_decompose_gg <- function(gg) {
  g_el <- ggplot2::ggplotGrob(gg)
  y <- c(
    panel = "panel",
    yaxis = "axis-l",
    xaxis = "axis-b",
    xlab = "xlab-b",
    ylab = "ylab-l",
    guide = "guide"
  )
  lapply(X = y, function(x) gtable::gtable_filter(g_el, x))
}

#' Helper: KM Layout
#'
#' @description `r lifecycle::badge("stable")`
#'
#' Prepares a (5 rows) x (2 cols) layout for the Kaplan-Meier curve.
#'
#' @inheritParams g_km
#' @inheritParams h_ggkm
#' @param g_el (`list` of `gtable`)\cr list as obtained by `h_decompose_gg()`.
#' @param annot_at_risk (`flag`)\cr compute and add the annotation table reporting the number of
#'   patient at risk matching the main grid of the Kaplan-Meier curve.
#'
#' @return A grid layout.
#'
#' @details The layout corresponds to a grid of two columns and five rows of unequal dimensions. Most of the
#'   dimension are fixed, only the curve is flexible and will accommodate with the remaining free space.
#'   * The left column gets the annotation of the `ggplot` (y-axis) and the names of the strata for the patient
#'     at risk tabulation. The main constraint is about the width of the columns which must allow the writing of
#'     the strata name.
#'   * The right column receive the `ggplot`, the legend, the x-axis and the patient at risk table.
#'
#' @examples
#' \donttest{
#' library(dplyr)
#' library(survival)
#' library(grid)
#'
#' fit_km <- tern_ex_adtte %>%
#'   filter(PARAMCD == "OS") %>%
#'   survfit(form = Surv(AVAL, 1 - CNSR) ~ ARMCD, data = .)
#' data_plot <- h_data_plot(fit_km = fit_km)
#' xticks <- h_xticks(data = data_plot)
#' gg <- h_ggkm(
#'   data = data_plot,
#'   censor_show = TRUE,
#'   xticks = xticks, xlab = "Days", ylab = "Survival Probability",
#'   title = "tt", footnotes = "ff", yval = "Survival"
#' )
#' g_el <- h_decompose_gg(gg)
#' lyt <- h_km_layout(data = data_plot, g_el = g_el, title = "t", footnotes = "f")
#' grid.show.layout(lyt)
#' }
#'
#' @export
h_km_layout <- function(data, g_el, title, footnotes, annot_at_risk = TRUE, annot_at_risk_title = TRUE) {
  txtlines <- levels(as.factor(data$strata))
  nlines <- nlevels(as.factor(data$strata))
  col_annot_width <- max(
    c(
      as.numeric(grid::convertX(g_el$yaxis$width + g_el$ylab$width, "pt")),
      as.numeric(
        grid::convertX(
          grid::stringWidth(txtlines) + grid::unit(7, "pt"), "pt"
        )
      )
    )
  )

  ttl_row <- as.numeric(!is.null(title))
  foot_row <- as.numeric(!is.null(footnotes))
  no_tbl_ind <- c()
  ht_x <- c()
  ht_units <- c()

  if (ttl_row == 1) {
    no_tbl_ind <- c(no_tbl_ind, TRUE)
    ht_x <- c(ht_x, 2)
    ht_units <- c(ht_units, "lines")
  }

  no_tbl_ind <- c(no_tbl_ind, rep(TRUE, 3), rep(FALSE, 2))
  ht_x <- c(
    ht_x,
    1,
    grid::convertX(with(g_el, xaxis$height + ylab$width), "pt") + grid::unit(5, "pt"),
    grid::convertX(g_el$guide$heights, "pt") + grid::unit(2, "pt"),
    1,
    nlines + 0.5,
    grid::convertX(with(g_el, xaxis$height + ylab$width), "pt")
  )
  ht_units <- c(
    ht_units,
    "null",
    "pt",
    "pt",
    "lines",
    "lines",
    "pt"
  )

  if (foot_row == 1) {
    no_tbl_ind <- c(no_tbl_ind, TRUE)
    ht_x <- c(ht_x, 1)
    ht_units <- c(ht_units, "lines")
  }
  if (annot_at_risk) {
    no_at_risk_tbl <- rep(TRUE, 6 + ttl_row + foot_row)
    if (!annot_at_risk_title) {
      no_at_risk_tbl[length(no_at_risk_tbl) - 2 - foot_row] <- FALSE
    }
  } else {
    no_at_risk_tbl <- no_tbl_ind
  }

  grid::grid.layout(
    nrow = sum(no_at_risk_tbl), ncol = 2,
    widths = grid::unit(c(col_annot_width, 1), c("pt", "null")),
    heights = grid::unit(
      x = ht_x[no_at_risk_tbl],
      units = ht_units[no_at_risk_tbl]
    )
  )
}

#' Helper: Patient-at-Risk Grobs
#'
#' @description `r lifecycle::badge("stable")`
#'
#' Two graphical objects are obtained, one corresponding to row labeling and the second to the table of
#' numbers of patients at risk. If `title = TRUE`, a third object corresponding to the table title is
#' also obtained.
#'
#' @inheritParams g_km
#' @inheritParams h_ggkm
#' @param annot_tbl (`data.frame`)\cr annotation as prepared by [survival::summary.survfit()] which
#'   includes the number of patients at risk at given time points.
#' @param xlim (`numeric`)\cr the maximum value on the x-axis (used to
#'   ensure the at risk table aligns with the KM graph).
#' @param title (`flag`)\cr whether the "Patients at Risk" title should be added above the `annot_at_risk`
#'   table. Has no effect if `annot_at_risk` is `FALSE`. Defaults to `TRUE`.
#'
#' @return A named `list` of two `gTree` objects if `title = FALSE`: `at_risk` and `label`, or three
#'   `gTree` objects if `title = TRUE`: `at_risk`, `label`, and `title`.
#'
#' @examples
#' \donttest{
#' library(dplyr)
#' library(survival)
#' library(grid)
#'
#' fit_km <- tern_ex_adtte %>%
#'   filter(PARAMCD == "OS") %>%
#'   survfit(form = Surv(AVAL, 1 - CNSR) ~ ARMCD, data = .)
#'
#' data_plot <- h_data_plot(fit_km = fit_km)
#'
#' xticks <- h_xticks(data = data_plot)
#'
#' gg <- h_ggkm(
#'   data = data_plot,
#'   censor_show = TRUE,
#'   xticks = xticks, xlab = "Days", ylab = "Survival Probability",
#'   title = "tt", footnotes = "ff", yval = "Survival"
#' )
#'
#' # The annotation table reports the patient at risk for a given strata and
#' # time (`xticks`).
#' annot_tbl <- summary(fit_km, time = xticks)
#' if (is.null(fit_km$strata)) {
#'   annot_tbl <- with(annot_tbl, data.frame(n.risk = n.risk, time = time, strata = "All"))
#' } else {
#'   strata_lst <- strsplit(sub("=", "equals", levels(annot_tbl$strata)), "equals")
#'   levels(annot_tbl$strata) <- matrix(unlist(strata_lst), ncol = 2, byrow = TRUE)[, 2]
#'   annot_tbl <- data.frame(
#'     n.risk = annot_tbl$n.risk,
#'     time = annot_tbl$time,
#'     strata = annot_tbl$strata
#'   )
#' }
#'
#' # The annotation table is transformed into a grob.
#' tbl <- h_grob_tbl_at_risk(data = data_plot, annot_tbl = annot_tbl, xlim = max(xticks))
#'
#' # For the representation, the layout is estimated for which the decomposition
#' # of the graphic element is necessary.
#' g_el <- h_decompose_gg(gg)
#' lyt <- h_km_layout(data = data_plot, g_el = g_el, title = "t", footnotes = "f")
#'
#' grid::grid.newpage()
#' pushViewport(viewport(layout = lyt, height = .95, width = .95))
#' grid.rect(gp = grid::gpar(lty = 1, col = "purple", fill = "gray85", lwd = 1))
#' pushViewport(viewport(layout.pos.row = 3:4, layout.pos.col = 2))
#' grid.rect(gp = grid::gpar(lty = 1, col = "orange", fill = "gray85", lwd = 1))
#' grid::grid.draw(tbl$at_risk)
#' popViewport()
#' pushViewport(viewport(layout.pos.row = 3:4, layout.pos.col = 1))
#' grid.rect(gp = grid::gpar(lty = 1, col = "green3", fill = "gray85", lwd = 1))
#' grid::grid.draw(tbl$label)
#' }
#'
#' @export
h_grob_tbl_at_risk <- function(data, annot_tbl, xlim, title = TRUE) {
  txtlines <- levels(as.factor(data$strata))
  nlines <- nlevels(as.factor(data$strata))
  y_int <- annot_tbl$time[2] - annot_tbl$time[1]
  annot_tbl <- expand.grid(
    time = seq(0, xlim, y_int),
    strata = unique(annot_tbl$strata)
  ) %>% dplyr::left_join(annot_tbl, by = c("time", "strata"))
  annot_tbl[is.na(annot_tbl)] <- 0
  y_str_unit <- as.numeric(annot_tbl$strata)
  vp_table <- grid::plotViewport(margins = grid::unit(c(0, 0, 0, 0), "lines"))
  if (title) {
    gb_table_title <- grid::gList(
      grid::textGrob(
        label = "Patients at Risk:",
        x = 1,
        y = grid::unit(0.2, "native"),
        gp = grid::gpar(fontface = "bold", fontsize = 10)
      )
    )
  }
  gb_table_left_annot <- grid::gList(
    grid::rectGrob(
      x = 0, y = grid::unit(c(1:nlines) - 1, "lines"),
      gp = grid::gpar(fill = c("gray95", "gray90"), alpha = 1, col = "white"),
      height = grid::unit(1, "lines"), just = "bottom", hjust = 0
    ),
    grid::textGrob(
      label = unique(annot_tbl$strata),
      x = 0.5,
      y = grid::unit(
        (max(unique(y_str_unit)) - unique(y_str_unit)) + 0.75,
        "native"
      ),
      gp = grid::gpar(fontface = "italic", fontsize = 10)
    )
  )
  gb_patient_at_risk <- grid::gList(
    grid::rectGrob(
      x = 0, y = grid::unit(c(1:nlines) - 1, "lines"),
      gp = grid::gpar(fill = c("gray95", "gray90"), alpha = 1, col = "white"),
      height = grid::unit(1, "lines"), just = "bottom", hjust = 0
    ),
    grid::textGrob(
      label = annot_tbl$n.risk,
      x = grid::unit(annot_tbl$time, "native"),
      y = grid::unit(
        (max(y_str_unit) - y_str_unit) + .5,
        "line"
      ) # maybe native
    )
  )

  ret <- list(
    at_risk = grid::gList(
      grid::gTree(
        vp = vp_table,
        children = grid::gList(
          grid::gTree(
            vp = grid::dataViewport(
              xscale = c(0, xlim) + c(-0.05, 0.05) * xlim,
              yscale = c(0, nlines + 1),
              extension = c(0.05, 0)
            ),
            children = grid::gList(gb_patient_at_risk)
          )
        )
      )
    ),
    label = grid::gList(
      grid::gTree(
        vp = grid::viewport(width = max(grid::stringWidth(txtlines))),
        children = grid::gList(
          grid::gTree(
            vp = grid::dataViewport(
              xscale = 0:1,
              yscale = c(0, nlines + 1),
              extension = c(0.0, 0)
            ),
            children = grid::gList(gb_table_left_annot)
          )
        )
      )
    )
  )

  if (title) {
    ret[["title"]] <- grid::gList(
      grid::gTree(
        vp = grid::viewport(width = max(grid::stringWidth(txtlines))),
        children = grid::gList(
          grid::gTree(
            vp = grid::dataViewport(
              xscale = 0:1,
              yscale = c(0, 1),
              extension = c(0, 0)
            ),
            children = grid::gList(gb_table_title)
          )
        )
      )
    )
  }

  ret
}

#' Helper Function: Survival Estimations
#'
#' @description `r lifecycle::badge("stable")`
#'
#' Transform a survival fit to a table with groups in rows characterized by N, median and confidence interval.
#'
#' @inheritParams h_data_plot
#'
#' @return A summary table with statistics `N`, `Median`, and `XX% CI` (`XX` taken from `fit_km`).
#'
#' @examples
#' \donttest{
#' library(dplyr)
#' library(survival)
#'
#' adtte <- tern_ex_adtte %>% filter(PARAMCD == "OS")
#' fit <- survfit(
#'   form = Surv(AVAL, 1 - CNSR) ~ ARMCD,
#'   data = adtte
#' )
#' h_tbl_median_surv(fit_km = fit)
#' }
#'
#' @export
h_tbl_median_surv <- function(fit_km, armval = "All") {
  y <- if (is.null(fit_km$strata)) {
    as.data.frame(t(summary(fit_km)$table), row.names = armval)
  } else {
    tbl <- summary(fit_km)$table
    rownames_lst <- strsplit(sub("=", "equals", rownames(tbl)), "equals")
    rownames(tbl) <- matrix(unlist(rownames_lst), ncol = 2, byrow = TRUE)[, 2]
    as.data.frame(tbl)
  }
  conf.int <- summary(fit_km)$conf.int # nolint
  y$records <- round(y$records)
  y$median <- signif(y$median, 4)
  y$`CI` <- paste0(
    "(", signif(y[[paste0(conf.int, "LCL")]], 4), ", ", signif(y[[paste0(conf.int, "UCL")]], 4), ")"
  )
  stats::setNames(
    y[c("records", "median", "CI")],
    c("N", "Median", f_conf_level(conf.int))
  )
}

#' Helper Function: Survival Estimation Grob
#'
#' @description `r lifecycle::badge("stable")`
#'
#' The survival fit is transformed in a grob containing a table with groups in
#' rows characterized by N, median and 95% confidence interval.
#'
#' @inheritParams g_km
#' @inheritParams h_data_plot
#' @param ttheme (`list`)\cr see [gridExtra::ttheme_default()].
#' @param x (`numeric`)\cr a value between 0 and 1 specifying x-location.
#' @param y (`numeric`)\cr a value between 0 and 1 specifying y-location.
#' @param width (`unit`)\cr width (as a unit) to use when printing the grob.
#'
#' @return A `grob` of a table containing statistics `N`, `Median`, and `XX% CI` (`XX` taken from `fit_km`).
#'
#' @examples
#' \donttest{
#' library(dplyr)
#' library(survival)
#' library(grid)
#'
#' grid::grid.newpage()
#' grid.rect(gp = grid::gpar(lty = 1, col = "pink", fill = "gray85", lwd = 1))
#' tern_ex_adtte %>%
#'   filter(PARAMCD == "OS") %>%
#'   survfit(form = Surv(AVAL, 1 - CNSR) ~ ARMCD, data = .) %>%
#'   h_grob_median_surv() %>%
#'   grid::grid.draw()
#' }
#'
#' @export
h_grob_median_surv <- function(fit_km,
                               armval = "All",
                               x = 0.9,
                               y = 0.9,
                               width = grid::unit(0.3, "npc"),
                               ttheme = gridExtra::ttheme_default()) {
  data <- h_tbl_median_surv(fit_km, armval = armval)

  width <- grid::convertUnit(width, "in")
  height <- width * (nrow(data) + 1) / 12

  w <- paste(" ", c(
    rownames(data)[which.max(nchar(rownames(data)))],
    sapply(names(data), function(x) c(x, data[[x]])[which.max(nchar(c(x, data[[x]])))])
  ))
  w_unit <- grid::convertWidth(grid::stringWidth(w), "in", valueOnly = TRUE)

  w_txt <- sapply(1:64, function(x) {
    graphics::par(ps = x)
    graphics::strwidth(w[4], units = "in")
  })
  f_size_w <- which.max(w_txt[w_txt < as.numeric((w_unit / sum(w_unit)) * width)[4]])

  h_txt <- sapply(1:64, function(x) {
    graphics::par(ps = x)
    graphics::strheight(grid::stringHeight("X"), units = "in")
  })
  f_size_h <- which.max(h_txt[h_txt < as.numeric(grid::unit(as.numeric(height) / 4, grid::unitType(height)))])

  if (ttheme$core$fg_params$fontsize == 12) {
    ttheme$core$fg_params$fontsize <- min(f_size_w, f_size_h)
    ttheme$colhead$fg_params$fontsize <- min(f_size_w, f_size_h)
    ttheme$rowhead$fg_params$fontsize <- min(f_size_w, f_size_h)
  }

  gt <- gridExtra::tableGrob(
    d = data,
    theme = ttheme
  )
  gt$widths <- ((w_unit / sum(w_unit)) * width)
  gt$heights <- rep(grid::unit(as.numeric(height) / 4, grid::unitType(height)), nrow(gt))

  vp <- grid::viewport(
    x = grid::unit(x, "npc") + grid::unit(1, "lines"),
    y = grid::unit(y, "npc") + grid::unit(1.5, "lines"),
    height = height,
    width = width,
    just = c("right", "top")
  )

  grid::gList(
    grid::gTree(
      vp = vp,
      children = grid::gList(gt)
    )
  )
}

#' Helper: Grid Object with y-axis Annotation
#'
#' @description `r lifecycle::badge("stable")`
#'
#' Build the y-axis annotation from a decomposed `ggplot`.
#'
#' @param ylab (`gtable`)\cr the y-lab as a graphical object derived from a `ggplot`.
#' @param yaxis (`gtable`)\cr the y-axis as a graphical object derived from a `ggplot`.
#'
#' @return a `gTree` object containing the y-axis annotation from a `ggplot`.
#'
#' @examples
#' \donttest{
#' library(dplyr)
#' library(survival)
#' library(grid)
#'
#' fit_km <- tern_ex_adtte %>%
#'   filter(PARAMCD == "OS") %>%
#'   survfit(form = Surv(AVAL, 1 - CNSR) ~ ARMCD, data = .)
#' data_plot <- h_data_plot(fit_km = fit_km)
#' xticks <- h_xticks(data = data_plot)
#' gg <- h_ggkm(
#'   data = data_plot,
#'   censor_show = TRUE,
#'   xticks = xticks, xlab = "Days", ylab = "Survival Probability",
#'   title = "title", footnotes = "footnotes", yval = "Survival"
#' )
#'
#' g_el <- h_decompose_gg(gg)
#'
#' grid::grid.newpage()
#' pvp <- grid::plotViewport(margins = c(5, 4, 2, 20))
#' pushViewport(pvp)
#' grid::grid.draw(h_grob_y_annot(ylab = g_el$ylab, yaxis = g_el$yaxis))
#' grid.rect(gp = grid::gpar(lty = 1, col = "gray35", fill = NA))
#' }
#'
#' @export
h_grob_y_annot <- function(ylab, yaxis) {
  grid::gList(
    grid::gTree(
      vp = grid::viewport(
        width = grid::convertX(yaxis$width + ylab$width, "pt"),
        x = grid::unit(1, "npc"),
        just = "right"
      ),
      children = grid::gList(cbind(ylab, yaxis))
    )
  )
}

#' Helper Function: Pairwise `CoxPH` table
#'
#' @description `r lifecycle::badge("stable")`
#'
#' Create a `data.frame` of pairwise stratified or unstratified `CoxPH` analysis results.
#'
#' @inheritParams g_km
#'
#' @return A `data.frame` containing statistics `HR`, `XX% CI` (`XX` taken from `control_coxph_pw`),
#'   and `p-value (log-rank)`.
#'
#' @examples
#' \donttest{
#' library(dplyr)
#'
#' adtte <- tern_ex_adtte %>%
#'   filter(PARAMCD == "OS") %>%
#'   mutate(is_event = CNSR == 0)
#'
#' h_tbl_coxph_pairwise(
#'   df = adtte,
#'   variables = list(tte = "AVAL", is_event = "is_event", arm = "ARM"),
#'   control_coxph_pw = control_coxph(conf_level = 0.9)
#' )
#' }
#'
#' @export
h_tbl_coxph_pairwise <- function(df,
                                 variables,
                                 control_coxph_pw = control_coxph()) {
  assert_df_with_variables(df, variables)
  arm <- variables$arm
  df[[arm]] <- factor(df[[arm]])
  ref_group <- levels(df[[arm]])[1]
  comp_group <- levels(df[[arm]])[-1]
  results <- Map(function(comp) {
    res <- s_coxph_pairwise(
      df = df[df[[arm]] == comp, , drop = FALSE],
      .ref_group = df[df[[arm]] == ref_group, , drop = FALSE],
      .in_ref_col = FALSE,
      .var = variables$tte,
      is_event = variables$is_event,
      strat = variables$strat,
      control = control_coxph_pw
    )
    res_df <- data.frame(
      hr = format(round(res$hr, 2), nsmall = 2),
      hr_ci = paste0(
        "(", format(round(res$hr_ci[1], 2), nsmall = 2), ", ",
        format(round(res$hr_ci[2], 2), nsmall = 2), ")"
      ),
      pvalue = if (res$pvalue < 0.0001) "<0.0001" else format(round(res$pvalue, 4), 4),
      stringsAsFactors = FALSE
    )
    colnames(res_df) <- c("HR", vapply(res[c("hr_ci", "pvalue")], obj_label, FUN.VALUE = "character"))
    row.names(res_df) <- comp
    res_df
  }, comp_group)
  do.call(rbind, results)
}

#' Helper Function: `CoxPH` Grob
#'
#' @description `r lifecycle::badge("stable")`
#'
#' Grob of `rtable` output from [h_tbl_coxph_pairwise()]
#'
#' @inheritParams h_grob_median_surv
#' @param ... arguments will be passed to [h_tbl_coxph_pairwise()].
#' @param x (`numeric`)\cr a value between 0 and 1 specifying x-location.
#' @param y (`numeric`)\cr a value between 0 and 1 specifying y-location.
#' @param width (`unit`)\cr width (as a unit) to use when printing the grob.
#'
#' @return A `grob` of a table containing statistics `HR`, `XX% CI` (`XX` taken from `control_coxph_pw`),
#'   and `p-value (log-rank)`.
#'
#' @examples
#' \donttest{
#' library(dplyr)
#' library(survival)
#' library(grid)
#'
#' grid::grid.newpage()
#' grid.rect(gp = grid::gpar(lty = 1, col = "pink", fill = "gray85", lwd = 1))
#' data <- tern_ex_adtte %>%
#'   filter(PARAMCD == "OS") %>%
#'   mutate(is_event = CNSR == 0)
#' tbl_grob <- h_grob_coxph(
#'   df = data,
#'   variables = list(tte = "AVAL", is_event = "is_event", arm = "ARMCD"),
#'   control_coxph_pw = control_coxph(conf_level = 0.9), x = 0.5, y = 0.5
#' )
#' grid::grid.draw(tbl_grob)
#' }
#'
#' @export
h_grob_coxph <- function(...,
                         x = 0,
                         y = 0,
                         width = grid::unit(0.4, "npc"),
                         ttheme = gridExtra::ttheme_default(
                           padding = grid::unit(c(1, .5), "lines"),
                           core = list(bg_params = list(fill = c("grey95", "grey90"), alpha = .5))
                         )) {
  data <- h_tbl_coxph_pairwise(...)

  width <- grid::convertUnit(width, "in")
  height <- width * (nrow(data) + 1) / 12

  w <- paste("    ", c(
    rownames(data)[which.max(nchar(rownames(data)))],
    sapply(names(data), function(x) c(x, data[[x]])[which.max(nchar(c(x, data[[x]])))])
  ))
  w_unit <- grid::convertWidth(grid::stringWidth(w), "in", valueOnly = TRUE)

  w_txt <- sapply(1:64, function(x) {
    graphics::par(ps = x)
    graphics::strwidth(w[4], units = "in")
  })
  f_size_w <- which.max(w_txt[w_txt < as.numeric((w_unit / sum(w_unit)) * width)[4]])

  h_txt <- sapply(1:64, function(x) {
    graphics::par(ps = x)
    graphics::strheight(grid::stringHeight("X"), units = "in")
  })
  f_size_h <- which.max(h_txt[h_txt < as.numeric(grid::unit(as.numeric(height) / 4, grid::unitType(height)))])

  if (ttheme$core$fg_params$fontsize == 12) {
    ttheme$core$fg_params$fontsize <- min(f_size_w, f_size_h)
    ttheme$colhead$fg_params$fontsize <- min(f_size_w, f_size_h)
    ttheme$rowhead$fg_params$fontsize <- min(f_size_w, f_size_h)
  }

  tryCatch(
    expr = {
      gt <- gridExtra::tableGrob(
        d = data,
        theme = ttheme
      ) # ERROR 'data' must be of a vector type, was 'NULL'
      gt$widths <- ((w_unit / sum(w_unit)) * width)
      gt$heights <- rep(grid::unit(as.numeric(height) / 4, grid::unitType(height)), nrow(gt))
      vp <- grid::viewport(
        x = grid::unit(x, "npc") + grid::unit(1, "lines"),
        y = grid::unit(y, "npc") + grid::unit(1.5, "lines"),
        height = height,
        width = width,
        just = c("left", "bottom")
      )
      grid::gList(
        grid::gTree(
          vp = vp,
          children = grid::gList(gt)
        )
      )
    },
    error = function(w) {
      message(paste(
        "Warning: Cox table will not be displayed as there is",
        "not any level to be compared in the arm variable."
      ))
      return(
        grid::gList(
          grid::gTree(
            vp = NULL,
            children = NULL
          )
        )
      )
    }
  )
}

#' Add Titles, Footnotes, Page Number, and a Bounding Box to a Grid Grob
#'
#' @description `r lifecycle::badge("stable")`
#'
#' This function is useful to label grid grobs (also `ggplot2`, and `lattice` plots)
#' with title, footnote, and page numbers.
#'
#' @inheritParams grid::grob
#' @param grob a grid grob object, optionally `NULL` if only a `grob` with the decoration should be shown.
#' @param titles vector of character strings. Vector elements are separated by a newline and strings are wrapped
#'   according to the page width.
#' @param footnotes vector of character string. Same rules as for `titles`.
#' @param page string with page numeration, if `NULL` then no page number is displayed.
#' @param width_titles unit object
#' @param width_footnotes unit object
#' @param border boolean, whether a a border should be drawn around the plot or not.
#' @param margins unit object of length 4
#' @param padding  unit object of length 4
#' @param outer_margins  unit object of length 4
#' @param gp_titles a `gpar` object
#' @param gp_footnotes a `gpar` object
#'
#' @return A grid grob (`gTree`).
#'
#' @details The titles and footnotes will be ragged, i.e. each title will be wrapped individually.
#'
#' @examples
#' library(grid)
#'
#' titles <- c(
#'   "Edgar Anderson's Iris Data",
#'   paste(
#'     "This famous (Fisher's or Anderson's) iris data set gives the measurements",
#'     "in centimeters of the variables sepal length and width and petal length",
#'     "and width, respectively, for 50 flowers from each of 3 species of iris."
#'   )
#' )
#'
#' footnotes <- c(
#'   "The species are Iris setosa, versicolor, and virginica.",
#'   paste(
#'     "iris is a data frame with 150 cases (rows) and 5 variables (columns) named",
#'     "Sepal.Length, Sepal.Width, Petal.Length, Petal.Width, and Species."
#'   )
#' )
#'
#' ## empty plot
#' grid.newpage()
#'
#' grid.draw(
#'   decorate_grob(
#'     NULL,
#'     titles = titles,
#'     footnotes = footnotes,
#'     page = "Page 4 of 10"
#'   )
#' )
#'
#' # grid
#' p <- gTree(
#'   children = gList(
#'     rectGrob(),
#'     xaxisGrob(),
#'     yaxisGrob(),
#'     textGrob("Sepal.Length", y = unit(-4, "lines")),
#'     textGrob("Petal.Length", x = unit(-3.5, "lines"), rot = 90),
#'     pointsGrob(iris$Sepal.Length, iris$Petal.Length, gp = gpar(col = iris$Species), pch = 16)
#'   ),
#'   vp = vpStack(plotViewport(), dataViewport(xData = iris$Sepal.Length, yData = iris$Petal.Length))
#' )
#' grid.newpage()
#' grid.draw(p)
#'
#' grid.newpage()
#' grid.draw(
#'   decorate_grob(
#'     grob = p,
#'     titles = titles,
#'     footnotes = footnotes,
#'     page = "Page 6 of 129"
#'   )
#' )
#'
#' ## with ggplot2
#' library(ggplot2)
#'
#' p_gg <- ggplot2::ggplot(iris, aes(Sepal.Length, Sepal.Width, col = Species)) +
#'   ggplot2::geom_point()
#' p_gg
#' p <- ggplotGrob(p_gg)
#' grid.newpage()
#' grid.draw(
#'   decorate_grob(
#'     grob = p,
#'     titles = titles,
#'     footnotes = footnotes,
#'     page = "Page 6 of 129"
#'   )
#' )
#'
#' ## with lattice
#' library(lattice)
#'
#' xyplot(Sepal.Length ~ Petal.Length, data = iris, col = iris$Species)
#' p <- grid.grab()
#' grid.newpage()
#' grid.draw(
#'   decorate_grob(
#'     grob = p,
#'     titles = titles,
#'     footnotes = footnotes,
#'     page = "Page 6 of 129"
#'   )
#' )
#'
#' # with gridExtra - no borders
#' library(gridExtra)
#' grid.newpage()
#' grid.draw(
#'   decorate_grob(
#'     tableGrob(
#'       head(mtcars)
#'     ),
#'     titles = "title",
#'     footnotes = "footnote",
#'     border = FALSE
#'   )
#' )
#'
#' @export
decorate_grob <- function(grob,
                          titles,
                          footnotes,
                          page = "",
                          width_titles = grid::unit(1, "npc") - grid::stringWidth(page),
                          width_footnotes = grid::unit(1, "npc") - grid::stringWidth(page),
                          border = TRUE,
                          margins = grid::unit(c(1, 0, 1, 0), "lines"),
                          padding = grid::unit(rep(1, 4), "lines"),
                          outer_margins = grid::unit(c(2, 1.5, 3, 1.5), "cm"),
                          gp_titles = grid::gpar(),
                          gp_footnotes = grid::gpar(fontsize = 8),
                          name = NULL,
                          gp = grid::gpar(),
                          vp = NULL) {
  st_titles <- split_text_grob(
    titles,
    x = 0, y = 1,
    just = c("left", "top"),
    width = width_titles,
    vp = grid::viewport(layout.pos.row = 1, layout.pos.col = 1),
    gp = gp_titles
  )

  st_footnotes <- split_text_grob(
    footnotes,
    x = 0, y = 1,
    just = c("left", "top"),
    width = width_footnotes,
    vp = grid::viewport(layout.pos.row = 3, layout.pos.col = 1),
    gp = gp_footnotes
  )

  grid::gTree(
    grob = grob,
    titles = titles,
    footnotes = footnotes,
    page = page,
    width_titles = width_titles,
    width_footnotes = width_footnotes,
    border = border,
    margins = margins,
    padding = padding,
    outer_margins = outer_margins,
    gp_titles = gp_titles,
    gp_footnotes = gp_footnotes,
    children = grid::gList(
      grid::gTree(
        children = grid::gList(
          st_titles,
          grid::gTree(
            children = grid::gList(
              if (border) grid::rectGrob(),
              grid::gTree(
                children = grid::gList(
                  grob
                ),
                vp = grid::plotViewport(margins = padding)
              )
            ),
            vp = grid::vpStack(
              grid::viewport(layout.pos.row = 2, layout.pos.col = 1),
              grid::plotViewport(margins = margins)
            )
          ),
          st_footnotes,
          grid::textGrob(
            page,
            x = 1, y = 0,
            just = c("right", "bottom"),
            vp = grid::viewport(layout.pos.row = 3, layout.pos.col = 1),
            gp = gp_footnotes
          )
        ),
        childrenvp = NULL,
        name = "titles_grob_footnotes",
        vp = grid::vpStack(
          grid::plotViewport(margins = outer_margins),
          grid::viewport(
            layout = grid::grid.layout(
              nrow = 3, ncol = 1,
              heights = grid::unit.c(
                grid::grobHeight(st_titles),
                grid::unit(1, "null"),
                grid::grobHeight(st_footnotes)
              )
            )
          )
        )
      )
    ),
    name = name,
    gp = gp,
    vp = vp,
    cl = "decoratedGrob"
  )
}

#' @importFrom grid validDetails
#' @noRd
validDetails.decoratedGrob <- function(x) {
  checkmate::assert_character(x$titles)
  checkmate::assert_character(x$footnotes)

  if (!is.null(x$grob)) {
    checkmate::assert_true(grid::is.grob(x$grob))
  }
  if (length(x$page) == 1) {
    checkmate::assert_character(x$page)
  }
  if (!grid::is.unit(x$outer_margins)) {
    checkmate::assert_vector(x$outer_margins, len = 4)
  }
  if (!grid::is.unit(x$margins)) {
    checkmate::assert_vector(x$margins, len = 4)
  }
  if (!grid::is.unit(x$padding)) {
    checkmate::assert_vector(x$padding, len = 4)
  }

  x
}

#' @importFrom grid widthDetails
#' @noRd
widthDetails.decoratedGrob <- function(x) {
  grid::unit(1, "null")
}

#' @importFrom grid heightDetails
#' @noRd
heightDetails.decoratedGrob <- function(x) {
  grid::unit(1, "null")
}

# Adapted from Paul Murell R Graphics 2nd Edition
# https://www.stat.auckland.ac.nz/~paul/RG2e/interactgrid-splittext.R
split_string <- function(text, width) {
  strings <- strsplit(text, " ")
  out_string <- NA
  for (string_i in seq_along(strings)) {
    newline_str <- strings[[string_i]]
    if (length(newline_str) == 0) newline_str <- ""
    if (is.na(out_string[string_i])) {
      out_string[string_i] <- newline_str[[1]][[1]]
      linewidth <- grid::stringWidth(out_string[string_i])
    }
    gapwidth <- grid::stringWidth(" ")
    availwidth <- as.numeric(width)
    if (length(newline_str) > 1) {
      for (i in seq(2, length(newline_str))) {
        width_i <- grid::stringWidth(newline_str[i])
        if (grid::convertWidth(linewidth + gapwidth + width_i, grid::unitType(width), valueOnly = TRUE) < availwidth) {
          sep <- " "
          linewidth <- linewidth + gapwidth + width_i
        } else {
          sep <- "\n"
          linewidth <- width_i
        }
        out_string[string_i] <- paste(out_string[string_i], newline_str[i], sep = sep)
      }
    }
  }
  paste(out_string, collapse = "\n")
}

#' Split Text According To Available Text Width
#'
#' Dynamically wrap text.
#'
#' @inheritParams grid::grid.text
#' @param text character string
#' @param width a unit object specifying max width of text
#'
#' @return A text grob.
#'
#' @details This code is taken from `R Graphics by Paul Murell, 2nd edition`
#'
#' @keywords internal
split_text_grob <- function(text,
                            x = grid::unit(0.5, "npc"),
                            y = grid::unit(0.5, "npc"),
                            width = grid::unit(1, "npc"),
                            just = "centre",
                            hjust = NULL,
                            vjust = NULL,
                            default.units = "npc", # nolint
                            name = NULL,
                            gp = grid::gpar(),
                            vp = NULL) {
  if (!grid::is.unit(x)) x <- grid::unit(x, default.units)
  if (!grid::is.unit(y)) y <- grid::unit(y, default.units)
  if (!grid::is.unit(width)) width <- grid::unit(width, default.units)
  if (grid::unitType(x) %in% c("sum", "min", "max")) x <- grid::convertUnit(x, default.units)
  if (grid::unitType(y) %in% c("sum", "min", "max")) y <- grid::convertUnit(y, default.units)
  if (grid::unitType(width) %in% c("sum", "min", "max")) width <- grid::convertUnit(width, default.units)

  ## if it is a fixed unit then we do not need to recalculate when viewport resized
  if (!inherits(width, "unit.arithmetic") &&
    !is.null(attr(width, "unit")) &&
    attr(width, "unit") %in% c("cm", "inches", "mm", "points", "picas", "bigpts", "dida", "cicero", "scaledpts")) {
    attr(text, "fixed_text") <- paste(vapply(text, split_string, character(1), width = width), collapse = "\n")
  }

  grid::grid.text(
    label = split_string(text, width),
    x = x, y = y,
    just = just,
    hjust = hjust,
    vjust = vjust,
    rot = 0,
    check.overlap = FALSE,
    name = name,
    gp = gp,
    vp = vp,
    draw = FALSE
  )
}

#' @importFrom grid validDetails
#' @noRd
validDetails.dynamicSplitText <- function(x) {
  checkmate::assert_character(x$text)
  checkmate::assert_true(grid::is.unit(x$width))
  checkmate::assert_vector(x$width, len = 1)
  x
}

#' @importFrom grid heightDetails
#' @noRd
heightDetails.dynamicSplitText <- function(x) {
  txt <- if (!is.null(attr(x$text, "fixed_text"))) {
    attr(x$text, "fixed_text")
  } else {
    paste(vapply(x$text, split_string, character(1), width = x$width), collapse = "\n")
  }
  grid::stringHeight(txt)
}

#' @importFrom grid widthDetails
#' @noRd
widthDetails.dynamicSplitText <- function(x) {
  x$width
}

#' @importFrom grid drawDetails
#' @noRd
drawDetails.dynamicSplitText <- function(x, recording) {
  txt <- if (!is.null(attr(x$text, "fixed_text"))) {
    attr(x$text, "fixed_text")
  } else {
    paste(vapply(x$text, split_string, character(1), width = x$width), collapse = "\n")
  }

  x$width <- NULL
  x$label <- txt
  x$text <- NULL
  class(x) <- c("text", class(x)[-1])

  grid::grid.draw(x)
}

#' Update Page Number
#'
#' Automatically updates page number.
#'
#' @param npages number of pages in total
#' @param ... passed on to [decorate_grob()]
#'
#' @return Closure that increments the page number.
#'
#' @keywords internal
decorate_grob_factory <- function(npages, ...) {
  current_page <- 0
  function(grob) {
    current_page <<- current_page + 1
    if (current_page > npages) {
      stop(paste("current page is", current_page, "but max.", npages, "specified."))
    }
    decorate_grob(grob = grob, page = paste("Page", current_page, "of", npages), ...)
  }
}

#' Decorate Set of `grobs` and Add Page Numbering
#'
#' @description `r lifecycle::badge("stable")`
#'
#' Note that this uses the [decorate_grob_factory()] function.
#'
#' @param grobs a list of grid grobs
#' @param ... arguments passed on to [decorate_grob()].
#'
#' @return A decorated grob.
#'
#' @examples
#' library(ggplot2)
#' library(grid)
#' g <- with(data = iris, {
#'   list(
#'     ggplot2::ggplotGrob(
#'       ggplot2::ggplot(mapping = aes(Sepal.Length, Sepal.Width, col = Species)) +
#'         ggplot2::geom_point()
#'     ),
#'     ggplot2::ggplotGrob(
#'       ggplot2::ggplot(mapping = aes(Sepal.Length, Petal.Length, col = Species)) +
#'         ggplot2::geom_point()
#'     ),
#'     ggplot2::ggplotGrob(
#'       ggplot2::ggplot(mapping = aes(Sepal.Length, Petal.Width, col = Species)) +
#'         ggplot2::geom_point()
#'     ),
#'     ggplot2::ggplotGrob(
#'       ggplot2::ggplot(mapping = aes(Sepal.Width, Petal.Length, col = Species)) +
#'         ggplot2::geom_point()
#'     ),
#'     ggplot2::ggplotGrob(
#'       ggplot2::ggplot(mapping = aes(Sepal.Width, Petal.Width, col = Species)) +
#'         ggplot2::geom_point()
#'     ),
#'     ggplot2::ggplotGrob(
#'       ggplot2::ggplot(mapping = aes(Petal.Length, Petal.Width, col = Species)) +
#'         ggplot2::geom_point()
#'     )
#'   )
#' })
#' lg <- decorate_grob_set(grobs = g, titles = "Hello\nOne\nTwo\nThree", footnotes = "")
#'
#' draw_grob(lg[[1]])
#' draw_grob(lg[[2]])
#' draw_grob(lg[[6]])
#'
#' @export
decorate_grob_set <- function(grobs, ...) {
  n <- length(grobs)
  lgf <- decorate_grob_factory(npages = n, ...)
  lapply(grobs, lgf)
}

#' Confidence Intervals for a Difference of Binomials
#'
#' @description `r lifecycle::badge("experimental")`
#'
#' Several confidence intervals for the difference between proportions.
#'
#' @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.
#'
#' @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,


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	22x	checkmate::assert_character(omit_columns, null.ok = TRUE, min.len = 1, any.missing = FALSE)
76	21x	checkmate::assert_data_frame(data)
77	20x	checkmate::assert_flag(char_as_factor)
78	19x	checkmate::assert_flag(logical_as_factor)
79	19x	checkmate::assert_string(na_level)
80
81	17x	target_vars <- if (is.null(omit_columns)) {
82	15x	names(data)
83		} else {
84	2x	setdiff(names(data), omit_columns) # May have duplicates.
85		}
86	17x	if (length(target_vars) == 0) {
87	1x	return(data)
88		}
89
90	16x	l_target_vars <- split(target_vars, target_vars)
91
92		# Makes sure target_vars exist in data and names are not duplicated.
93	16x	assert_df_with_variables(data, l_target_vars)
94
95	16x	for (x in target_vars) {
96	304x	xi <- data[[x]]
97	304x	xi_label <- obj_label(xi)
98
99		# Determine whether to convert character or logical input.
100	304x	do_char_conversion <- is.character(xi) && char_as_factor
101	304x	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	304x	if (do_logical_conversion) {
106	2x	xi <- as.character(xi)
107		}
108
109	304x	if (is.factor(xi) \|\| is.character(xi)) {
110		# Handle empty strings and NA values.
111	217x	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	217x	if (do_char_conversion \|\| do_logical_conversion) {
116	78x	levels_xi <- setdiff(sort(unique(xi)), na_level)
117	78x	if (na_level %in% unique(xi)) {
118	18x	levels_xi <- c(levels_xi, na_level)
119		}
120
121	78x	xi <- factor(xi, levels = levels_xi)
122		}
123
124	217x	data[, x] <- formatters::with_label(xi, label = xi_label)
125		}
126		}
127	16x	return(data)
128		}

1		#' Odds Ratio Estimation
2		#'
3		#' @description `r lifecycle::badge("stable")`
4		#'
5		#' Compares bivariate responses between two groups in terms of odds ratios
6		#' along with a confidence interval.
7		#'
8		#' @inheritParams argument_convention
9		#'
10		#' @details This function uses either logistic regression for unstratified
11		#' analyses, or conditional logistic regression for stratified analyses.
12		#' The Wald confidence interval with the specified confidence level is
13		#' calculated.
14		#'
15		#' @note For stratified analyses, there is currently no implementation for conditional
16		#' likelihood confidence intervals, therefore the likelihood confidence interval is not
17		#' yet available as an option. Besides, when `rsp` contains only responders or non-responders,
18		#' then the result values will be `NA`, because no odds ratio estimation is possible.
19		#'
20		#' @seealso Relevant helper function [h_odds_ratio()].
21		#'
22		#' @name odds_ratio
23		NULL
24
25		#' @describeIn odds_ratio Statistics function which estimates the odds ratio
26		#' between a treatment and a control. A `variables` list with `arm` and `strata`
27		#' variable names must be passed if a stratified analysis is required.
28		#'
29		#' @inheritParams split_cols_by_groups
30		#'
31		#' @return
32		#' * `s_odds_ratio()` returns a named list with the statistics `or_ci`
33		#' (containing `est`, `lcl`, and `ucl`) and `n_tot`.
34		#'
35		#' @examples
36		#' set.seed(12)
37		#' dta <- data.frame(
38		#' rsp = sample(c(TRUE, FALSE), 100, TRUE),
39		#' grp = factor(rep(c("A", "B"), each = 50), levels = c("B", "A")),
40		#' strata = factor(sample(c("C", "D"), 100, TRUE))
41		#' )
42		#'
43		#' # Unstratified analysis.
44		#' s_odds_ratio(
45		#' df = subset(dta, grp == "A"),
46		#' .var = "rsp",
47		#' .ref_group = subset(dta, grp == "B"),
48		#' .in_ref_col = FALSE,
49		#' .df_row = dta
50		#' )
51		#'
52		#' # Stratified analysis.
53		#' s_odds_ratio(
54		#' df = subset(dta, grp == "A"),
55		#' .var = "rsp",
56		#' .ref_group = subset(dta, grp == "B"),
57		#' .in_ref_col = FALSE,
58		#' .df_row = dta,
59		#' variables = list(arm = "grp", strata = "strata")
60		#' )
61		#'
62		#' @export
63		s_odds_ratio <- function(df,
64		.var,
65		.ref_group,
66		.in_ref_col,
67		.df_row,
68		variables = list(arm = NULL, strata = NULL),
69		conf_level = 0.95,
70		groups_list = NULL) {
71	65x	y <- list(or_ci = "", n_tot = "")
72
73	65x	if (!.in_ref_col) {
74	65x	assert_proportion_value(conf_level)
75	65x	assert_df_with_variables(df, list(rsp = .var))
76	65x	assert_df_with_variables(.ref_group, list(rsp = .var))
77
78	65x	if (is.null(variables$strata)) {
79	52x	data <- data.frame(
80	52x	rsp = c(.ref_group[[.var]], df[[.var]]),
81	52x	grp = factor(
82	52x	rep(c("ref", "Not-ref"), c(nrow(.ref_group), nrow(df))),
83	52x	levels = c("ref", "Not-ref")
84		)
85		)
86	52x	y <- or_glm(data, conf_level = conf_level)
87		} else {
88	13x	assert_df_with_variables(.df_row, c(list(rsp = .var), variables))
89
90		# The group variable prepared for clogit must be synchronised with combination groups definition.
91	13x	if (is.null(groups_list)) {
92	12x	ref_grp <- as.character(unique(.ref_group[[variables$arm]]))
93	12x	trt_grp <- as.character(unique(df[[variables$arm]]))
94	12x	grp <- stats::relevel(factor(.df_row[[variables$arm]]), ref = ref_grp)
95		} else {
96		# If more than one level in reference col.
97	1x	reference <- as.character(unique(.ref_group[[variables$arm]]))
98	1x	grp_ref_flag <- vapply(
99	1x	X = groups_list,
100	1x	FUN.VALUE = TRUE,
101	1x	FUN = function(x) all(reference %in% x)
102		)
103	1x	ref_grp <- names(groups_list)[grp_ref_flag]
104
105		# If more than one level in treatment col.
106	1x	treatment <- as.character(unique(df[[variables$arm]]))
107	1x	grp_trt_flag <- vapply(
108	1x	X = groups_list,
109	1x	FUN.VALUE = TRUE,
110	1x	FUN = function(x) all(treatment %in% x)
111		)
112	1x	trt_grp <- names(groups_list)[grp_trt_flag]
113
114	1x	grp <- combine_levels(.df_row[[variables$arm]], levels = reference, new_level = ref_grp)
115	1x	grp <- combine_levels(grp, levels = treatment, new_level = trt_grp)
116		}
117
118		# The reference level in `grp` must be the same as in the `rtables` column split.
119	13x	data <- data.frame(
120	13x	rsp = .df_row[[.var]],
121	13x	grp = grp,
122	13x	strata = interaction(.df_row[variables$strata])
123		)
124	13x	y_all <- or_clogit(data, conf_level = conf_level)
125	13x	checkmate::assert_string(trt_grp)
126	13x	checkmate::assert_subset(trt_grp, names(y_all$or_ci))
127	12x	y$or_ci <- y_all$or_ci[[trt_grp]]
128	12x	y$n_tot <- y_all$n_tot
129		}
130		}
131
132	64x	y$or_ci <- formatters::with_label(
133	64x	x = y$or_ci,
134	64x	label = paste0("Odds Ratio (", 100 * conf_level, "% CI)")
135		)
136
137	64x	y$n_tot <- formatters::with_label(
138	64x	x = y$n_tot,
139	64x	label = "Total n"
140		)
141
142	64x	y
143		}
144
145		#' @describeIn odds_ratio Formatted analysis function which is used as `afun` in `estimate_odds_ratio()`.
146		#'
147		#' @return
148		#' * `a_odds_ratio()` returns the corresponding list with formatted [rtables::CellValue()].
149		#'
150		#' @examples
151		#' a_odds_ratio(
152		#' df = subset(dta, grp == "A"),
153		#' .var = "rsp",
154		#' .ref_group = subset(dta, grp == "B"),
155		#' .in_ref_col = FALSE,
156		#' .df_row = dta
157		#' )
158		#'
159		#' @export
160		a_odds_ratio <- make_afun(
161		s_odds_ratio,
162		.formats = c(or_ci = "xx.xx (xx.xx - xx.xx)"),
163		.indent_mods = c(or_ci = 1L)
164		)
165
166		#' @describeIn odds_ratio Layout-creating function which can take statistics function arguments
167		#' and additional format arguments. This function is a wrapper for [rtables::analyze()].
168		#'
169		#' @param ... arguments passed to `s_odds_ratio()`.
170		#'
171		#' @return
172		#' * `estimate_odds_ratio()` returns a layout object suitable for passing to further layouting functions,
173		#' or to [rtables::build_table()]. Adding this function to an `rtable` layout will add formatted rows containing
174		#' the statistics from `s_odds_ratio()` to the table layout.
175		#'
176		#' @examples
177		#' dta <- data.frame(
178		#' rsp = sample(c(TRUE, FALSE), 100, TRUE),
179		#' grp = factor(rep(c("A", "B"), each = 50))
180		#' )
181		#'
182		#' l <- basic_table() %>%
183		#' split_cols_by(var = "grp", ref_group = "B") %>%
184		#' estimate_odds_ratio(vars = "rsp")
185		#'
186		#' build_table(l, df = dta)
187		#'
188		#' @export
189		estimate_odds_ratio <- function(lyt,
190		vars,
191		nested = TRUE,
192		...,
193		show_labels = "hidden",
194		table_names = vars,
195		.stats = "or_ci",
196		.formats = NULL,
197		.labels = NULL,
198		.indent_mods = NULL) {
199	3x	afun <- make_afun(
200	3x	a_odds_ratio,
201	3x	.stats = .stats,
202	3x	.formats = .formats,
203	3x	.labels = .labels,
204	3x	.indent_mods = .indent_mods
205		)
206
207	3x	analyze(
208	3x	lyt,
209	3x	vars,
210	3x	afun = afun,
211	3x	nested = nested,
212	3x	extra_args = list(...),
213	3x	show_labels = show_labels,
214	3x	table_names = table_names
215		)
216		}
217
218		#' Helper Functions for Odds Ratio Estimation
219		#'
220		#' @description `r lifecycle::badge("stable")`
221		#'
222		#' Functions to calculate odds ratios in [estimate_odds_ratio()].
223		#'
224		#' @inheritParams argument_convention
225		#' @param data (`data.frame`)\cr data frame containing at least the variables `rsp` and `grp`, and optionally
226		#' `strata` for [or_clogit()].
227		#'
228		#' @return A named `list` of elements `or_ci` and `n_tot`.
229		#'
230		#' @seealso [odds_ratio]
231		#'
232		#' @name h_odds_ratio
233		NULL
234
235		#' @describeIn h_odds_ratio Estimates the odds ratio based on [stats::glm()]. Note that there must be
236		#' exactly 2 groups in `data` as specified by the `grp` variable.
237		#'
238		#' @examples
239		#' # Data with 2 groups.
240		#' data <- data.frame(
241		#' rsp = as.logical(c(1, 1, 0, 1, 0, 0, 1, 1)),
242		#' grp = letters[c(1, 1, 1, 2, 2, 2, 1, 2)],
243		#' strata = letters[c(1, 2, 1, 2, 2, 2, 1, 2)],
244		#' stringsAsFactors = TRUE
245		#' )
246		#'
247		#' # Odds ratio based on glm.
248		#' or_glm(data, conf_level = 0.95)
249		#'
250		#' @export
251		or_glm <- function(data, conf_level) {
252	55x	checkmate::assert_logical(data$rsp)
253	55x	assert_proportion_value(conf_level)
254	55x	assert_df_with_variables(data, list(rsp = "rsp", grp = "grp"))
255	55x	checkmate::assert_multi_class(data$grp, classes = c("factor", "character"))
256
257	55x	data$grp <- as_factor_keep_attributes(data$grp)
258	55x	assert_df_with_factors(data, list(val = "grp"), min.levels = 2, max.levels = 2)
259	55x	formula <- stats::as.formula("rsp ~ grp")
260	55x	model_fit <- stats::glm(
261	55x	formula = formula, data = data,
262	55x	family = stats::binomial(link = "logit")
263		)
264
265		# Note that here we need to discard the intercept.
266	55x	or <- exp(stats::coef(model_fit)[-1])
267	55x	or_ci <- exp(
268	55x	stats::confint.default(model_fit, level = conf_level)[-1, , drop = FALSE]
269		)
270
271	55x	values <- stats::setNames(c(or, or_ci), c("est", "lcl", "ucl"))
272	55x	n_tot <- stats::setNames(nrow(model_fit$model), "n_tot")
273
274	55x	list(or_ci = values, n_tot = n_tot)
275		}
276
277		#' @describeIn h_odds_ratio estimates the odds ratio based on [survival::clogit()]. This is done for
278		#' the whole data set including all groups, since the results are not the same as when doing
279		#' pairwise comparisons between the groups.
280		#'
281		#' @examples
282		#' # Data with 3 groups.
283		#' data <- data.frame(
284		#' rsp = as.logical(c(1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0)),
285		#' grp = letters[c(1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3)],
286		#' strata = LETTERS[c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2)],
287		#' stringsAsFactors = TRUE
288		#' )
289		#'
290		#' # Odds ratio based on stratified estimation by conditional logistic regression.
291		#' or_clogit(data, conf_level = 0.95)
292		#'
293		#' @export
294		or_clogit <- function(data, conf_level) {
295	16x	checkmate::assert_logical(data$rsp)
296	16x	assert_proportion_value(conf_level)
297	16x	assert_df_with_variables(data, list(rsp = "rsp", grp = "grp", strata = "strata"))
298	16x	checkmate::assert_multi_class(data$grp, classes = c("factor", "character"))
299	16x	checkmate::assert_multi_class(data$strata, classes = c("factor", "character"))
300
301	16x	data$grp <- as_factor_keep_attributes(data$grp)
302	16x	data$strata <- as_factor_keep_attributes(data$strata)
303
304		# Deviation from convention: `survival::strata` must be simply `strata`.
305	16x	formula <- stats::as.formula("rsp ~ grp + strata(strata)")
306	16x	model_fit <- clogit_with_tryCatch(formula = formula, data = data)
307
308		# Create a list with one set of OR estimates and CI per coefficient, i.e.
309		# comparison of one group vs. the reference group.
310	16x	coef_est <- stats::coef(model_fit)
311	16x	ci_est <- stats::confint(model_fit, level = conf_level)
312	16x	or_ci <- list()
313	16x	for (coef_name in names(coef_est)) {
314	18x	grp_name <- gsub("^grp", "", x = coef_name)
315	18x	or_ci[[grp_name]] <- stats::setNames(
316	18x	object = exp(c(coef_est[coef_name], ci_est[coef_name, , drop = TRUE])),
317	18x	nm = c("est", "lcl", "ucl")
318		)
319		}
320	16x	list(or_ci = or_ci, n_tot = c(n_tot = model_fit$n))
321		}

1		#' Incidence Rate
2		#'
3		#' @description `r lifecycle::badge("stable")`
4		#'
5		#' Estimate the event rate adjusted for person-years at risk, otherwise known
6		#' as incidence rate. Primary analysis variable is the person-years at risk.
7		#'
8		#' @inheritParams argument_convention
9		#' @param control (`list`)\cr parameters for estimation details, specified by using
10		#' the helper function [control_incidence_rate()]. Possible parameter options are:
11		#' * `conf_level` (`proportion`)\cr confidence level for the estimated incidence rate.
12		#' * `conf_type` (`string`)\cr `normal` (default), `normal_log`, `exact`, or `byar`
13		#' for confidence interval type.
14		#' * `input_time_unit` (`string`)\cr `day`, `week`, `month`, or `year` (default)
15		#' indicating time unit for data input.
16		#' * `num_pt_year` (`numeric`)\cr time unit for desired output (in person-years).
17		#' @param n_events (`integer`)\cr number of events observed.
18		#'
19		#' @seealso [control_incidence_rate()] and helper functions [h_incidence_rate].
20		#'
21		#' @name incidence_rate
22		NULL
23
24		#' @describeIn incidence_rate Statistics function which estimates the incidence rate and the
25		#' associated confidence interval.
26		#'
27		#' @return
28		#' * `s_incidence_rate()` returns the following statistics:
29		#' - `person_years`: Total person-years at risk.
30		#' - `n_events`: Total number of events observed.
31		#' - `rate`: Estimated incidence rate.
32		#' - `rate_ci`: Confidence interval for the incidence rate.
33		#'
34		#' @examples
35		#' library(dplyr)
36		#'
37		#' df <- data.frame(
38		#' USUBJID = as.character(seq(6)),
39		#' CNSR = c(0, 1, 1, 0, 0, 0),
40		#' AVAL = c(10.1, 20.4, 15.3, 20.8, 18.7, 23.4),
41		#' ARM = factor(c("A", "A", "A", "B", "B", "B"))
42		#' ) %>%
43		#' mutate(is_event = CNSR == 0) %>%
44		#' mutate(n_events = as.integer(is_event))
45		#'
46		#' @keywords internal
47		s_incidence_rate <- function(df,
48		.var,
49		n_events,
50		is_event,
51		control = control_incidence_rate()) {
52	1x	if (!missing(is_event)) {
53	!	warning("argument is_event will be deprecated. Please use n_events.")
54
55	!	if (missing(n_events)) {
56	!	assert_df_with_variables(df, list(tte = .var, is_event = is_event))
57	!	checkmate::assert_string(.var)
58	!	checkmate::assert_logical(df[[is_event]], any.missing = FALSE)
59	!	checkmate::assert_numeric(df[[.var]], any.missing = FALSE)
60	!	n_events <- is_event
61		}
62		} else {
63	1x	assert_df_with_variables(df, list(tte = .var, n_events = n_events))
64	1x	checkmate::assert_string(.var)
65	1x	checkmate::assert_numeric(df[[.var]], any.missing = FALSE)
66	1x	checkmate::assert_integer(df[[n_events]], any.missing = FALSE)
67		}
68
69	1x	input_time_unit <- control$input_time_unit
70	1x	num_pt_year <- control$num_pt_year
71	1x	conf_level <- control$conf_level
72	1x	person_years <- sum(df[[.var]], na.rm = TRUE) * (
73	1x	1 * (input_time_unit == "year") +
74	1x	1 / 12 * (input_time_unit == "month") +
75	1x	1 / 52.14 * (input_time_unit == "week") +
76	1x	1 / 365.24 * (input_time_unit == "day")
77		)
78	1x	n_events <- sum(df[[n_events]], na.rm = TRUE)
79
80	1x	result <- h_incidence_rate(
81	1x	person_years,
82	1x	n_events,
83	1x	control
84		)
85	1x	list(
86	1x	person_years = formatters::with_label(person_years, "Total patient-years at risk"),
87	1x	n_events = formatters::with_label(n_events, "Number of adverse events observed"),
88	1x	rate = formatters::with_label(result$rate, paste("AE rate per", num_pt_year, "patient-years")),
89	1x	rate_ci = formatters::with_label(result$rate_ci, f_conf_level(conf_level))
90		)
91		}
92
93		#' @describeIn incidence_rate Formatted analysis function which is used as `afun`
94		#' in `estimate_incidence_rate()`.
95		#'
96		#' @return
97		#' * `a_incidence_rate()` returns the corresponding list with formatted [rtables::CellValue()].
98		#'
99		#'
100		#' @keywords internal
101		a_incidence_rate <- make_afun(
102		s_incidence_rate,
103		.formats = c(
104		"person_years" = "xx.x",
105		"n_events" = "xx",
106		"rate" = "xx.xx",
107		"rate_ci" = "(xx.xx, xx.xx)"
108		)
109		)
110
111		#' @describeIn incidence_rate Layout-creating function which can take statistics function arguments
112		#' and additional format arguments. This function is a wrapper for [rtables::analyze()].
113		#'
114		#' @return
115		#' * `estimate_incidence_rate()` returns a layout object suitable for passing to further layouting functions,
116		#' or to [rtables::build_table()]. Adding this function to an `rtable` layout will add formatted rows containing
117		#' the statistics from `s_incidence_rate()` to the table layout.
118		#'
119		#' @examples
120		#' basic_table() %>%
121		#' split_cols_by("ARM") %>%
122		#' add_colcounts() %>%
123		#' estimate_incidence_rate(
124		#' vars = "AVAL",
125		#' n_events = "n_events",
126		#' control = control_incidence_rate(
127		#' input_time_unit = "month",
128		#' num_pt_year = 100
129		#' )
130		#' ) %>%
131		#' build_table(df)
132		#'
133		#' @export
134		estimate_incidence_rate <- function(lyt,
135		vars,
136		nested = TRUE,
137		...,
138		show_labels = "hidden",
139		table_names = vars,
140		.stats = NULL,
141		.formats = NULL,
142		.labels = NULL,
143		.indent_mods = NULL) {
144	1x	afun <- make_afun(
145	1x	a_incidence_rate,
146	1x	.stats = .stats,
147	1x	.formats = .formats,
148	1x	.labels = .labels,
149	1x	.indent_mods = .indent_mods
150		)
151
152	1x	analyze(
153	1x	lyt,
154	1x	vars,
155	1x	show_labels = show_labels,
156	1x	table_names = table_names,
157	1x	afun = afun,
158	1x	nested = nested,
159	1x	extra_args = list(...)
160		)
161		}
162
163		#' Helper Functions for Incidence Rate
164		#'
165		#' @description `r lifecycle::badge("stable")`
166		#'
167		#' @param control (`list`)\cr parameters for estimation details, specified by using
168		#' the helper function [control_incidence_rate()]. Possible parameter options are:
169		#' * `conf_level`: (`proportion`)\cr confidence level for the estimated incidence rate.
170		#' * `conf_type`: (`string`)\cr `normal` (default), `normal_log`, `exact`, or `byar`
171		#' for confidence interval type.
172		#' * `input_time_unit`: (`string`)\cr `day`, `week`, `month`, or `year` (default)
173		#' indicating time unit for data input.
174		#' * `num_pt_year`: (`numeric`)\cr time unit for desired output (in person-years).
175		#' @param person_years (`numeric`)\cr total person-years at risk.
176		#' @param alpha (`numeric`)\cr two-sided alpha-level for confidence interval.
177		#' @param n_events (`integer`)\cr number of events observed.
178		#'
179		#' @return Estimated incidence rate `rate` and associated confidence interval `rate_ci`.
180		#'
181		#' @seealso [incidence_rate]
182		#'
183		#' @name h_incidence_rate
184		NULL
185
186		#' @describeIn h_incidence_rate Helper function to estimate the incidence rate and
187		#' associated confidence interval based on the normal approximation for the
188		#' incidence rate. Unit is one person-year.
189		#'
190		#' @examples
191		#' h_incidence_rate_normal(200, 2)
192		#'
193		#' @export
194		h_incidence_rate_normal <- function(person_years,
195		n_events,
196		alpha = 0.05) {
197	1x	checkmate::assert_number(person_years)
198	1x	checkmate::assert_number(n_events)
199	1x	assert_proportion_value(alpha)
200
201	1x	est <- n_events / person_years
202	1x	se <- sqrt(est / person_years)
203	1x	ci <- est + c(-1, 1) * stats::qnorm(1 - alpha / 2) * se
204
205	1x	list(rate = est, rate_ci = ci)
206		}
207
208		#' @describeIn h_incidence_rate Helper function to estimate the incidence rate and
209		#' associated confidence interval based on the normal approximation for the
210		#' logarithm of the incidence rate. Unit is one person-year.
211		#'
212		#' @examples
213		#' h_incidence_rate_normal_log(200, 2)
214		#'
215		#' @export
216		h_incidence_rate_normal_log <- function(person_years,
217		n_events,
218		alpha = 0.05) {
219	5x	checkmate::assert_number(person_years)
220	5x	checkmate::assert_number(n_events)
221	5x	assert_proportion_value(alpha)
222
223	5x	rate_est <- n_events / person_years
224	5x	rate_se <- sqrt(rate_est / person_years)
225	5x	lrate_est <- log(rate_est)
226	5x	lrate_se <- rate_se / rate_est
227	5x	ci <- exp(lrate_est + c(-1, 1) * stats::qnorm(1 - alpha / 2) * lrate_se)
228
229	5x	list(rate = rate_est, rate_ci = ci)
230		}
231
232		#' @describeIn h_incidence_rate Helper function to estimate the incidence rate and
233		#' associated exact confidence interval. Unit is one person-year.
234		#'
235		#' @examples
236		#' h_incidence_rate_exact(200, 2)
237		#'
238		#' @export
239		h_incidence_rate_exact <- function(person_years,
240		n_events,
241		alpha = 0.05) {
242	1x	checkmate::assert_number(person_years)
243	1x	checkmate::assert_number(n_events)
244	1x	assert_proportion_value(alpha)
245
246	1x	est <- n_events / person_years
247	1x	lcl <- stats::qchisq(p = (alpha) / 2, df = 2 * n_events) / (2 * person_years)
248	1x	ucl <- stats::qchisq(p = 1 - (alpha) / 2, df = 2 * n_events + 2) / (2 * person_years)
249
250	1x	list(rate = est, rate_ci = c(lcl, ucl))
251		}
252
253		#' @describeIn h_incidence_rate Helper function to estimate the incidence rate and
254		#' associated `Byar`'s confidence interval. Unit is one person-year.
255		#'
256		#' @examples
257		#' h_incidence_rate_byar(200, 2)
258		#'
259		#' @export
260		h_incidence_rate_byar <- function(person_years,
261		n_events,
262		alpha = 0.05) {
263	1x	checkmate::assert_number(person_years)
264	1x	checkmate::assert_number(n_events)
265	1x	assert_proportion_value(alpha)
266
267	1x	est <- n_events / person_years
268	1x	seg_1 <- n_events + 0.5
269	1x	seg_2 <- 1 - 1 / (9 * (n_events + 0.5))
270	1x	seg_3 <- stats::qnorm(1 - alpha / 2) * sqrt(1 / (n_events + 0.5)) / 3
271	1x	lcl <- seg_1 * ((seg_2 - seg_3)^3) / person_years
272	1x	ucl <- seg_1 * ((seg_2 + seg_3) ^ 3) / person_years # styler: off
273
274	1x	list(rate = est, rate_ci = c(lcl, ucl))
275		}
276
277		#' @describeIn h_incidence_rate Helper function to estimate the incidence rate and
278		#' associated confidence interval.
279		#'
280		#'
281		#' @keywords internal
282		h_incidence_rate <- function(person_years,
283		n_events,
284		control = control_incidence_rate()) {
285	4x	alpha <- 1 - control$conf_level
286	4x	est <- switch(control$conf_type,
287	4x	normal = h_incidence_rate_normal(person_years, n_events, alpha),
288	4x	normal_log = h_incidence_rate_normal_log(person_years, n_events, alpha),
289	4x	exact = h_incidence_rate_exact(person_years, n_events, alpha),
290	4x	byar = h_incidence_rate_byar(person_years, n_events, alpha)
291		)
292
293	4x	num_pt_year <- control$num_pt_year
294	4x	list(
295	4x	rate = est$rate * num_pt_year,
296	4x	rate_ci = est$rate_ci * num_pt_year
297		)
298		}

1		#' Add Titles, Footnotes, Page Number, and a Bounding Box to a Grid Grob
2		#'
3		#' @description `r lifecycle::badge("stable")`
4		#'
5		#' This function is useful to label grid grobs (also `ggplot2`, and `lattice` plots)
6		#' with title, footnote, and page numbers.
7		#'
8		#' @inheritParams grid::grob
9		#' @param grob a grid grob object, optionally `NULL` if only a `grob` with the decoration should be shown.
10		#' @param titles vector of character strings. Vector elements are separated by a newline and strings are wrapped
11		#' according to the page width.
12		#' @param footnotes vector of character string. Same rules as for `titles`.
13		#' @param page string with page numeration, if `NULL` then no page number is displayed.
14		#' @param width_titles unit object
15		#' @param width_footnotes unit object
16		#' @param border boolean, whether a a border should be drawn around the plot or not.
17		#' @param margins unit object of length 4
18		#' @param padding unit object of length 4
19		#' @param outer_margins unit object of length 4
20		#' @param gp_titles a `gpar` object
21		#' @param gp_footnotes a `gpar` object
22		#'
23		#' @return A grid grob (`gTree`).
24		#'
25		#' @details The titles and footnotes will be ragged, i.e. each title will be wrapped individually.
26		#'
27		#' @examples
28		#' library(grid)
29		#'
30		#' titles <- c(
31		#' "Edgar Anderson's Iris Data",
32		#' paste(
33		#' "This famous (Fisher's or Anderson's) iris data set gives the measurements",
34		#' "in centimeters of the variables sepal length and width and petal length",
35		#' "and width, respectively, for 50 flowers from each of 3 species of iris."
36		#' )
37		#' )
38		#'
39		#' footnotes <- c(
40		#' "The species are Iris setosa, versicolor, and virginica.",
41		#' paste(
42		#' "iris is a data frame with 150 cases (rows) and 5 variables (columns) named",
43		#' "Sepal.Length, Sepal.Width, Petal.Length, Petal.Width, and Species."
44		#' )
45		#' )
46		#'
47		#' ## empty plot
48		#' grid.newpage()
49		#'
50		#' grid.draw(
51		#' decorate_grob(
52		#' NULL,
53		#' titles = titles,
54		#' footnotes = footnotes,
55		#' page = "Page 4 of 10"
56		#' )
57		#' )
58		#'
59		#' # grid
60		#' p <- gTree(
61		#' children = gList(
62		#' rectGrob(),
63		#' xaxisGrob(),
64		#' yaxisGrob(),
65		#' textGrob("Sepal.Length", y = unit(-4, "lines")),
66		#' textGrob("Petal.Length", x = unit(-3.5, "lines"), rot = 90),
67		#' pointsGrob(iris$Sepal.Length, iris$Petal.Length, gp = gpar(col = iris$Species), pch = 16)
68		#' ),
69		#' vp = vpStack(plotViewport(), dataViewport(xData = iris$Sepal.Length, yData = iris$Petal.Length))
70		#' )
71		#' grid.newpage()
72		#' grid.draw(p)
73		#'
74		#' grid.newpage()
75		#' grid.draw(
76		#' decorate_grob(
77		#' grob = p,
78		#' titles = titles,
79		#' footnotes = footnotes,
80		#' page = "Page 6 of 129"
81		#' )
82		#' )
83		#'
84		#' ## with ggplot2
85		#' library(ggplot2)
86		#'
87		#' p_gg <- ggplot2::ggplot(iris, aes(Sepal.Length, Sepal.Width, col = Species)) +
88		#' ggplot2::geom_point()
89		#' p_gg
90		#' p <- ggplotGrob(p_gg)
91		#' grid.newpage()
92		#' grid.draw(
93		#' decorate_grob(
94		#' grob = p,
95		#' titles = titles,
96		#' footnotes = footnotes,
97		#' page = "Page 6 of 129"
98		#' )
99		#' )
100		#'
101		#' ## with lattice
102		#' library(lattice)
103		#'
104		#' xyplot(Sepal.Length ~ Petal.Length, data = iris, col = iris$Species)
105		#' p <- grid.grab()
106		#' grid.newpage()
107		#' grid.draw(
108		#' decorate_grob(
109		#' grob = p,
110		#' titles = titles,
111		#' footnotes = footnotes,
112		#' page = "Page 6 of 129"
113		#' )
114		#' )
115		#'
116		#' # with gridExtra - no borders
117		#' library(gridExtra)
118		#' grid.newpage()
119		#' grid.draw(
120		#' decorate_grob(
121		#' tableGrob(
122		#' head(mtcars)
123		#' ),
124		#' titles = "title",
125		#' footnotes = "footnote",
126		#' border = FALSE
127		#' )
128		#' )
129		#'
130		#' @export
131		decorate_grob <- function(grob,
132		titles,
133		footnotes,
134		page = "",
135		width_titles = grid::unit(1, "npc") - grid::stringWidth(page),
136		width_footnotes = grid::unit(1, "npc") - grid::stringWidth(page),
137		border = TRUE,
138		margins = grid::unit(c(1, 0, 1, 0), "lines"),
139		padding = grid::unit(rep(1, 4), "lines"),
140		outer_margins = grid::unit(c(2, 1.5, 3, 1.5), "cm"),
141		gp_titles = grid::gpar(),
142		gp_footnotes = grid::gpar(fontsize = 8),
143		name = NULL,
144		gp = grid::gpar(),
145		vp = NULL) {
146	8x	st_titles <- split_text_grob(
147	8x	titles,
148	8x	x = 0, y = 1,
149	8x	just = c("left", "top"),
150	8x	width = width_titles,
151	8x	vp = grid::viewport(layout.pos.row = 1, layout.pos.col = 1),
152	8x	gp = gp_titles
153		)
154
155	8x	st_footnotes <- split_text_grob(
156	8x	footnotes,
157	8x	x = 0, y = 1,
158	8x	just = c("left", "top"),
159	8x	width = width_footnotes,
160	8x	vp = grid::viewport(layout.pos.row = 3, layout.pos.col = 1),
161	8x	gp = gp_footnotes
162		)
163
164	8x	grid::gTree(
165	8x	grob = grob,
166	8x	titles = titles,
167	8x	footnotes = footnotes,
168	8x	page = page,
169	8x	width_titles = width_titles,
170	8x	width_footnotes = width_footnotes,
171	8x	border = border,
172	8x	margins = margins,
173	8x	padding = padding,
174	8x	outer_margins = outer_margins,
175	8x	gp_titles = gp_titles,
176	8x	gp_footnotes = gp_footnotes,
177	8x	children = grid::gList(
178	8x	grid::gTree(
179	8x	children = grid::gList(
180	8x	st_titles,
181	8x	grid::gTree(
182	8x	children = grid::gList(
183	8x	if (border) grid::rectGrob(),
184	8x	grid::gTree(
185	8x	children = grid::gList(
186	8x	grob
187		),
188	8x	vp = grid::plotViewport(margins = padding)
189		)
190		),
191	8x	vp = grid::vpStack(
192	8x	grid::viewport(layout.pos.row = 2, layout.pos.col = 1),
193	8x	grid::plotViewport(margins = margins)
194		)
195		),
196	8x	st_footnotes,
197	8x	grid::textGrob(
198	8x	page,
199	8x	x = 1, y = 0,
200	8x	just = c("right", "bottom"),
201	8x	vp = grid::viewport(layout.pos.row = 3, layout.pos.col = 1),
202	8x	gp = gp_footnotes
203		)
204		),
205	8x	childrenvp = NULL,
206	8x	name = "titles_grob_footnotes",
207	8x	vp = grid::vpStack(
208	8x	grid::plotViewport(margins = outer_margins),
209	8x	grid::viewport(
210	8x	layout = grid::grid.layout(
211	8x	nrow = 3, ncol = 1,
212	8x	heights = grid::unit.c(
213	8x	grid::grobHeight(st_titles),
214	8x	grid::unit(1, "null"),
215	8x	grid::grobHeight(st_footnotes)
216		)
217		)
218		)
219		)
220		)
221		),
222	8x	name = name,
223	8x	gp = gp,
224	8x	vp = vp,
225	8x	cl = "decoratedGrob"
226		)
227		}
228
229		#' @importFrom grid validDetails
230		#' @noRd
231		validDetails.decoratedGrob <- function(x) {
232	!	checkmate::assert_character(x$titles)
233	!	checkmate::assert_character(x$footnotes)
234
235	!	if (!is.null(x$grob)) {
236	!	checkmate::assert_true(grid::is.grob(x$grob))
237		}
238	!	if (length(x$page) == 1) {
239	!	checkmate::assert_character(x$page)
240		}
241	!	if (!grid::is.unit(x$outer_margins)) {
242	!	checkmate::assert_vector(x$outer_margins, len = 4)
243		}
244	!	if (!grid::is.unit(x$margins)) {
245	!	checkmate::assert_vector(x$margins, len = 4)
246		}
247	!	if (!grid::is.unit(x$padding)) {
248	!	checkmate::assert_vector(x$padding, len = 4)
249		}
250
251	!	x
252		}
253
254		#' @importFrom grid widthDetails
255		#' @noRd
256		widthDetails.decoratedGrob <- function(x) {
257	!	grid::unit(1, "null")
258		}
259
260		#' @importFrom grid heightDetails
261		#' @noRd
262		heightDetails.decoratedGrob <- function(x) {
263	!	grid::unit(1, "null")
264		}
265
266		# Adapted from Paul Murell R Graphics 2nd Edition
267		# https://www.stat.auckland.ac.nz/~paul/RG2e/interactgrid-splittext.R
268		split_string <- function(text, width) {
269	17x	strings <- strsplit(text, " ")
270	17x	out_string <- NA
271	17x	for (string_i in seq_along(strings)) {
272	17x	newline_str <- strings[[string_i]]
273	6x	if (length(newline_str) == 0) newline_str <- ""
274	17x	if (is.na(out_string[string_i])) {
275	17x	out_string[string_i] <- newline_str[[1]][[1]]
276	17x	linewidth <- grid::stringWidth(out_string[string_i])
277		}
278	17x	gapwidth <- grid::stringWidth(" ")
279	17x	availwidth <- as.numeric(width)
280	17x	if (length(newline_str) > 1) {
281	5x	for (i in seq(2, length(newline_str))) {
282	27x	width_i <- grid::stringWidth(newline_str[i])
283	27x	if (grid::convertWidth(linewidth + gapwidth + width_i, grid::unitType(width), valueOnly = TRUE) < availwidth) {
284	25x	sep <- " "
285	25x	linewidth <- linewidth + gapwidth + width_i
286		} else {
287	2x	sep <- "\n"
288	2x	linewidth <- width_i
289		}
290	27x	out_string[string_i] <- paste(out_string[string_i], newline_str[i], sep = sep)
291		}
292		}
293		}
294	17x	paste(out_string, collapse = "\n")
295		}
296
297		#' Split Text According To Available Text Width
298		#'
299		#' Dynamically wrap text.
300		#'
301		#' @inheritParams grid::grid.text
302		#' @param text character string
303		#' @param width a unit object specifying max width of text
304		#'
305		#' @return A text grob.
306		#'
307		#' @details This code is taken from `R Graphics by Paul Murell, 2nd edition`
308		#'
309		#' @keywords internal
310		split_text_grob <- function(text,
311		x = grid::unit(0.5, "npc"),
312		y = grid::unit(0.5, "npc"),
313		width = grid::unit(1, "npc"),
314		just = "centre",
315		hjust = NULL,
316		vjust = NULL,
317		default.units = "npc", # nolint
318		name = NULL,
319		gp = grid::gpar(),
320		vp = NULL) {
321	16x	if (!grid::is.unit(x)) x <- grid::unit(x, default.units)
322	16x	if (!grid::is.unit(y)) y <- grid::unit(y, default.units)
323	!	if (!grid::is.unit(width)) width <- grid::unit(width, default.units)
324	!	if (grid::unitType(x) %in% c("sum", "min", "max")) x <- grid::convertUnit(x, default.units)
325	!	if (grid::unitType(y) %in% c("sum", "min", "max")) y <- grid::convertUnit(y, default.units)
326	16x	if (grid::unitType(width) %in% c("sum", "min", "max")) width <- grid::convertUnit(width, default.units)
327
328		## if it is a fixed unit then we do not need to recalculate when viewport resized
329	16x	if (!inherits(width, "unit.arithmetic") &&
330	16x	!is.null(attr(width, "unit")) &&
331	16x	attr(width, "unit") %in% c("cm", "inches", "mm", "points", "picas", "bigpts", "dida", "cicero", "scaledpts")) {
332	!	attr(text, "fixed_text") <- paste(vapply(text, split_string, character(1), width = width), collapse = "\n")
333		}
334
335	16x	grid::grid.text(
336	16x	label = split_string(text, width),
337	16x	x = x, y = y,
338	16x	just = just,
339	16x	hjust = hjust,
340	16x	vjust = vjust,
341	16x	rot = 0,
342	16x	check.overlap = FALSE,
343	16x	name = name,
344	16x	gp = gp,
345	16x	vp = vp,
346	16x	draw = FALSE
347		)
348		}
349
350		#' @importFrom grid validDetails
351		#' @noRd
352		validDetails.dynamicSplitText <- function(x) {
353	!	checkmate::assert_character(x$text)
354	!	checkmate::assert_true(grid::is.unit(x$width))
355	!	checkmate::assert_vector(x$width, len = 1)
356	!	x
357		}
358
359		#' @importFrom grid heightDetails
360		#' @noRd
361		heightDetails.dynamicSplitText <- function(x) {
362	!	txt <- if (!is.null(attr(x$text, "fixed_text"))) {
363	!	attr(x$text, "fixed_text")
364		} else {
365	!	paste(vapply(x$text, split_string, character(1), width = x$width), collapse = "\n")
366		}
367	!	grid::stringHeight(txt)
368		}
369
370		#' @importFrom grid widthDetails
371		#' @noRd
372		widthDetails.dynamicSplitText <- function(x) {
373	!	x$width
374		}
375
376		#' @importFrom grid drawDetails
377		#' @noRd
378		drawDetails.dynamicSplitText <- function(x, recording) {
379	!	txt <- if (!is.null(attr(x$text, "fixed_text"))) {
380	!	attr(x$text, "fixed_text")
381		} else {
382	!	paste(vapply(x$text, split_string, character(1), width = x$width), collapse = "\n")
383		}
384
385	!	x$width <- NULL
386	!	x$label <- txt
387	!	x$text <- NULL
388	!	class(x) <- c("text", class(x)[-1])
389
390	!	grid::grid.draw(x)
391		}
392
393		#' Update Page Number
394		#'
395		#' Automatically updates page number.
396		#'
397		#' @param npages number of pages in total
398		#' @param ... passed on to [decorate_grob()]
399		#'
400		#' @return Closure that increments the page number.
401		#'
402		#' @keywords internal
403		decorate_grob_factory <- function(npages, ...) {
404	2x	current_page <- 0
405	2x	function(grob) {
406	7x	current_page <<- current_page + 1
407	7x	if (current_page > npages) {
408	1x	stop(paste("current page is", current_page, "but max.", npages, "specified."))
409		}
410	6x	decorate_grob(grob = grob, page = paste("Page", current_page, "of", npages), ...)
411		}
412		}
413
414		#' Decorate Set of `grobs` and Add Page Numbering
415		#'
416		#' @description `r lifecycle::badge("stable")`
417		#'
418		#' Note that this uses the [decorate_grob_factory()] function.
419		#'
420		#' @param grobs a list of grid grobs
421		#' @param ... arguments passed on to [decorate_grob()].
422		#'
423		#' @return A decorated grob.
424		#'
425		#' @examples
426		#' library(ggplot2)
427		#' library(grid)
428		#' g <- with(data = iris, {
429		#' list(
430		#' ggplot2::ggplotGrob(
431		#' ggplot2::ggplot(mapping = aes(Sepal.Length, Sepal.Width, col = Species)) +
432		#' ggplot2::geom_point()
433		#' ),
434		#' ggplot2::ggplotGrob(
435		#' ggplot2::ggplot(mapping = aes(Sepal.Length, Petal.Length, col = Species)) +
436		#' ggplot2::geom_point()
437		#' ),
438		#' ggplot2::ggplotGrob(
439		#' ggplot2::ggplot(mapping = aes(Sepal.Length, Petal.Width, col = Species)) +
440		#' ggplot2::geom_point()
441		#' ),
442		#' ggplot2::ggplotGrob(
443		#' ggplot2::ggplot(mapping = aes(Sepal.Width, Petal.Length, col = Species)) +
444		#' ggplot2::geom_point()
445		#' ),
446		#' ggplot2::ggplotGrob(
447		#' ggplot2::ggplot(mapping = aes(Sepal.Width, Petal.Width, col = Species)) +
448		#' ggplot2::geom_point()
449		#' ),
450		#' ggplot2::ggplotGrob(
451		#' ggplot2::ggplot(mapping = aes(Petal.Length, Petal.Width, col = Species)) +
452		#' ggplot2::geom_point()
453		#' )
454		#' )
455		#' })
456		#' lg <- decorate_grob_set(grobs = g, titles = "Hello\nOne\nTwo\nThree", footnotes = "")
457		#'
458		#' draw_grob(lg[[1]])
459		#' draw_grob(lg[[2]])
460		#' draw_grob(lg[[6]])
461		#'
462		#' @export
463		decorate_grob_set <- function(grobs, ...) {
464	1x	n <- length(grobs)
465	1x	lgf <- decorate_grob_factory(npages = n, ...)
466	1x	lapply(grobs, lgf)
467		}

1		#' Confidence Intervals for a Difference of Binomials
2		#'
3		#' @description `r lifecycle::badge("experimental")`
4		#'
5		#' Several confidence intervals for the difference between proportions.
6		#'
7		#' @name desctools_binom
8		NULL
9
10		#' Recycle List of Parameters
11		#'
12		#' This function recycles all supplied elements to the maximal dimension.
13		#'
14		#' @param ... (`any`)\cr Elements to recycle.
15		#'
16		#' @return A `list`.
17		#'
18		#' @keywords internal
19		#' @noRd
20		h_recycle <- function(...) {
21	60x	lst <- list(...)
22	60x	maxdim <- max(lengths(lst))
23	60x	res <- lapply(lst, rep, length.out = maxdim)
24	60x	attr(res, "maxdim") <- maxdim
25	60x	return(res)
26		}
27
28		#' @describeIn desctools_binom Several confidence intervals for the difference between proportions.
29		#'
30		#' @return A `matrix` of 3 values:
31		#' * `est`: estimate of proportion difference.
32		#' * `lwr.ci`: estimate of lower end of the confidence interval.
33		#' * `upr.ci`: estimate of upper end of the confidence interval.
34		#'
35		#' @keywords internal
36		desctools_binom <- function(x1, n1, x2, n2, conf.level = 0.95, sides = c( # nolint
37		"two.sided",
38		"left", "right"
39		), method = c(
40		"ac", "wald", "waldcc", "score",
41		"scorecc", "mn", "mee", "blj", "ha", "hal", "jp"
42		)) {
43	18x	if (missing(sides)) {
44	18x	sides <- match.arg(sides)
45		}
46	18x	if (missing(method)) {
47	1x	method <- match.arg(method)
48		}
49	18x	iBinomDiffCI <- function(x1, n1, x2, n2, conf.level, sides, # nolint
50	18x	method) {
51	18x	if (sides != "two.sided") {
52	!	conf.level <- 1 - 2 * (1 - conf.level) # nolint
53		}
54	18x	alpha <- 1 - conf.level
55	18x	kappa <- stats::qnorm(1 - alpha / 2)
56	18x	p1_hat <- x1 / n1
57	18x	p2_hat <- x2 / n2
58	18x	est <- p1_hat - p2_hat
59	18x	switch(method,
60	18x	wald = {
61	2x	vd <- p1_hat * (1 - p1_hat) / n1 + p2_hat * (1 - p2_hat) / n2
62	2x	term2 <- kappa * sqrt(vd)
63	2x	ci_lwr <- max(-1, est - term2)
64	2x	ci_upr <- min(1, est + term2)
65		},
66	18x	waldcc = {
67	2x	vd <- p1_hat * (1 - p1_hat) / n1 + p2_hat * (1 - p2_hat) / n2
68	2x	term2 <- kappa * sqrt(vd)
69	2x	term2 <- term2 + 0.5 * (1 / n1 + 1 / n2)
70	2x	ci_lwr <- max(-1, est - term2)
71	2x	ci_upr <- min(1, est + term2)
72		},
73	18x	ac = {
74	2x	n1 <- n1 + 2
75	2x	n2 <- n2 + 2
76	2x	x1 <- x1 + 1
77	2x	x2 <- x2 + 1
78	2x	p1_hat <- x1 / n1
79	2x	p2_hat <- x2 / n2
80	2x	est1 <- p1_hat - p2_hat
81	2x	vd <- p1_hat * (1 - p1_hat) / n1 + p2_hat * (1 - p2_hat) / n2
82	2x	term2 <- kappa * sqrt(vd)
83	2x	ci_lwr <- max(-1, est1 - term2)
84	2x	ci_upr <- min(1, est1 + term2)
85		},
86	18x	exact = {
87	!	ci_lwr <- NA
88	!	ci_upr <- NA
89		},
90	18x	score = {
91	2x	w1 <- desctools_binomci(
92	2x	x = x1, n = n1, conf.level = conf.level,
93	2x	method = "wilson"
94		)
95	2x	w2 <- desctools_binomci(
96	2x	x = x2, n = n2, conf.level = conf.level,
97	2x	method = "wilson"
98		)
99	2x	l1 <- w1[2]
100	2x	u1 <- w1[3]
101	2x	l2 <- w2[2]
102	2x	u2 <- w2[3]
103	2x	ci_lwr <- est - kappa * sqrt(l1 * (1 - l1) / n1 +
104	2x	u2 * (1 - u2) / n2)
105	2x	ci_upr <- est + kappa * sqrt(u1 * (1 - u1) / n1 +
106	2x	l2 * (1 - l2) / n2)
107		},
108	18x	scorecc = {
109	1x	w1 <- desctools_binomci(
110	1x	x = x1, n = n1, conf.level = conf.level,
111	1x	method = "wilsoncc"
112		)
113	1x	w2 <- desctools_binomci(
114	1x	x = x2, n = n2, conf.level = conf.level,
115	1x	method = "wilsoncc"
116		)
117	1x	l1 <- w1[2]
118	1x	u1 <- w1[3]
119	1x	l2 <- w2[2]
120	1x	u2 <- w2[3]
121	1x	ci_lwr <- max(-1, est - sqrt((p1_hat - l1)^2 +
122	1x	(u2 - p2_hat)^2))
123	1x	ci_upr <- min(1, est + sqrt((u1 - p1_hat)^2 + (p2_hat -
124	1x	l2)^2))
125		},
126	18x	mee = {
127	1x	.score <- function(p1, n1, p2, n2, dif) {
128	!	if (dif > 1) dif <- 1
129	!	if (dif < -1) dif <- -1
130	24x	diff <- p1 - p2 - dif
131	24x	if (abs(diff) == 0) {
132	!	res <- 0
133		} else {
134	24x	t <- n2 / n1
135	24x	a <- 1 + t
136	24x	b <- -(1 + t + p1 + t * p2 + dif * (t + 2))
137	24x	c <- dif * dif + dif * (2 * p1 + t + 1) + p1 +
138	24x	t * p2
139	24x	d <- -p1 * dif * (1 + dif)
140	24x	v <- (b / a / 3)^3 - b * c / (6 * a * a) + d / a / 2
141	24x	if (abs(v) < .Machine$double.eps) v <- 0
142	24x	s <- sqrt((b / a / 3)^2 - c / a / 3)
143	24x	u <- ifelse(v > 0, 1, -1) * s
144	24x	w <- (3.141592654 + acos(v / u^3)) / 3
145	24x	p1d <- 2 * u * cos(w) - b / a / 3
146	24x	p2d <- p1d - dif
147	24x	n <- n1 + n2
148	24x	res <- (p1d * (1 - p1d) / n1 + p2d * (1 - p2d) / n2)
149		}
150	24x	return(sqrt(res))
151		}
152	1x	pval <- function(delta) {
153	24x	z <- (est - delta) / .score(
154	24x	p1_hat, n1, p2_hat,
155	24x	n2, delta
156		)
157	24x	2 * min(stats::pnorm(z), 1 - stats::pnorm(z))
158		}
159	1x	ci_lwr <- max(-1, stats::uniroot(function(delta) {
160	12x	pval(delta) -
161	12x	alpha
162	1x	}, interval = c(-1 + 1e-06, est - 1e-06))$root)
163	1x	ci_upr <- min(1, stats::uniroot(function(delta) {
164	12x	pval(delta) -
165	12x	alpha
166	1x	}, interval = c(est + 1e-06, 1 - 1e-06))$root)
167		},
168	18x	blj = {
169	1x	p1_dash <- (x1 + 0.5) / (n1 + 1)
170	1x	p2_dash <- (x2 + 0.5) / (n2 + 1)
171	1x	vd <- p1_dash * (1 - p1_dash) / n1 + p2_dash * (1 -
172	1x	p2_dash) / n2
173	1x	term2 <- kappa * sqrt(vd)
174	1x	est_dash <- p1_dash - p2_dash
175	1x	ci_lwr <- max(-1, est_dash - term2)
176	1x	ci_upr <- min(1, est_dash + term2)
177		},
178	18x	ha = {
179	4x	term2 <- 1 / (2 * min(n1, n2)) + kappa * sqrt(p1_hat *
180	4x	(1 - p1_hat) / (n1 - 1) + p2_hat * (1 - p2_hat) / (n2 -
181	4x	1))
182	4x	ci_lwr <- max(-1, est - term2)
183	4x	ci_upr <- min(1, est + term2)
184		},
185	18x	mn = {
186	1x	.conf <- function(x1, n1, x2, n2, z, lower = FALSE) {
187	2x	p1 <- x1 / n1
188	2x	p2 <- x2 / n2
189	2x	p_hat <- p1 - p2
190	2x	dp <- 1 + ifelse(lower, 1, -1) * p_hat
191	2x	i <- 1
192	2x	while (i <= 50) {
193	46x	dp <- 0.5 * dp
194	46x	y <- p_hat + ifelse(lower, -1, 1) * dp
195	46x	score <- .score(p1, n1, p2, n2, y)
196	46x	if (score < z) {
197	20x	p_hat <- y
198		}
199	46x	if ((dp < 1e-07) \|\| (abs(z - score) < 1e-06)) {
200	2x	(break)()
201		} else {
202	44x	i <- i +
203	44x	1
204		}
205		}
206	2x	return(y)
207		}
208	1x	.score <- function(p1, n1, p2, n2, dif) {
209	46x	diff <- p1 - p2 - dif
210	46x	if (abs(diff) == 0) {
211	!	res <- 0
212		} else {
213	46x	t <- n2 / n1
214	46x	a <- 1 + t
215	46x	b <- -(1 + t + p1 + t * p2 + dif * (t + 2))
216	46x	c <- dif * dif + dif * (2 * p1 + t + 1) + p1 +
217	46x	t * p2
218	46x	d <- -p1 * dif * (1 + dif)
219	46x	v <- (b / a / 3)^3 - b * c / (6 * a * a) + d / a / 2
220	46x	s <- sqrt((b / a / 3)^2 - c / a / 3)
221	46x	u <- ifelse(v > 0, 1, -1) * s
222	46x	w <- (3.141592654 + acos(v / u^3)) / 3
223	46x	p1d <- 2 * u * cos(w) - b / a / 3
224	46x	p2d <- p1d - dif
225	46x	n <- n1 + n2
226	46x	var <- (p1d * (1 - p1d) / n1 + p2d * (1 - p2d) / n2) *
227	46x	n / (n - 1)
228	46x	res <- diff^2 / var
229		}
230	46x	return(res)
231		}
232	1x	z <- stats::qchisq(conf.level, 1)
233	1x	ci_lwr <- max(-1, .conf(x1, n1, x2, n2, z, TRUE))
234	1x	ci_upr <- min(1, .conf(x1, n1, x2, n2, z, FALSE))
235		},
236	18x	beal = {
237	!	a <- p1_hat + p2_hat
238	!	b <- p1_hat - p2_hat
239	!	u <- ((1 / n1) + (1 / n2)) / 4
240	!	v <- ((1 / n1) - (1 / n2)) / 4
241	!	V <- u * ((2 - a) * a - b^2) + 2 * v * (1 - a) * b # nolint
242	!	z <- stats::qchisq(p = 1 - alpha / 2, df = 1)
243	!	A <- sqrt(z * (V + z * u^2 * (2 - a) * a + z * v^2 * (1 - a)^2)) # nolint
244	!	B <- (b + z * v * (1 - a)) / (1 + z * u) # nolint
245	!	ci_lwr <- max(-1, B - A / (1 + z * u))
246	!	ci_upr <- min(1, B + A / (1 + z * u))
247		},
248	18x	hal = {
249	1x	psi <- (p1_hat + p2_hat) / 2
250	1x	u <- (1 / n1 + 1 / n2) / 4
251	1x	v <- (1 / n1 - 1 / n2) / 4
252	1x	z <- kappa
253	1x	theta <- ((p1_hat - p2_hat) + z^2 * v * (1 - 2 *
254	1x	psi)) / (1 + z^2 * u)
255	1x	w <- z / (1 + z^2 * u) * sqrt(u * (4 * psi * (1 - psi) -
256	1x	(p1_hat - p2_hat)^2) + 2 * v * (1 - 2 * psi) *
257	1x	(p1_hat - p2_hat) + 4 * z^2 * u^2 * (1 - psi) *
258	1x	psi + z^2 * v^2 * (1 - 2 * psi)^2)
259	1x	c(theta + w, theta - w)
260	1x	ci_lwr <- max(-1, theta - w)
261	1x	ci_upr <- min(1, theta + w)
262		},
263	18x	jp = {
264	1x	psi <- 0.5 * ((x1 + 0.5) / (n1 + 1) + (x2 + 0.5) / (n2 +
265	1x	1))
266	1x	u <- (1 / n1 + 1 / n2) / 4
267	1x	v <- (1 / n1 - 1 / n2) / 4
268	1x	z <- kappa
269	1x	theta <- ((p1_hat - p2_hat) + z^2 * v * (1 - 2 *
270	1x	psi)) / (1 + z^2 * u)
271	1x	w <- z / (1 + z^2 * u) * sqrt(u * (4 * psi * (1 - psi) -
272	1x	(p1_hat - p2_hat)^2) + 2 * v * (1 - 2 * psi) *
273	1x	(p1_hat - p2_hat) + 4 * z^2 * u^2 * (1 - psi) *
274	1x	psi + z^2 * v^2 * (1 - 2 * psi)^2)
275	1x	c(theta + w, theta - w)
276	1x	ci_lwr <- max(-1, theta - w)
277	1x	ci_upr <- min(1, theta + w)
278		},
279		)
280	18x	ci <- c(
281	18x	est = est, lwr.ci = min(ci_lwr, ci_upr),
282	18x	upr.ci = max(ci_lwr, ci_upr)
283		)
284	18x	if (sides == "left") {
285	!	ci[3] <- 1
286	18x	} else if (sides == "right") {
287	!	ci[2] <- -1
288		}
289	18x	return(ci)
290		}
291	18x	method <- match.arg(arg = method, several.ok = TRUE)
292	18x	sides <- match.arg(arg = sides, several.ok = TRUE)
293	18x	lst <- h_recycle(
294	18x	x1 = x1, n1 = n1, x2 = x2, n2 = n2, conf.level = conf.level,
295	18x	sides = sides, method = method
296		)
297	18x	res <- t(sapply(1:attr(lst, "maxdim"), function(i) {
298	18x	iBinomDiffCI(
299	18x	x1 = lst$x1[i],
300	18x	n1 = lst$n1[i], x2 = lst$x2[i], n2 = lst$n2[i], conf.level = lst$conf.level[i],
301	18x	sides = lst$sides[i], method = lst$method[i]
302		)
303		}))
304	18x	lgn <- h_recycle(x1 = if (is.null(names(x1))) {
305	18x	paste("x1", seq_along(x1), sep = ".")
306		} else {
307	!	names(x1)
308	18x	}, n1 = if (is.null(names(n1))) {
309	18x	paste("n1", seq_along(n1), sep = ".")
310		} else {
311	!	names(n1)
312	18x	}, x2 = if (is.null(names(x2))) {
313	18x	paste("x2", seq_along(x2), sep = ".")
314		} else {
315	!	names(x2)
316	18x	}, n2 = if (is.null(names(n2))) {
317	18x	paste("n2", seq_along(n2), sep = ".")
318		} else {
319	!	names(n2)
320	18x	}, conf.level = conf.level, sides = sides, method = method)
321	18x	xn <- apply(as.data.frame(lgn[sapply(lgn, function(x) {
322	126x	length(unique(x)) !=
323	126x	1
324	18x	})]), 1, paste, collapse = ":")
325	18x	rownames(res) <- xn
326	18x	return(res)
327		}
328
329		#' @describeIn desctools_binom Compute confidence intervals for binomial proportions.
330		#'
331		#' @param x (`count`)\cr number of successes
332		#' @param n (`count`)\cr number of trials
333		#' @param conf.level (`proportion`)\cr confidence level, defaults to 0.95.
334		#' @param sides (`character`)\cr side of the confidence interval to compute. Must be one of `"two-sided"` (default),
335		#' `"left"`, or `"right"`.
336		#' @param method (`character`)\cr method to use. Can be one out of: `"wald"`, `"wilson"`, `"wilsoncc"`,
337		#' `"agresti-coull"`, `"jeffreys"`, `"modified wilson"`, `"modified jeffreys"`, `"clopper-pearson"`, `"arcsine"`,
338		#' `"logit"`, `"witting"`, `"pratt"`, `"midp"`, `"lik"`, and `"blaker"`.
339		#'
340		#' @return A `matrix` with 3 columns containing:
341		#' * `est`: estimate of proportion difference.
342		#' * `lwr.ci`: lower end of the confidence interval.
343		#' * `upr.ci`: upper end of the confidence interval.
344		#'
345		#' @keywords internal
346		desctools_binomci <- function(x,
347		n,
348		conf.level = 0.95, # nolint
349		sides = c("two.sided", "left", "right"),
350		method = c(
351		"wilson", "wald", "waldcc", "agresti-coull",
352		"jeffreys", "modified wilson", "wilsoncc", "modified jeffreys",
353		"clopper-pearson", "arcsine", "logit", "witting", "pratt",
354		"midp", "lik", "blaker"
355		),
356		rand = 123,
357		tol = 1e-05) {
358	24x	if (missing(method)) {
359	1x	method <- "wilson"
360		}
361	24x	if (missing(sides)) {
362	23x	sides <- "two.sided"
363		}
364	24x	iBinomCI <- function(x, n, conf.level = 0.95, sides = c( # nolint
365	24x	"two.sided",
366	24x	"left", "right"
367	24x	), method = c(
368	24x	"wilson", "wilsoncc", "wald",
369	24x	"waldcc", "agresti-coull", "jeffreys", "modified wilson",
370	24x	"modified jeffreys", "clopper-pearson", "arcsine", "logit",
371	24x	"witting", "pratt", "midp", "lik", "blaker"
372	24x	), rand = 123,
373	24x	tol = 1e-05) {
374	24x	if (length(x) != 1) {
375	!	stop("'x' has to be of length 1 (number of successes)")
376		}
377	24x	if (length(n) != 1) {
378	!	stop("'n' has to be of length 1 (number of trials)")
379		}
380	24x	if (length(conf.level) != 1) {
381	!	stop("'conf.level' has to be of length 1 (confidence level)")
382		}
383	24x	if (conf.level < 0.5 \|\| conf.level > 1) {
384	!	stop("'conf.level' has to be in [0.5, 1]")
385		}
386	24x	sides <- match.arg(sides, choices = c(
387	24x	"two.sided", "left",
388	24x	"right"
389	24x	), several.ok = FALSE)
390	24x	if (sides != "two.sided") {
391	1x	conf.level <- 1 - 2 * (1 - conf.level) # nolint
392		}
393	24x	alpha <- 1 - conf.level
394	24x	kappa <- stats::qnorm(1 - alpha / 2)
395	24x	p_hat <- x / n
396	24x	q_hat <- 1 - p_hat
397	24x	est <- p_hat
398	24x	switch(match.arg(arg = method, choices = c(
399	24x	"wilson",
400	24x	"wald", "waldcc", "wilsoncc", "agresti-coull", "jeffreys",
401	24x	"modified wilson", "modified jeffreys", "clopper-pearson",
402	24x	"arcsine", "logit", "witting", "pratt", "midp", "lik",
403	24x	"blaker"
404		)),
405	24x	wald = {
406	1x	term2 <- kappa * sqrt(p_hat * q_hat) / sqrt(n)
407	1x	ci_lwr <- max(0, p_hat - term2)
408	1x	ci_upr <- min(1, p_hat + term2)
409		},
410	24x	waldcc = {
411	1x	term2 <- kappa * sqrt(p_hat * q_hat) / sqrt(n)
412	1x	term2 <- term2 + 1 / (2 * n)
413	1x	ci_lwr <- max(0, p_hat - term2)
414	1x	ci_upr <- min(1, p_hat + term2)
415		},
416	24x	wilson = {
417	6x	term1 <- (x + kappa^2 / 2) / (n + kappa^2)
418	6x	term2 <- kappa * sqrt(n) / (n + kappa^2) * sqrt(p_hat *
419	6x	q_hat + kappa^2 / (4 * n))
420	6x	ci_lwr <- max(0, term1 - term2)
421	6x	ci_upr <- min(1, term1 + term2)
422		},
423	24x	wilsoncc = {
424	3x	lci <- (2 * x + kappa^2 - 1 - kappa * sqrt(kappa^2 -
425	3x	2 - 1 / n + 4 * p_hat * (n * q_hat + 1))) / (2 *
426	3x	(n + kappa^2))
427	3x	uci <- (2 * x + kappa^2 + 1 + kappa * sqrt(kappa^2 +
428	3x	2 - 1 / n + 4 * p_hat * (n * q_hat - 1))) / (2 *
429	3x	(n + kappa^2))
430	3x	ci_lwr <- max(0, ifelse(p_hat == 0, 0, lci))
431	3x	ci_upr <- min(1, ifelse(p_hat == 1, 1, uci))
432		},
433	24x	`agresti-coull` = {
434	1x	x_tilde <- x + kappa^2 / 2
435	1x	n_tilde <- n + kappa^2
436	1x	p_tilde <- x_tilde / n_tilde
437	1x	q_tilde <- 1 - p_tilde
438	1x	est <- p_tilde
439	1x	term2 <- kappa * sqrt(p_tilde * q_tilde) / sqrt(n_tilde)
440	1x	ci_lwr <- max(0, p_tilde - term2)
441	1x	ci_upr <- min(1, p_tilde + term2)
442		},
443	24x	jeffreys = {
444	1x	if (x == 0) {
445	!	ci_lwr <- 0
446		} else {
447	1x	ci_lwr <- stats::qbeta(
448	1x	alpha / 2,
449	1x	x + 0.5, n - x + 0.5
450		)
451		}
452	1x	if (x == n) {
453	!	ci_upr <- 1
454		} else {
455	1x	ci_upr <- stats::qbeta(1 -
456	1x	alpha / 2, x + 0.5, n - x + 0.5)
457		}
458		},
459	24x	`modified wilson` = {
460	1x	term1 <- (x + kappa^2 / 2) / (n + kappa^2)
461	1x	term2 <- kappa * sqrt(n) / (n + kappa^2) * sqrt(p_hat *
462	1x	q_hat + kappa^2 / (4 * n))
463	1x	if ((n <= 50 & x %in% c(1, 2)) \| (n >= 51 & x %in%
464	1x	c(1:3))) {
465	!	ci_lwr <- 0.5 * stats::qchisq(alpha, 2 *
466	!	x) / n
467		} else {
468	1x	ci_lwr <- max(0, term1 - term2)
469		}
470	1x	if ((n <= 50 & x %in% c(n - 1, n - 2)) \| (n >= 51 &
471	1x	x %in% c(n - (1:3)))) {
472	!	ci_upr <- 1 - 0.5 * stats::qchisq(
473	!	alpha,
474	!	2 * (n - x)
475	!	) / n
476		} else {
477	1x	ci_upr <- min(1, term1 +
478	1x	term2)
479		}
480		},
481	24x	`modified jeffreys` = {
482	1x	if (x == n) {
483	!	ci_lwr <- (alpha / 2)^(1 / n)
484		} else {
485	1x	if (x <= 1) {
486	!	ci_lwr <- 0
487		} else {
488	1x	ci_lwr <- stats::qbeta(
489	1x	alpha / 2,
490	1x	x + 0.5, n - x + 0.5
491		)
492		}
493		}
494	1x	if (x == 0) {
495	!	ci_upr <- 1 - (alpha / 2)^(1 / n)
496		} else {
497	1x	if (x >= n - 1) {
498	!	ci_upr <- 1
499		} else {
500	1x	ci_upr <- stats::qbeta(1 -
501	1x	alpha / 2, x + 0.5, n - x + 0.5)
502		}
503		}
504		},
505	24x	`clopper-pearson` = {
506	1x	ci_lwr <- stats::qbeta(alpha / 2, x, n - x + 1)
507	1x	ci_upr <- stats::qbeta(1 - alpha / 2, x + 1, n - x)
508		},
509	24x	arcsine = {
510	1x	p_tilde <- (x + 0.375) / (n + 0.75)
511	1x	est <- p_tilde
512	1x	ci_lwr <- sin(asin(sqrt(p_tilde)) - 0.5 * kappa / sqrt(n))^2
513	1x	ci_upr <- sin(asin(sqrt(p_tilde)) + 0.5 * kappa / sqrt(n))^2
514		},
515	24x	logit = {
516	1x	lambda_hat <- log(x / (n - x))
517	1x	V_hat <- n / (x * (n - x)) # nolint
518	1x	lambda_lower <- lambda_hat - kappa * sqrt(V_hat)
519	1x	lambda_upper <- lambda_hat + kappa * sqrt(V_hat)
520	1x	ci_lwr <- exp(lambda_lower) / (1 + exp(lambda_lower))
521	1x	ci_upr <- exp(lambda_upper) / (1 + exp(lambda_upper))
522		},
523	24x	witting = {
524	1x	set.seed(rand)
525	1x	x_tilde <- x + stats::runif(1, min = 0, max = 1)
526	1x	pbinom_abscont <- function(q, size, prob) {
527	22x	v <- trunc(q)
528	22x	term1 <- stats::pbinom(v - 1, size = size, prob = prob)
529	22x	term2 <- (q - v) * stats::dbinom(v, size = size, prob = prob)
530	22x	return(term1 + term2)
531		}
532	1x	qbinom_abscont <- function(p, size, x) {
533	2x	fun <- function(prob, size, x, p) {
534	22x	pbinom_abscont(x, size, prob) - p
535		}
536	2x	stats::uniroot(fun,
537	2x	interval = c(0, 1), size = size,
538	2x	x = x, p = p
539	2x	)$root
540		}
541	1x	ci_lwr <- qbinom_abscont(1 - alpha, size = n, x = x_tilde)
542	1x	ci_upr <- qbinom_abscont(alpha, size = n, x = x_tilde)
543		},
544	24x	pratt = {
545	1x	if (x == 0) {
546	!	ci_lwr <- 0
547	!	ci_upr <- 1 - alpha^(1 / n)
548	1x	} else if (x == 1) {
549	!	ci_lwr <- 1 - (1 - alpha / 2)^(1 / n)
550	!	ci_upr <- 1 - (alpha / 2)^(1 / n)
551	1x	} else if (x == (n - 1)) {
552	!	ci_lwr <- (alpha / 2)^(1 / n)
553	!	ci_upr <- (1 - alpha / 2)^(1 / n)
554	1x	} else if (x == n) {
555	!	ci_lwr <- alpha^(1 / n)
556	!	ci_upr <- 1
557		} else {
558	1x	z <- stats::qnorm(1 - alpha / 2)
559	1x	A <- ((x + 1) / (n - x))^2 # nolint
560	1x	B <- 81 * (x + 1) * (n - x) - 9 * n - 8 # nolint
561	1x	C <- (0 - 3) * z * sqrt(9 * (x + 1) * (n - x) * (9 * n + 5 - z^2) + n + 1) # nolint
562	1x	D <- 81 * (x + 1)^2 - 9 * (x + 1) * (2 + z^2) + 1 # nolint
563	1x	E <- 1 + A * ((B + C) / D)^3 # nolint
564	1x	ci_upr <- 1 / E
565	1x	A <- (x / (n - x - 1))^2 # nolint
566	1x	B <- 81 * x * (n - x - 1) - 9 * n - 8 # nolint
567	1x	C <- 3 * z * sqrt(9 * x * (n - x - 1) * (9 * n + 5 - z^2) + n + 1) # nolint
568	1x	D <- 81 * x^2 - 9 * x * (2 + z^2) + 1 # nolint
569	1x	E <- 1 + A * ((B + C) / D)^3 # nolint
570	1x	ci_lwr <- 1 / E
571		}
572		},
573	24x	midp = {
574	1x	f_low <- function(pi, x, n) {
575	12x	1 / 2 * stats::dbinom(x, size = n, prob = pi) + stats::pbinom(x,
576	12x	size = n, prob = pi, lower.tail = FALSE
577		) -
578	12x	(1 - conf.level) / 2
579		}
580	1x	f_up <- function(pi, x, n) {
581	12x	1 / 2 * stats::dbinom(x, size = n, prob = pi) + stats::pbinom(x -
582	12x	1, size = n, prob = pi) - (1 - conf.level) / 2
583		}
584	1x	ci_lwr <- 0
585	1x	ci_upr <- 1
586	1x	if (x != 0) {
587	1x	ci_lwr <- stats::uniroot(f_low,
588	1x	interval = c(0, p_hat),
589	1x	x = x, n = n
590	1x	)$root
591		}
592	1x	if (x != n) {
593	1x	ci_upr <- stats::uniroot(f_up, interval = c(
594	1x	p_hat,
595	1x	1
596	1x	), x = x, n = n)$root
597		}
598		},
599	24x	lik = {
600	2x	ci_lwr <- 0
601	2x	ci_upr <- 1
602	2x	z <- stats::qnorm(1 - alpha * 0.5)
603	2x	tol <- .Machine$double.eps^0.5
604	2x	BinDev <- function(y, x, mu, wt, bound = 0, tol = .Machine$double.eps^0.5, # nolint
605		...) {
606	40x	ll_y <- ifelse(y %in% c(0, 1), 0, stats::dbinom(x, wt,
607	40x	y,
608	40x