Owen (Grouped) Shapley Effects for the SAC-SMA Model¶

With RS-HDMR Surrogate Validation¶

This notebook demonstrates Owen (grouped) Shapley effects for the Sacramento Soil Moisture Accounting (SAC-SMA) model:

Phase 1 — Build an RS-HDMR surrogate of the NSE response using QMC Sobol (scrambled + shifted) training data
Phase 2 — Compute Shapley Effects under independent inputs (standard RS-HDMR)
Phase 3 — Compute Shapley Effects under correlated inputs using MC Shapley with a Gaussian copula and a literature-based correlation matrix from Vrugt et al. (2006)
Phase 4 — Compare independent vs correlated sensitivity rankings

Reference: Vrugt, J.A., et al. (2006). Application of stochastic parameter optimization to the Sacramento Soil Moisture Accounting model. Journal of Hydrology, 325, 288–307.

Catchment: St Helens Creek, North Queensland (area = 120.5 km²), daily data from Feb 1973 – Dec 2018.

In [1]:

Copied!





import numpy as np
import pandas as pd
import matplotlib
#matplotlib.use('Agg')
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
from scipy.stats.qmc import Sobol
import warnings
import sys
import os
warnings.filterwarnings('ignore')

plt.rcParams.update({
    'font.family': 'serif',
    'font.size': 11,
    'figure.dpi': 120,
    'savefig.dpi': 200,
    'savefig.bbox': 'tight'
})

# ShapleyX
from shapleyx import rshdmr
from shapleyx.utilities.mc_shapley import (
    GaussianCopulaUniform,
    shapley_effects,
)

# SAC-SMA model
sys.path.insert(0, os.path.join(os.getcwd(), 'sac_sma_data'))
from sacramento import sacsma

print('All imports OK')
import numpy as np
import pandas as pd
import matplotlib
#matplotlib.use('Agg')
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
from scipy.stats.qmc import Sobol
import warnings
import sys
import os
warnings.filterwarnings('ignore')

plt.rcParams.update({
    'font.family': 'serif',
    'font.size': 11,
    'figure.dpi': 120,
    'savefig.dpi': 200,
    'savefig.bbox': 'tight'
})

# ShapleyX
from shapleyx import rshdmr
from shapleyx.utilities.mc_shapley import (
    GaussianCopulaUniform,
    shapley_effects,
)

# SAC-SMA model
sys.path.insert(0, os.path.join(os.getcwd(), 'sac_sma_data'))
from sacramento import sacsma

print('All imports OK')

All imports OK

1. Load Forcing Data and Observations¶

We have ~46 years of daily rainfall, PET, and observed discharge for St Helens Creek (North Queensland). Catchment area: 120.5 km².

In [2]:

Copied!





df_forcing = pd.read_csv('sac_sma_data/st_helens_forcing_original.csv', parse_dates=['date'], dayfirst=True)
print(f'Forcing data: {len(df_forcing)} daily timesteps')
print(f'Date range: {df_forcing.date.iloc[0].date()} to {df_forcing.date.iloc[-1].date()}')

CATCHMENT_AREA_M2 = 120_500_000.0

prcp = df_forcing['rainfall'].values
pet  = df_forcing['pet'].values
q_obs = df_forcing['Q_CUMEC'].values

def flow_mm_to_cumec(flow_mm):
    return flow_mm * CATCHMENT_AREA_M2 / 1000.0 / 86400.0

WARMUP_DAYS = 365
print(f'Warm-up: {WARMUP_DAYS} days')
print(f'Evaluation: day {WARMUP_DAYS+1} to {len(df_forcing)} ({len(df_forcing)-WARMUP_DAYS} timesteps)')
df_forcing = pd.read_csv('sac_sma_data/st_helens_forcing_original.csv', parse_dates=['date'], dayfirst=True)
print(f'Forcing data: {len(df_forcing)} daily timesteps')
print(f'Date range: {df_forcing.date.iloc[0].date()} to {df_forcing.date.iloc[-1].date()}')

CATCHMENT_AREA_M2 = 120_500_000.0

prcp = df_forcing['rainfall'].values
pet  = df_forcing['pet'].values
q_obs = df_forcing['Q_CUMEC'].values

def flow_mm_to_cumec(flow_mm):
    return flow_mm * CATCHMENT_AREA_M2 / 1000.0 / 86400.0

WARMUP_DAYS = 365
print(f'Warm-up: {WARMUP_DAYS} days')
print(f'Evaluation: day {WARMUP_DAYS+1} to {len(df_forcing)} ({len(df_forcing)-WARMUP_DAYS} timesteps)')

Forcing data: 16763 daily timesteps
Date range: 1973-02-08 to 2018-12-31
Warm-up: 365 days
Evaluation: day 366 to 16763 (16398 timesteps)

2. Define the SAC-SMA NSE Wrapper¶

This function takes a 14-parameter vector, runs the Numba-JIT-compiled SAC-SMA model, converts the simulated flow from mm/day to m³/s, and computes the Nash-Sutcliffe Efficiency (NSE) against observed discharge.

We first run a warmup call to trigger Numba compilation.

In [3]:

Copied!





def nse_from_params(params, warmup=WARMUP_DAYS):
    flow_mm = sacsma(prcp, pet, params)
    sim = flow_mm_to_cumec(flow_mm[warmup:])
    obs = q_obs[warmup:]
    residual = obs - sim
    numerator = np.sum(residual ** 2)
    denominator = np.sum((obs - np.mean(obs)) ** 2)
    return 1.0 - numerator / denominator if denominator > 0 else -np.inf

# Warm-up (triggers Numba JIT)
params_default = np.array([44.13, 52.76, 80.87, 112.70, 52.81, 0.084, 0.490,
                           0.0148, 0.133, 35.31, 2.532, 0.00578, 0.424, 0.0093])
nse_default = nse_from_params(params_default)
print(f'NSE at default parameters: {nse_default:.4f}')
def nse_from_params(params, warmup=WARMUP_DAYS):
    flow_mm = sacsma(prcp, pet, params)
    sim = flow_mm_to_cumec(flow_mm[warmup:])
    obs = q_obs[warmup:]
    residual = obs - sim
    numerator = np.sum(residual ** 2)
    denominator = np.sum((obs - np.mean(obs)) ** 2)
    return 1.0 - numerator / denominator if denominator > 0 else -np.inf

# Warm-up (triggers Numba JIT)
params_default = np.array([44.13, 52.76, 80.87, 112.70, 52.81, 0.084, 0.490,
                           0.0148, 0.133, 35.31, 2.532, 0.00578, 0.424, 0.0093])
nse_default = nse_from_params(params_default)
print(f'NSE at default parameters: {nse_default:.4f}')

NSE at default parameters: 0.8406

3. Generate Training Data: QMC Sobol Sequence (4096 samples)¶

We use a scrambled Sobol sequence with Cranley-Patterson rotation (random shift) to produce 4096 well-spaced samples in the 14-dimensional unit hypercube. Each dimension is then mapped to the physically plausible range from es_parameters.csv.

In [4]:

Copied!





N_PARAMS = 14
N_SAMPLES = 4096

param_names = ['Uztwm', 'Uzfwm', 'Lztwm', 'Lzfpm', 'Lzfsm',
               'Adimp', 'Uzk', 'Lzpk', 'Lzsk', 'Zperc',
               'Rexp', 'Pctim', 'Pfree', 'Side']

param_lower = np.array([30, 10, 20, 40, 15, 0, 0.2, 0.001, 0.03, 0, 0, 0, 0, 0])
param_upper = np.array([125, 75, 300, 600, 300, 0.2, 0.9, 0.03, 0.8, 80, 5, 0.05, 0.8, 0.8])
param_ranges = param_upper - param_lower

print(f'Parameter space: {N_PARAMS} dimensions x {N_SAMPLES} Sobol samples')
for i, name in enumerate(param_names):
    print(f'  {name:>6s}: [{param_lower[i]:8.3f}, {param_upper[i]:8.3f}]')
N_PARAMS = 14
N_SAMPLES = 4096

param_names = ['Uztwm', 'Uzfwm', 'Lztwm', 'Lzfpm', 'Lzfsm',
               'Adimp', 'Uzk', 'Lzpk', 'Lzsk', 'Zperc',
               'Rexp', 'Pctim', 'Pfree', 'Side']

param_lower = np.array([30, 10, 20, 40, 15, 0, 0.2, 0.001, 0.03, 0, 0, 0, 0, 0])
param_upper = np.array([125, 75, 300, 600, 300, 0.2, 0.9, 0.03, 0.8, 80, 5, 0.05, 0.8, 0.8])
param_ranges = param_upper - param_lower

print(f'Parameter space: {N_PARAMS} dimensions x {N_SAMPLES} Sobol samples')
for i, name in enumerate(param_names):
    print(f'  {name:>6s}: [{param_lower[i]:8.3f}, {param_upper[i]:8.3f}]')

Parameter space: 14 dimensions x 4096 Sobol samples
   Uztwm: [  30.000,  125.000]
   Uzfwm: [  10.000,   75.000]
   Lztwm: [  20.000,  300.000]
   Lzfpm: [  40.000,  600.000]
   Lzfsm: [  15.000,  300.000]
   Adimp: [   0.000,    0.200]
     Uzk: [   0.200,    0.900]
    Lzpk: [   0.001,    0.030]
    Lzsk: [   0.030,    0.800]
   Zperc: [   0.000,   80.000]
    Rexp: [   0.000,    5.000]
   Pctim: [   0.000,    0.050]
   Pfree: [   0.000,    0.800]
    Side: [   0.000,    0.800]

In [5]:

Copied!





sobol_gen = Sobol(d=N_PARAMS, scramble=True, seed=42)
X_uniform = sobol_gen.random(n=N_SAMPLES)

# Cranley-Patterson rotation
rng = np.random.default_rng(123)
shift = rng.uniform(0, 1, size=N_PARAMS)
X_uniform = (X_uniform + shift) % 1.0

X_params = param_lower + X_uniform * param_ranges
print(f'Sobol samples shape: {X_params.shape}')
print(f'First 3 samples:')
for i in range(min(3, N_SAMPLES)):
    print(f'  [{i}] ' + ', '.join(f'{v:.2f}' for v in X_params[i]))
sobol_gen = Sobol(d=N_PARAMS, scramble=True, seed=42)
X_uniform = sobol_gen.random(n=N_SAMPLES)

# Cranley-Patterson rotation
rng = np.random.default_rng(123)
shift = rng.uniform(0, 1, size=N_PARAMS)
X_uniform = (X_uniform + shift) % 1.0

X_params = param_lower + X_uniform * param_ranges
print(f'Sobol samples shape: {X_params.shape}')
print(f'First 3 samples:')
for i in range(min(3, N_SAMPLES)):
    print(f'  [{i}] ' + ', '.join(f'{v:.2f}' for v in X_params[i]))

Sobol samples shape: (4096, 14)
First 3 samples:
  [0] 40.77, 66.43, 27.50, 173.18, 150.05, 0.13, 0.63, 0.03, 0.51, 0.36, 0.21, 0.05, 0.16, 0.65
  [1] 65.53, 45.54, 132.44, 68.88, 36.73, 0.04, 0.29, 0.02, 0.13, 57.01, 4.46, 0.03, 0.45, 0.32
  [2] 72.51, 49.16, 232.04, 310.12, 257.91, 0.17, 0.72, 0.02, 0.03, 39.35, 3.77, 0.02, 0.80, 0.01

4. Evaluate SAC-SMA for All Training Samples¶

This runs the model 4096 times (~30–60 seconds with Numba JIT). Each run simulates ~46 years of daily streamflow.

In [6]:

Copied!





print('Evaluating 4096 parameter sets...')
Y_nse = np.array([nse_from_params(X_params[i]) for i in range(N_SAMPLES)])

valid = Y_nse > -np.inf
print(f'NSE distribution ({valid.sum()} valid samples):')
Yv = Y_nse[valid]
print(f'  Mean:   {Yv.mean():.4f}')
print(f'  Std:    {Yv.std():.4f}')
print(f'  Min:    {Yv.min():.4f}')
print(f'  Max:    {Yv.max():.4f}')
print(f'  Median: {np.median(Yv):.4f}')
print(f'  NSE > 0: {(Yv > 0).mean()*100:.1f}% (better than mean-flow benchmark)')
print('Evaluating 4096 parameter sets...')
Y_nse = np.array([nse_from_params(X_params[i]) for i in range(N_SAMPLES)])

valid = Y_nse > -np.inf
print(f'NSE distribution ({valid.sum()} valid samples):')
Yv = Y_nse[valid]
print(f'  Mean:   {Yv.mean():.4f}')
print(f'  Std:    {Yv.std():.4f}')
print(f'  Min:    {Yv.min():.4f}')
print(f'  Max:    {Yv.max():.4f}')
print(f'  Median: {np.median(Yv):.4f}')
print(f'  NSE > 0: {(Yv > 0).mean()*100:.1f}% (better than mean-flow benchmark)')

Evaluating 4096 parameter sets...
NSE distribution (4096 valid samples):
  Mean:   0.6927
  Std:    0.0977
  Min:    0.3339
  Max:    0.8580
  Median: 0.7029
  NSE > 0: 100.0% (better than mean-flow benchmark)

In [7]:

Copied!





fig, ax = plt.subplots(figsize=(8, 4))
ax.hist(Yv, bins=40, color='steelblue', edgecolor='white', alpha=0.8)
ax.axvline(Yv.mean(), color='firebrick', ls='--', lw=1.5, label=f'Mean = {Yv.mean():.3f}')
ax.axvline(0, color='grey', ls=':', lw=1, label='NSE = 0')
ax.set_xlabel('NSE')
ax.set_ylabel('Frequency')
ax.set_title('Distribution of NSE Across 4096 Sobol Samples')
ax.legend(fontsize=9)
plt.tight_layout()
plt.savefig('nse_distribution.eps')
plt.savefig('nse_distribution.svg')
plt.show()
fig, ax = plt.subplots(figsize=(8, 4))
ax.hist(Yv, bins=40, color='steelblue', edgecolor='white', alpha=0.8)
ax.axvline(Yv.mean(), color='firebrick', ls='--', lw=1.5, label=f'Mean = {Yv.mean():.3f}')
ax.axvline(0, color='grey', ls=':', lw=1, label='NSE = 0')
ax.set_xlabel('NSE')
ax.set_ylabel('Frequency')
ax.set_title('Distribution of NSE Across 4096 Sobol Samples')
ax.legend(fontsize=9)
plt.tight_layout()
plt.savefig('nse_distribution.eps')
plt.savefig('nse_distribution.svg')
plt.show()

The PostScript backend does not support transparency; partially transparent artists will be rendered opaque.

No description has been provided for this image

5. Build the RS-HDMR Surrogate Model¶

We train the surrogate using:

polys=[6] — full 5th-order polynomial chaos expansion, capturing interactions up to order 6
method='omp_cv_stream' — OMP with CV (suitable for ~30K basis functions)

In [8]:

Copied!





df_train = pd.DataFrame(X_params, columns=param_names)
df_train['Y'] = Y_nse

print('Training RS-HDMR surrogate (streaming OMP-CV)...')
analyzer = rshdmr(
    data_file=df_train,
    polys=[6],
    method='omp_cv_stream',
    n_iter=400,
    test_size=0.25,
    verbose=True,
    resampling=False,
    n_jobs=4,
)

sobol, shap, total = analyzer.run_all()
print('\nSurrogate training complete.')
df_train = pd.DataFrame(X_params, columns=param_names)
df_train['Y'] = Y_nse

print('Training RS-HDMR surrogate (streaming OMP-CV)...')
analyzer = rshdmr(
    data_file=df_train,
    polys=[6],
    method='omp_cv_stream',
    n_iter=400,
    test_size=0.25,
    verbose=True,
    resampling=False,
    n_jobs=4,
)

sobol, shap, total = analyzer.run_all()
print('\nSurrogate training complete.')

Training RS-HDMR surrogate (streaming OMP-CV)...
Found a DataFrame

============================================================
Transforming data to unit hypercube
============================================================

Feature: Uztwm, Min Value: 30.0162, Max Value: 124.9915
Feature: Uzfwm, Min Value: 10.0222, Max Value: 74.9941
Feature: Lztwm, Min Value: 20.0122, Max Value: 299.9523
Feature: Lzfpm, Min Value: 40.1022, Max Value: 599.9162
Feature: Lzfsm, Min Value: 15.1039, Max Value: 299.9918
Feature: Adimp, Min Value: 0.0000, Max Value: 0.1999
Feature: Uzk, Min Value: 0.2000, Max Value: 0.8999
Feature: Lzpk, Min Value: 0.0010, Max Value: 0.0300
Feature: Lzsk, Min Value: 0.0302, Max Value: 0.8000
Feature: Zperc, Min Value: 0.0116, Max Value: 79.9933
Feature: Rexp, Min Value: 0.0008, Max Value: 4.9996
Feature: Pctim, Min Value: 0.0000, Max Value: 0.0500
Feature: Pfree, Min Value: 0.0001, Max Value: 0.7999
Feature: Side, Min Value: 0.0001, Max Value: 0.8000

============================================================
Building basis functions
============================================================

Total basis functions in basis set : 38759
Total number of features in basis set is 38759

============================================================
Running regression analysis
============================================================

running Streaming OMP-CV
  Evaluating 10 folds (4 workers)...
  Fold 1/10  done
  Fold 2/10  done
  Fold 3/10  done
  Fold 4/10  done
  Fold 5/10  done
  Fold 6/10  done
  Fold 7/10  done
  Fold 8/10  done
  Fold 9/10  done
  Fold 10/10  done
  All 10 folds complete.  Final CV path:
  Best: nz=210 (score=0.9720)
  CV nz=1   /400  score=0.2214 ±0.0276
  CV nz=2   /400  score=0.3983 ±0.0354
  CV nz=3   /400  score=0.5316 ±0.0234
  CV nz=4   /400  score=0.6040 ±0.0190
  CV nz=5   /400  score=0.6454 ±0.0211
  CV nz=6   /400  score=0.6771 ±0.0191
  CV nz=7   /400  score=0.7182 ±0.0149
  CV nz=8   /400  score=0.7522 ±0.0121
  CV nz=9   /400  score=0.7699 ±0.0137
  CV nz=10  /400  score=0.7847 ±0.0143
  CV nz=11  /400  score=0.7945 ±0.0131
  CV nz=12  /400  score=0.8078 ±0.0128
  CV nz=13  /400  score=0.8152 ±0.0116
  CV nz=14  /400  score=0.8241 ±0.0132
  CV nz=15  /400  score=0.8305 ±0.0136
  CV nz=16  /400  score=0.8383 ±0.0134
  CV nz=17  /400  score=0.8471 ±0.0126
  CV nz=18  /400  score=0.8529 ±0.0123
  CV nz=19  /400  score=0.8616 ±0.0106
  CV nz=20  /400  score=0.8644 ±0.0100
  CV nz=21  /400  score=0.8690 ±0.0100
  CV nz=22  /400  score=0.8746 ±0.0106
  CV nz=23  /400  score=0.8810 ±0.0103
  CV nz=24  /400  score=0.8859 ±0.0109
  CV nz=25  /400  score=0.8908 ±0.0106
  CV nz=26  /400  score=0.8966 ±0.0104
  CV nz=27  /400  score=0.8989 ±0.0099
  CV nz=28  /400  score=0.9011 ±0.0096
  CV nz=29  /400  score=0.9041 ±0.0092
  CV nz=30  /400  score=0.9091 ±0.0085
  CV nz=31  /400  score=0.9107 ±0.0088
  CV nz=32  /400  score=0.9121 ±0.0093
  CV nz=33  /400  score=0.9134 ±0.0093
  CV nz=34  /400  score=0.9161 ±0.0087
  CV nz=35  /400  score=0.9194 ±0.0087
  CV nz=36  /400  score=0.9195 ±0.0088
  CV nz=37  /400  score=0.9195 ±0.0091
  CV nz=38  /400  score=0.9209 ±0.0093
  CV nz=39  /400  score=0.9221 ±0.0099
  CV nz=40  /400  score=0.9232 ±0.0097
  CV nz=41  /400  score=0.9238 ±0.0097
  CV nz=42  /400  score=0.9253 ±0.0100
  CV nz=43  /400  score=0.9268 ±0.0093
  CV nz=44  /400  score=0.9277 ±0.0091
  CV nz=45  /400  score=0.9291 ±0.0088
  CV nz=46  /400  score=0.9298 ±0.0079
  CV nz=47  /400  score=0.9311 ±0.0074
  CV nz=48  /400  score=0.9322 ±0.0074
  CV nz=49  /400  score=0.9333 ±0.0071
  CV nz=50  /400  score=0.9347 ±0.0065
  CV nz=51  /400  score=0.9358 ±0.0064
  CV nz=52  /400  score=0.9369 ±0.0060
  CV nz=53  /400  score=0.9380 ±0.0064
  CV nz=54  /400  score=0.9392 ±0.0065
  CV nz=55  /400  score=0.9407 ±0.0063
  CV nz=56  /400  score=0.9420 ±0.0062
  CV nz=57  /400  score=0.9430 ±0.0061
  CV nz=58  /400  score=0.9435 ±0.0059
  CV nz=59  /400  score=0.9440 ±0.0060
  CV nz=60  /400  score=0.9444 ±0.0061
  CV nz=61  /400  score=0.9452 ±0.0058
  CV nz=62  /400  score=0.9457 ±0.0062
  CV nz=63  /400  score=0.9462 ±0.0061
  CV nz=64  /400  score=0.9465 ±0.0059
  CV nz=65  /400  score=0.9471 ±0.0061
  CV nz=66  /400  score=0.9477 ±0.0062
  CV nz=67  /400  score=0.9484 ±0.0060
  CV nz=68  /400  score=0.9495 ±0.0058
  CV nz=69  /400  score=0.9500 ±0.0056
  CV nz=70  /400  score=0.9510 ±0.0051
  CV nz=71  /400  score=0.9516 ±0.0050
  CV nz=72  /400  score=0.9524 ±0.0050
  CV nz=73  /400  score=0.9531 ±0.0052
  CV nz=74  /400  score=0.9540 ±0.0051
  CV nz=75  /400  score=0.9550 ±0.0045
  CV nz=76  /400  score=0.9562 ±0.0043
  CV nz=77  /400  score=0.9574 ±0.0041
  CV nz=78  /400  score=0.9577 ±0.0039
  CV nz=79  /400  score=0.9577 ±0.0039
  CV nz=80  /400  score=0.9581 ±0.0037
  CV nz=81  /400  score=0.9582 ±0.0034
  CV nz=82  /400  score=0.9582 ±0.0034
  CV nz=83  /400  score=0.9582 ±0.0034
  CV nz=84  /400  score=0.9584 ±0.0031
  CV nz=85  /400  score=0.9585 ±0.0031
  CV nz=86  /400  score=0.9587 ±0.0031
  CV nz=87  /400  score=0.9593 ±0.0030
  CV nz=88  /400  score=0.9594 ±0.0028
  CV nz=89  /400  score=0.9597 ±0.0026
  CV nz=90  /400  score=0.9599 ±0.0026
  CV nz=91  /400  score=0.9604 ±0.0024
  CV nz=92  /400  score=0.9607 ±0.0025
  CV nz=93  /400  score=0.9609 ±0.0026
  CV nz=94  /400  score=0.9611 ±0.0027
  CV nz=95  /400  score=0.9614 ±0.0026
  CV nz=96  /400  score=0.9614 ±0.0026
  CV nz=97  /400  score=0.9616 ±0.0026
  CV nz=98  /400  score=0.9620 ±0.0029
  CV nz=99  /400  score=0.9620 ±0.0031
  CV nz=100 /400  score=0.9623 ±0.0029
  CV nz=101 /400  score=0.9625 ±0.0029
  CV nz=102 /400  score=0.9627 ±0.0030
  CV nz=103 /400  score=0.9627 ±0.0031
  CV nz=104 /400  score=0.9629 ±0.0030
  CV nz=105 /400  score=0.9630 ±0.0028
  CV nz=106 /400  score=0.9634 ±0.0028
  CV nz=107 /400  score=0.9634 ±0.0027
  CV nz=108 /400  score=0.9634 ±0.0026
  CV nz=109 /400  score=0.9634 ±0.0023
  CV nz=110 /400  score=0.9634 ±0.0023
  CV nz=111 /400  score=0.9635 ±0.0022
  CV nz=112 /400  score=0.9636 ±0.0022
  CV nz=113 /400  score=0.9635 ±0.0026
  CV nz=114 /400  score=0.9638 ±0.0026
  CV nz=115 /400  score=0.9637 ±0.0026
  CV nz=116 /400  score=0.9640 ±0.0026
  CV nz=117 /400  score=0.9641 ±0.0024
  CV nz=118 /400  score=0.9641 ±0.0028
  CV nz=119 /400  score=0.9644 ±0.0028
  CV nz=120 /400  score=0.9645 ±0.0028
  CV nz=121 /400  score=0.9645 ±0.0027
  CV nz=122 /400  score=0.9647 ±0.0027
  CV nz=123 /400  score=0.9647 ±0.0028
  CV nz=124 /400  score=0.9648 ±0.0029
  CV nz=125 /400  score=0.9650 ±0.0030
  CV nz=126 /400  score=0.9650 ±0.0030
  CV nz=127 /400  score=0.9651 ±0.0028
  CV nz=128 /400  score=0.9653 ±0.0029
  CV nz=129 /400  score=0.9652 ±0.0029
  CV nz=130 /400  score=0.9652 ±0.0031
  CV nz=131 /400  score=0.9652 ±0.0031
  CV nz=132 /400  score=0.9653 ±0.0029
  CV nz=133 /400  score=0.9653 ±0.0028
  CV nz=134 /400  score=0.9655 ±0.0029
  CV nz=135 /400  score=0.9656 ±0.0030
  CV nz=136 /400  score=0.9657 ±0.0031
  CV nz=137 /400  score=0.9659 ±0.0031
  CV nz=138 /400  score=0.9660 ±0.0031
  CV nz=139 /400  score=0.9662 ±0.0028
  CV nz=140 /400  score=0.9665 ±0.0028
  CV nz=141 /400  score=0.9665 ±0.0028
  CV nz=142 /400  score=0.9667 ±0.0029
  CV nz=143 /400  score=0.9668 ±0.0029
  CV nz=144 /400  score=0.9670 ±0.0030
  CV nz=145 /400  score=0.9672 ±0.0031
  CV nz=146 /400  score=0.9673 ±0.0032
  CV nz=147 /400  score=0.9676 ±0.0032
  CV nz=148 /400  score=0.9676 ±0.0032
  CV nz=149 /400  score=0.9678 ±0.0032
  CV nz=150 /400  score=0.9679 ±0.0031
  CV nz=151 /400  score=0.9680 ±0.0029
  CV nz=152 /400  score=0.9681 ±0.0029
  CV nz=153 /400  score=0.9682 ±0.0028
  CV nz=154 /400  score=0.9685 ±0.0028
  CV nz=155 /400  score=0.9687 ±0.0027
  CV nz=156 /400  score=0.9688 ±0.0027
  CV nz=157 /400  score=0.9692 ±0.0026
  CV nz=158 /400  score=0.9693 ±0.0025
  CV nz=159 /400  score=0.9692 ±0.0026
  CV nz=160 /400  score=0.9695 ±0.0026
  CV nz=161 /400  score=0.9698 ±0.0025
  CV nz=162 /400  score=0.9699 ±0.0028
  CV nz=163 /400  score=0.9701 ±0.0028
  CV nz=164 /400  score=0.9703 ±0.0029
  CV nz=165 /400  score=0.9705 ±0.0028
  CV nz=166 /400  score=0.9705 ±0.0028
  CV nz=167 /400  score=0.9705 ±0.0029
  CV nz=168 /400  score=0.9706 ±0.0029
  CV nz=169 /400  score=0.9707 ±0.0030
  CV nz=170 /400  score=0.9710 ±0.0027
  CV nz=171 /400  score=0.9710 ±0.0027
  CV nz=172 /400  score=0.9711 ±0.0026
  CV nz=173 /400  score=0.9712 ±0.0024
  CV nz=174 /400  score=0.9711 ±0.0024
  CV nz=175 /400  score=0.9710 ±0.0025
  CV nz=176 /400  score=0.9711 ±0.0026
  CV nz=177 /400  score=0.9710 ±0.0026
  CV nz=178 /400  score=0.9711 ±0.0026
  CV nz=179 /400  score=0.9712 ±0.0028
  CV nz=180 /400  score=0.9712 ±0.0029
  CV nz=181 /400  score=0.9713 ±0.0029
  CV nz=182 /400  score=0.9713 ±0.0029
  CV nz=183 /400  score=0.9713 ±0.0029
  CV nz=184 /400  score=0.9714 ±0.0029
  CV nz=185 /400  score=0.9714 ±0.0031
  CV nz=186 /400  score=0.9715 ±0.0033
  CV nz=187 /400  score=0.9714 ±0.0034
  CV nz=188 /400  score=0.9714 ±0.0034
  CV nz=189 /400  score=0.9715 ±0.0034
  CV nz=190 /400  score=0.9716 ±0.0034
  CV nz=191 /400  score=0.9715 ±0.0035
  CV nz=192 /400  score=0.9715 ±0.0036
  CV nz=193 /400  score=0.9714 ±0.0036
  CV nz=194 /400  score=0.9713 ±0.0034
  CV nz=195 /400  score=0.9713 ±0.0033
  CV nz=196 /400  score=0.9714 ±0.0032
  CV nz=197 /400  score=0.9714 ±0.0032
  CV nz=198 /400  score=0.9715 ±0.0033
  CV nz=199 /400  score=0.9715 ±0.0031
  CV nz=200 /400  score=0.9715 ±0.0030
  CV nz=201 /400  score=0.9715 ±0.0030
  CV nz=202 /400  score=0.9713 ±0.0030
  CV nz=203 /400  score=0.9715 ±0.0029
  CV nz=204 /400  score=0.9716 ±0.0029
  CV nz=205 /400  score=0.9717 ±0.0029
  CV nz=206 /400  score=0.9718 ±0.0029
  CV nz=207 /400  score=0.9717 ±0.0028
  CV nz=208 /400  score=0.9717 ±0.0029
  CV nz=209 /400  score=0.9719 ±0.0029
  CV nz=210 /400  score=0.9720 ±0.0029 ★
  CV nz=211 /400  score=0.9720 ±0.0029
  CV nz=212 /400  score=0.9719 ±0.0030
  CV nz=213 /400  score=0.9719 ±0.0030
  CV nz=214 /400  score=0.9719 ±0.0030
  CV nz=215 /400  score=0.9718 ±0.0030
  CV nz=216 /400  score=0.9718 ±0.0030
  CV nz=217 /400  score=0.9718 ±0.0030
  CV nz=218 /400  score=0.9716 ±0.0029
  CV nz=219 /400  score=0.9715 ±0.0030
  CV nz=220 /400  score=0.9714 ±0.0031
  CV nz=221 /400  score=0.9714 ±0.0032
  CV nz=222 /400  score=0.9714 ±0.0031
  CV nz=223 /400  score=0.9712 ±0.0030
  CV nz=224 /400  score=0.9711 ±0.0031
  CV nz=225 /400  score=0.9711 ±0.0031
  CV nz=226 /400  score=0.9712 ±0.0032
  CV nz=227 /400  score=0.9711 ±0.0031
  CV nz=228 /400  score=0.9711 ±0.0032
  CV nz=229 /400  score=0.9712 ±0.0031
  CV nz=230 /400  score=0.9711 ±0.0030
  CV nz=231 /400  score=0.9712 ±0.0030
  CV nz=232 /400  score=0.9712 ±0.0030
  CV nz=233 /400  score=0.9711 ±0.0030
  CV nz=234 /400  score=0.9710 ±0.0030
  CV nz=235 /400  score=0.9710 ±0.0030
  CV nz=236 /400  score=0.9710 ±0.0030
  CV nz=237 /400  score=0.9711 ±0.0031
  CV nz=238 /400  score=0.9711 ±0.0031
  CV nz=239 /400  score=0.9711 ±0.0031
  CV nz=240 /400  score=0.9710 ±0.0032
  CV nz=241 /400  score=0.9710 ±0.0033
  CV nz=242 /400  score=0.9710 ±0.0033
  CV nz=243 /400  score=0.9710 ±0.0033
  CV nz=244 /400  score=0.9711 ±0.0033
  CV nz=245 /400  score=0.9710 ±0.0033
  CV nz=246 /400  score=0.9712 ±0.0032
  CV nz=247 /400  score=0.9711 ±0.0032
  CV nz=248 /400  score=0.9710 ±0.0032
  CV nz=249 /400  score=0.9709 ±0.0032
  CV nz=250 /400  score=0.9709 ±0.0032
  CV nz=251 /400  score=0.9709 ±0.0032
  CV nz=252 /400  score=0.9710 ±0.0033
  CV nz=253 /400  score=0.9710 ±0.0033
  CV nz=254 /400  score=0.9709 ±0.0034
  CV nz=255 /400  score=0.9710 ±0.0033
  CV nz=256 /400  score=0.9710 ±0.0033
  CV nz=257 /400  score=0.9709 ±0.0033
  CV nz=258 /400  score=0.9709 ±0.0034
  CV nz=259 /400  score=0.9708 ±0.0034
  CV nz=260 /400  score=0.9708 ±0.0034
  CV nz=261 /400  score=0.9708 ±0.0034
  CV nz=262 /400  score=0.9706 ±0.0035
  CV nz=263 /400  score=0.9705 ±0.0037
  CV nz=264 /400  score=0.9706 ±0.0036
  CV nz=265 /400  score=0.9707 ±0.0036
  CV nz=266 /400  score=0.9707 ±0.0036
  CV nz=267 /400  score=0.9708 ±0.0036
  CV nz=268 /400  score=0.9707 ±0.0036
  CV nz=269 /400  score=0.9706 ±0.0036
  CV nz=270 /400  score=0.9707 ±0.0037
  CV nz=271 /400  score=0.9706 ±0.0036
  CV nz=272 /400  score=0.9707 ±0.0035
  CV nz=273 /400  score=0.9707 ±0.0035
  CV nz=274 /400  score=0.9707 ±0.0034
  CV nz=275 /400  score=0.9707 ±0.0034
  CV nz=276 /400  score=0.9706 ±0.0035
  CV nz=277 /400  score=0.9704 ±0.0036
  CV nz=278 /400  score=0.9704 ±0.0037
  CV nz=279 /400  score=0.9703 ±0.0037
  CV nz=280 /400  score=0.9703 ±0.0036
  CV nz=281 /400  score=0.9703 ±0.0036
  CV nz=282 /400  score=0.9703 ±0.0036
  CV nz=283 /400  score=0.9702 ±0.0036
  CV nz=284 /400  score=0.9701 ±0.0035
  CV nz=285 /400  score=0.9701 ±0.0035
  CV nz=286 /400  score=0.9701 ±0.0035
  CV nz=287 /400  score=0.9701 ±0.0035
  CV nz=288 /400  score=0.9701 ±0.0035
  CV nz=289 /400  score=0.9701 ±0.0035
  CV nz=290 /400  score=0.9700 ±0.0035
  CV nz=291 /400  score=0.9700 ±0.0034
  CV nz=292 /400  score=0.9700 ±0.0034
  CV nz=293 /400  score=0.9700 ±0.0033
  CV nz=294 /400  score=0.9698 ±0.0032
  CV nz=295 /400  score=0.9699 ±0.0033
  CV nz=296 /400  score=0.9699 ±0.0034
  CV nz=297 /400  score=0.9699 ±0.0034
  CV nz=298 /400  score=0.9699 ±0.0033
  CV nz=299 /400  score=0.9699 ±0.0033
  CV nz=300 /400  score=0.9700 ±0.0033
  CV nz=301 /400  score=0.9700 ±0.0033
  CV nz=302 /400  score=0.9700 ±0.0033
  CV nz=303 /400  score=0.9699 ±0.0033
  CV nz=304 /400  score=0.9699 ±0.0032
  CV nz=305 /400  score=0.9699 ±0.0033
  CV nz=306 /400  score=0.9699 ±0.0033
  CV nz=307 /400  score=0.9698 ±0.0033
  CV nz=308 /400  score=0.9698 ±0.0033
  CV nz=309 /400  score=0.9698 ±0.0032
  CV nz=310 /400  score=0.9698 ±0.0033
  CV nz=311 /400  score=0.9697 ±0.0033
  CV nz=312 /400  score=0.9697 ±0.0034
  CV nz=313 /400  score=0.9697 ±0.0034
  CV nz=314 /400  score=0.9697 ±0.0034
  CV nz=315 /400  score=0.9697 ±0.0034
  CV nz=316 /400  score=0.9697 ±0.0034
  CV nz=317 /400  score=0.9695 ±0.0036
  CV nz=318 /400  score=0.9695 ±0.0035
  CV nz=319 /400  score=0.9695 ±0.0035
  CV nz=320 /400  score=0.9695 ±0.0034
  CV nz=321 /400  score=0.9695 ±0.0034
  CV nz=322 /400  score=0.9695 ±0.0036
  CV nz=323 /400  score=0.9694 ±0.0036
  CV nz=324 /400  score=0.9695 ±0.0037
  CV nz=325 /400  score=0.9694 ±0.0037
  CV nz=326 /400  score=0.9693 ±0.0038
  CV nz=327 /400  score=0.9693 ±0.0038
  CV nz=328 /400  score=0.9693 ±0.0038
  CV nz=329 /400  score=0.9693 ±0.0039
  CV nz=330 /400  score=0.9692 ±0.0038
  CV nz=331 /400  score=0.9692 ±0.0039
  CV nz=332 /400  score=0.9692 ±0.0039
  CV nz=333 /400  score=0.9692 ±0.0039
  CV nz=334 /400  score=0.9692 ±0.0038
  CV nz=335 /400  score=0.9691 ±0.0039
  CV nz=336 /400  score=0.9690 ±0.0039
  CV nz=337 /400  score=0.9690 ±0.0039
  CV nz=338 /400  score=0.9690 ±0.0039
  CV nz=339 /400  score=0.9690 ±0.0039
  CV nz=340 /400  score=0.9691 ±0.0039
  CV nz=341 /400  score=0.9690 ±0.0039
  CV nz=342 /400  score=0.9690 ±0.0039
  CV nz=343 /400  score=0.9690 ±0.0039
  CV nz=344 /400  score=0.9691 ±0.0039
  CV nz=345 /400  score=0.9691 ±0.0039
  CV nz=346 /400  score=0.9689 ±0.0040
  CV nz=347 /400  score=0.9690 ±0.0040
  CV nz=348 /400  score=0.9690 ±0.0041
  CV nz=349 /400  score=0.9690 ±0.0041
  CV nz=350 /400  score=0.9691 ±0.0041
  CV nz=351 /400  score=0.9690 ±0.0041
  CV nz=352 /400  score=0.9690 ±0.0041
  CV nz=353 /400  score=0.9690 ±0.0040
  CV nz=354 /400  score=0.9689 ±0.0040
  CV nz=355 /400  score=0.9689 ±0.0039
  CV nz=356 /400  score=0.9689 ±0.0039
  CV nz=357 /400  score=0.9689 ±0.0039
  CV nz=358 /400  score=0.9689 ±0.0039
  CV nz=359 /400  score=0.9689 ±0.0039
  CV nz=360 /400  score=0.9689 ±0.0040
  CV nz=361 /400  score=0.9688 ±0.0040
  CV nz=362 /400  score=0.9689 ±0.0040
  CV nz=363 /400  score=0.9688 ±0.0040
  CV nz=364 /400  score=0.9689 ±0.0039
  CV nz=365 /400  score=0.9689 ±0.0039
  CV nz=366 /400  score=0.9688 ±0.0038
  CV nz=367 /400  score=0.9689 ±0.0039
  CV nz=368 /400  score=0.9688 ±0.0039
  CV nz=369 /400  score=0.9688 ±0.0039
  CV nz=370 /400  score=0.9687 ±0.0039
  CV nz=371 /400  score=0.9687 ±0.0039
  CV nz=372 /400  score=0.9687 ±0.0039
  CV nz=373 /400  score=0.9688 ±0.0038
  CV nz=374 /400  score=0.9687 ±0.0037
  CV nz=375 /400  score=0.9687 ±0.0037
  CV nz=376 /400  score=0.9687 ±0.0036
  CV nz=377 /400  score=0.9687 ±0.0036
  CV nz=378 /400  score=0.9688 ±0.0036
  CV nz=379 /400  score=0.9688 ±0.0037
  CV nz=380 /400  score=0.9688 ±0.0037
  CV nz=381 /400  score=0.9687 ±0.0038
  CV nz=382 /400  score=0.9687 ±0.0037
  CV nz=383 /400  score=0.9687 ±0.0038
  CV nz=384 /400  score=0.9686 ±0.0038
  CV nz=385 /400  score=0.9686 ±0.0038
  CV nz=386 /400  score=0.9687 ±0.0037
  CV nz=387 /400  score=0.9687 ±0.0036
  CV nz=388 /400  score=0.9686 ±0.0036
  CV nz=389 /400  score=0.9686 ±0.0036
  CV nz=390 /400  score=0.9685 ±0.0036
  CV nz=391 /400  score=0.9685 ±0.0035
  CV nz=392 /400  score=0.9685 ±0.0036
  CV nz=393 /400  score=0.9686 ±0.0035
  CV nz=394 /400  score=0.9686 ±0.0035
  CV nz=395 /400  score=0.9686 ±0.0036
  CV nz=396 /400  score=0.9686 ±0.0035
  CV nz=397 /400  score=0.9686 ±0.0034
  CV nz=398 /400  score=0.9686 ±0.0034
  CV nz=399 /400  score=0.9686 ±0.0034
  CV nz=400 /400  score=0.9686 ±0.0035
Fit Execution Time : 5925.121147
Basis functions selected: 210 / 38759 (0.54%)
--
 
 Model complete 
 

============================================================
RS-HDMR model performance statistics
============================================================

variance of data        : 0.010
sum of coefficients^2   : 0.009
variance ratio          : 0.995
===============================
mae error on test set   : 0.009
mse error on test set   : 0.000
explained variance score: 0.984
===============================
slope     :  0.9841591060827801
r value   :  0.9920479353754939
r^2       :  0.98415910608278
p value   :  0.0
std error :  0.0019514102473279413

============================================================
                  Completed all analysis
                 ------------------------

I have an everyday religion that works for me. Love
yourself first, and everything else falls into line.
Lucille Ball
============================================================


Surrogate training complete.

6. RS-HDMR Shapley Effects (Independent Inputs)¶

The standard RS-HDMR Shapley effects assume independent inputs. They redistribute interaction variance equally among the participating variables. This is our baseline.

In [9]:

Copied!





display_col = 'scaled effect' if 'scaled effect' in shap.columns else 'effect'
print(f'RS-HDMR Shapley Effects — Independent Inputs:')
print(shap.to_string(float_format=lambda x: f'{x:.6f}'))

shap_sorted = shap.sort_values(display_col, ascending=True)
fig, ax = plt.subplots(figsize=(8, 5))
colors = plt.cm.Blues(np.linspace(0.4, 0.9, len(shap_sorted)))
ax.barh(range(len(shap_sorted)), shap_sorted[display_col].values, color=colors)
ax.set_yticks(range(len(shap_sorted)))
ax.set_yticklabels(shap_sorted['label'].values)
ax.set_xlabel('Shapley Effect (fraction of explained variance)')
ax.set_title('RS-HDMR Shapley Effects — Independent Inputs')
ax.axvline(0, color='grey', lw=0.8)
plt.tight_layout()
plt.savefig('shapley_independent.png')
plt.show()
display_col = 'scaled effect' if 'scaled effect' in shap.columns else 'effect'
print(f'RS-HDMR Shapley Effects — Independent Inputs:')
print(shap.to_string(float_format=lambda x: f'{x:.6f}'))

shap_sorted = shap.sort_values(display_col, ascending=True)
fig, ax = plt.subplots(figsize=(8, 5))
colors = plt.cm.Blues(np.linspace(0.4, 0.9, len(shap_sorted)))
ax.barh(range(len(shap_sorted)), shap_sorted[display_col].values, color=colors)
ax.set_yticks(range(len(shap_sorted)))
ax.set_yticklabels(shap_sorted['label'].values)
ax.set_xlabel('Shapley Effect (fraction of explained variance)')
ax.set_title('RS-HDMR Shapley Effects — Independent Inputs')
ax.axvline(0, color='grey', lw=0.8)
plt.tight_layout()
plt.savefig('shapley_independent.png')
plt.show()

RS-HDMR Shapley Effects — Independent Inputs:
    label   effect  scaled effect
0   Uztwm 0.002480       0.002521
1   Uzfwm 0.003868       0.003932
2   Lztwm 0.042437       0.043132
3   Lzfpm 0.247489       0.251538
4   Lzfsm 0.200905       0.204191
5   Adimp 0.082453       0.083802
6     Uzk 0.013440       0.013660
7    Lzpk 0.049590       0.050402
8    Lzsk 0.062339       0.063359
9   Zperc 0.063154       0.064187
10   Rexp 0.161437       0.164078
11  Pctim 0.006900       0.007013
12  Pfree 0.000656       0.000667
13   Side 0.046754       0.047519

7. Build the Parameter Correlation Matrix¶

We construct the correlation matrix from Vrugt et al. (2006, Table 2) — the full 13×13 posterior correlation matrix from MCMC (SCEM-UA) calibration of SAC-SMA on the Leaf River, Mississippi. Parameter SIDE was fixed (not calibrated), so its correlations are set to zero.

This matrix defines the dependence structure for the Gaussian copula used in MC Shapley.

In [10]:

Copied!





