tern coverage - 94.72%

Files
Source

#' Line plot with optional table
#'
#' @description `r lifecycle::badge("stable")`
#'
#' Line plot with the optional table.
#'
#' @inheritParams argument_convention
#' @param alt_counts_df (`data.frame` or `NULL`)\cr data set that will be used (only)
#'   to counts objects in groups for stratification.
#' @param variables (named `character`) vector of variable names in `df` which should include:
#'   * `x` (`string`)\cr name of x-axis variable.
#'   * `y` (`string`)\cr name of y-axis variable.
#'   * `group_var` (`string` or `NULL`)\cr name of grouping variable (or strata), i.e. treatment arm.
#'     Can be `NA` to indicate lack of groups.
#'   * `subject_var` (`string` or `NULL`)\cr name of subject variable. Only applies if `group_var` is
#'      not NULL.
#'   * `paramcd` (`string` or `NA`)\cr name of the variable for parameter's code. Used for y-axis label and plot's
#'     subtitle. Can be `NA` if `paramcd` is not to be added to the y-axis label or subtitle.
#'   * `y_unit` (`string` or `NA`)\cr name of variable with units of `y`. Used for y-axis label and plot's subtitle.
#'     Can be `NA` if y unit is not to be added to the y-axis label or subtitle.
#'   * `facet_var` (`string` or `NA`)\cr name of the secondary grouping variable used for plot faceting, i.e. treatment
#'     arm. Can be `NA` to indicate lack of groups.
#' @param mid (`character` or `NULL`)\cr names of the statistics that will be plotted as midpoints.
#'   All the statistics indicated in `mid` variable must be present in the object returned by `sfun`,
#'   and be of a `double` or `numeric` type vector of length one.
#' @param interval (`character` or `NULL`)\cr names of the statistics that will be plotted as intervals.
#'   All the statistics indicated in `interval` variable must be present in the object returned by `sfun`,
#'   and be of a `double` or `numeric` type vector of length two. Set `interval = NULL` if intervals should not be
#'   added to the plot.
#' @param whiskers (`character`)\cr names of the interval whiskers that will be plotted. Names must match names
#'   of the list element `interval` that will be returned by `sfun` (e.g. `mean_ci_lwr` element of
#'   `sfun(x)[["mean_ci"]]`). It is possible to specify one whisker only, or to suppress all whiskers by setting
#'   `interval = NULL`.
#' @param table (`character` or `NULL`)\cr names of the statistics that will be displayed in the table below the plot.
#'   All the statistics indicated in `table` variable must be present in the object returned by `sfun`.
#' @param sfun (`function`)\cr the function to compute the values of required statistics. It must return a named `list`
#'   with atomic vectors. The names of the `list` elements refer to the names of the statistics and are used by `mid`,
#'   `interval`, `table`. It must be able to accept as input a vector with data for which statistics are computed.
#' @param ... optional arguments to `sfun`.
#' @param mid_type (`string`)\cr controls the type of the `mid` plot, it can be point (`"p"`), line (`"l"`),
#'   or point and line (`"pl"`).
#' @param mid_point_size (`numeric(1)`)\cr font size of the `mid` plot points.
#' @param position (`character` or `call`)\cr geom element position adjustment, either as a string, or the result of
#'   a call to a position adjustment function.
#' @param legend_title (`string`)\cr legend title.
#' @param legend_position (`string`)\cr the position of the plot legend (`"none"`, `"left"`, `"right"`, `"bottom"`,
#'   `"top"`, or a two-element numeric vector).
#' @param ggtheme (`theme`)\cr a graphical theme as provided by `ggplot2` to control styling of the plot.
#' @param xticks (`numeric` or `NULL`)\cr numeric vector of tick positions or a single number with spacing
#'   between ticks on the x-axis, for use when `variables$x` is numeric. If `NULL` (default), [labeling::extended()] is
#'   used to determine optimal tick positions on the x-axis. If `variables$x` is not numeric, this argument is ignored.
#' @param x_lab (`string` or `NULL`)\cr x-axis label. If `NULL` then no label will be added.
#' @param y_lab (`string` or `NULL`)\cr y-axis label. If `NULL` then no label will be added.
#' @param y_lab_add_paramcd (`flag`)\cr whether `paramcd`, i.e. `unique(df[[variables["paramcd"]]])` should be added
#'   to the y-axis label (`y_lab`).
#' @param y_lab_add_unit (`flag`)\cr whether y-axis unit, i.e. `unique(df[[variables["y_unit"]]])` should be added
#'   to the y-axis label (`y_lab`).
#' @param title (`string`)\cr plot title.
#' @param subtitle (`string`)\cr plot subtitle.
#' @param subtitle_add_paramcd (`flag`)\cr whether `paramcd`, i.e. `unique(df[[variables["paramcd"]]])` should be
#'   added to the plot's subtitle (`subtitle`).
#' @param subtitle_add_unit (`flag`)\cr whether the y-axis unit, i.e. `unique(df[[variables["y_unit"]]])` should be
#'   added to the plot's subtitle (`subtitle`).
#' @param caption (`string`)\cr optional caption below the plot.
#' @param table_format (named `character` or `NULL`)\cr format patterns for descriptive statistics used in the
#'   (optional) table appended to the plot. It is passed directly to the `h_format_row` function through the `format`
#'   parameter. Names of `table_format` must match the names of statistics returned by `sfun` function.
#' @param table_labels (named `character` or `NULL`)\cr labels for descriptive statistics used in the (optional) table
#'   appended to the plot. Names of `table_labels` must match the names of statistics returned by `sfun` function.
#' @param table_font_size (`numeric(1)`)\cr font size of the text in the table.
#' @param newpage `r lifecycle::badge("deprecated")` not used.
#' @param col (`character`)\cr color(s). See `?ggplot2::aes_colour_fill_alpha` for example values.
#' @param linetype (`character`)\cr line type(s). See `?ggplot2::aes_linetype_size_shape` for example values.
#' @param errorbar_width (`numeric(1)`)\cr width of the error bars.
#'
#' @return A `ggplot` line plot (and statistics table if applicable).
#'
#' @examples
#' library(nestcolor)
#'
#' adsl <- tern_ex_adsl
#' adlb <- tern_ex_adlb %>% dplyr::filter(ANL01FL == "Y", PARAMCD == "ALT", AVISIT != "SCREENING")
#' adlb$AVISIT <- droplevels(adlb$AVISIT)
#' adlb <- dplyr::mutate(adlb, AVISIT = forcats::fct_reorder(AVISIT, AVISITN, min))
#'
#' # Mean with CI
#' g_lineplot(adlb, adsl, subtitle = "Laboratory Test:")
#'
#' # Mean with CI, no stratification with group_var
#' g_lineplot(adlb, variables = control_lineplot_vars(group_var = NA))
#'
#' # Mean, upper whisker of CI, no group_var(strata) counts N
#' g_lineplot(
#'   adlb,
#'   whiskers = "mean_ci_upr",
#'   title = "Plot of Mean and Upper 95% Confidence Limit by Visit"
#' )
#'
#' # Median with CI
#' g_lineplot(
#'   adlb,
#'   adsl,
#'   mid = "median",
#'   interval = "median_ci",
#'   whiskers = c("median_ci_lwr", "median_ci_upr"),
#'   title = "Plot of Median and 95% Confidence Limits by Visit"
#' )
#'
#' # Mean, +/- SD
#' g_lineplot(adlb, adsl,
#'   interval = "mean_sdi",
#'   whiskers = c("mean_sdi_lwr", "mean_sdi_upr"),
#'   title = "Plot of Median +/- SD by Visit"
#' )
#'
#' # Mean with CI plot with stats table
#' g_lineplot(adlb, adsl, table = c("n", "mean", "mean_ci"))
#'
#' # Mean with CI, table and customized confidence level
#' g_lineplot(
#'   adlb,
#'   adsl,
#'   table = c("n", "mean", "mean_ci"),
#'   control = control_analyze_vars(conf_level = 0.80),
#'   title = "Plot of Mean and 80% Confidence Limits by Visit"
#' )
#'
#' # Mean with CI, table, filtered data
#' adlb_f <- dplyr::filter(adlb, ARMCD != "ARM A" | AVISIT == "BASELINE")
#' g_lineplot(adlb_f, table = c("n", "mean"))
#'
#' @export
g_lineplot <- function(df,
                       alt_counts_df = NULL,
                       variables = control_lineplot_vars(),
                       mid = "mean",
                       interval = "mean_ci",
                       whiskers = c("mean_ci_lwr", "mean_ci_upr"),
                       table = NULL,
                       sfun = s_summary,
                       ...,
                       mid_type = "pl",
                       mid_point_size = 2,
                       position = ggplot2::position_dodge(width = 0.4),
                       legend_title = NULL,
                       legend_position = "bottom",
                       ggtheme = nestcolor::theme_nest(),
                       xticks = NULL,
                       xlim = NULL,
                       ylim = NULL,
                       x_lab = obj_label(df[[variables[["x"]]]]),
                       y_lab = NULL,
                       y_lab_add_paramcd = TRUE,
                       y_lab_add_unit = TRUE,
                       title = "Plot of Mean and 95% Confidence Limits by Visit",
                       subtitle = "",
                       subtitle_add_paramcd = TRUE,
                       subtitle_add_unit = TRUE,
                       caption = NULL,
                       table_format = NULL,
                       table_labels = NULL,
                       table_font_size = 3,
                       errorbar_width = 0.45,
                       newpage = lifecycle::deprecated(),
                       col = NULL,
                       linetype = NULL) {
  checkmate::assert_character(variables, any.missing = TRUE)
  checkmate::assert_character(mid, null.ok = TRUE)
  checkmate::assert_character(interval, null.ok = TRUE)
  checkmate::assert_character(col, null.ok = TRUE)
  checkmate::assert_character(linetype, null.ok = TRUE)
  checkmate::assert_numeric(xticks, null.ok = TRUE)
  checkmate::assert_numeric(xlim, finite = TRUE, any.missing = FALSE, len = 2, sorted = TRUE, null.ok = TRUE)
  checkmate::assert_numeric(ylim, finite = TRUE, any.missing = FALSE, len = 2, sorted = TRUE, null.ok = TRUE)
  checkmate::assert_number(errorbar_width, lower = 0)
  checkmate::assert_string(title, null.ok = TRUE)
  checkmate::assert_string(subtitle, null.ok = TRUE)

  if (!is.null(table)) {
    table_format <- get_formats_from_stats(table)
    table_labels <- get_labels_from_stats(table)
  }

  extra_args <- list(...)
  if ("control" %in% names(extra_args)) {
    if (!is.null(table) && all(table_labels == get_labels_from_stats(table))) {
      table_labels <- table_labels %>% labels_use_control(extra_args[["control"]])
    }
  }

  if (is.character(interval)) {
    checkmate::assert_vector(whiskers, min.len = 0, max.len = 2)
  }

  if (length(whiskers) == 1) {
    checkmate::assert_character(mid)
  }

  if (is.character(mid)) {
    checkmate::assert_scalar(mid_type)
    checkmate::assert_subset(mid_type, c("pl", "p", "l"))
  }

  x <- variables[["x"]]
  y <- variables[["y"]]
  paramcd <- variables["paramcd"] # NA if paramcd == NA or it is not in variables
  y_unit <- variables["y_unit"] # NA if y_unit == NA or it is not in variables
  if (is.na(variables["group_var"])) {
    group_var <- NULL # NULL if group_var == NA or it is not in variables
  } else {
    group_var <- variables[["group_var"]]
    subject_var <- variables[["subject_var"]]
  }
  if (is.na(variables["facet_var"])) {
    facet_var <- NULL # NULL if facet_var == NA or it is not in variables
  } else {
    facet_var <- variables[["facet_var"]]
  }
  checkmate::assert_flag(y_lab_add_paramcd, null.ok = TRUE)
  checkmate::assert_flag(subtitle_add_paramcd, null.ok = TRUE)
  if ((!is.null(y_lab) && y_lab_add_paramcd) || (!is.null(subtitle) && subtitle_add_paramcd)) {
    checkmate::assert_false(is.na(paramcd))
    checkmate::assert_scalar(unique(df[[paramcd]]))
  }

  checkmate::assert_flag(y_lab_add_unit, null.ok = TRUE)
  checkmate::assert_flag(subtitle_add_unit, null.ok = TRUE)
  if ((!is.null(y_lab) && y_lab_add_unit) || (!is.null(subtitle) && subtitle_add_unit)) {
    checkmate::assert_false(is.na(y_unit))
    checkmate::assert_scalar(unique(df[[y_unit]]))
  }

  if (!is.null(group_var) && !is.null(alt_counts_df)) {
    checkmate::assert_set_equal(unique(alt_counts_df[[group_var]]), unique(df[[group_var]]))
  }

  ####################################### |
  # ---- Compute required statistics ----
  ####################################### |
  # Remove unused levels for x-axis
  if (is.factor(df[[x]])) {
    df[[x]] <- droplevels(df[[x]])
  }

  if (!is.null(facet_var) && !is.null(group_var)) {
    df_grp <- tidyr::expand(df, .data[[facet_var]], .data[[group_var]], .data[[x]]) # expand based on levels of factors
  } else if (!is.null(group_var)) {
    df_grp <- tidyr::expand(df, .data[[group_var]], .data[[x]]) # expand based on levels of factors
  } else {
    df_grp <- tidyr::expand(df, NULL, .data[[x]])
  }

  df_grp <- df_grp %>%
    dplyr::full_join(y = df[, c(facet_var, group_var, x, y)], by = c(facet_var, group_var, x), multiple = "all") %>%
    dplyr::group_by_at(c(facet_var, group_var, x))

  df_stats <- df_grp %>%
    dplyr::summarise(
      data.frame(t(do.call(c, unname(sfun(.data[[y]])[c(mid, interval)])))),
      .groups = "drop"
    )

  df_stats <- df_stats[!is.na(df_stats[[mid]]), ]

  # add number of objects N in group_var (strata)
  if (!is.null(group_var) && !is.null(alt_counts_df)) {
    strata_N <- paste0(group_var, "_N") # nolint

    df_N <- stats::aggregate(eval(parse(text = subject_var)) ~ eval(parse(text = group_var)), data = alt_counts_df, FUN = function(x) length(unique(x))) # nolint
    colnames(df_N) <- c(group_var, "N") # nolint
    df_N[[strata_N]] <- paste0(df_N[[group_var]], " (N = ", df_N$N, ")") # nolint

    # keep strata factor levels
    matches <- sapply(unique(df_N[[group_var]]), function(x) {
      regex_pattern <- gsub("([][(){}^$.|*+?\\\\])", "\\\\\\1", x)
      unique(df_N[[paste0(group_var, "_N")]])[grepl(
        paste0("^", regex_pattern),
        unique(df_N[[paste0(group_var, "_N")]])
      )]
    })
    df_N[[paste0(group_var, "_N")]] <- factor(df_N[[group_var]]) # nolint
    levels(df_N[[paste0(group_var, "_N")]]) <- unlist(matches) # nolint

    # strata_N should not be in colnames(df_stats)
    checkmate::assert_disjunct(strata_N, colnames(df_stats))

    df_stats <- merge(x = df_stats, y = df_N[, c(group_var, strata_N)], by = group_var)
  } else if (!is.null(group_var)) {
    strata_N <- group_var # nolint
  } else {
    strata_N <- NULL # nolint
  }

  ############################################### |
  # ---- Prepare certain plot's properties. ----
  ############################################### |
  # legend title
  if (is.null(legend_title) && !is.null(group_var) && legend_position != "none") {
    legend_title <- attr(df[[group_var]], "label")
  }

  # y label
  if (!is.null(y_lab)) {
    if (y_lab_add_paramcd) {
      y_lab <- paste(y_lab, unique(df[[paramcd]]))
    }

    if (y_lab_add_unit) {
      y_lab <- paste0(y_lab, " (", unique(df[[y_unit]]), ")")
    }

    y_lab <- trimws(y_lab)
  }

  # subtitle
  if (!is.null(subtitle)) {
    if (subtitle_add_paramcd) {
      subtitle <- paste(subtitle, unique(df[[paramcd]]))
    }

    if (subtitle_add_unit) {
      subtitle <- paste0(subtitle, " (", unique(df[[y_unit]]), ")")
    }

    subtitle <- trimws(subtitle)
  }

  ############################### |
  # ---- Build plot object. ----
  ############################### |
  p <- ggplot2::ggplot(
    data = df_stats,
    mapping = ggplot2::aes(
      x = .data[[x]], y = .data[[mid]],
      color = if (is.null(strata_N)) NULL else .data[[strata_N]],
      shape = if (is.null(strata_N)) NULL else .data[[strata_N]],
      lty = if (is.null(strata_N)) NULL else .data[[strata_N]],
      group = if (is.null(strata_N)) NULL else .data[[strata_N]]
    )
  )

  if (!is.null(group_var) && nlevels(df_stats[[strata_N]]) > 6) {
    p <- p +
      scale_shape_manual(values = seq(15, 15 + nlevels(df_stats[[strata_N]])))
  }

  if (!is.null(mid)) {
    # points
    if (grepl("p", mid_type, fixed = TRUE)) {
      p <- p + ggplot2::geom_point(position = position, size = mid_point_size, na.rm = TRUE)
    }

    # lines - plotted only if there is a strata grouping (group_var)
    if (grepl("l", mid_type, fixed = TRUE) && !is.null(strata_N)) { # nolint
      p <- p + ggplot2::geom_line(position = position, na.rm = TRUE)
    }
  }

  # interval
  if (!is.null(interval)) {
    p <- p +
      ggplot2::geom_errorbar(
        ggplot2::aes(ymin = .data[[whiskers[1]]], ymax = .data[[whiskers[max(1, length(whiskers))]]]),
        width = errorbar_width,
        position = position
      )

    if (length(whiskers) == 1) { # lwr or upr only; mid is then required
      # workaround as geom_errorbar does not provide single-direction whiskers
      p <- p +
        ggplot2::geom_linerange(
          data = df_stats[!is.na(df_stats[[whiskers]]), ], # as na.rm =TRUE does not suppress warnings
          ggplot2::aes(ymin = .data[[mid]], ymax = .data[[whiskers]]),
          position = position,
          na.rm = TRUE,
          show.legend = FALSE
        )
    }
  }

  if (is.numeric(df_stats[[x]])) {
    if (length(xticks) == 1) xticks <- seq(from = min(df_stats[[x]]), to = max(df_stats[[x]]), by = xticks)
    p <- p + ggplot2::scale_x_continuous(breaks = if (!is.null(xticks)) xticks else waiver(), limits = xlim)
  }

  p <- p +
    ggplot2::scale_y_continuous(labels = scales::comma, limits = ylim) +
    ggplot2::labs(
      title = title,
      subtitle = subtitle,
      caption = caption,
      color = legend_title,
      lty = legend_title,
      shape = legend_title,
      x = x_lab,
      y = y_lab
    )

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

  if (!is.null(facet_var)) {
    p <- p +
      facet_grid(cols = vars(df_stats[[facet_var]]))
  }

  if (!is.null(ggtheme)) {
    p <- p + ggtheme
  } else {
    p <- p +
      ggplot2::theme_bw() +
      ggplot2::theme(
        legend.key.width = grid::unit(1, "cm"),
        legend.position = legend_position,
        legend.direction = ifelse(
          legend_position %in% c("top", "bottom"),
          "horizontal",
          "vertical"
        )
      )
  }

  ############################################################# |
  # ---- Optionally, add table to the bottom of the plot. ----
  ############################################################# |
  if (!is.null(table)) {
    df_stats_table <- df_grp %>%
      dplyr::summarise(
        h_format_row(
          x = sfun(.data[[y]], ...)[table],
          format = table_format,
          labels = table_labels
        ),
        .groups = "drop"
      )

    stats_lev <- rev(setdiff(colnames(df_stats_table), c(group_var, x)))

    df_stats_table <- df_stats_table %>%
      tidyr::pivot_longer(
        cols = -dplyr::all_of(c(group_var, x)),
        names_to = "stat",
        values_to = "value",
        names_ptypes = list(stat = factor(levels = stats_lev))
      )

    tbl <- ggplot2::ggplot(
      df_stats_table,
      ggplot2::aes(x = .data[[x]], y = .data[["stat"]], label = .data[["value"]])
    ) +
      ggplot2::geom_text(size = table_font_size) +
      ggplot2::theme_bw() +
      ggplot2::theme(
        panel.border = ggplot2::element_blank(),
        panel.grid.major = ggplot2::element_blank(),
        panel.grid.minor = ggplot2::element_blank(),
        axis.ticks = ggplot2::element_blank(),
        axis.title = ggplot2::element_blank(),
        axis.text.x = ggplot2::element_blank(),
        axis.text.y = ggplot2::element_text(margin = ggplot2::margin(t = 0, r = 0, b = 0, l = 5)),
        strip.text = ggplot2::element_text(hjust = 0),
        strip.text.x = ggplot2::element_text(margin = ggplot2::margin(1.5, 0, 1.5, 0, "pt")),
        strip.background = ggplot2::element_rect(fill = "grey95", color = NA),
        legend.position = "none"
      )

    if (!is.null(group_var)) {
      tbl <- tbl + ggplot2::facet_wrap(facets = group_var, ncol = 1)
    }

    # align plot and table
    cowplot::plot_grid(p, tbl, ncol = 1, align = "v", axis = "tblr")
  } else {
    p
  }
}

