Review of the Zhang et al. (2026) GSA products

Scope: SpyDust/paper-products/AME_feature_analysis/GSA_and_surrogates/, the six GSA_{DC,HII,MC}[_rev].ipynb notebooks and obs_catalogue_vs_MonteCarlo_samples.ipynb.

Nothing in the SpyDust repository has been modified. This document reports findings only.

0. What the notebooks actually ran

All seven notebooks call gsa_pipeline with an identical configuration:

scaler="standard", kernel_kind="rbf", ard=True, length_scale_init=1.0,
N_sobol=1024, gp_random_state=0, drop_const_atol=0.0

length_scale_bounds is not passed anywhere — it was not settable before v0.2.0 and defaulted to (1e-2, 1e2).

The choice of scaler="standard" over the package default "minmax" turns out to be the single most consequential line in these notebooks. See §4.

1. Exact reproducibility

obs_catalogue_vs_MonteCarlo_samples.ipynb reads pre-computed inputs (data/{MC,DC,HII}_dist_ensem_{params,features}.txt, 5000 x 7 and 5000 x 2), so it re-runs exactly with no SpyDust dependency. All six published tables — three ensembles x two features — reproduce cell for cell:

table max abs. difference vs published row order
obs / MC / nu_peak 0.0005 identical
obs / MC / width 0.0005 identical
obs / DC / nu_peak 0.0005 identical
obs / DC / width 0.0005 identical
obs / HII / nu_peak 0.0005 identical
obs / HII / width 0.0005 identical

0.0005 is the rounding of the published LaTeX to three decimals, i.e. the agreement is exact. Every number quoted below is therefore measured on the same inputs that produced the paper’s tables, not on a re-simulation.

2. Item 1 — does the ARD upper-bound clamp change AggRank?

Short answer: not in a way that changes any conclusion, but the order is not completely stable either. The six obs_catalogue tables are unchanged under both a relaxed upper bound (1e-2, 1e4) and sklearn’s (1e-5, 1e5). Three of the six _rev tables — the ones that actually contain clamped length scales — do reorder, always as an adjacent swap among bottom-ranked parameters with ST ≈ 0. Detail below.

Two things are worth separating.

Where the clamp actually is. None of the six obs_catalogue tables contains a clamped length scale: the smallest 1/ARD_LS anywhere in them is 0.0234 (obs / DC / width, delta, i.e. ARD_LS = 42.7). The 1/ARD_LS = 0.010 entries — ARD_LS exactly 100, the upper bound — are in the main six GSA_{DC,HII,MC}[_rev].ipynb notebooks, on y, Chi, xC, xh and T.

What relaxing the bound does anyway. Re-fitting the obs_catalogue ensembles with wider bounds leaves five of the six tables numerically unchanged to four decimals. The sixth, obs / DC / width, does move:

parameter 1/ARD_LS (1e-2,1e2) (1e-2,1e4) ST (1e-2,1e2) ST (1e-2,1e4)
sigma 27.7875 30.3473 0.5679 0.6034
gamma 0.8309 0.9896 0.3981 0.4168
log_a 1.7440 1.8938 0.4366 0.4521
log_p 0.1408 0.2567 0.0054 0.0103
logp_width 0.1177 0.1741 0.0064 0.0092

AggRank shifts by at most 0.2 and the row order does not change. Note that obs / DC / width is nonetheless a fully converged fit — see §4.

The notebooks that actually contain a clamp

GSA_MC_rev.ipynb was rebuilt from source — data_revision/MC_params.txt plus MC_SED_tumbling.txt through SpyDust.SED_fit.fit_sed_ensemble, the peak- frequency mask 10 < nu_peak < 100, the first 5000 survivors, and log applied to the a column — reproducing the notebook’s (5000, 9) input with xh dropped as constant, exactly as published.

Both of its tables carry a clamped parameter (y, at ARD_LS = 100.000):

Naming. GSA_MC_rev.ipynb and obs_catalogue…ipynb both produce tables labelled “MC / width”, from different datasets: the former is 5000 SED fits over 9 environment parameters, the latter is a 5000 x 7 distribution-ensemble. They are written MC_rev / … and obs / MC / … below, because their conclusions differ.

