Nothing
tests/testthat/test-02-darfur.R
In sensemakr: Sensitivity Analysis Tools for Regression Models
context("Testing Darfur Example")
# runs regression model
data("darfur")
model2 <- lm(peacefactor ~ directlyharmed + age + farmer_dar + herder_dar +
pastvoted + hhsize_darfur + female + village,
data = within(darfur, directlyharmed <- directlyharmed*(-1)))
model <- lm(peacefactor ~ directlyharmed + age + farmer_dar + herder_dar +
pastvoted + hhsize_darfur + female + village, data = darfur)
test_that("testing darfur data", {
expect_equal(dim(darfur), c(1276, 14))
expect_equal(colnames(darfur),
c("wouldvote",
"peacefactor",
"peace_formerenemies",
"peace_jjindiv",
"peace_jjtribes",
"gos_soldier_execute",
"directlyharmed",
"age",
"farmer_dar",
"herder_dar",
"pastvoted",
"hhsize_darfur",
"village",
"female"))
})
test_that(desc = "testing darfur sensemakr",
{
darfur_out <- sensemakr(model, treatment = "directlyharmed", benchmark_covariates = "female", kd = 1:3)
plot(darfur_out)
ovb_contour_plot(darfur_out)
ovb_extreme_plot(darfur_out)
# info
expect_equal(darfur_out$info$formula, peacefactor ~ directlyharmed + age + farmer_dar + herder_dar +
pastvoted + hhsize_darfur + female + village)
expect_equal(darfur_out$info$treatment, "directlyharmed")
expect_equal(darfur_out$info$q, 1)
expect_equal(darfur_out$info$alpha, 0.05)
expect_equal(darfur_out$info$reduce, TRUE)
# sensitivity stats
expect_equal(darfur_out$sensitivity_stats$dof, 783)
expect_equal(darfur_out$sensitivity_stats$r2yd.x, 0.02187, tolerance = 1e-5)
expect_equivalent(c(darfur_out$sensitivity_stats$rv_q), 0.13878, tolerance = 1e-5)
expect_equivalent(c(darfur_out$sensitivity_stats$rv_qa), 0.07626, tolerance = 1e-5)
# bounds
check_bounds <- structure(list(bound_label = c("1x female", "2x female", "3x female"),
r2dz.x = c(0.00916428667504862, 0.0183285733500972, 0.0274928600251459),
r2yz.dx = c(0.12464092303637, 0.249324064199975, 0.374050471038094),
treatment = rep("directlyharmed", 3),
adjusted_estimate = c(0.0752202712144491, 0.0529151723844518, 0.0303960234641548),
adjusted_se = c(0.0218733277437572, 0.0203500620779637, 0.0186700648170924),
adjusted_t = c(3.43890386024675, 2.60024623913809, 1.62806202131271),
adjusted_lower_CI = c(0.032282966, 0.012968035,-0.006253282),
adjusted_upper_CI = c(0.11815758, 0.09286231, 0.06704533)
),
.Names = c("bound_label", "r2dz.x", "r2yz.dx", "treatment",
"adjusted_estimate", "adjusted_se", "adjusted_t",
"adjusted_lower_CI", "adjusted_upper_CI"),
row.names = c(NA, -3L), class = c("ovb_bounds", "data.frame"))
expect_equivalent(darfur_out$bounds, check_bounds)
out1 <- capture.output(darfur_out)
out1 <- capture.output(summary(darfur_out))
out3 <- capture.output(ovb_minimal_reporting(darfur_out))
darfur_out2 <- sensemakr(model, treatment = "directlyharmed")
plot(darfur_out2)
ovb_contour_plot(darfur_out2)
ovb_extreme_plot(darfur_out2)
out3 <- capture.output(ovb_minimal_reporting(darfur_out2))
})
test_that(desc = "testing darfur sensemakr but negative",
{
darfur_out <- sensemakr(model2, treatment = "directlyharmed", benchmark_covariates = "female", kd = 1:3)
plot(darfur_out)
ovb_contour_plot(darfur_out)
ovb_extreme_plot(darfur_out)
# info
expect_equal(darfur_out$info$formula, peacefactor ~ directlyharmed + age + farmer_dar + herder_dar +
pastvoted + hhsize_darfur + female + village)
expect_equal(darfur_out$info$treatment, "directlyharmed")
expect_equal(darfur_out$info$q, 1)
expect_equal(darfur_out$info$alpha, 0.05)
expect_equal(darfur_out$info$reduce, TRUE)
# sensitivity stats
expect_equal(darfur_out$sensitivity_stats$dof, 783)
expect_equal(darfur_out$sensitivity_stats$r2yd.x, 0.02187, tolerance = 1e-5)
expect_equivalent(c(darfur_out$sensitivity_stats$rv_q), 0.13878, tolerance = 1e-5)
expect_equivalent(c(darfur_out$sensitivity_stats$rv_qa), 0.07626, tolerance = 1e-5)
# bounds
check_bounds <- structure(list(bound_label = c("1x female", "2x female", "3x female"),
r2dz.x = c(0.00916428667504862, 0.0183285733500972, 0.0274928600251459),
r2yz.dx = c(0.12464092303637, 0.249324064199975, 0.374050471038094),
treatment = rep("directlyharmed", 3),
adjusted_estimate = c(-0.0752202712144491, -0.0529151723844518, -0.0303960234641548),
adjusted_se = c(0.0218733277437572, 0.0203500620779637, 0.0186700648170924),
adjusted_t = c(-3.43890386024675, -2.60024623913809, -1.62806202131271),
adjusted_lower_CI = -1*c(0.11815758, 0.09286231, 0.06704533),
adjusted_upper_CI = -1*c(0.032282966, 0.012968035,-0.006253282)
),
.Names = c("bound_label", "r2dz.x", "r2yz.dx", "treatment",
"adjusted_estimate", "adjusted_se", "adjusted_t",
"adjusted_lower_CI", "adjusted_upper_CI"),
row.names = c(NA, -3L), class = c("ovb_bounds", "data.frame"))
expect_equivalent(darfur_out$bounds, check_bounds)
out1 <- capture.output(darfur_out)
out1 <- capture.output(summary(darfur_out))
out3 <- capture.output(ovb_minimal_reporting(darfur_out))
darfur_out2 <- sensemakr(model, treatment = "directlyharmed")
plot(darfur_out2)
ovb_contour_plot(darfur_out2)
ovb_extreme_plot(darfur_out2)
out3 <- capture.output(ovb_minimal_reporting(darfur_out2))
})
test_that("testing darfur manual bounds",
{
sense.out <- sensemakr(model, treatment = "directlyharmed", benchmark_covariates = "female", r2dz.x = .1)
bounds.check <- sense.out$bounds
to_check <- bounds.check$adjusted_se[1]
true_check <- adjusted_se(model, treatment = "directlyharmed", r2dz.x = .1, r2yz.dx = .1)
expect_equal(to_check, unname(true_check))
}
)
test_that(desc = "testing darfur sensemakr with formula",
{
darfur_out <- sensemakr(formula = peacefactor ~ directlyharmed + age + farmer_dar + herder_dar +
pastvoted + hhsize_darfur + female + village, data = darfur, treatment = "directlyharmed", benchmark_covariates = "female", kd = 1:3)
# info
expect_equal(darfur_out$info$formula, peacefactor ~ directlyharmed + age + farmer_dar + herder_dar +
pastvoted + hhsize_darfur + female + village)
expect_equal(darfur_out$info$treatment, "directlyharmed")
expect_equal(darfur_out$info$q, 1)
expect_equal(darfur_out$info$alpha, 0.05)
expect_equal(darfur_out$info$reduce, TRUE)
# sensitivity stats
expect_equal(darfur_out$sensitivity_stats$dof, 783)
expect_equal(darfur_out$sensitivity_stats$r2yd.x, 0.02187, tolerance = 1e-5)
expect_equivalent(c(darfur_out$sensitivity_stats$rv_q), 0.13878, tolerance = 1e-5)
expect_equivalent(c(darfur_out$sensitivity_stats$rv_qa), 0.07626, tolerance = 1e-5)
# bounds
check_bounds <- structure(list(bound_label = c("1x female", "2x female", "3x female"),
r2dz.x = c(0.00916428667504862, 0.0183285733500972, 0.0274928600251459),
r2yz.dx = c(0.12464092303637, 0.249324064199975, 0.374050471038094),
treatment = rep("directlyharmed", 3),
adjusted_estimate = c(0.0752202712144491, 0.0529151723844518, 0.0303960234641548),
adjusted_se = c(0.0218733277437572, 0.0203500620779637, 0.0186700648170924),
adjusted_t = c(3.43890386024675, 2.60024623913809, 1.62806202131271),
adjusted_lower_CI = c(0.032282966, 0.012968035,-0.006253282),
adjusted_upper_CI = c(0.11815758, 0.09286231, 0.06704533)),
.Names = c("bound_label", "r2dz.x", "r2yz.dx", "treatment",
"adjusted_estimate", "adjusted_se", "adjusted_t",
"adjusted_lower_CI", "adjusted_upper_CI"),
row.names = c(NA, -3L), class = c("ovb_bounds", "data.frame"))
expect_equivalent(darfur_out$bounds, check_bounds)
})
test_that(desc = "testing darfur sensemakr manually",
{
model.treat <- lm(directlyharmed ~ age + farmer_dar + herder_dar +
pastvoted + hhsize_darfur + female + village, data = darfur)
darfur_out <- sensemakr(estimate = 0.09731582,
se = 0.02325654,
dof = 783,
treatment = "directlyharmed",
benchmark_covariates = "female",
r2dxj.x = partial_r2(model.treat, covariates = "female"),
r2yxj.dx = partial_r2(model, covariates = "female"),
kd = 1:3)
plot(darfur_out)
plot(darfur_out, type = "extreme")
plot(darfur_out, sensitivity.of = "t-value")
add_bound_to_contour(0.3,0.3, "test")
# info
expect_equal(darfur_out$info$formula, "Data provided manually")
expect_equal(darfur_out$info$treatment, "directlyharmed")
expect_equal(darfur_out$info$q, 1)
expect_equal(darfur_out$info$alpha, 0.05)
expect_equal(darfur_out$info$reduce, TRUE)
# sensitivity stats
expect_equal(darfur_out$sensitivity_stats$dof, 783)
expect_equal(darfur_out$sensitivity_stats$r2yd.x, 0.02187, tolerance = 1e-5)
expect_equivalent(c(darfur_out$sensitivity_stats$rv_q), 0.13878, tolerance = 1e-5)
expect_equivalent(c(darfur_out$sensitivity_stats$rv_qa), 0.07626, tolerance = 1e-5)
# bounds
check_bounds <-
structure(
list(
bound_label = c("1x female", "2x female", "3x female"),
r2dz.x = c(0.00916428667504862, 0.0183285733500972, 0.0274928600251459),
r2yz.dx = c(0.12464092303637, 0.249324064199975, 0.374050471038094),
adjusted_estimate = c(0.0752202698486415, 0.0529151689180575,
0.0303960178770157),
adjusted_se = c(0.0218733298036818, 0.0203500639944344,
0.0186700665753491),
adjusted_t = c(3.43890347394571, 2.60024582392121,
1.62806156873318),
adjusted_lower_CI = c(0.0322829603180086,
0.0129680276030601, -0.00625329133645187),
adjusted_upper_CI = c(0.118157579379274,
0.092862310233055, 0.0670453270904833)
),
row.names = c(NA, -3L),
class = "data.frame"
)
expect_equivalent(as.data.frame(darfur_out$bounds), as.data.frame(check_bounds))
})
test_that(desc = "testing darfur sensitivity stats",{
# checks RV
## RV q = 1
rv <- robustness_value(model, covariates = "directlyharmed")
expect_equivalent(c(rv), 0.138776, tolerance = 1e-5)
expect_equivalent(attributes(rv)$q, 1)
expect_equivalent(attributes(rv)$names, "directlyharmed")
## RV q = 1, alpha = 0.05
rv <- robustness_value(model, covariates = "directlyharmed", q = 1, alpha = 0.05)
expect_equivalent(c(rv), 0.07625797, tolerance = 1e-5)
expect_equivalent(attributes(rv)$q, 1)
expect_equivalent(attributes(rv)$alpha, 0.05)
expect_equivalent(attributes(rv)$names, "directlyharmed")
# checks partial R2
r2 <- partial_r2(model, covariates = "directlyharmed")
expect_equivalent(r2, 0.02187309, tolerance = 1e-5)
# checks partial f2
f2 <- partial_f2(model, covariates = "directlyharmed")
expect_equivalent(f2, 0.02236222, tolerance = 1e-5)
# sensitivity stats
sens_stats <- sensitivity_stats(model, treatment = "directlyharmed")
expect_equivalent(sens_stats$treatment, "directlyharmed")
expect_equivalent(sens_stats$estimate, 0.09731582, tolerance = 1e5)
expect_equivalent(sens_stats$se, 0.02325654, tolerance = 1e5)
expect_equivalent(sens_stats$t_statistic, 4.18445, tolerance = 1e5)
expect_equivalent(sens_stats$r2yd.x, 0.02187309, tolerance = 1e5)
expect_equivalent(sens_stats$rv_q , 0.1387764, tolerance = 1e5)
expect_equivalent(sens_stats$rv_qa , 0.07625797, tolerance = 1e5)
expect_equivalent(sens_stats$f2yd.x , 0.02236222, tolerance = 1e5)
expect_equivalent(sens_stats$dof , 783, tolerance = 1e5)
expect_equivalent(group_partial_r2(model, covariates = "directlyharmed"), partial_r2(model, covariates = "directlyharmed"))
expect_error(group_partial_r2(model))
expect_equal(group_partial_r2(model, covariates = c("directlyharmed", "female")), 0.1350435, tolerance = 1e-5)
})
test_that(desc = "testing darfur adjusted estimates",{
should_be_zero <- adjusted_estimate(model, treatment = "directlyharmed", r2yz.dx = 1, r2dz.x = partial_r2(model, covariates = "directlyharmed"))
expect_equivalent(should_be_zero, 0)
rv <- robustness_value(model, covariates = "directlyharmed")
should_be_zero <- adjusted_estimate(model, treatment = "directlyharmed", r2yz.dx = rv, r2dz.x = rv)
expect_equivalent(should_be_zero, 0)
rv <- robustness_value(model, covariates = "directlyharmed", alpha = 0.05)
thr <- qt(0.975, df = 783 - 1)
should_be_1.96 <- adjusted_t(model, treatment = "directlyharmed", r2yz.dx = rv, r2dz.x = rv)
expect_equivalent(c(should_be_1.96), thr)
should_be_estimate <- bias(model, treatment = "directlyharmed", r2yz.dx = 1, r2dz.x = partial_r2(model, covariates = "directlyharmed"))
expect_equivalent(should_be_estimate, coef(model)["directlyharmed"])
rv <- robustness_value(model, covariates = "directlyharmed")
should_be_estimate <- bias(model, treatment = "directlyharmed", r2yz.dx = rv, r2dz.x = rv)
expect_equivalent(should_be_estimate, coef(model)["directlyharmed"])
rv <- robustness_value(model, covariates = "directlyharmed", q = 0.5)
should_be_half_estimate <- bias(model, treatment = "directlyharmed", r2yz.dx = rv, r2dz.x = rv)
expect_equivalent(should_be_half_estimate, coef(model)["directlyharmed"]/2)
})
test_that(desc = "testing darfur plots",
{
# testing countour output with internal functions
## point estimate
contour_out <- ovb_contour_plot(model, treatment = "directlyharmed",
benchmark_covariates = "female",kd = 1:3)
add_bound_to_contour(model, treatment = "directlyharmed",
benchmark_covariates = "age", kd = 10)
add_bound_to_contour(model, treatment = "directlyharmed",
benchmark_covariates = "age", kd = 200, ky = 20)
adj_est <- adjusted_estimate(model, treatment = "directlyharmed",
r2dz.x = contour_out$r2dz.x[10],
r2yz.dx = contour_out$r2yz.dx[1])
expect_equivalent(adj_est, contour_out$value[10])
bounds <- ovb_bounds(model,
treatment = "directlyharmed",
benchmark_covariates = "female",
kd = 1:3)
expect_equivalent(contour_out$bounds, bounds[c(2,3,1)])
## t-value
contour_out <- ovb_contour_plot(model, treatment = "directlyharmed",
benchmark_covariates = "female",kd = 1:3,
sensitivity.of = "t-value")
add_bound_to_contour(model, treatment = "directlyharmed", benchmark_covariates = "age",
kd = 200, ky = 10, sensitivity.of = "t-value")
adj_t <- adjusted_t(model, treatment = "directlyharmed",
r2dz.x = contour_out$r2dz.x[10],
r2yz.dx = contour_out$r2yz.dx[1])
expect_equivalent(adj_t, contour_out$value[10])
bounds <- ovb_bounds(model,
treatment = "directlyharmed",
benchmark_covariates = "female",
kd = 1:3)
expect_equivalent(contour_out$bounds, bounds[c(2,3,1)])
# tests bounds numerically
expect_equivalent(bounds$adjusted_estimate, c(0.0752202712144491, 0.0529151723844518, 0.0303960234641548))
expect_equivalent(bounds$adjusted_t, c(3.43890386024675, 2.60024623913809, 1.62806202131271))
# test extreme scenario plot
extreme_out <- ovb_extreme_plot(model, treatment = "directlyharmed", kd = 1:3)
adj_est <- adjusted_estimate(model, treatment = "directlyharmed",
r2yz.dx = 1,
r2dz.x = extreme_out$scenario_r2yz.dx_1$r2dz.x[5])
expect_equivalent(adj_est, extreme_out$scenario_r2yz.dx_1$adjusted_estimate[5])
}
)
test_that("testing darfur print",
{
skip_on_cran()
darfur_out <- sensemakr(model, treatment = "directlyharmed", benchmark_covariates = "female", kd = 1:3)
darfur_out2 <- sensemakr(model, treatment = "directlyharmed")
print.sense <- "Sensitivity Analysis to Unobserved Confounding\n\nModel Formula: peacefactor ~ directlyharmed + age + farmer_dar + herder_dar + \n pastvoted + hhsize_darfur + female + village\n\nNull hypothesis: q = 1 and reduce = TRUE \n\nUnadjusted Estimates of ' directlyharmed ':\n Coef. estimate: 0.09732 \n Standard Error: 0.02326 \n t-value: 4.18445 \n\nSensitivity Statistics:\n Partial R2 of treatment with outcome: 0.02187 \n Robustness Value, q = 1 : 0.13878 \n Robustness Value, q = 1 alpha = 0.05 : 0.07626 \n\nFor more information, check summary."
compare <- capture_output(print(darfur_out))
expect_equal(compare, print.sense)
summary.sense <- "Sensitivity Analysis to Unobserved Confounding\n\nModel Formula: peacefactor ~ directlyharmed + age + farmer_dar + herder_dar + \n pastvoted + hhsize_darfur + female + village\n\nNull hypothesis: q = 1 and reduce = TRUE \n-- This means we are considering biases that reduce the absolute value of the current estimate.\n-- The null hypothesis deemed problematic is H0:tau = 0 \n\nUnadjusted Estimates of 'directlyharmed': \n Coef. estimate: 0.0973 \n Standard Error: 0.0233 \n t-value (H0:tau = 0): 4.1844 \n\nSensitivity Statistics:\n Partial R2 of treatment with outcome: 0.0219 \n Robustness Value, q = 1: 0.1388 \n Robustness Value, q = 1, alpha = 0.05: 0.0763 \n\nVerbal interpretation of sensitivity statistics:\n\n-- Partial R2 of the treatment with the outcome: an extreme confounder (orthogonal to the covariates) that explains 100% of the residual variance of the outcome, would need to explain at least 2.19% of the residual variance of the treatment to fully account for the observed estimated effect.\n\n-- Robustness Value, q = 1: unobserved confounders (orthogonal to the covariates) that explain more than 13.88% of the residual variance of both the treatment and the outcome are strong enough to bring the point estimate to 0 (a bias of 100% of the original estimate). Conversely, unobserved confounders that do not explain more than 13.88% of the residual variance of both the treatment and the outcome are not strong enough to bring the point estimate to 0.\n\n-- Robustness Value, q = 1, alpha = 0.05: unobserved confounders (orthogonal to the covariates) that explain more than 7.63% of the residual variance of both the treatment and the outcome are strong enough to bring the estimate to a range where it is no longer 'statistically different' from 0 (a bias of 100% of the original estimate), at the significance level of alpha = 0.05. Conversely, unobserved confounders that do not explain more than 7.63% of the residual variance of both the treatment and the outcome are not strong enough to bring the estimate to a range where it is no longer 'statistically different' from 0, at the significance level of alpha = 0.05.\n\nBounds on omitted variable bias:\n\n--The table below shows the maximum strength of unobserved confounders with association with the treatment and the outcome bounded by a multiple of the observed explanatory power of the chosen benchmark covariate(s).\n\n Bound Label R2dz.x R2yz.dx Treatment Adjusted Estimate Adjusted Se\n 1x female 0.0092 0.1246 directlyharmed 0.0752 0.0219\n 2x female 0.0183 0.2493 directlyharmed 0.0529 0.0204\n 3x female 0.0275 0.3741 directlyharmed 0.0304 0.0187\n Adjusted T Adjusted Lower CI Adjusted Upper CI\n 3.4389 0.0323 0.1182\n 2.6002 0.0130 0.0929\n 1.6281 -0.0063 0.0670"
compare <- capture_output(summary(darfur_out))
expect_equal(compare, summary.sense)
latex.table <- c("\\begin{table}[!h]", "\\centering", "\\begin{tabular}{lrrrrrr}",
"\\multicolumn{7}{c}{Outcome: \\textit{peacefactor}} \\\\", "\\hline \\hline ",
"Treatment: & Est. & S.E. & t-value & $R^2_{Y \\sim D |{\\bf X}}$ & $RV_{q = 1}$ & $RV_{q = 1, \\alpha = 0.05}$ \\\\ ",
"\\hline ", "\\textit{directlyharmed} & 0.097 & 0.023 & 4.184 & 2.2\\% & 13.9\\% & 7.6\\% \\\\ ",
"\\hline ", "df = 783 & & \\multicolumn{5}{r}{ \\small \\textit{Bound (1x female)}: $R^2_{Y\\sim Z| {\\bf X}, D}$ = 12.5\\%, $R^2_{D\\sim Z| {\\bf X} }$ = 0.9\\%} \\\\",
"\\end{tabular}", "\\end{table}")
expect_equal(capture.output(ovb_minimal_reporting(darfur_out)), latex.table)
latex2 <- c("\\begin{table}[!h]", "\\centering", "\\begin{tabular}{lrrrrrr}",
"\\multicolumn{7}{c}{Outcome: \\textit{peacefactor}} \\\\", "\\hline \\hline ",
"Treatment: & Est. & S.E. & t-value & $R^2_{Y \\sim D |{\\bf X}}$ & $RV_{q = 1}$ & $RV_{q = 1, \\alpha = 0.05}$ \\\\ ",
"\\hline ", "\\textit{directlyharmed} & 0.097 & 0.023 & 4.184 & 2.2\\% & 13.9\\% & 7.6\\% \\\\ ",
"\\hline ", "df = 783 & & \\multicolumn{5}{r}{ }\\end{tabular}",
"\\end{table}")
expect_equal(capture.output(ovb_minimal_reporting(darfur_out2)), latex2)
})
test_that("testing darfur different q",
{
darfur_out <- sensemakr(model, treatment = "directlyharmed", benchmark_covariates = "female", q = 2, kd = 1:3)
rvq <- darfur_out$sensitivity_stats$rv_q
rvqa <- darfur_out$sensitivity_stats$rv_qa
expect_equivalent(rvq, robustness_value(model, covariates = "directlyharmed", q = 2))
expect_equivalent(rvqa, robustness_value(model, covariates = "directlyharmed", q = 2,alpha = 0.05))
}
)
test_that("Darfur group benchmarks", {
village <- grep(pattern = "village", names(coef(model)), value = T)
sensitivity <- sensemakr(model, treatment = "directlyharmed",
benchmark_covariates = list(village = village),
kd = 0.3)
r2y <- group_partial_r2(model, covariates = village)
treat.model <- update(model, directlyharmed ~ .-directlyharmed)
r2d <- group_partial_r2(treat.model, covariates = village)
bounds.check <- ovb_partial_r2_bound(r2dxj.x = r2d, r2yxj.dx = r2y, kd = 0.3)
bounds <- sensitivity$bounds
expect_equal(bounds$r2dz.x, bounds.check$r2dz.x)
expect_equal(bounds$r2yz.dx, bounds.check$r2yz.dx)
}
)
Try the sensemakr package in your browser
Any scripts or data that you put into this service are public.
sensemakr documentation built on Oct. 8, 2021, 9:10 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
context("Testing Darfur Example")
# runs regression model
data("darfur")
model2 <- lm(peacefactor ~ directlyharmed + age + farmer_dar + herder_dar +
pastvoted + hhsize_darfur + female + village,
data = within(darfur, directlyharmed <- directlyharmed*(-1)))
model <- lm(peacefactor ~ directlyharmed + age + farmer_dar + herder_dar +
pastvoted + hhsize_darfur + female + village, data = darfur)
test_that("testing darfur data", {
expect_equal(dim(darfur), c(1276, 14))
expect_equal(colnames(darfur),
c("wouldvote",
"peacefactor",
"peace_formerenemies",
"peace_jjindiv",
"peace_jjtribes",
"gos_soldier_execute",
"directlyharmed",
"age",
"farmer_dar",
"herder_dar",
"pastvoted",
"hhsize_darfur",
"village",
"female"))
})
test_that(desc = "testing darfur sensemakr",
{
darfur_out <- sensemakr(model, treatment = "directlyharmed", benchmark_covariates = "female", kd = 1:3)
plot(darfur_out)
ovb_contour_plot(darfur_out)
ovb_extreme_plot(darfur_out)
# info
expect_equal(darfur_out$info$formula, peacefactor ~ directlyharmed + age + farmer_dar + herder_dar +
pastvoted + hhsize_darfur + female + village)
expect_equal(darfur_out$info$treatment, "directlyharmed")
expect_equal(darfur_out$info$q, 1)
expect_equal(darfur_out$info$alpha, 0.05)
expect_equal(darfur_out$info$reduce, TRUE)
# sensitivity stats
expect_equal(darfur_out$sensitivity_stats$dof, 783)
expect_equal(darfur_out$sensitivity_stats$r2yd.x, 0.02187, tolerance = 1e-5)
expect_equivalent(c(darfur_out$sensitivity_stats$rv_q), 0.13878, tolerance = 1e-5)
expect_equivalent(c(darfur_out$sensitivity_stats$rv_qa), 0.07626, tolerance = 1e-5)
# bounds
check_bounds <- structure(list(bound_label = c("1x female", "2x female", "3x female"),
r2dz.x = c(0.00916428667504862, 0.0183285733500972, 0.0274928600251459),
r2yz.dx = c(0.12464092303637, 0.249324064199975, 0.374050471038094),
treatment = rep("directlyharmed", 3),
adjusted_estimate = c(0.0752202712144491, 0.0529151723844518, 0.0303960234641548),
adjusted_se = c(0.0218733277437572, 0.0203500620779637, 0.0186700648170924),
adjusted_t = c(3.43890386024675, 2.60024623913809, 1.62806202131271),
adjusted_lower_CI = c(0.032282966, 0.012968035,-0.006253282),
adjusted_upper_CI = c(0.11815758, 0.09286231, 0.06704533)
),
.Names = c("bound_label", "r2dz.x", "r2yz.dx", "treatment",
"adjusted_estimate", "adjusted_se", "adjusted_t",
"adjusted_lower_CI", "adjusted_upper_CI"),
row.names = c(NA, -3L), class = c("ovb_bounds", "data.frame"))
expect_equivalent(darfur_out$bounds, check_bounds)
out1 <- capture.output(darfur_out)
out1 <- capture.output(summary(darfur_out))
out3 <- capture.output(ovb_minimal_reporting(darfur_out))
darfur_out2 <- sensemakr(model, treatment = "directlyharmed")
plot(darfur_out2)
ovb_contour_plot(darfur_out2)
ovb_extreme_plot(darfur_out2)
out3 <- capture.output(ovb_minimal_reporting(darfur_out2))
})
test_that(desc = "testing darfur sensemakr but negative",
{
darfur_out <- sensemakr(model2, treatment = "directlyharmed", benchmark_covariates = "female", kd = 1:3)
plot(darfur_out)
ovb_contour_plot(darfur_out)
ovb_extreme_plot(darfur_out)
# info
expect_equal(darfur_out$info$formula, peacefactor ~ directlyharmed + age + farmer_dar + herder_dar +
pastvoted + hhsize_darfur + female + village)
expect_equal(darfur_out$info$treatment, "directlyharmed")
expect_equal(darfur_out$info$q, 1)
expect_equal(darfur_out$info$alpha, 0.05)
expect_equal(darfur_out$info$reduce, TRUE)
# sensitivity stats
expect_equal(darfur_out$sensitivity_stats$dof, 783)
expect_equal(darfur_out$sensitivity_stats$r2yd.x, 0.02187, tolerance = 1e-5)
expect_equivalent(c(darfur_out$sensitivity_stats$rv_q), 0.13878, tolerance = 1e-5)
expect_equivalent(c(darfur_out$sensitivity_stats$rv_qa), 0.07626, tolerance = 1e-5)
# bounds
check_bounds <- structure(list(bound_label = c("1x female", "2x female", "3x female"),
r2dz.x = c(0.00916428667504862, 0.0183285733500972, 0.0274928600251459),
r2yz.dx = c(0.12464092303637, 0.249324064199975, 0.374050471038094),
treatment = rep("directlyharmed", 3),
adjusted_estimate = c(-0.0752202712144491, -0.0529151723844518, -0.0303960234641548),
adjusted_se = c(0.0218733277437572, 0.0203500620779637, 0.0186700648170924),
adjusted_t = c(-3.43890386024675, -2.60024623913809, -1.62806202131271),
adjusted_lower_CI = -1*c(0.11815758, 0.09286231, 0.06704533),
adjusted_upper_CI = -1*c(0.032282966, 0.012968035,-0.006253282)
),
.Names = c("bound_label", "r2dz.x", "r2yz.dx", "treatment",
"adjusted_estimate", "adjusted_se", "adjusted_t",
"adjusted_lower_CI", "adjusted_upper_CI"),
row.names = c(NA, -3L), class = c("ovb_bounds", "data.frame"))
expect_equivalent(darfur_out$bounds, check_bounds)
out1 <- capture.output(darfur_out)
out1 <- capture.output(summary(darfur_out))
out3 <- capture.output(ovb_minimal_reporting(darfur_out))
darfur_out2 <- sensemakr(model, treatment = "directlyharmed")
plot(darfur_out2)
ovb_contour_plot(darfur_out2)
ovb_extreme_plot(darfur_out2)
out3 <- capture.output(ovb_minimal_reporting(darfur_out2))
})
test_that("testing darfur manual bounds",
{
sense.out <- sensemakr(model, treatment = "directlyharmed", benchmark_covariates = "female", r2dz.x = .1)
bounds.check <- sense.out$bounds
to_check <- bounds.check$adjusted_se[1]
true_check <- adjusted_se(model, treatment = "directlyharmed", r2dz.x = .1, r2yz.dx = .1)
expect_equal(to_check, unname(true_check))
}
)
test_that(desc = "testing darfur sensemakr with formula",
{
darfur_out <- sensemakr(formula = peacefactor ~ directlyharmed + age + farmer_dar + herder_dar +
pastvoted + hhsize_darfur + female + village, data = darfur, treatment = "directlyharmed", benchmark_covariates = "female", kd = 1:3)
# info
expect_equal(darfur_out$info$formula, peacefactor ~ directlyharmed + age + farmer_dar + herder_dar +
pastvoted + hhsize_darfur + female + village)
expect_equal(darfur_out$info$treatment, "directlyharmed")
expect_equal(darfur_out$info$q, 1)
expect_equal(darfur_out$info$alpha, 0.05)
expect_equal(darfur_out$info$reduce, TRUE)
# sensitivity stats
expect_equal(darfur_out$sensitivity_stats$dof, 783)
expect_equal(darfur_out$sensitivity_stats$r2yd.x, 0.02187, tolerance = 1e-5)
expect_equivalent(c(darfur_out$sensitivity_stats$rv_q), 0.13878, tolerance = 1e-5)
expect_equivalent(c(darfur_out$sensitivity_stats$rv_qa), 0.07626, tolerance = 1e-5)
# bounds
check_bounds <- structure(list(bound_label = c("1x female", "2x female", "3x female"),
r2dz.x = c(0.00916428667504862, 0.0183285733500972, 0.0274928600251459),
r2yz.dx = c(0.12464092303637, 0.249324064199975, 0.374050471038094),
treatment = rep("directlyharmed", 3),
adjusted_estimate = c(0.0752202712144491, 0.0529151723844518, 0.0303960234641548),
adjusted_se = c(0.0218733277437572, 0.0203500620779637, 0.0186700648170924),
adjusted_t = c(3.43890386024675, 2.60024623913809, 1.62806202131271),
adjusted_lower_CI = c(0.032282966, 0.012968035,-0.006253282),
adjusted_upper_CI = c(0.11815758, 0.09286231, 0.06704533)),
.Names = c("bound_label", "r2dz.x", "r2yz.dx", "treatment",
"adjusted_estimate", "adjusted_se", "adjusted_t",
"adjusted_lower_CI", "adjusted_upper_CI"),
row.names = c(NA, -3L), class = c("ovb_bounds", "data.frame"))
expect_equivalent(darfur_out$bounds, check_bounds)
})
test_that(desc = "testing darfur sensemakr manually",
{
model.treat <- lm(directlyharmed ~ age + farmer_dar + herder_dar +
pastvoted + hhsize_darfur + female + village, data = darfur)
darfur_out <- sensemakr(estimate = 0.09731582,
se = 0.02325654,
dof = 783,
treatment = "directlyharmed",
benchmark_covariates = "female",
r2dxj.x = partial_r2(model.treat, covariates = "female"),
r2yxj.dx = partial_r2(model, covariates = "female"),
kd = 1:3)
plot(darfur_out)
plot(darfur_out, type = "extreme")
plot(darfur_out, sensitivity.of = "t-value")
add_bound_to_contour(0.3,0.3, "test")
# info
expect_equal(darfur_out$info$formula, "Data provided manually")
expect_equal(darfur_out$info$treatment, "directlyharmed")
expect_equal(darfur_out$info$q, 1)
expect_equal(darfur_out$info$alpha, 0.05)
expect_equal(darfur_out$info$reduce, TRUE)
# sensitivity stats
expect_equal(darfur_out$sensitivity_stats$dof, 783)
expect_equal(darfur_out$sensitivity_stats$r2yd.x, 0.02187, tolerance = 1e-5)
expect_equivalent(c(darfur_out$sensitivity_stats$rv_q), 0.13878, tolerance = 1e-5)
expect_equivalent(c(darfur_out$sensitivity_stats$rv_qa), 0.07626, tolerance = 1e-5)
# bounds
check_bounds <-
structure(
list(
bound_label = c("1x female", "2x female", "3x female"),
r2dz.x = c(0.00916428667504862, 0.0183285733500972, 0.0274928600251459),
r2yz.dx = c(0.12464092303637, 0.249324064199975, 0.374050471038094),
adjusted_estimate = c(0.0752202698486415, 0.0529151689180575,
0.0303960178770157),
adjusted_se = c(0.0218733298036818, 0.0203500639944344,
0.0186700665753491),
adjusted_t = c(3.43890347394571, 2.60024582392121,
1.62806156873318),
adjusted_lower_CI = c(0.0322829603180086,
0.0129680276030601, -0.00625329133645187),
adjusted_upper_CI = c(0.118157579379274,
0.092862310233055, 0.0670453270904833)
),
row.names = c(NA, -3L),
class = "data.frame"
)
expect_equivalent(as.data.frame(darfur_out$bounds), as.data.frame(check_bounds))
})
test_that(desc = "testing darfur sensitivity stats",{
# checks RV
## RV q = 1
rv <- robustness_value(model, covariates = "directlyharmed")
expect_equivalent(c(rv), 0.138776, tolerance = 1e-5)
expect_equivalent(attributes(rv)$q, 1)
expect_equivalent(attributes(rv)$names, "directlyharmed")
## RV q = 1, alpha = 0.05
rv <- robustness_value(model, covariates = "directlyharmed", q = 1, alpha = 0.05)
expect_equivalent(c(rv), 0.07625797, tolerance = 1e-5)
expect_equivalent(attributes(rv)$q, 1)
expect_equivalent(attributes(rv)$alpha, 0.05)
expect_equivalent(attributes(rv)$names, "directlyharmed")
# checks partial R2
r2 <- partial_r2(model, covariates = "directlyharmed")
expect_equivalent(r2, 0.02187309, tolerance = 1e-5)
# checks partial f2
f2 <- partial_f2(model, covariates = "directlyharmed")
expect_equivalent(f2, 0.02236222, tolerance = 1e-5)
# sensitivity stats
sens_stats <- sensitivity_stats(model, treatment = "directlyharmed")
expect_equivalent(sens_stats$treatment, "directlyharmed")
expect_equivalent(sens_stats$estimate, 0.09731582, tolerance = 1e5)
expect_equivalent(sens_stats$se, 0.02325654, tolerance = 1e5)
expect_equivalent(sens_stats$t_statistic, 4.18445, tolerance = 1e5)
expect_equivalent(sens_stats$r2yd.x, 0.02187309, tolerance = 1e5)
expect_equivalent(sens_stats$rv_q , 0.1387764, tolerance = 1e5)
expect_equivalent(sens_stats$rv_qa , 0.07625797, tolerance = 1e5)
expect_equivalent(sens_stats$f2yd.x , 0.02236222, tolerance = 1e5)
expect_equivalent(sens_stats$dof , 783, tolerance = 1e5)
expect_equivalent(group_partial_r2(model, covariates = "directlyharmed"), partial_r2(model, covariates = "directlyharmed"))
expect_error(group_partial_r2(model))
expect_equal(group_partial_r2(model, covariates = c("directlyharmed", "female")), 0.1350435, tolerance = 1e-5)
})
test_that(desc = "testing darfur adjusted estimates",{
should_be_zero <- adjusted_estimate(model, treatment = "directlyharmed", r2yz.dx = 1, r2dz.x = partial_r2(model, covariates = "directlyharmed"))
expect_equivalent(should_be_zero, 0)
rv <- robustness_value(model, covariates = "directlyharmed")
should_be_zero <- adjusted_estimate(model, treatment = "directlyharmed", r2yz.dx = rv, r2dz.x = rv)
expect_equivalent(should_be_zero, 0)
rv <- robustness_value(model, covariates = "directlyharmed", alpha = 0.05)
thr <- qt(0.975, df = 783 - 1)
should_be_1.96 <- adjusted_t(model, treatment = "directlyharmed", r2yz.dx = rv, r2dz.x = rv)
expect_equivalent(c(should_be_1.96), thr)
should_be_estimate <- bias(model, treatment = "directlyharmed", r2yz.dx = 1, r2dz.x = partial_r2(model, covariates = "directlyharmed"))
expect_equivalent(should_be_estimate, coef(model)["directlyharmed"])
rv <- robustness_value(model, covariates = "directlyharmed")
should_be_estimate <- bias(model, treatment = "directlyharmed", r2yz.dx = rv, r2dz.x = rv)
expect_equivalent(should_be_estimate, coef(model)["directlyharmed"])
rv <- robustness_value(model, covariates = "directlyharmed", q = 0.5)
should_be_half_estimate <- bias(model, treatment = "directlyharmed", r2yz.dx = rv, r2dz.x = rv)
expect_equivalent(should_be_half_estimate, coef(model)["directlyharmed"]/2)
})
test_that(desc = "testing darfur plots",
{
# testing countour output with internal functions
## point estimate
contour_out <- ovb_contour_plot(model, treatment = "directlyharmed",
benchmark_covariates = "female",kd = 1:3)
add_bound_to_contour(model, treatment = "directlyharmed",
benchmark_covariates = "age", kd = 10)
add_bound_to_contour(model, treatment = "directlyharmed",
benchmark_covariates = "age", kd = 200, ky = 20)
adj_est <- adjusted_estimate(model, treatment = "directlyharmed",
r2dz.x = contour_out$r2dz.x[10],
r2yz.dx = contour_out$r2yz.dx[1])
expect_equivalent(adj_est, contour_out$value[10])
bounds <- ovb_bounds(model,
treatment = "directlyharmed",
benchmark_covariates = "female",
kd = 1:3)
expect_equivalent(contour_out$bounds, bounds[c(2,3,1)])
## t-value
contour_out <- ovb_contour_plot(model, treatment = "directlyharmed",
benchmark_covariates = "female",kd = 1:3,
sensitivity.of = "t-value")
add_bound_to_contour(model, treatment = "directlyharmed", benchmark_covariates = "age",
kd = 200, ky = 10, sensitivity.of = "t-value")
adj_t <- adjusted_t(model, treatment = "directlyharmed",
r2dz.x = contour_out$r2dz.x[10],
r2yz.dx = contour_out$r2yz.dx[1])
expect_equivalent(adj_t, contour_out$value[10])
bounds <- ovb_bounds(model,
treatment = "directlyharmed",
benchmark_covariates = "female",
kd = 1:3)
expect_equivalent(contour_out$bounds, bounds[c(2,3,1)])
# tests bounds numerically
expect_equivalent(bounds$adjusted_estimate, c(0.0752202712144491, 0.0529151723844518, 0.0303960234641548))
expect_equivalent(bounds$adjusted_t, c(3.43890386024675, 2.60024623913809, 1.62806202131271))
# test extreme scenario plot
extreme_out <- ovb_extreme_plot(model, treatment = "directlyharmed", kd = 1:3)
adj_est <- adjusted_estimate(model, treatment = "directlyharmed",
r2yz.dx = 1,
r2dz.x = extreme_out$scenario_r2yz.dx_1$r2dz.x[5])
expect_equivalent(adj_est, extreme_out$scenario_r2yz.dx_1$adjusted_estimate[5])
}
)
test_that("testing darfur print",
{
skip_on_cran()
darfur_out <- sensemakr(model, treatment = "directlyharmed", benchmark_covariates = "female", kd = 1:3)
darfur_out2 <- sensemakr(model, treatment = "directlyharmed")
print.sense <- "Sensitivity Analysis to Unobserved Confounding\n\nModel Formula: peacefactor ~ directlyharmed + age + farmer_dar + herder_dar + \n pastvoted + hhsize_darfur + female + village\n\nNull hypothesis: q = 1 and reduce = TRUE \n\nUnadjusted Estimates of ' directlyharmed ':\n Coef. estimate: 0.09732 \n Standard Error: 0.02326 \n t-value: 4.18445 \n\nSensitivity Statistics:\n Partial R2 of treatment with outcome: 0.02187 \n Robustness Value, q = 1 : 0.13878 \n Robustness Value, q = 1 alpha = 0.05 : 0.07626 \n\nFor more information, check summary."
compare <- capture_output(print(darfur_out))
expect_equal(compare, print.sense)
summary.sense <- "Sensitivity Analysis to Unobserved Confounding\n\nModel Formula: peacefactor ~ directlyharmed + age + farmer_dar + herder_dar + \n pastvoted + hhsize_darfur + female + village\n\nNull hypothesis: q = 1 and reduce = TRUE \n-- This means we are considering biases that reduce the absolute value of the current estimate.\n-- The null hypothesis deemed problematic is H0:tau = 0 \n\nUnadjusted Estimates of 'directlyharmed': \n Coef. estimate: 0.0973 \n Standard Error: 0.0233 \n t-value (H0:tau = 0): 4.1844 \n\nSensitivity Statistics:\n Partial R2 of treatment with outcome: 0.0219 \n Robustness Value, q = 1: 0.1388 \n Robustness Value, q = 1, alpha = 0.05: 0.0763 \n\nVerbal interpretation of sensitivity statistics:\n\n-- Partial R2 of the treatment with the outcome: an extreme confounder (orthogonal to the covariates) that explains 100% of the residual variance of the outcome, would need to explain at least 2.19% of the residual variance of the treatment to fully account for the observed estimated effect.\n\n-- Robustness Value, q = 1: unobserved confounders (orthogonal to the covariates) that explain more than 13.88% of the residual variance of both the treatment and the outcome are strong enough to bring the point estimate to 0 (a bias of 100% of the original estimate). Conversely, unobserved confounders that do not explain more than 13.88% of the residual variance of both the treatment and the outcome are not strong enough to bring the point estimate to 0.\n\n-- Robustness Value, q = 1, alpha = 0.05: unobserved confounders (orthogonal to the covariates) that explain more than 7.63% of the residual variance of both the treatment and the outcome are strong enough to bring the estimate to a range where it is no longer 'statistically different' from 0 (a bias of 100% of the original estimate), at the significance level of alpha = 0.05. Conversely, unobserved confounders that do not explain more than 7.63% of the residual variance of both the treatment and the outcome are not strong enough to bring the estimate to a range where it is no longer 'statistically different' from 0, at the significance level of alpha = 0.05.\n\nBounds on omitted variable bias:\n\n--The table below shows the maximum strength of unobserved confounders with association with the treatment and the outcome bounded by a multiple of the observed explanatory power of the chosen benchmark covariate(s).\n\n Bound Label R2dz.x R2yz.dx Treatment Adjusted Estimate Adjusted Se\n 1x female 0.0092 0.1246 directlyharmed 0.0752 0.0219\n 2x female 0.0183 0.2493 directlyharmed 0.0529 0.0204\n 3x female 0.0275 0.3741 directlyharmed 0.0304 0.0187\n Adjusted T Adjusted Lower CI Adjusted Upper CI\n 3.4389 0.0323 0.1182\n 2.6002 0.0130 0.0929\n 1.6281 -0.0063 0.0670"
compare <- capture_output(summary(darfur_out))
expect_equal(compare, summary.sense)
latex.table <- c("\\begin{table}[!h]", "\\centering", "\\begin{tabular}{lrrrrrr}",
"\\multicolumn{7}{c}{Outcome: \\textit{peacefactor}} \\\\", "\\hline \\hline ",
"Treatment: & Est. & S.E. & t-value & $R^2_{Y \\sim D |{\\bf X}}$ & $RV_{q = 1}$ & $RV_{q = 1, \\alpha = 0.05}$ \\\\ ",
"\\hline ", "\\textit{directlyharmed} & 0.097 & 0.023 & 4.184 & 2.2\\% & 13.9\\% & 7.6\\% \\\\ ",
"\\hline ", "df = 783 & & \\multicolumn{5}{r}{ \\small \\textit{Bound (1x female)}: $R^2_{Y\\sim Z| {\\bf X}, D}$ = 12.5\\%, $R^2_{D\\sim Z| {\\bf X} }$ = 0.9\\%} \\\\",
"\\end{tabular}", "\\end{table}")
expect_equal(capture.output(ovb_minimal_reporting(darfur_out)), latex.table)
latex2 <- c("\\begin{table}[!h]", "\\centering", "\\begin{tabular}{lrrrrrr}",
"\\multicolumn{7}{c}{Outcome: \\textit{peacefactor}} \\\\", "\\hline \\hline ",
"Treatment: & Est. & S.E. & t-value & $R^2_{Y \\sim D |{\\bf X}}$ & $RV_{q = 1}$ & $RV_{q = 1, \\alpha = 0.05}$ \\\\ ",
"\\hline ", "\\textit{directlyharmed} & 0.097 & 0.023 & 4.184 & 2.2\\% & 13.9\\% & 7.6\\% \\\\ ",
"\\hline ", "df = 783 & & \\multicolumn{5}{r}{ }\\end{tabular}",
"\\end{table}")
expect_equal(capture.output(ovb_minimal_reporting(darfur_out2)), latex2)
})
test_that("testing darfur different q",
{
darfur_out <- sensemakr(model, treatment = "directlyharmed", benchmark_covariates = "female", q = 2, kd = 1:3)
rvq <- darfur_out$sensitivity_stats$rv_q
rvqa <- darfur_out$sensitivity_stats$rv_qa
expect_equivalent(rvq, robustness_value(model, covariates = "directlyharmed", q = 2))
expect_equivalent(rvqa, robustness_value(model, covariates = "directlyharmed", q = 2,alpha = 0.05))
}
)
test_that("Darfur group benchmarks", {
village <- grep(pattern = "village", names(coef(model)), value = T)
sensitivity <- sensemakr(model, treatment = "directlyharmed",
benchmark_covariates = list(village = village),
kd = 0.3)
r2y <- group_partial_r2(model, covariates = village)
treat.model <- update(model, directlyharmed ~ .-directlyharmed)
r2d <- group_partial_r2(treat.model, covariates = village)
bounds.check <- ovb_partial_r2_bound(r2dxj.x = r2d, r2yxj.dx = r2y, kd = 0.3)
bounds <- sensitivity$bounds
expect_equal(bounds$r2dz.x, bounds.check$r2dz.x)
expect_equal(bounds$r2yz.dx, bounds.check$r2yz.dx)
}
)
Try the sensemakr package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.