#' Helper function to format the optional `g_lineplot` table
#'
#' @description `r lifecycle::badge("stable")`
#'
#' @param x (named `list`)\cr list of numerical values to be formatted and optionally labeled.
#'   Elements of `x` must be `numeric` vectors.
#' @param format (named `character` or `NULL`)\cr format patterns for `x`. Names of the `format` must
#'   match the names of `x`. This parameter is passed directly to the `rtables::format_rcell`
#'   function through the `format` parameter.
#' @param labels (named `character` or `NULL`)\cr optional labels for `x`. Names of the `labels` must
#'   match the names of `x`. When a label is not specified for an element of `x`,
#'   then this function tries to use `label` or `names` (in this order) attribute of that element
#'   (depending on which one exists and it is not `NULL` or `NA` or `NaN`). If none of these attributes
#'   are attached to a given element of `x`, then the label is automatically generated.
#'
#' @return A single row `data.frame` object.
#'
#' @examples
#' mean_ci <- c(48, 51)
#' x <- list(mean = 50, mean_ci = mean_ci)
#' format <- c(mean = "xx.x", mean_ci = "(xx.xx, xx.xx)")
#' labels <- c(mean = "My Mean")
#' h_format_row(x, format, labels)
#'
#' attr(mean_ci, "label") <- "Mean 95% CI"
#' x <- list(mean = 50, mean_ci = mean_ci)
#' h_format_row(x, format, labels)
#'
#' @export
h_format_row <- function(x, format, labels = NULL) {
  # cell: one row, one column data.frame
  format_cell <- function(x, format, label = NULL) {
    fc <- format_rcell(x = x, format = unlist(format))
    if (is.na(fc)) {
      fc <- "NA"
    }
    x_label <- attr(x, "label")
    if (!is.null(label) && !is.na(label)) {
      names(fc) <- label
    } else if (!is.null(x_label) && !is.na(x_label)) {
      names(fc) <- x_label
    } else if (length(x) == length(fc)) {
      names(fc) <- names(x)
    }
    as.data.frame(t(fc))
  }

  row <- do.call(
    cbind,
    lapply(
      names(x), function(xn) format_cell(x[[xn]], format = format[xn], label = labels[xn])
    )
  )

  row
}

#' Control function for `g_lineplot()`
#'
#' @description `r lifecycle::badge("stable")`
#'
#' Default values for `variables` parameter in `g_lineplot` function.
#' A variable's default value can be overwritten for any variable.
#'
#' @param x (`string`)\cr x-variable name.
#' @param y (`string`)\cr y-variable name.
#' @param group_var (`string` or `NA`)\cr group variable name.
#' @param subject_var (`string` or `NA`)\cr subject variable name.
#' @param facet_var (`string` or `NA`)\cr faceting variable name.
#' @param paramcd (`string` or `NA`)\cr parameter code variable name.
#' @param y_unit (`string` or `NA`)\cr y-axis unit variable name.
#'
#' @return A named character vector of variable names.
#'
#' @examples
#' control_lineplot_vars()
#' control_lineplot_vars(group_var = NA)
#'
#' @export
control_lineplot_vars <- function(x = "AVISIT",
                                  y = "AVAL",
                                  group_var = "ARM",
                                  facet_var = NA,
                                  paramcd = "PARAMCD",
                                  y_unit = "AVALU",
                                  subject_var = "USUBJID") {
  checkmate::assert_string(x)
  checkmate::assert_string(y)
  checkmate::assert_string(group_var, na.ok = TRUE, null.ok = TRUE)
  checkmate::assert_string(facet_var, na.ok = TRUE, null.ok = TRUE)
  checkmate::assert_string(subject_var, na.ok = TRUE, null.ok = TRUE)
  checkmate::assert_string(paramcd, na.ok = TRUE, null.ok = TRUE)
  checkmate::assert_string(y_unit, na.ok = TRUE, null.ok = TRUE)

  variables <- c(
    x = x, y = y, group_var = group_var, paramcd = paramcd,
    y_unit = y_unit, subject_var = subject_var, facet_var = facet_var
  )
  return(variables)
}

#' Count patients with toxicity grades that have worsened from baseline by highest grade post-baseline
#'
#' @description `r lifecycle::badge("stable")`
#'
#' The analyze function [count_abnormal_lab_worsen_by_baseline()] creates a layout element to count patients with
#' analysis toxicity grades which have worsened from baseline, categorized by highest (worst) grade post-baseline.
#'
#' This function analyzes primary analysis variable `var` which indicates analysis toxicity grades. Additional
#' analysis variables that can be supplied as a list via the `variables` parameter are `id` (defaults to `USUBJID`),
#' a variable to indicate unique subject identifiers, `baseline_var` (defaults to `BTOXGR`), a variable to indicate
#' baseline toxicity grades, and `direction_var` (defaults to `GRADDIR`), a variable to indicate toxicity grade
#' directions of interest to include (e.g. `"H"` (high), `"L"` (low), or `"B"` (both)).
#'
#' For the direction(s) specified in `direction_var`, patient counts by worst grade for patients who have
#' worsened from baseline are calculated as follows:
#'   * `1` to `4`: The number of patients who have worsened from their baseline grades with worst
#'     grades 1-4, respectively.
#'   * `Any`: The total number of patients who have worsened from their baseline grades.
#'
#' Fractions are calculated by dividing the above counts by the number of patients who's analysis toxicity grades
#' have worsened from baseline toxicity grades during treatment.
#'
#' Prior to using this function in your table layout you must use [rtables::split_rows_by()] to create a row
#' split on variable `direction_var`.
#'
#' @inheritParams argument_convention
#' @param variables (named `list` of `string`)\cr list of additional analysis variables including:
#'   * `id` (`string`)\cr subject variable name.
#'   * `baseline_var` (`string`)\cr name of the data column containing baseline toxicity variable.
#'   * `direction_var` (`string`)\cr see `direction_var` for more details.
#' @param .stats (`character`)\cr statistics to select for the table. Run `get_stats("abnormal_by_worst_grade_worsen")`
#'   to see all available statistics.
#'
#' @seealso Relevant helper functions [h_adlb_worsen()] and [h_worsen_counter()] which are used within
#' [s_count_abnormal_lab_worsen_by_baseline()] to process input data.
#'
#' @name abnormal_by_worst_grade_worsen
#' @order 1
NULL

#' Helper function to prepare ADLB with worst labs
#'
#' @description `r lifecycle::badge("stable")`
#'
#' Helper function to prepare a `df` for generate the patient count shift table.
#'
#' @param adlb (`data.frame`)\cr ADLB data frame.
#' @param worst_flag_low (named `vector`)\cr worst low post-baseline lab grade flag variable. See how this is
#'   implemented in the following examples.
#' @param worst_flag_high (named `vector`)\cr worst high post-baseline lab grade flag variable. See how this is
#'   implemented in the following examples.
#' @param direction_var (`string`)\cr name of the direction variable specifying the direction of the shift table of
#'   interest. Only lab records flagged by `L`, `H` or `B` are included in the shift table.
#'   * `L`: low direction only
#'   * `H`: high direction only
#'   * `B`: both low and high directions
#'
#' @return `h_adlb_worsen()` returns the `adlb` `data.frame` containing only the
#'   worst labs specified according to `worst_flag_low` or `worst_flag_high` for the
#'   direction specified according to `direction_var`. For instance, for a lab that is
#'   needed for the low direction only, only records flagged by `worst_flag_low` are
#'   selected. For a lab that is needed for both low and high directions, the worst
#'   low records are selected for the low direction, and the worst high record are selected
#'   for the high direction.
#'
#' @seealso [abnormal_by_worst_grade_worsen]
#'
#' @examples
#' library(dplyr)
#'
#' # The direction variable, GRADDR, is based on metadata
#' adlb <- tern_ex_adlb %>%
#'   mutate(
#'     GRADDR = case_when(
#'       PARAMCD == "ALT" ~ "B",
#'       PARAMCD == "CRP" ~ "L",
#'       PARAMCD == "IGA" ~ "H"
#'     )
#'   ) %>%
#'   filter(SAFFL == "Y" & ONTRTFL == "Y" & GRADDR != "")
#'
#' df <- h_adlb_worsen(
#'   adlb,
#'   worst_flag_low = c("WGRLOFL" = "Y"),
#'   worst_flag_high = c("WGRHIFL" = "Y"),
#'   direction_var = "GRADDR"
#' )
#'
#' @export
h_adlb_worsen <- function(adlb,
                          worst_flag_low = NULL,
                          worst_flag_high = NULL,
                          direction_var) {
  checkmate::assert_string(direction_var)
  checkmate::assert_subset(as.character(unique(adlb[[direction_var]])), c("B", "L", "H"))
  assert_df_with_variables(adlb, list("Col" = direction_var))

  if (any(unique(adlb[[direction_var]]) == "H")) {
    assert_df_with_variables(adlb, list("High" = names(worst_flag_high)))
  }

  if (any(unique(adlb[[direction_var]]) == "L")) {
    assert_df_with_variables(adlb, list("Low" = names(worst_flag_low)))
  }

  if (any(unique(adlb[[direction_var]]) == "B")) {
    assert_df_with_variables(
      adlb,
      list(
        "Low" = names(worst_flag_low),
        "High" = names(worst_flag_high)
      )
    )
  }

  # extract patients with worst post-baseline lab, either low or high or both
  worst_flag <- c(worst_flag_low, worst_flag_high)
  col_names <- names(worst_flag)
  filter_values <- worst_flag
  temp <- Map(
    function(x, y) which(adlb[[x]] == y),
    col_names,
    filter_values
  )
  position_satisfy_filters <- Reduce(union, temp)

  # select variables of interest
  adlb_f <- adlb[position_satisfy_filters, ]

  # generate subsets for different directionality
  adlb_f_h <- adlb_f[which(adlb_f[[direction_var]] == "H"), ]
  adlb_f_l <- adlb_f[which(adlb_f[[direction_var]] == "L"), ]
  adlb_f_b <- adlb_f[which(adlb_f[[direction_var]] == "B"), ]

  # for labs requiring both high and low, data is duplicated and will be stacked on top of each other
  adlb_f_b_h <- adlb_f_b
  adlb_f_b_l <- adlb_f_b

  # extract data with worst lab
  if (!is.null(worst_flag_high) && !is.null(worst_flag_low)) {
    # change H to High, L to Low
    adlb_f_h[[direction_var]] <- rep("High", nrow(adlb_f_h))
    adlb_f_l[[direction_var]] <- rep("Low", nrow(adlb_f_l))

    # change, B to High and Low
    adlb_f_b_h[[direction_var]] <- rep("High", nrow(adlb_f_b_h))
    adlb_f_b_l[[direction_var]] <- rep("Low", nrow(adlb_f_b_l))

    adlb_out_h <- adlb_f_h[which(adlb_f_h[[names(worst_flag_high)]] == worst_flag_high), ]
    adlb_out_b_h <- adlb_f_b_h[which(adlb_f_b_h[[names(worst_flag_high)]] == worst_flag_high), ]
    adlb_out_l <- adlb_f_l[which(adlb_f_l[[names(worst_flag_low)]] == worst_flag_low), ]
    adlb_out_b_l <- adlb_f_b_l[which(adlb_f_b_l[[names(worst_flag_low)]] == worst_flag_low), ]

    out <- rbind(adlb_out_h, adlb_out_b_h, adlb_out_l, adlb_out_b_l)
  } else if (!is.null(worst_flag_high)) {
    adlb_f_h[[direction_var]] <- rep("High", nrow(adlb_f_h))
    adlb_f_b_h[[direction_var]] <- rep("High", nrow(adlb_f_b_h))

    adlb_out_h <- adlb_f_h[which(adlb_f_h[[names(worst_flag_high)]] == worst_flag_high), ]
    adlb_out_b_h <- adlb_f_b_h[which(adlb_f_b_h[[names(worst_flag_high)]] == worst_flag_high), ]

    out <- rbind(adlb_out_h, adlb_out_b_h)
  } else if (!is.null(worst_flag_low)) {
    adlb_f_l[[direction_var]] <- rep("Low", nrow(adlb_f_l))
    adlb_f_b_l[[direction_var]] <- rep("Low", nrow(adlb_f_b_l))

    adlb_out_l <- adlb_f_l[which(adlb_f_l[[names(worst_flag_low)]] == worst_flag_low), ]
    adlb_out_b_l <- adlb_f_b_l[which(adlb_f_b_l[[names(worst_flag_low)]] == worst_flag_low), ]

    out <- rbind(adlb_out_l, adlb_out_b_l)
  }

  # label
  formatters::var_labels(out) <- formatters::var_labels(adlb_f, fill = FALSE)
  # NA
  out
}