All six _rev tables were rebuilt and re-fitted. MC_rev and DC_rev drop xh as constant; HII_rev drops y.

table LML published LML, upper bound 1e4 LML, restarts=5 AggRank order vs upper 1e4
MC_rev / peak 5909.4 5941.2 5909.4 unchanged
MC_rev / width 6225.8 6237.1 6225.8 unchanged
DC_rev / peak 11908.1 12046.1 11908.1 unchanged
DC_rev / width 15818.3 15996.7 15818.3 changed
HII_rev / peak 6594.0 6728.4 6594.0 changed
HII_rev / width 4928.3 4978.0 4928.3 changed

restarts=5 reproduces the published likelihood exactly in all six. These notebooks are converged; the optimiser defect of §4 does not touch them.

The upper bound does change the reported AggRank order in three of six. Every such change is an adjacent swap between two parameters at the bottom of the ranking whose AggRank differs by at most 0.2:

table swap published AggRank variant AggRank their ST
DC_rev / width y ↔ Chi (7th, 8th) 7.4, 7.6 7.5, 7.5 (tied) 0.0000, 0.0000
HII_rev / peak Chi ↔ beta (4th, 5th) 4.8, 5.0 5.0, 4.8 0.0284, 0.0281
HII_rev / width T ↔ xh (7th, 8th) 7.1, 7.1 (already tied) 7.1, 7.0 0.0010, 0.0000

The mechanism is tie-breaking, not re-ranking:

  • In DC_rev / width, y and Chi are both clamped at 1/ARD_LS = 0.0100, so they tie exactly on that column. Relaxing the bound separates them (0.0001 vs 0.0022) and the tie resolves the other way. Under the variant their AggRank becomes exactly equal, so the printed order is then just the stable sort.
  • In HII_rev / width, xh is clamped and T is not, and their AggRank is already exactly tied at 7.1 in the published table — the published order between them was never determined by the data.
  • HII_rev / peak involves no clamped parameter at all (Chi and beta sit at 0.399 and 0.301). Their ST differs by 0.0003; the refit flips that ordering. This one is not a clamp effect but two indistinguishable parameters trading places.

In none of the three does a parameter move in or out of the important group. The top three rows are identical in every table, and every swapped parameter has ST ≤ 0.03 with S1, MI and dCor all ≈ 0.

What the clamp does mean is that 1/ARD_LS = 0.010 on those rows is a property of the bound, not a measurement — and, as DC_rev / width shows, it manufactures exact ties that then decide the printed row order. The conclusion drawn from those rows — that the parameters do not matter — is independently supported by their S1, ST, MI and dCor, so it stands. The AggRank position of the bottom two or three rows should not be read as meaningful.

3. Item 2 — how irreproducible are the published S1_conf / ST_conf?

Substantially. Roughly a quarter of their own value.

SALib guards its bootstrap seeding with if seed:. MCPost passes seed=gp_random_state, and the notebooks all use gp_random_state=0, which is falsy — so the confidence-interval bootstrap ran unseeded in every published run.

First, the scope. Re-running sobol.analyze on a fixed surrogate response with seeds {0, 1, 2, 42, 12345, None}:

quantity spread across seeds
S1 point estimates 0.000e+00
ST point estimates 0.000e+00
S1_conf 9.6e-03
ST_conf 7.9e-03

So the defect touches the two *_conf columns and nothing else. S1, ST, MI, dCor, PermMean, ARD_LS and AggRank are untouched by this defect (but see §4 for a different one that does touch them).

Second, the magnitude. Fifteen repeat analyses of the same surrogate — so every difference is purely the unseeded bootstrap:

table S1_conf relative spread ST_conf relative spread
obs / MC / nu_peak 16 – 37 % 19 – 35 %
obs / MC / width 18 – 29 % 20 – 35 %
obs / DC / nu_peak 16 – 26 % 23 – 47 %
obs / DC / width 22 – 31 % 13 – 36 %
obs / HII / nu_peak 21 – 35 % 19 – 37 %
obs / HII / width 21 – 42 % 15 – 31 %