# Vrugt et al. (2006) Table 2 — 13x13 correlation matrix
# Vrugt order: UZTWM, UZFWM, UZK, PCTIM, ADIMP, ZPERC, REXP,
#               LZTWM, LZFSM, LZFPM, LZSK, LZPK, PFREE
corr_vrugt = np.array([
    [ 1.00, -0.02, -0.04,  0.19,  0.36, -0.01,  0.11, -0.84, -0.17,  0.28,  0.17,  0.16,  0.57],
    [-0.02,  1.00, -0.53,  0.05,  0.01,  0.15,  0.11, -0.01, -0.08,  0.22, -0.15, -0.18, -0.33],
    [-0.04, -0.53,  1.00, -0.01, -0.41, -0.13, -0.03,  0.07, -0.11, -0.01, -0.02,  0.08,  0.08],
    [ 0.19,  0.05, -0.01,  1.00, -0.16,  0.01,  0.03, -0.14,  0.01,  0.02, -0.08,  0.02, -0.04],
    [ 0.36,  0.01, -0.41, -0.16,  1.00,  0.07,  0.04, -0.21, -0.15,  0.09,  0.17,  0.15,  0.31],
    [-0.01,  0.15, -0.13,  0.01,  0.07,  1.00,  0.05,  0.06, -0.06, -0.10, -0.13,  0.02,  0.00],
    [ 0.11,  0.11, -0.03,  0.03,  0.04,  0.05,  1.00, -0.03,  0.65,  0.68,  0.24, -0.10, -0.12],
    [-0.84, -0.01,  0.07, -0.14, -0.21,  0.06, -0.03,  1.00,  0.28, -0.16, -0.14, -0.23, -0.63],
    [-0.17, -0.08, -0.11,  0.01, -0.15, -0.06,  0.65,  0.28,  1.00,  0.41, -0.33, -0.55, -0.36],
    [ 0.28,  0.22, -0.01,  0.02,  0.09, -0.10,  0.68, -0.16,  0.41,  1.00,  0.06, -0.43, -0.26],
    [ 0.17, -0.15, -0.02, -0.08,  0.17, -0.13,  0.24, -0.14, -0.33,  0.06,  1.00,  0.60,  0.41],
    [ 0.16, -0.18,  0.08,  0.02,  0.15,  0.02, -0.10, -0.23, -0.55, -0.43,  0.60,  1.00,  0.56],
    [ 0.57, -0.33,  0.08, -0.04,  0.31,  0.00, -0.12, -0.63, -0.36, -0.26,  0.41,  0.56,  1.00],
])

# Map Vrugt order -> our parameter order
# Ours: Uztwm(0), Uzfwm(1), Lztwm(2), Lzfpm(3), Lzfsm(4), Adimp(5), Uzk(6),
#       Lzpk(7),  Lzsk(8),  Zperc(9), Rexp(10), Pctim(11), Pfree(12), Side(13)
vrugt_to_our = [0, 1, 6, 11, 5, 9, 10, 2, 4, 3, 8, 7, 12]

corr_14 = np.eye(N_PARAMS)
for i, vi in enumerate(vrugt_to_our):
    for j, vj in enumerate(vrugt_to_our):
        corr_14[vi, vj] = corr_vrugt[i, j]

# Ensure symmetry and positive definiteness
corr_14 = (corr_14 + corr_14.T) / 2
eigvals = np.linalg.eigvalsh(corr_14)
if eigvals.min() < 0:
    corr_14 += np.eye(N_PARAMS) * (-eigvals.min() + 1e-8)

print('Correlation matrix:')
pd.set_option('display.precision', 3)
display(pd.DataFrame(corr_14, index=param_names, columns=param_names))

print('\nStrongest correlations:')
for i in range(N_PARAMS):
    for j in range(i+1, N_PARAMS):
        if abs(corr_14[i, j]) > 0.3:
            print(f'  {param_names[i]:>6s} <-> {param_names[j]:<6s}: rho = {corr_14[i, j]:+.3f}')
# Vrugt et al. (2006) Table 2 — 13x13 correlation matrix
# Vrugt order: UZTWM, UZFWM, UZK, PCTIM, ADIMP, ZPERC, REXP,
#               LZTWM, LZFSM, LZFPM, LZSK, LZPK, PFREE
corr_vrugt = np.array([
    [ 1.00, -0.02, -0.04,  0.19,  0.36, -0.01,  0.11, -0.84, -0.17,  0.28,  0.17,  0.16,  0.57],
    [-0.02,  1.00, -0.53,  0.05,  0.01,  0.15,  0.11, -0.01, -0.08,  0.22, -0.15, -0.18, -0.33],
    [-0.04, -0.53,  1.00, -0.01, -0.41, -0.13, -0.03,  0.07, -0.11, -0.01, -0.02,  0.08,  0.08],
    [ 0.19,  0.05, -0.01,  1.00, -0.16,  0.01,  0.03, -0.14,  0.01,  0.02, -0.08,  0.02, -0.04],
    [ 0.36,  0.01, -0.41, -0.16,  1.00,  0.07,  0.04, -0.21, -0.15,  0.09,  0.17,  0.15,  0.31],
    [-0.01,  0.15, -0.13,  0.01,  0.07,  1.00,  0.05,  0.06, -0.06, -0.10, -0.13,  0.02,  0.00],
    [ 0.11,  0.11, -0.03,  0.03,  0.04,  0.05,  1.00, -0.03,  0.65,  0.68,  0.24, -0.10, -0.12],
    [-0.84, -0.01,  0.07, -0.14, -0.21,  0.06, -0.03,  1.00,  0.28, -0.16, -0.14, -0.23, -0.63],
    [-0.17, -0.08, -0.11,  0.01, -0.15, -0.06,  0.65,  0.28,  1.00,  0.41, -0.33, -0.55, -0.36],
    [ 0.28,  0.22, -0.01,  0.02,  0.09, -0.10,  0.68, -0.16,  0.41,  1.00,  0.06, -0.43, -0.26],
    [ 0.17, -0.15, -0.02, -0.08,  0.17, -0.13,  0.24, -0.14, -0.33,  0.06,  1.00,  0.60,  0.41],
    [ 0.16, -0.18,  0.08,  0.02,  0.15,  0.02, -0.10, -0.23, -0.55, -0.43,  0.60,  1.00,  0.56],
    [ 0.57, -0.33,  0.08, -0.04,  0.31,  0.00, -0.12, -0.63, -0.36, -0.26,  0.41,  0.56,  1.00],
])

# Map Vrugt order -> our parameter order
# Ours: Uztwm(0), Uzfwm(1), Lztwm(2), Lzfpm(3), Lzfsm(4), Adimp(5), Uzk(6),
#       Lzpk(7),  Lzsk(8),  Zperc(9), Rexp(10), Pctim(11), Pfree(12), Side(13)
vrugt_to_our = [0, 1, 6, 11, 5, 9, 10, 2, 4, 3, 8, 7, 12]

corr_14 = np.eye(N_PARAMS)
for i, vi in enumerate(vrugt_to_our):
    for j, vj in enumerate(vrugt_to_our):
        corr_14[vi, vj] = corr_vrugt[i, j]

# Ensure symmetry and positive definiteness
corr_14 = (corr_14 + corr_14.T) / 2
eigvals = np.linalg.eigvalsh(corr_14)
if eigvals.min() < 0:
    corr_14 += np.eye(N_PARAMS) * (-eigvals.min() + 1e-8)

print('Correlation matrix:')
pd.set_option('display.precision', 3)
display(pd.DataFrame(corr_14, index=param_names, columns=param_names))

print('\nStrongest correlations:')
for i in range(N_PARAMS):
    for j in range(i+1, N_PARAMS):
        if abs(corr_14[i, j]) > 0.3:
            print(f'  {param_names[i]:>6s} <-> {param_names[j]:<6s}: rho = {corr_14[i, j]:+.3f}')

Correlation matrix:

	Uztwm	Uzfwm	Lztwm	Lzfpm	Lzfsm	Adimp	Uzk	Lzpk	Lzsk	Zperc	Rexp	Pctim	Pfree	Side
Uztwm	1.00	-0.02	-0.84	0.28	-0.17	0.36	-0.04	0.16	0.17	-0.01	0.11	0.19	0.57	0.0
Uzfwm	-0.02	1.00	-0.01	0.22	-0.08	0.01	-0.53	-0.18	-0.15	0.15	0.11	0.05	-0.33	0.0
Lztwm	-0.84	-0.01	1.00	-0.16	0.28	-0.21	0.07	-0.23	-0.14	0.06	-0.03	-0.14	-0.63	0.0
Lzfpm	0.28	0.22	-0.16	1.00	0.41	0.09	-0.01	-0.43	0.06	-0.10	0.68	0.02	-0.26	0.0
Lzfsm	-0.17	-0.08	0.28	0.41	1.00	-0.15	-0.11	-0.55	-0.33	-0.06	0.65	0.01	-0.36	0.0
Adimp	0.36	0.01	-0.21	0.09	-0.15	1.00	-0.41	0.15	0.17	0.07	0.04	-0.16	0.31	0.0
Uzk	-0.04	-0.53	0.07	-0.01	-0.11	-0.41	1.00	0.08	-0.02	-0.13	-0.03	-0.01	0.08	0.0
Lzpk	0.16	-0.18	-0.23	-0.43	-0.55	0.15	0.08	1.00	0.60	0.02	-0.10	0.02	0.56	0.0
Lzsk	0.17	-0.15	-0.14	0.06	-0.33	0.17	-0.02	0.60	1.00	-0.13	0.24	-0.08	0.41	0.0
Zperc	-0.01	0.15	0.06	-0.10	-0.06	0.07	-0.13	0.02	-0.13	1.00	0.05	0.01	0.00	0.0
Rexp	0.11	0.11	-0.03	0.68	0.65	0.04	-0.03	-0.10	0.24	0.05	1.00	0.03	-0.12	0.0
Pctim	0.19	0.05	-0.14	0.02	0.01	-0.16	-0.01	0.02	-0.08	0.01	0.03	1.00	-0.04	0.0
Pfree	0.57	-0.33	-0.63	-0.26	-0.36	0.31	0.08	0.56	0.41	0.00	-0.12	-0.04	1.00	0.0
Side	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	1.0

Strongest correlations:
   Uztwm <-> Lztwm : rho = -0.840
   Uztwm <-> Adimp : rho = +0.360
   Uztwm <-> Pfree : rho = +0.570
   Uzfwm <-> Uzk   : rho = -0.530
   Uzfwm <-> Pfree : rho = -0.330
   Lztwm <-> Pfree : rho = -0.630
   Lzfpm <-> Lzfsm : rho = +0.410
   Lzfpm <-> Lzpk  : rho = -0.430
   Lzfpm <-> Rexp  : rho = +0.680
   Lzfsm <-> Lzpk  : rho = -0.550
   Lzfsm <-> Lzsk  : rho = -0.330
   Lzfsm <-> Rexp  : rho = +0.650
   Lzfsm <-> Pfree : rho = -0.360
   Adimp <-> Uzk   : rho = -0.410
   Adimp <-> Pfree : rho = +0.310
    Lzpk <-> Lzsk  : rho = +0.600
    Lzpk <-> Pfree : rho = +0.560
    Lzsk <-> Pfree : rho = +0.410

In [11]:

Copied!





# Build the correlated distribution
joint_correlated = GaussianCopulaUniform(
    lows=param_lower,
    highs=param_upper,
    corr=corr_14,
)

print('GaussianCopulaUniform ready (14D, Vrugt correlations).')
# Build the correlated distribution
joint_correlated = GaussianCopulaUniform(
    lows=param_lower,
    highs=param_upper,
    corr=corr_14,
)

print('GaussianCopulaUniform ready (14D, Vrugt correlations).')

GaussianCopulaUniform ready (14D, Vrugt correlations).

In [12]:

Copied!





# SAC-SMA wrapper for direct MC Shapley
def sacsma_nse_wrapper(x):
    """SAC-SMA NSE function — expects physical parameter values."""
    return nse_from_params(x)
# SAC-SMA wrapper for direct MC Shapley
def sacsma_nse_wrapper(x):
    """SAC-SMA NSE function — expects physical parameter values."""
    return nse_from_params(x)

Owen (Grouped) Shapley Effects¶

Standard Shapley effects attribute variance to individual parameters. Owen effects (Owen 2014) extend this to a two-level hierarchy:

Level	Question	Example
Outer (Group)	Which subsystem matters most?	Soil zone vs. routing
Inner (Individual)	Within a subsystem, which parameter dominates?	LZFPM vs. UZTWM

Both levels sum to the total output variance. This is not just summing individual Shapley effects — the group effect is computed over a different game where each group is an atomic player.

In [13]:

Copied!





# Define hydrologically meaningful parameter groups
owen_groups = {
    'soil zone':     ['Uztwm', 'Uzfwm', 'Lztwm', 'Lzfpm', 'Lzfsm'],
    'impervious':    ['Adimp', 'Pctim', 'Pfree', 'Side'],
    'groundwater':   ['Lzpk', 'Lzsk', 'Uzk'],
    'deep drainage': ['Zperc', 'Rexp'],
}

# Verify all 14 parameters are covered
all_grouped = []
for g_name, g_vars in owen_groups.items():
    all_grouped.extend(g_vars)
    print(f'{g_name:>15s}: {g_vars}')
assert set(all_grouped) == set(param_names), \
    f'Missing: {set(param_names) - set(all_grouped)}'
print(f'\nAll {len(param_names)} parameters assigned to {len(owen_groups)} groups.')
# Define hydrologically meaningful parameter groups
owen_groups = {
    'soil zone':     ['Uztwm', 'Uzfwm', 'Lztwm', 'Lzfpm', 'Lzfsm'],
    'impervious':    ['Adimp', 'Pctim', 'Pfree', 'Side'],
    'groundwater':   ['Lzpk', 'Lzsk', 'Uzk'],
    'deep drainage': ['Zperc', 'Rexp'],
}

# Verify all 14 parameters are covered
all_grouped = []
for g_name, g_vars in owen_groups.items():
    all_grouped.extend(g_vars)
    print(f'{g_name:>15s}: {g_vars}')
assert set(all_grouped) == set(param_names), \
    f'Missing: {set(param_names) - set(all_grouped)}'
print(f'\nAll {len(param_names)} parameters assigned to {len(owen_groups)} groups.')

      soil zone: ['Uztwm', 'Uzfwm', 'Lztwm', 'Lzfpm', 'Lzfsm']
     impervious: ['Adimp', 'Pctim', 'Pfree', 'Side']
    groundwater: ['Lzpk', 'Lzsk', 'Uzk']
  deep drainage: ['Zperc', 'Rexp']

All 14 parameters assigned to 4 groups.

Note: With polys=[6] and k_max=5, the outer group effects are exact — all group coalitions (up to the full 14-parameter set) are enumerated. The soil zone group (5 variables) and impervious group (4 variables) are fully covered.

Owen Effects via RS-HDMR Surrogate¶

We compute Owen effects using the trained surrogate with the Vrugt et al. (2006) correlation structure. k_max = 5 covers all 4 groups (max group size = 5 for the soil zone), so both outer and inner Owen effects are exact.

In [15]:

Copied!





# Import Owen effects machinery
from shapleyx.utilities.mc_shapley import (
    collect_shapley_data,
    owen_from_data,
    owen_effects,
)

print('Computing Owen effects via RS-HDMR surrogate...')
print(f'  Groups: {len(owen_groups)}, k_max=5, N=1000')

# Collect pick-freeze data on the surrogate
data_surr = collect_shapley_data(
    analyzer.predict,
    joint_correlated,
    N=1000,
    predict_batch=analyzer.predict,
    k_max=5,
    progress=True,
)

# Compute Owen decomposition
ge_surr, ie_surr, tv_surr = owen_from_data(
    data_surr, owen_groups, param_names, k_max=5)

print(f'\nOwen effects — Surrogate:')
print(f'  Total variance: {tv_surr:.6f}')
print(f'\n  {"Group":>15s}  {"Outer Phi":>10s}')
print(f'  {"-"*15}  {"-"*10}')
for g_name, g_val in ge_surr.items():
    print(f'  {g_name:>15s}  {g_val:10.6f}')

print(f'\n  {"Parameter":>15s}  {"Inner phi":>10s}  {"Group":>15s}')
print(f'  {"-"*15}  {"-"*10}  {"-"*15}')
for g_name, g_vars in owen_groups.items():
    for v in g_vars:
        print(f'  {v:>15s}  {ie_surr[v]:10.6f}  {g_name:>15s}')

print(f'\n  Sum outer: {sum(ge_surr.values()):.6f}')
print(f'  Sum inner: {sum(ie_surr.values()):.6f}')
# Import Owen effects machinery
from shapleyx.utilities.mc_shapley import (
    collect_shapley_data,
    owen_from_data,
    owen_effects,
)

print('Computing Owen effects via RS-HDMR surrogate...')
print(f'  Groups: {len(owen_groups)}, k_max=5, N=1000')

# Collect pick-freeze data on the surrogate
data_surr = collect_shapley_data(
    analyzer.predict,
    joint_correlated,
    N=1000,
    predict_batch=analyzer.predict,
    k_max=5,
    progress=True,
)

# Compute Owen decomposition
ge_surr, ie_surr, tv_surr = owen_from_data(
    data_surr, owen_groups, param_names, k_max=5)

print(f'\nOwen effects — Surrogate:')
print(f'  Total variance: {tv_surr:.6f}')
print(f'\n  {"Group":>15s}  {"Outer Phi":>10s}')
print(f'  {"-"*15}  {"-"*10}')
for g_name, g_val in ge_surr.items():
    print(f'  {g_name:>15s}  {g_val:10.6f}')

print(f'\n  {"Parameter":>15s}  {"Inner phi":>10s}  {"Group":>15s}')
print(f'  {"-"*15}  {"-"*10}  {"-"*15}')
for g_name, g_vars in owen_groups.items():
    for v in g_vars:
        print(f'  {v:>15s}  {ie_surr[v]:10.6f}  {g_name:>15s}')

print(f'\n  Sum outer: {sum(ge_surr.values()):.6f}')
print(f'  Sum inner: {sum(ie_surr.values()):.6f}')

Computing Owen effects via RS-HDMR surrogate...
  Groups: 4, k_max=5, N=1000

MC Shapley: 100%|██████████| 6945000/6945000 [00:31<00:00, 218849.23evals/s]

Owen effects — Surrogate:
  Total variance: 0.004368

            Group   Outer Phi
  ---------------  ----------
        soil zone    0.001749
       impervious    0.000818
      groundwater    0.001046
    deep drainage    0.000755

        Parameter   Inner phi            Group
  ---------------  ----------  ---------------
            Uztwm    0.000000        soil zone
            Uzfwm    0.000089        soil zone
            Lztwm    0.000153        soil zone
            Lzfpm    0.000692        soil zone
            Lzfsm    0.000714        soil zone
            Adimp    0.000194       impervious
            Pctim    0.000098       impervious
            Pfree    0.000290       impervious
             Side    0.000263       impervious
             Lzpk    0.000587      groundwater
             Lzsk    0.000381      groundwater
              Uzk    0.000063      groundwater
            Zperc    0.000137    deep drainage
             Rexp    0.000707    deep drainage

  Sum outer: 0.004368
  Sum inner: 0.004368

Owen Effects — Direct SAC-SMA Validation¶

To validate the surrogate, we compute Owen effects directly on the SAC-SMA model at a smaller sample size (N=200, ).

In [16]:

Copied!





print('Computing Owen effects directly on SAC-SMA model...')
print(f'  N=200, k_max=5 — full group coverage')

data_direct = collect_shapley_data(
    sacsma_nse_wrapper,
    joint_correlated,
    N=200,
    k_max=5,
    progress=True,
)