#' Helper function to analyze patients for `s_count_abnormal_lab_worsen_by_baseline()`
#'
#' @description `r lifecycle::badge("stable")`
#'
#' Helper function to count the number of patients and the fraction of patients according to
#' highest post-baseline lab grade variable `.var`, baseline lab grade variable `baseline_var`,
#' and the direction of interest specified in `direction_var`.
#'
#' @inheritParams argument_convention
#' @inheritParams h_adlb_worsen
#' @param baseline_var (`string`)\cr name of the baseline lab grade variable.
#'
#' @return The counts and fraction of patients
#'   whose worst post-baseline lab grades are worse than their baseline grades, for
#'   post-baseline worst grades "1", "2", "3", "4" and "Any".
#'
#' @seealso [abnormal_by_worst_grade_worsen]
#'
#' @examples
#' library(dplyr)
#'
#' # The direction variable, GRADDR, is based on metadata
#' adlb <- tern_ex_adlb %>%
#'   mutate(
#'     GRADDR = case_when(
#'       PARAMCD == "ALT" ~ "B",
#'       PARAMCD == "CRP" ~ "L",
#'       PARAMCD == "IGA" ~ "H"
#'     )
#'   ) %>%
#'   filter(SAFFL == "Y" & ONTRTFL == "Y" & GRADDR != "")
#'
#' df <- h_adlb_worsen(
#'   adlb,
#'   worst_flag_low = c("WGRLOFL" = "Y"),
#'   worst_flag_high = c("WGRHIFL" = "Y"),
#'   direction_var = "GRADDR"
#' )
#'
#' # `h_worsen_counter`
#' h_worsen_counter(
#'   df %>% filter(PARAMCD == "CRP" & GRADDR == "Low"),
#'   id = "USUBJID",
#'   .var = "ATOXGR",
#'   baseline_var = "BTOXGR",
#'   direction_var = "GRADDR"
#' )
#'
#' @export
h_worsen_counter <- function(df, id, .var, baseline_var, direction_var) {
  checkmate::assert_string(id)
  checkmate::assert_string(.var)
  checkmate::assert_string(baseline_var)
  checkmate::assert_scalar(unique(df[[direction_var]]))
  checkmate::assert_subset(unique(df[[direction_var]]), c("High", "Low"))
  assert_df_with_variables(df, list(val = c(id, .var, baseline_var, direction_var)))

  # remove post-baseline missing
  df <- df[df[[.var]] != "<Missing>", ]

  # obtain directionality
  direction <- unique(df[[direction_var]])

  if (direction == "Low") {
    grade <- -1:-4
    worst_grade <- -4
  } else if (direction == "High") {
    grade <- 1:4
    worst_grade <- 4
  }

  if (nrow(df) > 0) {
    by_grade <- lapply(grade, function(i) {
      # filter baseline values that is less than i or <Missing>
      df_temp <- df[df[[baseline_var]] %in% c((i + sign(i) * -1):(-1 * worst_grade), "<Missing>"), ]
      # num: number of patients with post-baseline worst lab equal to i
      num <- length(unique(df_temp[df_temp[[.var]] %in% i, id, drop = TRUE]))
      # denom: number of patients with baseline values less than i or <missing> and post-baseline in the same direction
      denom <- length(unique(df_temp[[id]]))
      rm(df_temp)
      c(num = num, denom = denom)
    })
  } else {
    by_grade <- lapply(1, function(i) {
      c(num = 0, denom = 0)
    })
  }

  names(by_grade) <- as.character(seq_along(by_grade))

  # baseline grade less 4 or missing
  df_temp <- df[!df[[baseline_var]] %in% worst_grade, ]

  # denom: number of patients with baseline values less than 4 or <missing> and post-baseline in the same direction
  denom <- length(unique(df_temp[, id, drop = TRUE]))

  # condition 1: missing baseline and in the direction of abnormality
  con1 <- which(df_temp[[baseline_var]] == "<Missing>" & df_temp[[.var]] %in% grade)
  df_temp_nm <- df_temp[which(df_temp[[baseline_var]] != "<Missing>" & df_temp[[.var]] %in% grade), ]

  # condition 2: if post-baseline values are present then post-baseline values must be worse than baseline
  if (direction == "Low") {
    con2 <- which(as.numeric(as.character(df_temp_nm[[.var]])) < as.numeric(as.character(df_temp_nm[[baseline_var]])))
  } else {
    con2 <- which(as.numeric(as.character(df_temp_nm[[.var]])) > as.numeric(as.character(df_temp_nm[[baseline_var]])))
  }

  # number of patients satisfy either conditions 1 or 2
  num <- length(unique(df_temp[union(con1, con2), id, drop = TRUE]))

  list(fraction = c(by_grade, list("Any" = c(num = num, denom = denom))))
}

#' @describeIn abnormal_by_worst_grade_worsen Statistics function for patients whose worst post-baseline
#'   lab grades are worse than their baseline grades.
#'
#' @return
#' * `s_count_abnormal_lab_worsen_by_baseline()` returns the counts and fraction of patients whose worst
#'   post-baseline lab grades are worse than their baseline grades, for post-baseline worst grades
#'   "1", "2", "3", "4" and "Any".
#'
#' @keywords internal
s_count_abnormal_lab_worsen_by_baseline <- function(df, # nolint
                                                    .var = "ATOXGR",
                                                    variables = list(
                                                      id = "USUBJID",
                                                      baseline_var = "BTOXGR",
                                                      direction_var = "GRADDR"
                                                    )) {
  checkmate::assert_string(.var)
  checkmate::assert_set_equal(names(variables), c("id", "baseline_var", "direction_var"))
  checkmate::assert_string(variables$id)
  checkmate::assert_string(variables$baseline_var)
  checkmate::assert_string(variables$direction_var)
  assert_df_with_variables(df, c(aval = .var, variables[1:3]))
  assert_list_of_variables(variables)

  h_worsen_counter(df, variables$id, .var, variables$baseline_var, variables$direction_var)
}

#' @describeIn abnormal_by_worst_grade_worsen Formatted analysis function which is used as `afun`
#'   in `count_abnormal_lab_worsen_by_baseline()`.
#'
#' @return
#' * `a_count_abnormal_lab_worsen_by_baseline()` returns the corresponding list with
#'   formatted [rtables::CellValue()].
#'
#' @keywords internal
a_count_abnormal_lab_worsen_by_baseline <- make_afun( # nolint
  s_count_abnormal_lab_worsen_by_baseline,
  .formats = c(fraction = format_fraction),
  .ungroup_stats = "fraction"
)

#' @describeIn abnormal_by_worst_grade_worsen Layout-creating function which can take statistics function
#'   arguments and additional format arguments. This function is a wrapper for [rtables::analyze()].
#'
#' @return
#' * `count_abnormal_lab_worsen_by_baseline()` 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_count_abnormal_lab_worsen_by_baseline()` to the table layout.
#'
#' @examples
#' library(dplyr)
#'
#' # The direction variable, GRADDR, is based on metadata
#' adlb <- tern_ex_adlb %>%
#'   mutate(
#'     GRADDR = case_when(
#'       PARAMCD == "ALT" ~ "B",
#'       PARAMCD == "CRP" ~ "L",
#'       PARAMCD == "IGA" ~ "H"
#'     )
#'   ) %>%
#'   filter(SAFFL == "Y" & ONTRTFL == "Y" & GRADDR != "")
#'
#' df <- h_adlb_worsen(
#'   adlb,
#'   worst_flag_low = c("WGRLOFL" = "Y"),
#'   worst_flag_high = c("WGRHIFL" = "Y"),
#'   direction_var = "GRADDR"
#' )
#'
#' basic_table() %>%
#'   split_cols_by("ARMCD") %>%
#'   add_colcounts() %>%
#'   split_rows_by("PARAMCD") %>%
#'   split_rows_by("GRADDR") %>%
#'   count_abnormal_lab_worsen_by_baseline(
#'     var = "ATOXGR",
#'     variables = list(
#'       id = "USUBJID",
#'       baseline_var = "BTOXGR",
#'       direction_var = "GRADDR"
#'     )
#'   ) %>%
#'   append_topleft("Direction of Abnormality") %>%
#'   build_table(df = df, alt_counts_df = tern_ex_adsl)
#'
#' @export
#' @order 2
count_abnormal_lab_worsen_by_baseline <- function(lyt, # nolint
                                                  var,
                                                  variables = list(
                                                    id = "USUBJID",
                                                    baseline_var = "BTOXGR",
                                                    direction_var = "GRADDR"
                                                  ),
                                                  na_str = default_na_str(),
                                                  nested = TRUE,
                                                  ...,
                                                  table_names = NULL,
                                                  .stats = NULL,
                                                  .formats = NULL,
                                                  .labels = NULL,
                                                  .indent_mods = NULL) {
  checkmate::assert_string(var)

  extra_args <- list(variables = variables, ...)

  afun <- make_afun(
    a_count_abnormal_lab_worsen_by_baseline,
    .stats = .stats,
    .formats = .formats,
    .labels = .labels,
    .indent_mods = .indent_mods
  )

  lyt <- analyze(
    lyt = lyt,
    vars = var,
    afun = afun,
    na_str = na_str,
    nested = nested,
    extra_args = extra_args,
    show_labels = "hidden"
  )

  lyt
}

#' Confidence intervals for a difference of binomials
#'
#' @description `r lifecycle::badge("experimental")`
#'
#' Several confidence intervals for the difference between proportions.
#'
#' @name desctools_binom
NULL

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

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

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

#' Create a forest plot from an `rtable`
#'
#' Given a [rtables::rtable()] object with at least one column with a single value and one column with 2
#' values, converts table to a [ggplot2::ggplot()] object and generates an accompanying forest plot. The
#' table and forest plot are printed side-by-side.
#'
#' @description `r lifecycle::badge("stable")`
#'
#' @inheritParams rtable2gg
#' @inheritParams argument_convention
#' @param tbl (`VTableTree`)\cr `rtables` table with at least one column with a single value and one column with 2
#'   values.
#' @param col_x (`integer(1)` or `NULL`)\cr column index with estimator. By default tries to get this from
#'   `tbl` attribute `col_x`, otherwise needs to be manually specified. If `NULL`, points will be excluded
#'   from forest plot.
#' @param col_ci (`integer(1)` or `NULL`)\cr column index with confidence intervals. By default tries to get this from
#'   `tbl` attribute `col_ci`, otherwise needs to be manually specified. If `NULL`, lines will be excluded
#'   from forest plot.
#' @param vline (`numeric(1)` or `NULL`)\cr x coordinate for vertical line, if `NULL` then the line is omitted.
#' @param forest_header (`character(2)`)\cr text displayed to the left and right of `vline`, respectively.
#'   If `vline = NULL` then `forest_header` is not printed. By default tries to get this from `tbl` attribute
#'   `forest_header`. If `NULL`, defaults will be extracted from the table if possible, and set to
#'   `"Comparison\nBetter"` and `"Treatment\nBetter"` if not.
#' @param xlim (`numeric(2)`)\cr limits for x axis.
#' @param logx (`flag`)\cr show the x-values on logarithm scale.
#' @param x_at (`numeric`)\cr x-tick locations, if `NULL`, `x_at` is set to `vline` and both `xlim` values.
#' @param width_row_names `r lifecycle::badge("deprecated")` Please use the `lbl_col_padding` argument instead.
#' @param width_columns (`numeric`)\cr a vector of column widths. Each element's position in
#'   `colwidths` corresponds to the column of `tbl` in the same position. If `NULL`, column widths are calculated
#'   according to maximum number of characters per column.
#' @param width_forest `r lifecycle::badge("deprecated")` Please use the `rel_width_forest` argument instead.
#' @param rel_width_forest (`proportion`)\cr proportion of total width to allocate to the forest plot. Relative
#'   width of table is then `1 - rel_width_forest`. If `as_list = TRUE`, this parameter is ignored.
#' @param font_size (`numeric(1)`)\cr font size.
#' @param col_symbol_size (`numeric` or `NULL`)\cr column index from `tbl` containing data to be used
#'   to determine relative size for estimator plot symbol. Typically, the symbol size is proportional
#'   to the sample size used to calculate the estimator. If `NULL`, the same symbol size is used for all subgroups.
#'   By default tries to get this from `tbl` attribute `col_symbol_size`, otherwise needs to be manually specified.
#' @param col (`character`)\cr color(s).
#' @param ggtheme (`theme`)\cr a graphical theme as provided by `ggplot2` to control styling of the plot.
#' @param as_list (`flag`)\cr whether the two `ggplot` objects should be returned as a list. If `TRUE`, a named list
#'   with two elements, `table` and `plot`, will be returned. If `FALSE` (default) the table and forest plot are
#'   printed side-by-side via [cowplot::plot_grid()].
#' @param gp `r lifecycle::badge("deprecated")` `g_forest` is now generated as a `ggplot` object. This argument
#'   is no longer used.
#' @param draw `r lifecycle::badge("deprecated")` `g_forest` is now generated as a `ggplot` object. This argument
#'   is no longer used.
#' @param newpage `r lifecycle::badge("deprecated")` `g_forest` is now generated as a `ggplot` object. This argument
#'   is no longer used.
#'
#' @return `ggplot` forest plot and table.
#'
#' @examples
#' library(dplyr)
#' library(forcats)
#' library(nestcolor)
#'
#' adrs <- tern_ex_adrs
#' n_records <- 20
#' adrs_labels <- formatters::var_labels(adrs, fill = TRUE)
#' adrs <- adrs %>%
#'   filter(PARAMCD == "BESRSPI") %>%
#'   filter(ARM %in% c("A: Drug X", "B: Placebo")) %>%
#'   slice(seq_len(n_records)) %>%
#'   droplevels() %>%
#'   mutate(
#'     # Reorder levels of factor to make the placebo group the reference arm.
#'     ARM = fct_relevel(ARM, "B: Placebo"),
#'     rsp = AVALC == "CR"
#'   )
#' formatters::var_labels(adrs) <- c(adrs_labels, "Response")
#' df <- extract_rsp_subgroups(
#'   variables = list(rsp = "rsp", arm = "ARM", subgroups = c("SEX", "STRATA2")),
#'   data = adrs
#' )
#' # Full commonly used response table.
#'
#' tbl <- basic_table() %>%
#'   tabulate_rsp_subgroups(df)
#' g_forest(tbl)
#'
#' # Odds ratio only table.
#'
#' tbl_or <- basic_table() %>%
#'   tabulate_rsp_subgroups(df, vars = c("n_tot", "or", "ci"))
#' g_forest(
#'   tbl_or,
#'   forest_header = c("Comparison\nBetter", "Treatment\nBetter")
#' )
#'
#' # Survival forest plot example.
#' adtte <- tern_ex_adtte
#' # Save variable labels before data processing steps.
#' adtte_labels <- formatters::var_labels(adtte, fill = TRUE)
#' adtte_f <- adtte %>%
#'   filter(
#'     PARAMCD == "OS",
#'     ARM %in% c("B: Placebo", "A: Drug X"),
#'     SEX %in% c("M", "F")
#'   ) %>%
#'   mutate(
#'     # Reorder levels of ARM to display reference arm before treatment arm.
#'     ARM = droplevels(fct_relevel(ARM, "B: Placebo")),
#'     SEX = droplevels(SEX),
#'     AVALU = as.character(AVALU),
#'     is_event = CNSR == 0
#'   )
#' labels <- list(
#'   "ARM" = adtte_labels["ARM"],
#'   "SEX" = adtte_labels["SEX"],
#'   "AVALU" = adtte_labels["AVALU"],
#'   "is_event" = "Event Flag"
#' )
#' formatters::var_labels(adtte_f)[names(labels)] <- as.character(labels)
#' df <- extract_survival_subgroups(
#'   variables = list(
#'     tte = "AVAL",
#'     is_event = "is_event",
#'     arm = "ARM", subgroups = c("SEX", "BMRKR2")
#'   ),
#'   data = adtte_f
#' )
#' table_hr <- basic_table() %>%
#'   tabulate_survival_subgroups(df, time_unit = adtte_f$AVALU[1])
#' g_forest(table_hr)
#'
#' # Works with any `rtable`.
#' tbl <- rtable(
#'   header = c("E", "CI", "N"),
#'   rrow("", 1, c(.8, 1.2), 200),
#'   rrow("", 1.2, c(1.1, 1.4), 50)
#' )
#' g_forest(
#'   tbl = tbl,
#'   col_x = 1,
#'   col_ci = 2,
#'   xlim = c(0.5, 2),
#'   x_at = c(0.5, 1, 2),
#'   col_symbol_size = 3
#' )
#'
#' tbl <- rtable(
#'   header = rheader(
#'     rrow("", rcell("A", colspan = 2)),
#'     rrow("", "c1", "c2")
#'   ),
#'   rrow("row 1", 1, c(.8, 1.2)),
#'   rrow("row 2", 1.2, c(1.1, 1.4))
#' )
#' g_forest(
#'   tbl = tbl,
#'   col_x = 1,
#'   col_ci = 2,
#'   xlim = c(0.5, 2),
#'   x_at = c(0.5, 1, 2),
#'   vline = 1,
#'   forest_header = c("Hello", "World")
#' )
#'
#' @export
g_forest <- function(tbl,
                     col_x = attr(tbl, "col_x"),
                     col_ci = attr(tbl, "col_ci"),
                     vline = 1,
                     forest_header = attr(tbl, "forest_header"),
                     xlim = c(0.1, 10),
                     logx = TRUE,
                     x_at = c(0.1, 1, 10),
                     width_row_names = lifecycle::deprecated(),
                     width_columns = NULL,
                     width_forest = lifecycle::deprecated(),
                     lbl_col_padding = 0,
                     rel_width_forest = 0.25,
                     font_size = 12,
                     col_symbol_size = attr(tbl, "col_symbol_size"),
                     col = getOption("ggplot2.discrete.colour")[1],
                     ggtheme = NULL,
                     as_list = FALSE,
                     gp = lifecycle::deprecated(),
                     draw = lifecycle::deprecated(),
                     newpage = lifecycle::deprecated()) {
  # Deprecated argument warnings
  if (lifecycle::is_present(width_row_names)) {
    lifecycle::deprecate_warn(
      "0.9.4", "g_forest(width_row_names)", "g_forest(lbl_col_padding)",
      details = "The width of the row label column can be adjusted via the `lbl_col_padding` parameter."
    )
  }
  if (lifecycle::is_present(width_forest)) {
    lifecycle::deprecate_warn(
      "0.9.4", "g_forest(width_forest)", "g_forest(rel_width_forest)",
      details = "Relative width of the forest plot (as a proportion) can be set via the `rel_width_forest` parameter."
    )
  }
  if (lifecycle::is_present(gp)) {
    lifecycle::deprecate_warn(
      "0.9.4", "g_forest(gp)", "g_forest(ggtheme)",
      details = paste(
        "`g_forest` is now generated as a `ggplot` object.",
        "Additional display settings should be supplied via the `ggtheme` parameter."
      )
    )
  }
  if (lifecycle::is_present(draw)) {
    lifecycle::deprecate_warn(
      "0.9.4", "g_forest(draw)",
      details = "`g_forest` now generates `ggplot` objects. This parameter has no effect."
    )
  }
  if (lifecycle::is_present(newpage)) {
    lifecycle::deprecate_warn(
      "0.9.4", "g_forest(newpage)",
      details = "`g_forest` now generates `ggplot` objects. This parameter has no effect."
    )
  }

  checkmate::assert_class(tbl, "VTableTree")
  checkmate::assert_number(col_x, lower = 0, upper = ncol(tbl), null.ok = TRUE)
  checkmate::assert_number(col_ci, lower = 0, upper = ncol(tbl), null.ok = TRUE)
  checkmate::assert_number(col_symbol_size, lower = 0, upper = ncol(tbl), null.ok = TRUE)
  checkmate::assert_number(font_size, lower = 0)
  checkmate::assert_character(col, null.ok = TRUE)
  checkmate::assert_true(is.null(col) | length(col) == 1 | length(col) == nrow(tbl))

  # Extract info from table
  mat <- matrix_form(tbl, indent_rownames = TRUE)
  mat_strings <- formatters::mf_strings(mat)
  nlines_hdr <- formatters::mf_nlheader(mat)
  nrows_body <- nrow(mat_strings) - nlines_hdr
  tbl_stats <- mat_strings[nlines_hdr, -1]

  # Generate and modify table as ggplot object
  gg_table <- rtable2gg(tbl, fontsize = font_size, colwidths = width_columns, lbl_col_padding = lbl_col_padding) +
    theme(plot.margin = margin(0, 0, 0, 0.025, "npc"))
  gg_table$scales$scales[[1]]$expand <- c(0.01, 0.01)
  gg_table$scales$scales[[2]]$limits[2] <- nrow(mat_strings) + 1
  if (nlines_hdr == 2) {
    gg_table$scales$scales[[2]]$expand <- c(0, 0)
    arms <- unique(mat_strings[1, ][nzchar(trimws(mat_strings[1, ]))])
  } else {
    arms <- NULL
  }

  tbl_df <- as_result_df(tbl)
  dat_cols <- seq(which(names(tbl_df) == "node_class") + 1, ncol(tbl_df))
  tbl_df <- tbl_df[, c(which(names(tbl_df) == "row_num"), dat_cols)]
  names(tbl_df) <- c("row_num", tbl_stats)

  # Check table data columns
  if (!is.null(col_ci)) {
    ci_col <- col_ci + 1
  } else {
    tbl_df[["empty_ci"]] <- rep(list(c(NA_real_, NA_real_)), nrow(tbl_df))
    ci_col <- which(names(tbl_df) == "empty_ci")
  }
  if (length(tbl_df[, ci_col][[1]]) != 2) stop("CI column must have two elements (lower and upper limits).")

  if (!is.null(col_x)) {
    x_col <- col_x + 1
  } else {
    tbl_df[["empty_x"]] <- NA_real_
    x_col <- which(names(tbl_df) == "empty_x")
  }
  if (!is.null(col_symbol_size)) {
    sym_size <- unlist(tbl_df[, col_symbol_size + 1])
  } else {
    sym_size <- rep(1, nrow(tbl_df))
  }

  tbl_df[, c("ci_lwr", "ci_upr")] <- t(sapply(tbl_df[, ci_col], unlist))
  x <- unlist(tbl_df[, x_col])
  lwr <- unlist(tbl_df[["ci_lwr"]])
  upr <- unlist(tbl_df[["ci_upr"]])
  row_num <- nrow(mat_strings) - tbl_df[["row_num"]] - as.numeric(nlines_hdr == 2)

  if (is.null(col)) col <- "#343cff"
  if (length(col) == 1) col <- rep(col, nrow(tbl_df))
  if (is.null(x_at)) x_at <- union(xlim, vline)
  x_labels <- x_at

  # Apply log transformation
  if (logx) {
    x_t <- log(x)
    lwr_t <- log(lwr)
    upr_t <- log(upr)
    xlim_t <- log(xlim)
  } else {
    x_t <- x
    lwr_t <- lwr
    upr_t <- upr
    xlim_t <- xlim
  }

  # Set up plot area
  gg_plt <- ggplot(data = tbl_df) +
    theme(
      panel.background = element_rect(fill = "transparent", color = NA_character_),
      plot.background = element_rect(fill = "transparent", color = NA_character_),
      panel.grid.major = element_blank(),
      panel.grid.minor = element_blank(),
      axis.title.x = element_blank(),
      axis.title.y = element_blank(),
      axis.line.x = element_line(),
      axis.text = element_text(size = font_size),
      legend.position = "none",
      plot.margin = margin(0, 0.1, 0.05, 0, "npc")
    ) +
    scale_x_continuous(
      transform = ifelse(logx, "log", "identity"),
      limits = xlim,
      breaks = x_at,
      labels = x_labels,
      expand = c(0.01, 0)
    ) +
    scale_y_continuous(
      limits = c(0, nrow(mat_strings) + 1),
      breaks = NULL,
      expand = c(0, 0)
    ) +
    coord_cartesian(clip = "off")

  if (is.null(ggtheme)) {
    gg_plt <- gg_plt + annotate(
      "rect",
      xmin = xlim[1],
      xmax = xlim[2],
      ymin = 0,
      ymax = nrows_body + 0.5,
      fill = "grey92"
    )
  }

  if (!is.null(vline)) {
    # Set default forest header
    if (is.null(forest_header)) {
      forest_header <- c(
        paste(if (length(arms) == 2) arms[1] else "Comparison", "Better", sep = "\n"),
        paste(if (length(arms) == 2) arms[2] else "Treatment", "Better", sep = "\n")
      )
    }

    # Add vline and forest header labels
    mid_pts <- if (logx) {
      c(exp(mean(log(c(xlim[1], vline)))), exp(mean(log(c(vline, xlim[2])))))
    } else {
      c(mean(c(xlim[1], vline)), mean(c(vline, xlim[2])))
    }
    gg_plt <- gg_plt +
      annotate(
        "segment",
        x = vline, xend = vline, y = 0, yend = nrows_body + 0.5
      ) +
      annotate(
        "text",
        x = mid_pts[1], y = nrows_body + 1.25,
        label = forest_header[1],
        size = font_size / .pt,
        lineheight = 0.9
      ) +
      annotate(
        "text",
        x = mid_pts[2], y = nrows_body + 1.25,
        label = forest_header[2],
        size = font_size / .pt,
        lineheight = 0.9
      )
  }

  # Add points to plot
  if (any(!is.na(x_t))) {
    x_t[x < xlim[1] | x > xlim[2]] <- NA
    gg_plt <- gg_plt + geom_point(
      x = x_t,
      y = row_num,
      color = col,
      aes(size = sym_size),
      na.rm = TRUE
    )
  }

  for (i in seq_len(nrow(tbl_df))) {
    # Determine which arrow(s) to add to CI lines
    which_arrow <- c(lwr_t[i] < xlim_t[1], upr_t[i] > xlim_t[2])
    which_arrow <- dplyr::case_when(
      all(which_arrow) ~ "both",
      which_arrow[1] ~ "first",
      which_arrow[2] ~ "last",
      TRUE ~ NA_character_
    )

    # Add CI lines
    gg_plt <- gg_plt +
      if (!is.na(which_arrow)) {
        annotate(
          "segment",
          x = if (!which_arrow %in% c("first", "both")) lwr[i] else xlim[1],
          xend = if (!which_arrow %in% c("last", "both")) upr[i] else xlim[2],
          y = row_num[i], yend = row_num[i],
          color = if (length(col) == 1) col else col[i],
          arrow = arrow(length = unit(0.05, "npc"), ends = which_arrow),
          na.rm = TRUE
        )
      } else {
        annotate(
          "segment",
          x = lwr[i], xend = upr[i],
          y = row_num[i], yend = row_num[i],
          color = if (length(col) == 1) col else col[i],
          na.rm = TRUE
        )
      }
  }

  # Apply custom ggtheme to plot
  if (!is.null(ggtheme)) gg_plt <- gg_plt + ggtheme

  if (as_list) {
    list(
      table = gg_table,
      plot = gg_plt
    )
  } else {
    cowplot::plot_grid(
      gg_table,
      gg_plt,
      align = "h",
      axis = "tblr",
      rel_widths = c(1 - rel_width_forest, rel_width_forest)
    )
  }
}