Typical spread ~25 % of the value. For several parameters the published value falls outside the range spanned by 15 re-runs — obs / MC / nu_peak’s gamma has ST_conf = 0.0336 published against an observed range of [0.0346, 0.0447].

Consequence. The published S1_conf and ST_conf are single draws from a distribution about a quarter as wide as the value itself. Quoting them to three decimals asserts precision that is not there. They should either be regenerated with a seeded bootstrap, or reported with an explicit statement that they are Monte-Carlo confidence estimates with ~25 % run-to-run variability. The scientific conclusions do not depend on them — no claim in these tables rests on a *_conf value.

4. A defect that does touch the published point estimates

§3 established that the unseeded bootstrap moves only S1_conf / ST_conf. A separate defect — the GP hyperparameter fit landing in a non-converged optimum — does reach S1, ST and ARD_LS, and it does so for some of these tables despite scaler="standard".

The synthetic n x d grid never shows a collapse under scaler="standard", so the initial read was that the paper’s configuration was entirely safe. On the paper’s actual data that is not uniformly true. Re-fitting each table with n_restarts_optimizer=5 and comparing log marginal likelihood:

table LML, published (restarts=0) LML, restarts=5 gain largest ST change
obs / MC / nu_peak 375.1 375.1 +0.0 0.0000
obs / MC / width −2860.3 −1904.3 +956 0.214
obs / DC / nu_peak −3635.0 −3436.7 +198 0.015
obs / DC / width −902.2 −902.2 +0.0 0.0000
obs / HII / nu_peak −200.9 −200.9 +0.0 0.0000
obs / HII / width −2875.8 −2804.6 +71 0.008

Three of the six were already at the optimum and are bit-identical. Two gain a modest amount of likelihood and move their indices by less than 0.015, which changes no ordering. One — obs / MC / width — is a genuine outlier.

“Not collapsed” is not the same as “converged”. obs / MC / width sits 956 nats below the optimum, and at the better optimum its total-effect indices move substantially:

parameter ST published ST at the better fit
gamma 0.6737 0.4600
sigma 0.2633 0.4671
logp_width 0.3041 0.2005
log_a 0.1716 0.2280

That reorders the top three: published gamma > logp_width > sigma, at the better fit sigma > gamma > logp_width. obs / MC / nu_peak by contrast is bit-for-bit identical — its fit was already at the optimum.

obs / DC / nu_peak (+198 nats, ST ≤0.015) and obs / HII / width (+71 nats, ST ≤0.008) move too little to change any ordering.

The lesson for future audits: a surrogate can pass every range invariant, avoid collapse entirely, and still be sitting in a bad local optimum. Log marginal likelihood is the diagnostic that catches it, which is why v0.2.0 reports it in extras["log_marginal_likelihood"].

5. Item 3 — obs_catalogue_vs_MonteCarlo_samples.ipynb

This notebook holds the three flagged numbers. All three are in the DC ensemble (notebook cell 12), and two of the three are in DC / width.

S1 = -0.030 (sigma, obs / DC / width)

Real, reproducible, and not usable as a sensitivity. A first-order Sobol index is a variance fraction and cannot be negative; a negative value means the estimator’s numerator is below its own noise floor. The honest reading is “sigma has no resolvable first-order effect”, i.e. S1 ≈ 0, with sigma’s influence entering through interactions — which its ST = 0.568 supports. The published -0.030 should be read as zero, not as a small negative sensitivity.

Note that v0.2.0 now warns on this: S1 < -0.05 trips the range invariant, and several parameters across these tables sit close to that threshold.

ST_conf = 0.382 and 0.333 (sigma and gamma, obs / DC / width)

Real, and the largest confidence intervals in the whole set — 67 % and 84 % of their own point estimates. Combined with §3’s finding that *_conf has ~25 % run-to-run variability, these two numbers carry essentially no information beyond “this is poorly determined”. DC / width is also the table with the largest 1/ARD_LS and the weakest fit, so the wide intervals are consistent with the rest of its diagnostics rather than an isolated anomaly.

1/ARD_LS = 27.787 (sigma, obs / DC / width)

