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

## ----setup--------------------------------------------------------------------
library(causaldef)
library(ggplot2) # Ensure ggplot2 is loaded

# DAG Helper (using ggplot2 only)
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 = "#8FBC8F", shape = 21, stroke = 1.5) + # DarkSeaGreen
    ggplot2::geom_text(ggplot2::aes(label = name), fontface = "bold", size = 2.5, color = "white") +
    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()
}

data(gene_perturbation)

head(gene_perturbation)

# Visualize the Structural Assumption
coords <- data.frame(
  name = c("Unobserved", "Knockout", "TargetExpr", "Housekp"),
  x = c(0, -1.5, 1.5, 0),
  y = c(1, 0, 0, -1)
)
edges <- data.frame(
  from = c("Unobserved", "Unobserved", "Unobserved", "Knockout"),
  to =   c("Knockout", "TargetExpr", "Housekp", "TargetExpr")
)
plot_dag(coords, edges, title = "Gene Perturbation DAG")

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

print(spec)

## ----deficiency---------------------------------------------------------------
# Compare unadjusted vs IPTW
results <- estimate_deficiency(
  spec,
  methods = c("unadjusted", "iptw"),
  n_boot = 50  # Low for vignette speed
)

print(results)
plot(results, type = "bar")

## ----diagnosis----------------------------------------------------------------
nc_test <- nc_diagnostic(spec, method = "iptw")
print(nc_test)

## ----decision-----------------------------------------------------------------
# Utility scale: 0 to 10 (e.g., normalized expression fold-change value)
bounds <- policy_regret_bound(results, utility_range = c(0, 10), method = "iptw")

print(bounds)

## ----effect-------------------------------------------------------------------
# Estimate ATE using the best performing method (IPTW)
effect <- estimate_effect(results, target_method = "iptw")
print(effect)