From a fitted model to a sensitivity report

Data and model

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.0000318234

Step 1: compute sensitivity paths

A 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.04486840

Both results are confoundsens objects, so everything downstream is the same regardless of framework.

Step 2: plot the paths

plot_robustness_curve(p_r2, points = FALSE)

The dotted line is the null. The band is a pointwise 95% interval for the bias-adjusted estimate.

plot_robustness_curve(it$path, points = FALSE)

On the ITCV path the dashed line marks the critical correlation; the path crosses it exactly at the ITCV.

Step 3: benchmark against observed covariates

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.0001094601

The 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.

plot_sensitivity_love(imp)

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 contour plot. Passing the itcv_lm() result draws the exact significance boundary (solid) next to the product reference \(|r_{YU} r_{DU}| = \text{ITCV}\) (dashed). The two agree only for confounders whose two correlations have equal magnitude; elsewhere the exact boundary decides. Pairs that no correlation matrix can produce are shaded as infeasible.

plot_sensitivity_contour(itcv = it,
                         benchmarks = imp_it[!is.na(imp_it$r_yu), ],
                         limits = c(-0.5, 0.5))

Step 4: a plain-language report

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.

Using results from other packages

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.07625797

k <- 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.0837785

For 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)

Matched designs: Rosenbaum bounds

For a matched-pair study, sens_path_pairs() computes the path under Rosenbaum’s sensitivity model, in which two matched units may differ in their odds of treatment by at most a factor \(\Gamma\). The path stores, for each \(\Gamma\), the deviate equivalent to the upper bound on the one-sided P-value, so the same crossing, plotting, and reporting functions apply. The werfel data in DOS2 compare DNA damage in 39 welders and 39 controls matched for age and smoking.

data("werfel", package = "DOS2")
d <- werfel$serpc_p - werfel$cerpc_p
p_gam <- sens_path_pairs(d, gamma = seq(1, 6, by = 0.01))
robustness_points(p_gam)
#>    level theta_start lambda_null lambda_sig        sig_rule
#> 1 pooled    4.535374          NA    3.86761 theta threshold
sens_report(p_gam)
#> Sensitivity report (framework: rosenbaum; source: pairs)
#> Matched pairs: 39 (39 with nonzero difference); Wilcoxon signed rank
#>   statistic = 715; sensitivity value Gamma = 3.87 (one-sided alpha =
#>   0.05).
#> Confounder strength (lambda): Gamma (odds of differential treatment
#>   assignment).
#> The upper-bound P-value reaches 0.05 at Gamma = 3.87; it stays below
#>   0.5 for Gamma <= 6.
#> 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, where available, observed
#>   benchmarks.
plot_robustness_curve(p_gam, points = FALSE)

The result agrees with DOS2::senWilcox(). For a binary outcome, use method = "binary", which applies the exact bound for McNemar’s statistic.

Two-parameter surfaces

Each framework describes an omitted confounder with two parameters, one for its association with treatment and one for its association with the outcome. A path is one line through that plane. sens_surface() returns the whole plane with the boundaries at which the conclusion changes, computed from the framework’s published formula. The partial \(R^2\) surface is the plane of sensemakr::ovb_contour_plot(); benchmarks from covariate_impacts() are placed on it automatically.

s_r2 <- sens_surface(p_r2)
s_r2
#> <confoundsurf>
#>   framework : partial_r2 
#>   a (x-axis): Partial R2 of U with treatment [0, 0.416329] 
#>   b (y-axis): Partial R2 of U with outcome [0, 0.416329] 
#>   value     : Adjusted estimate 
#>   boundaries: Null (point estimate = 0); Significance (alpha = 0.05); Significance, reversed direction (alpha = 0.05)
plot_sensitivity_surface(s_r2, benchmarks = head(imp[order(-imp$impact), ], 4))

For a matched-pair study, the surface shows the amplification of \(\Gamma\) (Rosenbaum and Silber, 2009): every pair \((\Lambda, \Delta)\) on the solid curve corresponds to the sensitivity value, as in sensitivitymv::amplify().

plot_sensitivity_surface(sens_surface(p_gam))

The size of a shaded zone depends on the axis ranges and carries no inferential meaning; the surface is read only through its boundaries.

Several analyses in one table

sens_table() collects the decision information of several paths, each in its own units. It is intended for reporting several analyses side by side, not for ranking them: a robustness value, an impact threshold, and a value of \(\Gamma\) measure different things.

sens_table(partial_r2 = p_r2, itcv = it$path, sensemakr = p_sm)
#> <confoundvis_table>
#>  path       framework  source   
#>  partial_r2 partial_r2 lm       
#>  itcv       itcv       lm       
#>  sensemakr  partial_r2 sensemakr
#>  parameter                                        summary_name      
#>  Confounder partial R2 with treatment and outcome RV (point null)   
#>  Confounder impact k = r(D,U) x r(Y,U)            ITCV (conditional)
#>  Confounder partial R2 with treatment and outcome RV (point null)   
#>  summary_value theta_start lambda_null lambda_sig sig_rule       
#>  0.13880       0.09732     0.1388      0.07626    t, df = 782    
#>  0.08378       0.14790     0.1479      0.08378    theta threshold
#>  0.13880       0.09732     0.1388      0.07626    t, df = 782    
#> 
#> Each row keeps its own sensitivity parameter (column `parameter`); lambda values are comparable only within a row.

What this package does and does not do

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:

A sensitivity path answers “how strong would confounding have to be?”, not “is there such a confounder?” Only substantive knowledge can judge that.
Covariate benchmarks describe observed variables. An omitted confounder can be stronger than every measured covariate, and benchmarks can themselves be distorted when the omitted confounder is correlated with them.
The partial \(R^2\) and ITCV paths are for linear regression coefficients. The E-value path concerns ratio measures; its bias factor is a bound, not an exact adjustment.
No display repairs a flawed identification strategy (for example, bad controls, selection, or measurement error in the treatment).
confoundvis adds no new sensitivity parameter or measure. Every number it reports is the published quantity of the framework that produced it.

References

Rosenbaum, P. R., & Silber, J. H. (2009). Amplification of sensitivity analysis in matched observational studies. Journal of the American Statistical Association, 104(488), 1398–1405.

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.

Rosenbaum, P. R. (2002). Observational Studies (2nd ed.). Springer.

VanderWeele, T. J., & Ding, P. (2017). Sensitivity analysis in observational research: Introducing the E-value. Annals of Internal Medicine, 167(4), 268–274.