tern coverage - 95.64%

Files
Source

#' Formatting functions
#'
#' See below for the list of formatting functions created in `tern` to work with `rtables`.
#'
#' Other available formats can be listed via [`formatters::list_valid_format_labels()`]. Additional
#' custom formats can be created via the [`formatters::sprintf_format()`] function.
#'
#' @family formatting functions
#' @name formatting_functions
NULL

#' Format fraction and percentage
#'
#' @description `r lifecycle::badge("stable")`
#'
#' Formats a fraction together with ratio in percent.
#'
#' @param x (named `integer`)\cr vector with elements `num` and `denom`.
#' @param ... not used. Required for `rtables` interface.
#'
#' @return A string in the format `num / denom (ratio %)`. If `num` is 0, the format is `num / denom`.
#'
#' @examples
#' format_fraction(x = c(num = 2L, denom = 3L))
#' format_fraction(x = c(num = 0L, denom = 3L))
#'
#' @family formatting functions
#' @export
format_fraction <- function(x, ...) {
  attr(x, "label") <- NULL

  checkmate::assert_vector(x)
  checkmate::assert_count(x["num"])
  checkmate::assert_count(x["denom"])

  result <- if (x["num"] == 0) {
    paste0(x["num"], "/", x["denom"])
  } else {
    paste0(
      x["num"], "/", x["denom"],
      " (", round(x["num"] / x["denom"] * 100, 1), "%)"
    )
  }

  return(result)
}

#' Format fraction and percentage with fixed single decimal place
#'
#' @description `r lifecycle::badge("stable")`
#'
#' Formats a fraction together with ratio in percent with fixed single decimal place.
#' Includes trailing zero in case of whole number percentages to always keep one decimal place.
#'
#' @inheritParams format_fraction
#'
#' @return A string in the format `num / denom (ratio %)`. If `num` is 0, the format is `num / denom`.
#'
#' @examples
#' format_fraction_fixed_dp(x = c(num = 1L, denom = 2L))
#' format_fraction_fixed_dp(x = c(num = 1L, denom = 4L))
#' format_fraction_fixed_dp(x = c(num = 0L, denom = 3L))
#'
#' @family formatting functions
#' @export
format_fraction_fixed_dp <- function(x, ...) {
  attr(x, "label") <- NULL
  checkmate::assert_vector(x)
  checkmate::assert_count(x["num"])
  checkmate::assert_count(x["denom"])

  result <- if (x["num"] == 0) {
    paste0(x["num"], "/", x["denom"])
  } else {
    paste0(
      x["num"], "/", x["denom"],
      " (", sprintf("%.1f", round(x["num"] / x["denom"] * 100, 1)), "%)"
    )
  }
  return(result)
}

#' Format count and fraction
#'
#' @description `r lifecycle::badge("stable")`
#'
#' Formats a count together with fraction with special consideration when count is `0`.
#'
#' @param x (`numeric(2)`)\cr vector of length 2 with count and fraction, respectively.
#' @param ... not used. Required for `rtables` interface.
#'
#' @return A string in the format `count (fraction %)`. If `count` is 0, the format is `0`.
#'
#' @examples
#' format_count_fraction(x = c(2, 0.6667))
#' format_count_fraction(x = c(0, 0))
#'
#' @family formatting functions
#' @export
format_count_fraction <- function(x, ...) {
  attr(x, "label") <- NULL

  if (any(is.na(x))) {
    return("NA")
  }

  checkmate::assert_vector(x)
  checkmate::assert_integerish(x[1])
  assert_proportion_value(x[2], include_boundaries = TRUE)

  result <- if (x[1] == 0) {
    "0"
  } else {
    paste0(x[1], " (", round(x[2] * 100, 1), "%)")
  }

  return(result)
}

#' Format count and percentage with fixed single decimal place
#'
#' @description `r lifecycle::badge("experimental")`
#'
#' Formats a count together with fraction with special consideration when count is `0`.
#'
#' @inheritParams format_count_fraction
#'
#' @return A string in the format `count (fraction %)`. If `count` is 0, the format is `0`.
#'
#' @examples
#' format_count_fraction_fixed_dp(x = c(2, 0.6667))
#' format_count_fraction_fixed_dp(x = c(2, 0.5))
#' format_count_fraction_fixed_dp(x = c(0, 0))
#'
#' @family formatting functions
#' @export
format_count_fraction_fixed_dp <- function(x, ...) {
  attr(x, "label") <- NULL

  if (any(is.na(x))) {
    return("NA")
  }

  checkmate::assert_vector(x)
  checkmate::assert_integerish(x[1])
  assert_proportion_value(x[2], include_boundaries = TRUE)

  result <- if (x[1] == 0) {
    "0"
  } else if (.is_equal_float(x[2], 1)) {
    sprintf("%d (100%%)", x[1])
  } else {
    sprintf("%d (%.1f%%)", x[1], x[2] * 100)
  }

  return(result)
}

#' Format count and fraction with special case for count < 10
#'
#' @description `r lifecycle::badge("stable")`
#'
#' Formats a count together with fraction with special consideration when count is less than 10.
#'
#' @inheritParams format_count_fraction
#'
#' @return A string in the format `count (fraction %)`. If `count` is less than 10, only `count` is printed.
#'
#' @examples
#' format_count_fraction_lt10(x = c(275, 0.9673))
#' format_count_fraction_lt10(x = c(2, 0.6667))
#' format_count_fraction_lt10(x = c(9, 1))
#'
#' @family formatting functions
#' @export
format_count_fraction_lt10 <- function(x, ...) {
  attr(x, "label") <- NULL

  if (any(is.na(x))) {
    return("NA")
  }

  checkmate::assert_vector(x)
  checkmate::assert_integerish(x[1])
  assert_proportion_value(x[2], include_boundaries = TRUE)

  result <- if (x[1] < 10) {
    paste0(x[1])
  } else {
    paste0(x[1], " (", round(x[2] * 100, 1), "%)")
  }

  return(result)
}

#' Format XX as a formatting function
#'
#' Translate a string where x and dots are interpreted as number place
#' holders, and others as formatting elements.
#'
#' @param str (`string`)\cr template.
#'
#' @return An `rtables` formatting function.
#'
#' @examples
#' test <- list(c(1.658, 0.5761), c(1e1, 785.6))
#'
#' z <- format_xx("xx (xx.x)")
#' sapply(test, z)
#'
#' z <- format_xx("xx.x - xx.x")
#' sapply(test, z)
#'
#' z <- format_xx("xx.x, incl. xx.x% NE")
#' sapply(test, z)
#'
#' @family formatting functions
#' @export
format_xx <- function(str) {
  # Find position in the string.
  positions <- gregexpr(pattern = "x+\\.x+|x+", text = str, perl = TRUE)
  x_positions <- regmatches(x = str, m = positions)[[1]]

  # Roundings depends on the number of x behind [.].
  roundings <- lapply(
    X = x_positions,
    function(x) {
      y <- strsplit(split = "\\.", x = x)[[1]]
      rounding <- function(x) {
        round(x, digits = ifelse(length(y) > 1, nchar(y[2]), 0))
      }
      return(rounding)
    }
  )

  rtable_format <- function(x, output) {
    values <- Map(y = x, fun = roundings, function(y, fun) fun(y))
    regmatches(x = str, m = positions)[[1]] <- values
    return(str)
  }

  return(rtable_format)
}

#' Format numeric values by significant figures
#'
#' Format numeric values to print with a specified number of significant figures.
#'
#' @param sigfig (`integer(1)`)\cr number of significant figures to display.
#' @param format (`string`)\cr the format label (string) to apply when printing the value. Decimal
#'   places in string are ignored in favor of formatting by significant figures. Formats options are:
#'   `"xx"`, `"xx / xx"`, `"(xx, xx)"`, `"xx - xx"`, and `"xx (xx)"`.
#' @param num_fmt (`string`)\cr numeric format modifiers to apply to the value. Defaults to `"fg"` for
#'   standard significant figures formatting - fixed (non-scientific notation) format (`"f"`)
#'   and `sigfig` equal to number of significant figures instead of decimal places (`"g"`). See the
#'   [formatC()] `format` argument for more options.
#'
#' @return An `rtables` formatting function.
#'
#' @examples
#' fmt_3sf <- format_sigfig(3)
#' fmt_3sf(1.658)
#' fmt_3sf(1e1)
#'
#' fmt_5sf <- format_sigfig(5)
#' fmt_5sf(0.57)
#' fmt_5sf(0.000025645)
#'
#' @family formatting functions
#' @export
format_sigfig <- function(sigfig, format = "xx", num_fmt = "fg") {
  checkmate::assert_integerish(sigfig)
  format <- gsub("xx\\.|xx\\.x+", "xx", format)
  checkmate::assert_choice(format, c("xx", "xx / xx", "(xx, xx)", "xx - xx", "xx (xx)"))
  function(x, ...) {
    if (!is.numeric(x)) stop("`format_sigfig` cannot be used for non-numeric values. Please choose another format.")
    num <- formatC(signif(x, digits = sigfig), digits = sigfig, format = num_fmt, flag = "#")
    num <- gsub("\\.$", "", num) # remove trailing "."

    format_value(num, format)
  }
}

#' Format fraction with lower threshold
#'
#' @description `r lifecycle::badge("stable")`
#'
#' Formats a fraction when the second element of the input `x` is the fraction. It applies
#' a lower threshold, below which it is just stated that the fraction is smaller than that.
#'
#' @param threshold (`proportion`)\cr lower threshold.
#'
#' @return An `rtables` formatting function that takes numeric input `x` where the second
#'   element is the fraction that is formatted. If the fraction is above or equal to the threshold,
#'   then it is displayed in percentage. If it is positive but below the threshold, it returns,
#'   e.g. "<1" if the threshold is `0.01`. If it is zero, then just "0" is returned.
#'
#' @examples
#' format_fun <- format_fraction_threshold(0.05)
#' format_fun(x = c(20, 0.1))
#' format_fun(x = c(2, 0.01))
#' format_fun(x = c(0, 0))
#'
#' @family formatting functions
#' @export
format_fraction_threshold <- function(threshold) {
  assert_proportion_value(threshold)
  string_below_threshold <- paste0("<", round(threshold * 100))
  function(x, ...) {
    assert_proportion_value(x[2], include_boundaries = TRUE)
    ifelse(
      x[2] > 0.01,
      round(x[2] * 100),
      ifelse(
        x[2] == 0,
        "0",
        string_below_threshold
      )
    )
  }
}

#' Format extreme values
#'
#' @description `r lifecycle::badge("stable")`
#'
#' `rtables` formatting functions that handle extreme values.
#'
#' @param digits (`integer(1)`)\cr number of decimal places to display.
#'
#' @details For each input, apply a format to the specified number of `digits`. If the value is
#'    below a threshold, it returns "<0.01" e.g. if the number of `digits` is 2. If the value is
#'    above a threshold, it returns ">999.99" e.g. if the number of `digits` is 2.
#'    If it is zero, then returns "0.00".
#'
#' @family formatting functions
#' @name extreme_format
NULL

#' @describeIn extreme_format Internal helper function to calculate the threshold and create formatted strings
#'  used in Formatting Functions. Returns a list with elements `threshold` and `format_string`.
#'
#' @return
#' * `h_get_format_threshold()` returns a `list` of 2 elements: `threshold`, with `low` and `high` thresholds,
#'   and `format_string`, with thresholds formatted as strings.
#'
#' @examples
#' h_get_format_threshold(2L)
#'
#' @export
h_get_format_threshold <- function(digits = 2L) {
  checkmate::assert_integerish(digits)

  low_threshold <- 1 / (10 ^ digits) # styler: off
  high_threshold <- 1000 - (1 / (10 ^ digits)) # styler: off

  string_below_threshold <- paste0("<", low_threshold)
  string_above_threshold <- paste0(">", high_threshold)

  list(
    "threshold" = c(low = low_threshold, high = high_threshold),
    "format_string" = c(low = string_below_threshold, high = string_above_threshold)
  )
}

#' @describeIn extreme_format Internal helper function to apply a threshold format to a value.
#'   Creates a formatted string to be used in Formatting Functions.
#'
#' @param x (`numeric(1)`)\cr value to format.
#'
#' @return
#' * `h_format_threshold()` returns the given value, or if the value is not within the digit threshold the relation
#'   of the given value to the digit threshold, as a formatted string.
#'
#' @examples
#' h_format_threshold(0.001)
#' h_format_threshold(1000)
#'
#' @export
h_format_threshold <- function(x, digits = 2L) {
  if (is.na(x)) {
    return(x)
  }

  checkmate::assert_numeric(x, lower = 0)

  l_fmt <- h_get_format_threshold(digits)

  result <- if (x < l_fmt$threshold["low"] && 0 < x) {
    l_fmt$format_string["low"]
  } else if (x > l_fmt$threshold["high"]) {
    l_fmt$format_string["high"]
  } else {
    sprintf(fmt = paste0("%.", digits, "f"), x)
  }

  unname(result)
}

#' Format a single extreme value
#'
#' @description `r lifecycle::badge("stable")`
#'
#' Create a formatting function for a single extreme value.
#'
#' @inheritParams extreme_format
#'
#' @return An `rtables` formatting function that uses threshold `digits` to return a formatted extreme value.
#'
#' @examples
#' format_fun <- format_extreme_values(2L)
#' format_fun(x = 0.127)
#' format_fun(x = Inf)
#' format_fun(x = 0)
#' format_fun(x = 0.009)
#'
#' @family formatting functions
#' @export
format_extreme_values <- function(digits = 2L) {
  function(x, ...) {
    checkmate::assert_scalar(x, na.ok = TRUE)

    h_format_threshold(x = x, digits = digits)
  }
}

#' Format extreme values part of a confidence interval
#'
#' @description `r lifecycle::badge("stable")`
#'
#' Formatting Function for extreme values part of a confidence interval. Values
#' are formatted as e.g. "(xx.xx, xx.xx)" if the number of `digits` is 2.
#'
#' @inheritParams extreme_format
#'
#' @return An `rtables` formatting function that uses threshold `digits` to return a formatted extreme
#'   values confidence interval.
#'
#' @examples
#' format_fun <- format_extreme_values_ci(2L)
#' format_fun(x = c(0.127, Inf))
#' format_fun(x = c(0, 0.009))
#'
#' @family formatting functions
#' @export
format_extreme_values_ci <- function(digits = 2L) {
  function(x, ...) {
    checkmate::assert_vector(x, len = 2)
    l_result <- h_format_threshold(x = x[1], digits = digits)
    h_result <- h_format_threshold(x = x[2], digits = digits)

    paste0("(", l_result, ", ", h_result, ")")
  }
}

