Skip to content

Add bootstrapping uncertainty and GPP/R constraints to metab.mle and metab.kalman - #154

Open
bmcafee wants to merge 15 commits into
GLEON:mainfrom
bmcafee:main
Open

Add bootstrapping uncertainty and GPP/R constraints to metab.mle and metab.kalman#154
bmcafee wants to merge 15 commits into
GLEON:mainfrom
bmcafee:main

Conversation

@bmcafee

@bmcafee bmcafee commented Jul 17, 2026

Copy link
Copy Markdown
Contributor

This updates adds new arguments to metab.mle and metab.kalman: constrain.sign and n.boot. The constrain.sign argument forces GPP and R to be positive and negative, respectively, via exponentiation. The n.boot argument adds new columns to the dataframe output describing the upper and lower 95% confidence intervals of the bootstrapped parameters as well as the coefficient of variation.

There is also a small change to the metab function adding the verbose argument. If TRUE, the verbose argument adds a progress bar to the console while analyzing a time series. Bootstrapping slows down the process significantly, so this is a way to confirm to the user that progress is being made.

Please see the added references for justification of each implements bootstrapping procedure. Of note, the uncertainty of GPP and R estimates when using metab.kalman with constrained parameters are much smaller than without constrained parameters. I am not certain what kind of interactions could cause this behavior. As far as I am aware, the procedure is accurate. I've added a warning message about the Kalman filter bootstrapping to let the user know to be wary.

Comment thread R/metab.kalman.R Outdated
Comment thread R/metab.kalman.R Outdated
Comment thread R/metab.R Outdated

@jzwart jzwart left a comment

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks @bmcafee! This is great. The bootstrapping approach is sound and nicely referenced, and the constrain.sign reparameterization is applied consistently across the NLL functions, the point estimates, and the bootstrap refits.

I have a few suggestions as in-line comments. It'd also might be nice to have some example code for CI output in demo/fig_metab.R, as well as some tests to cover the new functions. Let me know if you want help with either of those.

Comment thread R/metab.mle.R Outdated
n.obs <- length(do.obs)

if (ar1.resids){
ar1.lm <- lm(resids[1:(n.obs - 1)] ~ resids[2:n.obs] - 1)

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm pretty sure this is regressing the earlier residuals on the later ones (resids[t] ~ resids[t+1), but we'd want a forward AR(1), which is resids[t] ~ resids[t-1]

Suggested change
ar1.lm <- lm(resids[1:(n.obs - 1)] ~ resids[2:n.obs] - 1)
ar1.lm <- lm(resids[2:n.obs] ~ resids[1:(n.obs - 1)] - 1)

Comment thread R/metab.mle.R
if (ar1.resids) {
# Residual randomization preserving AR(1)
simRes <- rep(NA, n.obs)
simRes[1] <- sample(resids[!is.na(resids)], 1)

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here and elsewhere, the bootstraps call sample() / rnorm() with no seed handling, so the CIs change run-to-run and can't be tested deterministically. Could you add a seed argument threaded down from metab.mle / metab.kalman, and set it at the top of the bootstrap, or at minimum document that users should set.seed() beforehand if they want reproducible CI's?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've updated the documentation to mention set.seed in commit f081a57. Setting the seed in the global environment makes it work every time, and ensures that our function doesn't need to interact with that environment at all.

Comment thread R/metab.R Outdated
Comment thread R/metab.mle.R
Comment on lines +210 to +211
cv.R <- stats::sd(boot.results$R, na.rm=TRUE)/mean(boot.results$R, na.rm=TRUE) # coefficient of variation
cv.NEP <- stats::sd(boot.results$NEP, na.rm=TRUE)/mean(boot.results$NEP, na.rm=TRUE) # coefficient of variation

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

CV for NEP and R can be hard to interpret. NEP's mean can approach 0 and thus CV explodes, and R's mean is negative so the CV is negative. Maybe make a note in the return portion of the docs to let the users know what this output means.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've added these details to the documentation in commit f081a57.

@jzwart

jzwart commented Aug 4, 2026

Copy link
Copy Markdown
Member

re: the note about the constrained-Kalman intervals being surprisingly tight: I think that might be expected.

When coefficients are estimated in log‑space and then exponentiated, a wide spread in log‑space collapses into a narrow, right‑skewed range close to zero. As a result, the quantile‑based confidence interval on GPP (proportional to gppCoeff) becomes tighter. and in the unconstrained bootstrap, gppCoeff and rCoeff can cross zero across resamples, and that sign instability inflates the CI. constraining prevents sign changes, reducing variance and tightening the interval.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants