## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup--------------------------------------------------------------------
library(causaldef)
library(ggplot2)

# DAG Helper
plot_dag <- function(coords, edges, title = NULL) {
  edges_df <- merge(edges, coords, by.x = "from", by.y = "name")
  colnames(edges_df)[c(3,4)] <- c("x_start", "y_start")
  edges_df <- merge(edges_df, coords, by.x = "to", by.y = "name")
  colnames(edges_df)[c(5,6)] <- c("x_end", "y_end")
  
  ggplot2::ggplot(coords, ggplot2::aes(x = x, y = y)) +
    ggplot2::geom_segment(data = edges_df, ggplot2::aes(x = x_start, y = y_start, xend = x_end, yend = y_end),
                 arrow = ggplot2::arrow(length = ggplot2::unit(0.3, "cm"), type = "closed"), 
                 color = "gray40", size = 1, alpha = 0.8) +
    ggplot2::geom_point(size = 14, color = "white", fill = "#E7B800", shape = 21, stroke = 1.5) + # different color
    ggplot2::geom_text(ggplot2::aes(label = name), fontface = "bold", size = 3, color = "black") +
    ggplot2::ggtitle(title) + 
    ggplot2::theme_void(base_size = 14) +
    ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5, face = "bold", margin = ggplot2::margin(b = 10))) +
    ggplot2::coord_fixed()
}

## ----load_data----------------------------------------------------------------
data("gene_perturbation")
head(gene_perturbation)

## ----nc_spec------------------------------------------------------------------
spec_nc <- causal_spec(
  data = gene_perturbation,
  treatment = "knockout_status",
  outcome = "target_expression",
  covariates = c("library_size", "batch"),
  negative_control = "housekeeping_gene"
)

# Visualize the Structural Assumption
coords <- data.frame(
  name = c("Unobserved", "Treatment", "Target", "Housekp"),
  x = c(0, -1.5, 1.5, 0),
  y = c(1, 0, 0, -1)
)
edges <- data.frame(
  from = c("Unobserved", "Unobserved", "Unobserved", "Treatment"),
  to =   c("Treatment", "Target", "Housekp", "Target")
)
plot_dag(coords, edges, title = "Why Negative Controls Work:\nShared Unobserved Confounding")

## ----nc_run-------------------------------------------------------------------
set.seed(123)
# Check if IPTW approach effectively removes confounding using the negative control
nc_res <- nc_diagnostic(spec_nc, method = "iptw", n_boot = 100)
print(nc_res)

## ----sensitivity_setup--------------------------------------------------------
data("hct_outcomes")

# Convert treatment to numeric for the Gaussian model approximation
hct_outcomes$conditioning_numeric <- as.numeric(hct_outcomes$conditioning_intensity) - 1

# Simplified specification focusing on conditioning intensity
spec_sens <- causal_spec(
  data = hct_outcomes,
  treatment = "conditioning_numeric",
  outcome = "time_to_event", # Simplified for this example
  covariates = c("age", "disease_status")
)

## ----run_frontier-------------------------------------------------------------
frontier <- confounding_frontier(
  spec_sens,
  alpha_range = c(-1, 1),
  gamma_range = c(-1, 1),
  grid_size = 40
)

# The result contains a grid we can visualize
head(frontier$grid)

## ----plot_frontier, fig.width=6, fig.height=5---------------------------------
ggplot(frontier$grid, aes(x = alpha, y = gamma, fill = delta)) +
  geom_tile() +
  scale_fill_viridis_c(option = "magma", direction = -1) +
  labs(
    title = "Confounding Frontier",
    x = "Confounding Strength U -> A (alpha)",
    y = "Confounding Strength U -> Y (gamma)",
    fill = "Deficiency\n(Delta)"
  ) +
  theme_minimal() +
  geom_contour(aes(z = delta), color = "white", breaks = c(0.01, 0.05, 0.1))