#' Format automatically using data significant digits
#'
#' @description `r lifecycle::badge("stable")`
#'
#' Formatting function for the majority of default methods used in [analyze_vars()].
#' For non-derived values, the significant digits of data is used (e.g. range), while derived
#' values have one more digits (measure of location and dispersion like mean, standard deviation).
#' This function can be called internally with "auto" like, for example,
#' `.formats = c("mean" = "auto")`. See details to see how this works with the inner function.
#'
#' @param dt_var (`numeric`)\cr variable data the statistics were calculated from. Used only to
#'   find significant digits. In [analyze_vars] this comes from `.df_row` (see
#'   [rtables::additional_fun_params]), and it is the row data after the above row splits. No
#'   column split is considered.
#' @param x_stat (`string`)\cr string indicating the current statistical method used.
#'
#' @return A string that `rtables` prints in a table cell.
#'
#' @details
#' The internal function is needed to work with `rtables` default structure for
#' format functions, i.e. `function(x, ...)`, where is x are results from statistical evaluation.
#' It can be more than one element (e.g. for `.stats = "mean_sd"`).
#'
#' @examples
#' x_todo <- c(0.001, 0.2, 0.0011000, 3, 4)
#' res <- c(mean(x_todo[1:3]), sd(x_todo[1:3]))
#'
#' # x is the result coming into the formatting function -> res!!
#' format_auto(dt_var = x_todo, x_stat = "mean_sd")(x = res)
#' format_auto(x_todo, "range")(x = range(x_todo))
#' no_sc_x <- c(0.0000001, 1)
#' format_auto(no_sc_x, "range")(x = no_sc_x)
#'
#' @family formatting functions
#' @export
format_auto <- function(dt_var, x_stat) {
  function(x = "", ...) {
    checkmate::assert_numeric(x, min.len = 1)
    checkmate::assert_numeric(dt_var, min.len = 1)
    # Defaults - they may be a param in the future
    der_stats <- c(
      "mean", "sd", "se", "median", "geom_mean", "quantiles", "iqr",
      "mean_sd", "mean_se", "mean_se", "mean_ci", "mean_sei", "mean_sdi",
      "median_ci"
    )
    nonder_stats <- c("n", "range", "min", "max")

    # Safenet for miss-modifications
    stopifnot(length(intersect(der_stats, nonder_stats)) == 0) # nolint
    checkmate::assert_choice(x_stat, c(der_stats, nonder_stats))

    # Finds the max number of digits in data
    detect_dig <- vapply(dt_var, count_decimalplaces, FUN.VALUE = numeric(1)) %>%
      max()

    if (x_stat %in% der_stats) {
      detect_dig <- detect_dig + 1
    }

    # Render input
    str_vals <- formatC(x, digits = detect_dig, format = "f")
    def_fmt <- get_formats_from_stats(x_stat)[[x_stat]]
    str_fmt <- str_extract(def_fmt, invert = FALSE)[[1]]
    if (length(str_fmt) != length(str_vals)) {
      stop(
        "Number of inserted values as result (", length(str_vals),
        ") is not the same as there should be in the default tern formats for ",
        x_stat, " (-> ", def_fmt, " needs ", length(str_fmt), " values). ",
        "See tern_default_formats to check all of them."
      )
    }

    # Squashing them together
    inv_str_fmt <- str_extract(def_fmt, invert = TRUE)[[1]]
    stopifnot(length(inv_str_fmt) == length(str_vals) + 1) # nolint

    out <- vector("character", length = length(inv_str_fmt) + length(str_vals))
    is_even <- seq_along(out) %% 2 == 0
    out[is_even] <- str_vals
    out[!is_even] <- inv_str_fmt

    return(paste0(out, collapse = ""))
  }
}

# Utility function that could be useful in general
str_extract <- function(string, pattern = "xx|xx\\.|xx\\.x+", invert = FALSE) {
  regmatches(string, gregexpr(pattern, string), invert = invert)
}

# Helper function
count_decimalplaces <- function(dec) {
  if (is.na(dec)) {
    return(0)
  } else if (abs(dec - round(dec)) > .Machine$double.eps^0.5) { # For precision
    nchar(strsplit(format(dec, scientific = FALSE, trim = FALSE), ".", fixed = TRUE)[[1]][[2]])
  } else {
    return(0)
  }
}

#' Apply automatic formatting
#'
#' Checks if any of the listed formats in `.formats` are `"auto"`, and replaces `"auto"` with
#' the correct implementation of `format_auto` for the given statistics, data, and variable.
#'
#' @inheritParams argument_convention
#' @param x_stats (named `list`)\cr a named list of statistics where each element corresponds
#'   to an element in `.formats`, with matching names.
#'
#' @keywords internal
apply_auto_formatting <- function(.formats, x_stats, .df_row, .var) {
  is_auto_fmt <- vapply(.formats, function(ii) is.character(ii) && ii == "auto", logical(1))
  if (any(is_auto_fmt)) {
    auto_stats <- x_stats[is_auto_fmt]
    var_df <- .df_row[[.var]] # xxx this can be extended for the WHOLE data or single facets
    .formats[is_auto_fmt] <- lapply(names(auto_stats), format_auto, dt_var = var_df)
  }
  .formats
}

#' Horizontal waterfall plot
#'
#' @description `r lifecycle::badge("stable")`
#'
#' This basic waterfall plot visualizes a quantity `height` ordered by value with some markup.
#'
#' @param height (`numeric`)\cr vector containing values to be plotted as the waterfall bars.
#' @param id (`character`)\cr vector containing identifiers to use as the x-axis label for the waterfall bars.
#' @param col (`character`)\cr color(s).
#' @param col_var (`factor`, `character`, or `NULL`)\cr categorical variable for bar coloring. `NULL` by default.
#' @param xlab (`string`)\cr x label. Default is `"ID"`.
#' @param ylab (`string`)\cr y label. Default is `"Value"`.
#' @param title (`string`)\cr text to be displayed as plot title.
#' @param col_legend_title (`string`)\cr text to be displayed as legend title.
#'
#' @return A `ggplot` waterfall plot.
#'
#' @examples
#' library(dplyr)
#'
#' g_waterfall(height = c(3, 5, -1), id = letters[1:3])
#'
#' g_waterfall(
#'   height = c(3, 5, -1),
#'   id = letters[1:3],
#'   col_var = letters[1:3]
#' )
#'
#' adsl_f <- tern_ex_adsl %>%
#'   select(USUBJID, STUDYID, ARM, ARMCD, SEX)
#'
#' adrs_f <- tern_ex_adrs %>%
#'   filter(PARAMCD == "OVRINV") %>%
#'   mutate(pchg = rnorm(n(), 10, 50))
#'
#' adrs_f <- head(adrs_f, 30)
#' adrs_f <- adrs_f[!duplicated(adrs_f$USUBJID), ]
#' head(adrs_f)
#'
#' g_waterfall(
#'   height = adrs_f$pchg,
#'   id = adrs_f$USUBJID,
#'   col_var = adrs_f$AVALC
#' )
#'
#' g_waterfall(
#'   height = adrs_f$pchg,
#'   id = paste("asdfdsfdsfsd", adrs_f$USUBJID),
#'   col_var = adrs_f$SEX
#' )
#'
#' g_waterfall(
#'   height = adrs_f$pchg,
#'   id = paste("asdfdsfdsfsd", adrs_f$USUBJID),
#'   xlab = "ID",
#'   ylab = "Percentage Change",
#'   title = "Waterfall plot"
#' )
#'
#' @export
g_waterfall <- function(height,
                        id,
                        col_var = NULL,
                        col = getOption("ggplot2.discrete.colour"),
                        xlab = NULL,
                        ylab = NULL,
                        col_legend_title = NULL,
                        title = NULL) {
  if (!is.null(col_var)) {
    check_same_n(height = height, id = id, col_var = col_var)
  } else {
    check_same_n(height = height, id = id)
  }

  checkmate::assert_multi_class(col_var, c("character", "factor"), null.ok = TRUE)
  checkmate::assert_character(col, null.ok = TRUE)

  xlabel <- deparse(substitute(id))
  ylabel <- deparse(substitute(height))

  col_label <- if (!missing(col_var)) {
    deparse(substitute(col_var))
  }

  xlab <- if (is.null(xlab)) xlabel else xlab
  ylab <- if (is.null(ylab)) ylabel else ylab
  col_legend_title <- if (is.null(col_legend_title)) col_label else col_legend_title

  plot_data <- data.frame(
    height = height,
    id = as.character(id),
    col_var = if (is.null(col_var)) "x" else to_n(col_var, length(height)),
    stringsAsFactors = FALSE
  )

  plot_data_ord <- plot_data[order(plot_data$height, decreasing = TRUE), ]

  p <- ggplot2::ggplot(plot_data_ord, ggplot2::aes(x = factor(id, levels = id), y = height)) +
    ggplot2::geom_col() +
    ggplot2::geom_text(
      label = format(plot_data_ord$height, digits = 2),
      vjust = ifelse(plot_data_ord$height >= 0, -0.5, 1.5)
    ) +
    ggplot2::xlab(xlab) +
    ggplot2::ylab(ylab) +
    ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90, hjust = 0, vjust = .5))

  if (!is.null(col_var)) {
    p <- p +
      ggplot2::aes(fill = col_var) +
      ggplot2::labs(fill = col_legend_title) +
      ggplot2::theme(
        legend.position = "bottom",
        legend.background = ggplot2::element_blank(),
        legend.title = ggplot2::element_text(face = "bold"),
        legend.box.background = ggplot2::element_rect(colour = "black")
      )
  }

  if (!is.null(col)) {
    p <- p +
      ggplot2::scale_fill_manual(values = col)
  }

  if (!is.null(title)) {
    p <- p +
      ggplot2::labs(title = title) +
      ggplot2::theme(plot.title = ggplot2::element_text(face = "bold"))
  }

  p
}

#' Control functions for Kaplan-Meier plot annotation tables
#'
#' @description `r lifecycle::badge("stable")`
#'
#' Auxiliary functions for controlling arguments for formatting the annotation tables that can be added to plots
#' generated via [g_km()].
#'
#' @param x (`proportion`)\cr x-coordinate for center of annotation table.
#' @param y (`proportion`)\cr y-coordinate for center of annotation table.
#' @param w (`proportion`)\cr relative width of the annotation table.
#' @param h (`proportion`)\cr relative height of the annotation table.
#' @param fill (`flag` or `character`)\cr whether the annotation table should have a background fill color.
#'   Can also be a color code to use as the background fill color. If `TRUE`, color code defaults to `"#00000020"`.
#'
#' @return A list of components with the same names as the arguments.
#'
#' @seealso [g_km()]
#'
#' @name control_annot
NULL

#' @describeIn control_annot Control function for formatting the median survival time annotation table. This annotation
#'   table can be added in [g_km()] by setting `annot_surv_med=TRUE`, and can be configured using the
#'   `control_surv_med_annot()` function by setting it as the `control_annot_surv_med` argument.
#'
#' @examples
#' control_surv_med_annot()
#'
#' @export
control_surv_med_annot <- function(x = 0.8, y = 0.85, w = 0.32, h = 0.16, fill = TRUE) {
  assert_proportion_value(x)
  assert_proportion_value(y)
  assert_proportion_value(w)
  assert_proportion_value(h)

  list(x = x, y = y, w = w, h = h, fill = fill)
}

#' @describeIn control_annot Control function for formatting the Cox-PH annotation table. This annotation table can be
#'   added in [g_km()] by setting `annot_coxph=TRUE`, and can be configured using the `control_coxph_annot()` function
#'   by setting it as the `control_annot_coxph` argument.
#'
#' @param ref_lbls (`flag`)\cr whether the reference group should be explicitly printed in labels for the
#'   annotation table. If `FALSE` (default), only comparison groups will be printed in the table labels.
#'
#' @examples
#' control_coxph_annot()
#'
#' @export
control_coxph_annot <- function(x = 0.29, y = 0.51, w = 0.4, h = 0.125, fill = TRUE, ref_lbls = FALSE) {
  checkmate::assert_logical(ref_lbls, any.missing = FALSE)

  res <- c(control_surv_med_annot(x = x, y = y, w = w, h = h), list(ref_lbls = ref_lbls))
  res
}

#' Helper function to calculate 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
#' library(dplyr)
#' library(survival)
#'
#' data <- tern_ex_adtte %>%
#'   filter(PARAMCD == "OS") %>%
#'   survfit(formula = 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 for 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
#' library(dplyr)
#' library(survival)
#'
#' adtte <- tern_ex_adtte %>% filter(PARAMCD == "OS")
#' fit <- survfit(
#'   formula = 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 for generating a pairwise Cox-PH table
#'
#' @description `r lifecycle::badge("stable")`
#'
#' Create a `data.frame` of pairwise stratified or unstratified Cox-PH analysis results.
#'
#' @inheritParams g_km
#' @param annot_coxph_ref_lbls (`flag`)\cr whether the reference group should be explicitly printed in labels for the
#'   `annot_coxph` table. If `FALSE` (default), only comparison groups will be printed in `annot_coxph` table labels.
#'
#' @return A `data.frame` containing statistics `HR`, `XX% CI` (`XX` taken from `control_coxph_pw`),
#'   and `p-value (log-rank)`.
#'
#' @examples
#' 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,
                                 ref_group_coxph = NULL,
                                 control_coxph_pw = control_coxph(),
                                 annot_coxph_ref_lbls = FALSE) {
  if ("strat" %in% names(variables)) {
    warning(
      "Warning: the `strat` element name of the `variables` list argument to `h_tbl_coxph_pairwise() ",
      "was deprecated in tern 0.9.4.\n  ",
      "Please use the name `strata` instead of `strat` in the `variables` argument."
    )
    variables[["strata"]] <- variables[["strat"]]
  }

  assert_df_with_variables(df, variables)
  checkmate::assert_choice(ref_group_coxph, levels(df[[variables$arm]]), null.ok = TRUE)
  checkmate::assert_flag(annot_coxph_ref_lbls)

  arm <- variables$arm
  df[[arm]] <- factor(df[[arm]])

  ref_group <- if (!is.null(ref_group_coxph)) ref_group_coxph else levels(df[[variables$arm]])[1]
  comp_group <- setdiff(levels(df[[arm]]), ref_group)

  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,
      strata = variables$strata,
      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)
  if (annot_coxph_ref_lbls) names(results) <- paste(comp_group, "vs.", ref_group)

  do.call(rbind, results)
}

#' Helper function to tidy survival fit data
#'
#' @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
#' library(dplyr)
#' library(survival)
#'
#' # Test with multiple arms
#' tern_ex_adtte %>%
#'   filter(PARAMCD == "OS") %>%
#'   survfit(formula = Surv(AVAL, 1 - CNSR) ~ ARMCD, data = .) %>%
#'   h_data_plot()
#'
#' # Test with single arm
#' tern_ex_adtte %>%
#'   filter(PARAMCD == "OS", ARMCD == "ARM B") %>%
#'   survfit(formula = 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
}

## Deprecated Functions ----