#' Forest plot grob
#'
#' @description `r lifecycle::badge("deprecated")`
#'
#' @inheritParams g_forest
#' @param tbl (`VTableTree`)\cr `rtables` table object.
#' @param x (`numeric`)\cr coordinate of point.
#' @param lower,upper (`numeric`)\cr lower/upper bound of the confidence interval.
#' @param symbol_size (`numeric`)\cr vector with relative size for plot symbol.
#'   If `NULL`, the same symbol size is used.
#'
#' @details
#' The heights get automatically determined.
#'
#' @examples
#' tbl <- rtable(
#'   header = rheader(
#'     rrow("", "E", rcell("CI", colspan = 2), "N"),
#'     rrow("", "A", "B", "C", "D")
#'   ),
#'   rrow("row 1", 1, 0.8, 1.1, 16),
#'   rrow("row 2", 1.4, 0.8, 1.6, 25),
#'   rrow("row 3", 1.2, 0.8, 1.6, 36)
#' )
#'
#' x <- c(1, 1.4, 1.2)
#' lower <- c(0.8, 0.8, 0.8)
#' upper <- c(1.1, 1.6, 1.6)
#' # numeric vector with multiplication factor to scale each circle radius
#' # default radius is 1/3.5 lines
#' symbol_scale <- c(1, 1.25, 1.5)
#'
#' # Internal function - forest_grob
#' \donttest{
#' p <- forest_grob(tbl, x, lower, upper,
#'   vline = 1, forest_header = c("A", "B"),
#'   x_at = c(.1, 1, 10), xlim = c(0.1, 10), logx = TRUE, symbol_size = symbol_scale,
#'   vp = grid::plotViewport(margins = c(1, 1, 1, 1))
#' )
#'
#' draw_grob(p)
#' }
#'
#' @noRd
#' @keywords internal
forest_grob <- function(tbl,
                        x,
                        lower,
                        upper,
                        vline,
                        forest_header,
                        xlim = NULL,
                        logx = FALSE,
                        x_at = NULL,
                        width_row_names = NULL,
                        width_columns = NULL,
                        width_forest = grid::unit(1, "null"),
                        symbol_size = NULL,
                        col = "blue",
                        name = NULL,
                        gp = NULL,
                        vp = NULL) {
  lifecycle::deprecate_warn(
    "0.9.4", "forest_grob()",
    details = "`g_forest` now generates `ggplot` objects. This function is no longer used within `tern`."
  )

  nr <- nrow(tbl)
  if (is.null(vline)) {
    checkmate::assert_true(is.null(forest_header))
  } else {
    checkmate::assert_number(vline)
    checkmate::assert_character(forest_header, len = 2, null.ok = TRUE)
  }

  checkmate::assert_numeric(x, len = nr)
  checkmate::assert_numeric(lower, len = nr)
  checkmate::assert_numeric(upper, len = nr)
  checkmate::assert_numeric(symbol_size, len = nr, null.ok = TRUE)
  checkmate::assert_character(col)

  if (is.null(symbol_size)) {
    symbol_size <- rep(1, nr)
  }

  if (is.null(xlim)) {
    r <- range(c(x, lower, upper), na.rm = TRUE)
    xlim <- r + c(-0.05, 0.05) * diff(r)
  }

  if (logx) {
    if (is.null(x_at)) {
      x_at <- pretty(log(stats::na.omit(c(x, lower, upper))))
      x_labels <- exp(x_at)
    } else {
      x_labels <- x_at
      x_at <- log(x_at)
    }
    xlim <- log(xlim)
    x <- log(x)
    lower <- log(lower)
    upper <- log(upper)
    if (!is.null(vline)) {
      vline <- log(vline)
    }
  } else {
    x_labels <- TRUE
  }

  data_forest_vp <- grid::dataViewport(xlim, c(0, 1))

  # Get table content as matrix form.
  mf <- matrix_form(tbl)

  # Use `rtables` indent_string eventually.
  mf$strings[, 1] <- paste0(
    strrep("    ", c(rep(0, attr(mf, "nrow_header")), mf$row_info$indent)),
    mf$strings[, 1]
  )

  n_header <- attr(mf, "nrow_header")

  if (any(mf$display[, 1] == FALSE)) stop("row names need to be always displayed")

  # Pre-process the data to be used in lapply and cell_in_rows.
  to_args_for_cell_in_rows_fun <- function(part = c("body", "header"),
                                           underline_colspan = FALSE) {
    part <- match.arg(part)
    if (part == "body") {
      mat_row_indices <- seq_len(nrow(tbl)) + n_header
      row_ind_offset <- -n_header
    } else {
      mat_row_indices <- seq_len(n_header)
      row_ind_offset <- 0
    }

    lapply(mat_row_indices, function(i) {
      disp <- mf$display[i, -1]
      list(
        row_name = mf$strings[i, 1],
        cells = mf$strings[i, -1][disp],
        cell_spans = mf$spans[i, -1][disp],
        row_index = i + row_ind_offset,
        underline_colspan = underline_colspan
      )
    })
  }

  args_header <- to_args_for_cell_in_rows_fun("header", underline_colspan = TRUE)
  args_body <- to_args_for_cell_in_rows_fun("body", underline_colspan = FALSE)

  grid::gTree(
    name = name,
    children = grid::gList(
      grid::gTree(
        children = do.call(grid::gList, lapply(args_header, do.call, what = cell_in_rows)),
        vp = grid::vpPath("vp_table_layout", "vp_header")
      ),
      grid::gTree(
        children = do.call(grid::gList, lapply(args_body, do.call, what = cell_in_rows)),
        vp = grid::vpPath("vp_table_layout", "vp_body")
      ),
      grid::linesGrob(
        grid::unit(c(0, 1), "npc"),
        y = grid::unit(c(.5, .5), "npc"),
        vp = grid::vpPath("vp_table_layout", "vp_spacer")
      ),
      # forest part
      if (is.null(vline)) {
        NULL
      } else {
        grid::gTree(
          children = grid::gList(
            grid::gTree(
              children = grid::gList(
                # this may overflow, to fix, look here
                # https://stackoverflow.com/questions/33623169/add-multi-line-footnote-to-tablegrob-while-using-gridextra-in-r # nolint
                grid::textGrob(
                  forest_header[1],
                  x = grid::unit(vline, "native") - grid::unit(1, "lines"),
                  just = c("right", "center")
                ),
                grid::textGrob(
                  forest_header[2],
                  x = grid::unit(vline, "native") + grid::unit(1, "lines"),
                  just = c("left", "center")
                )
              ),
              vp = grid::vpStack(grid::viewport(layout.pos.col = ncol(tbl) + 2), data_forest_vp)
            )
          ),
          vp = grid::vpPath("vp_table_layout", "vp_header")
        )
      },
      grid::gTree(
        children = grid::gList(
          grid::gTree(
            children = grid::gList(
              grid::rectGrob(gp = grid::gpar(col = "gray90", fill = "gray90")),
              if (is.null(vline)) {
                NULL
              } else {
                grid::linesGrob(
                  x = grid::unit(rep(vline, 2), "native"),
                  y = grid::unit(c(0, 1), "npc"),
                  gp = grid::gpar(lwd = 2),
                  vp = data_forest_vp
                )
              },
              grid::xaxisGrob(at = x_at, label = x_labels, vp = data_forest_vp)
            ),
            vp = grid::viewport(layout.pos.col = ncol(tbl) + 2)
          )
        ),
        vp = grid::vpPath("vp_table_layout", "vp_body")
      ),
      grid::gTree(
        children = do.call(
          grid::gList,
          Map(
            function(xi, li, ui, row_index, size_i, col) {
              forest_dot_line(
                xi,
                li,
                ui,
                row_index,
                xlim,
                symbol_size = size_i,
                col = col,
                datavp = data_forest_vp
              )
            },
            x,
            lower,
            upper,
            seq_along(x),
            symbol_size,
            col,
            USE.NAMES = FALSE
          )
        ),
        vp = grid::vpPath("vp_table_layout", "vp_body")
      )
    ),
    childrenvp = forest_viewport(tbl, width_row_names, width_columns, width_forest),
    vp = vp,
    gp = gp
  )
}

cell_in_rows <- function(row_name,
                         cells,
                         cell_spans,
                         row_index,
                         underline_colspan = FALSE) {
  checkmate::assert_string(row_name)
  checkmate::assert_character(cells, min.len = 1, any.missing = FALSE)
  checkmate::assert_numeric(cell_spans, len = length(cells), any.missing = FALSE)
  checkmate::assert_number(row_index)
  checkmate::assert_flag(underline_colspan)

  vp_name_rn <- paste0("rowname-", row_index)
  g_rowname <- if (!is.null(row_name) && row_name != "") {
    grid::textGrob(
      name = vp_name_rn,
      label = row_name,
      x = grid::unit(0, "npc"),
      just = c("left", "center"),
      vp = grid::vpPath(paste0("rowname-", row_index))
    )
  } else {
    NULL
  }

  gl_cols <- if (!(length(cells) > 0)) {
    list(NULL)
  } else {
    j <- 1 # column index of cell

    lapply(seq_along(cells), function(k) {
      cell_ascii <- cells[[k]]
      cs <- cell_spans[[k]]

      if (is.na(cell_ascii) || is.null(cell_ascii)) {
        cell_ascii <- "NA"
      }

      cell_name <- paste0("g-cell-", row_index, "-", j)

      cell_grobs <- if (identical(cell_ascii, "")) {
        NULL
      } else {
        if (cs == 1) {
          grid::textGrob(
            label = cell_ascii,
            name = cell_name,
            vp = grid::vpPath(paste0("cell-", row_index, "-", j))
          )
        } else {
          # +1 because of rowname
          vp_joined_cols <- grid::viewport(layout.pos.row = row_index, layout.pos.col = seq(j + 1, j + cs))

          lab <- grid::textGrob(
            label = cell_ascii,
            name = cell_name,
            vp = vp_joined_cols
          )

          if (!underline_colspan || grepl("^[[:space:]]*$", cell_ascii)) {
            lab
          } else {
            grid::gList(
              lab,
              grid::linesGrob(
                x = grid::unit.c(grid::unit(.2, "lines"), grid::unit(1, "npc") - grid::unit(.2, "lines")),
                y = grid::unit(c(0, 0), "npc"),
                vp = vp_joined_cols
              )
            )
          }
        }
      }
      j <<- j + cs

      cell_grobs
    })
  }

  grid::gList(
    g_rowname,
    do.call(grid::gList, gl_cols)
  )
}

