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

## ----load-package-------------------------------------------------------------
library(FPScausal)

## ----sim-scalar---------------------------------------------------------------
set.seed(42)
dat <- simulate_fps_data(
  n                      = 200,
  setting                = "LL",
  outcome_type           = "scalar",
  include_functional_cov = FALSE,
  seed                   = 42
)

cat("Treatment X:", nrow(dat$X), "x", ncol(dat$X), "\n")
cat("Outcome Y:   length", length(dat$Y), "\n")
cat("Scalar C:   ", nrow(dat$C), "x", ncol(dat$C), "\n")

## ----plot-true-mu, fig.cap="True causal effect function."---------------------
plot(dat$t_grid, dat$true_beta, type = "l", lwd = 2, col = "black",
     xlab = "s", ylab = expression(mu(s)), main = "True causal effect")
abline(h = 0, lty = 2, col = "grey")

## ----weight-scalar------------------------------------------------------------
w_obj <- fps_weighting(
  treatment  = dat$X,
  treat_grid = dat$t_grid,
  domain_name = "s",
  pve        = 0.95,
  covariates = dat$C
)
print(w_obj)

## ----plot-weights-scalar, fig.cap="Distribution of FPS weights."--------------
plot(w_obj, type = "weights")

## ----plot-fpca-treatment, fig.cap="Treatment FPCA.", fig.width=9--------------
plot(w_obj, type = "fpca_treatment")

## ----plot-balance-scalar, fig.cap="Covariate balance for scalar outcome."-----
plot(w_obj, type = "balance")

## ----effect-scalar------------------------------------------------------------
eff <- fps_effect_estimation(
  outcome    = dat$Y,
  fps_object = w_obj,
  true_beta  = dat$true_beta
)
print(eff)

## ----plot-comparison-scalar, fig.cap="Weighted vs unweighted causal effect."----
plot(eff, type = "comparison")

## ----effect-scalar-boot-------------------------------------------------------
eff_boot <- fps_effect_estimation(
  outcome    = dat$Y,
  fps_object = w_obj,
  bootstrap  = TRUE,
  B          = 200,
  alpha      = 0.05,
  true_beta  = dat$true_beta,
  seed       = 123
)

## ----plot-effect-scalar, fig.cap="Causal effect with 95% bootstrap CI."-------
plot(eff_boot, type = "effect")

## ----plot-sig-scalar, fig.cap="Significant regions at alpha = 0.05."----------
plot(eff_boot, type = "significance")

## ----effect-binary------------------------------------------------------------
Y_bin   <- as.integer(dat$Y > median(dat$Y))
eff_bin <- fps_effect_estimation(Y_bin, w_obj)
print(eff_bin)

## ----plot-comparison-binary, fig.cap="Binary outcome: weighted vs unweighted."----
plot(eff_bin, type = "comparison")

## ----sim-scalar-funcov--------------------------------------------------------
dat_fc <- simulate_fps_data(
  n                      = 2000,
  setting                = "LL",
  outcome_type           = "scalar",
  include_functional_cov = TRUE,
  seed                   = 7
)

## ----weight-scalar-funcov-----------------------------------------------------
w_fc <- fps_weighting(
  treatment   = dat_fc$X,
  treat_grid  = dat_fc$t_grid,
  domain_name = "s",
  pve         = 0.95,
  covariates  = list(
    scalar     = dat_fc$C,
    functional = list(dat_fc$D)
  ),
  cov_grids   = list(dat_fc$t_grid),
  cov_pve     = 0.95
)
print(w_fc)

## ----plot-balance-funcov, fig.cap="Balance with functional covariate."--------
plot(w_fc, type = "balance")

## ----effect-scalar-funcov-----------------------------------------------------
eff_fc <- fps_effect_estimation(
  outcome    = dat_fc$Y,
  fps_object = w_fc,
  true_beta  = dat_fc$true_beta
)
plot(eff_fc, type = "comparison")

## ----all-settings, results='asis'---------------------------------------------
settings    <- c("LL", "LN", "NL", "NN")
results_tbl <- lapply(settings, function(s) {
  d   <- simulate_fps_data(200, setting = s, outcome_type = "scalar",
                            include_functional_cov = FALSE, seed = 1)
  w   <- fps_weighting(d$X, treat_grid = d$t_grid,
                        covariates = d$C)
  eff <- fps_effect_estimation(d$Y, w, true_beta = d$true_beta)
  data.frame(
    Setting        = s,
    ISE_weighted   = round(mean((eff$beta - d$true_beta)^2),  4),
    ISE_unweighted = round(mean((eff$beta_unweighted - d$true_beta)^2), 4),
    ISB_weighted   = round(mean(eff$beta - d$true_beta)^2,   6),
    ISB_unweighted = round(mean(eff$beta_unweighted - d$true_beta)^2, 6)
  )
})

