## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>",
    fig.width = 7,
    fig.height = 5
)

## ----example-setup------------------------------------------------------------
library(causaldef)
set.seed(42)

n <- 500
U <- rnorm(n) # Unmeasured confounder
W <- 0.7 * U + rnorm(n, sd = 0.5) # Observed covariate
A <- rbinom(n, 1, plogis(0.5 * U)) # Treatment
Y <- 2 * A + 1.5 * U + rnorm(n) # Outcome (true effect = 2)

df <- data.frame(W = W, A = A, Y = Y)

spec <- causal_spec(df, "A", "Y", "W")

## ----deficiency-calc----------------------------------------------------------
def_results <- estimate_deficiency(
    spec,
    methods = c("unadjusted", "iptw"),
    n_boot = 100
)

print(def_results)

## ----frontier-plot, fig.height=5----------------------------------------------
frontier <- confounding_frontier(
    spec,
    alpha_range = c(-3, 3),
    gamma_range = c(-3, 3),
    grid_size = 40
)

plot(frontier)

## ----policy-bound-------------------------------------------------------------
bounds <- policy_regret_bound(def_results, utility_range = c(0, 10), method = "iptw")
print(bounds)

## ----evalue-calculation-------------------------------------------------------
# Effect estimate
effect <- estimate_effect(def_results, target_method = "iptw")
print(effect)

# For E-value calculation, we need to convert to risk ratio scale
# This is an approximation; exact E-values require binary outcomes
# Here we use the standardized effect size

# Assuming effect is on continuous scale, we can compute
# a pseudo-risk ratio via effect size transformation
effect_se <- 1 # Approximate SE (would be from bootstrap in practice)
effect_est <- effect$estimate
t_stat <- effect_est / effect_se

# Approximate conversion to OR (for conceptual illustration)
# Using Chinn's (2000) formula for approximate OR from mean difference
approx_or <- exp(effect_est / 1.81)

# E-value formula
if (approx_or > 1) {
    e_value <- approx_or + sqrt(approx_or * (approx_or - 1))
} else {
    e_value <- 1 / approx_or + sqrt(1 / approx_or * (1 / approx_or - 1))
}

cat(sprintf("
Approximate E-value: %.2f

Interpretation: To explain away the observed effect, an unmeasured
confounder would need a risk ratio of at least %.2f with both
treatment and outcome.
", e_value, e_value))

## ----benchmark----------------------------------------------------------------
# Add more covariates for benchmarking
df$W2 <- 0.3 * U + rnorm(n, sd = 0.8)
df$W3 <- 0.9 * U + rnorm(n, sd = 0.3)

spec_multi <- causal_spec(df, "A", "Y", c("W", "W2", "W3"))

frontier_bench <- confounding_frontier(
    spec_multi,
    grid_size = 40
)

# Access benchmarks
if (!is.null(frontier_bench$benchmarks)) {
    print(frontier_bench$benchmarks)
}

## ----full-analysis------------------------------------------------------------
# Add negative control
df$Y_nc <- U + rnorm(n, sd = 0.5) # Affected by U, not by A

spec_full <- causal_spec(
    df, "A", "Y", c("W", "W2", "W3"),
    negative_control = "Y_nc"
)

# Complete analysis
def_full <- estimate_deficiency(
    spec_full,
    methods = c("unadjusted", "iptw"),
    n_boot = 100
)

nc_full <- nc_diagnostic(spec_full, method = "iptw", n_boot = 100)

print(def_full)
print(nc_full)