#' Graphic object: forest dot line
#'
#' @description `r lifecycle::badge("deprecated")`
#'
#' Calculate the `grob` corresponding to the dot line within the forest plot.
#'
#' @noRd
#' @keywords internal
forest_dot_line <- function(x,
                            lower,
                            upper,
                            row_index,
                            xlim,
                            symbol_size = 1,
                            col = "blue",
                            datavp) {
  lifecycle::deprecate_warn(
    "0.9.4", "forest_dot_line()",
    details = "`g_forest` now generates `ggplot` objects. This function is no longer used within `tern`."
  )

  ci <- c(lower, upper)
  if (any(!is.na(c(x, ci)))) {
    # line
    y <- grid::unit(c(0.5, 0.5), "npc")

    g_line <- if (all(!is.na(ci)) && ci[2] > xlim[1] && ci[1] < xlim[2]) {
      # -
      if (ci[1] >= xlim[1] && ci[2] <= xlim[2]) {
        grid::linesGrob(x = grid::unit(c(ci[1], ci[2]), "native"), y = y)
      } else if (ci[1] < xlim[1] && ci[2] > xlim[2]) {
        # <->
        grid::linesGrob(
          x = grid::unit(xlim, "native"),
          y = y,
          arrow = grid::arrow(angle = 30, length = grid::unit(0.5, "lines"), ends = "both")
        )
      } else if (ci[1] < xlim[1] && ci[2] <= xlim[2]) {
        # <-
        grid::linesGrob(
          x = grid::unit(c(xlim[1], ci[2]), "native"),
          y = y,
          arrow = grid::arrow(angle = 30, length = grid::unit(0.5, "lines"), ends = "first")
        )
      } else if (ci[1] >= xlim[1] && ci[2] > xlim[2]) {
        # ->
        grid::linesGrob(
          x = grid::unit(c(ci[1], xlim[2]), "native"),
          y = y,
          arrow = grid::arrow(angle = 30, length = grid::unit(0.5, "lines"), ends = "last")
        )
      }
    } else {
      NULL
    }

    g_circle <- if (!is.na(x) && x >= xlim[1] && x <= xlim[2]) {
      grid::circleGrob(
        x = grid::unit(x, "native"),
        y = y,
        r = grid::unit(1 / 3.5 * symbol_size, "lines"),
        name = "point"
      )
    } else {
      NULL
    }

    grid::gTree(
      children = grid::gList(
        grid::gTree(
          children = grid::gList(
            grid::gList(
              g_line,
              g_circle
            )
          ),
          vp = datavp,
          gp = grid::gpar(col = col, fill = col)
        )
      ),
      vp = grid::vpPath(paste0("forest-", row_index))
    )
  } else {
    NULL
  }
}

#' Create a viewport tree for the forest plot
#'
#' @description `r lifecycle::badge("deprecated")`
#'
#' @param tbl (`VTableTree`)\cr `rtables` table object.
#' @param width_row_names (`grid::unit`)\cr width of row names.
#' @param width_columns (`grid::unit`)\cr width of column spans.
#' @param width_forest (`grid::unit`)\cr width of the forest plot.
#' @param gap_column (`grid::unit`)\cr gap width between the columns.
#' @param gap_header (`grid::unit`)\cr gap width between the header.
#' @param mat_form (`MatrixPrintForm`)\cr matrix print form of the table.
#'
#' @return A viewport tree.
#'
#' @examples
#' library(grid)
#'
#' tbl <- rtable(
#'   header = rheader(
#'     rrow("", "E", rcell("CI", colspan = 2)),
#'     rrow("", "A", "B", "C")
#'   ),
#'   rrow("row 1", 1, 0.8, 1.1),
#'   rrow("row 2", 1.4, 0.8, 1.6),
#'   rrow("row 3", 1.2, 0.8, 1.2)
#' )
#'
#' \donttest{
#' v <- forest_viewport(tbl)
#'
#' grid::grid.newpage()
#' showViewport(v)
#' }
#'
#' @export
forest_viewport <- function(tbl,
                            width_row_names = NULL,
                            width_columns = NULL,
                            width_forest = grid::unit(1, "null"),
                            gap_column = grid::unit(1, "lines"),
                            gap_header = grid::unit(1, "lines"),
                            mat_form = NULL) {
  lifecycle::deprecate_warn(
    "0.9.4",
    "forest_viewport()",
    details = "`g_forest` now generates `ggplot` objects. This function is no longer used within `tern`."
  )

  checkmate::assert_class(tbl, "VTableTree")
  checkmate::assert_true(grid::is.unit(width_forest))
  if (!is.null(width_row_names)) {
    checkmate::assert_true(grid::is.unit(width_row_names))
  }
  if (!is.null(width_columns)) {
    checkmate::assert_true(grid::is.unit(width_columns))
  }

  if (is.null(mat_form)) mat_form <- matrix_form(tbl)

  mat_form$strings[!mat_form$display] <- ""

  nr <- nrow(tbl)
  nc <- ncol(tbl)
  nr_h <- attr(mat_form, "nrow_header")

  if (is.null(width_row_names) || is.null(width_columns)) {
    tbl_widths <- formatters::propose_column_widths(mat_form)
    strs_with_width <- strrep("x", tbl_widths) # that works for mono spaced fonts
    if (is.null(width_row_names)) width_row_names <- grid::stringWidth(strs_with_width[1])
    if (is.null(width_columns)) width_columns <- grid::stringWidth(strs_with_width[-1])
  }

  # Widths for row name, cols, forest.
  widths <- grid::unit.c(
    width_row_names + gap_column,
    width_columns + gap_column,
    width_forest
  )

  n_lines_per_row <- apply(
    X = mat_form$strings,
    MARGIN = 1,
    FUN = function(row) {
      tmp <- vapply(
        gregexpr("\n", row, fixed = TRUE),
        attr, numeric(1),
        "match.length"
      ) + 1
      max(c(tmp, 1))
    }
  )

  i_header <- seq_len(nr_h)

  height_body_rows <- grid::unit(n_lines_per_row[-i_header] * 1.2, "lines")
  height_header_rows <- grid::unit(n_lines_per_row[i_header] * 1.2, "lines")

  height_body <- grid::unit(sum(n_lines_per_row[-i_header]) * 1.2, "lines")
  height_header <- grid::unit(sum(n_lines_per_row[i_header]) * 1.2, "lines")

  nc_g <- nc + 2 # number of columns incl. row names and forest

  vp_tbl <- grid::vpTree(
    parent = grid::viewport(
      name = "vp_table_layout",
      layout = grid::grid.layout(
        nrow = 3, ncol = 1,
        heights = grid::unit.c(height_header, gap_header, height_body)
      )
    ),
    children = grid::vpList(
      vp_forest_table_part(nr_h, nc_g, 1, 1, widths, height_header_rows, "vp_header"),
      vp_forest_table_part(nr, nc_g, 3, 1, widths, height_body_rows, "vp_body"),
      grid::viewport(name = "vp_spacer", layout.pos.row = 2, layout.pos.col = 1)
    )
  )
  vp_tbl
}

#' Viewport forest plot: table part
#'
#' @description `r lifecycle::badge("deprecated")`
#'
#' Prepares a viewport for the table included in the forest plot.
#'
#' @noRd
#' @keywords internal
vp_forest_table_part <- function(nrow,
                                 ncol,
                                 l_row,
                                 l_col,
                                 widths,
                                 heights,
                                 name) {
  lifecycle::deprecate_warn(
    "0.9.4", "vp_forest_table_part()",
    details = "`g_forest` now generates `ggplot` objects. This function is no longer used within `tern`."
  )

  grid::vpTree(
    grid::viewport(
      name = name,
      layout.pos.row = l_row,
      layout.pos.col = l_col,
      layout = grid::grid.layout(nrow = nrow, ncol = ncol, widths = widths, heights = heights)
    ),
    children = grid::vpList(
      do.call(
        grid::vpList,
        lapply(
          seq_len(nrow), function(i) {
            grid::viewport(layout.pos.row = i, layout.pos.col = 1, name = paste0("rowname-", i))
          }
        )
      ),
      do.call(
        grid::vpList,
        apply(
          expand.grid(seq_len(nrow), seq_len(ncol - 2)),
          1,
          function(x) {
            i <- x[1]
            j <- x[2]
            grid::viewport(layout.pos.row = i, layout.pos.col = j + 1, name = paste0("cell-", i, "-", j))
          }
        )
      ),
      do.call(
        grid::vpList,
        lapply(
          seq_len(nrow),
          function(i) {
            grid::viewport(layout.pos.row = i, layout.pos.col = ncol, name = paste0("forest-", i))
          }
        )
      )
    )
  )
}

#' Forest rendering
#'
#' @description `r lifecycle::badge("deprecated")`
#'
#' Renders the forest grob.
#'
#' @noRd
#' @keywords internal
grid.forest <- function(...) { # nolint
  lifecycle::deprecate_warn(
    "0.9.4", "grid.forest()",
    details = "`g_forest` now generates `ggplot` objects. This function is no longer used within `tern`."
  )

  grid::grid.draw(forest_grob(...))
}

#' Univariate formula special term
#'
#' @description `r lifecycle::badge("stable")`
#'
#' The special term `univariate` indicate that the model should be fitted individually for
#' every variable included in univariate.
#'
#' @param x (`character`)\cr a vector of variable names separated by commas.
#'
#' @return When used within a model formula, produces univariate models for each variable provided.
#'
#' @details
#' If provided alongside with pairwise specification, the model
#' `y ~ ARM + univariate(SEX, AGE, RACE)` lead to the study and comparison of the models
#' + `y ~ ARM`
#' + `y ~ ARM + SEX`
#' + `y ~ ARM + AGE`
#' + `y ~ ARM + RACE`
#'
#' @export
univariate <- function(x) {
  structure(x, varname = deparse(substitute(x)))
}

# Get the right-hand-term of a formula
rht <- function(x) {
  checkmate::assert_formula(x)
  y <- as.character(rev(x)[[1]])
  return(y)
}

#' Hazard ratio estimation in interactions
#'
#' This function estimates the hazard ratios between arms when an interaction variable is given with
#' specific values.
#'
#' @param variable,given (`character(2)`)\cr names of the two variables in the interaction. We seek the estimation of
#'   the levels of `variable` given the levels of `given`.
#' @param lvl_var,lvl_given (`character`)\cr corresponding levels given by [levels()].
#' @param mmat (named `numeric`) a vector filled with `0`s used as a template to obtain the design matrix.
#' @param coef (`numeric`)\cr vector of estimated coefficients.
#' @param vcov (`matrix`)\cr variance-covariance matrix of underlying model.
#' @param conf_level (`proportion`)\cr confidence level of estimate intervals.
#'
#' @details Given the cox regression investigating the effect of Arm (A, B, C; reference A)
#'   and Sex (F, M; reference Female). The model is abbreviated: y ~ Arm + Sex + Arm x Sex.
#'   The cox regression estimates the coefficients along with a variance-covariance matrix for:
#'
#'   - b1 (arm b), b2 (arm c)
#'   - b3 (sex m)
#'   - b4 (arm b: sex m), b5 (arm c: sex m)
#'
#'   Given that I want an estimation of the Hazard Ratio for arm C/sex M, the estimation
#'   will be given in reference to arm A/Sex M by exp(b2 + b3 + b5)/ exp(b3) = exp(b2 + b5),
#'   therefore the interaction coefficient is given by b2 + b5 while the standard error is obtained
#'   as $1.96 * sqrt(Var b2 + Var b5 + 2 * covariance (b2,b5))$ for a confidence level of 0.95.
#'
#' @return A list of matrices (one per level of variable) with rows corresponding to the combinations of
#'   `variable` and `given`, with columns:
#'   * `coef_hat`: Estimation of the coefficient.
#'   * `coef_se`: Standard error of the estimation.
#'   * `hr`: Hazard ratio.
#'   * `lcl, ucl`: Lower/upper confidence limit of the hazard ratio.
#'
#' @seealso [s_cox_multivariate()].
#'
#' @examples
#' library(dplyr)
#' library(survival)
#'
#' ADSL <- tern_ex_adsl %>%
#'   filter(SEX %in% c("F", "M"))
#'
#' adtte <- tern_ex_adtte %>% filter(PARAMCD == "PFS")
#' adtte$ARMCD <- droplevels(adtte$ARMCD)
#' adtte$SEX <- droplevels(adtte$SEX)
#'
#' mod <- coxph(
#'   formula = Surv(time = AVAL, event = 1 - CNSR) ~ (SEX + ARMCD)^2,
#'   data = adtte
#' )
#'
#' mmat <- stats::model.matrix(mod)[1, ]
#' mmat[!mmat == 0] <- 0
#'
#' @keywords internal
estimate_coef <- function(variable, given,
                          lvl_var, lvl_given,
                          coef,
                          mmat,
                          vcov,
                          conf_level = 0.95) {
  var_lvl <- paste0(variable, lvl_var[-1]) # [-1]: reference level
  giv_lvl <- paste0(given, lvl_given)

  design_mat <- expand.grid(variable = var_lvl, given = giv_lvl)
  design_mat <- design_mat[order(design_mat$variable, design_mat$given), ]
  design_mat <- within(
    data = design_mat,
    expr = {
      inter <- paste0(variable, ":", given)
      rev_inter <- paste0(given, ":", variable)
    }
  )

  split_by_variable <- design_mat$variable
  interaction_names <- paste(design_mat$variable, design_mat$given, sep = "/")

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

  betas <- as.matrix(coef)

  coef_hat <- t(design_mat) %*% betas
  dimnames(coef_hat)[2] <- "coef"

  coef_se <- apply(design_mat, 2, function(x) {
    vcov_el <- as.logical(x)
    y <- vcov[vcov_el, vcov_el]
    y <- sum(y)
    y <- sqrt(y)
    return(y)
  })

  q_norm <- stats::qnorm((1 + conf_level) / 2)
  y <- cbind(coef_hat, `se(coef)` = coef_se)

  y <- apply(y, 1, function(x) {
    x["hr"] <- exp(x["coef"])
    x["lcl"] <- exp(x["coef"] - q_norm * x["se(coef)"])
    x["ucl"] <- exp(x["coef"] + q_norm * x["se(coef)"])

    return(x)
  })

  y <- t(y)
  y <- by(y, split_by_variable, identity)
  y <- lapply(y, as.matrix)

  attr(y, "details") <- paste0(
    "Estimations of ", variable,
    " hazard ratio given the level of ", given, " compared to ",
    variable, " level ", lvl_var[1], "."
  )
  return(y)
}

#' `tryCatch` around `car::Anova`
#'
#' Captures warnings when executing [car::Anova].
#'
#' @inheritParams car::Anova
#'
#' @return A list with item `aov` for the result of the model and `error_text` for the captured warnings.
#'
#' @examples
#' # `car::Anova` on cox regression model including strata and expected
#' # a likelihood ratio test triggers a warning as only Wald method is
#' # accepted.
#'
#' library(survival)
#'
#' mod <- coxph(
#'   formula = Surv(time = futime, event = fustat) ~ factor(rx) + strata(ecog.ps),
#'   data = ovarian
#' )
#'
#' @keywords internal
try_car_anova <- function(mod,
                          test.statistic) { # nolint
  y <- tryCatch(
    withCallingHandlers(
      expr = {
        warn_text <- c()
        list(
          aov = car::Anova(
            mod,
            test.statistic = test.statistic,
            type = "III"
          ),
          warn_text = warn_text
        )
      },
      warning = function(w) {
        # If a warning is detected it is handled as "w".
        warn_text <<- trimws(paste0("Warning in `try_car_anova`: ", w))

        # A warning is sometimes expected, then, we want to restart
        # the execution while ignoring the warning.
        invokeRestart("muffleWarning")
      }
    ),
    finally = {
    }
  )

  return(y)
}

#' Fit a Cox regression model and ANOVA
#'
#' The functions derives the effect p-values using [car::Anova()] from [survival::coxph()] results.
#'
#' @inheritParams t_coxreg
#'
#' @return A list with items `mod` (results of [survival::coxph()]), `msum` (result of `summary`) and
#'   `aov` (result of [car::Anova()]).
#'
#' @noRd
fit_n_aov <- function(formula,
                      data = data,
                      conf_level = conf_level,
                      pval_method = c("wald", "likelihood"),
                      ...) {
  pval_method <- match.arg(pval_method)

  environment(formula) <- environment()
  suppressWarnings({
    # We expect some warnings due to coxph which fails strict programming.
    mod <- survival::coxph(formula, data = data, ...)
    msum <- summary(mod, conf.int = conf_level)
  })

  aov <- try_car_anova(
    mod,
    test.statistic = switch(pval_method,
      "wald" = "Wald",
      "likelihood" = "LR"
    )
  )

  warn_attr <- aov$warn_text
  if (!is.null(aov$warn_text)) message(warn_attr)

  aov <- aov$aov
  y <- list(mod = mod, msum = msum, aov = aov)
  attr(y, "message") <- warn_attr

  return(y)
}

# argument_checks
check_formula <- function(formula) {
  if (!(inherits(formula, "formula"))) {
    stop("Check `formula`. A formula should resemble `Surv(time = AVAL, event = 1 - CNSR) ~ study_arm(ARMCD)`.")
  }

  invisible()
}

check_covariate_formulas <- function(covariates) {
  if (!all(vapply(X = covariates, FUN = inherits, what = "formula", FUN.VALUE = TRUE)) || is.null(covariates)) {
    stop("Check `covariates`, it should be a list of right-hand-term formulas, e.g. list(Age = ~AGE).")
  }

  invisible()
}

name_covariate_names <- function(covariates) {
  miss_names <- names(covariates) == ""
  no_names <- is.null(names(covariates))
  if (any(miss_names)) names(covariates)[miss_names] <- vapply(covariates[miss_names], FUN = rht, FUN.VALUE = "name")
  if (no_names) names(covariates) <- vapply(covariates, FUN = rht, FUN.VALUE = "name")
  return(covariates)
}

check_increments <- function(increments, covariates) {
  if (!is.null(increments)) {
    covariates <- vapply(covariates, FUN = rht, FUN.VALUE = "name")
    lapply(
      X = names(increments), FUN = function(x) {
        if (!x %in% covariates) {
          warning(
            paste(
              "Check `increments`, the `increment` for ", x,
              "doesn't match any names in investigated covariate(s)."
            )
          )
        }
      }
    )
  }

  invisible()
}

