gsaot: Optimal Transport GSA in R

• gsaot is an R package that uses Optimal Transport to perform global sensitivity analysis from existing Monte Carlo data without requiring bespoke experimental designs.
• It implements classical OT, Wasserstein-Bures, and entropic OT indices to handle multivariate outputs and correlated inputs within a unified, model-agnostic framework.
• The package offers efficient numerical solvers, visualization tools, and bootstrapping capabilities that extend traditional variance-based sensitivity methods.

gsaot is an R package for Optimal Transport-based global sensitivity analysis (GSA). It is designed to estimate sensitivity indices from a dataset only, using any Monte Carlo sample of inputs and outputs as a post-processing step, rather than requiring bespoke experimental designs. Within the GSA landscape, it is positioned as model-agnostic, able to handle correlated inputs without additional tricks, and able to handle multivariate outputs such as vectors, time series, and spatial fields through a unified Optimal Transport (OT) framework. The package implements classical OT-based indices, Wasserstein-Bures (WB) indices, and entropic OT indices, together with visualization functions and multiple numerical solvers including network simplex and Sinkhorn-Knopp (Chiani et al., 24 Jul 2025).

1. Position in global sensitivity analysis

gsaot is situated against a background of traditional GSA methods such as Sobol’, Delta, and PAWN, which often require special sampling designs, assume independent inputs, or focus on scalar outputs. In contrast, gsaot is described as taking as input any Monte Carlo sample of inputs and outputs and performing sensitivity analysis as a post-processing step. It also accommodates correlated inputs and multivariate outputs within a single framework, and it decomposes the OT-based index into variance-based (“advective”), covariance-based (“diffusive”), and higher-order components (Chiani et al., 24 Jul 2025).

The package is also located explicitly within the R ecosystem. The source compares it with packages including sensitivity, sensobol, and multisensi, and with broader software ecosystems such as SAFE and SALib. The authors state that gsaot is, to their knowledge, the first R package to offer OT-based sensitivity indices that satisfy “zero-independence” and “max-functionality” while allowing correlated inputs and multivariate outputs. This suggests that its intended niche is not to replace classical variance-based GSA, but to extend sensitivity analysis to settings in which dependence structures and non-scalar outputs are intrinsic rather than exceptional.

A common misconception in this area is that generalized sensitivity analysis for structured outputs necessarily requires redesigning the simulation campaign. The package is presented specifically to counter that assumption: it estimates indices based on an existing dataset only, which makes it suitable for retrospective analysis of simulation archives or empirical datasets.

2. OT-based sensitivity indices and their properties

The package is built around a formulation in which sensitivity is measured through the expected separation between the unconditional output law and a conditional output law. Let $X \in \mathbb{R}^d$ be the input vector and $Y \in \mathbb{R}^k$ the output. If $\zeta(\cdot,\cdot)$ is a separation measure on probability distributions with $\zeta(P,P)=0$, then the local separation is

$\zeta(x_i) = \zeta(P_1, P_2(x_i)),$

where $P_1 = \mathrm{Law}(Y)$ and $P_2(x_i) = \mathrm{Law}(Y \mid X_i=x_i)$. The corresponding global importance is

$$\xi^\zeta(Y,X_i) = E_{X_i}\big[\zeta(P_1,P_2(X_i))\big].$$