#' Helper function to create a KM plot
#'
#' @description `r lifecycle::badge("deprecated")`
#'
#' 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(formula = 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()) {
  lifecycle::deprecate_warn(
    "0.9.4",
    "h_ggkm()",
    details = "`g_km` now generates `ggplot` objects. This function is no longer used within `tern`."
  )
  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("deprecated")`
#'
#' 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(formula = 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) {
  lifecycle::deprecate_warn(
    "0.9.4",
    "h_decompose_gg()",
    details = "`g_km` now generates `ggplot` objects. This function is no longer used within `tern`."
  )
  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 function to prepare a KM layout
#'
#' @description `r lifecycle::badge("deprecated")`
#'
#' 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(formula = 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) {
  lifecycle::deprecate_warn(
    "0.9.4",
    "h_km_layout()",
    details = "`g_km` now generates `ggplot` objects. This function is no longer used within `tern`."
  )
  txtlines <- levels(as.factor(data$strata))
  nlines <- nlevels(as.factor(data$strata))
  col_annot_width <- max(
    c(
      as.numeric(grid::convertX(g_el$yaxis$widths + g_el$ylab$widths, "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$heights + ylab$widths), "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$heights + ylab$widths), "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 function to create patient-at-risk grobs
#'
#' @description `r lifecycle::badge("deprecated")`
#'
#' 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(1)`)\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(formula = 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
#' # times (`xticks`).
#' annot_tbl <- summary(fit_km, times = 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) {
  lifecycle::deprecate_warn(
    "0.9.4",
    "h_grob_tbl_at_risk()",
    details = "`g_km` now generates `ggplot` objects. This function is no longer used within `tern`."
  )
  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 to create survival estimation grobs
#'
#' @description `r lifecycle::badge("deprecated")`
#'
#' 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 (`proportion`)\cr a value between 0 and 1 specifying x-location.
#' @param y (`proportion`)\cr a value between 0 and 1 specifying y-location.
#' @param width (`grid::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(formula = 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()) {
  lifecycle::deprecate_warn(
    "0.9.4",
    "h_grob_median_surv()",
    details = "`g_km` now generates `ggplot` objects. This function is no longer used within `tern`."
  )
  data <- h_tbl_median_surv(fit_km, armval = armval)

  width <- grid::convertUnit(grid::unit(as.numeric(width), grid::unitType(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 function to create grid object with y-axis annotation
#'
#' @description `r lifecycle::badge("deprecated")`
#'
#' 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(formula = 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) {
  lifecycle::deprecate_warn(
    "0.9.4",
    "h_grob_y_annot()",
    details = "`g_km` now generates `ggplot` objects. This function is no longer used within `tern`."
  )
  grid::gList(
    grid::gTree(
      vp = grid::viewport(
        width = grid::convertX(yaxis$widths + ylab$widths, "pt"),
        x = grid::unit(1, "npc"),
        just = "right"
      ),
      children = grid::gList(cbind(ylab, yaxis))
    )
  )
}

#' Helper function to create Cox-PH grobs
#'
#' @description `r lifecycle::badge("deprecated")`
#'
#' Grob of `rtable` output from [h_tbl_coxph_pairwise()]
#'
#' @inheritParams h_grob_median_surv
#' @param ... arguments to pass to [h_tbl_coxph_pairwise()].
#' @param x (`proportion`)\cr a value between 0 and 1 specifying x-location.
#' @param y (`proportion`)\cr a value between 0 and 1 specifying y-location.
#' @param width (`grid::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))
                         )) {
  lifecycle::deprecate_warn(
    "0.9.4",
    "h_grob_coxph()",
    details = "`g_km` now generates `ggplot` objects. This function is no longer used within `tern`."
  )
  data <- h_tbl_coxph_pairwise(...)

  width <- grid::convertUnit(grid::unit(as.numeric(width), grid::unitType(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
          )
        )
      )
    }
  )
}

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  checkmate::assert_string(x)

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

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

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

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

#' @describeIn h_logistic_regression Helper function to tabulate the interaction
#'   results of a logistic regression model. This basically is a wrapper for
#'   [h_or_interaction()] and [h_glm_simple_term_extract()] which puts the results
#'   in the right data frame format.
#'
#' @return A `data.frame` of tabulated interaction term results from a logistic regression model.
#'
#' @examples
#' h_glm_inter_term_extract("AGE", "ARMCD", mod2)
#'
#' @export
h_glm_inter_term_extract <- function(odds_ratio_var,
                                     interaction_var,
                                     fit_glm,
                                     ...) {
  # First obtain the main effects.
  main_stats <- h_glm_simple_term_extract(odds_ratio_var, fit_glm)
  main_stats$is_reference_summary <- FALSE
  main_stats$odds_ratio <- NA
  main_stats$lcl <- NA
  main_stats$ucl <- NA

  # Then we get the odds ratio estimates and put into df form.
  or_numbers <- h_or_interaction(odds_ratio_var, interaction_var, fit_glm, ...)
  is_num_or_var <- attr(fit_glm$terms, "dataClasses")[odds_ratio_var] == "numeric"

  if (is_num_or_var) {
    # Numeric OR variable case.
    references <- names(or_numbers)
    n_ref <- length(references)

    extract_from_list <- function(l, name, pos = 1) {
      unname(unlist(
        lapply(or_numbers, function(x) {
          x[[name]][pos]
        })
      ))
    }
    or_stats <- data.frame(
      variable = odds_ratio_var,
      variable_label = unname(formatters::var_labels(fit_glm$data[odds_ratio_var], fill = TRUE)),
      term = odds_ratio_var,
      term_label = unname(formatters::var_labels(fit_glm$data[odds_ratio_var], fill = TRUE)),
      interaction = interaction_var,
      interaction_label = unname(formatters::var_labels(fit_glm$data[interaction_var], fill = TRUE)),
      reference = references,
      reference_label = references,
      estimate = NA,
      std_error = NA,
      odds_ratio = extract_from_list(or_numbers, "or"),
      lcl = extract_from_list(or_numbers, "ci", pos = "lcl"),
      ucl = extract_from_list(or_numbers, "ci", pos = "ucl"),
      df = NA,
      pvalue = NA,
      is_variable_summary = FALSE,
      is_term_summary = FALSE,
      is_reference_summary = TRUE
    )
  } else {
    # Categorical OR variable case.
    references <- names(or_numbers[[1]])
    n_ref <- length(references)

    extract_from_list <- function(l, name, pos = 1) {
      unname(unlist(
        lapply(or_numbers, function(x) {
          lapply(x, function(y) y[[name]][pos])
        })
      ))
    }
    or_stats <- data.frame(
      variable = odds_ratio_var,
      variable_label = unname(formatters::var_labels(fit_glm$data[odds_ratio_var], fill = TRUE)),
      term = rep(names(or_numbers), each = n_ref),
      term_label = h_simple_term_labels(rep(names(or_numbers), each = n_ref), table(fit_glm$data[[odds_ratio_var]])),
      interaction = interaction_var,
      interaction_label = unname(formatters::var_labels(fit_glm$data[interaction_var], fill = TRUE)),
      reference = unlist(lapply(or_numbers, names)),
      reference_label = unlist(lapply(or_numbers, names)),
      estimate = NA,
      std_error = NA,
      odds_ratio = extract_from_list(or_numbers, "or"),
      lcl = extract_from_list(or_numbers, "ci", pos = "lcl"),
      ucl = extract_from_list(or_numbers, "ci", pos = "ucl"),
      df = NA,
      pvalue = NA,
      is_variable_summary = FALSE,
      is_term_summary = FALSE,
      is_reference_summary = TRUE
    )
  }

  df <- rbind(
    main_stats[, names(or_stats)],
    or_stats
  )
  df[order(-df$is_variable_summary, df$term, -df$is_term_summary, df$reference), ]
}

#' @describeIn h_logistic_regression Helper function to tabulate the results including
#'   odds ratios and confidence intervals of simple terms.
#'
#' @return Tabulated statistics for the given variable(s) from the logistic regression model.
#'
#' @examples
#' h_logistic_simple_terms("AGE", mod1)
#'
#' @export
h_logistic_simple_terms <- function(x, fit_glm, conf_level = 0.95) {
  checkmate::assert_multi_class(fit_glm, c("glm", "clogit"))
  if (inherits(fit_glm, "glm")) {
    checkmate::assert_set_equal(fit_glm$family$family, "binomial")
  }
  terms_name <- attr(stats::terms(fit_glm), "term.labels")
  xs_class <- attr(fit_glm$terms, "dataClasses")
  interaction <- terms_name[which(!terms_name %in% names(xs_class))]
  checkmate::assert_subset(x, terms_name)
  if (length(interaction) != 0) {
    # Make sure any item in x is not part of interaction term
    checkmate::assert_disjunct(x, unlist(strsplit(interaction, ":")))
  }
  x_stats <- lapply(x, h_glm_simple_term_extract, fit_glm)
  x_stats <- do.call(rbind, x_stats)
  q_norm <- stats::qnorm((1 + conf_level) / 2)
  x_stats$odds_ratio <- lapply(x_stats$estimate, exp)
  x_stats$lcl <- Map(function(or, se) exp(log(or) - q_norm * se), x_stats$odds_ratio, x_stats$std_error)
  x_stats$ucl <- Map(function(or, se) exp(log(or) + q_norm * se), x_stats$odds_ratio, x_stats$std_error)
  x_stats$ci <- Map(function(lcl, ucl) c(lcl, ucl), lcl = x_stats$lcl, ucl = x_stats$ucl)
  x_stats
}

#' @describeIn h_logistic_regression Helper function to tabulate the results including
#'   odds ratios and confidence intervals of interaction terms.
#'
#' @return Tabulated statistics for the given variable(s) from the logistic regression model.
#'
#' @examples
#' h_logistic_inter_terms(c("RACE", "AGE", "ARMCD", "AGE:ARMCD"), mod2)
#'
#' @export
h_logistic_inter_terms <- function(x,
                                   fit_glm,
                                   conf_level = 0.95,
                                   at = NULL) {
  # Find out the interaction variables and interaction term.
  inter_vars <- h_get_interaction_vars(fit_glm)
  checkmate::assert_vector(inter_vars, len = 2)


  inter_term_index <- intersect(grep(inter_vars[1], x), grep(inter_vars[2], x))
  inter_term <- x[inter_term_index]

  # For the non-interaction vars we need the standard stuff.
  normal_terms <- setdiff(x, union(inter_vars, inter_term))

  x_stats <- lapply(normal_terms, h_glm_simple_term_extract, fit_glm)
  x_stats <- do.call(rbind, x_stats)
  q_norm <- stats::qnorm((1 + conf_level) / 2)
  x_stats$odds_ratio <- lapply(x_stats$estimate, exp)
  x_stats$lcl <- Map(function(or, se) exp(log(or) - q_norm * se), x_stats$odds_ratio, x_stats$std_error)
  x_stats$ucl <- Map(function(or, se) exp(log(or) + q_norm * se), x_stats$odds_ratio, x_stats$std_error)
  normal_stats <- x_stats
  normal_stats$is_reference_summary <- FALSE

  # Now the interaction term itself.
  inter_term_stats <- h_glm_interaction_extract(inter_term, fit_glm)
  inter_term_stats$odds_ratio <- NA
  inter_term_stats$lcl <- NA
  inter_term_stats$ucl <- NA
  inter_term_stats$is_reference_summary <- FALSE

  is_intervar1_numeric <- attr(fit_glm$terms, "dataClasses")[inter_vars[1]] == "numeric"

  # Interaction stuff.
  inter_stats_one <- h_glm_inter_term_extract(
    inter_vars[1],
    inter_vars[2],
    fit_glm,
    conf_level = conf_level,
    at = `if`(is_intervar1_numeric, NULL, at)
  )
  inter_stats_two <- h_glm_inter_term_extract(
    inter_vars[2],
    inter_vars[1],
    fit_glm,
    conf_level = conf_level,
    at = `if`(is_intervar1_numeric, at, NULL)
  )

  # Now just combine everything in one data frame.
  col_names <- c(
    "variable",
    "variable_label",
    "term",
    "term_label",
    "interaction",
    "interaction_label",
    "reference",
    "reference_label",
    "estimate",
    "std_error",
    "df",
    "pvalue",
    "odds_ratio",
    "lcl",
    "ucl",
    "is_variable_summary",
    "is_term_summary",
    "is_reference_summary"
  )
  df <- rbind(
    inter_stats_one[, col_names],
    inter_stats_two[, col_names],
    inter_term_stats[, col_names]
  )
  if (length(normal_terms) > 0) {
    df <- rbind(
      normal_stats[, col_names],
      df
    )
  }
  df$ci <- combine_vectors(df$lcl, df$ucl)
  df
}

#' Create a STEP graph
#'
#' @description `r lifecycle::badge("stable")`
#'
#' Based on the STEP results, creates a `ggplot` graph showing the estimated HR or OR
#' along the continuous biomarker value subgroups.
#'
#' @param df (`tibble`)\cr result of [tidy.step()].
#' @param use_percentile (`flag`)\cr whether to use percentiles for the x axis or actual
#'   biomarker values.
#' @param est (named `list`)\cr `col` and `lty` settings for estimate line.
#' @param ci_ribbon (named `list` or `NULL`)\cr `fill` and `alpha` settings for the confidence interval
#'   ribbon area, or `NULL` to not plot a CI ribbon.
#' @param col (`character`)\cr color(s).
#'
#' @return A `ggplot` STEP graph.
#'
#' @seealso Custom tidy method [tidy.step()].
#'
#' @examples
#' library(survival)
#' lung$sex <- factor(lung$sex)
#'
#' # Survival example.
#' vars <- list(
#'   time = "time",
#'   event = "status",
#'   arm = "sex",
#'   biomarker = "age"
#' )
#'
#' step_matrix <- fit_survival_step(
#'   variables = vars,
#'   data = lung,
#'   control = c(control_coxph(), control_step(num_points = 10, degree = 2))
#' )
#' step_data <- broom::tidy(step_matrix)
#'
#' # Default plot.
#' g_step(step_data)
#'
#' # Add the reference 1 horizontal line.
#' library(ggplot2)
#' g_step(step_data) +
#'   ggplot2::geom_hline(ggplot2::aes(yintercept = 1), linetype = 2)
#'
#' # Use actual values instead of percentiles, different color for estimate and no CI,
#' # use log scale for y axis.
#' g_step(
#'   step_data,
#'   use_percentile = FALSE,
#'   est = list(col = "blue", lty = 1),
#'   ci_ribbon = NULL
#' ) + scale_y_log10()
#'
#' # Adding another curve based on additional column.
#' step_data$extra <- exp(step_data$`Percentile Center`)
#' g_step(step_data) +
#'   ggplot2::geom_line(ggplot2::aes(y = extra), linetype = 2, color = "green")
#'
#' # Response example.
#' vars <- list(
#'   response = "status",
#'   arm = "sex",
#'   biomarker = "age"
#' )
#'
#' step_matrix <- fit_rsp_step(
#'   variables = vars,
#'   data = lung,
#'   control = c(
#'     control_logistic(response_definition = "I(response == 2)"),
#'     control_step()
#'   )
#' )
#' step_data <- broom::tidy(step_matrix)
#' g_step(step_data)
#'
#' @export
g_step <- function(df,
                   use_percentile = "Percentile Center" %in% names(df),
                   est = list(col = "blue", lty = 1),
                   ci_ribbon = list(fill = getOption("ggplot2.discrete.colour")[1], alpha = 0.5),
                   col = getOption("ggplot2.discrete.colour")) {
  checkmate::assert_tibble(df)
  checkmate::assert_flag(use_percentile)
  checkmate::assert_character(col, null.ok = TRUE)
  checkmate::assert_list(est, names = "named")
  checkmate::assert_list(ci_ribbon, names = "named", null.ok = TRUE)

  x_var <- ifelse(use_percentile, "Percentile Center", "Interval Center")
  df$x <- df[[x_var]]
  attrs <- attributes(df)
  df$y <- df[[attrs$estimate]]

  # Set legend names. To be modified also at call level
  legend_names <- c("Estimate", "CI 95%")

  p <- ggplot2::ggplot(df, ggplot2::aes(x = .data[["x"]], y = .data[["y"]]))

  if (!is.null(col)) {
    p <- p +
      ggplot2::scale_color_manual(values = col)
  }

  if (!is.null(ci_ribbon)) {
    if (is.null(ci_ribbon$fill)) {
      ci_ribbon$fill <- "lightblue"
    }
    p <- p + ggplot2::geom_ribbon(
      ggplot2::aes(
        ymin = .data[["ci_lower"]], ymax = .data[["ci_upper"]],
        fill = legend_names[2]
      ),
      alpha = ci_ribbon$alpha
    ) +
      scale_fill_manual(
        name = "", values = c("CI 95%" = ci_ribbon$fill)
      )
  }
  suppressMessages(p <- p +
    ggplot2::geom_line(
      ggplot2::aes(y = .data[["y"]], color = legend_names[1]),
      linetype = est$lty
    ) +
    scale_colour_manual(
      name = "", values = c("Estimate" = "blue")
    ))

  p <- p + ggplot2::labs(x = attrs$biomarker, y = attrs$estimate)
  if (use_percentile) {
    p <- p + ggplot2::scale_x_continuous(labels = scales::percent)
  }
  p
}

