R/tip.R

In LucyMcGowan/tipr: Tipping Point Analyses

#' Tip a result with a normally distributed confounder.
#'
#' choose one of the following, and the other will be estimated:
#' * `exposure_confounder_effect`
#' * `confounder_outcome_effect`
#'
#' @param effect_observed Numeric positive value. Observed exposure - outcome effect
#'    (assumed to be the exponentiated coefficient, so a risk ratio, odds
#'    ratio, or hazard ratio). This can be the point estimate, lower confidence
#'    bound, or upper confidence bound.
#' @param exposure_confounder_effect Numeric. Estimated difference in scaled means between the
#'    unmeasured confounder in the exposed population and unexposed population
#' @param confounder_outcome_effect Numeric positive value. Estimated relationship
#'    between the unmeasured confounder and the outcome
#' @param verbose Logical. Indicates whether to print informative message.
#'    Default: `TRUE`
#' @param correction_factor Character string. Options are "none", "hr", "or".
#'   For common outcomes (>15%), the odds ratio or hazard ratio is not a good
#'   estimate for the risk ratio. In these cases, we can apply a correction
#'   factor. If you are supplying a hazard ratio for a common outcome, set
#'   this to "hr"; if you are supplying an odds ratio for a common outcome, set
#'   this to "or"; if you are supplying a risk ratio or your outcome is rare,
#'   set this to "none" (default).
#'
#'
#' @return Data frame.
#'
#' @examples
#' ## to estimate the relationship between an unmeasured confounder and outcome
#' ## needed to tip analysis
#' tip(1.2, exposure_confounder_effect = -2)
#'
#' ## to estimate the number of unmeasured confounders specified needed to tip
#' ## the analysis
#' tip(1.2, exposure_confounder_effect = -2, confounder_outcome_effect = .99)
#'
#' ## Example with broom
#' if (requireNamespace("broom", quietly = TRUE) &&
#'     requireNamespace("dplyr", quietly = TRUE)) {
#'   glm(am ~ mpg, data = mtcars, family = "binomial") %>%
#'    broom::tidy(conf.int = TRUE, exponentiate = TRUE) %>%
#'    dplyr::filter(term == "mpg") %>%
#'    dplyr::pull(conf.low) %>%
#'    tip(confounder_outcome_effect = 2.5)
#'}
#' @export
tip <- function(effect_observed, exposure_confounder_effect = NULL, confounder_outcome_effect = NULL,
                verbose = getOption("tipr.verbose", TRUE), correction_factor = "none") {

  check_arguments(
    "tip()",
    exposure_confounder_effect,
    confounder_outcome_effect
  )

  exposure_confounder_effect <- exposure_confounder_effect %||% list(NULL)
  confounder_outcome_effect <- confounder_outcome_effect %||% list(NULL)

  o <- purrr::pmap(
    list(b = effect_observed,
         exposure_confounder_effect = exposure_confounder_effect,
         confounder_outcome_effect = confounder_outcome_effect,
         verbose = verbose,
         correction_factor = correction_factor),
    tip_one
  )
  do.call(rbind, o)
}

tip_one <- function(b, exposure_confounder_effect, confounder_outcome_effect, verbose, correction_factor) {
  check_effect(b)
  check_gamma(confounder_outcome_effect)

  correction <- ""
  if (correction_factor == "hr") {
    b <- hr_transform(b)
    confounder_outcome_effect <- hr_transform(confounder_outcome_effect)
    correction <- 'You opted to use the hazard ratio correction to convert your hazard ratios to approximate risk ratios.\nThis is a good idea if the outcome is common (>15%).'
  }

  if (correction_factor == "or") {
    b <- or_transform(b)
    confounder_outcome_effect <- or_transform(confounder_outcome_effect)
    correction <- 'You opted to use the odds ratio correction to convert your odds ratios to approximate risk ratios.\nThis is a good idea if the outcome is common (>15%).'
  }

  n_unmeasured_confounders <- 1

  if (is.null(confounder_outcome_effect)) {
    confounder_outcome_effect <- b ^ (1 / exposure_confounder_effect)
  } else if (is.null(exposure_confounder_effect)) {
    exposure_confounder_effect <- log(b) / log(confounder_outcome_effect)
  } else {
    n_unmeasured_confounders <-
      log(b) / (exposure_confounder_effect * log(confounder_outcome_effect))
    if (any(n_unmeasured_confounders < 0)) {
      if (length(exposure_confounder_effect) > 1) {
        exposure_confounder_effects <- exposure_confounder_effect[n_unmeasured_confounders < 0]
      } else {
        exposure_confounder_effects <- exposure_confounder_effect
      }
      if (length(confounder_outcome_effect) > 1) {
        confounder_outcome_effects <- confounder_outcome_effect[n_unmeasured_confounders < 0]
      } else {
        confounder_outcome_effects <- confounder_outcome_effect
      }

      warning_cli(c(
        "!" = "The observed effect {b} would not tip with the unmeasured \\
        confounder given:",
        "*" = "`exposure_confounder_effect`: {exposure_confounder_effects}",
        "*" = "`confounder_outcome_effect`: {confounder_outcome_effects}"
      ))
      n_unmeasured_confounders <- max(0, n_unmeasured_confounders)
    }
    too_small <-
      n_unmeasured_confounders < 1 & n_unmeasured_confounders > 0
    if (any(too_small)) {
      exposure_confounder_effects <- ifelse(length(exposure_confounder_effect) > 1, exposure_confounder_effect[too_small], exposure_confounder_effect)
      confounder_outcome_effects <-
        ifelse(length(confounder_outcome_effect) > 1,
               confounder_outcome_effect[too_small],
               confounder_outcome_effect)
      warning_cli(c(
        "!" = "The observed effect {b} would tip with < 1 of the given \\
        unmeasured confounders:",
        "*" = "`exposure_confounder_effect`: {exposure_confounder_effects}",
        "*" = "`confounder_outcome_effect`: {confounder_outcome_effects}"
      ))
    }
  }
  o <- tibble::tibble(
    effect_adjusted = 1,
    effect_observed = b,
    exposure_confounder_effect = exposure_confounder_effect,
    confounder_outcome_effect = confounder_outcome_effect,
    n_unmeasured_confounders = n_unmeasured_confounders
  )
  if (verbose) {
    if (all(o$n_unmeasured_confounders == 0)) {
      o_notip <- o[o$n_unmeasured_confounders == 0,]
      message_cli(c(
        "i" = "The observed effect ({round(o_notip$effect_observed, 2)}) \\
        cannot be tipped by an unmeasured confounder\nwith the \\
        following specifications:",
        "*" = "estimated difference in scaled means between the \\
        unmeasured confounder\n     in the exposed population and \\
        unexposed population: {round(o_notip$exposure_confounder_effect, 2)}",
        "*" = "estimated relationship between the unmeasured confounder and \\
        the outcome: {round(o_notip$confounder_outcome_effect, 2)}"
      ))

      if (correction != "") message_cli(c("i" = correction))

    } else if (any(o$n_unmeasured_confounders == 0)) {
      o_notip <- o[o$n_unmeasured_confounders == 0,]
      message_cli(c(
        "i" = "The observed effect ({round(o_notip$effect_observed, 2)}) \\
        cannot be tipped by an unmeasured confounder with the \\
        following specifications:",
        "*" = "estimated difference in scaled means between the \\
        unmeasured confounder in the exposed population and \\
        unexposed population: {round(o_notip$exposure_confounder_effect, 2)}",
        "*" = "estimated relationship between the unmeasured confounder and \\
        the outcome: {round(o_notip$confounder_outcome_effect, 2)}"
      ))

      if (correction != "") message_cli(c("i" = correction))

      o_tip <- o[o$n_unmeasured_confounders != 0,]
      message_cli(c(
        "i" = "The observed effect ({round(o_tip$effect_observed, 2)}) WOULD \\
        be tipped by {round(o$n_unmeasured_confounders)} \\
        unmeasured confounder{ifelse(o_tip$n_unmeasured_confounders > 1, 's', '')} \\
        with the following specifications:",
        "*" = "estimated difference in scaled means between the \\
        unmeasured confounder in the exposed population and \\
        unexposed population: {round(o_tip$exposure_confounder_effect, 2)}",
        "*" = "estimated relationship between the unmeasured confounder and \\
        the outcome: {round(o_tip$confounder_outcome_effect, 2)}"
      ))

      if (correction != "") message_cli(c("i" = correction))

    } else {
      message_cli(c(
        "i" = "The observed effect ({round(o$effect_observed, 2)}) WOULD \\
        be tipped by {round(o$n_unmeasured_confounders)} \\
        unmeasured confounder{ifelse(o$n_unmeasured_confounders > 1, 's', '')} \\
        with the following specifications:",
        "*" = "estimated difference in scaled means between the \\
        unmeasured confounder in the exposed population and \\
        unexposed population: {round(o$exposure_confounder_effect, 2)}",
        "*" = "estimated relationship between the unmeasured confounder and \\
        the outcome: {round(o$confounder_outcome_effect, 2)}"
      ))

      if (correction != "") message_cli(c("i" = correction))
    }
  }
  o
}

#' Tip an observed risk ratio with a normally distributed confounder.
#'
#' choose one of the following, and the other will be estimated:
#' * `exposure_confounder_effect`
#' * `confounder_outcome_effect`
#'
#' @param effect_observed Numeric positive value. Observed exposure - outcome
#'    risk ratio. This can be the point estimate, lower confidence bound,
#'    or upper confidence bound.
#' @param exposure_confounder_effect Numeric. Estimated difference in scaled means between the
#'    unmeasured confounder in the exposed population and unexposed population
#' @param confounder_outcome_effect Numeric positive value. Estimated relationship
#'    between the unmeasured confounder and the outcome
#' @param verbose Logical. Indicates whether to print informative message.
#'    Default: `TRUE`
#'
#'
#' @return Data frame.
#'
#' @examples
#' ## to estimate the relationship between an unmeasured confounder and outcome
#' ## needed to tip analysis
#' tip_rr(1.2, exposure_confounder_effect = -2)
#'
#' ## to estimate the number of unmeasured confounders specified needed to tip
#' ## the analysis
#' tip_rr(1.2, exposure_confounder_effect = -2, confounder_outcome_effect = .99)
#'
#' @export
tip_rr <- function(effect_observed, exposure_confounder_effect = NULL, confounder_outcome_effect = NULL, verbose = getOption("tipr.verbose", TRUE)) {
  check_arguments(
    "tip_rr()",
    exposure_confounder_effect,
    confounder_outcome_effect
  )

  tip(
    effect_observed,
    exposure_confounder_effect = exposure_confounder_effect,
    confounder_outcome_effect = confounder_outcome_effect,
    verbose = verbose
  )
}


#' Tip an observed hazard ratio with a normally distributed confounder.
#'
#' choose one of the following, and the other will be estimated:
#' * `exposure_confounder_effect`
#' * `confounder_outcome_effect`
#'
#' @param effect_observed Numeric positive value. Observed exposure - outcome hazard ratio.
#'    This can be the point estimate, lower confidence bound, or upper
#'    confidence bound.
#' @param exposure_confounder_effect Numeric. Estimated difference in scaled means between the
#'    unmeasured confounder in the exposed population and unexposed population
#' @param confounder_outcome_effect Numeric positive value. Estimated relationship
#'    between the unmeasured confounder and the outcome
#' @param verbose Logical. Indicates whether to print informative message.
#'    Default: `TRUE`
#' @param hr_correction Logical. Indicates whether to use a correction factor.
#'    The methods used for this function are based on risk ratios. For rare
#'    outcomes, a hazard ratio approximates a risk ratio. For common outcomes,
#'    a correction factor is needed. If you have a common outcome (>15%),
#'    set this to `TRUE`. Default: `FALSE`.
#'
#' @return Data frame.
#'
#' @examples
#' ## to estimate the relationship between an unmeasured confounder and outcome
#' ## needed to tip analysis
#' tip_hr(1.2, exposure_confounder_effect = -2)
#'
#' ## to estimate the number of unmeasured confounders specified needed to tip
#' ## the analysis
#' tip_hr(1.2, exposure_confounder_effect = -2, confounder_outcome_effect = .99)
#'
#' @export
tip_hr <- function(effect_observed, exposure_confounder_effect = NULL, confounder_outcome_effect = NULL, verbose = getOption("tipr.verbose", TRUE), hr_correction = FALSE) {
  check_arguments(
    "tip_hr()",
    exposure_confounder_effect,
    confounder_outcome_effect
  )
  correction_factor <- ifelse(hr_correction, "hr", "none")
  tip(
    effect_observed,
    exposure_confounder_effect = exposure_confounder_effect,
    confounder_outcome_effect = confounder_outcome_effect,
    verbose = verbose,
    correction_factor = correction_factor
  )
}

