| Type: | Package |
| Title: | Functional Propensity Score for Causal Inference |
| Version: | 0.1.0 |
| Description: | Implements functional propensity score (FPS) weighting for causal inference with functional treatments. Weights are estimated by maximising the empirical likelihood subject to covariate-balancing constraints and solving the resulting dual problem via the BFGS quasi-Newton algorithm, following Ciardulli, S. and Fontana, N. (2026). The package supports scalar, binary, and functional outcomes, as well as functional covariates. |
| License: | MIT + file LICENSE |
| Encoding: | UTF-8 |
| Depends: | R (≥ 4.1.0) |
| Imports: | fda (≥ 6.0.0), ggplot2 (≥ 3.4.0), tidyr (≥ 1.2.0), MASS (≥ 7.3-0), wCorr, patchwork (≥ 1.1.0), progress (≥ 1.2.0), stats, utils |
| Suggests: | testthat (≥ 3.0.0), knitr, rmarkdown |
| VignetteBuilder: | knitr |
| RoxygenNote: | 7.3.1 |
| Config/testthat/edition: | 3 |
| NeedsCompilation: | no |
| Packaged: | 2026-08-05 07:31:09 UTC; nicole.fontana |
| Author: | Nicole Fontana [aut, cre], Simone Ciardulli [aut] |
| Maintainer: | Nicole Fontana <nicole.fontana@polimi.it> |
| Repository: | CRAN |
| Date/Publication: | 2026-08-09 08:00:02 UTC |
FPScausal: Causal Inference with Functional Propensity Score Weighting
Description
FPScausal implements the functional propensity score (FPS) weighting methodology for causal inference with functional treatments. Given a functional treatment X(t), a set of scalar and/or functional confounders, and an outcome (continuous, binary, or functional), the package:
1. Decomposes X(t) into Functional Principal Component (FPC) scores via
fda.
2. Estimates covariate-balancing weights by maximising the empirical
likelihood subject to balancing constraints and solving the resulting
dual problem via the BFGS quasi-Newton algorithm.
3. Estimates the causal effect function \mu(t) (scalar/binary
outcome) or causal effect surface \mu(s,t) (functional outcome)
via weighted least squares.
4. Provides pointwise bootstrap confidence intervals and a rich set of
diagnostic and result plots.
Main functions
fps_weightingEstimate FPS weights.
fps_effect_estimationEstimate causal effects using FPS weights.
simulate_fps_dataGenerate synthetic datasets for testing and vignette illustrations.
S3 methods
Both fps_weighting and fps_effect_estimation objects support
print, summary, and plot with a type argument.
References
Ciardulli, S. and Fontana, N. (2026). Functional propensity score weighting for causal inference with functional treatments, covariates, and outcomes. Preprint.
Author(s)
Maintainer: Nicole Fontana nicole.fontana@polimi.it
Authors:
Simone Ciardulli simone.ciardulli@mail.polimi.it
Auto-select number of B-spline basis functions
Description
Auto-select number of B-spline basis functions
Usage
.auto_nbasis(t_grid)
Arguments
t_grid |
Numeric vector of observation time points. |
Value
Integer: suggested nbasis.
Check for domain overlap between treatment and functional outcome
Description
Emits a warning when the treatment and outcome share the same domain name and the domains numerically overlap, because the historical constraint (non-anticipativity) may be needed but is not enforced.
Usage
.check_domains_overlap(domain_t, name_t, domain_o, name_o)
Arguments
domain_t |
Numeric c(a, b) treatment domain. |
name_t |
Character treatment domain name. |
domain_o |
Numeric c(a, b) outcome domain. |
name_o |
Character outcome domain name. |
Check that an object is of class fps_weighting
Description
Check that an object is of class fps_weighting
Usage
.check_fps_weighting(fps_object)
Compute MISE, AISE, and ISB between estimated and true beta
Description
Compute MISE, AISE, and ISB between estimated and true beta
Usage
.compute_error_metrics(beta_hat, beta_true)
Arguments
beta_hat |
Numeric vector (or matrix for functional case). |
beta_true |
Numeric vector (or matrix) of the same shape. |
Value
Named numeric vector with MISE, AISE, ISB (for single estimate, MISE=AISE=ISB=ISE).
Find contiguous intervals where a logical mask is TRUE
Description
Returns a data frame with 'start' and 'end' columns (values from t_grid) for each run of TRUE values in 'mask'.
Usage
.contiguous_intervals(t_grid, mask)
Arguments
t_grid |
Numeric vector. |
mask |
Logical vector of same length as t_grid. |
Value
data.frame with columns start, end (or 0-row frame if no TRUE).
Detect the type of an outcome variable
Description
Detect the type of an outcome variable
Usage
.detect_outcome_type(outcome)
Arguments
outcome |
A numeric vector or matrix. |
Value
One of '"scalar"', '"binary"', or '"functional"'.
Estimate beta (effect function or surface)
Description
Estimate beta (effect function or surface)
Usage
.estimate_beta(outcome_type, outcome, A, efX, w, fpca_out = NULL)
Fit weighted regression model for a given outcome type
Description
Fit weighted regression model for a given outcome type
Usage
.fit_weighted_model(outcome_type, outcome, A, w)
Internal dual formulation weight estimator
Description
Estimates functional propensity score (FPS) weights by maximising the empirical likelihood subject to covariate-balancing constraints and solving the resulting dual problem. The dual reduces to minimising log(sum_i exp(-theta' g_i)) over the unconstrained parameter vector theta, where g_i = [A_i, C_i, vec(A_i C_i')] stacks the balancing moments. The optimisation is solved via the BFGS quasi-Newton algorithm and uses the log-sum-exp trick for numerical stability.
Usage
.fps_compute_weights(
treat_scores,
conf_matrix,
normalize = TRUE,
tol = 1e-08,
maxit = 1000
)
Arguments
treat_scores |
n x L numeric matrix of FPC scores for the treatment. |
conf_matrix |
n x p numeric matrix of confounders (scalar + FPC scores of functional covariates). |
normalize |
Logical. If TRUE (default), standardise A and C before building g_i. |
tol |
Relative convergence tolerance for |
maxit |
Maximum number of BFGS iterations (default 1000). |
Value
A named list:
weights |
Numeric vector of length n. Positive, sums to 1. |
theta |
Optimal dual parameter vector. |
convergence |
Convergence code from |
value |
Optimal dual objective value. |
Internal FPCA helper
Description
Performs Functional Principal Component Analysis on a matrix of functional
observations using a B-spline basis expansion, then retains enough components
to explain at least 'pve' of the total variance. Accepts either an n x T
numeric matrix or an fd object from the fda package.
Usage
.fps_fpca(X, pve = 0.95, t_grid = NULL, domain = NULL, nbasis = NULL)
Arguments
X |
n x T numeric matrix of observed functional data (rows = subjects),
or an |
pve |
Proportion of variance explained threshold (default 0.95). |
t_grid |
Numeric vector of observation time points (length T).
Inferred from |
domain |
Numeric vector c(a, b) giving the domain of the functions.
Inferred from |
nbasis |
Number of B-spline basis functions. If NULL, set automatically. |
Value
A named list with components:
scr |
n x L matrix of FPC scores. |
efn |
T x L matrix of eigenfunction values evaluated on t_grid. |
mean |
Numeric vector of length T: mean function evaluated on t_grid. |
eval |
Numeric vector of eigenvalues (all, not just L). |
varprop |
Numeric vector of variance proportions for retained components. |
perc |
Cumulative variance proportions (all components). |
pca_fd |
The raw |
L |
Number of retained components. |
t_grid |
The t_grid used. |
domain |
The domain used. |
Extract treatment matrix from fps_weighting or fps_effect_estimation object
Description
Extract treatment matrix from fps_weighting or fps_effect_estimation object
Usage
.get_fpca_treatment(fps_obj)
Arguments
fps_obj |
An 'fps_weighting' or 'fps_effect_estimation' object. |
Value
The 'fpca_treatment' list element.
Balance plot: point-line chart of absolute correlations
Description
Balance plot: point-line chart of absolute correlations
Usage
.plot_balance(x)
Bootstrap slice plot for functional outcome (single point) with legend
Description
Bootstrap slice plot for functional outcome (single point) with legend
Usage
.plot_bootstrap_slice(x, point, which_domain)
Comparison plot: weighted vs unweighted, both with CI
Description
Comparison plot: weighted vs unweighted, both with CI
Usage
.plot_comparison(x)
Effect plot: mu(t) or mu(s,t) with CI
Description
Effect plot: mu(t) or mu(s,t) with CI
Usage
.plot_effect(x)
1-D effect plot with CI ribbon and legend
Description
1-D effect plot with CI ribbon and legend
Usage
.plot_effect_1d(x, use_weighted = TRUE, show_true = TRUE)
Surface heatmap for functional outcome
Description
Surface heatmap for functional outcome
Usage
.plot_effect_surface(
x,
beta,
ci_lower = NULL,
ci_upper = NULL,
title = "Causal Effect Surface"
)
FPCA panel: scree + eigenfunctions
Description
FPCA panel: scree + eigenfunctions
Usage
.plot_fpca_panel(fpca_res, title = "FPCA", max_efn = 4, domain_name = "t")
Significance plot
Description
Significance plot
Usage
.plot_significance(x, alpha)
Weights boxplot
Description
Weights boxplot
Usage
.plot_weights(x)
Reflected-percentile bootstrap confidence interval
Description
Reflected-percentile bootstrap confidence interval
Usage
.reflected_ci(obs_val, boot_vals, alpha = 0.05)
Arguments
obs_val |
Observed statistic (scalar). |
boot_vals |
Numeric vector of bootstrap replicates. |
alpha |
Significance level. |
Value
Named numeric vector c(lwr, obs, upr).
Estimate causal effect of a functional treatment
Description
Given the FPS weights produced by fps_weighting, estimates the
causal effect function \hat\mu(t) (scalar/binary outcome) or the
causal effect surface \hat\mu(s,t) (functional outcome) via weighted
least squares. Optional bootstrap inference is available.
Usage
fps_effect_estimation(
outcome,
fps_object,
outcome_t_grid = NULL,
outcome_domain = NULL,
outcome_domain_name = "t",
outcome_nbasis = NULL,
outcome_pve = 0.95,
treatment_pve = NULL,
bootstrap = FALSE,
B = 1000,
alpha = 0.05,
true_beta = NULL,
seed = NULL
)
Arguments
outcome |
Numeric vector (scalar/binary, length n) or n x T matrix (functional outcome). |
fps_object |
Object of class |
outcome_t_grid |
Numeric vector. Observation grid for functional
outcome. Required when |
outcome_domain |
Numeric |
outcome_domain_name |
Character. Name of the outcome domain (default
|
outcome_nbasis |
Integer or |
outcome_pve |
Numeric in (0, 1]. PVE threshold for outcome FPCA (default 0.95). |
treatment_pve |
Numeric or |
bootstrap |
Logical. If |
B |
Integer. Number of bootstrap resamples (default 1000). |
alpha |
Numeric. Significance level for bootstrap CIs (default 0.05). |
true_beta |
Optional. Numeric vector (scalar/binary) or matrix
(functional) containing the true causal effect, used for visual comparison
in plots and error metrics in |
seed |
Integer or |
Details
Scalar and binary outcomes.
The treatment FPC scores A_i are regressed on the outcome using
lm (scalar or binary, linear probability model) with the FPS
weights. The estimated effect function is then reconstructed as
\hat\mu(t) = \sum_k \hat\mu_k \phi_k(t).
Functional outcome.
For each outcome FPC component j, the regression
c_{ij} \sim A_i is solved with the FPS weights. The causal surface
is reconstructed as
\hat\mu(s,t) = \Phi_X \hat{B} \Phi_Y^\top
where \hat{B} collects the regression coefficients.
Bootstrap CIs. Scalar/binary: residual bootstrap, B resamples. Functional: pairs bootstrap, B resamples. Pointwise reflected-percentile confidence intervals are returned.
Value
An object of class "fps_effect_estimation", a named list with:
- outcome_type
Character:
"scalar","binary", or"functional".- beta
Estimated causal effect, evaluated on
t_grid(numeric vector for scalar/binary) or on thet_gridxoutcome_t_gridgrid (matrix for functional).- beta_unweighted
Same as
betabut from unweighted regression.- fpca_treatment
FPCA of the treatment used in estimation.
- fpca_outcome
NULLfor scalar/binary; FPCA list for functional outcome.- ci_lower, ci_upper
NULLifbootstrap = FALSE; otherwise lower and upper bootstrap CI bounds (same shape asbeta).- alpha
Significance level used.
- t_grid
Treatment domain grid.
- outcome_t_grid
NULLfor scalar/binary; outcome grid for functional.- domain_name
Treatment domain name.
- outcome_domain_name
Outcome domain name.
- true_beta
Passed through unchanged.
- fps_object
The input
fps_weightingobject.- call
The matched call.
See Also
fps_weighting, simulate_fps_data
Examples
dat <- simulate_fps_data(n = 2000, setting = "LL", seed = 1)
w <- fps_weighting(dat$X, dat$t_grid, c(0, 1), covariates = dat$C)
# Scalar outcome, no bootstrap
eff <- fps_effect_estimation(dat$Y, w, true_beta = dat$true_beta)
plot(eff, type = "effect")
plot(eff, type = "comparison")
# With bootstrap (small B for illustration)
eff_boot <- fps_effect_estimation(dat$Y, w, bootstrap = TRUE, B = 100,
true_beta = dat$true_beta, seed = 42)
plot(eff_boot, type = "significance")
Estimate functional propensity score weights
Description
Computes covariate-balancing weights for a functional treatment using the empirical-likelihood balancing framework of Ciardulli, S. and Fontana, N. (2026). Treatment is represented via FPCA (Karhunen–Loeve expansion); The treatment is first represented via Functional Principal Component Analysis (FPCA) through its Karhunen–Loeve expansion truncated at rank L; the resulting FPC scores and observed confounders are balanced by solving the dual of the empirical-likelihood problem via the BFGS quasi-Newton algorithm. Functional covariates enter the balancing step through their own FPC scores.
Usage
fps_weighting(
treatment,
treat_grid = NULL,
treat_domain = NULL,
domain_name = "s",
nbasis = NULL,
pve = 0.95,
covariates,
cov_grids = NULL,
cov_domains = NULL,
cov_nbasis = NULL,
cov_pve = 0.95,
normalize = TRUE,
tol = 1e-08,
maxit = 1000
)
Arguments
treatment |
n x T numeric matrix of observed treatment trajectories, or
an |
treat_grid |
Numeric vector of length T giving the observation grid
of the treatment. Required when |
treat_domain |
Numeric vector |
domain_name |
Character string naming the domain variable (default
|
nbasis |
Integer. Number of B-spline basis functions used for the
treatment FPCA. If |
pve |
Numeric in (0, 1]. Proportion of variance explained threshold for the treatment FPCA (default 0.95). |
covariates |
Either (a) an n x p numeric matrix of scalar covariates,
or (b) a named list with elements |
cov_grids |
A list of numeric vectors (one per functional covariate)
giving the observation grids. Required if |
cov_domains |
A list of numeric vectors |
cov_nbasis |
A list of integers (or |
cov_pve |
Numeric in (0, 1]. PVE threshold for functional covariate FPCA (default 0.95). |
normalize |
Logical. If |
tol |
Relative convergence tolerance for the BFGS optimiser (default 1e-8). |
maxit |
Maximum number of BFGS iterations (default 1000). |
Value
An object of class "fps_weighting", which is a named list
with the following components:
- weights
Numeric vector of length n. Positive weights summing to 1.
- fpca_treatment
List returned by the internal FPCA routine, containing FPC scores (
scr), eigenfunctions (efn), mean function (mean), eigenvalues (eval), variance proportions (varprop), cumulative PVE (perc), rawpca.fdobject (pca_fd), number of components retained (L), and thet_gridanddomainused.- fpca_covariates
List of FPCA results for functional covariates, or
NULLif none were supplied.- scalar_covariates
The n x p scalar covariate matrix used.
- conf_matrix
Full augmented confounder matrix fed to the optimiser (scalar covariates column-bound with FPC scores of functional covariates).
- convergence
Convergence code from
optim(0 = success).- domain_name
The domain name passed via
domain_name.- call
The matched call.
See Also
fps_effect_estimation, simulate_fps_data
Examples
dat <- simulate_fps_data(n = 2000, setting = "LL", seed = 1)
# Scalar covariates only (treat_domain inferred from treat_grid)
w <- fps_weighting(
treatment = dat$X,
treat_grid = dat$t_grid,
covariates = dat$C
)
print(w)
plot(w, type = "balance")
# Include one functional covariate
w2 <- fps_weighting(
treatment = dat$X,
treat_grid = dat$t_grid,
treat_domain = c(0, 1),
covariates = list(scalar = dat$C, functional = list(dat$D)),
cov_grids = list(dat$t_grid)
)
plot(w2, type = "balance")
Plot diagnostics for fps_effect_estimation objects
Description
Plot diagnostics for fps_effect_estimation objects
Usage
## S3 method for class 'fps_effect_estimation'
plot(
x,
type = "effect",
point = NULL,
which_domain = "treatment",
alpha = NULL,
max_efn = 4,
...
)
Arguments
x |
An object of class |
type |
Character. One of:
|
point |
Numeric scalar or vector. Time point(s) at which to slice the
effect surface (for |
which_domain |
Character. Either |
alpha |
Numeric. Significance level; defaults to |
max_efn |
Integer. Maximum eigenfunctions shown in FPCA panels. |
... |
Ignored. |
Value
A ggplot2 object.
Plot diagnostics for fps_weighting objects
Description
Plot diagnostics for fps_weighting objects
Usage
## S3 method for class 'fps_weighting'
plot(x, type = "balance", max_efn = 4, ...)
Arguments
x |
An object of class |
type |
Character. One of:
|
max_efn |
Integer. Maximum number of eigenfunctions to show in FPCA panels (default 4). |
... |
Ignored. |
Value
A ggplot2 object (or a list of them for
type = "fpca_covariates").
Print method for fps_effect_estimation objects
Description
Print method for fps_effect_estimation objects
Usage
## S3 method for class 'fps_effect_estimation'
print(x, ...)
Arguments
x |
An object of class |
... |
Ignored. |
Value
Invisibly returns x.
Print method for fps_weighting objects
Description
Print method for fps_weighting objects
Usage
## S3 method for class 'fps_weighting'
print(x, ...)
Arguments
x |
An object of class |
... |
Ignored. |
Value
Invisibly returns x.
Simulate functional propensity score data
Description
Generates a synthetic dataset for testing and illustrating the FPScausal workflow.
Usage
simulate_fps_data(
n = 200,
setting = c("LL", "LN", "NL", "NN"),
outcome_type = c("scalar", "functional"),
p_scalar = 3,
include_functional_cov = TRUE,
domain = c(0, 1),
seed = NULL
)
Arguments
n |
Integer. Number of subjects (default 200). |
setting |
Character. One of '"LL"', '"LN"', '"NL"', '"NN"', where the first letter controls the treatment-confounder relationship and the second controls the confounder-outcome relationship. Default '"LL"'. |
outcome_type |
Character. Either '"scalar"' or '"functional"'. Default '"scalar"'. |
p_scalar |
Integer. Number of scalar confounders. Default 3. |
include_functional_cov |
Logical. If 'TRUE' (default), include one functional covariate D(t) in the returned list. |
domain |
Numeric vector |
seed |
Integer or NULL. Random seed for reproducibility. |
Details
**Treatment** X_i(t) is built from 6 Fourier eigenfunctions with eigenvalues (16, 12, 8, 4, 1, 0.5). **Scalar confounders** C_i are 3-dimensional vectors whose relationship to X's FPC scores is either linear or quadratic. An optional **functional covariate** D_i(t) is generated from 4 Fourier components. The **scalar outcome** is Y_i = 1 + integral(beta(t) * X_i(t)) + g(C_i) + N(0,25), and the **functional outcome** is Y_i(t) = mu0(t) + integral(mu(s,t) * X_i(s) ds) + h(D_i) + GP_error.
The four settings ("LL", "LN", "NL", "NN") vary whether the treatment-confounder ("L"inear / "N"onlinear) and confounder-outcome ("L"inear / "N"onlinear) relationships are linear or quadratic.
Value
A named list with:
- X
n x 51 matrix. Observed treatment trajectories on [0,1].
- Y
If 'outcome_type = "scalar"': numeric vector of length n. If 'outcome_type = "functional"': n x 51 matrix.
- C
n x p_scalar matrix. Scalar confounders.
- D
n x 51 matrix. Functional covariate (if 'include_functional_cov = TRUE', else 'NULL').
- t_grid
Numeric vector of 51 equally-spaced points on [0,1].
- true_beta
True causal effect. For scalar outcome: numeric vector of length 51. For functional outcome: 51 x 51 matrix mu(s,t).
- setting
The 'setting' argument used.
- outcome_type
The 'outcome_type' argument used.
Examples
dat <- simulate_fps_data(n = 100, setting = "LL", outcome_type = "scalar",
seed = 42)
str(dat)
Summary method for fps_effect_estimation objects
Description
Summary method for fps_effect_estimation objects
Usage
## S3 method for class 'fps_effect_estimation'
summary(object, ...)
Arguments
object |
An object of class |
... |
Ignored. |
Value
Invisibly returns object.
Summary method for fps_weighting objects
Description
Summary method for fps_weighting objects
Usage
## S3 method for class 'fps_weighting'
summary(object, ...)
Arguments
object |
An object of class |
... |
Ignored. |
Value
Invisibly returns object.