The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

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

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.