Implement multi-site local sensitivity analysis workflow

[divine7022]

Summary:

Implements a complete multi-site local sensitivity analysis (SA) workflow for SIPNET as specified in #151. This PR delivers aggregated sensitivity results across design points, environmental gradient analysis, and a Quarto-based report summarizing parameter importance patterns.

Addresses Issue #151 Requirements

[dlebauer]

BTW this looks really nice, well done. Though I'd like to work through the steps as I review.

Could you please add a readme with instructions on how to use the repository?

[divine7022]

thanks for reviewing, updated readme in #2

[infotroph]

Just realized I left these comments unposted on my partial readthrough the other day, so adding now. Would it be useful for me to finish a full review now or wait for updates first?

[infotroph]

Suggested change

	#----------------------------------
	# read.settings() uses xmlToList loses duplicate <variable> tags, so read XML directly
	library(XML)
	xml_doc <- XML::xmlParse(file.path("output", "pecan.CONFIGS.xml"))
	sa_variables <- unique(XML::xpathSApply(
	xml_doc,
	"//sensitivity.analysis//variable",
	XML::xmlValue
	))
	XML::free(xml_doc)
	sa_variables <- settings$sensitivity.analysis |>
	(\(x) x[names(x) == "variable"])() |>
	unlist() |>
	unique()

[divine7022]

👍

[infotroph]

Does this site selection process differ from what David already implemented in the downscaling repo? Seems preferable to use an existing tool rather than add that complexity to the UA process.

[divine7022]

I have refactored this script to completely remove the clustering logic. Instead it now acts purely as a consumer of the shared design_points.csv artifact (currently the david's manually selected 198 sites) to perform the necessary PFT mapping and sub-sampling for UA.

We will treat this as the interface for now, and can revisit or refactor the integration logic once the broader coordination strategy between the Downscaling and Uncertainty workflows is finalized.

see comment #1 (comment)

[infotroph]

I'm not sure what "local" means here

[divine7022]

good catch regarding the ambiguity

In this context, 'Local' refers to the mathematical method used by pecan's variance decomposition one-at-a-time (OAT) SA, where parameters are varied individually around their median to calculate partial derivatives (elasticities). This is distinct from the global sensitivity analysis we perform later.

I have updated the function documentation to explicitly state this aggregates one-at-a-time (OAT) parameter sensitivity results to avoid confusion with geographical locality.

[infotroph]

I recommend documenting what columns are expected in this df, and whether unexpected columns are ignored or cause errors

[divine7022]

👍 documented

[infotroph]

This will be much easier to read and understand if you define it as a function with an informative name -- yes, even if it's only called once as all_results <- purrr::map_dfr(sa_files, name_of_function)

[divine7022]

I have extracted the file processing logic into a standalone helper function process_sa_file

[divine7022]

Would it be useful for me to finish a full review now or wait for updates first?

Thanks for the comments,
@infotroph let me address these first, you can continue the full review once I push the next update

[divine7022]

.

[dlebauer]

Overall, this is very nice work. The code is clear and easy to understand.

I've read through the scripts and R functions. My next step will be to review the report. I've made a few comments. Not all need to be done now, but please write tickets or todos to capture future work.

[dlebauer]

The goal with this file was to handle settings that aren't in the pecan xml. Is there a reason that these are included here rather than the template.xml?

Minimizing config options here, and providing sensible defaults in the template.xml could make it more clear to end users.

[divine7022]

good point, template.xml already holds most pecan specific defaults. The items in 000-config.yml fall into two categories,
1 ) project level paths (ccmmf_dir, design_points, pecan_outdir) -- these are system specific and shouldn't live in template.xml since can change per machine

2 ) run configs (n_ensemble, n_met, start_date, end_date) -- moved these here from being hardcoded in
002_build_xml.R so users can configure them in one place. They feed into the XML building step

I have minimized config, let me know if you still have any specifics

[dlebauer]

Once we have ccmmf_dir, can we change this to use that as a variable? That way it is only necessary to change the system-specific path once. (also applies to master_design points etc. )

Suggested change

	pecan_outdir: "/projectnb2/dietzelab/ccmmf/modelout/ccmmf_phase_2b_mixed_pfts_20250701"
	pecan_outdir: "$ccmmf_dir/modelout/ccmmf_phase_2b_mixed_pfts_20250701"

[divine7022]

I looked into this, unfortunately yaml doesn't support variable substitution like $ccmmf_dir/modelout/...
If you'd prefer we could switch to an R based config

[dlebauer]

nit: these are just design points ...

Suggested change

	master_design_points: "/projectnb2/dietzelab/ccmmf/data/design_points.csv"
	design_points: "/projectnb2/dietzelab/ccmmf/data/design_points.csv"

[divine7022]

rename 👍

[dlebauer]

These may seem like trivial naming requests, but clarity and specificity is important.

• why is this 'master_design_points' ... that implies that there are other sets of design points.
• I would call this object design_points_csv instead of master_file
• instead of master_data, I would call that object design_points.

