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

Nicole Fontana, Simone Ciardulli

2026-08-05

1 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

estimated via weighted least squares.

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

library(FPScausal)

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


3 Part 1: Scalar continuous outcome

3.1 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 3

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

True causal effect function.

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

3.2.1 Diagnostic plots

Weight distribution:

plot(w_obj, type = "weights")

Distribution of FPS weights.

Treatment FPCA: scree and eigenfunctions:

plot(w_obj, type = "fpca_treatment")

Treatment FPCA.

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

plot(w_obj, type = "balance")

Covariate balance for scalar outcome.

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

plot(eff, type = "comparison")

Weighted vs unweighted causal effect.

3.4 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")

Causal effect with 95% bootstrap CI.

Significant time points (CI excludes 0):

plot(eff_boot, type = "significance")

Significant regions at alpha = 0.05.

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

plot(eff_bin, type = "comparison")

Binary outcome: weighted vs unweighted.


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

plot(w_fc, type = "balance")

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


3.7 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"
)
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

4 Part 2: Functional outcome

4.1 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 51

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

True causal effect surface mu(s,t).

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

Covariate balance for functional outcome.

4.3 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.4% variance explained)
#> Bootstrap CIs    : no
#> Beta surface range: [-3.7939, 3.9133]

Outcome FPCA:

plot(eff_fn, type = "fpca_outcome")

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

plot(eff_fn, type = "effect")

Estimated causal effect surface.

Weighted vs unweighted comparison:

plot(eff_fn, type = "comparison")

Weighted vs unweighted surface.

4.4 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")

4.5 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"
)
Surface ISE and ISB for weighted vs unweighted estimate
Setting ISE_weighted ISE_unweighted ISB_weighted ISB_unweighted
LL 1.9529 1.8064 0.000155 0.000140
LN 2.1111 2.0170 0.000005 0.000077
NL 1.6413 1.9482 0.000024 0.000131
NN 2.0524 2.1111 0.000092 0.000072

5 Session info

sessionInfo()
#> R version 4.3.0 (2023-04-21)
#> Platform: aarch64-apple-darwin20 (64-bit)
#> Running under: macOS Ventura 13.6
#> 
#> Matrix products: default
#> BLAS:   /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib 
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0
#> 
#> locale:
#> [1] C/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#> 
#> time zone: Europe/Rome
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] FPScausal_0.1.0
#> 
#> loaded via a namespace (and not attached):
#>  [1] ks_1.14.1          tidyr_1.3.0        sass_0.4.8         utf8_1.2.4        
#>  [5] generics_0.1.3     bitops_1.0-7       KernSmooth_2.23-20 lattice_0.21-8    
#>  [9] wCorr_1.9.8        hdrcde_3.4         hms_1.1.3          pracma_2.4.4      
#> [13] digest_0.6.33      magrittr_2.0.3     evaluate_0.23      grid_4.3.0        
#> [17] mvtnorm_1.2-4      fastmap_1.1.1      jsonlite_1.8.8     Matrix_1.6-4      
#> [21] progress_1.2.3     deSolve_1.40       mclust_6.0.1       purrr_1.0.2       
#> [25] fansi_1.0.6        scales_1.3.0       isoband_0.2.7      jquerylib_0.1.4   
#> [29] mnormt_2.1.1       cli_3.6.2          crayon_1.5.3       rlang_1.1.6       
#> [33] fda_6.2.0          munsell_0.5.0      splines_4.3.0      withr_2.5.2       
#> [37] cachem_1.0.8       yaml_2.3.8         tools_4.3.0        minqa_1.2.6       
#> [41] dplyr_1.1.4        colorspace_2.1-0   ggplot2_3.5.2      rainbow_3.7       
#> [45] vctrs_0.6.5        R6_2.5.1           lifecycle_1.0.4    fds_1.8           
#> [49] MASS_7.3-58.4      pcaPP_2.0-4        cluster_2.1.4      pkgconfig_2.0.3   
#> [53] bslib_0.6.1        pillar_1.9.0       gtable_0.3.6       Rcpp_1.1.0        
#> [57] glue_1.6.2         xfun_0.41          tibble_3.2.1       tidyselect_1.2.0  
#> [61] highr_0.10         rstudioapi_0.15.0  knitr_1.45         farver_2.1.1      
#> [65] patchwork_1.2.0    htmltools_0.5.7    rmarkdown_2.25     labeling_0.4.3    
#> [69] compiler_4.3.0     prettyunits_1.2.0  RCurl_1.98-1.13