knitr::kable(
  do.call(rbind, results_tbl),
  caption = "ISE and ISB for weighted vs unweighted estimate across settings"
)

## ----sim-functional-----------------------------------------------------------
dat_fn <- simulate_fps_data(
  n                      = 200,
  setting                = "LL",
  outcome_type           = "functional",
  include_functional_cov = FALSE,
  seed                   = 99
)

cat("Treatment X:", nrow(dat_fn$X), "x", ncol(dat_fn$X), "\n")
cat("Outcome Y:  ", nrow(dat_fn$Y), "x", ncol(dat_fn$Y), "\n")

## ----plot-true-surface, fig.cap="True causal effect surface mu(s,t)."---------
image(dat_fn$t_grid, dat_fn$t_grid, dat_fn$true_beta,
      xlab = "s (treatment)", ylab = "t (outcome)",
      main = expression(paste("True  ", mu, "(s,t)")),
      col  = hcl.colors(50, "Blue-Red 3"))

## ----weight-functional--------------------------------------------------------
w_fn <- fps_weighting(
  treatment   = dat_fn$X,
  treat_grid  = dat_fn$t_grid,
  treat_domain = c(0, 1),
  domain_name = "s",
  pve         = 0.95,
  covariates  = dat_fn$C
)
print(w_fn)

## ----plot-balance-functional, fig.cap="Covariate balance for functional outcome."----
plot(w_fn, type = "balance")

## ----effect-functional--------------------------------------------------------
eff_fn <- fps_effect_estimation(
  outcome             = dat_fn$Y,
  fps_object          = w_fn,
  outcome_t_grid      = dat_fn$t_grid,
  outcome_domain      = c(0, 1),
  outcome_domain_name = "t",
  outcome_pve         = 0.95,
  true_beta           = dat_fn$true_beta
)
print(eff_fn)

## ----plot-fpca-outcome, fig.cap="Outcome FPCA.", fig.width=9------------------
plot(eff_fn, type = "fpca_outcome")

## ----plot-surface, fig.cap="Estimated causal effect surface."-----------------
plot(eff_fn, type = "effect")

## ----plot-comparison-fn, fig.cap="Weighted vs unweighted surface.", fig.width=9----
plot(eff_fn, type = "comparison")

## ----effect-functional-boot---------------------------------------------------
eff_fn_boot <- fps_effect_estimation(
  outcome             = dat_fn$Y,
  fps_object          = w_fn,
  outcome_t_grid      = dat_fn$t_grid,
  outcome_domain      = c(0, 1),
  outcome_domain_name = "t",
  outcome_pve         = 0.95,
  bootstrap           = TRUE,
  B                   = 200,
  alpha               = 0.05,
  true_beta           = dat_fn$true_beta,
  seed                = 42
)

## ----plot-slice-outcome-------------------------------------------------------
plot(eff_fn_boot, type = "bootstrap_slice",
     point = 0.5, which_domain = "outcome")

## ----plot-slice-treatment-----------------------------------------------------
plot(eff_fn_boot, type = "bootstrap_slice",
     point = 0.5, which_domain = "treatment")

## ----plot-sig-fn--------------------------------------------------------------
plot(eff_fn_boot, type = "significance")

## ----all-settings-fn, results='asis'------------------------------------------
settings   <- c("LL", "LN", "NL", "NN")
results_fn <- lapply(settings, function(s) {
  d   <- simulate_fps_data(200, setting = s, outcome_type = "functional",
                            include_functional_cov = FALSE, seed = 2)
  w   <- fps_weighting(d$X, treat_grid = d$t_grid,
                        domain_name = "s",
                        covariates  = d$C)
  eff <- fps_effect_estimation(d$Y, w,
                                outcome_t_grid      = d$t_grid,
                                outcome_domain      = c(0, 1),
                                outcome_domain_name = "t",
                                true_beta           = d$true_beta)
  data.frame(
    Setting        = s,
    ISE_weighted   = round(mean((eff$beta - d$true_beta)^2),  4),
    ISE_unweighted = round(mean((eff$beta_unweighted - d$true_beta)^2), 4),
    ISB_weighted   = round(mean(eff$beta - d$true_beta)^2,   6),
    ISB_unweighted = round(mean(eff$beta_unweighted - d$true_beta)^2, 6)
  )
})

knitr::kable(
  do.call(rbind, results_fn),
  caption = "Surface ISE and ISB for weighted vs unweighted estimate"
)

## ----session-info-------------------------------------------------------------
sessionInfo()