README.md
In gsaot: Compute Global Sensitivity Analysis Indices Using Optimal Transport

gsaot

The package gsaot provides a set of tools to compute and plot Optimal Transport (OT) based sensitivity indices. The core functions of the package are:

ot_indices(): compute OT indices for multivariate outputs using different solvers for OT (network simplex, Sinkhorn, and so on).
ot_indices_wb(): compute OT indices for univariate or multivariate outputs using the Wasserstein-Bures semi-metric.
ot_indices_1d(): compute OT indices for univariate outputs using OT solution in one dimension. The package provides also functions to plot the resulting indices and the inner statistics.

install.packages("gsaot")

You can install the development version of gsaot from GitHub with:

# install.packages("devtools")
devtools::install_github("pietrocipolla/gsaot")

The sinkhorn and sinkhorn_log solvers in gsaot greatly benefit from optimization in compilation. To add this option (before package installation), edit your .R/Makevars file with the desired flags. Even though different compilers have different options, a common flag to enable a safe level of optimization is

CXXFLAGS+=-O2

More detailed information on how to customize the R packages compilation can be found in the R guide.

We can use a gaussian toy model with three outputs as an example:

library(gsaot)

N <- 1000

mx <- c(1, 1, 1)
Sigmax <- matrix(data = c(1, 0.5, 0.5, 0.5, 1, 0.5, 0.5, 0.5, 1), nrow = 3)

x1 <- rnorm(N)
x2 <- rnorm(N)
x3 <- rnorm(N)

x <- cbind(x1, x2, x3)
x <- mx + x %*% chol(Sigmax)

A <- matrix(data = c(4, -2, 1, 2, 5, -1), nrow = 2, byrow = TRUE)
y <- t(A %*% t(x))

x <- data.frame(x)

After having defined the number of partitions, we compute the sensitivity indices using different solvers. First, Sinkhorn solver and default parameters:

M <- 25

sensitivity_indices <- ot_indices(x, y, M)
#> Using default values for solver sinkhorn
sensitivity_indices
#> Method: sinkhorn 
#> 
#> Indices:
#>        X1        X2        X3 
#> 0.5656021 0.6131907 0.2594603 
#> 
#> Upper bound: 97.74799

Second, Network Simplex solver:

sensitivity_indices <- ot_indices(x, y, M, solver = "transport")
#> Using default values for solver transport
sensitivity_indices
#> Method: transport 
#> 
#> Indices:
#>        X1        X2        X3 
#> 0.4878673 0.5265048 0.1709062 
#> 
#> Upper bound: 97.74799

Third, Wasserstein-Bures solver, with bootstrap:

sensitivity_indices <- ot_indices_wb(x, y, M, boot = TRUE, R = 100)
sensitivity_indices
#> Method: wasserstein-bures 
#> 
#> Indices:
#>        X1        X2        X3 
#> 0.4591200 0.4921479 0.1073590 
#> 
#> Advective component:
#>         X1         X2         X3 
#> 0.28001929 0.31345309 0.09888083 
#> 
#> Diffusive component:
#>          X1          X2          X3 
#> 0.179100703 0.178694835 0.008478178 
#> 
#> Type of confidence interval: norm 
#> Number of replicates: 100 
#> Confidence level: 0.95 
#> Indices confidence intervals:
#>   Inputs     Index      low.ci    high.ci
#> 1     X1        WB 0.439721916 0.47851807
#> 2     X2        WB 0.475448734 0.50884711
#> 3     X3        WB 0.087350066 0.12736794
#> 4     X1 Advective 0.267679667 0.29235891
#> 5     X2 Advective 0.303581126 0.32332505
#> 6     X3 Advective 0.081559278 0.11620237
#> 7     X1 Diffusive 0.170346138 0.18785527
#> 8     X2 Diffusive 0.170261391 0.18712828
#> 9     X3 Diffusive 0.003887014 0.01306934
#> 
#> Upper bound: 97.85844

Fourth, we can use the package to compute the sensitivity map on the output:

sensitivity_indices <- ot_indices_smap(x, y, M)
sensitivity_indices
#>             X1         X2        X3
#> [1,] 0.5744058 0.04464331 0.1685419
#> [2,] 0.2909957 0.71343703 0.1274479