---
title: "FPScausal: Functional propensity score weighting for causal inference with functional treatments, covariates, and outcomes"
author: "Nicole Fontana, Simone Ciardulli"
date: "`r Sys.Date()`"
output:
  rmarkdown::html_vignette:
    toc: true
    toc_depth: 3
    number_sections: true
vignette: >
  %\VignetteIndexEntry{FPScausal: Functional propensity score weighting for causal inference with functional treatments, covariates, and outcomes}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# Introduction

**FPScausal** implements the *Functional Propensity Score (FPS) weighting*
methodology for causal inference with functional treatments
(Ciardulli, S. and Fontana, N., 2026).

The core idea is to represent a functional treatment $X(s)$ through its
Functional Principal Component (FPC) scores $\mathbf{A} \in \mathbb{R}^L$,
and then to estimate covariate-balancing weights $\{w_i\}$ by maximising the
empirical likelihood subject to the balancing constraints
$$\frac{1}{n}\sum_{i=1}^n w_i \mathbf{g}_i = \mathbf{0}, \quad
  \sum_{i=1}^n w_i = 1,$$
where $\mathbf{g}_i = [\mathbf{A}_i^\top, \mathbf{C}_i^\top,
\mathrm{vec}(\mathbf{A}_i\mathbf{C}_i^\top)^\top]^\top$ stacks the
balancing moments for unit $i$.  Here $\mathbf{C}_i$ denotes the vector
of **confounders** observed for unit $i$ (e.g.\ demographic variables
or baseline measurements).
The dual of this empirical-likelihood problem reduces to the smooth,
unconstrained minimisation
$$\min_{\boldsymbol\theta} \log\!\Bigl(\sum_{i=1}^n
  e^{-\boldsymbol\theta^\top \mathbf{g}_i}\Bigr),$$
solved via the BFGS quasi-Newton algorithm.  The weights are recovered as
the softmax transformation
$w_i = e^{-\boldsymbol\theta^{*\top}\mathbf{g}_i}/
       \sum_j e^{-\boldsymbol\theta^{*\top}\mathbf{g}_j}$.

Once the weights are obtained, the causal estimand is

- the **causal effect function** $\mu(s)$ (scalar or binary outcome), or
- the **causal effect surface** $\mu(s,t)$ (functional outcome)

estimated via weighted least squares.

This vignette walks through the full workflow on simulated data for two
outcome types.

```{r load-package}
library(FPScausal)
```

---

# Simulation settings

`simulate_fps_data()` implements the data-generating process from the
simulation study in the paper.  The four settings ("LL", "LN", "NL", "NN")
control whether the treatment-confounder and the confounder-outcome
relationships are **L**inear or **N**onlinear:

| Setting | Treatment–Confounder | Confounder–Outcome |
|---------|---------------------|--------------------|
| LL      | Linear              | Linear             |
| LN      | Linear              | Nonlinear          |
| NL      | Nonlinear           | Linear             |
| NN      | Nonlinear           | Nonlinear          |

The treatment $X(s)$ is built from six Fourier eigenfunctions; the scalar
confounders $\mathbf{C}$ are 3-dimensional; one functional covariate $D(s)$
(4 Fourier components) is optionally included.

---

# Part 1: Scalar continuous outcome

## Data generation

We simulate $n = 200$ subjects under setting **"LL"** with **scalar
covariates only** (no functional covariate) and a scalar continuous outcome.

```{r 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")
```

The true causal effect function is:

$$\mu(s) = 2\sqrt{2}\sin(2\pi s) + \sqrt{2}\cos(2\pi s)
           + \tfrac{\sqrt{2}}{2}\sin(4\pi s)
           + \tfrac{\sqrt{2}}{2}\cos(4\pi s)$$

```{r 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 estimation

The treatment domain `treat_domain` is inferred automatically from
`treat_grid` when omitted:

```{r 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)
```

### Diagnostic plots

**Weight distribution:**

```{r plot-weights-scalar, fig.cap="Distribution of FPS weights."}
plot(w_obj, type = "weights")
```

**Treatment FPCA:** scree and eigenfunctions:

```{r plot-fpca-treatment, fig.cap="Treatment FPCA.", fig.width=9}
plot(w_obj, type = "fpca_treatment")
```

**Covariate balance:** absolute Pearson correlations before (red) and after
(blue) weighting. Dashed line at 0.1:

```{r plot-balance-scalar, fig.cap="Covariate balance for scalar outcome."}
plot(w_obj, type = "balance")
```

## Effect estimation (analytical CI)

```{r effect-scalar}
eff <- fps_effect_estimation(
  outcome    = dat$Y,
  fps_object = w_obj,
  true_beta  = dat$true_beta
)
print(eff)
```

**Weighted vs unweighted comparison (with analytical CI):**

```{r plot-comparison-scalar, fig.cap="Weighted vs unweighted causal effect."}
plot(eff, type = "comparison")
```

## Effect estimation with bootstrap CI

```{r 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
)
```

**Effect with 95% bootstrap CI:**

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

**Significant time points** (CI excludes 0):

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

## Binary outcome

When the outcome is binary (0/1), `fps_effect_estimation` automatically
detects it and fits a **linear probability model** (weighted least squares),
returning the average treatment effect on the probability scale.

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

**Weighted vs unweighted comparison (with analytical CI):**

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

---

## Scalar outcome with a functional covariate

When a functional covariate $D(s)$ is available, it enters the balancing step
through its own FPC scores.  We use $n = 2000$ to ensure a stable weight
solution (the constraint dimension grows with the number of FPCs).

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

```{r 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)
```

The FPC scores of the functional covariate are automatically named
`Func_Cov1_FPC1`, `Func_Cov1_FPC2`, ... in the balance plot:

```{r plot-balance-funcov, fig.cap="Balance with functional covariate."}
plot(w_fc, type = "balance")
```

```{r 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 four simulation settings

The table below shows the Integrated Squared Error (ISE) and Integrated
Squared Bias (ISB) of the weighted vs unweighted estimate across settings.

```{r 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"
)
```

---

# Part 2: Functional outcome

## Data generation

Now we simulate with a **functional outcome** $Y(t)$, so the causal estimand
is the bivariate effect surface $\mu(s,t)$.
We use $n = 200$ with scalar covariates only for this illustration.

```{r 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")
```

The true surface is:

$$\mu(s,t) = 2\sqrt{2}\sin(2\pi s)\cos(2\pi t)
            + 2\sqrt{2}\sin(2\pi t)\cos(2\pi s)
            + \sqrt{2}\cos(4\pi t)\sin(4\pi s)
            + \sqrt{2}\cos(4\pi s)\sin(4\pi t)$$

```{r 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 estimation

```{r 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)
```

```{r plot-balance-functional, fig.cap="Covariate balance for functional outcome."}
plot(w_fn, type = "balance")
```

## Effect estimation (no bootstrap)

```{r 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)
```

**Outcome FPCA:**

```{r plot-fpca-outcome, fig.cap="Outcome FPCA.", fig.width=9}
plot(eff_fn, type = "fpca_outcome")
```

**Estimated effect surface (weighted):**

The dashed black contour lines overlay the true surface $\mu(s,t)$ for
reference — they appear because `true_beta` was passed to
`fps_effect_estimation()`.

```{r plot-surface, fig.cap="Estimated causal effect surface."}
plot(eff_fn, type = "effect")
```

**Weighted vs unweighted comparison:**

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

## Effect estimation with bootstrap

```{r 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
)
```

**1-D slice of the effect surface — fixing outcome time t = 0.5:**

```{r plot-slice-outcome}
plot(eff_fn_boot, type = "bootstrap_slice",
     point = 0.5, which_domain = "outcome")
```

**1-D slice — fixing exposure time s = 0.5:**

```{r plot-slice-treatment}
plot(eff_fn_boot, type = "bootstrap_slice",
     point = 0.5, which_domain = "treatment")
```

**Significance map:**

```{r plot-sig-fn}
plot(eff_fn_boot, type = "significance")
```

## All four simulation settings

```{r 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

```{r session-info}
sessionInfo()
```