R/plot.simBA.R

In sim.BA: Simulation-Based Bias Analysis for Observational Studies

#' Plot the results of a simulation
#'
#' `plot()` plots the output of a call to `simBA()`. The plot can contain either the estimated hazard ratios or standardized mean differences across simulations, each in a set of box plots.
#'
#' @param x a `simBA` object; the output of a call to [simBA()].
#' @param type the type of plot to produce; allowable options include `"balance"` (default), `"hr"`, and `"bias"`. Abbreviations allowed. See Details.
#' @param ... further arguments passed to [ggplot2::geom_boxplot()].
#'
#' @returns
#' A `ggplot` object, which can be modified using `ggplot2` syntax.
#'
#' @details
#' The balance plot plots absolute standardized mean differences. Vertical lines are placed at 0 (solid) and .1 (dashed). The hazard ratio (HR) plot plots hazard ratios on a log scale for the x-axis. Vertical lines are placed at 1 (solid) and the true marginal HR (dashed). The bias plot plots the relative error in the HR with a vertical line at 0% (indicating no error).
#'
#' @seealso [simBA()] for performing the simulation.
#'
#' @examples
#' # See help("simBA") for examples.
#'

#' @exportS3Method plot simBA
plot.simBA <- function(x, type = "balance", ...) {

  chk::chk_string(type)
  type <- tolower(type)
  type <- match_arg(type, c("balance", "hr", "bias"))
  adj <- attr(x, "adj")

  if (type == "balance") {
    bal <- do.call("rbind", lapply(x$sim_out, `[[`, "all_balance_tables"))

    proxy_vars <- sort(setdiff(unique(bal$variable), attr(x, "unmeasured_conf")))
    new_proxy_vars <- sub("p", "Proxy ", proxy_vars, fixed = TRUE)

    ## prepare the the SMD table to output
    bal$SMD[is.na(bal$SMD)] <- 0

    bal$variable <- factor(bal$variable,
                            levels = c(attr(x, "unmeasured_conf"), rev(proxy_vars)),
                            labels = c("Unmeasured\nConfounder", rev(new_proxy_vars)))
    bal$adjustment <- factor(bal$adjustment, levels = c("crude", "L1", "L2"),
                             labels = c("Unadjusted",
                                        paste("Level 1", firstup(adj)),
                                        paste("Level 2", firstup(adj))))

    bal$scenario <- factor(paste(bal$variable, bal$adjustment),
                           nmax = nlevels(bal$variable) * nlevels(bal$adjustment ))

    p <- ggplot(bal) +
      geom_boxplot(aes(x = abs(.data$SMD),
                       y = .data$variable,
                       fill = .data$variable),
                   ...) +
      geom_vline(xintercept = 0) +
      geom_vline(xintercept = .1, linetype = "dashed") +
      geom_vline(xintercept = mean(abs(bal$SMD[bal$variable == "Unmeasured\nConfounder" &
                                                   bal$adjustment == "Unadjusted"])),
                 linetype = "dotted") +
      facet_grid(~adjustment, switch = "y") +
      theme_bw() +
      theme(strip.placement = "outside") +
      guides(fill = "none") +
      labs(x = "Absolute Standardized\nMean Difference", y = NULL)

    p

  }
  else if (type == "hr") {
    est <- do.call("rbind", lapply(x$sim_out, `[[`, "all_HRs"))

    est <- est[est$adjustment != "true",]

    est$adjustment <- factor(est$adjustment, levels = c("L2", "L1", "crude"),
                             labels = c(paste("Level 2", firstup(adj)),
                                        paste("Level 1", firstup(adj)),
                                        "Unadjusted"))

    p <- ggplot(est) +
      geom_boxplot(aes(x = exp(.data$HR),
                       y = .data$adjustment,
                       fill = .data$adjustment),
                   ...) +
      geom_vline(xintercept = 1) +
      geom_vline(xintercept = x$HR_table$HR[x$HR_table$adjustment == "true"],
                 linetype = "dashed") +
      scale_x_log10() +
      scale_fill_grey(start = .55, end = .95) +
      theme_bw() +
      guides(fill = "none") +
      labs(x = "Hazard Ratio", y = NULL)

    p
  }
  else {
    est <- do.call("rbind", lapply(x$sim_out, `[[`, "all_HRs"))

    est <- est[est$adjustment != "true",]

    true <- x$HR_table$HR[x$HR_table$adjustment == "true"]

    est$adjustment <- factor(est$adjustment, levels = c("L2", "L1", "crude"),
                             labels = c(paste("Level 2", firstup(adj)),
                                        paste("Level 1", firstup(adj)),
                                        "Unadjusted"))

    p <- ggplot(est) +
      geom_boxplot(aes(x = (exp(.data$HR) - true) / true,
                       y = .data$adjustment,
                       fill = .data$adjustment),
                   ...) +
      geom_vline(xintercept = 0) +
      scale_x_continuous(labels = scales::label_percent()) +
      scale_fill_grey(start = .55, end = .95) +
      theme_bw() +
      guides(fill = "none") +
      labs(x = "Hazard Ratio Relative Bias", y = NULL)

    p
  }
}