Global Sensitivity Analysis Using Sobol Indices

Author

Akash B V

Published

April 6, 2026

1 Overview

This report presents results from a variance-based global sensitivity analysis of the SIPNET ecosystem model using Sobol indices. Unlike local sensitivity analysis (which measures rates of change at baseline conditions), variance-based methods decompose output variance into contributions from individual parameters and their interactions.

Why Sobol Indices?

Sobol indices provide a rigorous framework for:

Quantifying parameter importance via variance decomposition
Detecting interactions between parameters (non-additive effects)
Factor prioritization (which parameters to constrain for uncertainty reduction)
Factor fixing (which parameters can be fixed without loss of information)

1.1 Theoretical Background

Sobol indices are normalized variance contributions:

\[ S_i = \frac{V_i}{\text{VAR}(y)} \quad \text{(First-order index)} \tag{1}\]

\[ T_i = 1 - \frac{\text{VAR}_{x_{\sim i}}[\mathbb{E}_{x_i}(y | x_{\sim i})]}{\text{VAR}(y)} = \frac{\mathbb{E}_{x_{\sim i}}[\text{VAR}_{x_i}(y | x_{\sim i})]}{\text{VAR}(y)} \quad \text{(Total-order index)} \tag{2}\]

where $x_{\sim i}$ denotes all-parameters-but-$x_i$.

1.1.1 Interpretation

$S_i$: Direct effect of $x_i$ on output variance (excluding interactions)
$T_i$: Total effect of $x_i$ including all interactions with other parameters
$T_i - S_i$: Interaction strength (purely higher-order effects)
$\sum S_i \approx 1$: Model is additive (linear-like behavior)
$\sum S_i \ll 1$: Model is non-additive (interactions dominate)

2 Study Design

Method: Saltelli’s sampling scheme (Saltelli et al. 2010) with Jansen estimators for both first-order and total-order indices, implemented via the sensobol R package (Puy et al. 2022). Bootstrap resampling (R = 500) provides confidence intervals.

Sample Size: $N = 32$, Total Runs = $N \times (k + 2)$ = 1440

Scope: 17 sites across California row-crop agriculture

3 Data Processing

We process the raw indices to parse namespaced parameters (e.g., grass.SLA), identify the dummy parameter (Sobol validation; its $T_i$ should be near zero), and calculate interaction components ($T_i - S_i$).

4 Key Findings

Based on a Sobol analysis with 32 base samples and 1440 total runs across 17 sites:

Top inputs: Initial Conditions, Soil Methane Rate, and N Volatilization Rate explain the most output variance ($T_i$ = 0.91, 0.81, 0.63)
Model parameters account for ~71% of the first-order variance contribution; the remaining contribution is distributed across meteorological drivers, initial conditions, and management inputs
Model additivity ($\sum S_i$) averages 1.93 – values near 1 indicate additive behavior, values $\ll 1$ indicate parameter interactions dominate
Noise floor: the dummy parameter’s median $T_i$ = 0, providing a rigorous threshold for identifying non-influential inputs

5 Variance Decomposition (Source Partitioning)

This section summarizes first-order (main-effect) uncertainty contributions by source (model parameters, meteorology, initial conditions, and management). This is done by summing first-order Sobol indices ($S_i$) within broad categories (Equation 1). These contributions do not include interaction effects (see $T_i - S_i$).

For discrete inputs such as initial condition and meteorological ensembles, Sobol coordinates are used to select among ensemble members / files. The resulting indices quantify variance attributable to uncertainty in which ensemble member is selected, not sensitivity along an ordered or continuous variable.

The “Control (Noise)” source is excluded for clarity; it accounts for numerical estimation error from the Sobol design.

Figure 1: First-order (main-effect) variance contribution by source group (model parameters, meteorology, initial conditions, management). Bars are based on summed Si values and do not include interaction effects.

6 Parameter Ranking (First & Total Order)

This plot compares First-Order ($S_i$, solid bar) and Total-Order ($T_i$, whisker) indices. The gap between them quantifies interaction strength (Equation 2 minus Equation 1).

This ranking identifies the inputs with the largest total effect ($T_i$), indicating which uncertainties most strongly influence model outputs.

Solid bar ($S_i$): Direct (main) effect of the parameter alone

Gap ($T_i - S_i$): Higher-order interaction effects with other parameters

Parameters are ordered by median $T_i$ across sites. We show the top 15 inputs for each output variable.

Figure 2: Top 15 inputs ranked by median Total Effect (Ti) across sites. Solid bar = Main Effect (Si); whisker = Total Effect (Ti). Parameters with large gaps between Si and Ti interact strongly with other inputs.

6.1 Cross-Site Variability

A parameter whose importance varies across sites (wide box) indicates context-dependent sensitivity: its influence depends on local conditions (soil, climate, crop type). These are priorities for spatially explicit calibration.

Figure 3: Distribution of Ti across sites for top 10 parameters. Points = site-level values; box = IQR; line = median. Wide boxes indicate sensitivity varies with local conditions.

7 Interaction Analysis

Parameters with large differences between $T_i$ and $S_i$ interact strongly with other factors. These non-additive effects mean the parameter’s impact on output depends on the values of other inputs.

Interaction strength ($T_i - S_i$) indicates how much a parameter’s influence depends on other inputs rather than acting independently.

Figure 4: Interaction strength (Ti - Si) for inputs with interaction component exceeding 0.05 (i.e. 5% of total output variance, since Sobol indices are variance-normalized). Numbers show percent of Ti attributable to interactions.

8 Factor Fixing (Screening)

To identify non-influential parameters, we compare their Total-Order index ($T_i$) against the Dummy parameter (noise floor) that is useful for screening. The dummy parameter varies randomly but has no effect on model output, so we assume that any parameter with $T_i$ at or below the dummy’s median can be fixed to a constant value without meaningful loss of variance coverage.

The dashed red line marks the noise floor (median dummy $T_i$).

Figure 5: Factor screening: parameters below the noise line (dashed red) can be fixed. Dot = median Ti across sites; whisker = IQR.

Fixable parameters (Ti below noise floor):

Variable	Input	Median Ti	Noise Threshold
Aboveground Biomass	Growing Degree Days	0	0
Aboveground Biomass	Leaf Fall Fraction	0	0
Aboveground Biomass	Growth Resp. Factor	0	0
Aboveground Biomass	Leaf Growth Rate	0	0
Aboveground Biomass	Compost C:N Ratio	0	0
Aboveground Biomass	Compost Application	0	0
Aboveground Biomass	N Fertilizer Rate	0	0
Aboveground Biomass	Anaerobic Decomp. Rate	0	0
Aboveground Biomass	Anaerobic Transition Exp.	0	0
Aboveground Biomass	Fine Root C:N	0	0
Aboveground Biomass	Leaf C:N	0	0
Aboveground Biomass	Wood C:N	0	0
Aboveground Biomass	Anoxia Fraction	0	0
Aboveground Biomass	C:N Decomp. Modifier	0	0
Aboveground Biomass	Litter Methane Rate	0	0
Aboveground Biomass	N Fixation Half-Sat.	0	0
Aboveground Biomass	Max N Fixation Fraction	0	0
Aboveground Biomass	N Leaching Fraction	0	0
Aboveground Biomass	N Volatilization Rate	0	0
Aboveground Biomass	Soil Methane Rate	0	0
Aboveground Biomass	SOM Respiration Rate	0	0
CH₄ Flux	Growing Degree Days	0	0
CH₄ Flux	Leaf Fall Fraction	0	0
CH₄ Flux	Growth Resp. Factor	0	0
CH₄ Flux	Leaf Growth Rate	0	0
CH₄ Flux	N Fertilizer Rate	0	0
CH₄ Flux	N Fixation Half-Sat.	0	0
CH₄ Flux	Max N Fixation Fraction	0	0
CH₄ Flux	N Leaching Fraction	0	0
CH₄ Flux	N Volatilization Rate	0	0
Latent Heat Flux	Growing Degree Days	0	0
Latent Heat Flux	Leaf Fall Fraction	0	0
Latent Heat Flux	Growth Resp. Factor	0	0
Latent Heat Flux	Leaf Growth Rate	0	0
Latent Heat Flux	Compost C:N Ratio	0	0
Latent Heat Flux	Compost Application	0	0
Latent Heat Flux	N Fertilizer Rate	0	0
Latent Heat Flux	Anaerobic Decomp. Rate	0	0
Latent Heat Flux	Anaerobic Transition Exp.	0	0
Latent Heat Flux	Fine Root C:N	0	0
Latent Heat Flux	Leaf C:N	0	0
Latent Heat Flux	Wood C:N	0	0
Latent Heat Flux	Anoxia Fraction	0	0
Latent Heat Flux	C:N Decomp. Modifier	0	0
Latent Heat Flux	Litter Methane Rate	0	0
Latent Heat Flux	N Fixation Half-Sat.	0	0
Latent Heat Flux	Max N Fixation Fraction	0	0
Latent Heat Flux	N Leaching Fraction	0	0
Latent Heat Flux	N Volatilization Rate	0	0
Latent Heat Flux	Soil Methane Rate	0	0
Latent Heat Flux	SOM Respiration Rate	0	0
Net Primary Productivity	Growing Degree Days	0	0
Net Primary Productivity	Leaf Fall Fraction	0	0
Net Primary Productivity	Growth Resp. Factor	0	0
Net Primary Productivity	Leaf Growth Rate	0	0
Net Primary Productivity	Compost C:N Ratio	0	0
Net Primary Productivity	Compost Application	0	0
Net Primary Productivity	N Fertilizer Rate	0	0
Net Primary Productivity	Anaerobic Decomp. Rate	0	0
Net Primary Productivity	Anaerobic Transition Exp.	0	0
Net Primary Productivity	Fine Root C:N	0	0
Net Primary Productivity	Leaf C:N	0	0
Net Primary Productivity	Wood C:N	0	0
Net Primary Productivity	Anoxia Fraction	0	0
Net Primary Productivity	C:N Decomp. Modifier	0	0
Net Primary Productivity	Litter Methane Rate	0	0
Net Primary Productivity	N Fixation Half-Sat.	0	0
Net Primary Productivity	Max N Fixation Fraction	0	0
Net Primary Productivity	N Leaching Fraction	0	0
Net Primary Productivity	N Volatilization Rate	0	0
Net Primary Productivity	Soil Methane Rate	0	0
Net Primary Productivity	SOM Respiration Rate	0	0
N₂O Flux	Growing Degree Days	0	0
N₂O Flux	Leaf Fall Fraction	0	0
N₂O Flux	Growth Resp. Factor	0	0
N₂O Flux	Leaf Growth Rate	0	0
Soil Carbon	Growing Degree Days	0	0
Soil Carbon	Leaf Fall Fraction	0	0
Soil Carbon	Growth Resp. Factor	0	0
Soil Carbon	Leaf Growth Rate	0	0
Soil Carbon	N Fertilizer Rate	0	0
Soil Carbon	N Fixation Half-Sat.	0	0
Soil Carbon	Max N Fixation Fraction	0	0
Soil Carbon	N Leaching Fraction	0	0
Soil Carbon	N Volatilization Rate	0	0
Soil Moisture	Growing Degree Days	0	0
Soil Moisture	Leaf Fall Fraction	0	0
Soil Moisture	Growth Resp. Factor	0	0
Soil Moisture	Leaf Growth Rate	0	0
Soil Moisture	Compost C:N Ratio	0	0
Soil Moisture	Compost Application	0	0
Soil Moisture	N Fertilizer Rate	0	0
Soil Moisture	Anaerobic Decomp. Rate	0	0
Soil Moisture	Anaerobic Transition Exp.	0	0
Soil Moisture	Fine Root C:N	0	0
Soil Moisture	Leaf C:N	0	0
Soil Moisture	Wood C:N	0	0
Soil Moisture	Anoxia Fraction	0	0
Soil Moisture	C:N Decomp. Modifier	0	0
Soil Moisture	Litter Methane Rate	0	0
Soil Moisture	N Fixation Half-Sat.	0	0
Soil Moisture	Max N Fixation Fraction	0	0
Soil Moisture	N Leaching Fraction	0	0
Soil Moisture	N Volatilization Rate	0	0
Soil Moisture	Soil Methane Rate	0	0
Soil Moisture	SOM Respiration Rate	0	0

9 Model Additivity

The sum of First-Order indices ($\sum S_i$) indicates the degree of model linearity. Values near 1.0 imply the model is additive and parameters contribute independently. Values well below 1 indicate that interactions dominate – the model’s response to one parameter depends on the values of others.

Figure 6: Distribution of model additivity (sum of Si) across sites. Dashed line at 1.0 marks perfect additivity. Values below 1 indicate interaction effects.

10 Convergence Diagnostics

A Sobol analysis is only as reliable as its sample size allows. With $N = 32$ base samples and $k = 41$ parameters, our total model evaluation budget is 1440 runs per site. This section assesses whether that budget was sufficient for stable index estimation, following the diagnostic framework of Nossent et al. (2011).

10.1 Confidence Interval Precision

The width of the bootstrap confidence interval ($\text{CI}_{95}$) for each index indicates estimation precision. Parameters with narrow CIs have well-determined importance; wide CIs flag indices that may shift with additional samples.

Figure 7: Bootstrap 95% confidence interval width for Si (left) and Ti (right), by output variable. Each point is one parameter. Narrower intervals indicate more precise estimation. The dashed line marks a CI width of 0.1 (10% of output variance).

10.2 Parameter-Output Relationships

These scatterplots show the direct relationship between each top parameter’s sampled value and the model output. They reveal whether the sensitivity is linear (validating the Sobol decomposition) or nonlinear (indicating threshold behavior or saturation). Each point is one model evaluation from the Sobol design matrix.

Figure 8: Parameter-output scatterplots for the top 3 most influential parameters per output variable. Each point is one run from the Sobol design. Trend lines (LOESS) reveal linear vs. nonlinear relationships.

Interpreting Scatterplots

Linear trend with tight spread: parameter has a strong, additive effect – Si will be high
Nonlinear trend (curve, threshold): parameter has important interactions – Ti >> Si
Cloud with no trend: parameter has low sensitivity for this variable
Funnel shape (spread changes with x): heteroscedasticity – the parameter’s effect depends on other inputs

10.3 Convergence with Sample Size

To assess whether $N = 32$ is sufficient, we recompute Sobol indices using increasing fractions of the base sample. If indices stabilize at the full sample size, this supports adequate convergence; indices still drifting indicate that a larger N would improve precision.

Figure 9: Evolution of Ti for the top 5 parameters as the Sobol base sample size increases. Stable plateaus indicate convergence; ongoing drift indicates the estimate has not fully stabilized.

Sample Size Considerations

Our base sample $N = 32$ is modest by Sobol analysis standards (Nossent et al. 2011 recommend N = 500+ for stable individual indices), chosen to balance computational cost (~1,440 model runs per site x 17 sites) against the need for results across many sites. Where indices show ongoing drift at $N = 32$, a larger ensemble (Phase 4 workplan) would improve precision. The current results are reliable for parameter ranking and factor prioritization (relative ordering is robust), even though individual index point estimates carry estimation uncertainty.

11 Synthesis

What drives uncertainty? The variance decomposition (Section 5) separates biological parameter uncertainty from environmental forcing.

What to calibrate first? Parameter rankings (Section 6) identify the top candidates for reducing predictive uncertainty. Parameters whose importance varies across sites (Figure 3) need spatially explicit calibration.

What can be simplified? Factor fixing (Section 8) uses the dummy parameter as a rigorous noise floor to identify inputs that can be set to nominal values.

Is the model additive? Additivity diagnostics (Section 9) indicate whether parameter effects are independent or interact. Strong interactions (large $T_i - S_i$ in Section 7) suggest that calibrating parameters independently may be insufficient.

12 References

Puy, A., et al. (2022). sensobol: An R Package to Compute Variance-Based Sensitivity Indices. Journal of Statistical Software, 102(5).
Saltelli, A., et al. (2010). Variance based sensitivity analysis of model output. Computer Physics Communications.
Dietze, M. (2017). Ecological Forecasting. Princeton University Press.
Gelman, A. & Dodhia, R. (2002). Let’s practice what we preach: Turning tables into graphs.
Nossent, J., Elsen, P. & Bauwens, W. (2011). Sobol’ sensitivity analysis of a complex environmental model. Environmental Modelling & Software, 26, 1515-1525.