#' Multivariate Cox model - summarized results
#'
#' Analyses based on multivariate Cox model are usually not performed for the Controlled Substance Reporting or
#' regulatory documents but serve exploratory purposes only (e.g., for publication). In practice, the model usually
#' includes only the main effects (without interaction terms). It produces the hazard ratio estimates for each of the
#' covariates included in the model.
#' The analysis follows the same principles (e.g., stratified vs. unstratified analysis and tie handling) as the
#' usual Cox model analysis. Since there is usually no pre-specified hypothesis testing for such analysis,
#' the p.values need to be interpreted with caution. (**Statistical Analysis of Clinical Trials Data with R**,
#' `NEST's bookdown`)
#'
#' @param formula (`formula`)\cr a formula corresponding to the investigated [survival::Surv()] survival model
#'   including covariates.
#' @param data (`data.frame`)\cr a data frame which includes the variable in formula and covariates.
#' @param conf_level (`proportion`)\cr the confidence level for the hazard ratio interval estimations. Default is 0.95.
#' @param pval_method (`string`)\cr the method used for the estimation of p-values, should be one of
#'   `"wald"` (default) or `"likelihood"`.
#' @param ... optional parameters passed to [survival::coxph()]. Can include `ties`, a character string specifying the
#'   method for tie handling, one of `exact` (default), `efron`, `breslow`.
#'
#' @return A `list` with elements `mod`, `msum`, `aov`, and `coef_inter`.
#'
#' @details The output is limited to single effect terms. Work in ongoing for estimation of interaction terms
#'   but is out of scope as defined by the  Global Data Standards Repository
#'   (**`GDS_Standard_TLG_Specs_Tables_2.doc`**).
#'
#' @seealso [estimate_coef()].
#'
#' @examples
#' library(dplyr)
#'
#' adtte <- tern_ex_adtte
#' adtte_f <- subset(adtte, PARAMCD == "OS") # _f: filtered
#' adtte_f <- filter(
#'   adtte_f,
#'   PARAMCD == "OS" &
#'     SEX %in% c("F", "M") &
#'     RACE %in% c("ASIAN", "BLACK OR AFRICAN AMERICAN", "WHITE")
#' )
#' adtte_f$SEX <- droplevels(adtte_f$SEX)
#' adtte_f$RACE <- droplevels(adtte_f$RACE)
#'
#' @keywords internal
s_cox_multivariate <- function(formula, data,
                               conf_level = 0.95,
                               pval_method = c("wald", "likelihood"),
                               ...) {
  tf <- stats::terms(formula, specials = c("strata"))
  covariates <- rownames(attr(tf, "factors"))[-c(1, unlist(attr(tf, "specials")))]
  lapply(
    X = covariates,
    FUN = function(x) {
      if (is.character(data[[x]])) {
        data[[x]] <<- as.factor(data[[x]])
      }
      invisible()
    }
  )
  pval_method <- match.arg(pval_method)

  # Results directly exported from environment(fit_n_aov) to environment(s_function_draft)
  y <- fit_n_aov(
    formula = formula,
    data = data,
    conf_level = conf_level,
    pval_method = pval_method,
    ...
  )
  mod <- y$mod
  aov <- y$aov
  msum <- y$msum
  list2env(as.list(y), environment())

  all_term_labs <- attr(mod$terms, "term.labels")
  term_labs <- all_term_labs[which(attr(mod$terms, "order") == 1)]
  names(term_labs) <- term_labs

  coef_inter <- NULL
  if (any(attr(mod$terms, "order") > 1)) {
    for_inter <- all_term_labs[attr(mod$terms, "order") > 1]
    names(for_inter) <- for_inter
    mmat <- stats::model.matrix(mod)[1, ]
    mmat[!mmat == 0] <- 0
    mcoef <- stats::coef(mod)
    mvcov <- stats::vcov(mod)

    estimate_coef_local <- function(variable, given) {
      estimate_coef(
        variable, given,
        coef = mcoef, mmat = mmat, vcov = mvcov, conf_level = conf_level,
        lvl_var = levels(data[[variable]]), lvl_given = levels(data[[given]])
      )
    }

    coef_inter <- lapply(
      for_inter, function(x) {
        y <- attr(mod$terms, "factors")[, x]
        y <- names(y[y > 0])
        Map(estimate_coef_local, variable = y, given = rev(y))
      }
    )
  }

  list(mod = mod, msum = msum, aov = aov, coef_inter = coef_inter)
}

#' Count patients with marked laboratory abnormalities
#'
#' @description `r lifecycle::badge("stable")`
#'
#' The analyze function [count_abnormal_by_marked()] creates a layout element to count patients with marked laboratory
#' abnormalities for each direction of abnormality, categorized by parameter value.
#'
#' This function analyzes primary analysis variable `var` which indicates whether a single, replicated,
#' or last marked laboratory abnormality was observed. Levels of `var` to include for each marked lab
#' abnormality (`single` and `last_replicated`) can be supplied via the `category` parameter. Additional
#' analysis variables that can be supplied as a list via the `variables` parameter are `id` (defaults
#' to `USUBJID`), a variable to indicate unique subject identifiers, `param` (defaults to `PARAM`), a
#' variable to indicate parameter values, and `direction` (defaults to `abn_dir`), a variable to indicate
#' abnormality directions.
#'
#' For each combination of `param` and `direction` levels, marked lab abnormality counts are calculated
#' as follows:
#'   * `Single, not last` & `Last or replicated`: The number of patients with `Single, not last`
#'     and `Last or replicated` values, respectively.
#'   * `Any`: The number of patients with either single or replicated marked abnormalities.
#'
#' Fractions are calculated by dividing the above counts by the number of patients with at least one
#' valid measurement recorded during the analysis.
#'
#' Prior to using this function in your table layout you must use [rtables::split_rows_by()] to create two
#' row splits, one on variable `param` and one on variable `direction`.
#'
#' @inheritParams argument_convention
#' @param category (`list`)\cr a list with different marked category names for single
#'   and last or replicated.
#' @param .stats (`character`)\cr statistics to select for the table. Run `get_stats("abnormal_by_marked")`
#'   to see available statistics for this function.
#'
#' @note `Single, not last` and `Last or replicated` levels are mutually exclusive. If a patient has
#'   abnormalities that meet both the `Single, not last` and `Last or replicated` criteria, then the
#'   patient will be counted only under the `Last or replicated` category.
#'
#' @name abnormal_by_marked
#' @order 1
NULL

#' @describeIn abnormal_by_marked Statistics function for patients with marked lab abnormalities.
#'
#' @return
#' * `s_count_abnormal_by_marked()` returns statistic `count_fraction` with `Single, not last`,
#'   `Last or replicated`, and `Any` results.
#'
#' @keywords internal
s_count_abnormal_by_marked <- function(df,
                                       .var = "AVALCAT1",
                                       .spl_context,
                                       category = list(single = "SINGLE", last_replicated = c("LAST", "REPLICATED")),
                                       variables = list(id = "USUBJID", param = "PARAM", direction = "abn_dir")) {
  checkmate::assert_string(.var)
  checkmate::assert_list(variables)
  checkmate::assert_list(category)
  checkmate::assert_subset(names(category), c("single", "last_replicated"))
  checkmate::assert_subset(names(variables), c("id", "param", "direction"))
  checkmate::assert_vector(unique(df[[variables$direction]]), max.len = 1)

  assert_df_with_variables(df, c(aval = .var, variables))
  checkmate::assert_multi_class(df[[.var]], classes = c("factor", "character"))
  checkmate::assert_multi_class(df[[variables$id]], classes = c("factor", "character"))


  first_row <- .spl_context[.spl_context$split == variables[["param"]], ]
  # Patients in the denominator have at least one post-baseline visit.
  subj <- first_row$full_parent_df[[1]][[variables[["id"]]]]
  subj_cur_col <- subj[first_row$cur_col_subset[[1]]]
  # Some subjects may have a record for high and low directions but
  # should be counted only once.
  denom <- length(unique(subj_cur_col))

  if (denom != 0) {
    subjects_last_replicated <- unique(
      df[df[[.var]] %in% category[["last_replicated"]], variables$id, drop = TRUE]
    )
    subjects_single <- unique(
      df[df[[.var]] %in% category[["single"]], variables$id, drop = TRUE]
    )
    # Subjects who have both single and last/replicated abnormalities are counted in only the last/replicated group.
    subjects_single <- setdiff(subjects_single, subjects_last_replicated)
    n_single <- length(subjects_single)
    n_last_replicated <- length(subjects_last_replicated)
    n_any <- n_single + n_last_replicated
    result <- list(count_fraction = list(
      "Single, not last" = c(n_single, n_single / denom),
      "Last or replicated" = c(n_last_replicated, n_last_replicated / denom),
      "Any Abnormality" = c(n_any, n_any / denom)
    ))
  } else {
    result <- list(count_fraction = list(
      "Single, not last" = c(0, 0),
      "Last or replicated" = c(0, 0),
      "Any Abnormality" = c(0, 0)
    ))
  }

  result
}

#' @describeIn abnormal_by_marked Formatted analysis function which is used as `afun`
#'   in `count_abnormal_by_marked()`.
#'
#' @return
#' * `a_count_abnormal_by_marked()` returns the corresponding list with formatted [rtables::CellValue()].
#'
#' @keywords internal
a_count_abnormal_by_marked <- make_afun(
  s_count_abnormal_by_marked,
  .formats = c(count_fraction = format_count_fraction)
)

#' @describeIn abnormal_by_marked Layout-creating function which can take statistics function arguments
#'   and additional format arguments. This function is a wrapper for [rtables::analyze()].
#'
#' @return
#' * `count_abnormal_by_marked()` 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_count_abnormal_by_marked()` to the table layout.
#'
#' @examples
#' library(dplyr)
#'
#' df <- data.frame(
#'   USUBJID = as.character(c(rep(1, 5), rep(2, 5), rep(1, 5), rep(2, 5))),
#'   ARMCD = factor(c(rep("ARM A", 5), rep("ARM B", 5), rep("ARM A", 5), rep("ARM B", 5))),
#'   ANRIND = factor(c(
#'     "NORMAL", "HIGH", "HIGH", "HIGH HIGH", "HIGH",
#'     "HIGH", "HIGH", "HIGH HIGH", "NORMAL", "HIGH HIGH", "NORMAL", "LOW", "LOW", "LOW LOW", "LOW",
#'     "LOW", "LOW", "LOW LOW", "NORMAL", "LOW LOW"
#'   )),
#'   ONTRTFL = rep(c("", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y"), 2),
#'   PARAMCD = factor(c(rep("CRP", 10), rep("ALT", 10))),
#'   AVALCAT1 = factor(rep(c("", "", "", "SINGLE", "REPLICATED", "", "", "LAST", "", "SINGLE"), 2)),
#'   stringsAsFactors = FALSE
#' )
#'
#' df <- df %>%