[dlebauer]

why are we sampling here?

[divine7022]

this i just to sample few site to test.
added a development/production mode

[dlebauer]

these were not manually selected, they were from the site selection

[divine7022]

updated comment to reference clustering

[dlebauer]

The fact that this script exists makes me wonder whether the clustering workflow that generates design_points.csv should be updated - that will help ensure consistency in workflows that consume it.

[divine7022]

yes, put mapping step (annual crop -> grass, woody perennial crop -> temperate.deciduous) should ideally be part of the clustering output; added a TODO to coordinate with downscaling

[dlebauer]

it is unclear - why are configs hard coded here if they are already in the config file and/or template.xml?

[divine7022]

moved n_ens, n_met, start_date, end_date to the run: section in 000-config.yml. The script now reads all configurable values from config. Paths like ic_dir and met_dir are derived from ccmmf_dir in the script

[dlebauer]

consider - what is the minimum amount of information that we want the end users to be able to configure, and why?

In addition to file paths, I can see n_sample, start_date, end_date being useful for being able to trigger 'development' mode.

[divine7022]

current user facing config:

• flags.production -- true/false (controls site count)
• sites.n_sample -- number of sites in dev mode
• run.start_date / run.end_date
• run.n_ensemble / run.n_met
• paths.ccmmf_dir
  everything else is derived

[dlebauer]

is it worth the effort and overhead to manage a STATUS file? My recollection is that the STATUS file was developed primarily for the PHP web interface to display. It seems that dropping the status file would simplify the code a lot. What benefit does it provide, if any?

This is a context where targets might be useful.

[divine7022]

keeping it for now, STATUS file enables the --continue flag for workflow resumption, which is valuable when running large site SA jobs; these checks lets us skip completed steps on restart without re running everything

agree that targets would be a cleaner approach for this. Added to TODO list for future work

[dlebauer]

Overall the report is a nice analysis. A few changes could make the results easier to interpret.

• Read and implement "turning tables into graphs"
• Lead with a 3–5 bullet summary of the top findings (what parameters are most important to calibrate, how do these vary w/ environmental drivers).
• display all of the elasticity plots at once using faceting.
• show variance across sites in elasticity plot
• Show sites on a map similar to the one in showing the design points in the downscaling repo
• Replace site table with a rank‑distribution plot (median + IQR) and link to downloadable CSV.
  Make gradient results a small multiple of partial dependence plots for top 3 parameters, rather than a long table.

[dlebauer]

It would be helpful to be able to view all of these plots at the same time, either as panels or different color bars.

By default, I prefer to show medians, especially for non-normally distributed values. If they exist, consider showing the distribution of values (e.g., across sites)

[divine7022]

done, replaced the panel-tabset with a single faceted plot showing all response variables simultaneously. points show median signed elasticity across sites; whiskers show IQR

[dlebauer]

The site-by-site table is difficult to interpret and doesn't add much once across site variance is shown in the elasticity plot. Lets remove it. Andrew Gelman’s “Let’s practice what we preach: turning tables into graphs” makes a good case for visual summaries when tables get this dense: https://sites.stat.columbia.edu/gelman/research/published/dodhia.pdf

(The reason that I put a dense table of county aggregated means in the downscaling report is that those values are the numbers that the client wants)

A map would be a nice way to show the distribution of sites. It could eventually (not required for MVP) be interesting to map elasticity of the first, say, five parameters.

A few small usability notes for the style of table:

• when rendered as HTML, the table scrolls horizontally but headers don’t stay aligned.
• Two significant figures likely be enough; additional precision is noise.
• The cell color bars are helpful for relative magnitude.

[dlebauer]

Consider something like a heat map. it isn't necessary to show every site, but perhaps summarize by environmental cluster?

See Dietze et al 2017 https://agupubs.onlinelibrary.wiley.com/doi/10.1002/2013JG002392

[divine7022]

replaced dense DT table with:

• CA site map (points colored by pff)
• boxplots showing the cross site distribution of elasticity for the top 10 parameters

CSV data remains available download

[dlebauer]

link to this section above, or perhaps it could be in the beginning.

[divine7022]

added cross reference link to @sec-technical in the overview callout note where sensitivity metrics are first mentioned

[dlebauer]

its not clear why there are only two unique values of r2 in the table

[divine7022]

duplicate values occurred because the elasticity and variance_explained regression targets were producing identical results (same underlying data, same gradient variables). Restructured the gradient section to focus on signed elasticity only, should eliminates the duplicate

[dlebauer]

this plot has two rows with the same values

[dlebauer]

one way these could be displayed is using facets

• each row a different parameter (for top params / greatest slopes)
• each column a different environmental covariate
• each plot has x = env covariate and y = elasticity. Could use color to plot multiple outputs (SOC, NPP, etc) on the same plot, or create separate plots for each output.