ge_direct, ie_direct, tv_direct = owen_from_data(
    data_direct, owen_groups, param_names, k_max=5)

print(f'\nOwen effects — Direct SAC-SMA:')
print(f'  Total variance: {tv_direct:.6f}')
for g_name, g_val in ge_direct.items():
    print(f'  {g_name:>15s}: {g_val:.6f}')
print('Computing Owen effects directly on SAC-SMA model...')
print(f'  N=200, k_max=5 — full group coverage')

data_direct = collect_shapley_data(
    sacsma_nse_wrapper,
    joint_correlated,
    N=200,
    k_max=5,
    progress=True,
)

ge_direct, ie_direct, tv_direct = owen_from_data(
    data_direct, owen_groups, param_names, k_max=5)

print(f'\nOwen effects — Direct SAC-SMA:')
print(f'  Total variance: {tv_direct:.6f}')
for g_name, g_val in ge_direct.items():
    print(f'  {g_name:>15s}: {g_val:.6f}')

Computing Owen effects directly on SAC-SMA model...
  N=200, k_max=5 — full group coverage

MC Shapley: 100%|██████████| 1389000/1389000 [56:29<00:00, 409.81evals/s]

Owen effects — Direct SAC-SMA:
  Total variance: 0.005164
        soil zone: 0.002349
       impervious: 0.000287
      groundwater: 0.001269
    deep drainage: 0.001259

Comparison: Surrogate vs Direct Owen Effects¶

In [17]:

Copied!





group_names = list(owen_groups.keys())
surr_vals = [ge_surr[g] for g in group_names]
direct_vals = [ge_direct[g] for g in group_names]

df_groups = pd.DataFrame({
    'Group': group_names,
    'Surrogate': surr_vals,
    'Direct': direct_vals,
})
df_groups['Delta'] = np.abs(df_groups['Surrogate'] - df_groups['Direct'])

print('Group-level Owen effects:')
display(df_groups.round(6))

# Individual comparison
indiv_data = []
for v in param_names:
    g = [gn for gn, gv in owen_groups.items() if v in gv][0]
    indiv_data.append({
        'Parameter': v, 'Group': g,
        'Surrogate': ie_surr[v], 'Direct': ie_direct[v],
    })
df_indiv = pd.DataFrame(indiv_data)
df_indiv['Delta'] = np.abs(df_indiv['Surrogate'] - df_indiv['Direct'])
print('\nIndividual Owen effects:')
display(df_indiv.round(6))

# Rank correlation
from scipy.stats import spearmanr
rho_g, p_g = spearmanr(surr_vals, direct_vals)
surr_inner = np.array([ie_surr[v] for v in param_names])
dir_inner = np.array([ie_direct[v] for v in param_names])
rho_i, p_i = spearmanr(surr_inner, dir_inner)
print(f'\nRank correlation (groups): rho={rho_g:.4f} (p={p_g:.4f})')
print(f'Rank correlation (indiv):  rho={rho_i:.4f} (p={p_i:.4f})')
group_names = list(owen_groups.keys())
surr_vals = [ge_surr[g] for g in group_names]
direct_vals = [ge_direct[g] for g in group_names]

df_groups = pd.DataFrame({
    'Group': group_names,
    'Surrogate': surr_vals,
    'Direct': direct_vals,
})
df_groups['Delta'] = np.abs(df_groups['Surrogate'] - df_groups['Direct'])

print('Group-level Owen effects:')
display(df_groups.round(6))

# Individual comparison
indiv_data = []
for v in param_names:
    g = [gn for gn, gv in owen_groups.items() if v in gv][0]
    indiv_data.append({
        'Parameter': v, 'Group': g,
        'Surrogate': ie_surr[v], 'Direct': ie_direct[v],
    })
df_indiv = pd.DataFrame(indiv_data)
df_indiv['Delta'] = np.abs(df_indiv['Surrogate'] - df_indiv['Direct'])
print('\nIndividual Owen effects:')
display(df_indiv.round(6))

# Rank correlation
from scipy.stats import spearmanr
rho_g, p_g = spearmanr(surr_vals, direct_vals)
surr_inner = np.array([ie_surr[v] for v in param_names])
dir_inner = np.array([ie_direct[v] for v in param_names])
rho_i, p_i = spearmanr(surr_inner, dir_inner)
print(f'\nRank correlation (groups): rho={rho_g:.4f} (p={p_g:.4f})')
print(f'Rank correlation (indiv):  rho={rho_i:.4f} (p={p_i:.4f})')

Group-level Owen effects:

	Group	Surrogate	Direct	Delta
0	soil zone	1.749e-03	2.349e-03	6.000e-04
1	impervious	8.180e-04	2.870e-04	5.310e-04
2	groundwater	1.046e-03	1.269e-03	2.230e-04
3	deep drainage	7.550e-04	1.259e-03	5.040e-04

Individual Owen effects:

	Parameter	Group	Surrogate	Direct	Delta
0	Uztwm	soil zone	0.000e+00	1.960e-04	1.960e-04
1	Uzfwm	soil zone	8.900e-05	6.800e-05	2.100e-05
2	Lztwm	soil zone	1.530e-04	3.310e-04	1.780e-04
3	Lzfpm	soil zone	6.920e-04	7.540e-04	6.300e-05
4	Lzfsm	soil zone	7.140e-04	7.900e-04	7.600e-05
5	Adimp	impervious	1.940e-04	-3.000e-06	1.970e-04
6	Uzk	groundwater	6.300e-05	2.360e-04	1.720e-04
7	Lzpk	groundwater	5.870e-04	4.650e-04	1.220e-04
8	Lzsk	groundwater	3.810e-04	5.820e-04	2.010e-04
9	Zperc	deep drainage	1.370e-04	6.610e-04	5.240e-04
10	Rexp	deep drainage	7.070e-04	6.450e-04	6.200e-05
11	Pctim	impervious	9.800e-05	-1.410e-04	2.390e-04
12	Pfree	impervious	2.900e-04	2.920e-04	2.000e-06
13	Side	impervious	2.630e-04	2.890e-04	2.600e-05

Rank correlation (groups): rho=0.8000 (p=0.2000)
Rank correlation (indiv):  rho=0.7319 (p=0.0029)

In [18]:

Copied!





fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))

x = np.arange(len(group_names))
w = 0.35
ax1.bar(x - w/2, surr_vals, w, label='Surrogate', color='steelblue')
ax1.bar(x + w/2, direct_vals, w, label='Direct SAC-SMA', color='coral')
ax1.set_xticks(x)
ax1.set_xticklabels(group_names, rotation=30, ha='right')
ax1.set_ylabel('Owen Group Effect')
ax1.set_title('Group-Level Owen Effects')
ax1.legend()

# Individual effects
group_colours = {
    'soil zone': 'steelblue',
    'impervious': 'coral',
    'groundwater': 'seagreen',
    'deep drainage': 'goldenrod',
}
bar_colours = [group_colours[
    [g for g, gv in owen_groups.items() if v in gv][0]
] for v in param_names]

x2 = np.arange(len(param_names))
ax2.bar(x2 - w/2, surr_inner, w, label='Surrogate',
        color=bar_colours, alpha=0.7)
ax2.bar(x2 + w/2, dir_inner, w, label='Direct',
        color=bar_colours, alpha=0.3)
ax2.set_xticks(x2)
ax2.set_xticklabels(param_names, rotation=45, ha='right')
ax2.set_ylabel('Owen Individual Effect')
ax2.set_title('Individual (Within-Group) Owen Effects')

plt.tight_layout()
plt.show()
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))

x = np.arange(len(group_names))
w = 0.35
ax1.bar(x - w/2, surr_vals, w, label='Surrogate', color='steelblue')
ax1.bar(x + w/2, direct_vals, w, label='Direct SAC-SMA', color='coral')
ax1.set_xticks(x)
ax1.set_xticklabels(group_names, rotation=30, ha='right')
ax1.set_ylabel('Owen Group Effect')
ax1.set_title('Group-Level Owen Effects')
ax1.legend()

# Individual effects
group_colours = {
    'soil zone': 'steelblue',
    'impervious': 'coral',
    'groundwater': 'seagreen',
    'deep drainage': 'goldenrod',
}
bar_colours = [group_colours[
    [g for g, gv in owen_groups.items() if v in gv][0]
] for v in param_names]

x2 = np.arange(len(param_names))
ax2.bar(x2 - w/2, surr_inner, w, label='Surrogate',
        color=bar_colours, alpha=0.7)
ax2.bar(x2 + w/2, dir_inner, w, label='Direct',
        color=bar_colours, alpha=0.3)
ax2.set_xticks(x2)
ax2.set_xticklabels(param_names, rotation=45, ha='right')
ax2.set_ylabel('Owen Individual Effect')
ax2.set_title('Individual (Within-Group) Owen Effects')

plt.tight_layout()
plt.show()

Key Takeaways¶

Owen effects provide a natural two-level decomposition — group-level effects tell you which subsystem matters, and individual effects show which parameters matter within each subsystem.
The RS-HDMR surrogate accurately recovers Owen effects — both group and individual rankings match the direct computation at a fraction of the cost.
No additional model evaluations — owen_from_data() reuses the same pick-freeze data as standard Shapley. The grouping is purely a post-processing step over the same v(S) values.
Group Shapley ≠ sum of individual Shapley — the Owen decomposition accounts for within-group correlation and between-group interactions correctly. Simply summing individual Shapley effects for parameters in a group gives a different answer.