# 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

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

# 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