In gsaot, $\zeta$ is taken from OT. Given a cost function $c(y,y') \ge 0$ on $Y \in \mathbb{R}^k$0, the Kantorovich OT cost is

$Y \in \mathbb{R}^k$1

If $Y \in \mathbb{R}^k$2 for a metric $Y \in \mathbb{R}^k$3 and $Y \in \mathbb{R}^k$4, then $Y \in \mathbb{R}^k$5 is the $Y \in \mathbb{R}^k$6-Wasserstein distance. Using OT as the separation measure yields the local and global quantities

$Y \in \mathbb{R}^k$7

$Y \in \mathbb{R}^k$8

The normalized OT-based sensitivity index is

$Y \in \mathbb{R}^k$9

with

$\zeta(\cdot,\cdot)$0

where $\zeta(\cdot,\cdot)$1 is an independent copy of $\zeta(\cdot,\cdot)$2 (Chiani et al., 24 Jul 2025).

The package emphasizes three formal properties of this normalized index. First, $\zeta(\cdot,\cdot)$3 if and only if $\zeta(\cdot,\cdot)$4, termed “zero-independence.” Second, $\zeta(\cdot,\cdot)$5 if and only if $\zeta(\cdot,\cdot)$6 is a deterministic function of $\zeta(\cdot,\cdot)$7, termed “max-functionality.” Third, the index is normalized to the interval $\zeta(\cdot,\cdot)$8. These properties are central because they define sensitivity in distributional rather than purely moment-based terms.

For the specific squared-Euclidean cost $\zeta(\cdot,\cdot)$9, the OT cost admits the Gelbrich decomposition

$\zeta(P,P)=0$0

where $\zeta(P,P)=0$1 and $\zeta(P,P)=0$2 are the means and covariances of $\zeta(P,P)=0$3, and $\zeta(P,P)=0$4 is a nonnegative remainder. This yields the decomposition

$\zeta(P,P)=0$5

with an “advective” term

$\zeta(P,P)=0$6

and a “diffusive” term

$\zeta(P,P)=0$7

The sum $\zeta(P,P)=0$8 is the Wasserstein-Bures index $\zeta(P,P)=0$9; equality holds for elliptical laws with a common generator. This structure gives the package a bridge between classical variance-based reasoning and broader distributional sensitivity analysis.

3. Entropic regularization and index variants

gsaot implements both classical OT-based indices, corresponding to $\zeta(x_i) = \zeta(P_1, P_2(x_i)),$0, and entropic OT indices, corresponding to $\zeta(x_i) = \zeta(P_1, P_2(x_i)),$1. The entropic formulation augments the OT problem with a Kullback-Leibler penalty:

$\zeta(x_i) = \zeta(P_1, P_2(x_i)),$2

The source characterizes entropic OT as faster to compute via Sinkhorn-Knopp and as converging to the classical indices as $\zeta(x_i) = \zeta(P_1, P_2(x_i)),$3 (Chiani et al., 24 Jul 2025).

This distinction matters both theoretically and computationally. Classical OT preserves the exact “zero-independence” and “max-functionality” properties stated above. Entropic indices retain normalization and max-functionality, but only an approximate zero-independence. That caveat is important: regularization improves numerical tractability, but it perturbs the exact sensitivity semantics. A plausible implication is that the choice between classical and entropic OT in gsaot is not merely an implementation detail; it determines whether exact distributional separation or a regularized approximation is being used as the sensitivity criterion.

The package therefore supports multiple variants of OT-based sensitivity analysis under one interface: exact classical OT, the moment-based WB special case, and entropic OT. This unification is one of its defining features.

4. Estimation from sampled data

The estimation strategy in gsaot is explicitly designed to avoid double Monte Carlo loops. Given a sample of $\zeta(x_i) = \zeta(P_1, P_2(x_i)),$4 realizations $\zeta(x_i) = \zeta(P_1, P_2(x_i)),$5 of $\zeta(x_i) = \zeta(P_1, P_2(x_i)),$6, the package uses a “given-data” U-statistic approach. For each input $\zeta(x_i) = \zeta(P_1, P_2(x_i)),$7, its empirical support is partitioned into $\zeta(x_i) = \zeta(P_1, P_2(x_i)),$8 bins $\zeta(x_i) = \zeta(P_1, P_2(x_i)),$9. Let $P_1 = \mathrm{Law}(Y)$0 denote the empirical law of $P_1 = \mathrm{Law}(Y)$1 for the $P_1 = \mathrm{Law}(Y)$2 points in bin $P_1 = \mathrm{Law}(Y)$3, and let $P_1 = \mathrm{Law}(Y)$4 denote the unconditional empirical law. The estimator of the numerator is

$P_1 = \mathrm{Law}(Y)$5

The denominator $P_1 = \mathrm{Law}(Y)$6 is estimated by the U-statistic

$P_1 = \mathrm{Law}(Y)$7

The final estimator is then $P_1 = \mathrm{Law}(Y)$8 (Chiani et al., 24 Jul 2025).

Two special-case accelerations are implemented. For $P_1 = \mathrm{Law}(Y)$9 and $P_2(x_i) = \mathrm{Law}(Y \mid X_i=x_i)$0, one sorts the unconditional and conditional samples, takes $P_2(x_i) = \mathrm{Law}(Y \mid X_i=x_i)$1-th powers of the differences of matching quantiles, and averages. For elliptical families with $P_2(x_i) = \mathrm{Law}(Y \mid X_i=x_i)$2, closed-form expressions based on sample means and covariances are available via the Gelbrich formula.

The partition parameter is operationally important. The package documentation recommends choosing the number of partitions M so that there are at least 100 points per bin; for $P_2(x_i) = \mathrm{Law}(Y \mid X_i=x_i)$3 from 200 to 10,000, M=N/100 is given as a safe rule, and for very large $P_2(x_i) = \mathrm{Law}(Y \mid X_i=x_i)$4, results plateau for $P_2(x_i) = \mathrm{Law}(Y \mid X_i=x_i)$5. These recommendations are not abstract tuning advice; they are part of the estimator design, because the conditional empirical laws are defined through binning.

5. Solvers, computational trade-offs, and software interface

gsaot exposes a set of OT solvers for different regimes, using the transport package for $P_2(x_i) = \mathrm{Law}(Y \mid X_i=x_i)$6 and C++ implementations for $P_2(x_i) = \mathrm{Law}(Y \mid X_i=x_i)$7 (Chiani et al., 24 Jul 2025).

Solver	OT type	Package
1d-quantile	classic	gsaot
WB-closed	classic	gsaot
network simplex	classic	transport
primal-dual	classic	transport
revised simplex	classic	transport
shortlist	classic	transport
Sinkhorn	entropic	gsaot
Sinkhorn-stable	entropic	gsaot

The source gives specific computational characterizations. Network simplex and related simplex variants solve the linear program exactly. Their worst-case complexity is roughly $P_2(x_i) = \mathrm{Law}(Y \mid X_i=x_i)$8 for $P_2(x_i) = \mathrm{Law}(Y \mid X_i=x_i)$9 support points, though the implementation in transport is described as highly optimized and able to handle up to a few thousand points per partition. Sinkhorn-Knopp performs matrix-scaling updates for the $$\xi^\zeta(Y,X_i) = E_{X_i}\big[\zeta(P_1,P_2(X_i))\big].$$0-regularized problem. Each iteration costs $$\xi^\zeta(Y,X_i) = E_{X_i}\big[\zeta(P_1,P_2(X_i))\big].$$1 to update two scaling vectors, and the theoretical convergence rate is $$\xi^\zeta(Y,X_i) = E_{X_i}\big[\zeta(P_1,P_2(X_i))\big].$$2 iterations; in practice, a few hundred to a few thousand iterations suffice for $$\xi^\zeta(Y,X_i) = E_{X_i}\big[\zeta(P_1,P_2(X_i))\big].$$3. A stabilized variant is provided for numerical robustness at small $$\xi^\zeta(Y,X_i) = E_{X_i}\big[\zeta(P_1,P_2(X_i))\big].$$4.

The package also provides a concise user-facing workflow. After installation from CRAN and loading with library(gsaot), the basic steps are to generate or load a Monte Carlo sample x and y, choose the number of partitions M, and call one of three front-end functions:

• ot_indices_1d(x, y, M, ...) for scalar y and $c(y,y&#39;)=|y-y&#39;|^p`</li> <li><code>ot_indices_wb(x, y, M, ...)</code> for multivariate <code>y</code> with squared-Euclidean cost</li> <li><code>ot_indices(x, y, M, solver=..., solver_optns=...)</code> for general classical or entropic OT</li> </ul> <p>The package then supports <code>print()</code>, <code>plot()</code>, <code>plot_comparison()</code>, and <code>plot_separations()</code> for inspection and visualization. Recommended solver settings include <code>solver_optns=list(method=&quot;networkflow&quot;)</code> for network simplex, and for entropic OT, defaults such as <code>epsilon=1e-2</code> or <code>1e-3</code>, <code>maxErr=1e-3</code>, and <code>numIterations=1e5</code>. Bootstrapped <a href="https://www.emergentmind.com/topics/confidence" title="" rel="nofollow" data-turbo="false" class="assistant-link" x-data x-tooltip.raw="">confidence</a> intervals are supported with <code>boot=TRUE, R=1000, conf=0.95, type=&quot;norm&quot;</code>, and parallelization is available via <code>parallel=&quot;multicore&quot;, ncpus=...</code>.</p> <h2 class='paper-heading' id='worked-examples-and-empirical-usage'>6. Worked examples and empirical usage</h2> <p>The package paper presents three worked examples that illustrate how the framework behaves in distinct modeling contexts (<a href="/papers/2507.18588" title="" rel="nofollow" data-turbo="false" class="assistant-link" x-data x-tooltip.raw="">Chiani et al., 24 Jul 2025</a>).</p> <p>The first example is a linear Gaussian toy model with $\xi^\zeta(Y,X_i) = E_{X_i}\big[\zeta(P_1,P_2(X_i))\big].$5, where $\xi^\zeta(Y,X_i) = E_{X_i}\big[\zeta(P_1,P_2(X_i))\big].$6 is a $\xi^\zeta(Y,X_i) = E_{X_i}\big[\zeta(P_1,P_2(X_i))\big].$7 matrix and $\xi^\zeta(Y,X_i) = E_{X_i}\big[\zeta(P_1,P_2(X_i))\big].$8. Using $\xi^\zeta(Y,X_i) = E_{X_i}\big[\zeta(P_1,P_2(X_i))\big].$9 points and $\zeta$$0, the package is applied through <code>ot_indices_wb</code>, <code>ot_indices</code> with network simplex, and <code>ot_indices</code> with Sinkhorn. The reported interpretation is that all three methods recover the analytic ranking$$\zeta$1. The example also serves to calibrate the regularization parameter: Sinkhorn with large $\zeta$2 such as $\zeta$3 overestimates the indices, whereas small $\zeta$4 matches the exact ones.

  The second example is a spruce budworm and forest ODE model with a 3-dimensional dynamic state—budworm density $\zeta$5, size $\zeta$6, and energy $\zeta$7—and 10 inputs. The system is solved over $\zeta$8 with deSolve, using $\zeta$9 parameter sets and parallelization via future and furrr. OT indices are computed for output curves, and the interpretation reported in the source is that $c(y,y') \ge 0$0 (max budworm density) and $c(y,y') \ge 0$1 (branch growth rate) dominate $c(y,y') \ge 0$$2, while many other parameters fall below the dummy threshold and can be considered negligible. The function <code>plot_separations(idx_B, ranking=3)</code> is used to visualize$$c(y,y') \ge 0$3 across the support of the most influential inputs, showing which regions of the input support drive most of the change.

  The third example is a simple DICE climate module based on a discrete-time carbon-cycle and temperature anomaly model from DICE2016, with 9 uncertain inputs and an endogenous emission path from 2015 to 2100. The analysis focuses on $c(y,y') \ge 0$4 trajectories from $c(y,y') \ge 0$5 runs with $c(y,y') \ge 0$6, and uses a custom Minkowski $c(y,y') \ge 0$7 cost $c(y,y') \ge 0$8 together with entropic OT and bootstrapping. The interpretation reported is that $c(y,y') \ge 0$9, $Y \in \mathbb{R}^k$00, $Y \in \mathbb{R}^k$01, and $Y \in \mathbb{R}^k$02 are the top levers, while other inputs remain near the zero-line and may be fixed.

  These examples collectively show that gsaot is intended for scalar, multivariate, dynamic, and trajectory-valued outputs. They also demonstrate that the package is not restricted to squared-Euclidean costs: custom costs can be supplied for domain-specific notions of output discrepancy.

  7. Scope, interpretation, and limitations

  The package is presented as a flexible toolkit that complements moment-based and distribution-based methods by handling complex multivariate outputs, allowing correlated inputs, unifying classical OT, Wasserstein-Bures, and entropic OT, and supporting both closed-form and numerical solvers, as well as bootstrapping and parallelization (Chiani et al., 24 Jul 2025).

  Its interpretation framework is correspondingly broader than conventional scalar-output GSA. The advective term is described as akin to variance-based Sobol’, the diffusive term captures covariance-level changes, and the higher-order remainder covers distributional effects beyond first and second moments. This suggests that the package is especially suited to applications where sensitivity is expressed through changes in the geometry of the output distribution rather than through scalar dispersion alone.

  Several practical cautions are explicit in the source. Users are advised to tune M to the sample size, to tune $Y \in \mathbb{R}^k$03 and tolerances for Sinkhorn, and to use dummy thresholds to guard against numerical noise. Another important clarification is that entropic indices provide only approximate zero-independence, so a near-zero value under regularization should not be interpreted identically to the exact classical OT case.

  In this sense, gsaot occupies a specific methodological position. It does not merely add another set of scalar sensitivity estimators; it formulates GSA as a problem of comparing unconditional and conditional output laws through Optimal Transport, and then packages that formulation into a usable R interface for both exact and regularized computation.