#' Tip an observed odds ratio with a normally distributed confounder.
#'
#' choose one of the following, and the other will be estimated:
#' * `exposure_confounder_effect`
#' * `confounder_outcome_effect`
#'
#' @param effect_observed Numeric positive value. Observed exposure - outcome odds ratio.
#'    This can be the point estimate, lower confidence bound, or upper
#'    confidence bound.
#' @param exposure_confounder_effect Numeric. Estimated difference in scaled means between the
#'    unmeasured confounder in the exposed population and unexposed population
#' @param confounder_outcome_effect Numeric positive value. Estimated relationship
#'    between the unmeasured confounder and the outcome
#' @param verbose Logical. Indicates whether to print informative message.
#'    Default: `TRUE`
#' @param or_correction Logical. Indicates whether to use a correction factor.
#'    The methods used for this function are based on risk ratios. For rare
#'    outcomes, an odds ratio approximates a risk ratio. For common outcomes,
#'    a correction factor is needed. If you have a common outcome (>15%),
#'    set this to `TRUE`. Default: `FALSE`.
#'
#' @return Data frame.
#'
#' @examples
#' ## to estimate the relationship between an unmeasured confounder and outcome
#' ## needed to tip analysis
#' tip_or(1.2, exposure_confounder_effect = -2)
#'
#' ## to estimate the number of unmeasured confounders specified needed to tip
#' ## the analysis
#' tip_or(1.2, exposure_confounder_effect = -2, confounder_outcome_effect = .99)
#'
#' ## Example with broom
#' if (requireNamespace("broom", quietly = TRUE) &&
#'     requireNamespace("dplyr", quietly = TRUE)) {
#'   glm(am ~ mpg, data = mtcars, family = "binomial") %>%
#'    broom::tidy(conf.int = TRUE, exponentiate = TRUE) %>%
#'    dplyr::filter(term == "mpg") %>%
#'    dplyr::pull(conf.low) %>%
#'    tip_or(confounder_outcome_effect = 2.5, or_correction = TRUE)
#'}
#' @export
tip_or <- function(effect_observed, exposure_confounder_effect = NULL, confounder_outcome_effect = NULL, verbose = getOption("tipr.verbose", TRUE), or_correction = FALSE) {
  check_arguments(
    "tip_or()",
    exposure_confounder_effect,
    confounder_outcome_effect
  )

  correction_factor <- ifelse(or_correction, "or", "none")

  tip(
    effect_observed,
    exposure_confounder_effect = exposure_confounder_effect,
    confounder_outcome_effect = confounder_outcome_effect,
    verbose = verbose,
    correction_factor = correction_factor
  )
}

#' @rdname tip_rr
#' @export
tip_rr_with_continuous <- tip_rr

#' @rdname tip_hr
#' @export
tip_hr_with_continuous <- tip_hr

#' @rdname tip_or
#' @export
tip_or_with_continuous <- tip_or

#' @rdname tip
#' @export
tip_with_continuous <- tip

#' @rdname tip
#' @export
tip_c <- tip