#' Custom tidy method for STEP results
#'
#' @description `r lifecycle::badge("stable")`
#'
#' Tidy the STEP results into a `tibble` format ready for plotting.
#'
#' @param x (`matrix`)\cr results from [fit_survival_step()].
#' @param ... not used.
#'
#' @return A `tibble` with one row per STEP subgroup. The estimates and CIs are on the HR or OR scale,
#'   respectively. Additional attributes carry metadata also used for plotting.
#'
#' @seealso [g_step()] which consumes the result from this function.
#'
#' @method tidy step
#'
#' @examples
#' library(survival)
#' lung$sex <- factor(lung$sex)
#' vars <- list(
#'   time = "time",
#'   event = "status",
#'   arm = "sex",
#'   biomarker = "age"
#' )
#' step_matrix <- fit_survival_step(
#'   variables = vars,
#'   data = lung,
#'   control = c(control_coxph(), control_step(num_points = 10, degree = 2))
#' )
#' broom::tidy(step_matrix)
#'
#' @export
tidy.step <- function(x, ...) { # nolint
  checkmate::assert_class(x, "step")
  dat <- as.data.frame(x)
  nams <- names(dat)
  is_surv <- "loghr" %in% names(dat)
  est_var <- ifelse(is_surv, "loghr", "logor")
  new_est_var <- ifelse(is_surv, "Hazard Ratio", "Odds Ratio")
  new_y_vars <- c(new_est_var, c("ci_lower", "ci_upper"))
  names(dat)[match(est_var, nams)] <- new_est_var
  dat[, new_y_vars] <- exp(dat[, new_y_vars])
  any_is_na <- any(is.na(dat[, new_y_vars]))
  any_is_very_large <- any(abs(dat[, new_y_vars]) > 1e10, na.rm = TRUE)
  if (any_is_na) {
    warning(paste(
      "Missing values in the point estimate or CI columns,",
      "this will lead to holes in the `g_step()` plot"
    ))
  }
  if (any_is_very_large) {
    warning(paste(
      "Very large absolute values in the point estimate or CI columns,",
      "consider adding `scale_y_log10()` to the `g_step()` result for plotting"
    ))
  }
  if (any_is_na || any_is_very_large) {
    warning("Consider using larger `bandwidth`, less `num_points` in `control_step()` settings for fitting")
  }
  structure(
    tibble::as_tibble(dat),
    estimate = new_est_var,
    biomarker = attr(x, "variables")$biomarker,
    ci = f_conf_level(attr(x, "control")$conf_level)
  )
}

# summarize_glm_count ----------------------------------------------------------
#' Summarize Poisson negative binomial regression
#'
#' @description `r lifecycle::badge("experimental")`
#'
#' Summarize results of a Poisson negative binomial regression.
#' This can be used to analyze count and/or frequency data using a linear model.
#' It is specifically useful for analyzing count data (using the Poisson or Negative
#' Binomial distribution) that is result of a generalized linear model of one (e.g. arm) or more
#' covariates.
#'
#' @inheritParams h_glm_count
#' @inheritParams argument_convention
#' @param rate_mean_method (`character(1)`)\cr method used to estimate the mean odds ratio. Defaults to `emmeans`.
#'   see details for more information.
#' @param scale (`numeric(1)`)\cr linear scaling factor for rate and confidence intervals. Defaults to `1`.
#' @param .stats (`character`)\cr statistics to select for the table.
#'
#'   Options are: ``r shQuote(get_stats("summarize_glm_count"))``
#'
#' @details
#' `summarize_glm_count()` uses `s_glm_count()` to calculate the statistics for the table. This
#' analysis function uses [h_glm_count()] to estimate the GLM with [stats::glm()] for Poisson and Quasi-Poisson
#' distributions or [MASS::glm.nb()] for Negative Binomial distribution. All methods assume a
#' logarithmic link function.
#'
#' At this point, rates and confidence intervals are estimated from the model using
#' either [emmeans::emmeans()] when `rate_mean_method = "emmeans"` or [h_ppmeans()]
#' when `rate_mean_method = "ppmeans"`.
#'
#' If a reference group is specified while building the table with `split_cols_by(ref_group)`,
#' no rate ratio or `p-value` are calculated. Otherwise, we use [emmeans::contrast()] to
#' calculate the rate ratio and `p-value` for the reference group. Values are always estimated
#' with `method = "trt.vs.ctrl"` and `ref` equal to the first `arm` value.
#'
#' @name summarize_glm_count
NULL

#' @describeIn summarize_glm_count Layout-creating function which can take statistics function arguments
#'   and additional format arguments. This function is a wrapper for [rtables::analyze()].
#'
#' @return
#' * `summarize_glm_count()` 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_glm_count()` to the table layout.
#'
#' @examples
#' library(dplyr)
#'
#' anl <- tern_ex_adtte %>% filter(PARAMCD == "TNE")
#' anl$AVAL_f <- as.factor(anl$AVAL)
#'
#' lyt <- basic_table() %>%
#'   split_cols_by("ARM", ref_group = "B: Placebo") %>%
#'   add_colcounts() %>%
#'   analyze_vars(
#'     "AVAL_f",
#'     var_labels = "Number of exacerbations per patient",
#'     .stats = c("count_fraction"),
#'     .formats = c("count_fraction" = "xx (xx.xx%)"),
#'     .labels = c("Number of exacerbations per patient")
#'   ) %>%
#'   summarize_glm_count(
#'     vars = "AVAL",
#'     variables = list(arm = "ARM", offset = "lgTMATRSK", covariates = NULL),
#'     conf_level = 0.95,
#'     distribution = "poisson",
#'     rate_mean_method = "emmeans",
#'     var_labels = "Adjusted (P) exacerbation rate (per year)",
#'     table_names = "adjP",
#'     .stats = c("rate"),
#'     .labels = c(rate = "Rate")
#'   ) %>%
#'   summarize_glm_count(
#'     vars = "AVAL",
#'     variables = list(arm = "ARM", offset = "lgTMATRSK", covariates = c("REGION1")),
#'     conf_level = 0.95,
#'     distribution = "quasipoisson",
#'     rate_mean_method = "ppmeans",
#'     var_labels = "Adjusted (QP) exacerbation rate (per year)",
#'     table_names = "adjQP",
#'     .stats = c("rate", "rate_ci", "rate_ratio", "rate_ratio_ci", "pval"),
#'     .labels = c(
#'       rate = "Rate", rate_ci = "Rate CI", rate_ratio = "Rate Ratio",
#'       rate_ratio_ci = "Rate Ratio CI", pval = "p value"
#'     )
#'   )
#'
#' build_table(lyt = lyt, df = anl)
#'
#' @export
summarize_glm_count <- function(lyt,
                                vars,
                                variables,
                                distribution,
                                conf_level,
                                rate_mean_method = c("emmeans", "ppmeans")[1],
                                weights = stats::weights,
                                scale = 1,
                                var_labels,
                                na_str = default_na_str(),
                                nested = TRUE,
                                ...,
                                show_labels = "visible",
                                table_names = vars,
                                .stats = get_stats("summarize_glm_count"),
                                .formats = NULL,
                                .labels = NULL,
                                .indent_mods = c(
                                  "n" = 0L,
                                  "rate" = 0L,
                                  "rate_ci" = 1L,
                                  "rate_ratio" = 0L,
                                  "rate_ratio_ci" = 1L,
                                  "pval" = 1L
                                )) {
  checkmate::assert_choice(rate_mean_method, c("emmeans", "ppmeans"))

  extra_args <- list(
    variables = variables, distribution = distribution, conf_level = conf_level,
    rate_mean_method = rate_mean_method, weights = weights, scale = scale, ...
  )

  # Selecting parameters following the statistics
  .formats <- get_formats_from_stats(.stats, formats_in = .formats)
  .labels <- get_labels_from_stats(.stats, labels_in = .labels)
  .indent_mods <- get_indents_from_stats(.stats, indents_in = .indent_mods)

  afun <- make_afun(
    s_glm_count,
    .stats = .stats,
    .formats = .formats,
    .labels = .labels,
    .indent_mods = .indent_mods,
    .null_ref_cells = FALSE
  )

  analyze(
    lyt,
    vars,
    var_labels = var_labels,
    show_labels = show_labels,
    table_names = table_names,
    afun = afun,
    na_str = na_str,
    nested = nested,
    extra_args = extra_args
  )
}

#' @describeIn summarize_glm_count Statistics function that produces a named list of results
#'   of the investigated Poisson model.
#'
#' @return
#' * `s_glm_count()` returns a named `list` of 5 statistics:
#'   * `n`: Count of complete sample size for the group.
#'   * `rate`: Estimated event rate per follow-up time.
#'   * `rate_ci`: Confidence level for estimated rate per follow-up time.
#'   * `rate_ratio`: Ratio of event rates in each treatment arm to the reference arm.
#'   * `rate_ratio_ci`: Confidence level for the rate ratio.
#'   * `pval`: p-value.
#'
#' @keywords internal
s_glm_count <- function(df,
                        .var,
                        .df_row,
                        variables,
                        .ref_group,
                        .in_ref_col,
                        distribution,
                        conf_level,
                        rate_mean_method,
                        weights,
                        scale = 1) {
  arm <- variables$arm

  y <- df[[.var]]
  smry_level <- as.character(unique(df[[arm]]))

  # ensure there is only 1 value
  checkmate::assert_scalar(smry_level)

  results <- h_glm_count(
    .var = .var,
    .df_row = .df_row,
    variables = variables,
    distribution = distribution,
    weights
  )

  if (rate_mean_method == "emmeans") {
    emmeans_smry <- summary(results$emmeans_fit, level = conf_level)
  } else if (rate_mean_method == "ppmeans") {
    emmeans_smry <- h_ppmeans(results$glm_fit, .df_row, arm, conf_level)
  }

  emmeans_smry_level <- emmeans_smry[emmeans_smry[[arm]] == smry_level, ]

  # This happens if there is a reference col. No Ratio is calculated?
  if (.in_ref_col) {
    list(
      n = length(y[!is.na(y)]),
      rate = formatters::with_label(
        ifelse(distribution == "negbin", emmeans_smry_level$response * scale, emmeans_smry_level$rate * scale),
        "Adjusted Rate"
      ),
      rate_ci = formatters::with_label(
        c(emmeans_smry_level$asymp.LCL * scale, emmeans_smry_level$asymp.UCL * scale),
        f_conf_level(conf_level)
      ),
      rate_ratio = formatters::with_label(character(), "Adjusted Rate Ratio"),
      rate_ratio_ci = formatters::with_label(character(), f_conf_level(conf_level)),
      pval = formatters::with_label(character(), "p-value")
    )
  } else {
    emmeans_contrasts <- emmeans::contrast(
      results$emmeans_fit,
      method = "trt.vs.ctrl",
      ref = grep(
        as.character(unique(.ref_group[[arm]])),
        as.data.frame(results$emmeans_fit)[[arm]]
      )
    )

    contrasts_smry <- summary(
      emmeans_contrasts,
      infer = TRUE,
      adjust = "none"
    )

    smry_contrasts_level <- contrasts_smry[grepl(smry_level, contrasts_smry$contrast), ]

    list(
      n = length(y[!is.na(y)]),
      rate = formatters::with_label(
        ifelse(distribution == "negbin",
          emmeans_smry_level$response * scale,
          emmeans_smry_level$rate * scale
        ),
        "Adjusted Rate"
      ),
      rate_ci = formatters::with_label(
        c(emmeans_smry_level$asymp.LCL * scale, emmeans_smry_level$asymp.UCL * scale),
        f_conf_level(conf_level)
      ),
      rate_ratio = formatters::with_label(
        smry_contrasts_level$ratio,
        "Adjusted Rate Ratio"
      ),
      rate_ratio_ci = formatters::with_label(
        c(smry_contrasts_level$asymp.LCL, smry_contrasts_level$asymp.UCL),
        f_conf_level(conf_level)
      ),
      pval = formatters::with_label(
        smry_contrasts_level$p.value,
        "p-value"
      )
    )
  }
}
# h_glm_count ------------------------------------------------------------------
#' Helper functions for Poisson models
#'
#' @description `r lifecycle::badge("experimental")`
#'
#' Helper functions that returns the results of [stats::glm()] when Poisson or Quasi-Poisson
#' distributions are needed (see `family` parameter), or [MASS::glm.nb()] for Negative Binomial
#' distributions. Link function for the GLM is `log`.
#'
#' @inheritParams argument_convention
#'
#' @seealso [summarize_glm_count]
#'
#' @name h_glm_count
NULL

#' @describeIn h_glm_count Helper function to return the results of the
#'   selected model (Poisson, Quasi-Poisson, negative binomial).
#'
#' @param .df_row (`data.frame`)\cr dataset that includes all the variables that are called
#'   in `.var` and `variables`.
#' @param variables (named `list` of `string`)\cr list of additional analysis variables, with
#'   expected elements:
#'   * `arm` (`string`)\cr group variable, for which the covariate adjusted means of multiple
#'     groups will be summarized. Specifically, the first level of `arm` variable is taken as the
#'     reference group.
#'   * `covariates` (`character`)\cr a vector that can contain single variable names (such as
#'     `"X1"`), and/or interaction terms indicated by `"X1 * X2"`.
#'   * `offset` (`numeric`)\cr a numeric vector or scalar adding an offset.
#' @param distribution (`character`)\cr a character value specifying the distribution
#'   used in the regression (Poisson, Quasi-Poisson, negative binomial).
#' @param weights (`character`)\cr a character vector specifying weights used
#'   in averaging predictions. Number of weights must equal the number of levels included in the covariates.
#'   Weights option passed to [emmeans::emmeans()].
#'
#' @return
#' * `h_glm_count()` returns the results of the selected model.
#'
#' @keywords internal
h_glm_count <- function(.var,
                        .df_row,
                        variables,
                        distribution,
                        weights) {
  checkmate::assert_subset(distribution, c("poisson", "quasipoisson", "negbin"), empty.ok = FALSE)
  switch(distribution,
    poisson = h_glm_poisson(.var, .df_row, variables, weights),
    quasipoisson = h_glm_quasipoisson(.var, .df_row, variables, weights),
    negbin = h_glm_negbin(.var, .df_row, variables, weights)
  )
}

#' @describeIn h_glm_count Helper function to return results of a Poisson model.
#'
#' @return
#' * `h_glm_poisson()` returns the results of a Poisson model.
#'
#' @keywords internal
h_glm_poisson <- function(.var,
                          .df_row,
                          variables,
                          weights) {
  arm <- variables$arm
  covariates <- variables$covariates
  offset <- .df_row[[variables$offset]]

  formula <- stats::as.formula(paste0(
    .var, " ~ ",
    " + ",
    paste(covariates, collapse = " + "),
    " + ",
    arm
  ))

  glm_fit <- stats::glm(
    formula = formula,
    offset = offset,
    data = .df_row,
    family = stats::poisson(link = "log")
  )

  emmeans_fit <- emmeans::emmeans(
    glm_fit,
    specs = arm,
    data = .df_row,
    type = "response",
    offset = 0,
    weights = weights
  )

  list(
    glm_fit = glm_fit,
    emmeans_fit = emmeans_fit
  )
}

#' @describeIn h_glm_count Helper function to return results of a Quasi-Poisson model.
#'
#' @return
#' * `h_glm_quasipoisson()` returns the results of a Quasi-Poisson model.
#'
#' @keywords internal
h_glm_quasipoisson <- function(.var,
                               .df_row,
                               variables,
                               weights) {
  arm <- variables$arm
  covariates <- variables$covariates
  offset <- .df_row[[variables$offset]]

  formula <- stats::as.formula(paste0(
    .var, " ~ ",
    " + ",
    paste(covariates, collapse = " + "),
    " + ",
    arm
  ))

  glm_fit <- stats::glm(
    formula = formula,
    offset = offset,
    data = .df_row,
    family = stats::quasipoisson(link = "log")
  )

  emmeans_fit <- emmeans::emmeans(
    glm_fit,
    specs = arm,
    data = .df_row,
    type = "response",
    offset = 0,
    weights = weights
  )

  list(
    glm_fit = glm_fit,
    emmeans_fit = emmeans_fit
  )
}