Another note on plots - by default, keep the axes (and other scales) the same size even across different variations of a plot (e.g. when all elasticity should have the same range. Sometimes it is necessary to change this role - e.g. to show interesting trends that would otherwise be obscured. This makes it easier to compare.

[divine7022]

Implemented as faceted scatter plots. Each panel shows one parameter response gradient combination with an OLS fit and R2 annotation. Filtered to the top 8 relationships by R2; and axes are consistent where possible

[dlebauer]

taking the absolute value of elasticity makes sense when ranking parameter importance, but for plotting and regression, we should use untransformed elasticity.

There are lots of ways to plot this information but this is from LeBauer et al 2013 (shown here mostly to note that the direction of elasticity is also important information).

[divine7022]

good point, added signed (untransformed) elasticity summary stats to parameter_rankings:

• mean_elasticity, median_elasticity (signed means for plotting)
• q25_elasticity, q75_elasticity (for IQR whiskers)

absolute elasticity (mean_abs_elasticity) is retained for ranking and priority scoring; all plots now use signed elasticity

[dlebauer]

Additional items for the "TODO" list

• split diagnostics into a separate doc from sensitivity results; report should focus on sensitivity metrics that are useful for model development, parameterization
• add plots for each output variable of parameter value vs output (i.e. continuous sensitivity, these should be automatically generated by PEcAn)
• consider grouping sites by cluster (from clustering workflow)
• consider hierarchical bayesian model using brms::brm with outputs, parameters, and sites as random effects.
• test assumptions of regression model
• do not report P-values. They are okay for filtering but not hypothesis testing, and the P-values here would need to be corrected for the fact that we are running lots of tests (a simple correction is p * number of times lm is called

[dlebauer]

I haven't completed my review, but am going to start focusing on GSA (#2) per discussion.

These are the comments/suggestions I have so far.

[dlebauer]

these are not manually selected - these were generated by an earlier version of the downscaling workflow.

To ensure that we are using the correct version of an external dataset (in this case, design_points.csv), we should not check a copy into this repository.

We should have a shared understanding of where to write data and

[dlebauer]

I think it would be more clear to put n_sample and any other values that will be used in dev mode into a separate 'development' chunk

[dlebauer]

This is a small point to improve clarity & maintainability:

• add select to initial changes to design_points
  • only needs to be done once
• then
  • if(production) just sends a message
  • if(dev) only subsets + sends message

Suggested change

	dplyr::mutate(
	pft = dplyr::case_when(
	pft == "annual crop" ~ "grass",
	pft == "woody perennial crop" ~ "temperate.deciduous",
	TRUE ~ NA_character_
	)
	) |>
	dplyr::filter(!is.na(pft))
	if (production) {
	# production: use all sites
	final_data <- design_points |>
	dplyr::select(site_id, lat, lon, pft)
	PEcAn.logger::logger.info(
	"Production mode: using all ", nrow(final_data), " sites"
	)
	} else {
	# development: random subset for fast testing
	n_sample <- cfg$sites$n_sample %||% 10L
	set.seed(42)
	final_data <- design_points |>
	dplyr::slice_sample(n = min(n_sample, nrow(design_points))) |>
	dplyr::select(site_id, lat, lon, pft)
	PEcAn.logger::logger.info(
	"Development mode: sampled ", nrow(final_data), " of ",
	nrow(design_points), " sites"
	)
	}
	readr::write_csv(final_data, out_file)
	dplyr::select(site_id, lat, lon, pft) |>
	dplyr::mutate(
	pft = dplyr::case_when(
	pft == "annual crop" ~ "grass",
	pft == "woody perennial crop" ~ "temperate.deciduous",
	TRUE ~ NA_character_
	)
	) |>
	dplyr::filter(!is.na(pft))
	num_sites <- nrow(design_points)
	if (production) { # use all sites
	PEcAn.logger::logger.info(
	"Production mode: using all ", num_sites, " sites"
	)
	} else { # random subset for fast testing
	n_sample <- cfg$sites$n_sample %||% 10L
	set.seed(42)
	design_points <- design_points |>
	dplyr::slice_sample(n = min(n_sample, nrow(design_points)))
	PEcAn.logger::logger.info(
	"Development mode: sampled ", nrow(final_data), " of ",
	num_sites, " sites"
	)
	}
	readr::write_csv(design_points, out_file)

[dlebauer]

Suggested change

	We use the PEcAn variance decomposition workflow, which proceeds in three steps:
	We use the PEcAn variance decomposition workflow (LeBauer et al, 2013), which proceeds in three steps:

[dlebauer]

Suggested change

	2. **Sensitivity analysis.** Each parameter is perturbed one at a time to the quantile equivalents of $\pm 1\sigma$ and $\pm 2\sigma$ of its posterior distribution, and the model is re-run. The resulting response curves are fit with splines, and elasticity is computed as the normalized sensitivity:
	2. **Sensitivity analysis.** The model is run first with all parameters are set to their posterior median values, and then with each parameter perturbed one at a time to the quantile equivalents of $\pm 1\sigma$ and $\pm 2\sigma$ of its posterior distribution, and the model is re-run. The resulting response curves are fit with splines, and elasticity is computed as the normalized sensitivity: