FPScausal: Functional propensity score weighting for causal inference with functional treatments, covariates, and outcomes
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.
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")
#> Treatment X: 200 x 51
cat("Outcome Y: length", length(dat$Y), "\n")
#> Outcome Y: length 200
cat("Scalar C: ", nrow(dat$C), "x", ncol(dat$C), "\n")
#> Scalar C: 200 x 3The 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")
True causal effect function.
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)
#> Functional Propensity Score Weighting
#> --------------------------------------
#> Observations : 200
#> Treatment FPCs (L) : 4 (96.9% variance explained)
#> Convergence code : 0 (converged)
#> Weights : min = 0.0000 max = 0.0315 mean = 0.0050Diagnostic plots
Weight distribution:
Distribution of FPS weights.
Treatment FPCA: scree and eigenfunctions:
Treatment FPCA.
Covariate balance: absolute Pearson correlations before (red) and after (blue) weighting. Dashed line at 0.1:
Covariate balance for scalar outcome.
Effect estimation (analytical CI)
eff <- fps_effect_estimation(
outcome = dat$Y,
fps_object = w_obj,
true_beta = dat$true_beta
)
print(eff)
#> Functional Propensity Score Effect Estimation
#> ----------------------------------------------
#> Outcome type : scalar
#> Treatment FPCs : 4 (96.9% variance explained)
#> Bootstrap CIs : no
#> Beta range : [-3.4264, 3.8265]Weighted vs unweighted comparison (with analytical CI):
Weighted vs unweighted causal effect.
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:
Causal effect with 95% bootstrap CI.
Significant time points (CI excludes 0):
Significant regions at alpha = 0.05.
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)
#> Functional Propensity Score Effect Estimation
#> ----------------------------------------------
#> Outcome type : binary
#> Treatment FPCs : 4 (96.9% variance explained)
#> Bootstrap CIs : no
#> Beta range : [-0.1644, 0.1944]Weighted vs unweighted comparison (with analytical CI):
Binary outcome: weighted vs unweighted.
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)
#> Functional Propensity Score Weighting
#> --------------------------------------
#> Observations : 2000
#> Treatment FPCs (L) : 4 (96.5% variance explained)
#> Convergence code : 0 (converged)
#> Weights : min = 0.0000 max = 0.0288 mean = 0.0005
#> Functional covariates : 1 (FPCs retained: 4)The FPC scores of the functional covariate are automatically named
Func_Cov1_FPC1, Func_Cov1_FPC2, … in the
balance plot:
Balance with functional covariate.
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"
)| Setting | ISE_weighted | ISE_unweighted | ISB_weighted | ISB_unweighted |
|---|---|---|---|---|
| LL | 0.1431 | 0.4140 | 0e+00 | 1e-06 |
| LN | 0.1471 | 0.3976 | 0e+00 | 3e-06 |
| NL | 0.0666 | 0.3082 | 2e-06 | 5e-06 |
| NN | 0.0702 | 0.2902 | 2e-06 | 3e-06 |
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")
#> Treatment X: 200 x 51
cat("Outcome Y: ", nrow(dat_fn$Y), "x", ncol(dat_fn$Y), "\n")
#> Outcome Y: 200 x 51The 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"))
True causal effect surface mu(s,t).
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)
#> Functional Propensity Score Weighting
#> --------------------------------------
#> Observations : 200
#> Treatment FPCs (L) : 4 (96.2% variance explained)
#> Convergence code : 0 (converged)
#> Weights : min = 0.0000 max = 0.0662 mean = 0.0050
Covariate balance for functional outcome.
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)
#> Functional Propensity Score Effect Estimation
#> ----------------------------------------------
#> Outcome type : functional
#> Treatment FPCs : 4 (96.2% variance explained)
#> Outcome FPCs : 3 (95.5% variance explained)
#> Bootstrap CIs : no
#> Beta surface range: [-3.7391, 3.9210]Outcome FPCA:
Outcome FPCA.
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().
Estimated causal effect surface.
Weighted vs unweighted comparison:
Weighted vs unweighted surface.
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:
1-D slice — fixing exposure time s = 0.5:
Significance map:
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"
)| Setting | ISE_weighted | ISE_unweighted | ISB_weighted | ISB_unweighted |
|---|---|---|---|---|
| LL | 1.9535 | 1.8060 | 0.000158 | 0.000142 |
| LN | 2.1097 | 2.0167 | 0.000004 | 0.000078 |
| NL | 1.6394 | 1.9476 | 0.000026 | 0.000135 |
| NN | 2.0504 | 2.1107 | 0.000090 | 0.000074 |
Session info
sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#>
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so; LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] FPScausal_0.1.1 rmarkdown_2.32
#>
#> loaded via a namespace (and not attached):
#> [1] gtable_0.3.6 xfun_0.61 bslib_0.12.0 ggplot2_4.0.3
#> [5] ks_1.15.3 lattice_0.23-1 vctrs_0.7.3 tools_4.6.1
#> [9] bitops_1.1-0 generics_0.1.4 tibble_3.3.1 cluster_2.1.8.3
#> [13] pkgconfig_2.0.3 fds_1.9 Matrix_1.7-6 KernSmooth_2.23-27
#> [17] RColorBrewer_1.1-3 S7_0.2.2 lifecycle_1.0.5 compiler_4.6.1
#> [21] farver_2.1.2 progress_1.2.3 mnormt_2.1.2 htmltools_0.5.9
#> [25] sys_3.4.3 buildtools_1.0.0 rainbow_3.8 sass_0.4.10
#> [29] RCurl_1.98-1.20 yaml_2.3.12 pracma_2.4.6 crayon_1.5.3
#> [33] pillar_1.11.1 jquerylib_0.1.4 tidyr_1.3.2 MASS_7.3-66
#> [37] cachem_1.1.0 mclust_6.1.3 tidyselect_1.2.1 digest_0.6.39
#> [41] mvtnorm_1.4-2 dplyr_1.2.1 purrr_1.2.2 maketools_1.3.2
#> [45] labeling_0.4.3 splines_4.6.1 pcaPP_2.0-5 fastmap_1.2.0
#> [49] grid_4.6.1 colorspace_2.1-3 cli_3.6.6 magrittr_2.0.5
#> [53] patchwork_1.3.2 withr_3.0.3 prettyunits_1.2.0 scales_1.4.0
#> [57] otel_0.2.0 deSolve_1.42 hms_1.1.4 evaluate_1.0.5
#> [61] knitr_1.52 rlang_1.3.0 isoband_0.3.0 Rcpp_1.1.2
#> [65] glue_1.8.1 fda_6.3.0 hdrcde_3.5.0 wCorr_1.9.8
#> [69] minqa_1.2.8 jsonlite_2.0.0 R6_2.6.1