#' @describeIn h_glm_count Helper function to return results of a negative binomial model.
#'
#' @return
#' * `h_glm_negbin()` returns the results of a negative binomial model.
#'
#' @keywords internal
h_glm_negbin <- function(.var,
                         .df_row,
                         variables,
                         weights) {
  arm <- variables$arm
  covariates <- variables$covariates

  formula <- stats::as.formula(paste0(
    .var, " ~ ",
    " + ",
    paste(covariates, collapse = " + "),
    " + ",
    arm
  ))

  glm_fit <- MASS::glm.nb(
    formula = formula,
    data = .df_row,
    link = "log"
  )

  emmeans_fit <- emmeans::emmeans(
    glm_fit,
    specs = arm,
    data = .df_row,
    type = "response",
    offset = 0,
    weights = weights
  )

  list(
    glm_fit = glm_fit,
    emmeans_fit = emmeans_fit
  )
}

# h_ppmeans --------------------------------------------------------------------
#' Function to return the estimated means using predicted probabilities
#'
#' @description
#' For each arm level, the predicted mean rate is calculated using the fitted model object, with `newdata`
#' set to the result of `stats::model.frame`, a reconstructed data or the original data, depending on the
#' object formula (coming from the fit). The confidence interval is derived using the `conf_level` parameter.
#'
#' @param obj (`glm.fit`)\cr fitted model object used to derive the mean rate estimates in each treatment arm.
#' @param .df_row (`data.frame`)\cr dataset that includes all the variables that are called in `.var` and `variables`.
#' @param arm (`string`)\cr group variable, for which the covariate adjusted means of multiple groups will be
#'   summarized. Specifically, the first level of `arm` variable is taken as the reference group.
#' @param conf_level (`proportion`)\cr value used to derive the confidence interval for the rate.
#'
#' @return
#' * `h_ppmeans()` returns the estimated means.
#'
#' @seealso [summarize_glm_count()].
#'
#' @export
h_ppmeans <- function(obj, .df_row, arm, conf_level) {
  alpha <- 1 - conf_level
  p <- 1 - alpha / 2

  arm_levels <- levels(.df_row[[arm]])

  out <- lapply(arm_levels, function(lev) {
    temp <- .df_row
    temp[[arm]] <- factor(lev, levels = arm_levels)

    mf <- stats::model.frame(obj$formula, data = temp)
    X <- stats::model.matrix(obj$formula, data = mf) # nolint

    rate <- stats::predict(obj, newdata = mf, type = "response")
    rate_hat <- mean(rate)

    zz <- colMeans(rate * X)
    se <- sqrt(as.numeric(t(zz) %*% stats::vcov(obj) %*% zz))
    rate_lwr <- rate_hat * exp(-stats::qnorm(p) * se / rate_hat)
    rate_upr <- rate_hat * exp(stats::qnorm(p) * se / rate_hat)

    c(rate_hat, rate_lwr, rate_upr)
  })

  names(out) <- arm_levels
  out <- do.call(rbind, out)
  if ("negbin" %in% class(obj)) {
    colnames(out) <- c("response", "asymp.LCL", "asymp.UCL")
  } else {
    colnames(out) <- c("rate", "asymp.LCL", "asymp.UCL")
  }
  out <- as.data.frame(out)
  out[[arm]] <- rownames(out)
  out
}

#' Compare variables between groups
#'
#' @description `r lifecycle::badge("stable")`
#'
#' The analyze function [compare_vars()] creates a layout element to summarize and compare one or more variables, using
#' the S3 generic function [s_summary()] to calculate a list of summary statistics. A list of all available statistics
#' for numeric variables can be viewed by running `get_stats("analyze_vars_numeric", add_pval = TRUE)` and for
#' non-numeric variables by running `get_stats("analyze_vars_counts", add_pval = TRUE)`. Use the `.stats` parameter to
#' specify the statistics to include in your output summary table.
#'
#' Prior to using this function in your table layout you must use [rtables::split_cols_by()] to create a column
#' split on the variable to be used in comparisons, and specify a reference group via the `ref_group` parameter.
#' Comparisons can be performed for each group (column) against the specified reference group by including the p-value
#' statistic.
#'
#' @inheritParams argument_convention
#' @param .stats (`character`)\cr statistics to select for the table.
#'
#'   Options for numeric variables are: ``r shQuote(get_stats("analyze_vars_numeric", add_pval = TRUE))``
#'
#'   Options for non-numeric variables are: ``r shQuote(get_stats("analyze_vars_counts", add_pval = TRUE))``
#'


