{"id":"a64870ff-ec5d-43b9-b0c3-6901ac126130","arxiv_id":"2507.18588","paper_version":1,"verdict":"CONDITIONAL","confidence":"HIGH","novelty_score":5.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":3,"one_line_summary":"gsaot is a new R package implementing Optimal Transport-based sensitivity indices for multivariate outputs and correlated inputs using given-data estimators.","lead":"This paper introduces gsaot, an R package that computes Optimal Transport-based global sensitivity indices from input-output datasets. It lets analysts rank which model inputs drive multivariate or time-dependent outputs, even when inputs are correlated.","discovery_kind":"extension","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Eq. (16)'s binning estimator is downward biased by convexity of the OT cost; with no consistency or bias analysis, the package's finite-sample indices—and hence the central claim—rest on an unvalidated estimator.","rationale":"The reader's weakest assumption identifies exactly the most load-bearing point: the given-data estimator in Eq. (16) is the algorithmic heart of every index the package reports, and its statistical behavior is not established. I agree with that identification. The concern is not merely 'more theory would be nice': because the OT cost is convex in its marginals, the binning step is provably biased downward, and the bias is input-dependent. This can change the very rankings the package is designed to produce. The Gaussian example in §3.1 is intended as a validation, but the mismatch between the textual A matrix and the code's A matrix means the reported indices are not compared to the stated ground truth; the closeness of the numbers is coincidental or at least unverified. A synthetic test with a known index and engineered within-bin heterogeneity would directly measure the bias and show whether it is practically significant. Until such evidence is provided, a conditional acceptance is appropriate; the package may be correct, but the core estimator's reliability is an open question. The paper's novelty claim (first package with these features) is plausible and not contradicted by my review, but it is secondary: even a first package must have trustworthy indices.","tokens_in":20369,"tokens_out":6619,"duration_ms":67271,"concrete_test":"Use a model with known pointwise conditional distributions and strong within-bin heterogeneity, e.g., X1~U(0,1), Y = X1 + 2 sin(2π X1) + ε with ε~N(0,0.01). Compute the true ξ_K by numerical integration over 500 fine bins. Then run gsaot::ot_indices with the default networkflow solver on N=2000 samples for H=10,20,50,100, replicating 100 times. If the median estimated index falls below the true value by more than the bootstrap CI width for every H, the downward bias in Eq. (16) is confirmed, and the paper's recommendation to use the package without bias correction is not supported.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The central claim (Section 4) is that gsaot is the first package handling multivariate outputs and correlated inputs from given data. The correctness of this claim depends on the estimator in Eq. (16): it replaces the conditioning event Xi=xi with a partition cell Xi∈X^h_i and averages empirical OT costs over H bins. The paper cites Borgonovo et al. (2016) for given-data estimation but provides no bias, variance, or consistency analysis for OT-based indices. This matters because the OT cost is convex in each marginal (the Kantorovich value function), so K(PY, P_Y|Xi∈bin) ≤ (1/P(bin))∫_{bin} K(PY, P_Y|Xi=xi) dP(xi). Binning therefore biases the numerator downward, with the bias depending on within-bin heterogeneity of the conditional distributions. Different inputs will generally have different within-bin heterogeneity, so rankings can change. For fixed H, the estimator converges to the biased bin-average, not to ξ_K, as N→∞. The heuristics offered (at least 100 points per bin; Strong and Oakley plateau; dummy threshold) do not quantify this error. The §3.1 Gaussian validation cannot rescue the argument because the code's matrix A (c(4,-2,1,2,5,-1)) differs from the text's A ([[4,2,3],[2,5,-1]]), so the comparison against Table 2 is not a clean accuracy check. This leaves the package's headline capability unsupported for finite N and H.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The manuscript presents gsaot, an R package for global sensitivity analysis based on Optimal Transport (OT) distances. It defines OT-based sensitivity indices from Borgonovo et al. (2024), including the Wasserstein-Bures decomposition into advective, diffusive, and higher-moment components, and implements given-data estimators that partition each input into H bins and average OT costs between the marginal output distribution and conditional output distributions. The package supports one-dimensional outputs, Gaussian/Bures closed-form estimation, classical OT via the transport package, and entropic OT via custom C++ Sinkhorn solvers, together with bootstrap confidence intervals, plotting, and an irrelevance-threshold tool based on dummy inputs. Three applications are shown: a linear Gaussian model with known analytical values, a spruce budworm ODE model with time-dependent output, and a DICE-based climate model with a custom ground cost.","tokens_in":20705,"tokens_out":9899,"duration_ms":97754,"significance":"If the estimates are reliable, the package is a useful contribution: it is model-agnostic, works directly from a dataset, supports correlated inputs and multivariate outputs, and implements indices with zero-independence and max-functionality. The paper also ships reproducible R code, makes the package available on CRAN, and checks one analytical benchmark. The main weakness is that the reliability of the central given-data estimator is not established beyond heuristics, and the benchmark contains a matrix inconsistency, so the validation currently does not close that gap.","major_comments":[{"comment":"The given-data estimator in Eq. (16) is the backbone of the package, yet the paper provides no bias, variance, or consistency analysis for it. For fixed H, by convexity of the Kantorovich cost in the second marginal, K(P_Y, P_{Y|X_i∈X^h_i}) ≤ (1/P(X^h_i)) ∫_{X^h_i} K(P_Y, P_{Y|X_i=x_i}) dP(x_i), so the binning estimator is downward biased and converges as N→∞ to a bin-average, not to ξ_K(Y,X_i). Because different inputs generally have different within-bin heterogeneity, rankings can change even in large samples. The heuristics in the text (at least 100 points per bin, the Strong-Oakley plateau, the dummy threshold) do not quantify this error. The Gaussian example in Section 3.1 cannot settle the question because the benchmark itself is compromised (see next comment). I ask the authors to add either a consistency result with H=H(N)→∞, or a systematic numerical study showing that the bias is negligible and rankings are stable for representative nonlinear and non-Gaussian models.","section":"Section 3.1"},{"comment":"The validation against Table 2 is not a clean check because the matrix used in the R code differs from the matrix stated in the text. The text defines A = [[4,2,3],[2,5,-1]], but the code 'A <- matrix(data = c(4,-2,1,2,5,-1), nrow = 2, byrow = TRUE)' produces A = [[4,-2,1],[2,5,-1]]. As a result, the column for X3 (and the correlated-input structure) differs between the claimed analytical model and the simulated data. Please reconcile the two definitions and report which matrix was actually used to generate the results; the comparison between ot_indices_wb, networkflow, and Sinkhorn remains informative only after this correction.","section":"Section 3.1"},{"comment":"Equation (10) is missing the normalization by M_K[Y] in the remainder term. Since the first two terms are defined in Eqs. (11)–(12) as normalized indices, the decomposition should read ι_K = ι_V + ι_Σ + M_K[Y]^{-1} E[Γ(P_Y,P_{Y|X_i})]; as written, the right-hand side mixes normalized and unnormalized quantities. This is a local fix, but the equation is part of the theoretical basis of the package's decomposition property and should be corrected.","section":"Equation (10)"}],"minor_comments":[{"comment":"The one-dimensional formula in Eq. (17) is written as a sum over j=1,...,N of the p-th differences of empirical quantile functions at j/N, but the conditional sample has N_h < N points; please state how the conditional quantile function is evaluated at the finer grid and confirm that the expression equals the L_p^p Wasserstein cost between the two empirical measures.","section":"Eq. (17)"},{"comment":"The number of partitions is denoted H in Eq. (16) and in the statement 'partitioned into H subsets,' but the code and the later discussion use M throughout; please unify the notation.","section":"Section 2.4 / code examples"},{"comment":"The text says the example uses 'the fourth power of Minkowski distance of order 4 (L3^3)', while the code custom_metric uses dist(..., method='minkowski', p=3)^3; please correct the verbal description to match the implemented L_3^3 cost.","section":"Section 3.3"},{"comment":"The third panel in Figure 4 is labeled K_e, and the text says the plot shows the three most important inputs; the subsequent sentence about local separations for T_e should refer to a variable that is actually plotted, or the figure/ranking should be changed.","section":"Section 3.2, Figure 4"},{"comment":"There are several spelling and style slips, e.g., 'inizialize' should be 'initialize', and the sentence in Section 3.3 beginning 'We inizialize' should be reworded; a careful proofread is recommended.","section":"Throughout"}],"recommendation":"major_revision","confidential_remarks":"The package appears genuinely useful and the authors are credible, but the estimator bias issue plus the matrix mismatch in the headline benchmark are enough that I cannot accept as is. I would send back for major revision rather than reject: the issues are fixable with an added numerical study and corrections."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Short take: gsaot is a solid, useful R package that fills a real gap—no other R package computes OT-based sensitivity indices for multivariate outputs with correlated inputs from raw data. The paper deserves a serious referee, but there are a couple of documentation errors that need fixing, and the estimation section should be honest about the binning bias.\n\nWhat's new and good: The package is genuinely new as a software artifact. It wraps the transport package's network simplex and other solvers, implements a C++ Sinkhorn, includes the Wasserstein-Bures closed form and the 1D quantile formula, and provides plotting, bootstrap CIs, and a dummy-variable threshold. The three examples cover a Gaussian test case, a dynamic model, and a climate model. The theory recap is clean, and the code is on CRAN.\n\nSoft spots: Section 3.1 has a clear mismatch—the matrix A in the text ([[4,2,3],[2,5,-1]]) does not match the A in the R code (c(4,-2,1,2,5,-1)). So the printed estimates are not a clean comparison against Table 2. This must be fixed, either by changing the text or the code. Equation (10) also drops the M_K[Y] normalization on the higher-moment term. Minor typo. More substantively, the given-data estimator in Eq. (16) replaces the conditioning event Xi=xi with a bin. Since the OT cost is convex in each marginal, Jensen's inequality gives K(PY, P_Y|Xi in bin) <= the bin average of K(PY, P_Y|Xi=xi). So for fixed H, the estimator is downward biased as N grows. The paper's heuristics (100 points per bin, the plateau, the dummy threshold) do not quantify this bias. This is not fatal for a software paper—it is an inherited property of the given-data approach—but it should be acknowledged, preferably with a short consistency or bias discussion. Also, in the climate example, \"fourth power of Minkowski distance of order 4 (L3_3)\" is contradictory; the code uses p=3 and raises to the third power.\n\nOverall: the central claim about being the first package of this kind holds up. The package will be useful to UQ practitioners, and the paper deserves peer review after minor revision and a corrected validation. I would cite it if I work on sensitivity analysis.","headline":"A genuinely useful R package for OT-based sensitivity analysis, blighted by a few fixable reproducibility mistakes and an unacknowledged binning bias.","tokens_in":21201,"tokens_out":5138,"would_cite":true,"duration_ms":52430,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["62H20","62-04","49Q22"],"pacs":[],"model":"deepseek-v4-flash","headline":"Optimal transport ranks model inputs straight from data","keywords":["global sensitivity analysis","optimal transport","Wasserstein distance","multivariate output","correlated inputs","given-data estimation","entropic regularization","R package"],"falsifier":"Take a model with known true OT indices (for example, a Gaussian model as in the paper) and add a small-amplitude, high-frequency oscillation to the conditional output as a function of $X_i$, so that within every chosen bin the conditional distribution varies appreciably. Then compute the package's estimate at increasing $N$ with $H$ fixed; if the estimate stays biased away from the truth as $N$ grows while an exact-conditioning Monte Carlo estimate converges, the binning approximation—and with it the package's reliability—is refuted.","tokens_in":20202,"feed_emoji":"📦","tokens_out":6305,"duration_ms":60328,"temperature":0.7,"pith_summary":"The paper introduces gsaot, an R package that computes Optimal Transport-based global sensitivity indices directly from an input-output dataset, without requiring a special sampling design or access to the simulation model. Its central claim is that gsaot is the first package able to handle both multivariate outputs and correlated inputs in this given-data setting. The indices it computes are normalized measures of statistical association with the properties of zero-independence and max-functionality, and they decompose an input's effect into components acting on the mean, the covariance, and higher-order moments of the output distribution. If the estimators work as claimed, any stored collection of model runs—from a climate simulator, an ODE model, or field data—can be post-processed into input rankings and local influence curves.","feed_headline":"Optimal transport ranks model inputs straight from data","feed_subtitle":"gsaot post-processes any input-output sample into normalized sensitivity indices, even with multivariate outputs and correlated inputs.","key_machinery":"The load-bearing object is the OT-based sensitivity index $\\iota_K(Y,X_i) = E[K(P_Y, P_{Y|X_i})] / M_K[Y]$, where $K$ is the optimal transport cost under a user-chosen ground metric and $M_K[Y] = E[c(Y, Y')]$ is the expected cost between independent output replicas. The estimator replaces the conditioning event $X_i = x_i$ with membership in one of $H$ bins covering the input support, then averages the empirical OT costs over bins. This single plug-in device connects three computational routes: a closed-form quantile solution in one dimension, the Wasserstein-Bures closed form for elliptical distributions, and general empirical OT solved by network-flow or Sinkhorn algorithms.","core_discovery":"The paper's central claim is that the Wasserstein distance between the marginal output distribution and the output distribution conditional on an input defines a global sensitivity index that can be estimated from data alone, and that this index preserves the strong axiomatic properties of genuine measures of statistical association. Specifically, fixing input $X_i$ at $x_i$, the local separation $\\zeta_K(x_i) = K(P_Y, P_{Y|X_i=x_i})$ is averaged over $X_i$ and divided by the maximum attainable cost $M_K[Y]$, yielding $\\iota_K(Y,X_i)$ in $[0,1]$. The index is zero exactly when $Y$ and $X_i$ are independent and one exactly when $Y$ is a deterministic function of $X_i$. The paper implements this for one-dimensional outputs through a quantile formula, for elliptical or Gaussian outputs through the closed-form Wasserstein-Bures metric, and for general multivariate outputs through numerical optimal transport solvers, including entropic Sinkhorn regularization as a fast approximation.","pith_inferences":["If the binning estimator is unbiased in practice, gsaot could be used as a screening tool on archival datasets in fields where experimenter-controlled sampling is impossible, such as epidemiology or observational climate analysis.","The mean/covariance decomposition suggests a diagnostic not developed in the paper: comparing advective with diffusive components across inputs can reveal which inputs drive heteroscedasticity, something variance-based indices miss.","Because the population index satisfies max-functionality, the package could double as an independence test for multivariate dependence, but the paper does not provide distributional results or calibration for such testing.","The dummy-variable irrelevance threshold is heuristic; a formal false-discovery or calibration analysis would strengthen its use for variable screening, but the paper does not propose one."],"forward_implications":["Analysts can rank input importance for any black-box model from a stored dataset, removing the need for purpose-built sampling designs or model reruns.","Multivariate and time-series outputs are treated as whole distributions, so an input that changes the shape or covariance of the output is detected even when marginal means do not move.","Correlated inputs do not invalidate the index; the ranking and its zero and one extremes remain well defined.","The index decomposes into advective (mean-shift), diffusive (covariance-shift), and higher-order components, allowing analysts to say not only how much but in what way an input affects the output.","Entropic regularization provides a fast, lower-bounded approximation that converges to the classical index as the regularization parameter goes to zero, making large-$N$ analyses feasible."],"supporting_citations":[{"why":"Introduces the OT-based sensitivity indices, their normalization, decomposition, and axiomatic properties that gsaot implements.","marker":"Borgonovo et al. (2024)"},{"why":"Provides the closed-form $L^2$ Wasserstein metric between elliptical distributions used in the Wasserstein-Bures estimator.","marker":"Gelbrich (1990)"},{"why":"Supplies the one-dimensional quantile formula and the general computational optimal transport toolkit underlying the solvers.","marker":"Peyré and Cuturi (2020)"},{"why":"Introduces entropic regularization and the Sinkhorn-Knopp algorithm used for fast entropic indices.","marker":"Cuturi (2013)"},{"why":"Establishes Wasserstein-distance-based measures of association with zero-independence and max-functionality, the properties claimed for gsaot's indices.","marker":"Wiesel (2022)"},{"why":"Documents the plateau effect used as heuristic guidance for choosing the number of bins $H$.","marker":"Strong and Oakley (2013)"},{"why":"Provides the common rationale for global sensitivity measures and the given-data estimation approach used in the estimator.","marker":"Borgonovo et al. (2016)"},{"why":"Introduces the dummy-variable practice that gsaot turns into an irrelevance threshold.","marker":"Noacco et al. (2019)"}],"fun_headline_variants":["gsaot: optimal transport sensitivity indices in R","Wasserstein distance ranks input influence","Model-agnostic sensitivity via optimal transport","Data-only sensitivity indices via optimal transport","gsaot: measure input importance via Wasserstein"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The whole method rests on replacing \"fix $X_i$ at a precise value\" with \"$X_i$ falls in one of a few broad bins\" and averaging the results over those bins; if this replacement is biased at finite sample sizes, the computed indices and rankings are unreliable.","fun_headline_variants_meta":{"raw":{"variants":["gsaot: optimal transport sensitivity indices in R","Wasserstein distance ranks input influence","Model-agnostic sensitivity via optimal transport","Data-only sensitivity indices via optimal transport","gsaot: measure input importance via Wasserstein"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.001165,"raw_usage":{"total_tokens":4765,"prompt_tokens":832,"completion_tokens":3933,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":448,"completion_tokens_details":{"reasoning_tokens":3865}},"tokens_in":448,"tokens_out":3933,"duration_ms":27705,"temperature":1.0,"reasoning_tokens":3865,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-15T18:10:06.695039+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Take a model with known true OT indices (for example, a Gaussian model as in the paper) and add a small-amplitude, high-frequency oscillation to the conditional output as a function of $X_i$, so that within every chosen bin the conditional distribution varies appreciably. Then compute the package's estimate at increasing $N$ with $H$ fixed; if the estimate stays biased away from the truth as $N$ grows while an exact-conditioning Monte Carlo estimate converges, the binning approximation—and with it the package's reliability—is refuted.","supporting_citations":[{"cited_title":"Global Sensitivity Analysis via Optimal Transport","cited_arxiv_id":null,"evidence_quote":"Introduces the OT-based sensitivity indices, their normalization, decomposition, and axiomatic properties that gsaot implements."},{"cited_title":"Sinkhorn Distances : Lightspeed Computation of Optimal Transport","cited_arxiv_id":null,"evidence_quote":"Introduces entropic regularization and the Sinkhorn-Knopp algorithm used for fast entropic indices."},{"cited_title":"Measuring Association with Wasserstein Distances","cited_arxiv_id":null,"evidence_quote":"Establishes Wasserstein-distance-based measures of association with zero-independence and max-functionality, the properties claimed for gsaot's indices."},{"cited_title":"An efficient method for computing partial expected value of perfect information for correlated inputs","cited_arxiv_id":null,"evidence_quote":"Documents the plateau effect used as heuristic guidance for choosing the number of bins $H$."},{"cited_title":"A Common Rationale for Global Sensitivity Measures and Their Estimation","cited_arxiv_id":null,"evidence_quote":"Provides the common rationale for global sensitivity measures and the given-data estimation approach used in the estimator."},{"cited_title":"MATLAB / R workflows to assess critical choices in Global Sensitivity Analysis using the SAFE toolbox","cited_arxiv_id":null,"evidence_quote":"Introduces the dummy-variable practice that gsaot turns into an irrelevance threshold."}],"review_version":2}