update to add bias functions

NAMESPACE

@@ -1,6 +1,10 @@
 # Generated by roxygen2: do not edit by hand
 
 export("%>%")
+export(observed_bias_order)
+export(observed_bias_tbl)
+export(observed_bias_tip)
+export(observed_covariate_e_value)
 export(tip)
 export(tip_b)
 export(tip_c)

R/observed-bias-plot-helpers.R

@@ -0,0 +1,104 @@
+get_y <- function(m) {
+  deparse(stats::formula(m)[[2]])
+}
+
+parse_formula <- function(m) {
+  as.character(
+    attr(m$terms, "variables")
+  )[-c(1,2)]
+
+}
+
+create_covariate_lists <- function(ps_mod, outcome_mod) {
+  exposure <- get_y(ps_mod)
+
+  ps_covariates <- parse_formula(ps_mod)
+  outcome_covariates <- parse_formula(outcome_mod)
+
+  ps_covariates_clean <- unique(clean_covariate(ps_covariates))
+  outcome_covariates_clean <- unique(clean_covariate(outcome_covariates))
+  outcome_covariates_clean <- outcome_covariates_clean[
+    !(outcome_covariates_clean %in% exposure)
+  ]
+  list(exposure = exposure,
+       ps_covariates = ps_covariates,
+       ps_covariates_clean = ps_covariates_clean,
+       outcome_covariates = outcome_covariates,
+       outcome_covariates_clean = outcome_covariates_clean
+  )
+}
+
+drop_one_mod_tbl <- function(cov, names, covariate_lists) {
+  ps_covariates <- covariate_lists[["ps_covariates"]]
+  outcome_covariates <- covariate_lists[["outcome_covariates"]]
+
+  cov_ps <- cov[cov %in% covariate_lists[["ps_covariates_clean"]]]
+  cov_outcome <- cov[cov %in% covariate_lists[["outcome_covariates_clean"]]]
+  if (all(clean_covariate(ps_covariates) %in% cov_ps)) {
+    new_ps = 1
+  } else{
+    new_ps = ps_covariates[
+      !(clean_covariate(ps_covariates) %in% cov_ps)
+    ]
+  }
+  tibble::tibble(
+    dropped = names,
+    new_ps = list(new_ps),
+    new_outcome = list(
+      outcome_covariates[
+        !(clean_covariate(outcome_covariates) %in% cov_outcome)
+      ])
+  )
+}
+
+
+create_individual_covariate_list <- function(covariate_lists) {
+  covs <- as.list(unique(c(covariate_lists[["ps_covariates_clean"]],
+                           covariate_lists[["outcome_covariates_clean"]])))
+  names(covs) <- covs
+  covs
+}
+
+drop_tbl <- function(covs, covariate_lists) {
+
+  g <- purrr::map2(covs, names(covs), drop_one_mod_tbl, covariate_lists)
+  g <- do.call(rbind, g)
+  g$type <- ifelse(purrr::map(covs, length) == 1, "covariate", "group")
+  g
+}
+
+
+add_formula <- function(d, exposure, outcome) {
+  tibble::add_column(
+    d,
+    ps_form = purrr::map(d$new_ps, build_formula, y = exposure),
+    outcome_form = purrr::map(d$new_outcome, build_formula, y = outcome)
+  )
+}
+
+clean_covariate <- function(x) {
+  gsub(".*\\(|\\).*|\\^.*|,.*$", "", x)
+}
+
+build_formula <- function(y, x) {
+  covs <- glue::glue_collapse(x, sep = "+")
+  stats::as.formula(
+    glue::glue("{y} ~ {covs}")
+  )
+}
+
+check_drop_list <- function(l) {
+  if (!is.null(l)) {
+    n <- names(l)
+    if (length(n) != length(l)) {
+      stop_glue("`drop_list` must be a named list.")
+    }
+    c <- purrr::map_lgl(l, is.character)
+    if (!all(c)) {
+      stop_glue("`drop_list` must be a named list of character vectors.")
+    }
+  }
+}
+
+
+

R/observed_bias_order.R

@@ -0,0 +1,21 @@
+#' Order observed bias data frame for plotting
+#'
+#' @param d Observed bias data frame. Must have columns `dropped` and `type`
+#' @param by Character. Variable in `d` to order by.
+#'
+#' @return Data frame in the correct order
+#' @export
+observed_bias_order <- function(d, by) {
+  grps_ <- d[d$type == "group" & !grepl("Hypothetical", d$dropped), ]
+  grps <- which(d$type == "group" & !grepl("Hypothetical", d$dropped))
+  grps <- grps[order(grps_[[by]], decreasing = TRUE)]
+
+  hypo_ <- d[d$type == "tip", ]
+  hypo <- which(d$type == "tip")
+  hypo <- hypo[order(hypo_[[by]])]
+
+  d <- d[c(hypo, grps, order(d[[by]][d$type == "covariate"], decreasing = TRUE)), ]
+  d$dropped <- factor(d$dropped,
+                      levels = d$dropped)
+  d
+}

R/observed_bias_tbl.R