1		#' Formatting functions
2		#'
3		#' See below for the list of formatting functions created in `tern` to work with `rtables`.
4		#'
5		#' Other available formats can be listed via [`formatters::list_valid_format_labels()`]. Additional
6		#' custom formats can be created via the [`formatters::sprintf_format()`] function.
7		#'
8		#' @family formatting functions
9		#' @name formatting_functions
10		NULL
11
12		#' Format fraction and percentage
13		#'
14		#' @description `r lifecycle::badge("stable")`
15		#'
16		#' Formats a fraction together with ratio in percent.
17		#'
18		#' @param x (named `integer`)\cr vector with elements `num` and `denom`.
19		#' @param ... not used. Required for `rtables` interface.
20		#'
21		#' @return A string in the format `num / denom (ratio %)`. If `num` is 0, the format is `num / denom`.
22		#'
23		#' @examples
24		#' format_fraction(x = c(num = 2L, denom = 3L))
25		#' format_fraction(x = c(num = 0L, denom = 3L))
26		#'
27		#' @family formatting functions
28		#' @export
29		format_fraction <- function(x, ...) {
30	4x	attr(x, "label") <- NULL
31
32	4x	checkmate::assert_vector(x)
33	4x	checkmate::assert_count(x["num"])
34	2x	checkmate::assert_count(x["denom"])
35
36	2x	result <- if (x["num"] == 0) {
37	1x	paste0(x["num"], "/", x["denom"])
38		} else {
39	1x	paste0(
40	1x	x["num"], "/", x["denom"],
41	1x	" (", round(x["num"] / x["denom"] * 100, 1), "%)"
42		)
43		}
44
45	2x	return(result)
46		}
47
48		#' Format fraction and percentage with fixed single decimal place
49		#'
50		#' @description `r lifecycle::badge("stable")`
51		#'
52		#' Formats a fraction together with ratio in percent with fixed single decimal place.
53		#' Includes trailing zero in case of whole number percentages to always keep one decimal place.
54		#'
55		#' @inheritParams format_fraction
56		#'
57		#' @return A string in the format `num / denom (ratio %)`. If `num` is 0, the format is `num / denom`.
58		#'
59		#' @examples
60		#' format_fraction_fixed_dp(x = c(num = 1L, denom = 2L))
61		#' format_fraction_fixed_dp(x = c(num = 1L, denom = 4L))
62		#' format_fraction_fixed_dp(x = c(num = 0L, denom = 3L))
63		#'
64		#' @family formatting functions
65		#' @export
66		format_fraction_fixed_dp <- function(x, ...) {
67	3x	attr(x, "label") <- NULL
68	3x	checkmate::assert_vector(x)
69	3x	checkmate::assert_count(x["num"])
70	3x	checkmate::assert_count(x["denom"])
71
72	3x	result <- if (x["num"] == 0) {
73	1x	paste0(x["num"], "/", x["denom"])
74		} else {
75	2x	paste0(
76	2x	x["num"], "/", x["denom"],
77	2x	" (", sprintf("%.1f", round(x["num"] / x["denom"] * 100, 1)), "%)"
78		)
79		}
80	3x	return(result)
81		}
82
83		#' Format count and fraction
84		#'
85		#' @description `r lifecycle::badge("stable")`
86		#'
87		#' Formats a count together with fraction with special consideration when count is `0`.
88		#'
89		#' @param x (`numeric(2)`)\cr vector of length 2 with count and fraction, respectively.
90		#' @param ... not used. Required for `rtables` interface.
91		#'
92		#' @return A string in the format `count (fraction %)`. If `count` is 0, the format is `0`.
93		#'
94		#' @examples
95		#' format_count_fraction(x = c(2, 0.6667))
96		#' format_count_fraction(x = c(0, 0))
97		#'
98		#' @family formatting functions
99		#' @export
100		format_count_fraction <- function(x, ...) {
101	3x	attr(x, "label") <- NULL
102
103	3x	if (any(is.na(x))) {
104	1x	return("NA")
105		}
106
107	2x	checkmate::assert_vector(x)
108	2x	checkmate::assert_integerish(x[1])
109	2x	assert_proportion_value(x[2], include_boundaries = TRUE)
110
111	2x	result <- if (x[1] == 0) {
112	1x	"0"
113		} else {
114	1x	paste0(x[1], " (", round(x[2] * 100, 1), "%)")
115		}
116
117	2x	return(result)
118		}
119
120		#' Format count and percentage with fixed single decimal place
121		#'
122		#' @description `r lifecycle::badge("experimental")`
123		#'
124		#' Formats a count together with fraction with special consideration when count is `0`.
125		#'
126		#' @inheritParams format_count_fraction
127		#'
128		#' @return A string in the format `count (fraction %)`. If `count` is 0, the format is `0`.
129		#'
130		#' @examples
131		#' format_count_fraction_fixed_dp(x = c(2, 0.6667))
132		#' format_count_fraction_fixed_dp(x = c(2, 0.5))
133		#' format_count_fraction_fixed_dp(x = c(0, 0))
134		#'
135		#' @family formatting functions
136		#' @export
137		format_count_fraction_fixed_dp <- function(x, ...) {
138	1093x	attr(x, "label") <- NULL
139
140	1093x	if (any(is.na(x))) {
141	!	return("NA")
142		}
143
144	1093x	checkmate::assert_vector(x)
145	1093x	checkmate::assert_integerish(x[1])
146	1093x	assert_proportion_value(x[2], include_boundaries = TRUE)
147
148	1093x	result <- if (x[1] == 0) {
149	177x	"0"
150	1093x	} else if (.is_equal_float(x[2], 1)) {
151	543x	sprintf("%d (100%%)", x[1])
152		} else {
153	373x	sprintf("%d (%.1f%%)", x[1], x[2] * 100)
154		}
155
156	1093x	return(result)
157		}
158
159		#' Format count and fraction with special case for count < 10
160		#'
161		#' @description `r lifecycle::badge("stable")`
162		#'
163		#' Formats a count together with fraction with special consideration when count is less than 10.
164		#'
165		#' @inheritParams format_count_fraction
166		#'
167		#' @return A string in the format `count (fraction %)`. If `count` is less than 10, only `count` is printed.
168		#'
169		#' @examples
170		#' format_count_fraction_lt10(x = c(275, 0.9673))
171		#' format_count_fraction_lt10(x = c(2, 0.6667))
172		#' format_count_fraction_lt10(x = c(9, 1))
173		#'
174		#' @family formatting functions
175		#' @export
176		format_count_fraction_lt10 <- function(x, ...) {
177	7x	attr(x, "label") <- NULL
178
179	7x	if (any(is.na(x))) {
180	1x	return("NA")
181		}
182
183	6x	checkmate::assert_vector(x)
184	6x	checkmate::assert_integerish(x[1])
185	6x	assert_proportion_value(x[2], include_boundaries = TRUE)
186
187	6x	result <- if (x[1] < 10) {
188	3x	paste0(x[1])
189		} else {
190	3x	paste0(x[1], " (", round(x[2] * 100, 1), "%)")
191		}
192
193	6x	return(result)
194		}
195
196		#' Format XX as a formatting function
197		#'
198		#' Translate a string where x and dots are interpreted as number place
199		#' holders, and others as formatting elements.
200		#'
201		#' @param str (`string`)\cr template.
202		#'
203		#' @return An `rtables` formatting function.
204		#'
205		#' @examples
206		#' test <- list(c(1.658, 0.5761), c(1e1, 785.6))
207		#'
208		#' z <- format_xx("xx (xx.x)")
209		#' sapply(test, z)
210		#'
211		#' z <- format_xx("xx.x - xx.x")
212		#' sapply(test, z)
213		#'
214		#' z <- format_xx("xx.x, incl. xx.x% NE")
215		#' sapply(test, z)
216		#'
217		#' @family formatting functions
218		#' @export
219		format_xx <- function(str) {
220		# Find position in the string.
221	1x	positions <- gregexpr(pattern = "x+\\.x+\|x+", text = str, perl = TRUE)
222	1x	x_positions <- regmatches(x = str, m = positions)[[1]]
223
224		# Roundings depends on the number of x behind [.].
225	1x	roundings <- lapply(
226	1x	X = x_positions,
227	1x	function(x) {
228	2x	y <- strsplit(split = "\\.", x = x)[[1]]
229	2x	rounding <- function(x) {
230	4x	round(x, digits = ifelse(length(y) > 1, nchar(y[2]), 0))
231		}
232	2x	return(rounding)
233		}
234		)
235
236	1x	rtable_format <- function(x, output) {
237	2x	values <- Map(y = x, fun = roundings, function(y, fun) fun(y))
238	2x	regmatches(x = str, m = positions)[[1]] <- values
239	2x	return(str)
240		}
241
242	1x	return(rtable_format)
243		}
244
245		#' Format numeric values by significant figures
246		#'
247		#' Format numeric values to print with a specified number of significant figures.
248		#'
249		#' @param sigfig (`integer(1)`)\cr number of significant figures to display.
250		#' @param format (`string`)\cr the format label (string) to apply when printing the value. Decimal
251		#' places in string are ignored in favor of formatting by significant figures. Formats options are:
252		#' `"xx"`, `"xx / xx"`, `"(xx, xx)"`, `"xx - xx"`, and `"xx (xx)"`.
253		#' @param num_fmt (`string`)\cr numeric format modifiers to apply to the value. Defaults to `"fg"` for
254		#' standard significant figures formatting - fixed (non-scientific notation) format (`"f"`)
255		#' and `sigfig` equal to number of significant figures instead of decimal places (`"g"`). See the
256		#' [formatC()] `format` argument for more options.
257		#'
258		#' @return An `rtables` formatting function.
259		#'
260		#' @examples
261		#' fmt_3sf <- format_sigfig(3)
262		#' fmt_3sf(1.658)
263		#' fmt_3sf(1e1)
264		#'
265		#' fmt_5sf <- format_sigfig(5)
266		#' fmt_5sf(0.57)
267		#' fmt_5sf(0.000025645)
268		#'
269		#' @family formatting functions
270		#' @export
271		format_sigfig <- function(sigfig, format = "xx", num_fmt = "fg") {
272	3x	checkmate::assert_integerish(sigfig)
273	3x	format <- gsub("xx\\.\|xx\\.x+", "xx", format)
274	3x	checkmate::assert_choice(format, c("xx", "xx / xx", "(xx, xx)", "xx - xx", "xx (xx)"))
275	3x	function(x, ...) {
276	!	if (!is.numeric(x)) stop("`format_sigfig` cannot be used for non-numeric values. Please choose another format.")
277	12x	num <- formatC(signif(x, digits = sigfig), digits = sigfig, format = num_fmt, flag = "#")
278	12x	num <- gsub("\\.$", "", num) # remove trailing "."
279
280	12x	format_value(num, format)
281		}
282		}
283
284		#' Format fraction with lower threshold
285		#'
286		#' @description `r lifecycle::badge("stable")`
287		#'
288		#' Formats a fraction when the second element of the input `x` is the fraction. It applies
289		#' a lower threshold, below which it is just stated that the fraction is smaller than that.
290		#'
291		#' @param threshold (`proportion`)\cr lower threshold.
292		#'
293		#' @return An `rtables` formatting function that takes numeric input `x` where the second
294		#' element is the fraction that is formatted. If the fraction is above or equal to the threshold,
295		#' then it is displayed in percentage. If it is positive but below the threshold, it returns,
296		#' e.g. "<1" if the threshold is `0.01`. If it is zero, then just "0" is returned.
297		#'
298		#' @examples
299		#' format_fun <- format_fraction_threshold(0.05)
300		#' format_fun(x = c(20, 0.1))
301		#' format_fun(x = c(2, 0.01))
302		#' format_fun(x = c(0, 0))
303		#'
304		#' @family formatting functions
305		#' @export
306		format_fraction_threshold <- function(threshold) {
307	1x	assert_proportion_value(threshold)
308	1x	string_below_threshold <- paste0("<", round(threshold * 100))
309	1x	function(x, ...) {
310	3x	assert_proportion_value(x[2], include_boundaries = TRUE)
311	3x	ifelse(
312	3x	x[2] > 0.01,
313	3x	round(x[2] * 100),
314	3x	ifelse(
315	3x	x[2] == 0,
316	3x	"0",
317	3x	string_below_threshold
318		)
319		)
320		}
321		}
322
323		#' Format extreme values
324		#'
325		#' @description `r lifecycle::badge("stable")`
326		#'
327		#' `rtables` formatting functions that handle extreme values.
328		#'
329		#' @param digits (`integer(1)`)\cr number of decimal places to display.
330		#'
331		#' @details For each input, apply a format to the specified number of `digits`. If the value is
332		#' below a threshold, it returns "<0.01" e.g. if the number of `digits` is 2. If the value is
333		#' above a threshold, it returns ">999.99" e.g. if the number of `digits` is 2.
334		#' If it is zero, then returns "0.00".
335		#'
336		#' @family formatting functions
337		#' @name extreme_format
338		NULL
339
340		#' @describeIn extreme_format Internal helper function to calculate the threshold and create formatted strings
341		#' used in Formatting Functions. Returns a list with elements `threshold` and `format_string`.
342		#'
343		#' @return
344		#' * `h_get_format_threshold()` returns a `list` of 2 elements: `threshold`, with `low` and `high` thresholds,
345		#' and `format_string`, with thresholds formatted as strings.
346		#'
347		#' @examples
348		#' h_get_format_threshold(2L)
349		#'
350		#' @export
351		h_get_format_threshold <- function(digits = 2L) {
352	2113x	checkmate::assert_integerish(digits)
353
354	2113x	low_threshold <- 1 / (10 ^ digits) # styler: off
355	2113x	high_threshold <- 1000 - (1 / (10 ^ digits)) # styler: off
356
357	2113x	string_below_threshold <- paste0("<", low_threshold)
358	2113x	string_above_threshold <- paste0(">", high_threshold)
359
360	2113x	list(
361	2113x	"threshold" = c(low = low_threshold, high = high_threshold),
362	2113x	"format_string" = c(low = string_below_threshold, high = string_above_threshold)
363		)
364		}
365
366		#' @describeIn extreme_format Internal helper function to apply a threshold format to a value.
367		#' Creates a formatted string to be used in Formatting Functions.
368		#'
369		#' @param x (`numeric(1)`)\cr value to format.
370		#'
371		#' @return
372		#' * `h_format_threshold()` returns the given value, or if the value is not within the digit threshold the relation
373		#' of the given value to the digit threshold, as a formatted string.
374		#'
375		#' @examples
376		#' h_format_threshold(0.001)
377		#' h_format_threshold(1000)
378		#'
379		#' @export
380		h_format_threshold <- function(x, digits = 2L) {
381	2115x	if (is.na(x)) {
382	4x	return(x)
383		}
384
385	2111x	checkmate::assert_numeric(x, lower = 0)
386
387	2111x	l_fmt <- h_get_format_threshold(digits)
388
389	2111x	result <- if (x < l_fmt$threshold["low"] && 0 < x) {
390	44x	l_fmt$format_string["low"]
391	2111x	} else if (x > l_fmt$threshold["high"]) {
392	99x	l_fmt$format_string["high"]
393		} else {
394	1968x	sprintf(fmt = paste0("%.", digits, "f"), x)
395		}
396
397	2111x	unname(result)
398		}
399
400		#' Format a single extreme value
401		#'
402		#' @description `r lifecycle::badge("stable")`
403		#'
404		#' Create a formatting function for a single extreme value.
405		#'
406		#' @inheritParams extreme_format
407		#'
408		#' @return An `rtables` formatting function that uses threshold `digits` to return a formatted extreme value.
409		#'
410		#' @examples
411		#' format_fun <- format_extreme_values(2L)
412		#' format_fun(x = 0.127)
413		#' format_fun(x = Inf)
414		#' format_fun(x = 0)
415		#' format_fun(x = 0.009)
416		#'
417		#' @family formatting functions
418		#' @export
419		format_extreme_values <- function(digits = 2L) {
420	63x	function(x, ...) {
421	657x	checkmate::assert_scalar(x, na.ok = TRUE)
422
423	657x	h_format_threshold(x = x, digits = digits)
424		}
425		}
426
427		#' Format extreme values part of a confidence interval
428		#'
429		#' @description `r lifecycle::badge("stable")`
430		#'
431		#' Formatting Function for extreme values part of a confidence interval. Values
432		#' are formatted as e.g. "(xx.xx, xx.xx)" if the number of `digits` is 2.
433		#'
434		#' @inheritParams extreme_format
435		#'
436		#' @return An `rtables` formatting function that uses threshold `digits` to return a formatted extreme
437		#' values confidence interval.
438		#'
439		#' @examples
440		#' format_fun <- format_extreme_values_ci(2L)
441		#' format_fun(x = c(0.127, Inf))
442		#' format_fun(x = c(0, 0.009))
443		#'
444		#' @family formatting functions
445		#' @export
446		format_extreme_values_ci <- function(digits = 2L) {
447	71x	function(x, ...) {
448	726x	checkmate::assert_vector(x, len = 2)
449	726x	l_result <- h_format_threshold(x = x[1], digits = digits)
450	726x	h_result <- h_format_threshold(x = x[2], digits = digits)
451
452	726x	paste0("(", l_result, ", ", h_result, ")")
453		}
454		}
455
456		#' Format automatically using data significant digits
457		#'
458		#' @description `r lifecycle::badge("stable")`
459		#'
460		#' Formatting function for the majority of default methods used in [analyze_vars()].
461		#' For non-derived values, the significant digits of data is used (e.g. range), while derived
462		#' values have one more digits (measure of location and dispersion like mean, standard deviation).
463		#' This function can be called internally with "auto" like, for example,
464		#' `.formats = c("mean" = "auto")`. See details to see how this works with the inner function.
465		#'
466		#' @param dt_var (`numeric`)\cr variable data the statistics were calculated from. Used only to
467		#' find significant digits. In [analyze_vars] this comes from `.df_row` (see
468		#' [rtables::additional_fun_params]), and it is the row data after the above row splits. No
469		#' column split is considered.
470		#' @param x_stat (`string`)\cr string indicating the current statistical method used.
471		#'
472		#' @return A string that `rtables` prints in a table cell.
473		#'
474		#' @details
475		#' The internal function is needed to work with `rtables` default structure for
476		#' format functions, i.e. `function(x, ...)`, where is x are results from statistical evaluation.
477		#' It can be more than one element (e.g. for `.stats = "mean_sd"`).
478		#'
479		#' @examples
480		#' x_todo <- c(0.001, 0.2, 0.0011000, 3, 4)
481		#' res <- c(mean(x_todo[1:3]), sd(x_todo[1:3]))
482		#'
483		#' # x is the result coming into the formatting function -> res!!
484		#' format_auto(dt_var = x_todo, x_stat = "mean_sd")(x = res)
485		#' format_auto(x_todo, "range")(x = range(x_todo))
486		#' no_sc_x <- c(0.0000001, 1)
487		#' format_auto(no_sc_x, "range")(x = no_sc_x)
488		#'
489		#' @family formatting functions
490		#' @export
491		format_auto <- function(dt_var, x_stat) {
492	10x	function(x = "", ...) {
493	18x	checkmate::assert_numeric(x, min.len = 1)
494	18x	checkmate::assert_numeric(dt_var, min.len = 1)
495		# Defaults - they may be a param in the future
496	18x	der_stats <- c(
497	18x	"mean", "sd", "se", "median", "geom_mean", "quantiles", "iqr",
498	18x	"mean_sd", "mean_se", "mean_se", "mean_ci", "mean_sei", "mean_sdi",
499	18x	"median_ci"
500		)
501	18x	nonder_stats <- c("n", "range", "min", "max")
502
503		# Safenet for miss-modifications
504	18x	stopifnot(length(intersect(der_stats, nonder_stats)) == 0) # nolint
505	18x	checkmate::assert_choice(x_stat, c(der_stats, nonder_stats))
506
507		# Finds the max number of digits in data
508	18x	detect_dig <- vapply(dt_var, count_decimalplaces, FUN.VALUE = numeric(1)) %>%
509	18x	max()
510
511	18x	if (x_stat %in% der_stats) {
512	8x	detect_dig <- detect_dig + 1
513		}
514
515		# Render input
516	18x	str_vals <- formatC(x, digits = detect_dig, format = "f")
517	18x	def_fmt <- get_formats_from_stats(x_stat)[[x_stat]]
518	18x	str_fmt <- str_extract(def_fmt, invert = FALSE)[[1]]
519	18x	if (length(str_fmt) != length(str_vals)) {
520	2x	stop(
521	2x	"Number of inserted values as result (", length(str_vals),
522	2x	") is not the same as there should be in the default tern formats for ",
523	2x	x_stat, " (-> ", def_fmt, " needs ", length(str_fmt), " values). ",
524	2x	"See tern_default_formats to check all of them."
525		)
526		}
527
528		# Squashing them together
529	16x	inv_str_fmt <- str_extract(def_fmt, invert = TRUE)[[1]]
530	16x	stopifnot(length(inv_str_fmt) == length(str_vals) + 1) # nolint
531
532	16x	out <- vector("character", length = length(inv_str_fmt) + length(str_vals))
533	16x	is_even <- seq_along(out) %% 2 == 0
534	16x	out[is_even] <- str_vals
535	16x	out[!is_even] <- inv_str_fmt
536
537	16x	return(paste0(out, collapse = ""))
538		}
539		}
540
541		# Utility function that could be useful in general
542		str_extract <- function(string, pattern = "xx\|xx\\.\|xx\\.x+", invert = FALSE) {
543	34x	regmatches(string, gregexpr(pattern, string), invert = invert)
544		}
545
546		# Helper function
547		count_decimalplaces <- function(dec) {
548	161x	if (is.na(dec)) {
549	6x	return(0)
550	155x	} else if (abs(dec - round(dec)) > .Machine$double.eps^0.5) { # For precision
551	122x	nchar(strsplit(format(dec, scientific = FALSE, trim = FALSE), ".", fixed = TRUE)[[1]][[2]])
552		} else {
553	33x	return(0)
554		}
555		}
556
557		#' Apply automatic formatting
558		#'
559		#' Checks if any of the listed formats in `.formats` are `"auto"`, and replaces `"auto"` with
560		#' the correct implementation of `format_auto` for the given statistics, data, and variable.
561		#'
562		#' @inheritParams argument_convention
563		#' @param x_stats (named `list`)\cr a named list of statistics where each element corresponds
564		#' to an element in `.formats`, with matching names.
565		#'
566		#' @keywords internal
567		apply_auto_formatting <- function(.formats, x_stats, .df_row, .var) {
568	476x	is_auto_fmt <- vapply(.formats, function(ii) is.character(ii) && ii == "auto", logical(1))
569	476x	if (any(is_auto_fmt)) {
570	3x	auto_stats <- x_stats[is_auto_fmt]
571	3x	var_df <- .df_row[[.var]] # xxx this can be extended for the WHOLE data or single facets
572	3x	.formats[is_auto_fmt] <- lapply(names(auto_stats), format_auto, dt_var = var_df)
573		}
574	476x	.formats
575		}

1		#' Horizontal waterfall plot
2		#'
3		#' @description `r lifecycle::badge("stable")`
4		#'
5		#' This basic waterfall plot visualizes a quantity `height` ordered by value with some markup.
6		#'
7		#' @param height (`numeric`)\cr vector containing values to be plotted as the waterfall bars.
8		#' @param id (`character`)\cr vector containing identifiers to use as the x-axis label for the waterfall bars.
9		#' @param col (`character`)\cr color(s).
10		#' @param col_var (`factor`, `character`, or `NULL`)\cr categorical variable for bar coloring. `NULL` by default.
11		#' @param xlab (`string`)\cr x label. Default is `"ID"`.
12		#' @param ylab (`string`)\cr y label. Default is `"Value"`.
13		#' @param title (`string`)\cr text to be displayed as plot title.
14		#' @param col_legend_title (`string`)\cr text to be displayed as legend title.
15		#'
16		#' @return A `ggplot` waterfall plot.
17		#'
18		#' @examples
19		#' library(dplyr)
20		#'
21		#' g_waterfall(height = c(3, 5, -1), id = letters[1:3])
22		#'
23		#' g_waterfall(
24		#' height = c(3, 5, -1),
25		#' id = letters[1:3],
26		#' col_var = letters[1:3]
27		#' )
28		#'
29		#' adsl_f <- tern_ex_adsl %>%
30		#' select(USUBJID, STUDYID, ARM, ARMCD, SEX)
31		#'
32		#' adrs_f <- tern_ex_adrs %>%
33		#' filter(PARAMCD == "OVRINV") %>%
34		#' mutate(pchg = rnorm(n(), 10, 50))
35		#'
36		#' adrs_f <- head(adrs_f, 30)
37		#' adrs_f <- adrs_f[!duplicated(adrs_f$USUBJID), ]
38		#' head(adrs_f)
39		#'
40		#' g_waterfall(
41		#' height = adrs_f$pchg,
42		#' id = adrs_f$USUBJID,
43		#' col_var = adrs_f$AVALC
44		#' )
45		#'
46		#' g_waterfall(
47		#' height = adrs_f$pchg,
48		#' id = paste("asdfdsfdsfsd", adrs_f$USUBJID),
49		#' col_var = adrs_f$SEX
50		#' )
51		#'
52		#' g_waterfall(
53		#' height = adrs_f$pchg,
54		#' id = paste("asdfdsfdsfsd", adrs_f$USUBJID),
55		#' xlab = "ID",
56		#' ylab = "Percentage Change",
57		#' title = "Waterfall plot"
58		#' )
59		#'
60		#' @export
61		g_waterfall <- function(height,
62		id,
63		col_var = NULL,
64		col = getOption("ggplot2.discrete.colour"),
65		xlab = NULL,
66		ylab = NULL,
67		col_legend_title = NULL,
68		title = NULL) {
69	2x	if (!is.null(col_var)) {
70	1x	check_same_n(height = height, id = id, col_var = col_var)
71		} else {
72	1x	check_same_n(height = height, id = id)
73		}
74
75	2x	checkmate::assert_multi_class(col_var, c("character", "factor"), null.ok = TRUE)
76	2x	checkmate::assert_character(col, null.ok = TRUE)
77
78	2x	xlabel <- deparse(substitute(id))
79	2x	ylabel <- deparse(substitute(height))
80
81	2x	col_label <- if (!missing(col_var)) {
82	1x	deparse(substitute(col_var))
83		}
84
85	2x	xlab <- if (is.null(xlab)) xlabel else xlab
86	2x	ylab <- if (is.null(ylab)) ylabel else ylab
87	2x	col_legend_title <- if (is.null(col_legend_title)) col_label else col_legend_title
88
89	2x	plot_data <- data.frame(
90	2x	height = height,
91	2x	id = as.character(id),
92	2x	col_var = if (is.null(col_var)) "x" else to_n(col_var, length(height)),
93	2x	stringsAsFactors = FALSE
94		)
95
96	2x	plot_data_ord <- plot_data[order(plot_data$height, decreasing = TRUE), ]
97
98	2x	p <- ggplot2::ggplot(plot_data_ord, ggplot2::aes(x = factor(id, levels = id), y = height)) +
99	2x	ggplot2::geom_col() +
100	2x	ggplot2::geom_text(
101	2x	label = format(plot_data_ord$height, digits = 2),
102	2x	vjust = ifelse(plot_data_ord$height >= 0, -0.5, 1.5)
103		) +
104	2x	ggplot2::xlab(xlab) +
105	2x	ggplot2::ylab(ylab) +
106	2x	ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90, hjust = 0, vjust = .5))
107
108	2x	if (!is.null(col_var)) {
109	1x	p <- p +
110	1x	ggplot2::aes(fill = col_var) +
111	1x	ggplot2::labs(fill = col_legend_title) +
112	1x	ggplot2::theme(
113	1x	legend.position = "bottom",
114	1x	legend.background = ggplot2::element_blank(),
115	1x	legend.title = ggplot2::element_text(face = "bold"),
116	1x	legend.box.background = ggplot2::element_rect(colour = "black")
117		)
118		}
119
120	2x	if (!is.null(col)) {
121	1x	p <- p +
122	1x	ggplot2::scale_fill_manual(values = col)
123		}
124
125	2x	if (!is.null(title)) {
126	1x	p <- p +
127	1x	ggplot2::labs(title = title) +
128	1x	ggplot2::theme(plot.title = ggplot2::element_text(face = "bold"))
129		}
130
131	2x	p
132		}

