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

## ----setup--------------------------------------------------------------------
library(CausalSpline)

## ----simulate-----------------------------------------------------------------
set.seed(42)
dat <- simulate_dose_response(n = 600, dgp = "threshold", confounding = 0.6)
head(dat)

## ----true-curve, fig.cap="True vs observed relationship"----------------------
plot(dat$T, dat$Y, pch = 16, col = rgb(0, 0, 0, 0.2),
     xlab = "Treatment T", ylab = "Outcome Y",
     main = "Observed data (confounded)")
lines(sort(dat$T), dat$true_effect[order(dat$T)],
      col = "red", lwd = 2)
legend("topleft", legend = "True f(T)", col = "red", lwd = 2)

## ----fit-ipw------------------------------------------------------------------
fit_ipw <- causal_spline(
  Y ~ T | X1 + X2 + X3,
  data       = dat,
  method     = "ipw",
  df_exposure = 5,
  eval_grid  = 100
)
summary(fit_ipw)

## ----plot-ipw, fig.cap="IPW estimated dose-response with 95% CI"--------------
# Build true curve data frame for overlay
truth_df <- data.frame(
  t           = dat$T,
  true_effect = dat$true_effect
)
plot(fit_ipw, truth = truth_df)

## ----fit-gcomp----------------------------------------------------------------
fit_gc <- causal_spline(
  Y ~ T | X1 + X2 + X3,
  data        = dat,
  method      = "gcomp",
  df_exposure = 5
)
plot(fit_gc, truth = truth_df,
     title = "G-computation — Threshold DGP")

## ----overlap------------------------------------------------------------------
ov <- check_overlap(dat$T, fit_ipw$weights, plot = TRUE)
cat("ESS:", round(ov$ess), "/ n =", nrow(dat), "\n")
ov$plot

## ----compare-dgps, fig.height=8, fig.width=7----------------------------------
dgps <- c("threshold", "diminishing", "nonmonotone", "sinusoidal")
plots <- lapply(dgps, function(d) {
  dat_d <- simulate_dose_response(500, dgp = d, seed = 1)
  fit_d <- causal_spline(Y ~ T | X1 + X2 + X3, data = dat_d,
                          method = "ipw", df_exposure = 5,
                          verbose = FALSE)
  truth_d <- data.frame(t = dat_d$T, true_effect = dat_d$true_effect)
  plot(fit_d, truth = truth_d,
       title = paste("DGP:", d),
       rug = FALSE)
})

# Combine with patchwork (if available) or print individually
for (p in plots) print(p)