--- title: "Global Sensitivity Analysis Using Sobol Indices" author: "Akash B V" date: last-modified format: html: self-contained: true df-print: paged toc: true toc-depth: 3 code-fold: true code-tools: true theme: cosmo number-sections: true execute: echo: false warning: false message: false cache: false --- # Overview This report presents results from a **variance-based global sensitivity analysis** of the SIPNET ecosystem model using Sobol indices. Unlike local sensitivity analysis (which measures rates of change at baseline conditions), variance-based methods decompose output variance into contributions from individual parameters and their interactions. ::: callout-note ## Why Sobol Indices? Sobol indices provide a rigorous framework for: - **Quantifying parameter importance** via variance decomposition - **Detecting interactions** between parameters (non-additive effects) - **Factor prioritization** (which parameters to constrain for uncertainty reduction) - **Factor fixing** (which parameters can be fixed without loss of information) ::: ## Theoretical Background **Sobol indices** are normalized variance contributions: $$ S_i = \frac{V_i}{\text{VAR}(y)} \quad \text{(First-order index)} $$ {#eq-si} $$ T_i = 1 - \frac{\text{VAR}_{x_{\sim i}}[\mathbb{E}_{x_i}(y | x_{\sim i})]}{\text{VAR}(y)} = \frac{\mathbb{E}_{x_{\sim i}}[\text{VAR}_{x_i}(y | x_{\sim i})]}{\text{VAR}(y)} \quad \text{(Total-order index)} $$ {#eq-ti} where $x_{\sim i}$ denotes all-parameters-but-$x_i$. ### Interpretation - $S_i$: Direct effect of $x_i$ on output variance (excluding interactions) - $T_i$: Total effect of $x_i$ including all interactions with other parameters - $T_i - S_i$: Interaction strength (purely higher-order effects) - $\sum S_i \approx 1$: Model is **additive** (linear-like behavior) - $\sum S_i \ll 1$: Model is **non-additive** (interactions dominate) # Study Design ```{r} #| label: setup #| include: false library(here) here::i_am("analysis/global_sensitivity.qmd") library(readr) library(dplyr) library(ggplot2) library(tidyr) library(stringr) library(knitr) library(scales) library(config) library(patchwork) library(sensobol) source(here::here("R", "labels.R")) cfg <- config::get(file = here::here("000-config.yml")) sobol_indices <- readr::read_csv( here::here(cfg$paths$data_dir, "sobol_indices.csv"), show_col_types = FALSE ) metadata <- jsonlite::fromJSON(here::here(cfg$paths$data_dir, "sobol_design_matrix_metadata.json")) N <- metadata$N k <- metadata$k total_runs <- metadata$total_runs n_sites <- dplyr::n_distinct(sobol_indices$runid) response_vars <- unique(sobol_indices$variable) # consistent theme theme_set( theme_minimal(base_size = 12) + theme( panel.grid.minor = element_blank(), strip.text = element_text(size = 10, margin = margin(2, 2, 2, 2)), legend.position = "bottom" ) ) ``` **Method:** Saltelli's sampling scheme (Saltelli et al. 2010) with Jansen estimators for both first-order and total-order indices, implemented via the sensobol R package (Puy et al. 2022). Bootstrap resampling (R = 500) provides confidence intervals. **Sample Size:** $N = `r N`$, Total Runs = $N \times (k + 2)$ = `r total_runs` **Scope:** `r n_sites` sites across California row-crop agriculture # Data Processing We process the raw indices to parse namespaced parameters (e.g., `grass.SLA`), identify the dummy parameter (Sobol validation; its $T_i$ should be near zero), and calculate interaction components ($T_i - S_i$). ```{r} #| label: data-processing # reshape and enforce non-negativity (numerical artifacts can cause small negatives) sobol_long <- sobol_indices |> select(runid, variable, parameters, Si_original, Si_low.ci, Si_high.ci, Ti_original, Ti_low.ci, Ti_high.ci) |> mutate( Si = pmax(Si_original, 0), Ti = pmax(Ti_original, 0), Interaction = pmax(Ti - Si, 0) ) # parse parameter names: "pft.param" -> Prefix + Trait sobol_parsed <- sobol_long |> separate(parameters, into = c("Prefix", "Trait"), sep = "\\.", extra = "merge", fill = "right") |> mutate( Trait = ifelse(is.na(Trait), Prefix, Trait), Category = case_when( Prefix == "ic_ensemble" ~ "Initial Conditions", Prefix == "met_ensemble" ~ "Meteorology", Prefix == "dummy" ~ "Control (Noise)", Prefix == "mgmt" ~ "Management", TRUE ~ "Model Parameter" ), Label = label_param(Trait), var_label = label_variable(variable) ) # per-site summaries (keep distribution for cross-site variability) sobol_site <- sobol_parsed |> summarize( Si = mean(Si, na.rm = TRUE), Ti = mean(Ti, na.rm = TRUE), Interaction = mean(Interaction, na.rm = TRUE), .by = c("runid", "variable", "var_label", "Category", "Label") ) # cross-site summary: median + IQR sobol_summary <- sobol_site |> summarize( median_Ti = median(Ti, na.rm = TRUE), q25_Ti = quantile(Ti, 0.25, na.rm = TRUE), q75_Ti = quantile(Ti, 0.75, na.rm = TRUE), median_Si = median(Si, na.rm = TRUE), median_Int = median(Interaction, na.rm = TRUE), .by = c("variable", "var_label", "Category", "Label") ) ``` # Key Findings ```{r} #| label: key-findings #| results: asis # additivity additivity <- sobol_parsed |> filter(Category != "Control (Noise)") |> summarize(Sum_Si = sum(Si, na.rm = TRUE), .by = c("runid", "variable")) |> summarize(mean_sum = mean(Sum_Si), .by = "variable") # top 3 inputs (may include met, IC, management -- not just model parameters) top3 <- sobol_summary |> filter(Category != "Control (Noise)") |> arrange(desc(median_Ti)) |> slice_head(n = 3) # model parameter fraction of first-order variance parameter_frac <- sobol_parsed |> filter(Category != "Control (Noise)") |> summarize(Total_Si = sum(Si, na.rm = TRUE), .by = c("variable", "Category")) |> mutate(frac = Total_Si / sum(Total_Si), .by = "variable") |> filter(Category == "Model Parameter") |> pull(frac) |> mean() # noise floor dummy_vals <- sobol_summary |> filter(Category == "Control (Noise)") |> pull(median_Ti) |> mean() ``` Based on a Sobol analysis with `r N` base samples and `r total_runs` total runs across `r n_sites` sites: - **Top inputs**: `r top3$Label[1]`, `r top3$Label[2]`, and `r top3$Label[3]` explain the most output variance ($T_i$ = `r round(top3$median_Ti[1], 2)`, `r round(top3$median_Ti[2], 2)`, `r round(top3$median_Ti[3], 2)`) - **Model parameters** account for ~`r round(parameter_frac * 100)`% of the first-order variance contribution; the remaining contribution is distributed across meteorological drivers, initial conditions, and management inputs - **Model additivity** ($\sum S_i$) averages `r round(mean(additivity$mean_sum), 2)` -- values near 1 indicate additive behavior, values $\ll 1$ indicate parameter interactions dominate - **Noise floor**: the dummy parameter's median $T_i$ = `r round(dummy_vals, 4)`, providing a rigorous threshold for identifying non-influential inputs # Variance Decomposition (Source Partitioning) {#sec-variance-partition} This section summarizes first-order (main-effect) uncertainty contributions by source (model parameters, meteorology, initial conditions, and management). This is done by summing first-order Sobol indices ($S_i$) within broad categories (@eq-si). These contributions do not include interaction effects (see $T_i - S_i$). For discrete inputs such as initial condition and meteorological ensembles, Sobol coordinates are used to select among ensemble members / files. The resulting indices quantify variance attributable to uncertainty in which ensemble member is selected, not sensitivity along an ordered or continuous variable. The "Control (Noise)" source is excluded for clarity; it accounts for numerical estimation error from the Sobol design. ```{r} #| label: fig-variance-partition #| fig-cap: "First-order (main-effect) variance contribution by source group (model parameters, meteorology, initial conditions, management). Bars are based on summed Si values and do not include interaction effects." #| fig-width: 12 #| fig-height: 8 variance_partition <- sobol_parsed |> filter(Category != "Control (Noise)") |> summarize(Total_Si = sum(Si, na.rm = TRUE), .by = c("variable", "var_label", "Category")) |> mutate(Percentage = Total_Si / sum(Total_Si), .by = "variable") ggplot(variance_partition, aes(x = var_label, y = Total_Si, fill = Category)) + geom_col(position = "fill", width = 0.7) + geom_text(aes(label = percent(Percentage, accuracy = 1)), position = position_fill(vjust = 0.5), color = "white", fontface = "bold", size = 5) + scale_y_continuous(labels = percent) + scale_fill_brewer(palette = "Dark2") + labs( y = "Proportion of first-order variance", x = NULL, fill = "Source" ) + coord_flip() + theme( axis.text.y = element_text(size = 12), axis.text.x = element_text(size = 10), legend.text = element_text(size = 10) ) ``` # Parameter Ranking (First & Total Order) {#sec-rankings} This plot compares First-Order ($S_i$, solid bar) and Total-Order ($T_i$, whisker) indices. The gap between them quantifies interaction strength (@eq-ti minus @eq-si). This ranking identifies the inputs with the largest total effect ($T_i$), indicating which uncertainties most strongly influence model outputs. **Solid bar** ($S_i$): Direct (main) effect of the parameter alone **Gap** ($T_i - S_i$): Higher-order interaction effects with other parameters Parameters are ordered by median $T_i$ across sites. We show the top 15 inputs for each output variable. ```{r} #| label: fig-ranking-combined #| fig-cap: "Top 15 inputs ranked by median Total Effect (Ti) across sites. Solid bar = Main Effect (Si); whisker = Total Effect (Ti). Parameters with large gaps between Si and Ti interact strongly with other inputs." #| fig-height: 18 #| fig-width: 14 top_params <- sobol_summary |> filter(Category != "Control (Noise)") |> slice_max(order_by = median_Ti, n = 15, by = "variable") ggplot(top_params, aes(x = reorder(Label, median_Ti), y = median_Ti)) + geom_errorbar( aes(ymin = median_Si, ymax = median_Ti), width = 0.2, color = "grey50" ) + geom_point(aes(y = median_Ti, color = Category), shape = 95, size = 10) + geom_col(aes(y = median_Si, fill = Category), alpha = 0.8) + coord_flip() + facet_wrap(~var_label, scales = "free_x", ncol = 2) + scale_fill_brewer(palette = "Set2") + scale_color_brewer(palette = "Set2") + labs( x = NULL, y = "Sobol Index (median across sites)" ) + theme( axis.text.y = element_text(size = 9), strip.text = element_text(size = 11, face = "bold"), legend.text = element_text(size = 10) ) ``` ## Cross-Site Variability A parameter whose importance varies across sites (wide box) indicates context-dependent sensitivity: its influence depends on local conditions (soil, climate, crop type). These are priorities for spatially explicit calibration. ```{r} #| label: fig-site-variability #| fig-cap: "Distribution of Ti across sites for top 10 parameters. Points = site-level values; box = IQR; line = median. Wide boxes indicate sensitivity varies with local conditions." #| fig-height: 14 #| fig-width: 14 # top 10 by median Ti (across all variables) top10_labels <- sobol_summary |> filter(Category != "Control (Noise)") |> slice_max(order_by = median_Ti, n = 10) |> pull(Label) |> unique() sobol_site |> filter(Label %in% top10_labels, Category != "Control (Noise)") |> ggplot(aes(x = reorder(Label, Ti, FUN = median), y = Ti)) + geom_boxplot(aes(fill = Category), outlier.size = 0.8, alpha = 0.7) + coord_flip() + facet_wrap(~var_label, scales = "free_x", ncol = 2) + scale_fill_brewer(palette = "Set2") + labs( x = NULL, y = "Total-order index (Ti)", fill = "Source" ) + theme( axis.text.y = element_text(size = 9), strip.text = element_text(size = 11, face = "bold"), legend.text = element_text(size = 10) ) ``` # Interaction Analysis {#sec-interactions} Parameters with large differences between $T_i$ and $S_i$ interact strongly with other factors. These non-additive effects mean the parameter's impact on output depends on the values of other inputs. Interaction strength ($T_i - S_i$) indicates how much a parameter's influence depends on other inputs rather than acting independently. ```{r} #| label: fig-interactions #| fig-cap: "Interaction strength (Ti - Si) for inputs with interaction component exceeding 0.05 (i.e. 5% of total output variance, since Sobol indices are variance-normalized). Numbers show percent of Ti attributable to interactions." #| fig-height: 12 #| fig-width: 12 high_int <- sobol_summary |> filter(Category != "Control (Noise)", median_Int > 0.05) |> mutate(Int_pct = (median_Int / median_Ti) * 100) if (nrow(high_int) > 0) { ggplot(high_int, aes(x = reorder(Label, median_Int), y = median_Int)) + geom_col(aes(fill = Category), alpha = 0.8) + geom_text( aes(label = paste0(round(Int_pct, 0), "%")), hjust = -0.1, size = 4 ) + coord_flip() + facet_wrap(~var_label, scales = "free_x", ncol = 2) + scale_fill_brewer(palette = "Set2") + labs( x = NULL, y = "Interaction strength (Ti - Si)", fill = "Source" ) + theme( axis.text.y = element_text(size = 9), strip.text = element_text(size = 11, face = "bold"), legend.text = element_text(size = 10) ) } else { cat("No parameters exceeded the 5% interaction threshold.") } ``` # Factor Fixing (Screening) {#sec-screening} To identify non-influential parameters, we compare their Total-Order index ($T_i$) against the **Dummy parameter** (noise floor) that is useful for screening. The dummy parameter varies randomly but has no effect on model output, so we assume that any parameter with $T_i$ at or below the dummy's median can be fixed to a constant value without meaningful loss of variance coverage. The dashed red line marks the noise floor (median dummy $T_i$). ```{r} #| label: fig-screening #| fig-cap: "Factor screening: parameters below the noise line (dashed red) can be fixed. Dot = median Ti across sites; whisker = IQR." #| fig-height: 16 #| fig-width: 12 dummy_thresholds <- sobol_summary |> filter(Category == "Control (Noise)") |> select(variable, var_label, noise_Ti = median_Ti) screening <- sobol_summary |> filter(Category != "Control (Noise)") |> left_join(dummy_thresholds, by = c("variable", "var_label")) ggplot(screening, aes(x = reorder(Label, median_Ti), y = median_Ti)) + geom_point(aes(color = Category), size = 2.5) + geom_errorbar( aes(ymin = q25_Ti, ymax = q75_Ti, color = Category), width = 0.3, alpha = 0.6 ) + geom_hline( data = dummy_thresholds, aes(yintercept = noise_Ti), linetype = "dashed", color = "red", alpha = 0.7 ) + coord_flip() + facet_wrap(~var_label, scales = "free_x") + scale_color_brewer(palette = "Set2") + labs( x = NULL, y = "Total-order index (Ti)", color = "Source" ) + theme(axis.text.y = element_text(size = 8)) ``` ```{r} #| label: fixable-params #| results: asis fixable <- screening |> filter(median_Ti <= noise_Ti) |> mutate(Variable = label_variable(variable)) |> select(Variable, Input = Label, `Median Ti` = median_Ti, `Noise Threshold` = noise_Ti) |> arrange(Variable, `Median Ti`) if (nrow(fixable) > 0) { cat("\n**Fixable parameters** (Ti below noise floor):\n\n") knitr::kable(fixable, digits = 4) } else { cat("\n**No parameters met the screening criteria.** All inputs exceed the numerical noise floor.\n") } ``` # Model Additivity {#sec-additivity} The sum of First-Order indices ($\sum S_i$) indicates the degree of model linearity. Values near 1.0 imply the model is **additive** and parameters contribute independently. Values well below 1 indicate that **interactions dominate** -- the model's response to one parameter depends on the values of others. ```{r} #| label: fig-additivity #| fig-cap: "Distribution of model additivity (sum of Si) across sites. Dashed line at 1.0 marks perfect additivity. Values below 1 indicate interaction effects." #| fig-height: 4 #| fig-width: 7 additivity_site <- sobol_parsed |> filter(Category != "Control (Noise)") |> summarize(Sum_Si = sum(Si, na.rm = TRUE), .by = c("runid", "variable", "var_label")) ggplot(additivity_site, aes(x = var_label, y = Sum_Si)) + geom_boxplot(fill = "steelblue", alpha = 0.5, outlier.size = 0.8) + geom_hline(yintercept = 1, linetype = "dashed", color = "grey50") + labs(x = NULL, y = expression(sum(S[i]))) + coord_flip() ``` # Convergence Diagnostics {#sec-convergence} A Sobol analysis is only as reliable as its sample size allows. With $N = `r N`$ base samples and $k = `r k`$ parameters, our total model evaluation budget is `r total_runs` runs per site. This section assesses whether that budget was sufficient for stable index estimation, following the diagnostic framework of Nossent et al. (2011). ## Confidence Interval Precision The width of the bootstrap confidence interval ($\text{CI}_{95}$) for each index indicates estimation precision. Parameters with narrow CIs have well-determined importance; wide CIs flag indices that may shift with additional samples. ```{r} #| label: fig-ci-precision #| fig-cap: "Bootstrap 95% confidence interval width for Si (left) and Ti (right), by output variable. Each point is one parameter. Narrower intervals indicate more precise estimation. The dashed line marks a CI width of 0.1 (10% of output variance)." #| fig-height: 8 #| fig-width: 12 ci_data <- sobol_indices |> mutate( Si_ci_width = Si_high.ci - Si_low.ci, Ti_ci_width = Ti_high.ci - Ti_low.ci, var_label = label_variable(variable) ) p_si_ci <- ggplot(ci_data, aes(x = reorder(label_param(parameters), Si_ci_width, FUN = median), y = Si_ci_width)) + geom_boxplot(fill = "steelblue", alpha = 0.5, outlier.size = 0.5) + geom_hline(yintercept = 0.1, linetype = "dashed", color = "red", alpha = 0.5) + coord_flip() + labs(x = NULL, y = "Si CI width (95%)", title = "First-Order Index") + theme(axis.text.y = element_text(size = 7)) p_ti_ci <- ggplot(ci_data, aes(x = reorder(label_param(parameters), Ti_ci_width, FUN = median), y = Ti_ci_width)) + geom_boxplot(fill = "#D55E00", alpha = 0.5, outlier.size = 0.5) + geom_hline(yintercept = 0.1, linetype = "dashed", color = "red", alpha = 0.5) + coord_flip() + labs(x = NULL, y = "Ti CI width (95%)", title = "Total-Order Index") + theme(axis.text.y = element_text(size = 7)) p_si_ci + p_ti_ci ``` ## Parameter-Output Relationships {#sec-scatterplots} These scatterplots show the direct relationship between each top parameter's sampled value and the model output. They reveal whether the sensitivity is linear (validating the Sobol decomposition) or nonlinear (indicating threshold behavior or saturation). Each point is one model evaluation from the Sobol design matrix. ```{r} #| label: fig-param-output-scatter #| fig-cap: "Parameter-output scatterplots for the top 3 most influential parameters per output variable. Each point is one run from the Sobol design. Trend lines (LOESS) reveal linear vs. nonlinear relationships." #| fig-height: 14 #| fig-width: 12 #| warning: false # load design matrix design_matrix <- readr::read_csv( here::here(cfg$paths$data_dir, "sobol_design_matrix.csv"), show_col_types = FALSE ) # top 3 params per variable (by median Ti, excluding dummy/drivers) top3_per_var <- sobol_summary |> filter(!Category %in% c("Control (Noise)", "Initial Conditions", "Meteorology")) |> slice_max(order_by = median_Ti, n = 3, by = "variable") # load Y vectors for one representative site output_files <- list.files(here::here("output"), "ensemble\\.output.*\\.Rdata$", full.names = TRUE) representative_rid <- unique(sobol_indices$runid)[1] rid_files <- output_files[grepl(representative_rid, output_files)] scatter_list <- list() for (i in seq_len(nrow(top3_per_var))) { var_code <- top3_per_var$variable[i] var_pretty <- top3_per_var$var_label[i] param_label <- top3_per_var$Label[i] # find the raw parameter name in the design matrix # Label was derived from Trait, need to map back param_row <- sobol_parsed |> filter(variable == var_code, Label == param_label) |> slice_head(n = 1) if (nrow(param_row) == 0) next raw_param <- paste0(param_row$Prefix, ".", param_row$Trait) if (param_row$Prefix == param_row$Trait) raw_param <- param_row$Prefix if (!raw_param %in% names(design_matrix)) next # load Y yfile <- rid_files[grepl(paste0("\\.", var_code, "\\."), rid_files)] if (length(yfile) != 1) next env <- new.env(parent = emptyenv()) load(yfile, envir = env) Y <- as.numeric(unlist(env$ensemble.output)) if (length(Y) != nrow(design_matrix)) next scatter_list[[paste0(var_code, "_", raw_param)]] <- tibble::tibble( x = design_matrix[[raw_param]], y = Y, param = param_label, variable = var_pretty ) } if (length(scatter_list) > 0) { scatter_df <- dplyr::bind_rows(scatter_list) |> mutate(facet = paste0(variable, " ~ ", param)) ggplot(scatter_df, aes(x = x, y = y)) + geom_point(alpha = 0.15, size = 0.8, color = "steelblue") + geom_smooth(method = "loess", se = FALSE, color = "#D55E00", linewidth = 0.8) + facet_wrap(~facet, scales = "free", ncol = 3) + labs(x = "Parameter value (scaled)", y = "Model output") + theme(strip.text = element_text(size = 9)) } else { cat("Could not load design matrix or output files for scatterplots.\n") } ``` ::: {.callout-note} ## Interpreting Scatterplots - **Linear trend with tight spread**: parameter has a strong, additive effect -- Si will be high - **Nonlinear trend (curve, threshold)**: parameter has important interactions -- Ti >> Si - **Cloud with no trend**: parameter has low sensitivity for this variable - **Funnel shape (spread changes with x)**: heteroscedasticity -- the parameter's effect depends on other inputs ::: ## Convergence with Sample Size {#sec-convergence-n} To assess whether $N = `r N`$ is sufficient, we recompute Sobol indices using increasing fractions of the base sample. If indices stabilize at the full sample size, this supports adequate convergence; indices still drifting indicate that a larger N would improve precision. ```{r} #| label: fig-convergence #| fig-cap: "Evolution of Ti for the top 5 parameters as the Sobol base sample size increases. Stable plateaus indicate convergence; ongoing drift indicates the estimate has not fully stabilized." #| fig-height: 6 #| fig-width: 10 #| warning: false # subsample convergence for one representative site and variable conv_var <- "TotSoilCarb" conv_file <- rid_files[grepl(paste0("\\.", conv_var, "\\."), rid_files)] if (length(conv_file) == 1) { env <- new.env(parent = emptyenv()) load(conv_file, envir = env) Y_full <- as.numeric(unlist(env$ensemble.output)) # full parameter list including drivers (ic_ensemble, met_ensemble) dm_cols <- names(readr::read_csv( here::here(cfg$paths$data_dir, "sobol_design_matrix.csv"), n_max = 0, show_col_types = FALSE )) param_names <- dm_cols[dm_cols != "sample_id"] k_params <- length(param_names) # test at subset sizes n_vals <- c(8, 12, 16, 20, 24, 28, 32) n_vals <- n_vals[n_vals <= N] conv_results <- list() for (n_sub in n_vals) { n_total <- n_sub * (k_params + 2) if (n_total > length(Y_full)) next tryCatch({ idx <- sensobol::sobol_indices( Y = Y_full[1:n_total], N = n_sub, params = param_names, boot = FALSE ) df <- tibble::as_tibble(idx$results) |> filter(sensitivity == "Ti") |> mutate(N_sub = n_sub) conv_results[[as.character(n_sub)]] <- df }, error = function(e) NULL) } if (length(conv_results) > 0) { conv_df <- dplyr::bind_rows(conv_results) |> mutate(param_label = label_param(parameters)) # top 5 by final Ti top5_conv <- conv_df |> filter(N_sub == max(N_sub)) |> slice_max(order_by = original, n = 5) |> pull(parameters) conv_df |> filter(parameters %in% top5_conv) |> ggplot(aes(x = N_sub, y = original, color = param_label)) + geom_line(linewidth = 0.8) + geom_point(size = 2) + labs( x = "Base sample size (N)", y = paste("Ti for", label_variable(conv_var)), color = "Parameter" ) + scale_x_continuous(breaks = n_vals) + theme(legend.position = "bottom") } else { cat("Could not compute convergence subsets.\n") } } else { cat("Output file not found for convergence analysis.\n") } ``` ::: {.callout-note} ## Sample Size Considerations Our base sample $N = `r N`$ is modest by Sobol analysis standards (Nossent et al. 2011 recommend N = 500+ for stable individual indices), chosen to balance computational cost (~1,440 model runs per site x `r n_sites` sites) against the need for results across many sites. Where indices show ongoing drift at $N = `r N`$, a larger ensemble (Phase 4 workplan) would improve precision. The current results are reliable for parameter ranking and factor prioritization (relative ordering is robust), even though individual index point estimates carry estimation uncertainty. ::: ------------------------------------------------------------------------ # Synthesis **What drives uncertainty?** The variance decomposition (@sec-variance-partition) separates biological parameter uncertainty from environmental forcing. **What to calibrate first?** Parameter rankings (@sec-rankings) identify the top candidates for reducing predictive uncertainty. Parameters whose importance varies across sites (@fig-site-variability) need spatially explicit calibration. **What can be simplified?** Factor fixing (@sec-screening) uses the dummy parameter as a rigorous noise floor to identify inputs that can be set to nominal values. **Is the model additive?** Additivity diagnostics (@sec-additivity) indicate whether parameter effects are independent or interact. Strong interactions (large $T_i - S_i$ in @sec-interactions) suggest that calibrating parameters independently may be insufficient. # References 1. Puy, A., et al. (2022). sensobol: An R Package to Compute Variance-Based Sensitivity Indices. *Journal of Statistical Software*, 102(5). 2. Saltelli, A., et al. (2010). Variance based sensitivity analysis of model output. *Computer Physics Communications*. 3. Dietze, M. (2017). *Ecological Forecasting*. Princeton University Press. 4. Gelman, A. & Dodhia, R. (2002). Let's practice what we preach: Turning tables into graphs. 5. Nossent, J., Elsen, P. & Bauwens, W. (2011). Sobol' sensitivity analysis of a complex environmental model. *Environmental Modelling & Software*, 26, 1515-1525.