1		#' Count patients with toxicity grades that have worsened from baseline by highest grade post-baseline
2		#'
3		#' @description `r lifecycle::badge("stable")`
4		#'
5		#' The analyze function [count_abnormal_lab_worsen_by_baseline()] creates a layout element to count patients with
6		#' analysis toxicity grades which have worsened from baseline, categorized by highest (worst) grade post-baseline.
7		#'
8		#' This function analyzes primary analysis variable `var` which indicates analysis toxicity grades. Additional
9		#' analysis variables that can be supplied as a list via the `variables` parameter are `id` (defaults to `USUBJID`),
10		#' a variable to indicate unique subject identifiers, `baseline_var` (defaults to `BTOXGR`), a variable to indicate
11		#' baseline toxicity grades, and `direction_var` (defaults to `GRADDIR`), a variable to indicate toxicity grade
12		#' directions of interest to include (e.g. `"H"` (high), `"L"` (low), or `"B"` (both)).
13		#'
14		#' For the direction(s) specified in `direction_var`, patient counts by worst grade for patients who have
15		#' worsened from baseline are calculated as follows:
16		#' * `1` to `4`: The number of patients who have worsened from their baseline grades with worst
17		#' grades 1-4, respectively.
18		#' * `Any`: The total number of patients who have worsened from their baseline grades.
19		#'
20		#' Fractions are calculated by dividing the above counts by the number of patients who's analysis toxicity grades
21		#' have worsened from baseline toxicity grades during treatment.
22		#'
23		#' Prior to using this function in your table layout you must use [rtables::split_rows_by()] to create a row
24		#' split on variable `direction_var`.
25		#'
26		#' @inheritParams argument_convention
27		#' @param variables (named `list` of `string`)\cr list of additional analysis variables including:
28		#' * `id` (`string`)\cr subject variable name.
29		#' * `baseline_var` (`string`)\cr name of the data column containing baseline toxicity variable.
30		#' * `direction_var` (`string`)\cr see `direction_var` for more details.
31		#' @param .stats (`character`)\cr statistics to select for the table. Run `get_stats("abnormal_by_worst_grade_worsen")`
32		#' to see all available statistics.
33		#'
34		#' @seealso Relevant helper functions [h_adlb_worsen()] and [h_worsen_counter()] which are used within
35		#' [s_count_abnormal_lab_worsen_by_baseline()] to process input data.
36		#'
37		#' @name abnormal_by_worst_grade_worsen
38		#' @order 1
39		NULL
40
41		#' Helper function to prepare ADLB with worst labs
42		#'
43		#' @description `r lifecycle::badge("stable")`
44		#'
45		#' Helper function to prepare a `df` for generate the patient count shift table.
46		#'
47		#' @param adlb (`data.frame`)\cr ADLB data frame.
48		#' @param worst_flag_low (named `vector`)\cr worst low post-baseline lab grade flag variable. See how this is
49		#' implemented in the following examples.
50		#' @param worst_flag_high (named `vector`)\cr worst high post-baseline lab grade flag variable. See how this is
51		#' implemented in the following examples.
52		#' @param direction_var (`string`)\cr name of the direction variable specifying the direction of the shift table of
53		#' interest. Only lab records flagged by `L`, `H` or `B` are included in the shift table.
54		#' * `L`: low direction only
55		#' * `H`: high direction only
56		#' * `B`: both low and high directions
57		#'
58		#' @return `h_adlb_worsen()` returns the `adlb` `data.frame` containing only the
59		#' worst labs specified according to `worst_flag_low` or `worst_flag_high` for the
60		#' direction specified according to `direction_var`. For instance, for a lab that is
61		#' needed for the low direction only, only records flagged by `worst_flag_low` are
62		#' selected. For a lab that is needed for both low and high directions, the worst
63		#' low records are selected for the low direction, and the worst high record are selected
64		#' for the high direction.
65		#'
66		#' @seealso [abnormal_by_worst_grade_worsen]
67		#'
68		#' @examples
69		#' library(dplyr)
70		#'
71		#' # The direction variable, GRADDR, is based on metadata
72		#' adlb <- tern_ex_adlb %>%
73		#' mutate(
74		#' GRADDR = case_when(
75		#' PARAMCD == "ALT" ~ "B",
76		#' PARAMCD == "CRP" ~ "L",
77		#' PARAMCD == "IGA" ~ "H"
78		#' )
79		#' ) %>%
80		#' filter(SAFFL == "Y" & ONTRTFL == "Y" & GRADDR != "")
81		#'
82		#' df <- h_adlb_worsen(
83		#' adlb,
84		#' worst_flag_low = c("WGRLOFL" = "Y"),
85		#' worst_flag_high = c("WGRHIFL" = "Y"),
86		#' direction_var = "GRADDR"
87		#' )
88		#'
89		#' @export
90		h_adlb_worsen <- function(adlb,
91		worst_flag_low = NULL,
92		worst_flag_high = NULL,
93		direction_var) {
94	5x	checkmate::assert_string(direction_var)
95	5x	checkmate::assert_subset(as.character(unique(adlb[[direction_var]])), c("B", "L", "H"))
96	5x	assert_df_with_variables(adlb, list("Col" = direction_var))
97
98	5x	if (any(unique(adlb[[direction_var]]) == "H")) {
99	4x	assert_df_with_variables(adlb, list("High" = names(worst_flag_high)))
100		}
101
102	5x	if (any(unique(adlb[[direction_var]]) == "L")) {
103	4x	assert_df_with_variables(adlb, list("Low" = names(worst_flag_low)))
104		}
105
106	5x	if (any(unique(adlb[[direction_var]]) == "B")) {
107	3x	assert_df_with_variables(
108	3x	adlb,
109	3x	list(
110	3x	"Low" = names(worst_flag_low),
111	3x	"High" = names(worst_flag_high)
112		)
113		)
114		}
115
116		# extract patients with worst post-baseline lab, either low or high or both
117	5x	worst_flag <- c(worst_flag_low, worst_flag_high)
118	5x	col_names <- names(worst_flag)
119	5x	filter_values <- worst_flag
120	5x	temp <- Map(
121	5x	function(x, y) which(adlb[[x]] == y),
122	5x	col_names,
123	5x	filter_values
124		)
125	5x	position_satisfy_filters <- Reduce(union, temp)
126
127		# select variables of interest
128	5x	adlb_f <- adlb[position_satisfy_filters, ]
129
130		# generate subsets for different directionality
131	5x	adlb_f_h <- adlb_f[which(adlb_f[[direction_var]] == "H"), ]
132	5x	adlb_f_l <- adlb_f[which(adlb_f[[direction_var]] == "L"), ]
133	5x	adlb_f_b <- adlb_f[which(adlb_f[[direction_var]] == "B"), ]
134
135		# for labs requiring both high and low, data is duplicated and will be stacked on top of each other
136	5x	adlb_f_b_h <- adlb_f_b
137	5x	adlb_f_b_l <- adlb_f_b
138
139		# extract data with worst lab
140	5x	if (!is.null(worst_flag_high) && !is.null(worst_flag_low)) {
141		# change H to High, L to Low
142	3x	adlb_f_h[[direction_var]] <- rep("High", nrow(adlb_f_h))
143	3x	adlb_f_l[[direction_var]] <- rep("Low", nrow(adlb_f_l))
144
145		# change, B to High and Low
146	3x	adlb_f_b_h[[direction_var]] <- rep("High", nrow(adlb_f_b_h))
147	3x	adlb_f_b_l[[direction_var]] <- rep("Low", nrow(adlb_f_b_l))
148
149	3x	adlb_out_h <- adlb_f_h[which(adlb_f_h[[names(worst_flag_high)]] == worst_flag_high), ]
150	3x	adlb_out_b_h <- adlb_f_b_h[which(adlb_f_b_h[[names(worst_flag_high)]] == worst_flag_high), ]
151	3x	adlb_out_l <- adlb_f_l[which(adlb_f_l[[names(worst_flag_low)]] == worst_flag_low), ]
152	3x	adlb_out_b_l <- adlb_f_b_l[which(adlb_f_b_l[[names(worst_flag_low)]] == worst_flag_low), ]
153
154	3x	out <- rbind(adlb_out_h, adlb_out_b_h, adlb_out_l, adlb_out_b_l)
155	2x	} else if (!is.null(worst_flag_high)) {
156	1x	adlb_f_h[[direction_var]] <- rep("High", nrow(adlb_f_h))
157	1x	adlb_f_b_h[[direction_var]] <- rep("High", nrow(adlb_f_b_h))
158
159	1x	adlb_out_h <- adlb_f_h[which(adlb_f_h[[names(worst_flag_high)]] == worst_flag_high), ]
160	1x	adlb_out_b_h <- adlb_f_b_h[which(adlb_f_b_h[[names(worst_flag_high)]] == worst_flag_high), ]
161
162	1x	out <- rbind(adlb_out_h, adlb_out_b_h)
163	1x	} else if (!is.null(worst_flag_low)) {
164	1x	adlb_f_l[[direction_var]] <- rep("Low", nrow(adlb_f_l))
165	1x	adlb_f_b_l[[direction_var]] <- rep("Low", nrow(adlb_f_b_l))
166
167	1x	adlb_out_l <- adlb_f_l[which(adlb_f_l[[names(worst_flag_low)]] == worst_flag_low), ]
168	1x	adlb_out_b_l <- adlb_f_b_l[which(adlb_f_b_l[[names(worst_flag_low)]] == worst_flag_low), ]
169
170	1x	out <- rbind(adlb_out_l, adlb_out_b_l)
171		}
172
173		# label
174	5x	formatters::var_labels(out) <- formatters::var_labels(adlb_f, fill = FALSE)
175		# NA
176	5x	out
177		}
178
179		#' Helper function to analyze patients for `s_count_abnormal_lab_worsen_by_baseline()`
180		#'
181		#' @description `r lifecycle::badge("stable")`
182		#'
183		#' Helper function to count the number of patients and the fraction of patients according to
184		#' highest post-baseline lab grade variable `.var`, baseline lab grade variable `baseline_var`,
185		#' and the direction of interest specified in `direction_var`.
186		#'
187		#' @inheritParams argument_convention
188		#' @inheritParams h_adlb_worsen
189		#' @param baseline_var (`string`)\cr name of the baseline lab grade variable.
190		#'
191		#' @return The counts and fraction of patients
192		#' whose worst post-baseline lab grades are worse than their baseline grades, for
193		#' post-baseline worst grades "1", "2", "3", "4" and "Any".
194		#'
195		#' @seealso [abnormal_by_worst_grade_worsen]
196		#'
197		#' @examples
198		#' library(dplyr)
199		#'
200		#' # The direction variable, GRADDR, is based on metadata
201		#' adlb <- tern_ex_adlb %>%
202		#' mutate(
203		#' GRADDR = case_when(
204		#' PARAMCD == "ALT" ~ "B",
205		#' PARAMCD == "CRP" ~ "L",
206		#' PARAMCD == "IGA" ~ "H"
207		#' )
208		#' ) %>%
209		#' filter(SAFFL == "Y" & ONTRTFL == "Y" & GRADDR != "")
210		#'
211		#' df <- h_adlb_worsen(
212		#' adlb,
213		#' worst_flag_low = c("WGRLOFL" = "Y"),
214		#' worst_flag_high = c("WGRHIFL" = "Y"),
215		#' direction_var = "GRADDR"
216		#' )
217		#'
218		#' # `h_worsen_counter`
219		#' h_worsen_counter(
220		#' df %>% filter(PARAMCD == "CRP" & GRADDR == "Low"),
221		#' id = "USUBJID",
222		#' .var = "ATOXGR",
223		#' baseline_var = "BTOXGR",
224		#' direction_var = "GRADDR"
225		#' )
226		#'
227		#' @export
228		h_worsen_counter <- function(df, id, .var, baseline_var, direction_var) {
229	17x	checkmate::assert_string(id)
230	17x	checkmate::assert_string(.var)
231	17x	checkmate::assert_string(baseline_var)
232	17x	checkmate::assert_scalar(unique(df[[direction_var]]))
233	17x	checkmate::assert_subset(unique(df[[direction_var]]), c("High", "Low"))
234	17x	assert_df_with_variables(df, list(val = c(id, .var, baseline_var, direction_var)))
235
236		# remove post-baseline missing
237	17x	df <- df[df[[.var]] != "<Missing>", ]
238
239		# obtain directionality
240	17x	direction <- unique(df[[direction_var]])
241
242	17x	if (direction == "Low") {
243	10x	grade <- -1:-4
244	10x	worst_grade <- -4
245	7x	} else if (direction == "High") {
246	7x	grade <- 1:4
247	7x	worst_grade <- 4
248		}
249
250	17x	if (nrow(df) > 0) {
251	17x	by_grade <- lapply(grade, function(i) {
252		# filter baseline values that is less than i or <Missing>
253	68x	df_temp <- df[df[[baseline_var]] %in% c((i + sign(i) * -1):(-1 * worst_grade), "<Missing>"), ]
254		# num: number of patients with post-baseline worst lab equal to i
255	68x	num <- length(unique(df_temp[df_temp[[.var]] %in% i, id, drop = TRUE]))
256		# denom: number of patients with baseline values less than i or <missing> and post-baseline in the same direction
257	68x	denom <- length(unique(df_temp[[id]]))
258	68x	rm(df_temp)
259	68x	c(num = num, denom = denom)
260		})
261		} else {
262	!	by_grade <- lapply(1, function(i) {
263	!	c(num = 0, denom = 0)
264		})
265		}
266
267	17x	names(by_grade) <- as.character(seq_along(by_grade))
268
269		# baseline grade less 4 or missing
270	17x	df_temp <- df[!df[[baseline_var]] %in% worst_grade, ]
271
272		# denom: number of patients with baseline values less than 4 or <missing> and post-baseline in the same direction
273	17x	denom <- length(unique(df_temp[, id, drop = TRUE]))
274
275		# condition 1: missing baseline and in the direction of abnormality
276	17x	con1 <- which(df_temp[[baseline_var]] == "<Missing>" & df_temp[[.var]] %in% grade)
277	17x	df_temp_nm <- df_temp[which(df_temp[[baseline_var]] != "<Missing>" & df_temp[[.var]] %in% grade), ]
278
279		# condition 2: if post-baseline values are present then post-baseline values must be worse than baseline
280	17x	if (direction == "Low") {
281	10x	con2 <- which(as.numeric(as.character(df_temp_nm[[.var]])) < as.numeric(as.character(df_temp_nm[[baseline_var]])))
282		} else {
283	7x	con2 <- which(as.numeric(as.character(df_temp_nm[[.var]])) > as.numeric(as.character(df_temp_nm[[baseline_var]])))
284		}
285
286		# number of patients satisfy either conditions 1 or 2
287	17x	num <- length(unique(df_temp[union(con1, con2), id, drop = TRUE]))
288
289	17x	list(fraction = c(by_grade, list("Any" = c(num = num, denom = denom))))
290		}
291
292		#' @describeIn abnormal_by_worst_grade_worsen Statistics function for patients whose worst post-baseline
293		#' lab grades are worse than their baseline grades.
294		#'
295		#' @return
296		#' * `s_count_abnormal_lab_worsen_by_baseline()` returns the counts and fraction of patients whose worst
297		#' post-baseline lab grades are worse than their baseline grades, for post-baseline worst grades
298		#' "1", "2", "3", "4" and "Any".
299		#'
300		#' @keywords internal
301		s_count_abnormal_lab_worsen_by_baseline <- function(df, # nolint
302		.var = "ATOXGR",
303		variables = list(
304		id = "USUBJID",
305		baseline_var = "BTOXGR",
306		direction_var = "GRADDR"
307		)) {
308	1x	checkmate::assert_string(.var)
309	1x	checkmate::assert_set_equal(names(variables), c("id", "baseline_var", "direction_var"))
310	1x	checkmate::assert_string(variables$id)
311	1x	checkmate::assert_string(variables$baseline_var)
312	1x	checkmate::assert_string(variables$direction_var)
313	1x	assert_df_with_variables(df, c(aval = .var, variables[1:3]))
314	1x	assert_list_of_variables(variables)
315
316	1x	h_worsen_counter(df, variables$id, .var, variables$baseline_var, variables$direction_var)
317		}
318
319		#' @describeIn abnormal_by_worst_grade_worsen Formatted analysis function which is used as `afun`
320		#' in `count_abnormal_lab_worsen_by_baseline()`.
321		#'
322		#' @return
323		#' * `a_count_abnormal_lab_worsen_by_baseline()` returns the corresponding list with
324		#' formatted [rtables::CellValue()].
325		#'
326		#' @keywords internal
327		a_count_abnormal_lab_worsen_by_baseline <- make_afun( # nolint
328		s_count_abnormal_lab_worsen_by_baseline,
329		.formats = c(fraction = format_fraction),
330		.ungroup_stats = "fraction"
331		)
332
333		#' @describeIn abnormal_by_worst_grade_worsen Layout-creating function which can take statistics function
334		#' arguments and additional format arguments. This function is a wrapper for [rtables::analyze()].
335		#'
336		#' @return
337		#' * `count_abnormal_lab_worsen_by_baseline()` returns a layout object suitable for passing to further layouting
338		#' functions, or to [rtables::build_table()]. Adding this function to an `rtable` layout will add formatted
339		#' rows containing the statistics from `s_count_abnormal_lab_worsen_by_baseline()` to the table layout.
340		#'
341		#' @examples
342		#' library(dplyr)
343		#'
344		#' # The direction variable, GRADDR, is based on metadata
345		#' adlb <- tern_ex_adlb %>%
346		#' mutate(
347		#' GRADDR = case_when(
348		#' PARAMCD == "ALT" ~ "B",
349		#' PARAMCD == "CRP" ~ "L",
350		#' PARAMCD == "IGA" ~ "H"
351		#' )
352		#' ) %>%
353		#' filter(SAFFL == "Y" & ONTRTFL == "Y" & GRADDR != "")
354		#'
355		#' df <- h_adlb_worsen(
356		#' adlb,
357		#' worst_flag_low = c("WGRLOFL" = "Y"),
358		#' worst_flag_high = c("WGRHIFL" = "Y"),
359		#' direction_var = "GRADDR"
360		#' )
361		#'
362		#' basic_table() %>%
363		#' split_cols_by("ARMCD") %>%
364		#' add_colcounts() %>%
365		#' split_rows_by("PARAMCD") %>%
366		#' split_rows_by("GRADDR") %>%
367		#' count_abnormal_lab_worsen_by_baseline(
368		#' var = "ATOXGR",
369		#' variables = list(
370		#' id = "USUBJID",
371		#' baseline_var = "BTOXGR",
372		#' direction_var = "GRADDR"
373		#' )
374		#' ) %>%
375		#' append_topleft("Direction of Abnormality") %>%
376		#' build_table(df = df, alt_counts_df = tern_ex_adsl)
377		#'
378		#' @export
379		#' @order 2
380		count_abnormal_lab_worsen_by_baseline <- function(lyt, # nolint
381		var,
382		variables = list(
383		id = "USUBJID",
384		baseline_var = "BTOXGR",
385		direction_var = "GRADDR"
386		),
387		na_str = default_na_str(),
388		nested = TRUE,
389		...,
390		table_names = NULL,
391		.stats = NULL,
392		.formats = NULL,
393		.labels = NULL,
394		.indent_mods = NULL) {
395	1x	checkmate::assert_string(var)
396
397	1x	extra_args <- list(variables = variables, ...)
398
399	1x	afun <- make_afun(
400	1x	a_count_abnormal_lab_worsen_by_baseline,
401	1x	.stats = .stats,
402	1x	.formats = .formats,
403	1x	.labels = .labels,
404	1x	.indent_mods = .indent_mods
405		)
406
407	1x	lyt <- analyze(
408	1x	lyt = lyt,
409	1x	vars = var,
410	1x	afun = afun,
411	1x	na_str = na_str,
412	1x	nested = nested,
413	1x	extra_args = extra_args,
414	1x	show_labels = "hidden"
415		)
416
417	1x	lyt
418		}

