# graphics ----------------------------------------------------------------
##' @title Sensitivity plots
##' @description Several sensitivity plots.
##'
##' @param x sensemade object
##' @param type type of the plot
##' @param ... extra arguments
##'
##'
##' @return a ggplot object with the plot.
##' @export
plot.sensemade <- function(x, type = c("contour", "worst-case"), ...){
type <- match.arg(type)
switch(type,
contour = contourplot(x, ...),
"worst-case" = worstcaseplot(x, ...)
)
}
# A method must have all the arguments of the generic, including … if the generic does.
##' @export
contourplot <- function(x,
contour = c("estimate","t-value", "lower bound", "upper bound"),
nlevels = 15,
pch = 20,
cex = 1,
lim = NULL,
xlab = "Hypothetical partial R2 of unobserved confounder with the treatment",
ylab = "Hypothetical partial R2 of unobserved confounder with the outcome",
main = paste("Sensitivity of", contour, "to unobserved confounder\nContours of adjusted estimates"),
top = 3,
x.label = NULL,
y.label = NULL){
contour <- match.arg(contour)
#### benchmark data ####
benchmarks <- x$benchmarks$benchmark_R2
top <- min(c(nrow(benchmarks), top))
benchmarks <- benchmarks[1:top, ]
r2y <- benchmarks$r2y
r2d <- benchmarks$r2d
benchmarks_group <- x$benchmarks$benchmark_R2_group
r2y_group = benchmarks_group$r2y
r2d_group = benchmarks_group$r2d
####
# NOTE: if-else masking of design matrix columns tied to overall grouping
# design choice, make user mask these low level points
# via 'showvars' argument
# eg, we will always plot the points tied to 'grouping' benchmarks
# but we will not mask the points tied to the levels of the 'grouping' benchmarks
# one way to 'mask' is line 178 of https://github.com/statsccpr/ovb/blob/master/R/ovb.R
#### contour data (treatment effect data) ####
estimate <- x$treat.stats$estimate
se <- x$treat.stats$se
df <- x$treat.stats$df
t <- estimate/se
#### computing contours ####
## sequence of R2y and R2d
if (is.null(lim)) lim <- max(c(r2y, r2d, r2y_group, r2d_group), na.rm = TRUE) + 0.1
s <- seq(0, lim, by = 0.001)
## level curves
if (contour == "estimate") {
# estimate level curve
z <- adjust_estimate(estimate, outer(s, s, get_bias, se = se, df = df))
labels <- paste0(benchmarks$covariate, "\n", "(",round(benchmarks$adj_est_r2, 3),")")
lev <- 0
} else if (contour == "t-value") {
# t-value level curve
z <- outer(s, s, get_t, t = t, df = df)
labels <- paste0(benchmarks$covariate, "\n", "(",round(benchmarks$adj_t_r2, 3),")")
lev <- 2
} else if (contour == "lower bound" | contour == "upper bound" ) {
# CI level curves
new_estimate <- adjust_estimate(estimate, outer(s, s, get_bias, se = se, df = df))
new_se <- outer(s, s, get_se, se = se, df = df)
if (contour == "lower bound") {
# CI lower bound
z <- new_estimate - 1.96*new_se
benchmarks$adj_lw_r2 <- benchmarks$adj_est_r2 - 1.96*benchmarks$adj_se_r2
labs <- benchmarks$adj_lw_r2
} else {
# CI upper bound
z <- new_estimate + 1.96*new_se
benchmarks$adj_up_r2 <- benchmarks$adj_est_r2 + 1.96*benchmarks$adj_se_r2
labs <- benchmarks$adj_up_r2
}
labels <- paste0(benchmarks$covariate, "\n", "(",round(labs, 3),")")
lev <- 0
}
####
## plot contours with benchmarks
contour(s, s, z, nlevels = nlevels,
xlab = xlab,
ylab = ylab,
main = main)
points(r2d, r2y, pch = 23, col = "black", bg = "red", cex = cex)
contour(s, s, z = z, level = lev, add = TRUE, col = "red", lwd = 2, lty = 2)
points(r2d_group, r2y_group, pch = 23, col = "black", bg = "cyan", cex = cex)
#### add labels ####
# todo: try to create a better function for positioning labels and substitute jitter
if (is.null(x.label))
r2dl <- jitter(r2d, factor = 20)
if (is.null(y.label))
r2yl <- jitter(r2y, factor = 20)
text(r2dl, r2yl, labels = labels, cex = 0.7)
if(nrow(benchmarks_group)>0){
labels_group <- benchmarks_group$covariate
text(x=benchmarks_group$r2d,
y=benchmarks_group$r2y,
labels=labels_group,
cex=0.7,
pos=4)
}
labels <- data.frame(labels = labels,x = r2dl, y = r2yl, stringsAsFactors = FALSE)
####
#### data for reproducibility ####
rownames(z) <- colnames(z) <- s
out <- list(plot_type = paste(contour,"contours"),
contours = z,
benchmarks = benchmarks,
benchmarks_group = benchmarks_group,
labels = labels)
####
invisible(out)
}
contourplot
##' @export
worstcaseplot <- function(x,
lim = NULL,
scenarios = c(1, 0.3),
cex.legend = 0.6,
index = NULL,
xlab = "Hypothetical partial R2 of unobserved confounder(s) with treatment",
ylab = "Adjusted estimate",
main = "Sensitivity of estimate to unobserved confounder(s)\n\"Worst-case\" scenarios of partial R2 with outcome"){
benchmarks <- x$benchmarks$benchmark_R2
r2d <- benchmarks$r2d
if (!is.null(index)) r2d <- r2d[index]
estimate <- x$treat.stats$estimate
se <- x$treat.stats$se
df <- x$treat.stats$df
if (is.null(lim)) lim <- max(r2d, na.rm = TRUE) + 0.1
s <- seq(0, lim, by = 0.001)
mr2y <- x$benchmarks$benchmark_all_vars$r2y_all
y <- adjust_estimate(estimate, get_bias(se = se, df = df, r2y = scenarios[1], r2d = s))
out <- data.frame(r2d = s,
r2y = scenarios[1],
adj_est = y)
plot(s, y,
type = "l", bty = "L",
xlab = xlab,
ylab = ylab,
main = main)
abline(h = 0, col = "red", lty = 5)
scenarios2 <- scenarios[-1]
for (i in seq_along(scenarios2)) {
y <- adjust_estimate(estimate, get_bias(se = se, df = df, r2y = scenarios2[i], r2d = s))
out2 <- data.frame(r2d = s,
r2y = scenarios2[i],
adj_est = y)
out <- rbind(out, out2)
lines(s, y, lty = i + 1)
}
p <- length(scenarios)
# max observed R2y
y <- adjust_estimate(estimate, get_bias(se = se, df = df, r2y = mr2y, r2d = s))
lines(s, y, lty = p + 1, col = "red")
out2 <- data.frame(r2d = s,
r2y = mr2y,
adj_est = y)
out <- rbind(out, out2)
out <- list(plot_type = "worst case plot",
adjusted_estimates = out)
legend("topright",
lty = c((1:p), p + 1),
col = c(rep("black", p), "red"),
legend = paste0(c(scenarios*100, round(mr2y*100)), "%"),
ncol = p + 1,
title = "Hypothetical partial R2 of unobserved confounder(s) with outcome",
cex = cex.legend)
rug(x = r2d, col = "red")
invisible(out)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.