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 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)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.
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.
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.0050Weight 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.
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.
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.
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.
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")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 |
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).
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.0050plot(w_fn, type = "balance")Covariate balance for functional outcome.
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.
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")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.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 |
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.13These 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.