Step 2: plot the paths
The dotted line is the null. The band is a pointwise 95% interval for the bias-adjusted estimate.
On the ITCV path the dashed line marks the critical correlation; the path crosses it exactly at the ITCV.
This vignette walks through the complete workflow:
fitted model -> sensitivity path -> confoundsens object -> plots -> reportusing a real, publicly available dataset.
The darfur data (Hazlett, 2020), distributed with the
sensemakr package, come from a survey of Darfurian
refugees in eastern Chad. The question is whether being directly harmed
in the violence (directlyharmed) changed attitudes toward
peace (peacefactor). Following Cinelli and Hazlett (2020),
we adjust for age, occupation, past voting, household size, sex, and
village fixed effects.
library(confoundvis)
data("darfur", package = "sensemakr")
fit <- lm(peacefactor ~ directlyharmed + age + farmer_dar + herder_dar +
pastvoted + hhsize_darfur + female + village, data = darfur)
coef(summary(fit))["directlyharmed", ]
#> Estimate Std. Error t value Pr(>|t|)
#> 0.0973158193 0.0232565378 4.1844499848 0.0000318234A sensitivity path is the treatment estimate as a function
of the strength \(\lambda\) of a
hypothetical omitted confounder. confoundvis computes paths
for two regression frameworks directly from the fitted model.
Partial \(R^2\) (Cinelli and Hazlett, 2020). \(\lambda\) is the partial \(R^2\) of the confounder with both treatment and outcome:
p_r2 <- sens_path_lm(fit, treatment = "directlyharmed",
lambda = seq(0, 0.3, by = 0.001))
p_r2
#> <confoundsens>
#> n : 301
#> lambda : [0, 0.3]
#> theta : [-0.136029, 0.0973158]
#> level : <none>
#> se : yes
#> t : yes
#> frame : partial_r2 (from lm)Impact threshold (Frank, 2000). \(\lambda\) is the impact \(k = r_{DU}\,r_{YU}\); the path is the adjusted partial correlation:
it <- itcv_lm(fit, treatment = "directlyharmed",
lambda = seq(0, 0.3, by = 0.001))
c(r = it$r, r_crit = it$r_crit, ITCV = it$itcv,
ITCV_unconditional = it$itcv_unconditional)
#> r r_crit ITCV ITCV_unconditional
#> 0.14789555 0.06997985 0.08377850 0.04486840Both results are confoundsens objects, so everything
downstream is the same regardless of framework.
The dotted line is the null. The band is a pointwise 95% interval for the bias-adjusted estimate.
On the ITCV path the dashed line marks the critical correlation; the path crosses it exactly at the ITCV.
A threshold alone does not say whether an omitted confounder that
strong is plausible. The usual aid is to compare it with observed
covariates. covariate_impacts() computes, for every
covariate term, its partial associations with treatment and outcome.
imp <- covariate_impacts(fit, "directlyharmed", metric = "partial_r2")
imp[order(-imp$impact), c("covariate", "term_df", "r2dz.x", "r2yz.dx",
"bias_index")]
#> covariate term_df r2dz.x r2yz.dx bias_index
#> 7 village 485 4.169266e-01 0.4242888431 0.5508064275
#> 6 female 1 9.081065e-03 0.1090339154 0.0316104106
#> 1 age 1 1.115766e-03 0.0080242581 0.0029938602
#> 4 pastvoted 1 1.489741e-03 0.0050684125 0.0027498882
#> 3 herder_dar 1 7.898770e-03 0.0002521486 0.0014168700
#> 5 hhsize_darfur 1 1.740140e-04 0.0005195193 0.0003006981
#> 2 farmer_dar 1 4.917499e-06 0.0024364946 0.0001094601The sensitivity Love plot places each covariate on the chosen metric
and draws the threshold. With metric = "partial_r2", a
covariate to the right of the line would, if an omitted confounder were
equally strong, bring the point estimate to zero.
The ITCV metric uses signed products of partial correlations and is
defined for single-column terms; the multi-column village
term is dropped with a message.
imp_it <- covariate_impacts(fit, "directlyharmed", metric = "itcv")
#> Multi-column terms have no ITCV impact (NA); use metric = "partial_r2" to include them.
plot_sensitivity_love(imp_it)
The same object supplies benchmarks for the ITCV contour:
sens_report(p_r2)
#> Sensitivity report (framework: partial_r2; source: lm)
#> Treatment: directlyharmed
#> Estimate 0.0973 (SE 0.0233, df 783); robustness value RV = 0.139.
#> Confounder strength (lambda): Confounder partial R2 with treatment and
#> outcome.
#> Starting at 0.0973, the point estimate reaches the null at lambda =
#> 0.139; significance (alpha = 0.05) is lost at lambda = 0.0763.
#> These values describe how strong an omitted confounder would need to
#> be. They do not show whether such a confounder exists, and they rest
#> on the framework's own modelling assumptions. Judge plausibility
#> against substantive knowledge and observed benchmarks (see
#> covariate_impacts()).
sens_report(it$path)
#> Sensitivity report (framework: itcv; source: lm)
#> Treatment: directlyharmed
#> Partial correlation r = 0.148; critical r = 0.07; ITCV = 0.0838.
#> Confounder strength (lambda): Confounder impact k = r(D,U) x r(Y,U).
#> Starting at 0.148, the point estimate reaches the null at lambda =
#> 0.148; significance (alpha = 0.05) is lost at lambda = 0.0838.
#> These values describe how strong an omitted confounder would need to
#> be. They do not show whether such a confounder exists, and they rest
#> on the framework's own modelling assumptions. Judge plausibility
#> against substantive knowledge and observed benchmarks (see
#> covariate_impacts()).robustness_points() returns the same numbers as a data
frame for tables. On the partial \(R^2\) path, the null crossing is the
robustness value and the significance crossing is the robustness value
at \(\alpha = 0.05\) reported by
sensemakr.
If an analysis was run with another package, convert its result instead of recomputing it.
s <- sensemakr::sensemakr(fit, treatment = "directlyharmed",
benchmark_covariates = "female", kd = 1:3)
p_sm <- from_sensemakr(s, lambda = seq(0, 0.3, by = 0.001))
robustness_points(p_sm)
#> level theta_start lambda_null lambda_sig sig_rule
#> 1 pooled 0.09731582 0.1387762 0.07625786 t, df = 782
unlist(s$sensitivity_stats[c("rv_q", "rv_qa")])
#> rv_q rv_qa
#> 0.13877635 0.07625797k <- suppressWarnings(suppressMessages(
konfound::konfound(fit, directlyharmed, to_return = "raw_output")))
p_kf <- from_konfound(k, lambda = seq(0, 0.3, by = 0.001))
c(confoundvis = it$itcv, konfound_itcvGz = k$itcvGz)
#> confoundvis konfound_itcvGz
#> 0.0837785 0.0837785For binary or ratio-scale outcomes, E-value results convert to a path on the log risk-ratio scale that crosses the null at the E-value:
e <- EValue::evalues.RR(est = 1.8, lo = 1.3, hi = 2.5)
p_ev <- from_evalue(e, lambda = seq(1, 4, by = 0.01))
plot_robustness_curve(p_ev, points = FALSE)
confoundvis is a presentation layer. Its computations
reproduce the published formulas of each framework (the test suite
checks them against sensemakr,
konfound, and EValue), and its
contribution is a shared object and a shared set of displays. It
inherits every assumption of the framework that generated a path. In
particular:
plot_reversal_cone(),
plot_taylor_panels()) are conceptual tools built on
stylized models and do not analyse data.Cinelli, C., & Hazlett, C. (2020). Making sense of sensitivity: Extending omitted variable bias. Journal of the Royal Statistical Society: Series B, 82(1), 39–67.
Frank, K. A. (2000). Impact of a confounding variable on a regression coefficient. Sociological Methods & Research, 29(2), 147–194.
Frank, K. A., Maroulis, S. J., Duong, M. Q., & Kelcey, B. M. (2013). What would it take to change an inference? Educational Evaluation and Policy Analysis, 35(4), 437–460.
Hazlett, C. (2020). Angry or weary? How violence impacts attitudes toward peace among Darfurian refugees. Journal of Conflict Resolution, 64(5), 844–870.
VanderWeele, T. J., & Ding, P. (2017). Sensitivity analysis in observational research: Introducing the E-value. Annals of Internal Medicine, 167(4), 268–274.