1		#' Univariate formula special term
2		#'
3		#' @description `r lifecycle::badge("stable")`
4		#'
5		#' The special term `univariate` indicate that the model should be fitted individually for
6		#' every variable included in univariate.
7		#'
8		#' @param x (`character`)\cr a vector of variable names separated by commas.
9		#'
10		#' @return When used within a model formula, produces univariate models for each variable provided.
11		#'
12		#' @details
13		#' If provided alongside with pairwise specification, the model
14		#' `y ~ ARM + univariate(SEX, AGE, RACE)` lead to the study and comparison of the models
15		#' + `y ~ ARM`
16		#' + `y ~ ARM + SEX`
17		#' + `y ~ ARM + AGE`
18		#' + `y ~ ARM + RACE`
19		#'
20		#' @export
21		univariate <- function(x) {
22	1x	structure(x, varname = deparse(substitute(x)))
23		}
24
25		# Get the right-hand-term of a formula
26		rht <- function(x) {
27	4x	checkmate::assert_formula(x)
28	4x	y <- as.character(rev(x)[[1]])
29	4x	return(y)
30		}
31
32		#' Hazard ratio estimation in interactions
33		#'
34		#' This function estimates the hazard ratios between arms when an interaction variable is given with
35		#' specific values.
36		#'
37		#' @param variable,given (`character(2)`)\cr names of the two variables in the interaction. We seek the estimation of
38		#' the levels of `variable` given the levels of `given`.
39		#' @param lvl_var,lvl_given (`character`)\cr corresponding levels given by [levels()].
40		#' @param mmat (named `numeric`) a vector filled with `0`s used as a template to obtain the design matrix.
41		#' @param coef (`numeric`)\cr vector of estimated coefficients.
42		#' @param vcov (`matrix`)\cr variance-covariance matrix of underlying model.
43		#' @param conf_level (`proportion`)\cr confidence level of estimate intervals.
44		#'
45		#' @details Given the cox regression investigating the effect of Arm (A, B, C; reference A)
46		#' and Sex (F, M; reference Female). The model is abbreviated: y ~ Arm + Sex + Arm x Sex.
47		#' The cox regression estimates the coefficients along with a variance-covariance matrix for:
48		#'
49		#' - b1 (arm b), b2 (arm c)
50		#' - b3 (sex m)
51		#' - b4 (arm b: sex m), b5 (arm c: sex m)
52		#'
53		#' Given that I want an estimation of the Hazard Ratio for arm C/sex M, the estimation
54		#' will be given in reference to arm A/Sex M by exp(b2 + b3 + b5)/ exp(b3) = exp(b2 + b5),
55		#' therefore the interaction coefficient is given by b2 + b5 while the standard error is obtained
56		#' as $1.96 * sqrt(Var b2 + Var b5 + 2 * covariance (b2,b5))$ for a confidence level of 0.95.
57		#'
58		#' @return A list of matrices (one per level of variable) with rows corresponding to the combinations of
59		#' `variable` and `given`, with columns:
60		#' * `coef_hat`: Estimation of the coefficient.
61		#' * `coef_se`: Standard error of the estimation.
62		#' * `hr`: Hazard ratio.
63		#' * `lcl, ucl`: Lower/upper confidence limit of the hazard ratio.
64		#'
65		#' @seealso [s_cox_multivariate()].
66		#'
67		#' @examples
68		#' library(dplyr)
69		#' library(survival)
70		#'
71		#' ADSL <- tern_ex_adsl %>%
72		#' filter(SEX %in% c("F", "M"))
73		#'
74		#' adtte <- tern_ex_adtte %>% filter(PARAMCD == "PFS")
75		#' adtte$ARMCD <- droplevels(adtte$ARMCD)
76		#' adtte$SEX <- droplevels(adtte$SEX)
77		#'
78		#' mod <- coxph(
79		#' formula = Surv(time = AVAL, event = 1 - CNSR) ~ (SEX + ARMCD)^2,
80		#' data = adtte
81		#' )
82		#'
83		#' mmat <- stats::model.matrix(mod)[1, ]
84		#' mmat[!mmat == 0] <- 0
85		#'
86		#' @keywords internal
87		estimate_coef <- function(variable, given,
88		lvl_var, lvl_given,
89		coef,
90		mmat,
91		vcov,
92		conf_level = 0.95) {
93	8x	var_lvl <- paste0(variable, lvl_var[-1]) # [-1]: reference level
94	8x	giv_lvl <- paste0(given, lvl_given)
95
96	8x	design_mat <- expand.grid(variable = var_lvl, given = giv_lvl)
97	8x	design_mat <- design_mat[order(design_mat$variable, design_mat$given), ]
98	8x	design_mat <- within(
99	8x	data = design_mat,
100	8x	expr = {
101	8x	inter <- paste0(variable, ":", given)
102	8x	rev_inter <- paste0(given, ":", variable)
103		}
104		)
105
106	8x	split_by_variable <- design_mat$variable
107	8x	interaction_names <- paste(design_mat$variable, design_mat$given, sep = "/")
108
109	8x	design_mat <- apply(
110	8x	X = design_mat, MARGIN = 1, FUN = function(x) {
111	27x	mmat[names(mmat) %in% x[-which(names(x) == "given")]] <- 1
112	27x	return(mmat)
113		}
114		)
115	8x	colnames(design_mat) <- interaction_names
116
117	8x	betas <- as.matrix(coef)
118
119	8x	coef_hat <- t(design_mat) %*% betas
120	8x	dimnames(coef_hat)[2] <- "coef"
121
122	8x	coef_se <- apply(design_mat, 2, function(x) {
123	27x	vcov_el <- as.logical(x)
124	27x	y <- vcov[vcov_el, vcov_el]
125	27x	y <- sum(y)
126	27x	y <- sqrt(y)
127	27x	return(y)
128		})
129
130	8x	q_norm <- stats::qnorm((1 + conf_level) / 2)
131	8x	y <- cbind(coef_hat, `se(coef)` = coef_se)
132
133	8x	y <- apply(y, 1, function(x) {
134	27x	x["hr"] <- exp(x["coef"])
135	27x	x["lcl"] <- exp(x["coef"] - q_norm * x["se(coef)"])
136	27x	x["ucl"] <- exp(x["coef"] + q_norm * x["se(coef)"])
137
138	27x	return(x)
139		})
140
141	8x	y <- t(y)
142	8x	y <- by(y, split_by_variable, identity)
143	8x	y <- lapply(y, as.matrix)
144
145	8x	attr(y, "details") <- paste0(
146	8x	"Estimations of ", variable,
147	8x	" hazard ratio given the level of ", given, " compared to ",
148	8x	variable, " level ", lvl_var[1], "."
149		)
150	8x	return(y)
151		}
152
153		#' `tryCatch` around `car::Anova`
154		#'
155		#' Captures warnings when executing [car::Anova].
156		#'
157		#' @inheritParams car::Anova
158		#'
159		#' @return A list with item `aov` for the result of the model and `error_text` for the captured warnings.
160		#'
161		#' @examples
162		#' # `car::Anova` on cox regression model including strata and expected
163		#' # a likelihood ratio test triggers a warning as only Wald method is
164		#' # accepted.
165		#'
166		#' library(survival)
167		#'
168		#' mod <- coxph(
169		#' formula = Surv(time = futime, event = fustat) ~ factor(rx) + strata(ecog.ps),
170		#' data = ovarian
171		#' )
172		#'
173		#' @keywords internal
174		try_car_anova <- function(mod,
175		test.statistic) { # nolint
176	2x	y <- tryCatch(
177	2x	withCallingHandlers(
178	2x	expr = {
179	2x	warn_text <- c()
180	2x	list(
181	2x	aov = car::Anova(
182	2x	mod,
183	2x	test.statistic = test.statistic,
184	2x	type = "III"
185		),
186	2x	warn_text = warn_text
187		)
188		},
189	2x	warning = function(w) {
190		# If a warning is detected it is handled as "w".
191	!	warn_text <<- trimws(paste0("Warning in `try_car_anova`: ", w))
192
193		# A warning is sometimes expected, then, we want to restart
194		# the execution while ignoring the warning.
195	!	invokeRestart("muffleWarning")
196		}
197		),
198	2x	finally = {
199		}
200		)
201
202	2x	return(y)
203		}
204
205		#' Fit a Cox regression model and ANOVA
206		#'
207		#' The functions derives the effect p-values using [car::Anova()] from [survival::coxph()] results.
208		#'
209		#' @inheritParams t_coxreg
210		#'
211		#' @return A list with items `mod` (results of [survival::coxph()]), `msum` (result of `summary`) and
212		#' `aov` (result of [car::Anova()]).
213		#'
214		#' @noRd
215		fit_n_aov <- function(formula,
216		data = data,
217		conf_level = conf_level,
218		pval_method = c("wald", "likelihood"),
219		...) {
220	1x	pval_method <- match.arg(pval_method)
221
222	1x	environment(formula) <- environment()
223	1x	suppressWarnings({
224		# We expect some warnings due to coxph which fails strict programming.
225	1x	mod <- survival::coxph(formula, data = data, ...)
226	1x	msum <- summary(mod, conf.int = conf_level)
227		})
228
229	1x	aov <- try_car_anova(
230	1x	mod,
231	1x	test.statistic = switch(pval_method,
232	1x	"wald" = "Wald",
233	1x	"likelihood" = "LR"
234		)
235		)
236
237	1x	warn_attr <- aov$warn_text
238	!	if (!is.null(aov$warn_text)) message(warn_attr)
239
240	1x	aov <- aov$aov
241	1x	y <- list(mod = mod, msum = msum, aov = aov)
242	1x	attr(y, "message") <- warn_attr
243
244	1x	return(y)
245		}
246
247		# argument_checks
248		check_formula <- function(formula) {
249	1x	if (!(inherits(formula, "formula"))) {
250	1x	stop("Check `formula`. A formula should resemble `Surv(time = AVAL, event = 1 - CNSR) ~ study_arm(ARMCD)`.")
251		}
252
253	!	invisible()
254		}
255
256		check_covariate_formulas <- function(covariates) {
257	1x	if (!all(vapply(X = covariates, FUN = inherits, what = "formula", FUN.VALUE = TRUE)) \|\| is.null(covariates)) {
258	1x	stop("Check `covariates`, it should be a list of right-hand-term formulas, e.g. list(Age = ~AGE).")
259		}
260
261	!	invisible()
262		}
263
264		name_covariate_names <- function(covariates) {
265	1x	miss_names <- names(covariates) == ""
266	1x	no_names <- is.null(names(covariates))
267	!	if (any(miss_names)) names(covariates)[miss_names] <- vapply(covariates[miss_names], FUN = rht, FUN.VALUE = "name")
268	!	if (no_names) names(covariates) <- vapply(covariates, FUN = rht, FUN.VALUE = "name")
269	1x	return(covariates)
270		}
271
272		check_increments <- function(increments, covariates) {
273	1x	if (!is.null(increments)) {
274	1x	covariates <- vapply(covariates, FUN = rht, FUN.VALUE = "name")
275	1x	lapply(
276	1x	X = names(increments), FUN = function(x) {
277	3x	if (!x %in% covariates) {
278	1x	warning(
279	1x	paste(
280	1x	"Check `increments`, the `increment` for ", x,
281	1x	"doesn't match any names in investigated covariate(s)."
282		)
283		)
284		}
285		}
286		)
287		}
288
289	1x	invisible()
290		}
291
292		#' Multivariate Cox model - summarized results
293		#'
294		#' Analyses based on multivariate Cox model are usually not performed for the Controlled Substance Reporting or
295		#' regulatory documents but serve exploratory purposes only (e.g., for publication). In practice, the model usually
296		#' includes only the main effects (without interaction terms). It produces the hazard ratio estimates for each of the
297		#' covariates included in the model.
298		#' The analysis follows the same principles (e.g., stratified vs. unstratified analysis and tie handling) as the
299		#' usual Cox model analysis. Since there is usually no pre-specified hypothesis testing for such analysis,
300		#' the p.values need to be interpreted with caution. (Statistical Analysis of Clinical Trials Data with R,
301		#' `NEST's bookdown`)
302		#'
303		#' @param formula (`formula`)\cr a formula corresponding to the investigated [survival::Surv()] survival model
304		#' including covariates.
305		#' @param data (`data.frame`)\cr a data frame which includes the variable in formula and covariates.
306		#' @param conf_level (`proportion`)\cr the confidence level for the hazard ratio interval estimations. Default is 0.95.
307		#' @param pval_method (`string`)\cr the method used for the estimation of p-values, should be one of
308		#' `"wald"` (default) or `"likelihood"`.
309		#' @param ... optional parameters passed to [survival::coxph()]. Can include `ties`, a character string specifying the
310		#' method for tie handling, one of `exact` (default), `efron`, `breslow`.
311		#'
312		#' @return A `list` with elements `mod`, `msum`, `aov`, and `coef_inter`.
313		#'
314		#' @details The output is limited to single effect terms. Work in ongoing for estimation of interaction terms
315		#' but is out of scope as defined by the Global Data Standards Repository
316		#' (`GDS_Standard_TLG_Specs_Tables_2.doc`).
317		#'
318		#' @seealso [estimate_coef()].
319		#'
320		#' @examples
321		#' library(dplyr)
322		#'
323		#' adtte <- tern_ex_adtte
324		#' adtte_f <- subset(adtte, PARAMCD == "OS") # _f: filtered
325		#' adtte_f <- filter(
326		#' adtte_f,
327		#' PARAMCD == "OS" &
328		#' SEX %in% c("F", "M") &
329		#' RACE %in% c("ASIAN", "BLACK OR AFRICAN AMERICAN", "WHITE")
330		#' )
331		#' adtte_f$SEX <- droplevels(adtte_f$SEX)
332		#' adtte_f$RACE <- droplevels(adtte_f$RACE)
333		#'
334		#' @keywords internal
335		s_cox_multivariate <- function(formula, data,
336		conf_level = 0.95,
337		pval_method = c("wald", "likelihood"),
338		...) {
339	1x	tf <- stats::terms(formula, specials = c("strata"))
340	1x	covariates <- rownames(attr(tf, "factors"))[-c(1, unlist(attr(tf, "specials")))]
341	1x	lapply(
342	1x	X = covariates,
343	1x	FUN = function(x) {
344	3x	if (is.character(data[[x]])) {
345	1x	data[[x]] <<- as.factor(data[[x]])
346		}
347	3x	invisible()
348		}
349		)
350	1x	pval_method <- match.arg(pval_method)
351
352		# Results directly exported from environment(fit_n_aov) to environment(s_function_draft)
353	1x	y <- fit_n_aov(
354	1x	formula = formula,
355	1x	data = data,
356	1x	conf_level = conf_level,
357	1x	pval_method = pval_method,
358		...
359		)
360	1x	mod <- y$mod
361	1x	aov <- y$aov
362	1x	msum <- y$msum
363	1x	list2env(as.list(y), environment())
364
365	1x	all_term_labs <- attr(mod$terms, "term.labels")
366	1x	term_labs <- all_term_labs[which(attr(mod$terms, "order") == 1)]
367	1x	names(term_labs) <- term_labs
368
369	1x	coef_inter <- NULL
370	1x	if (any(attr(mod$terms, "order") > 1)) {
371	1x	for_inter <- all_term_labs[attr(mod$terms, "order") > 1]
372	1x	names(for_inter) <- for_inter
373	1x	mmat <- stats::model.matrix(mod)[1, ]
374	1x	mmat[!mmat == 0] <- 0
375	1x	mcoef <- stats::coef(mod)
376	1x	mvcov <- stats::vcov(mod)
377
378	1x	estimate_coef_local <- function(variable, given) {
379	6x	estimate_coef(
380	6x	variable, given,
381	6x	coef = mcoef, mmat = mmat, vcov = mvcov, conf_level = conf_level,
382	6x	lvl_var = levels(data[[variable]]), lvl_given = levels(data[[given]])
383		)
384		}
385
386	1x	coef_inter <- lapply(
387	1x	for_inter, function(x) {
388	3x	y <- attr(mod$terms, "factors")[, x]
389	3x	y <- names(y[y > 0])
390	3x	Map(estimate_coef_local, variable = y, given = rev(y))
391		}
392		)
393		}
394
395	1x	list(mod = mod, msum = msum, aov = aov, coef_inter = coef_inter)
396		}

1		#' Count patients with marked laboratory abnormalities
2		#'
3		#' @description `r lifecycle::badge("stable")`
4		#'
5		#' The analyze function [count_abnormal_by_marked()] creates a layout element to count patients with marked laboratory
6		#' abnormalities for each direction of abnormality, categorized by parameter value.
7		#'
8		#' This function analyzes primary analysis variable `var` which indicates whether a single, replicated,
9		#' or last marked laboratory abnormality was observed. Levels of `var` to include for each marked lab
10		#' abnormality (`single` and `last_replicated`) can be supplied via the `category` parameter. Additional
11		#' analysis variables that can be supplied as a list via the `variables` parameter are `id` (defaults
12		#' to `USUBJID`), a variable to indicate unique subject identifiers, `param` (defaults to `PARAM`), a
13		#' variable to indicate parameter values, and `direction` (defaults to `abn_dir`), a variable to indicate
14		#' abnormality directions.
15		#'
16		#' For each combination of `param` and `direction` levels, marked lab abnormality counts are calculated
17		#' as follows:
18		#' * `Single, not last` & `Last or replicated`: The number of patients with `Single, not last`
19		#' and `Last or replicated` values, respectively.
20		#' * `Any`: The number of patients with either single or replicated marked abnormalities.
21		#'
22		#' Fractions are calculated by dividing the above counts by the number of patients with at least one
23		#' valid measurement recorded during the analysis.
24		#'
25		#' Prior to using this function in your table layout you must use [rtables::split_rows_by()] to create two
26		#' row splits, one on variable `param` and one on variable `direction`.
27		#'
28		#' @inheritParams argument_convention
29		#' @param category (`list`)\cr a list with different marked category names for single
30		#' and last or replicated.
31		#' @param .stats (`character`)\cr statistics to select for the table. Run `get_stats("abnormal_by_marked")`
32		#' to see available statistics for this function.
33		#'
34		#' @note `Single, not last` and `Last or replicated` levels are mutually exclusive. If a patient has
35		#' abnormalities that meet both the `Single, not last` and `Last or replicated` criteria, then the
36		#' patient will be counted only under the `Last or replicated` category.
37		#'
38		#' @name abnormal_by_marked
39		#' @order 1
40		NULL
41
42		#' @describeIn abnormal_by_marked Statistics function for patients with marked lab abnormalities.
43		#'
44		#' @return
45		#' * `s_count_abnormal_by_marked()` returns statistic `count_fraction` with `Single, not last`,
46		#' `Last or replicated`, and `Any` results.
47		#'
48		#' @keywords internal
49		s_count_abnormal_by_marked <- function(df,
50		.var = "AVALCAT1",
51		.spl_context,
52		category = list(single = "SINGLE", last_replicated = c("LAST", "REPLICATED")),
53		variables = list(id = "USUBJID", param = "PARAM", direction = "abn_dir")) {
54	3x	checkmate::assert_string(.var)
55	3x	checkmate::assert_list(variables)
56	3x	checkmate::assert_list(category)
57	3x	checkmate::assert_subset(names(category), c("single", "last_replicated"))
58	3x	checkmate::assert_subset(names(variables), c("id", "param", "direction"))
59	3x	checkmate::assert_vector(unique(df[[variables$direction]]), max.len = 1)
60
61	2x	assert_df_with_variables(df, c(aval = .var, variables))
62	2x	checkmate::assert_multi_class(df[[.var]], classes = c("factor", "character"))
63	2x	checkmate::assert_multi_class(df[[variables$id]], classes = c("factor", "character"))
64
65
66	2x	first_row <- .spl_context[.spl_context$split == variables[["param"]], ]
67		# Patients in the denominator have at least one post-baseline visit.
68	2x	subj <- first_row$full_parent_df[[1]][[variables[["id"]]]]
69	2x	subj_cur_col <- subj[first_row$cur_col_subset[[1]]]
70		# Some subjects may have a record for high and low directions but
71		# should be counted only once.
72	2x	denom <- length(unique(subj_cur_col))
73
74	2x	if (denom != 0) {
75	2x	subjects_last_replicated <- unique(
76	2x	df[df[[.var]] %in% category[["last_replicated"]], variables$id, drop = TRUE]
77		)
78	2x	subjects_single <- unique(
79	2x	df[df[[.var]] %in% category[["single"]], variables$id, drop = TRUE]
80		)
81		# Subjects who have both single and last/replicated abnormalities are counted in only the last/replicated group.
82	2x	subjects_single <- setdiff(subjects_single, subjects_last_replicated)
83	2x	n_single <- length(subjects_single)
84	2x	n_last_replicated <- length(subjects_last_replicated)
85	2x	n_any <- n_single + n_last_replicated
86	2x	result <- list(count_fraction = list(
87	2x	"Single, not last" = c(n_single, n_single / denom),
88	2x	"Last or replicated" = c(n_last_replicated, n_last_replicated / denom),
89	2x	"Any Abnormality" = c(n_any, n_any / denom)
90		))
91		} else {
92	!	result <- list(count_fraction = list(
93	!	"Single, not last" = c(0, 0),
94	!	"Last or replicated" = c(0, 0),
95	!	"Any Abnormality" = c(0, 0)
96		))
97		}
98
99	2x	result
100		}
101
102		#' @describeIn abnormal_by_marked Formatted analysis function which is used as `afun`
103		#' in `count_abnormal_by_marked()`.
104		#'
105		#' @return
106		#' * `a_count_abnormal_by_marked()` returns the corresponding list with formatted [rtables::CellValue()].
107		#'
108		#' @keywords internal
109		a_count_abnormal_by_marked <- make_afun(
110		s_count_abnormal_by_marked,
111		.formats = c(count_fraction = format_count_fraction)
112		)
113
114		#' @describeIn abnormal_by_marked Layout-creating function which can take statistics function arguments
115		#' and additional format arguments. This function is a wrapper for [rtables::analyze()].
116		#'
117		#' @return
118		#' * `count_abnormal_by_marked()` returns a layout object suitable for passing to further layouting functions,
119		#' or to [rtables::build_table()]. Adding this function to an `rtable` layout will add formatted rows containing
120		#' the statistics from `s_count_abnormal_by_marked()` to the table layout.
121		#'
122		#' @examples
123		#' library(dplyr)
124		#'
125		#' df <- data.frame(
126		#' USUBJID = as.character(c(rep(1, 5), rep(2, 5), rep(1, 5), rep(2, 5))),
127		#' ARMCD = factor(c(rep("ARM A", 5), rep("ARM B", 5), rep("ARM A", 5), rep("ARM B", 5))),
128		#' ANRIND = factor(c(
129		#' "NORMAL", "HIGH", "HIGH", "HIGH HIGH", "HIGH",
130		#' "HIGH", "HIGH", "HIGH HIGH", "NORMAL", "HIGH HIGH", "NORMAL", "LOW", "LOW", "LOW LOW", "LOW",
131		#' "LOW", "LOW", "LOW LOW", "NORMAL", "LOW LOW"
132		#' )),
133		#' ONTRTFL = rep(c("", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y", "Y"), 2),
134		#' PARAMCD = factor(c(rep("CRP", 10), rep("ALT", 10))),
135		#' AVALCAT1 = factor(rep(c("", "", "", "SINGLE", "REPLICATED", "", "", "LAST", "", "SINGLE"), 2)),
136		#' stringsAsFactors = FALSE
137		#' )
138		#'
139		#' df <- df %>%
140