ARD_LS = 0.036, i.e. 3.6x the lower bound of 0.01 — close to it, but not pinned, so it is a measurement rather than an artefact. It is nonetheless the signature of a GP that is fitting sigma as very short-correlation-length, which is what a surrogate does when it cannot explain a direction smoothly. Re-fitting with a relaxed upper bound moves it to 30.35 and ST from 0.568 to 0.603, so the number is sensitive to fit settings at the ~10 % level.

Did these enter the paper? They appear in the notebook’s stored LaTeX output for the DC ensemble. Whether that specific table was carried into the manuscript is a question for the authors — the notebook emits LaTeX for all three ensembles, and this review cannot tell which were typeset.

6. Summary — which numbers to trust

Quantity Verdict
Reproducibility of all six obs_catalogue tables Exact — cell for cell
MI, dCor, PermMean, PermStd Trustworthy. Untouched by every defect found.
AggRank ordering, top rows Trustworthy. The leading parameters are identical under every variant tested.
AggRank ordering, bottom 2-3 rows Do not read as ranked. In 3 of 6 _rev tables the clamp creates exact 1/ARD_LS ties that decide the printed order; relaxing the bound swaps them. All have ST ≤ 0.03.
S1, ST for obs / MC / nu_peak Trustworthy. Bit-identical at the converged fit.
S1, ST for obs / DC / nu_peak, obs / HII / width Trustworthy. Move by ≤0.015 and ≤0.008; ordering unchanged.
S1, ST for obs / DC / width, obs / HII / nu_peak Trustworthy. Bit-identical at the converged fit.
S1, ST for obs / MC / width Re-run before quoting. The published fit is 956 nats below the optimum and the top-three ST ordering changes.
S1, ST, AggRank for MC_rev, DC_rev / peak Trustworthy. Fits already converged (restarts=5 reproduces the published likelihood exactly) and insensitive to the ARD upper bound.
S1_conf, ST_conf, all tables Do not quote to three decimals. Unseeded bootstrap; ~25 % run-to-run spread. Re-report with a seed, or state the variability.
S1 = -0.030 (obs / DC / width, sigma) Report as ~0. Negative variance fractions are estimator noise.
ST_conf = 0.382 / 0.333 (obs / DC / width) Real but uninformative. 67 % and 84 % of their point estimates.
1/ARD_LS = 0.010 rows (main six notebooks) Boundary artefact, not a measurement. ARD_LS = 100 = the bound exactly. The associated “this parameter is unimportant” conclusion is independently supported by S1, ST, MI and dCor all being ~0.
1/ARD_LS = 27.787 (obs / DC / width, sigma) Measurement, but fit-sensitive at the ~10 % level.

None of the paper’s qualitative conclusions — which parameters drive the AME peak frequency and width — is affected by any defect found here. The parameter rankings are stable under every variation tested. What is affected is the precision claimed for individual numbers, and one table’s total-effect values.

7. Scope and limits of this review

What was done:

  • All three obs_catalogue ensembles re-run exactly, from the notebook’s own pre-computed inputs, and verified cell-for-cell against the published LaTeX.
  • All three _rev notebooks rebuilt end to end from data_revision/ through fit_sed_ensemble, and re-fitted under three bound/optimiser settings each (18 GP fits on 5000 x 8).
  • The bootstrap seed’s scope established directly, by re-analysing a fixed surrogate response across six seeds.
  • Cross-run dispersion of S1_conf / ST_conf measured over 15 repeats per table.

What was not done:

  • The three non-_rev notebooks (GSA_{DC,HII,MC}.ipynb, which read data/ rather than data_revision/) were not re-fitted. All six _rev tables and all six obs_catalogue tables were.
  • SED_feature_emulation.ipynb, SEDfeature_emulator.ipynb, moment_expansion.ipynb and everything under case_study/ were not examined.
  • No attempt was made to determine which of the notebooks’ stored LaTeX tables were actually typeset into the manuscript.

No file in the SpyDust repository was read-modified or re-executed. Every re-run wrote only to a scratch directory, and every fit used a pristine checkout of MCPost v0.1.4 rather than the patched working tree, so the numbers above describe the code the paper ran.