@@ -0,0 +1,55 @@
+#' Create a data frame to assist with creating an observed bias plot
+#'
+#' @param ps_mod Model object for the propensity score model
+#' @param outcome_mod Model object for the outcome model
+#' @param drop_list Named list of covariates or groups of covariates to drop if
+#'   `NULL`, will default to dropping each covariate one at a time.
+#'
+#' @return Data frame with the following columns:
+#'   * `dropped`: The covariate or group of covariates that were dropped
+#'   * `type`: Explanation of `dropped`, whether it refers to a single covariate (`covariate`) or a group of covariates (`group`)
+#'   * `ps_formula`: The new formula for the updated propensity score model
+#'   * `outcome_formula`: The new formula for the updated outcome model
+#'   * `ps_model`: The new model object for the updated propensity score model
+#'   * `p`: The updated propensity score
+#' @export
+#'
+#' @examples
+#' ps_mod <- glm(am ~ mpg + cyl + I(hp^2), data = mtcars)
+#' outcome_mod <- lm(qsec ~ am + hp + disp + wt, data = mtcars)
+#' observed_bias_tbl(
+#'  ps_mod,
+#'  outcome_mod,
+#'  drop_list = list(
+#'    group_one = c("mpg", "hp"),
+#'    group_two = c("cyl", "wt")
+#'    )
+#'  )
+
+observed_bias_tbl <- function(ps_mod, outcome_mod, drop_list = NULL) {
+  c <- create_covariate_lists(ps_mod, outcome_mod)
+
+  if (is.null(drop_list)) {
+    drop_list <- create_individual_covariate_list(c)
+  }
+
+  check_drop_list(drop_list)
+  outcome <- get_y(outcome_mod)
+
+  g <- drop_tbl(drop_list, c)
+  d <- add_formula(g, c[["exposure"]], outcome)
+
+  observed_bias_tbl <- tibble::tibble(
+    dropped = d$dropped,
+    type = d$type,
+    ps_formula = d$ps_form,
+    outcome_formula = d$outcome_form,
+    ps_model = purrr::map(d$ps_form, ~ stats::update(ps_mod, .x))
+  )
+
+  tibble::add_column(
+    observed_bias_tbl,
+    p = purrr::map(observed_bias_tbl$ps_model,
+                   stats::predict, type = "response")
+  )
+}

R/observed_bias_tip.R

@@ -0,0 +1,26 @@
+#' Create a data frame to combine with an observed bias data frame demonstrating a hypothetical unmeasured confounder
+#'
+#' @param tip Numeric. Value you would like to tip to.
+#' @param point_estimate Numeric. Result estimate from the full model.
+#' @param lb Numeric. Result lower bound from the full model.
+#' @param ub Numeric. Result upper bound from the full model.
+#' @param tip_desc Character. A description of the tipping point.
+#'
+#' @return A data frame with five columns:
+#' * `dropped`: the input from `tip_desc`
+#' * `type`: Explanation of `dropped`, here `tip` to clarify that this was calculated as a tipping point.
+#' * `point_estimate`: the shifted point estimate
+#' * `lb`: the shifted lower bound
+#' * `ub`: the shifted upper bound
+#' @export
+#'
+observed_bias_tip <- function(tip, point_estimate, lb, ub, tip_desc = "Hypothetical unmeasured confounder") {
+  shift <- 1 - tip
+  tibble::tibble(
+    dropped = tip_desc,
+    type = "tip",
+    point_estimate = point_estimate + shift,
+    lb = lb + shift,
+    ub = ub + shift
+  )
+}

R/observed_covariate_e_value.R

@@ -0,0 +1,43 @@
+#' Calculate the Observed Covariate E-value
+#'
+#' @param lb Numeric. The lower bound of the full model
+#' @param ub Numeric. The upper bound of the full model
+#' @param lb_adj Numeric. The lower bound of the adjusted model
+#' @param ub_adj Numeric. The upper bound of the adjusted model
+#' @param transform Character. If your effect is an odds ratio or hazard ratio, this will
+#' perform the transformation suggested by VanderWeele and Ding. Allowed values are:
+#'  * "OR"
+#'  * "HR"
+#'
+#' @return The Observed Covariate E-value
+#' @export
+observed_covariate_e_value <- function(lb, ub, lb_adj, ub_adj, transform = NULL) {
+  if (!is.null(transform)) {
+    if (!transform %in% c("OR", "HR")) {
+      stop_glue("You input\n * `transform`: {transform}\n ",
+                "The only valid `transform` inputs are\n * 'HR'\n * 'OR'")
+    }
+    if (transform == "OR") {
+      lb <- sqrt(lb)
+      ub <- sqrt(ub)
+      lb_adj <- sqrt(lb_adj)
+      ub_adj <- sqrt(ub_adj)
+    }
+    if (transform == "HR") {
+      lb <- hr_transform(lb)
+      ub <- hr_transform(ub)
+      lb_adj <- hr_transform(lb_adj)
+      ub_adj <- hr_transform(ub_adj)
+    }
+  }
+  b <- get_limiting_bound(lb, ub)
+  b_adj <- get_limiting_bound_adj(b, lb_adj, ub_adj)
+  if (b < 1) {
+    b <- 1 / b
+    b_adj <- 1 / b_adj
+  }
+  if (b < b_adj) {
+    return((b_adj / b) + sqrt((b_adj / b) * ((b_adj / b) - 1)))
+  }
+  (b / b_adj) + sqrt((b / b_adj) * ((b / b_adj) - 1))
+}

R/tip-helpers.R

@@ -163,9 +163,9 @@ tip_n <- function(p0, p1, gamma, b) {
   as.numeric(n)
 }
 
-# e_value <- function(lb, ub) {
-#   observed_covariate_e_value(lb, ub, 1, 1)
-# }
+e_value <- function(lb, ub) {
+  observed_covariate_e_value(lb, ub, 1, 1)
+}
 
 hr_transform <- function(hr) {
   (1 - (0.5^sqrt(hr))) / (1 - (0.5^sqrt(1 / hr)))