FPScausal: Functional propensity score weighting for causal inference with functional treatments, covariates, and outcomes

In FPScausal: Functional Propensity Score for Causal Inference

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.

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 Linear or Nonlinear:

| 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.

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)$$

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:

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:

plot(w_obj, type = "weights")

Treatment FPCA: scree and eigenfunctions:

plot(w_obj, type = "fpca_treatment")

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

plot(w_obj, type = "balance")

Effect estimation (analytical CI)

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):

plot(eff, type = "comparison")

Effect estimation with bootstrap CI

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:

plot(eff_boot, type = "effect")

Significant time points (CI excludes 0):

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.

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):

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).

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

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:

plot(w_fc, type = "balance")

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.

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.

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)$$

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

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(w_fn, type = "balance")

Effect estimation (no bootstrap)

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:

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().

plot(eff_fn, type = "effect")

Weighted vs unweighted comparison:

plot(eff_fn, type = "comparison")

Effect estimation with bootstrap

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:

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

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

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

Significance map:

plot(eff_fn_boot, type = "significance")

All four simulation settings

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()

FPScausal documentation built on Aug. 9, 2026, 9:07 a.m.