1		#' Create a STEP graph
2		#'
3		#' @description `r lifecycle::badge("stable")`
4		#'
5		#' Based on the STEP results, creates a `ggplot` graph showing the estimated HR or OR
6		#' along the continuous biomarker value subgroups.
7		#'
8		#' @param df (`tibble`)\cr result of [tidy.step()].
9		#' @param use_percentile (`flag`)\cr whether to use percentiles for the x axis or actual
10		#' biomarker values.
11		#' @param est (named `list`)\cr `col` and `lty` settings for estimate line.
12		#' @param ci_ribbon (named `list` or `NULL`)\cr `fill` and `alpha` settings for the confidence interval
13		#' ribbon area, or `NULL` to not plot a CI ribbon.
14		#' @param col (`character`)\cr color(s).
15		#'
16		#' @return A `ggplot` STEP graph.
17		#'
18		#' @seealso Custom tidy method [tidy.step()].
19		#'
20		#' @examples
21		#' library(survival)
22		#' lung$sex <- factor(lung$sex)
23		#'
24		#' # Survival example.
25		#' vars <- list(
26		#' time = "time",
27		#' event = "status",
28		#' arm = "sex",
29		#' biomarker = "age"
30		#' )
31		#'
32		#' step_matrix <- fit_survival_step(
33		#' variables = vars,
34		#' data = lung,
35		#' control = c(control_coxph(), control_step(num_points = 10, degree = 2))
36		#' )
37		#' step_data <- broom::tidy(step_matrix)
38		#'
39		#' # Default plot.
40		#' g_step(step_data)
41		#'
42		#' # Add the reference 1 horizontal line.
43		#' library(ggplot2)
44		#' g_step(step_data) +
45		#' ggplot2::geom_hline(ggplot2::aes(yintercept = 1), linetype = 2)
46		#'
47		#' # Use actual values instead of percentiles, different color for estimate and no CI,
48		#' # use log scale for y axis.
49		#' g_step(
50		#' step_data,
51		#' use_percentile = FALSE,
52		#' est = list(col = "blue", lty = 1),
53		#' ci_ribbon = NULL
54		#' ) + scale_y_log10()
55		#'
56		#' # Adding another curve based on additional column.
57		#' step_data$extra <- exp(step_data$`Percentile Center`)
58		#' g_step(step_data) +
59		#' ggplot2::geom_line(ggplot2::aes(y = extra), linetype = 2, color = "green")
60		#'
61		#' # Response example.
62		#' vars <- list(
63		#' response = "status",
64		#' arm = "sex",
65		#' biomarker = "age"
66		#' )
67		#'
68		#' step_matrix <- fit_rsp_step(
69		#' variables = vars,
70		#' data = lung,
71		#' control = c(
72		#' control_logistic(response_definition = "I(response == 2)"),
73		#' control_step()
74		#' )
75		#' )
76		#' step_data <- broom::tidy(step_matrix)
77		#' g_step(step_data)
78		#'
79		#' @export
80		g_step <- function(df,
81		use_percentile = "Percentile Center" %in% names(df),
82		est = list(col = "blue", lty = 1),
83		ci_ribbon = list(fill = getOption("ggplot2.discrete.colour")[1], alpha = 0.5),
84		col = getOption("ggplot2.discrete.colour")) {
85	2x	checkmate::assert_tibble(df)
86	2x	checkmate::assert_flag(use_percentile)
87	2x	checkmate::assert_character(col, null.ok = TRUE)
88	2x	checkmate::assert_list(est, names = "named")
89	2x	checkmate::assert_list(ci_ribbon, names = "named", null.ok = TRUE)
90
91	2x	x_var <- ifelse(use_percentile, "Percentile Center", "Interval Center")
92	2x	df$x <- df[[x_var]]
93	2x	attrs <- attributes(df)
94	2x	df$y <- df[[attrs$estimate]]
95
96		# Set legend names. To be modified also at call level
97	2x	legend_names <- c("Estimate", "CI 95%")
98
99	2x	p <- ggplot2::ggplot(df, ggplot2::aes(x = .data[["x"]], y = .data[["y"]]))
100
101	2x	if (!is.null(col)) {
102	2x	p <- p +
103	2x	ggplot2::scale_color_manual(values = col)
104		}
105
106	2x	if (!is.null(ci_ribbon)) {
107	1x	if (is.null(ci_ribbon$fill)) {
108	!	ci_ribbon$fill <- "lightblue"
109		}
110	1x	p <- p + ggplot2::geom_ribbon(
111	1x	ggplot2::aes(
112	1x	ymin = .data[["ci_lower"]], ymax = .data[["ci_upper"]],
113	1x	fill = legend_names[2]
114		),
115	1x	alpha = ci_ribbon$alpha
116		) +
117	1x	scale_fill_manual(
118	1x	name = "", values = c("CI 95%" = ci_ribbon$fill)
119		)
120		}
121	2x	suppressMessages(p <- p +
122	2x	ggplot2::geom_line(
123	2x	ggplot2::aes(y = .data[["y"]], color = legend_names[1]),
124	2x	linetype = est$lty
125		) +
126	2x	scale_colour_manual(
127	2x	name = "", values = c("Estimate" = "blue")
128		))
129
130	2x	p <- p + ggplot2::labs(x = attrs$biomarker, y = attrs$estimate)
131	2x	if (use_percentile) {
132	1x	p <- p + ggplot2::scale_x_continuous(labels = scales::percent)
133		}
134	2x	p
135		}
136
137		#' Custom tidy method for STEP results
138		#'
139		#' @description `r lifecycle::badge("stable")`
140		#'
141		#' Tidy the STEP results into a `tibble` format ready for plotting.
142		#'
143		#' @param x (`matrix`)\cr results from [fit_survival_step()].
144		#' @param ... not used.
145		#'
146		#' @return A `tibble` with one row per STEP subgroup. The estimates and CIs are on the HR or OR scale,
147		#' respectively. Additional attributes carry metadata also used for plotting.
148		#'
149		#' @seealso [g_step()] which consumes the result from this function.
150		#'
151		#' @method tidy step
152		#'
153		#' @examples
154		#' library(survival)
155		#' lung$sex <- factor(lung$sex)
156		#' vars <- list(
157		#' time = "time",
158		#' event = "status",
159		#' arm = "sex",
160		#' biomarker = "age"
161		#' )
162		#' step_matrix <- fit_survival_step(
163		#' variables = vars,
164		#' data = lung,
165		#' control = c(control_coxph(), control_step(num_points = 10, degree = 2))
166		#' )
167		#' broom::tidy(step_matrix)
168		#'
169		#' @export
170		tidy.step <- function(x, ...) { # nolint
171	7x	checkmate::assert_class(x, "step")
172	7x	dat <- as.data.frame(x)
173	7x	nams <- names(dat)
174	7x	is_surv <- "loghr" %in% names(dat)
175	7x	est_var <- ifelse(is_surv, "loghr", "logor")
176	7x	new_est_var <- ifelse(is_surv, "Hazard Ratio", "Odds Ratio")
177	7x	new_y_vars <- c(new_est_var, c("ci_lower", "ci_upper"))
178	7x	names(dat)[match(est_var, nams)] <- new_est_var
179	7x	dat[, new_y_vars] <- exp(dat[, new_y_vars])
180	7x	any_is_na <- any(is.na(dat[, new_y_vars]))
181	7x	any_is_very_large <- any(abs(dat[, new_y_vars]) > 1e10, na.rm = TRUE)
182	7x	if (any_is_na) {
183	2x	warning(paste(
184	2x	"Missing values in the point estimate or CI columns,",
185	2x	"this will lead to holes in the `g_step()` plot"
186		))
187		}
188	7x	if (any_is_very_large) {
189	2x	warning(paste(
190	2x	"Very large absolute values in the point estimate or CI columns,",
191	2x	"consider adding `scale_y_log10()` to the `g_step()` result for plotting"
192		))
193		}
194	7x	if (any_is_na \|\| any_is_very_large) {
195	4x	warning("Consider using larger `bandwidth`, less `num_points` in `control_step()` settings for fitting")
196		}
197	7x	structure(
198	7x	tibble::as_tibble(dat),
199	7x	estimate = new_est_var,
200	7x	biomarker = attr(x, "variables")$biomarker,
201	7x	ci = f_conf_level(attr(x, "control")$conf_level)
202		)
203		}

1		# summarize_glm_count ----------------------------------------------------------
2		#' Summarize Poisson negative binomial regression
3		#'
4		#' @description `r lifecycle::badge("experimental")`
5		#'
6		#' Summarize results of a Poisson negative binomial regression.
7		#' This can be used to analyze count and/or frequency data using a linear model.
8		#' It is specifically useful for analyzing count data (using the Poisson or Negative
9		#' Binomial distribution) that is result of a generalized linear model of one (e.g. arm) or more
10		#' covariates.
11		#'
12		#' @inheritParams h_glm_count
13		#' @inheritParams argument_convention
14		#' @param rate_mean_method (`character(1)`)\cr method used to estimate the mean odds ratio. Defaults to `emmeans`.
15		#' see details for more information.
16		#' @param scale (`numeric(1)`)\cr linear scaling factor for rate and confidence intervals. Defaults to `1`.
17		#' @param .stats (`character`)\cr statistics to select for the table.
18		#'
19		#' Options are: ``r shQuote(get_stats("summarize_glm_count"))``
20		#'
21		#' @details
22		#' `summarize_glm_count()` uses `s_glm_count()` to calculate the statistics for the table. This
23		#' analysis function uses [h_glm_count()] to estimate the GLM with [stats::glm()] for Poisson and Quasi-Poisson
24		#' distributions or [MASS::glm.nb()] for Negative Binomial distribution. All methods assume a
25		#' logarithmic link function.
26		#'
27		#' At this point, rates and confidence intervals are estimated from the model using
28		#' either [emmeans::emmeans()] when `rate_mean_method = "emmeans"` or [h_ppmeans()]
29		#' when `rate_mean_method = "ppmeans"`.
30		#'
31		#' If a reference group is specified while building the table with `split_cols_by(ref_group)`,
32		#' no rate ratio or `p-value` are calculated. Otherwise, we use [emmeans::contrast()] to
33		#' calculate the rate ratio and `p-value` for the reference group. Values are always estimated
34		#' with `method = "trt.vs.ctrl"` and `ref` equal to the first `arm` value.
35		#'
36		#' @name summarize_glm_count
37		NULL
38
39		#' @describeIn summarize_glm_count Layout-creating function which can take statistics function arguments
40		#' and additional format arguments. This function is a wrapper for [rtables::analyze()].
41		#'
42		#' @return
43		#' * `summarize_glm_count()` returns a layout object suitable for passing to further layouting functions,
44		#' or to [rtables::build_table()]. Adding this function to an `rtable` layout will add formatted rows containing
45		#' the statistics from `s_glm_count()` to the table layout.
46		#'
47		#' @examples
48		#' library(dplyr)
49		#'
50		#' anl <- tern_ex_adtte %>% filter(PARAMCD == "TNE")
51		#' anl$AVAL_f <- as.factor(anl$AVAL)
52		#'
53		#' lyt <- basic_table() %>%
54		#' split_cols_by("ARM", ref_group = "B: Placebo") %>%
55		#' add_colcounts() %>%
56		#' analyze_vars(
57		#' "AVAL_f",
58		#' var_labels = "Number of exacerbations per patient",
59		#' .stats = c("count_fraction"),
60		#' .formats = c("count_fraction" = "xx (xx.xx%)"),
61		#' .labels = c("Number of exacerbations per patient")
62		#' ) %>%
63		#' summarize_glm_count(
64		#' vars = "AVAL",
65		#' variables = list(arm = "ARM", offset = "lgTMATRSK", covariates = NULL),
66		#' conf_level = 0.95,
67		#' distribution = "poisson",
68		#' rate_mean_method = "emmeans",
69		#' var_labels = "Adjusted (P) exacerbation rate (per year)",
70		#' table_names = "adjP",
71		#' .stats = c("rate"),
72		#' .labels = c(rate = "Rate")
73		#' ) %>%
74		#' summarize_glm_count(
75		#' vars = "AVAL",
76		#' variables = list(arm = "ARM", offset = "lgTMATRSK", covariates = c("REGION1")),
77		#' conf_level = 0.95,
78		#' distribution = "quasipoisson",
79		#' rate_mean_method = "ppmeans",
80		#' var_labels = "Adjusted (QP) exacerbation rate (per year)",
81		#' table_names = "adjQP",
82		#' .stats = c("rate", "rate_ci", "rate_ratio", "rate_ratio_ci", "pval"),
83		#' .labels = c(
84		#' rate = "Rate", rate_ci = "Rate CI", rate_ratio = "Rate Ratio",
85		#' rate_ratio_ci = "Rate Ratio CI", pval = "p value"
86		#' )
87		#' )
88		#'
89		#' build_table(lyt = lyt, df = anl)
90		#'
91		#' @export
92		summarize_glm_count <- function(lyt,
93		vars,
94		variables,
95		distribution,
96		conf_level,
97		rate_mean_method = c("emmeans", "ppmeans")[1],
98		weights = stats::weights,
99		scale = 1,
100		var_labels,
101		na_str = default_na_str(),
102		nested = TRUE,
103		...,
104		show_labels = "visible",
105		table_names = vars,
106		.stats = get_stats("summarize_glm_count"),
107		.formats = NULL,
108		.labels = NULL,
109		.indent_mods = c(
110		"n" = 0L,
111		"rate" = 0L,
112		"rate_ci" = 1L,
113		"rate_ratio" = 0L,
114		"rate_ratio_ci" = 1L,
115		"pval" = 1L
116		)) {
117	3x	checkmate::assert_choice(rate_mean_method, c("emmeans", "ppmeans"))
118
119	3x	extra_args <- list(
120	3x	variables = variables, distribution = distribution, conf_level = conf_level,
121	3x	rate_mean_method = rate_mean_method, weights = weights, scale = scale, ...
122		)
123
124		# Selecting parameters following the statistics
125	3x	.formats <- get_formats_from_stats(.stats, formats_in = .formats)
126	3x	.labels <- get_labels_from_stats(.stats, labels_in = .labels)
127	3x	.indent_mods <- get_indents_from_stats(.stats, indents_in = .indent_mods)
128
129	3x	afun <- make_afun(
130	3x	s_glm_count,
131	3x	.stats = .stats,
132	3x	.formats = .formats,
133	3x	.labels = .labels,
134	3x	.indent_mods = .indent_mods,
135	3x	.null_ref_cells = FALSE
136		)
137
138	3x	analyze(
139	3x	lyt,
140	3x	vars,
141	3x	var_labels = var_labels,
142	3x	show_labels = show_labels,
143	3x	table_names = table_names,
144	3x	afun = afun,
145	3x	na_str = na_str,
146	3x	nested = nested,
147	3x	extra_args = extra_args
148		)
149		}
150
151		#' @describeIn summarize_glm_count Statistics function that produces a named list of results
152		#' of the investigated Poisson model.
153		#'
154		#' @return
155		#' * `s_glm_count()` returns a named `list` of 5 statistics:
156		#' * `n`: Count of complete sample size for the group.
157		#' * `rate`: Estimated event rate per follow-up time.
158		#' * `rate_ci`: Confidence level for estimated rate per follow-up time.
159		#' * `rate_ratio`: Ratio of event rates in each treatment arm to the reference arm.
160		#' * `rate_ratio_ci`: Confidence level for the rate ratio.
161		#' * `pval`: p-value.
162		#'
163		#' @keywords internal
164		s_glm_count <- function(df,
165		.var,
166		.df_row,
167		variables,
168		.ref_group,
169		.in_ref_col,
170		distribution,
171		conf_level,
172		rate_mean_method,
173		weights,
174		scale = 1) {
175	14x	arm <- variables$arm
176
177	14x	y <- df[[.var]]
178	13x	smry_level <- as.character(unique(df[[arm]]))
179
180		# ensure there is only 1 value
181	13x	checkmate::assert_scalar(smry_level)
182
183	13x	results <- h_glm_count(
184	13x	.var = .var,
185	13x	.df_row = .df_row,
186	13x	variables = variables,
187	13x	distribution = distribution,
188	13x	weights
189		)
190
191	13x	if (rate_mean_method == "emmeans") {
192	13x	emmeans_smry <- summary(results$emmeans_fit, level = conf_level)
193	!	} else if (rate_mean_method == "ppmeans") {
194	!	emmeans_smry <- h_ppmeans(results$glm_fit, .df_row, arm, conf_level)
195		}
196
197	13x	emmeans_smry_level <- emmeans_smry[emmeans_smry[[arm]] == smry_level, ]
198
199		# This happens if there is a reference col. No Ratio is calculated?
200	13x	if (.in_ref_col) {
201	5x	list(
202	5x	n = length(y[!is.na(y)]),
203	5x	rate = formatters::with_label(
204	5x	ifelse(distribution == "negbin", emmeans_smry_level$response * scale, emmeans_smry_level$rate * scale),
205	5x	"Adjusted Rate"
206		),
207	5x	rate_ci = formatters::with_label(
208	5x	c(emmeans_smry_level$asymp.LCL * scale, emmeans_smry_level$asymp.UCL * scale),
209	5x	f_conf_level(conf_level)
210		),
211	5x	rate_ratio = formatters::with_label(character(), "Adjusted Rate Ratio"),
212	5x	rate_ratio_ci = formatters::with_label(character(), f_conf_level(conf_level)),
213	5x	pval = formatters::with_label(character(), "p-value")
214		)
215		} else {
216	8x	emmeans_contrasts <- emmeans::contrast(
217	8x	results$emmeans_fit,
218	8x	method = "trt.vs.ctrl",
219	8x	ref = grep(
220	8x	as.character(unique(.ref_group[[arm]])),
221	8x	as.data.frame(results$emmeans_fit)[[arm]]
222		)
223		)
224
225	8x	contrasts_smry <- summary(
226	8x	emmeans_contrasts,
227	8x	infer = TRUE,
228	8x	adjust = "none"
229		)
230
231	8x	smry_contrasts_level <- contrasts_smry[grepl(smry_level, contrasts_smry$contrast), ]
232
233	8x	list(
234	8x	n = length(y[!is.na(y)]),
235	8x	rate = formatters::with_label(
236	8x	ifelse(distribution == "negbin",
237	8x	emmeans_smry_level$response * scale,
238	8x	emmeans_smry_level$rate * scale
239		),
240	8x	"Adjusted Rate"
241		),
242	8x	rate_ci = formatters::with_label(
243	8x	c(emmeans_smry_level$asymp.LCL * scale, emmeans_smry_level$asymp.UCL * scale),
244	8x	f_conf_level(conf_level)
245		),
246	8x	rate_ratio = formatters::with_label(
247	8x	smry_contrasts_level$ratio,
248	8x	"Adjusted Rate Ratio"
249		),
250	8x	rate_ratio_ci = formatters::with_label(
251	8x	c(smry_contrasts_level$asymp.LCL, smry_contrasts_level$asymp.UCL),
252	8x	f_conf_level(conf_level)
253		),
254	8x	pval = formatters::with_label(
255	8x	smry_contrasts_level$p.value,
256	8x	"p-value"
257		)
258		)
259		}
260		}
261		# h_glm_count ------------------------------------------------------------------
262		#' Helper functions for Poisson models
263		#'
264		#' @description `r lifecycle::badge("experimental")`
265		#'
266		#' Helper functions that returns the results of [stats::glm()] when Poisson or Quasi-Poisson
267		#' distributions are needed (see `family` parameter), or [MASS::glm.nb()] for Negative Binomial
268		#' distributions. Link function for the GLM is `log`.
269		#'
270		#' @inheritParams argument_convention
271		#'
272		#' @seealso [summarize_glm_count]
273		#'
274		#' @name h_glm_count
275		NULL
276
277		#' @describeIn h_glm_count Helper function to return the results of the
278		#' selected model (Poisson, Quasi-Poisson, negative binomial).
279		#'
280		#' @param .df_row (`data.frame`)\cr dataset that includes all the variables that are called
281		#' in `.var` and `variables`.
282		#' @param variables (named `list` of `string`)\cr list of additional analysis variables, with
283		#' expected elements:
284		#' * `arm` (`string`)\cr group variable, for which the covariate adjusted means of multiple
285		#' groups will be summarized. Specifically, the first level of `arm` variable is taken as the
286		#' reference group.
287		#' * `covariates` (`character`)\cr a vector that can contain single variable names (such as
288		#' `"X1"`), and/or interaction terms indicated by `"X1 * X2"`.
289		#' * `offset` (`numeric`)\cr a numeric vector or scalar adding an offset.
290		#' @param distribution (`character`)\cr a character value specifying the distribution
291		#' used in the regression (Poisson, Quasi-Poisson, negative binomial).
292		#' @param weights (`character`)\cr a character vector specifying weights used
293		#' in averaging predictions. Number of weights must equal the number of levels included in the covariates.
294		#' Weights option passed to [emmeans::emmeans()].
295		#'
296		#' @return
297		#' * `h_glm_count()` returns the results of the selected model.
298		#'
299		#' @keywords internal
300		h_glm_count <- function(.var,
301		.df_row,
302		variables,
303		distribution,
304		weights) {
305	21x	checkmate::assert_subset(distribution, c("poisson", "quasipoisson", "negbin"), empty.ok = FALSE)
306	19x	switch(distribution,
307	13x	poisson = h_glm_poisson(.var, .df_row, variables, weights),
308	1x	quasipoisson = h_glm_quasipoisson(.var, .df_row, variables, weights),
309	5x	negbin = h_glm_negbin(.var, .df_row, variables, weights)
310		)
311		}
312
313		#' @describeIn h_glm_count Helper function to return results of a Poisson model.
314		#'
315		#' @return
316		#' * `h_glm_poisson()` returns the results of a Poisson model.
317		#'
318		#' @keywords internal
319		h_glm_poisson <- function(.var,
320		.df_row,
321		variables,
322		weights) {
323	17x	arm <- variables$arm
324	17x	covariates <- variables$covariates
325	17x	offset <- .df_row[[variables$offset]]
326
327	15x	formula <- stats::as.formula(paste0(
328	15x	.var, " ~ ",
329		" + ",
330	15x	paste(covariates, collapse = " + "),
331		" + ",
332	15x	arm
333		))
334
335	15x	glm_fit <- stats::glm(
336	15x	formula = formula,
337	15x	offset = offset,
338	15x	data = .df_row,
339	15x	family = stats::poisson(link = "log")
340		)
341
342	15x	emmeans_fit <- emmeans::emmeans(
343	15x	glm_fit,
344	15x	specs = arm,
345	15x	data = .df_row,
346	15x	type = "response",
347	15x	offset = 0,
348	15x	weights = weights
349		)
350
351	15x	list(
352	15x	glm_fit = glm_fit,
353	15x	emmeans_fit = emmeans_fit
354		)
355		}
356
357		#' @describeIn h_glm_count Helper function to return results of a Quasi-Poisson model.
358		#'
359		#' @return
360		#' * `h_glm_quasipoisson()` returns the results of a Quasi-Poisson model.
361		#'
362		#' @keywords internal
363		h_glm_quasipoisson <- function(.var,
364		.df_row,
365		variables,
366		weights) {
367	5x	arm <- variables$arm
368	5x	covariates <- variables$covariates
369	5x	offset <- .df_row[[variables$offset]]
370
371	3x	formula <- stats::as.formula(paste0(
372	3x	.var, " ~ ",
373		" + ",
374	3x	paste(covariates, collapse = " + "),
375		" + ",
376	3x	arm
377		))
378
379	3x	glm_fit <- stats::glm(
380	3x	formula = formula,
381	3x	offset = offset,
382	3x	data = .df_row,
383	3x	family = stats::quasipoisson(link = "log")
384		)
385
386	3x	emmeans_fit <- emmeans::emmeans(
387	3x	glm_fit,
388	3x	specs = arm,
389	3x	data = .df_row,
390	3x	type = "response",
391	3x	offset = 0,
392	3x	weights = weights
393		)
394
395	3x	list(
396	3x	glm_fit = glm_fit,
397	3x	emmeans_fit = emmeans_fit
398		)
399		}
400
401		#' @describeIn h_glm_count Helper function to return results of a negative binomial model.
402		#'
403		#' @return
404		#' * `h_glm_negbin()` returns the results of a negative binomial model.
405		#'
406		#' @keywords internal
407		h_glm_negbin <- function(.var,
408		.df_row,
409		variables,
410		weights) {
411	9x	arm <- variables$arm
412	9x	covariates <- variables$covariates
413
414	9x	formula <- stats::as.formula(paste0(
415	9x	.var, " ~ ",
416		" + ",
417	9x	paste(covariates, collapse = " + "),
418		" + ",
419	9x	arm
420		))
421
422	9x	glm_fit <- MASS::glm.nb(
423	9x	formula = formula,
424	9x	data = .df_row,
425	9x	link = "log"
426		)
427
428	7x	emmeans_fit <- emmeans::emmeans(
429	7x	glm_fit,
430	7x	specs = arm,
431	7x	data = .df_row,
432	7x	type = "response",
433	7x	offset = 0,
434	7x	weights = weights
435		)
436
437	7x	list(
438	7x	glm_fit = glm_fit,
439	7x	emmeans_fit = emmeans_fit
440		)
441		}
442
443		# h_ppmeans --------------------------------------------------------------------
444		#' Function to return the estimated means using predicted probabilities
445		#'
446		#' @description
447		#' For each arm level, the predicted mean rate is calculated using the fitted model object, with `newdata`
448		#' set to the result of `stats::model.frame`, a reconstructed data or the original data, depending on the
449		#' object formula (coming from the fit). The confidence interval is derived using the `conf_level` parameter.
450		#'
451		#' @param obj (`glm.fit`)\cr fitted model object used to derive the mean rate estimates in each treatment arm.
452		#' @param .df_row (`data.frame`)\cr dataset that includes all the variables that are called in `.var` and `variables`.
453		#' @param arm (`string`)\cr group variable, for which the covariate adjusted means of multiple groups will be
454		#' summarized. Specifically, the first level of `arm` variable is taken as the reference group.
455		#' @param conf_level (`proportion`)\cr value used to derive the confidence interval for the rate.
456		#'
457		#' @return
458		#' * `h_ppmeans()` returns the estimated means.
459		#'
460		#' @seealso [summarize_glm_count()].
461		#'
462		#' @export
463		h_ppmeans <- function(obj, .df_row, arm, conf_level) {
464	1x	alpha <- 1 - conf_level
465	1x	p <- 1 - alpha / 2
466
467	1x	arm_levels <- levels(.df_row[[arm]])
468
469	1x	out <- lapply(arm_levels, function(lev) {
470	3x	temp <- .df_row
471	3x	temp[[arm]] <- factor(lev, levels = arm_levels)
472
473	3x	mf <- stats::model.frame(obj$formula, data = temp)
474	3x	X <- stats::model.matrix(obj$formula, data = mf) # nolint
475
476	3x	rate <- stats::predict(obj, newdata = mf, type = "response")
477	3x	rate_hat <- mean(rate)
478
479	3x	zz <- colMeans(rate * X)
480	3x	se <- sqrt(as.numeric(t(zz) %% stats::vcov(obj) %% zz))
481	3x	rate_lwr <- rate_hat * exp(-stats::qnorm(p) * se / rate_hat)
482	3x	rate_upr <- rate_hat * exp(stats::qnorm(p) * se / rate_hat)
483
484	3x	c(rate_hat, rate_lwr, rate_upr)
485		})
486
487	1x	names(out) <- arm_levels
488	1x	out <- do.call(rbind, out)
489	1x	if ("negbin" %in% class(obj)) {
490	!	colnames(out) <- c("response", "asymp.LCL", "asymp.UCL")
491		} else {
492	1x	colnames(out) <- c("rate", "asymp.LCL", "asymp.UCL")
493		}
494	1x	out <- as.data.frame(out)
495	1x	out[[arm]] <- rownames(out)
496	1x	out
497		}

1		#' Compare variables between groups
2		#'
3		#' @description `r lifecycle::badge("stable")`
4		#'
5		#' The analyze function [compare_vars()] creates a layout element to summarize and compare one or more variables, using
6		#' the S3 generic function [s_summary()] to calculate a list of summary statistics. A list of all available statistics
7		#' for numeric variables can be viewed by running `get_stats("analyze_vars_numeric", add_pval = TRUE)` and for
8		#' non-numeric variables by running `get_stats("analyze_vars_counts", add_pval = TRUE)`. Use the `.stats` parameter to
9		#' specify the statistics to include in your output summary table.
10		#'
11		#' Prior to using this function in your table layout you must use [rtables::split_cols_by()] to create a column
12		#' split on the variable to be used in comparisons, and specify a reference group via the `ref_group` parameter.
13		#' Comparisons can be performed for each group (column) against the specified reference group by including the p-value
14		#' statistic.
15		#'
16		#' @inheritParams argument_convention
17		#' @param .stats (`character`)\cr statistics to select for the table.
18		#'
19		#' Options for numeric variables are: ``r shQuote(get_stats("analyze_vars_numeric", add_pval = TRUE))``
20		#'
21		#' Options for non-numeric variables are: ``r shQuote(get_stats("analyze_vars_counts", add_pval = TRUE))``
22		#'
23