From 4339554046721a7568796e9abd8cb33b4ddafaa5 Mon Sep 17 00:00:00 2001 From: Bennett McAfee <105509324+bmcafee@users.noreply.github.com> Date: Fri, 12 Jun 2026 16:14:48 -0400 Subject: [PATCH 01/13] Add bootstrapping to metab.mle --- NAMESPACE | 2 + R/metab.R | 15 ++++- R/metab.mle.R | 160 +++++++++++++++++++++++++++++++++++++++++++---- man/metab.Rd | 6 +- man/metab.mle.Rd | 59 +++++++++++++---- 5 files changed, 212 insertions(+), 30 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index 9c78e66..86fceee 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -59,6 +59,8 @@ importFrom(stats,lm) importFrom(stats,median) importFrom(stats,model.matrix) importFrom(stats,optim) +importFrom(stats,quantile) +importFrom(stats,rnorm) importFrom(stats,sd) importFrom(utils,flush.console) useDynLib(LakeMetabolizer, .registration = TRUE) diff --git a/R/metab.R b/R/metab.R index d7585a5..aea9bc3 100644 --- a/R/metab.R +++ b/R/metab.R @@ -7,7 +7,7 @@ #' #' #'@usage -#'metab(data, method, wtr.name="wtr", irr.name="irr", do.obs.name="do.obs", ...) +#'metab(data, method, verbose = TRUE, wtr.name="wtr", irr.name="irr", do.obs.name="do.obs", ...) #' #'@param data #' a data.frame whose columns are @@ -31,6 +31,7 @@ #'@param method #' a character string specifying one of the 5 statistical methods #'(bayesian, bookkeep, kalman, ols, mle) +#'@param verbose logical. If TRUE, a progress bar will be displayed. Defaults to TRUE #'@param wtr.name the name of the column containing temperature at the depth of do.obs (predictor variable for R) #'@param irr.name the name of the column containing irradiance (predictor variable for GPP) #'@param do.obs.name the name of the column in data containing the DO observations (in mg/L) to be used as the response variable @@ -59,7 +60,7 @@ #' #' To calculate k.gas: \link{k600.2.kGAS} #' -#' To calculate k600 values for k.gas: \link{k.cole}, \link{k.crusius}, \link{k.macIntyre}, \link{k.read} +#' To calculate k600 values for k.gas: \link{k.cole}, \link{k.crusius}, \link{k.macIntyre}, \link{k.read}, \link{k.klaus} #' #' #'@examples @@ -107,7 +108,7 @@ #' #'@export -metab <- function(data, method = NULL, wtr.name="wtr", irr.name="irr", do.obs.name="do.obs", ...){ +metab <- function(data, method = NULL, verbose = TRUE, wtr.name="wtr", irr.name="irr", do.obs.name="do.obs", ...){ m.args <- list(...) @@ -163,6 +164,9 @@ metab <- function(data, method = NULL, wtr.name="wtr", irr.name="irr", do.obs.na # ================================== # = Apply metab to subsets of data = # ================================== + if (verbose == TRUE){ + pb <- utils::txtProgressBar(min = 0, max = nid, initial = 0, style = 3) + } for(i in unique(ids)){ poss.args <- c("do.obs","do.sat","k.gas","z.mix", "irr", "wtr", "datetime") # data2 columns that could correspond to arguments @@ -173,7 +177,12 @@ metab <- function(data, method = NULL, wtr.name="wtr", irr.name="irr", do.obs.na # print(paste("Analyzing day #", i)); flush.console(); # Is this annoying? I'm commenting-out results[[i]] <- do.call(mtdCall, largs) # this is where all of the work happens + + if (verbose == TRUE){ + utils::setTxtProgressBar(pb,i) + } } + if (verbose == TRUE){close(pb)} answer0 <- conquerList(results, naming=data.frame("year"=data2[!duplicated(ids),"year"], "doy"=trunc(data2[!duplicated(ids),"doy"]))) diff --git a/R/metab.mle.R b/R/metab.mle.R index 754befa..8743a8f 100644 --- a/R/metab.mle.R +++ b/R/metab.mle.R @@ -1,6 +1,6 @@ #'@title Metabolism calculated from the maximum likelihood estimates of the parameters in a standard linear regression model #'@description Process-error-only model with parameters fitted via maximum likelihood estimation (MLE). This function runs the maximum likelihood metabolism model on the supplied gas concentration and other supporting data. -#'@param do.obs Vector of dissolved oxygen concentration observations, \eqn{mg O[2] L^{-1}}{mg O2 / L} +#'@param do.obs Vector of dissolved oxygen concentration observations, \eqn{mg O_2 L^{-1}}{mg O2 / L} #'@param do.sat Vector of dissolved oxygen saturation values based on water temperature. Calculate using \link{o2.at.sat} #'@param k.gas Vector of kGAS values calculated from any of the gas flux models #'(e.g., \link{k.cole}) and converted to kGAS using \link{k600.2.kGAS} @@ -8,6 +8,8 @@ #'@param irr Vector of photosynthetically active radiation in \eqn{\mu mol\ m^{-2} s^{-1}}{micro mols / m^2 / s} #'@param wtr Vector of water temperatures in \eqn{^{\circ}C}{degrees C}. Used in scaling respiration with temperature #'@param error.type Option specifying if model should assume pure Process Error 'PE' or Observation Error 'OE'. Defaults to observation error 'OE'. +#'@param n.boot Numeric of the number of bootstrap iterations. Set to 0 to bypass uncertainty quantification (default). Set to 2 or greater to quantify uncertainty via bootstrapping. +#'@param ar1.resids Logical. If TRUE, the AR(1) structure of the residuals will be retained during randomization when bootstrapping. #'@param ... additional arguments; currently "datetime" is the only recognized argument passed through \code{...} #'@return #'A data.frame with columns corresponding to components of metabolism @@ -16,6 +18,18 @@ #'\item{R}{numeric estimate of Respiration, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} #'\item{NEP}{numeric estimate of Net Ecosystem production, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} #'} +#'If n.boot is 2 or greater, the data frame will also contain the following columns: +#'\describe{ +#'\item{GPP.lci}{lower 95\% confidence interval of the estimate of Gross Primary Production, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} +#'\item{GPP.uci}{upper 95\% confidence interval of the estimate of Gross Primary Production, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} +#'\item{GPP.cv}{coefficient of variation of the bootstrapped estimates of Gross Primary Production} +#'\item{R.lci}{lower 95\% confidence interval of the estimate of Respiration, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} +#'\item{R.uci}{upper 95\% confidence interval of the estimate of Respiration, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} +#'\item{R.cv}{coefficient of variation of the bootstrapped estimates of Respiration} +#'\item{NEP.lci}{lower 95\% confidence interval of the estimate of Net Ecosystem Production, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} +#'\item{NEP.uci}{upper 95\% confidence interval of the estimate of Net Ecosystem Production, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} +#'\item{NEP.cv}{coefficient of variation of the bootstrapped estimates of Net Ecosystem Production} +#'} #' The maximum likelihood estimates of model parameters can be accessed via \code{attributes(metab.mle(...))[["params"]]} #' #'@details @@ -46,19 +60,27 @@ #'NEP is the sum of GPP and R. #' #'@note Currently, missing values in any arguments will result in an error, so freq must always equal nobs. -#'@author Luke A Winslow, Ryan Batt, GLEON Fellows +#'@author Luke A Winslow, Ryan Batt, Bennett McAfee, GLEON Fellows #'@references -#'Hanson, PC, SR Carpenter, N Kimura, C Wu, SP Cornelius, TK Kratz. 2008 -#'\emph{Evaluation of metabolism models for free-water dissolved oxygen in lakes}. -#'Limnology and Oceanography: Methods 6: 454:465. +#'Hanson, Paul C., Stephen R. Carpenter, Nobuaki Kimura, Chin Wu, Sean P. Cornelius, +#'and Timothy K. Kratz. 2008. \emph{Evaluation of Metabolism Models for Free-Water Dissolved +#'Oxygen Methods in Lakes}. Limnology and Oceanography: Methods 6 (9): 454–65. +#'\url{https://doi.org/10.4319/lom.2008.6.454}. +#' +#'Solomon, Christopher T., Denise A. Bruesewitz, David C. Richardson, et al. 2013. +#'\emph{Ecosystem Respiration: Drivers of Daily Variability and Background Respiration in Lakes around the Globe}. +#'Limnology and Oceanography 58 (3): 849–66. \url{https://doi.org/10.4319/lo.2013.58.3.0849}. +#' +#'Dugan, Hilary A., R. Iestyn Woolway, Arianto B. Santoso, et al. 2016. +#'\emph{Consequences of Gas Flux Model Choice on the Interpretation of Metabolic +#'Balance across 15 Lakes}. Inland Waters 6 (4): 581–92. \url{https://doi.org/10.1080/IW-6.4.836}. #' -#'Solomon CT, DA Bruesewitz, DC Richardson, KC Rose, MC Van de Bogert, PC Hanson, TK Kratz, B Larget, -#'R Adrian, B Leroux Babin, CY Chiu, DP Hamilton, EE Gaiser, S Hendricks, V Istvanovics, A Laas, DM O'Donnell, -#'ML Pace, E Ryder, PA Staehr, T Torgersen, MJ Vanni, KC Weathers, G Zhuw. 2013. -#'\emph{Ecosystem Respiration: Drivers of Daily Variability and Background Respiration in Lakes around the Globe}. -#'Limnology and Oceanography 58 (3): 849:866. doi:10.4319/lo.2013.58.3.0849. +#'Corman, Jessica R., Jacob A. Zwart, Jennifer Klug, et al. 2023. \emph{Response of +#'Lake Metabolism to Catchment Inputs Inferred Using High-Frequency Lake and +#'Stream Data from across the Northern Hemisphere}. Limnology and Oceanography +#'68 (12): 2617–31. \url{https://doi.org/10.1002/lno.12449}. #' -#'@importFrom stats dnorm optim +#'@importFrom stats dnorm optim rnorm quantile #' #'@seealso #'\link{metab}, \link{metab.bookkeep}, \link{metab.ols}, \link{metab.kalman}, \link{metab.bayesian} @@ -87,7 +109,7 @@ #' #'metab.mle(doobs[,2], do.sat, k.gas, z.mix[,2], irr[,2], wtr[,3]) #'@export -metab.mle <- function(do.obs, do.sat, k.gas, z.mix, irr, wtr, error.type="OE", ...){ +metab.mle <- function(do.obs, do.sat, k.gas, z.mix, irr, wtr, error.type="OE", n.boot = 0, ar1.resids = TRUE, ...){ complete.inputs(do.obs=do.obs, do.sat=do.sat, k.gas=k.gas, z.mix=z.mix, irr=irr, wtr=wtr, error=TRUE) @@ -159,7 +181,29 @@ metab.mle <- function(do.obs, do.sat, k.gas, z.mix, irr, wtr, error.type="OE", . GPP <- mean(pars[1]*irr, na.rm=TRUE) * freq R <- mean(pars[2]*log(wtr), na.rm=TRUE) * freq - return(list("params"=pars, "metab"=c("GPP"=GPP,"R"=R,"NEP"=GPP+R))) + # ================================= + # = Bootstrapping for uncertainty = + # ================================= + if (n.boot == 0){ + results <- list("params"=pars, "metab"=c("GPP"=GPP,"R"=R,"NEP"=GPP+R)) + } else if (n.boot >= 2){ + boot.results <- bootstrap.mle(n.boot, do.obs, pars, guesses, do.sat, k.gas, freq, z.mix, irr, wtr, error.type, ar1.resids) + ci.GPP <- stats::quantile(boot.results$GPP, c(0.025, 0.975), na.rm=TRUE) # 95% confidence interval + ci.R <- stats::quantile(boot.results$R, c(0.025, 0.975), na.rm=TRUE) # 95% confidence interval + ci.NEP <- stats::quantile(boot.results$NEP, c(0.025, 0.975), na.rm=TRUE) # 95% confidence interval + cv.GPP <- stats::sd(boot.results$GPP, na.rm=TRUE)/mean(boot.results$GPP, na.rm=TRUE) # coefficient of variation + 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 + results <- list("params"=pars, + "metab"=c("GPP"=GPP,"R"=R,"NEP"=GPP+R, + "GPP.lci" = unname(ci.GPP[1]), "GPP.uci" = unname(ci.GPP[2]), "GPP.cv" = cv.GPP, + "R.lci" = unname(ci.R[1]), "R.uci" = unname(ci.R[2]), "R.cv" = cv.R, + "NEP.lci" = unname(ci.NEP[1]), "NEP.uci" = unname(ci.NEP[2]), "NEP.cv" = cv.NEP)) + } else { + stop("n.boot must be 0 (no uncertainty quantification) or greater than 2 (uncertainty quantification via bootstrapping).") + } + + return(results) } # ============================================ @@ -226,3 +270,93 @@ mleNllOE <- function(Params, do.obs, do.sat, k.gas, z.mix, irr, wtr, error.type) return(-sum(dnorm(do.obs, alpha, sd=sqrt(Q), log=TRUE), na.rm=TRUE)) }#End function + +# ========================================= +# = Function to bootstrap for uncertainty = +# ========================================= +bootstrap.mle <- function(n.boot, do.obs, pars, guesses, do.sat, k.gas, freq, z.mix, irr, wtr, error.type = "OE", ar1.resids = FALSE){ + + ## Calculate fitted values and residuals + if (error.type == "OE"){ + alpha <- mleLoopOE(alpha = c(do.obs[1], rep(0, length(do.obs) - 1)), + doobs = do.obs, + c1 = pars[1], + c2 = pars[2], + beta = exp(-(k.gas / freq / z.mix)), + irr = irr, + wtr = wtr, + kz = k.gas / freq / z.mix, + dosat = do.sat) + } else if (error.type == "PE"){ + alpha <- mleLoopPE(alpha = c(do.obs[1], rep(0, length(do.obs) - 1)), + doobs = do.obs, + c1 = pars[1], + c2 = pars[2], + beta = exp(-(k.gas / freq / z.mix)), + irr = irr, + wtr = wtr, + kz = k.gas / freq / z.mix, + dosat = do.sat) + } else { + stop("error.type must be either 'OE' or 'PE', Observation Error or Process Error respectively.") + } + + + resids <- do.obs - alpha + n.obs <- length(do.obs) + + if (ar1.resids){ + ar1.lm <- lm(resids[1:(n.obs - 1)] ~ resids[2:n.obs] - 1) + ar1.coeff <- as.numeric(ar1.lm$coefficients) + ar1.sd <- sd(ar1.lm$residuals) + } + + ## Initialize results object + boot.results <- data.frame(boot.iter = seq_len(n.boot), + gppCoeff = rep(NA, n.boot), + rCoeff = rep(NA, n.boot), + Q = rep(NA, n.boot), + doInit = rep(NA, n.boot), + convergence = rep(NA, n.boot), + nll = rep(NA, n.boot), + GPP = rep(NA, n.boot), + R = rep(NA, n.boot), + NEP = rep(NA, n.boot)) + + ## Bootstrap loop + for (i in boot.results$boot.iter){ + + if (ar1.resids) { + # Residual randomization preserving AR(1) + simRes <- rep(NA, n.obs) + simRes[1] <- sample(resids[!is.na(resids)], 1) + for (j in 2:n.obs){simRes[j] <- ar1.coeff * simRes[j - 1] + stats::rnorm(1, sd = ar1.sd)} + } else { + # Residual randomization without replacement + simRes <- sample(resids[!is.na(resids)], length(resids), replace = FALSE) + } + + doSim <- alpha + simRes + + if (error.type == "OE"){ + simFit <- optim(guesses, fn = mleNllOE, do.obs = doSim, do.sat = do.sat, k.gas = k.gas / freq, z.mix = z.mix, irr = irr, wtr = wtr) + } else if (error.type == "PE"){ + simFit <- optim(guesses, fn = mleNllPE, do.obs = doSim, do.sat = do.sat, k.gas = k.gas / freq, z.mix = z.mix, irr = irr, wtr = wtr) + } + simFitPar <- simFit$par + simFitGPP <- mean(exp(simFitPar[1]) * irr, na.rm = TRUE) * freq + simFitR <- mean(-exp(simFitPar[2]) * log(wtr), na.rm = TRUE) * freq + + boot.results[i, "gppCoeff"] <- exp(simFitPar[1]) + boot.results[i, "rCoeff"] <- -exp(simFitPar[2]) + boot.results[i, "Q"] <- exp(simFitPar[3]) + boot.results[i, "doInit"] <- simFitPar[4] + boot.results[i, "convergence"] <- simFit$convergence + boot.results[i, "nll"] <- simFit$value + boot.results[i, "GPP"] <- simFitGPP + boot.results[i, "R"] <- simFitR + boot.results[i, "NEP"] <- simFitGPP + simFitR + } + + return(boot.results) +} diff --git a/man/metab.Rd b/man/metab.Rd index 27012ad..6a8c256 100644 --- a/man/metab.Rd +++ b/man/metab.Rd @@ -4,7 +4,7 @@ \alias{metab} \title{Calculate metabolism} \usage{ -metab(data, method, wtr.name="wtr", irr.name="irr", do.obs.name="do.obs", ...) +metab(data, method, verbose = TRUE, wtr.name="wtr", irr.name="irr", do.obs.name="do.obs", ...) } \arguments{ \item{data}{a data.frame whose columns are @@ -28,6 +28,8 @@ Columns that are not used by a particular statistical method do not need to be s \item{method}{a character string specifying one of the 5 statistical methods (bayesian, bookkeep, kalman, ols, mle)} +\item{verbose}{logical. If TRUE, a progress bar will be displayed. Defaults to TRUE} + \item{wtr.name}{the name of the column containing temperature at the depth of do.obs (predictor variable for R)} \item{irr.name}{the name of the column containing irradiance (predictor variable for GPP)} @@ -101,7 +103,7 @@ To calculate do.sat: \link{o2.at.sat} To calculate k.gas: \link{k600.2.kGAS} -To calculate k600 values for k.gas: \link{k.cole}, \link{k.crusius}, \link{k.macIntyre}, \link{k.read} +To calculate k600 values for k.gas: \link{k.cole}, \link{k.crusius}, \link{k.macIntyre}, \link{k.read}, \link{k.klaus} } \author{ Ryan D. Batt diff --git a/man/metab.mle.Rd b/man/metab.mle.Rd index 22fb799..c50bebb 100644 --- a/man/metab.mle.Rd +++ b/man/metab.mle.Rd @@ -4,10 +4,21 @@ \alias{metab.mle} \title{Metabolism calculated from the maximum likelihood estimates of the parameters in a standard linear regression model} \usage{ -metab.mle(do.obs, do.sat, k.gas, z.mix, irr, wtr, error.type = "OE", ...) +metab.mle( + do.obs, + do.sat, + k.gas, + z.mix, + irr, + wtr, + error.type = "OE", + n.boot = 0, + ar1.resids = TRUE, + ... +) } \arguments{ -\item{do.obs}{Vector of dissolved oxygen concentration observations, \eqn{mg O[2] L^{-1}}{mg O2 / L}} +\item{do.obs}{Vector of dissolved oxygen concentration observations, \eqn{mg O_2 L^{-1}}{mg O2 / L}} \item{do.sat}{Vector of dissolved oxygen saturation values based on water temperature. Calculate using \link{o2.at.sat}} @@ -22,6 +33,10 @@ metab.mle(do.obs, do.sat, k.gas, z.mix, irr, wtr, error.type = "OE", ...) \item{error.type}{Option specifying if model should assume pure Process Error 'PE' or Observation Error 'OE'. Defaults to observation error 'OE'.} +\item{n.boot}{Numeric of the number of bootstrap iterations. Set to 0 to bypass uncertainty quantification (default). Set to 2 or greater to quantify uncertainty via bootstrapping.} + +\item{ar1.resids}{Logical. If TRUE, the AR(1) structure of the residuals will be retained during randomization when bootstrapping.} + \item{...}{additional arguments; currently "datetime" is the only recognized argument passed through \code{...}} } \value{ @@ -31,6 +46,18 @@ A data.frame with columns corresponding to components of metabolism \item{R}{numeric estimate of Respiration, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} \item{NEP}{numeric estimate of Net Ecosystem production, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} } +If n.boot is 2 or greater, the data frame will also contain the following columns: +\describe{ +\item{GPP.lci}{lower 95\% confidence interval of the estimate of Gross Primary Production, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} +\item{GPP.uci}{upper 95\% confidence interval of the estimate of Gross Primary Production, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} +\item{GPP.cv}{coefficient of variation of the bootstrapped estimates of Gross Primary Production} +\item{R.lci}{lower 95\% confidence interval of the estimate of Respiration, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} +\item{R.uci}{upper 95\% confidence interval of the estimate of Respiration, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} +\item{R.cv}{coefficient of variation of the bootstrapped estimates of Respiration} +\item{NEP.lci}{lower 95\% confidence interval of the estimate of Net Ecosystem Production, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} +\item{NEP.uci}{upper 95\% confidence interval of the estimate of Net Ecosystem Production, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} +\item{NEP.cv}{coefficient of variation of the bootstrapped estimates of Net Ecosystem Production} +} The maximum likelihood estimates of model parameters can be accessed via \code{attributes(metab.mle(...))[["params"]]} } \description{ @@ -92,19 +119,27 @@ do.sat = o2.at.sat.base(wtr[,3], altitude=300) metab.mle(doobs[,2], do.sat, k.gas, z.mix[,2], irr[,2], wtr[,3]) } \references{ -Hanson, PC, SR Carpenter, N Kimura, C Wu, SP Cornelius, TK Kratz. 2008 -\emph{Evaluation of metabolism models for free-water dissolved oxygen in lakes}. -Limnology and Oceanography: Methods 6: 454:465. - -Solomon CT, DA Bruesewitz, DC Richardson, KC Rose, MC Van de Bogert, PC Hanson, TK Kratz, B Larget, -R Adrian, B Leroux Babin, CY Chiu, DP Hamilton, EE Gaiser, S Hendricks, V Istvanovics, A Laas, DM O'Donnell, -ML Pace, E Ryder, PA Staehr, T Torgersen, MJ Vanni, KC Weathers, G Zhuw. 2013. -\emph{Ecosystem Respiration: Drivers of Daily Variability and Background Respiration in Lakes around the Globe}. -Limnology and Oceanography 58 (3): 849:866. doi:10.4319/lo.2013.58.3.0849. +Hanson, Paul C., Stephen R. Carpenter, Nobuaki Kimura, Chin Wu, Sean P. Cornelius, +and Timothy K. Kratz. 2008. \emph{Evaluation of Metabolism Models for Free-Water Dissolved +Oxygen Methods in Lakes}. Limnology and Oceanography: Methods 6 (9): 454–65. +\url{https://doi.org/10.4319/lom.2008.6.454}. + +Solomon, Christopher T., Denise A. Bruesewitz, David C. Richardson, et al. 2013. +\emph{Ecosystem Respiration: Drivers of Daily Variability and Background Respiration in Lakes around the Globe}. +Limnology and Oceanography 58 (3): 849–66. \url{https://doi.org/10.4319/lo.2013.58.3.0849}. + +Dugan, Hilary A., R. Iestyn Woolway, Arianto B. Santoso, et al. 2016. +\emph{Consequences of Gas Flux Model Choice on the Interpretation of Metabolic +Balance across 15 Lakes}. Inland Waters 6 (4): 581–92. \url{https://doi.org/10.1080/IW-6.4.836}. + +Corman, Jessica R., Jacob A. Zwart, Jennifer Klug, et al. 2023. \emph{Response of +Lake Metabolism to Catchment Inputs Inferred Using High-Frequency Lake and +Stream Data from across the Northern Hemisphere}. Limnology and Oceanography +68 (12): 2617–31. \url{https://doi.org/10.1002/lno.12449}. } \seealso{ \link{metab}, \link{metab.bookkeep}, \link{metab.ols}, \link{metab.kalman}, \link{metab.bayesian} } \author{ -Luke A Winslow, Ryan Batt, GLEON Fellows +Luke A Winslow, Ryan Batt, Bennett McAfee, GLEON Fellows } From 9c766a28a334e365e80d22ac6c10ea1f6bf497d2 Mon Sep 17 00:00:00 2001 From: Bennett McAfee <105509324+bmcafee@users.noreply.github.com> Date: Tue, 30 Jun 2026 13:50:33 -0400 Subject: [PATCH 02/13] Add bootstrap.kalman --- R/metab.kalman.R | 559 +++++++++++++++++++++++++++----------------- man/metab.kalman.Rd | 21 +- 2 files changed, 366 insertions(+), 214 deletions(-) diff --git a/R/metab.kalman.R b/R/metab.kalman.R index 6eab269..2565a12 100644 --- a/R/metab.kalman.R +++ b/R/metab.kalman.R @@ -8,6 +8,7 @@ #'@param z.mix Vector of mixed-layer depths in meters. To calculate, see \link[rLakeAnalyzer]{ts.meta.depths} #'@param irr Vector of photosynthetically active radiation in \eqn{\mu mol\ m^{-2} s^{-1}}{micro mols / m^2 / s} #'@param wtr Vector of water temperatures in \eqn{^{\circ}C}{degrees C}. Used in scaling respiration with temperature +#'@param n.boot Numeric of the number of bootstrap iterations. Set to 0 to bypass uncertainty quantification (default). Set to 2 or greater to quantify uncertainty via bootstrapping. #'@param ... additional arguments; currently "datetime" is the only recognized argument passed through \code{...} #'@return #'A data.frame with columns corresponding to components of metabolism @@ -16,6 +17,18 @@ #'\item{R}{numeric estimate of Respiration, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} #'\item{NEP}{numeric estimate of Net Ecosystem production, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} #'} +#'If n.boot is 2 or greater, the data frame will also contain the following columns: +#'\describe{ +#'\item{GPP.lci}{lower 95\% confidence interval of the estimate of Gross Primary Production, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} +#'\item{GPP.uci}{upper 95\% confidence interval of the estimate of Gross Primary Production, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} +#'\item{GPP.cv}{coefficient of variation of the bootstrapped estimates of Gross Primary Production} +#'\item{R.lci}{lower 95\% confidence interval of the estimate of Respiration, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} +#'\item{R.uci}{upper 95\% confidence interval of the estimate of Respiration, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} +#'\item{R.cv}{coefficient of variation of the bootstrapped estimates of Respiration} +#'\item{NEP.lci}{lower 95\% confidence interval of the estimate of Net Ecosystem Production, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} +#'\item{NEP.uci}{upper 95\% confidence interval of the estimate of Net Ecosystem Production, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} +#'\item{NEP.cv}{coefficient of variation of the bootstrapped estimates of Net Ecosystem Production} +#'} #' #' #' Use \link{attributes} to access more model output: @@ -51,10 +64,13 @@ #' #'@references #'Batt, Ryan D. and Stephen R. Carpenter. 2012. \emph{Free-water lake metabolism: -#'addressing noisy time series with a Kalman filter}. Limnology and Oceanography: Methods 10: 20-30. doi: 10.4319/lom.2012.10.20 +#'addressing noisy time series with a Kalman filter}. Limnology and Oceanography: Methods 10: 20-30. \url{https://doi.org/10.4319/lom.2012.10.20} +#' +#'Stoffer, David S., and Kent D. Wall. 1991. \emph{Bootstrapping State-Space Models: Gaussian Maximum Likelihood Estimation and the Kalman Filter}. +#'Journal of the American Statistical Association 86 (416): 1024–33. \url{https://doi.org/10.2307/2290521}. #'@seealso #'\link{temp.kalman}, \link{watts.in}, \link{metab}, \link{metab.bookkeep}, \link{metab.ols}, \link{metab.mle}, \link{metab.bayesian} -#'@author Ryan Batt, Luke A. Winslow +#'@author Ryan Batt, Luke A. Winslow, Bennett McAfee #'@note If observation error is substantial, consider applying a Kalman filter to the water temperature time series by supplying #' \code{wtr} as the output from \link{temp.kalman} #'@examples @@ -84,22 +100,22 @@ #' do.sat=do.sat, wtr=wtr[,2], #' k.gas=k.gas, do.obs=doobs[,2]) #'@export -metab.kalman <- function(do.obs, do.sat, k.gas, z.mix, irr, wtr, ...){ - +metab.kalman <- function(do.obs, do.sat, k.gas, z.mix, irr, wtr, n.boot = 0, ...){ + complete.inputs(do.obs=do.obs, do.sat=do.sat, k.gas=k.gas, z.mix=z.mix, irr=irr, wtr=wtr, error=TRUE) - + nobs <- length(do.obs) - + mm.args <- list(...) - + if(any(z.mix <= 0)){ stop("z.mix must be greater than zero.") } if(any(wtr <= 0)){ stop("all wtr must be positive.") } - + if("datetime"%in%names(mm.args)){ # check to see if datetime is in the ... args datetime <- mm.args$datetime # extract datetime freq <- calc.freq(datetime) # calculate sampling frequency from datetime @@ -113,22 +129,46 @@ metab.kalman <- function(do.obs, do.sat, k.gas, z.mix, irr, wtr, ...){ # warning will only be seen through direct use of metab.bookkeep when datettime is not supplied freq <- nobs } - + # Filter and fit guesses <- c(1E-4,1E-4,log(5),log(5)) fit <- optim(guesses, fn=KFnllDO, do.obs=do.obs, do.sat=do.sat, k.gas=(k.gas/freq), z.mix=z.mix, irr=irr, wtr=wtr) pars0 <- fit$par pars <- c("gppCoeff"=pars0[1], "rCoeff"=pars0[2], "Q"=exp(pars0[3]), "H"=exp(pars0[4])) - + # Smooth - smoothDO <- KFsmoothDO(pars, do.obs=do.obs, do.sat=do.sat, k.gas=(k.gas/freq), z.mix=z.mix, irr=irr, wtr=wtr) - + KFresults <- KFsmoothDO(pars, do.obs=do.obs, do.sat=do.sat, k.gas=(k.gas/freq), z.mix=z.mix, irr=irr, wtr=wtr) + smoothDO <- KFresults$smoothDO + # Use fits to calculate metabolism - GPP <- mean(pars[1]*irr, na.rm=TRUE) * freq R <- mean(pars[2]*log(wtr), na.rm=TRUE) * freq - - return(list("smoothDO"=smoothDO,"params"=pars, "metab"=c("GPP"=GPP,"R"=R, "NEP"=GPP+R))) + + if (n.boot == 0){ + + return(list("smoothDO"=smoothDO,"params"=pars, "metab"=c("GPP"=GPP,"R"=R, "NEP"=GPP+R))) + + } else if (n.boot >= 2){ + + boot.results <- bootstrap.kalman(n.boot, Params = pars, KFresults, guesses, do.obs, do.sat, k.gas =(k.gas/freq), freq, z.mix, irr, wtr) + ci.GPP <- stats::quantile(boot.results$GPP, c(0.025, 0.975), na.rm=TRUE) # 95% confidence interval + ci.R <- stats::quantile(boot.results$R, c(0.025, 0.975), na.rm=TRUE) # 95% confidence interval + ci.NEP <- stats::quantile(boot.results$NEP, c(0.025, 0.975), na.rm=TRUE) # 95% confidence interval + cv.GPP <- stats::sd(boot.results$GPP, na.rm=TRUE)/mean(boot.results$GPP, na.rm=TRUE) # coefficient of variation + 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 + results <- list("smoothDO"=smoothDO, + "params"=pars, + "metab"=c("GPP"=GPP,"R"=R,"NEP"=GPP+R, + "GPP.lci" = unname(ci.GPP[1]), "GPP.uci" = unname(ci.GPP[2]), "GPP.cv" = cv.GPP, + "R.lci" = unname(ci.R[1]), "R.uci" = unname(ci.R[2]), "R.cv" = cv.R, + "NEP.lci" = unname(ci.NEP[1]), "NEP.uci" = unname(ci.NEP[2]), "NEP.cv" = cv.NEP)) + + } else { + stop("n.boot must be 0 (no uncertainty quantification) or greater than 2 (uncertainty quantification via bootstrapping).") + } + + } @@ -139,79 +179,79 @@ metab.kalman <- function(do.obs, do.sat, k.gas, z.mix, irr, wtr, ...){ # ====================== # Main recursion written in C KFnllDO <- function(Params, do.obs, do.sat, k.gas, z.mix, irr, wtr){ - - # =========================== - # = Unpack and set initials = - # =========================== - #!Pseudocode #1: Initial guesses for B, C, and Q t - c1 <- Params[1] #PAR coeff - c2 <- Params[2] #log(Temp) coeff - Q <- exp(Params[3]) # Variance of the process error - H <- exp(Params[4]) # Variance of observation error - - # See KalmanDO_smooth.R comments for explanation of beta - kz <- k.gas/z.mix # K and Zmix are both vector of length nobs - # beta <- 1-kz # beta is a vector of length nobs (this beta is for difference equation form) - beta <- exp(-kz) # This beta is for using the differential equation form - - # Set first true value equal to first observation - alpha <- do.obs[1]#Let's give this model some starting values - - # Set process covariance, P, equal to Q - P <- Q #starting value - - # Empty vector for nll's - nlls <- rep(0,length(do.obs)) - - # ================== - # = Main Recursion = - # ================== - # =============== - # = Predictions = - # =============== - # Equations for Predictions from Harvey - # a[t|t-1] = T[t]*a[t-1] + c[t] Harvey pg 105 eq. 3.2.2a - # P[t|t-1] = T[t]*P[t-1]*T'[t] + R[t]*Q[t]*R'[t] Harvey pg 106 eq. 3.2.2b - - # Predictions where gas flux not split into beta etc. (I'm pretty sure this is wrong): - # Uk <- K[i-1]*(do.sat[i-1] - alpha)/Zmix[i-1] - # alpha <- alpha + c1*irr[i-1] + c2*log(wtr[i-1]) + Uk - # P <- (Uk*P*Uk) + Q - - # Predictions where gas flux is split into beta: - # Difference Equation Version (original): - # alpha <- beta[i-1]*alpha + c1*irr[i-1] + c2*log(wtr[i-1]) + kz[i-1]*do.sat[i-1] - - # Differential Equation Version: - # Gordon's code (from BayesMetabSS_indivProcessErr_072012.txt): - # alpha[f,j] <- P[f,j] - rho[f,j] + KO2[f,j]* DOSat[f,j]; - # DOHat[f+1,j] <- (alpha[f,j]-(alpha[f,j]-KO2[f,j]*DOTrue[f,j])*exp(-KO2[f,j]))/KO2[f,j]; - - # Define Gordon's "alpha" as "a1": - # a1 = c1*irr + c2*log(wtr) + kz*do.sat # my a1 is his alpha - - # Coefficients in front of alpha are defined above as "beta", and play a particular role in propagating uncertainty - # Need to algebraically regarrange Gordon's equation for DOHat (my "alpha" is his "DOHat") - # Note that alpha and a1 are rewritten each iteration of the loop... - # Therefore, alpha on the left side is "alpha[t]", and alpha on the right side is "alpha[t-1]" - - # alpha = (a1 - (a1 - kz*alpha)*exp(-kz))/kz # my alpha is his DOHat - # alpha = (a1 + (-a1 + kz*alpha)*exp(-kz))/kz # redistribute -1 - # alpha = (a1 + exp(-kz)*(-a1 + kz*alpha))/kz # regarrange - # alpha = (a1 + -exp(-kz)*a1 + exp(-kz)*kz*alpha)/kz # multiply "exp(kz)" through - # alpha = a1/kz + -exp(-kz)*a1/kz + exp(-kz)*alpha # multiply "/kz" through - - # ====================== - # = Updating Equations = - # ====================== - # Updating Equations from Harvey - # a[t] = a[t|t-1] + P[t|t-1]*Z'[t]*F[t]^-1(y[t] - Z[t]*a[t|t-1] - d[t]). Harvey, page 106, 3.2.3a - # P[t] = P[t|t-1] - P[t|t-1]*Z'[t]*F[t]^-1*Z[t]*P[t|t-1] Harvey, page 106, eq. 3.2.3b - # F[t] = Z[t]*P[t|t-1]*Z'[t] + H[t] Harvey, page 106, eq. 3.2.3c - # kalmanLoopC(double *alpha, double *doobs, double *c1, double *c2, double *P, double *Q, double *H, double *beta, double *irr, double *wtr, double *kz, double *dosat, int *nobs) - nlls <- kalmanLoopR(nlls=nlls, alpha=alpha, doobs=do.obs, c1=c1, c2=c2, P=P, Q=Q, H=H, beta=beta, irr=irr, wtr=wtr, kz=kz, dosat=do.sat) - - return(sum(nlls)) # return the sum of nll's + + # =========================== + # = Unpack and set initials = + # =========================== + #!Pseudocode #1: Initial guesses for B, C, and Q t + c1 <- Params[1] #PAR coeff + c2 <- Params[2] #log(Temp) coeff + Q <- exp(Params[3]) # Variance of the process error + H <- exp(Params[4]) # Variance of observation error + + # See KalmanDO_smooth.R comments for explanation of beta + kz <- k.gas/z.mix # K and Zmix are both vector of length nobs + # beta <- 1-kz # beta is a vector of length nobs (this beta is for difference equation form) + beta <- exp(-kz) # This beta is for using the differential equation form + + # Set first true value equal to first observation + alpha <- do.obs[1]#Let's give this model some starting values + + # Set process covariance, P, equal to Q + P <- Q #starting value + + # Empty vector for nll's + nlls <- rep(0,length(do.obs)) + + # ================== + # = Main Recursion = + # ================== + # =============== + # = Predictions = + # =============== + # Equations for Predictions from Harvey + # a[t|t-1] = T[t]*a[t-1] + c[t] Harvey pg 105 eq. 3.2.2a + # P[t|t-1] = T[t]*P[t-1]*T'[t] + R[t]*Q[t]*R'[t] Harvey pg 106 eq. 3.2.2b + + # Predictions where gas flux not split into beta etc. (I'm pretty sure this is wrong): + # Uk <- K[i-1]*(do.sat[i-1] - alpha)/Zmix[i-1] + # alpha <- alpha + c1*irr[i-1] + c2*log(wtr[i-1]) + Uk + # P <- (Uk*P*Uk) + Q + + # Predictions where gas flux is split into beta: + # Difference Equation Version (original): + # alpha <- beta[i-1]*alpha + c1*irr[i-1] + c2*log(wtr[i-1]) + kz[i-1]*do.sat[i-1] + + # Differential Equation Version: + # Gordon's code (from BayesMetabSS_indivProcessErr_072012.txt): + # alpha[f,j] <- P[f,j] - rho[f,j] + KO2[f,j]* DOSat[f,j]; + # DOHat[f+1,j] <- (alpha[f,j]-(alpha[f,j]-KO2[f,j]*DOTrue[f,j])*exp(-KO2[f,j]))/KO2[f,j]; + + # Define Gordon's "alpha" as "a1": + # a1 = c1*irr + c2*log(wtr) + kz*do.sat # my a1 is his alpha + + # Coefficients in front of alpha are defined above as "beta", and play a particular role in propagating uncertainty + # Need to algebraically regarrange Gordon's equation for DOHat (my "alpha" is his "DOHat") + # Note that alpha and a1 are rewritten each iteration of the loop... + # Therefore, alpha on the left side is "alpha[t]", and alpha on the right side is "alpha[t-1]" + + # alpha = (a1 - (a1 - kz*alpha)*exp(-kz))/kz # my alpha is his DOHat + # alpha = (a1 + (-a1 + kz*alpha)*exp(-kz))/kz # redistribute -1 + # alpha = (a1 + exp(-kz)*(-a1 + kz*alpha))/kz # regarrange + # alpha = (a1 + -exp(-kz)*a1 + exp(-kz)*kz*alpha)/kz # multiply "exp(kz)" through + # alpha = a1/kz + -exp(-kz)*a1/kz + exp(-kz)*alpha # multiply "/kz" through + + # ====================== + # = Updating Equations = + # ====================== + # Updating Equations from Harvey + # a[t] = a[t|t-1] + P[t|t-1]*Z'[t]*F[t]^-1(y[t] - Z[t]*a[t|t-1] - d[t]). Harvey, page 106, 3.2.3a + # P[t] = P[t|t-1] - P[t|t-1]*Z'[t]*F[t]^-1*Z[t]*P[t|t-1] Harvey, page 106, eq. 3.2.3b + # F[t] = Z[t]*P[t|t-1]*Z'[t] + H[t] Harvey, page 106, eq. 3.2.3c + # kalmanLoopC(double *alpha, double *doobs, double *c1, double *c2, double *P, double *Q, double *H, double *beta, double *irr, double *wtr, double *kz, double *dosat, int *nobs) + nlls <- kalmanLoopR(nlls=nlls, alpha=alpha, doobs=do.obs, c1=c1, c2=c2, P=P, Q=Q, H=H, beta=beta, irr=irr, wtr=wtr, kz=kz, dosat=do.sat) + + return(sum(nlls)) # return the sum of nll's }#End function @@ -220,127 +260,133 @@ KFnllDO <- function(Params, do.obs, do.sat, k.gas, z.mix, irr, wtr){ # = Kalman Smoother = # =================== KFsmoothDO <- function(Params, do.obs, do.sat, k.gas, z.mix, irr, wtr, Hfac=NULL){ - nobs <- length(do.obs) - d0 <- double(nobs-1) - # beta <- 1-KO2zmix #do.obs_t = 1*do.obs_t-1 + -KO2zmix*do.obs_t-1 + Sea%*%Ewe + eta === (1-KO2zmix)*do.obs_t-1..... - - # Unpack parameters (these were previously fitted) - c1 <- Params[1] # irr Coeff - c2 <- Params[2] # log(wtr) Coeff - Q <- Params[3] # Variance of the Process Error - if(is.null(Hfac)){ - H <- Params[4] - }else{ - H <- Params[4]*Hfac - } - # Variance of Observation Error - - # Need to define portion of K multiplied by state variable (DO) - # Gas flux = K[t-1](do.sat[t-1] - alpha[t-1])/Zmix[t-1] - # Gas flux = (K[t-1]/Zmix[t-1])*(do.sat[t-1]) - (K[t-1]/Zmix[t-1])*(alpha[t-1]) - # Note that K/Z is essentially a coefficient sitting in front of alpha, the estimate of DO - # Therefore, - # alpha[t] = 1*alpha[t-1] + c1*irr[t-1] + c2*log(wtr[t-1]) + K[t-1](do.sat[t-1] - alpha[t-1])/Zmix[t-1] - # Becomes: - # kz[t] = k[t]/Zmix[t] - # alpha[t] = 1*alpha[t-1] + -kz[t-1]*alpha[t-1] + c1*irr[t-1] + c2*log(wtr[t-1]) + kz[t-1]*do.sat[t-1] - # Or, - # alpha[t] = (1-kz[t-1])*alpha[t-1] + c1*irr[t-1] + c2*log(wtr[t-1]) + kz[t-1]*do.sat[t-1] - # Defining kz and redefining (1-kz[t]) as beta[t]: - kz <- k.gas/z.mix # K and Zmix are both vector of length nobs - # beta <- 1-kz # beta is a vector of length nobs (this beta is for difference equation form) - beta <- exp(-kz) # This beta is for using the differential equation form - - # Set first true value equal to first observation - alpha <- do.obs[1] - - # Set process covariance, P, equal to Q - P <- Q # starting value - - # Initial values - aHat <- c(alpha, d0) # aHat[t] == "a[t|t-1]" (estimate of a before updating) - pHat <- c(P, d0) # pHat[t] == "p[t|t-1]" (estimate of a before updating) - aVec <- aHat # aVec[t] == "a[t|t]" or "a[t]" (aVec is the "updated" version of aHat) - pVec <- pHat # pVec[t] == "P[t|t]" or "P[t]" (pVec is the "updated" version of pHat) - etaVec <- double(nobs) - - for(i in 2:nobs){ - # =============== - # = Predictions = - # =============== - # Equations for Predictions from Harvey - # a[t|t-1] = T[t]*a[t-1] + c[t] Harvey pg 105 eq. 3.2.2a - # P[t|t-1] = T[t]*P[t-1]*T'[t] + R[t]*Q[t]*R'[t] Harvey pg 106 eq. 3.2.2b - - # Predictions where gas flux not split into beta etc.: - # Uk <- K[i-1]*(do.sat[i-1] - alpha)/Zmix[i-1] - # alpha <- alpha + c1*irr[i-1] + c2*log(wtr[i-1]) + Uk - # aHat[i] <- alpha - # P <- (Uk*P*Uk) + Q - # pHat[i] <- P - - # Predictions where gas flux is split into beta (see explanation above): - - # Difference Equation Version: - # alpha <- beta[i-1]*alpha + c1*irr[i-1] + c2*log(wtr[i-1]) + kz[i-1]*do.sat[i-1] - - # Differential Equation Version (see kalmanDO_nll.R for explanation): - if(is.finite(1/kz[i-1])){ - - a1 <- c1*irr[i-1] + c2*log(wtr[i-1]) + kz[i-1]*do.sat[i-1] - alpha <- a1/kz[i-1] + -beta[i-1]*a1/kz[i-1] + beta[i-1]*alpha # NOTE: beta==exp(-kz); kz=K/Zmix - - }else{ - - alpha <- c1*irr[i-1] + c2*log(wtr[i-1]) - - } - - - aHat[i] <- alpha - P <- (beta[i-1]*P*beta[i-1]) + Q - pHat[i] <- P - - # ====================== - # = Updating Equations = - # ====================== - # Updating Equations from Harvey - # a[t] = a[t|t-1] + P[t|t-1]*Z'[t]*F[t]^-1(y[t] - Z[t]*a[t|t-1] - d[t]). Harvey, page 106, 3.2.3a - # P[t] = P[t|t-1] - P[t|t-1]*Z'[t]*F[t]^-1*Z[t]*P[t|t-1] Harvey, page 106, eq. 3.2.3b - # F[t] = Z[t]*P[t|t-1]*Z'[t] + H[t] Harvey, page 106, eq. 3.2.3c - - eta <- do.obs[i] - alpha - Eff <- P + H - alpha <- alpha + P/Eff*eta - P <- P - P*P/Eff - - aVec[i] <- alpha - pVec[i] <- P - etaVec[i] <- eta - } - - #Kalman Smoother - aSmooth <- rep(NA,nobs) - Psmooth <- rep(NA,nobs) - aSmooth[nobs] <- aVec[nobs] # "starting" value for smoother (smoother starts at end and works backwards) - # pSmooth[nobs] <- pVec[nobs] - - # Filtering is informed by past information - # Smoothing includes the information from filtering (estimates of parameters), but also future information. - # "The aim of filtering is to find the expected value of the state vector, alpha[t], conditional on the information available at time t, that is E(alpha[t]|Y[t]). The aim of smoothing is to take account of the information made available after time t. The mean of the distribution of alpha[t], conditional on all the sample, may be written as E(alpha[t]|Y[T]) and is known as the smoothed estimate. THe corresponding estimator is called the SMOOTHER. Since the smoother is based on more information than the filtered estimator, it will have a MSE which, in general, is smaller than that of the filtered estimator; it cannot be greater." ~ Harvey 1989, pgs 149-150. - #a[t|T] = a[t] + Pstar[t]*(a[t+1|T] - T[t+1]*a[t]) - # P[t|T] = P[t] + Pstar[t]*(P[t+1|T] - P[t+1|t])*Pstar[t] - # Pstar[t] = P[t]*T[t+1]/P[t+1|t] - # t is current time step, T is last time step (when in []), T contains AR parameters (when NOT in []) - - for(i in length(d0):1){ - pStar <- pVec[i]*beta[i+1]/pHat[i+1] - aSmooth[i] <- aVec[i] + pStar*(aSmooth[i+1] - aHat[i+1]) - - # CAN ALSO SMOOTH P, WHICH GIVES THE SMOOTHED COVARIANCE MATRIX (not a matrix for univariate; gives estimate of accuracy of state estimate) - # pSmooth[i] <- pVec[i] + pStar*(pSmooth[i+1] - pHat[i+1])*pStar - } - return(aSmooth) # return smoothed DO time series + nobs <- length(do.obs) + d0 <- double(nobs-1) + # beta <- 1-KO2zmix #do.obs_t = 1*do.obs_t-1 + -KO2zmix*do.obs_t-1 + Sea%*%Ewe + eta === (1-KO2zmix)*do.obs_t-1..... + + # Unpack parameters (these were previously fitted) + c1 <- Params[1] # irr Coeff + c2 <- Params[2] # log(wtr) Coeff + Q <- Params[3] # Variance of the Process Error + if(is.null(Hfac)){ + H <- Params[4] + }else{ + H <- Params[4]*Hfac + } + # Variance of Observation Error + + # Need to define portion of K multiplied by state variable (DO) + # Gas flux = K[t-1](do.sat[t-1] - alpha[t-1])/Zmix[t-1] + # Gas flux = (K[t-1]/Zmix[t-1])*(do.sat[t-1]) - (K[t-1]/Zmix[t-1])*(alpha[t-1]) + # Note that K/Z is essentially a coefficient sitting in front of alpha, the estimate of DO + # Therefore, + # alpha[t] = 1*alpha[t-1] + c1*irr[t-1] + c2*log(wtr[t-1]) + K[t-1](do.sat[t-1] - alpha[t-1])/Zmix[t-1] + # Becomes: + # kz[t] = k[t]/Zmix[t] + # alpha[t] = 1*alpha[t-1] + -kz[t-1]*alpha[t-1] + c1*irr[t-1] + c2*log(wtr[t-1]) + kz[t-1]*do.sat[t-1] + # Or, + # alpha[t] = (1-kz[t-1])*alpha[t-1] + c1*irr[t-1] + c2*log(wtr[t-1]) + kz[t-1]*do.sat[t-1] + # Defining kz and redefining (1-kz[t]) as beta[t]: + kz <- k.gas/z.mix # K and Zmix are both vector of length nobs + # beta <- 1-kz # beta is a vector of length nobs (this beta is for difference equation form) + beta <- exp(-kz) # This beta is for using the differential equation form + + # Set first true value equal to first observation + alpha <- do.obs[1] + + # Set process covariance, P, equal to Q + P <- Q # starting value + + # Initial values + aHat <- c(alpha, d0) # aHat[t] == "a[t|t-1]" (estimate of a before updating) + pHat <- c(P, d0) # pHat[t] == "p[t|t-1]" (estimate of a before updating) + aVec <- aHat # aVec[t] == "a[t|t]" or "a[t]" (aVec is the "updated" version of aHat) + pVec <- pHat # pVec[t] == "P[t|t]" or "P[t]" (pVec is the "updated" version of pHat) + etaVec <- double(nobs) + + for(i in 2:nobs){ + # =============== + # = Predictions = + # =============== + # Equations for Predictions from Harvey + # a[t|t-1] = T[t]*a[t-1] + c[t] Harvey pg 105 eq. 3.2.2a + # P[t|t-1] = T[t]*P[t-1]*T'[t] + R[t]*Q[t]*R'[t] Harvey pg 106 eq. 3.2.2b + + # Predictions where gas flux not split into beta etc.: + # Uk <- K[i-1]*(do.sat[i-1] - alpha)/Zmix[i-1] + # alpha <- alpha + c1*irr[i-1] + c2*log(wtr[i-1]) + Uk + # aHat[i] <- alpha + # P <- (Uk*P*Uk) + Q + # pHat[i] <- P + + # Predictions where gas flux is split into beta (see explanation above): + + # Difference Equation Version: + # alpha <- beta[i-1]*alpha + c1*irr[i-1] + c2*log(wtr[i-1]) + kz[i-1]*do.sat[i-1] + + # Differential Equation Version (see kalmanDO_nll.R for explanation): + if(is.finite(1/kz[i-1])){ + + a1 <- c1*irr[i-1] + c2*log(wtr[i-1]) + kz[i-1]*do.sat[i-1] + alpha <- a1/kz[i-1] + -beta[i-1]*a1/kz[i-1] + beta[i-1]*alpha # NOTE: beta==exp(-kz); kz=K/Zmix + + }else{ + + alpha <- c1*irr[i-1] + c2*log(wtr[i-1]) + + } + + + aHat[i] <- alpha + P <- (beta[i-1]*P*beta[i-1]) + Q + pHat[i] <- P + + # ====================== + # = Updating Equations = + # ====================== + # Updating Equations from Harvey + # a[t] = a[t|t-1] + P[t|t-1]*Z'[t]*F[t]^-1(y[t] - Z[t]*a[t|t-1] - d[t]). Harvey, page 106, 3.2.3a + # P[t] = P[t|t-1] - P[t|t-1]*Z'[t]*F[t]^-1*Z[t]*P[t|t-1] Harvey, page 106, eq. 3.2.3b + # F[t] = Z[t]*P[t|t-1]*Z'[t] + H[t] Harvey, page 106, eq. 3.2.3c + + eta <- do.obs[i] - alpha + Eff <- P + H + alpha <- alpha + P/Eff*eta + P <- P - P*P/Eff + + aVec[i] <- alpha + pVec[i] <- P + etaVec[i] <- eta + } + + #Kalman Smoother + aSmooth <- rep(NA,nobs) + Psmooth <- rep(NA,nobs) + aSmooth[nobs] <- aVec[nobs] # "starting" value for smoother (smoother starts at end and works backwards) + # pSmooth[nobs] <- pVec[nobs] + + # Filtering is informed by past information + # Smoothing includes the information from filtering (estimates of parameters), but also future information. + # "The aim of filtering is to find the expected value of the state vector, alpha[t], conditional on the information available at time t, that is E(alpha[t]|Y[t]). The aim of smoothing is to take account of the information made available after time t. The mean of the distribution of alpha[t], conditional on all the sample, may be written as E(alpha[t]|Y[T]) and is known as the smoothed estimate. THe corresponding estimator is called the SMOOTHER. Since the smoother is based on more information than the filtered estimator, it will have a MSE which, in general, is smaller than that of the filtered estimator; it cannot be greater." ~ Harvey 1989, pgs 149-150. + #a[t|T] = a[t] + Pstar[t]*(a[t+1|T] - T[t+1]*a[t]) + # P[t|T] = P[t] + Pstar[t]*(P[t+1|T] - P[t+1|t])*Pstar[t] + # Pstar[t] = P[t]*T[t+1]/P[t+1|t] + # t is current time step, T is last time step (when in []), T contains AR parameters (when NOT in []) + + for(i in length(d0):1){ + pStar <- pVec[i]*beta[i+1]/pHat[i+1] + aSmooth[i] <- aVec[i] + pStar*(aSmooth[i+1] - aHat[i+1]) + + # CAN ALSO SMOOTH P, WHICH GIVES THE SMOOTHED COVARIANCE MATRIX (not a matrix for univariate; gives estimate of accuracy of state estimate) + # pSmooth[i] <- pVec[i] + pStar*(pSmooth[i+1] - pHat[i+1])*pStar + } + return(list(smoothDO = aSmooth, # return smoothed DO time series + filteredDO = aVec, + predictedDO = aHat, + innovations = etaVec, + innovation.var = pHat + H, + innovation.std = etaVec / sqrt(pHat + H), + state.var = pVec)) } # =========================================== @@ -348,8 +394,97 @@ KFsmoothDO <- function(Params, do.obs, do.sat, k.gas, z.mix, irr, wtr, Hfac=NULL # =========================================== # kalmanLoopC(double *alpha, double *doobs, double *c1, double *c2, double *P, double *Q, double *H, double *beta, double *irr, double *wtr, double *kz, double *dosat, int *nobs) kalmanLoopR <- function(nlls, alpha, doobs, c1, c2, P, Q, H, beta, irr, wtr, kz, dosat){ - nobs <- length(doobs) - a.loop <- .C("kalmanLoopC", nlls=as.double(nlls), as.double(alpha), as.double(doobs), as.double(c1), as.double(c2), as.double(P), as.double(Q), as.double(H), as.double(beta), as.double(irr), as.double(wtr), as.double(kz), as.double(dosat), as.integer(nobs), PACKAGE="LakeMetabolizer") - return(a.loop[["nlls"]]) + nobs <- length(doobs) + a.loop <- .C("kalmanLoopC", nlls=as.double(nlls), as.double(alpha), as.double(doobs), as.double(c1), as.double(c2), as.double(P), as.double(Q), as.double(H), as.double(beta), as.double(irr), as.double(wtr), as.double(kz), as.double(dosat), as.integer(nobs), PACKAGE="LakeMetabolizer") + return(a.loop[["nlls"]]) } +# ========================================= +# = Function to bootstrap for uncertainty = +# ========================================= +bootstrap.kalman <- function(n.boot, Params, KFresults, guesses, do.obs, do.sat, k.gas, freq, z.mix, irr, wtr, Hfac = NULL){ + + nobs <- length(do.obs) + #d0 <- double(nobs-1) + # beta <- 1-KO2zmix #do.obs_t = 1*do.obs_t-1 + -KO2zmix*do.obs_t-1 + Sea%*%Ewe + eta === (1-KO2zmix)*do.obs_t-1..... + + # Unpack parameters (these were previously fitted) + c1 <- Params[1] # irr Coeff + c2 <- Params[2] # log(wtr) Coeff + Q <- Params[3] # Variance of the Process Error + if(is.null(Hfac)){ + H <- Params[4] + }else{ + H <- Params[4]*Hfac + } + # Variance of Observation Error + + kz <- k.gas / z.mix + beta <- exp(-kz) + + ## Initialize results object + boot.results <- data.frame(boot.iter = seq_len(n.boot), + gppCoeff = rep(NA, n.boot), + rCoeff = rep(NA, n.boot), + Q = rep(NA, n.boot), + H = rep(NA, n.boot), + convergence = rep(NA, n.boot), + nll = rep(NA, n.boot), + GPP = rep(NA, n.boot), + R = rep(NA, n.boot), + NEP = rep(NA, n.boot)) + + # KFresults for reference + # list(smoothDO = aSmooth, # return smoothed DO time series + # filteredDO = aVec, + # predictedDO = aHat, + # innovations = etaVec, + # innovation.var = pHat + H, + # innovation.std = etaVec / sqrt(pHat + H), + # state.var = pVec) + + for (j in boot.results$boot.iter){ + + # Resample innovations + innovation.star <- sample(KFresults$innovation.std, replace = TRUE) + eta.star <- innovation.star * sqrt(KFresults$innovation.var) + + do.sim <- numeric(nobs) + do.sim[1] <- do.obs[1] + alpha <- do.obs[1] + + # Generate Simulated DO time series from innovations + for(i in 2:nobs){ + if(is.finite(1/kz[i-1])){ + a1 <- c1*irr[i-1] + c2*log(wtr[i-1]) + kz[i-1]*do.sat[i-1] + alpha <- a1/kz[i-1] + -beta[i-1]*a1/kz[i-1] + beta[i-1]*alpha # NOTE: beta==exp(-kz); kz=K/Zmix + }else{ + alpha <- c1*irr[i-1] + c2*log(wtr[i-1]) + } + do.sim[i] <- alpha + eta.star[i] + } + + # MLE + simFit <- optim(guesses, fn=KFnllDO, do.obs=do.sim, do.sat=do.sat, k.gas=k.gas, z.mix=z.mix, irr=irr, wtr=wtr) + simFitPar <- simFit$par + simFitGPP <- mean(simFitPar[1] * irr, na.rm = TRUE) * freq + simFitR <- mean(simFitPar[2] * log(wtr), na.rm = TRUE) * freq + + # Writing results to table + boot.results[j, "gppCoeff"] <- simFitPar[1] + boot.results[j, "rCoeff"] <- simFitPar[2] + boot.results[j, "Q"] <- exp(simFitPar[3]) + boot.results[j, "H"] <- exp(simFitPar[4]) + boot.results[j, "convergence"] <- simFit$convergence + boot.results[j, "nll"] <- simFit$value + boot.results[j, "GPP"] <- simFitGPP + boot.results[j, "R"] <- simFitR + boot.results[j, "NEP"] <- simFitGPP + simFitR + + } + + return(boot.results) + +} + + diff --git a/man/metab.kalman.Rd b/man/metab.kalman.Rd index b28d3fb..5dd4963 100644 --- a/man/metab.kalman.Rd +++ b/man/metab.kalman.Rd @@ -4,7 +4,7 @@ \alias{metab.kalman} \title{Metabolism calculated from parameters estimated using a Kalman filter} \usage{ -metab.kalman(do.obs, do.sat, k.gas, z.mix, irr, wtr, ...) +metab.kalman(do.obs, do.sat, k.gas, z.mix, irr, wtr, n.boot = 0, ...) } \arguments{ \item{do.obs}{Vector of dissovled oxygen concentration observations, \eqn{mg O[2] L^{-1}}{mg O2 / L}} @@ -20,6 +20,8 @@ metab.kalman(do.obs, do.sat, k.gas, z.mix, irr, wtr, ...) \item{wtr}{Vector of water temperatures in \eqn{^{\circ}C}{degrees C}. Used in scaling respiration with temperature} +\item{n.boot}{Numeric of the number of bootstrap iterations. Set to 0 to bypass uncertainty quantification (default). Set to 2 or greater to quantify uncertainty via bootstrapping.} + \item{...}{additional arguments; currently "datetime" is the only recognized argument passed through \code{...}} } \value{ @@ -29,6 +31,18 @@ A data.frame with columns corresponding to components of metabolism \item{R}{numeric estimate of Respiration, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} \item{NEP}{numeric estimate of Net Ecosystem production, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} } +If n.boot is 2 or greater, the data frame will also contain the following columns: +\describe{ +\item{GPP.lci}{lower 95\% confidence interval of the estimate of Gross Primary Production, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} +\item{GPP.uci}{upper 95\% confidence interval of the estimate of Gross Primary Production, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} +\item{GPP.cv}{coefficient of variation of the bootstrapped estimates of Gross Primary Production} +\item{R.lci}{lower 95\% confidence interval of the estimate of Respiration, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} +\item{R.uci}{upper 95\% confidence interval of the estimate of Respiration, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} +\item{R.cv}{coefficient of variation of the bootstrapped estimates of Respiration} +\item{NEP.lci}{lower 95\% confidence interval of the estimate of Net Ecosystem Production, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} +\item{NEP.uci}{upper 95\% confidence interval of the estimate of Net Ecosystem Production, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} +\item{NEP.cv}{coefficient of variation of the bootstrapped estimates of Net Ecosystem Production} +} Use \link{attributes} to access more model output: @@ -99,7 +113,10 @@ metab.kalman(irr=irr[,2], z.mix=rep(1, length(k.gas)), } \references{ Batt, Ryan D. and Stephen R. Carpenter. 2012. \emph{Free-water lake metabolism: -addressing noisy time series with a Kalman filter}. Limnology and Oceanography: Methods 10: 20-30. doi: 10.4319/lom.2012.10.20 +addressing noisy time series with a Kalman filter}. Limnology and Oceanography: Methods 10: 20-30. \url{https://doi.org/10.4319/lom.2012.10.20} + +Stoffer, David S., and Kent D. Wall. 1991. \emph{Bootstrapping State-Space Models: Gaussian Maximum Likelihood Estimation and the Kalman Filter}. +Journal of the American Statistical Association 86 (416): 1024–33. \url{https://doi.org/10.2307/2290521}. } \seealso{ \link{temp.kalman}, \link{watts.in}, \link{metab}, \link{metab.bookkeep}, \link{metab.ols}, \link{metab.mle}, \link{metab.bayesian} From 1c19c618df242e522e17e14ac630b912bd6d2183 Mon Sep 17 00:00:00 2001 From: Bennett McAfee <105509324+bmcafee@users.noreply.github.com> Date: Thu, 2 Jul 2026 13:45:53 -0400 Subject: [PATCH 03/13] Constrain GPP and R to correct sign in metab.kalman --- R/metab.kalman.R | 16 ++++++++-------- man/metab.kalman.Rd | 2 +- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/R/metab.kalman.R b/R/metab.kalman.R index 2565a12..3ce46e9 100644 --- a/R/metab.kalman.R +++ b/R/metab.kalman.R @@ -131,10 +131,10 @@ metab.kalman <- function(do.obs, do.sat, k.gas, z.mix, irr, wtr, n.boot = 0, ... } # Filter and fit - guesses <- c(1E-4,1E-4,log(5),log(5)) + guesses <- c(log(1E-4),log(1E-4),log(5),log(5)) fit <- optim(guesses, fn=KFnllDO, do.obs=do.obs, do.sat=do.sat, k.gas=(k.gas/freq), z.mix=z.mix, irr=irr, wtr=wtr) pars0 <- fit$par - pars <- c("gppCoeff"=pars0[1], "rCoeff"=pars0[2], "Q"=exp(pars0[3]), "H"=exp(pars0[4])) + pars <- c("gppCoeff"=exp(pars0[1]), "rCoeff"=-exp(pars0[2]), "Q"=exp(pars0[3]), "H"=exp(pars0[4])) # Smooth KFresults <- KFsmoothDO(pars, do.obs=do.obs, do.sat=do.sat, k.gas=(k.gas/freq), z.mix=z.mix, irr=irr, wtr=wtr) @@ -184,8 +184,8 @@ KFnllDO <- function(Params, do.obs, do.sat, k.gas, z.mix, irr, wtr){ # = Unpack and set initials = # =========================== #!Pseudocode #1: Initial guesses for B, C, and Q t - c1 <- Params[1] #PAR coeff - c2 <- Params[2] #log(Temp) coeff + c1 <- exp(Params[1]) #PAR coeff + c2 <- -exp(Params[2]) #log(Temp) coeff Q <- exp(Params[3]) # Variance of the process error H <- exp(Params[4]) # Variance of observation error @@ -467,12 +467,12 @@ bootstrap.kalman <- function(n.boot, Params, KFresults, guesses, do.obs, do.sat, # MLE simFit <- optim(guesses, fn=KFnllDO, do.obs=do.sim, do.sat=do.sat, k.gas=k.gas, z.mix=z.mix, irr=irr, wtr=wtr) simFitPar <- simFit$par - simFitGPP <- mean(simFitPar[1] * irr, na.rm = TRUE) * freq - simFitR <- mean(simFitPar[2] * log(wtr), na.rm = TRUE) * freq + simFitGPP <- mean(exp(simFitPar[1]) * irr, na.rm = TRUE) * freq + simFitR <- mean(-exp(simFitPar[2]) * log(wtr), na.rm = TRUE) * freq # Writing results to table - boot.results[j, "gppCoeff"] <- simFitPar[1] - boot.results[j, "rCoeff"] <- simFitPar[2] + boot.results[j, "gppCoeff"] <- exp(simFitPar[1]) + boot.results[j, "rCoeff"] <- -exp(simFitPar[2]) boot.results[j, "Q"] <- exp(simFitPar[3]) boot.results[j, "H"] <- exp(simFitPar[4]) boot.results[j, "convergence"] <- simFit$convergence diff --git a/man/metab.kalman.Rd b/man/metab.kalman.Rd index 5dd4963..c2ba52c 100644 --- a/man/metab.kalman.Rd +++ b/man/metab.kalman.Rd @@ -122,5 +122,5 @@ Journal of the American Statistical Association 86 (416): 1024–33. \url{https: \link{temp.kalman}, \link{watts.in}, \link{metab}, \link{metab.bookkeep}, \link{metab.ols}, \link{metab.mle}, \link{metab.bayesian} } \author{ -Ryan Batt, Luke A. Winslow +Ryan Batt, Luke A. Winslow, Bennett McAfee } From efc6ef4bbca0b743392d58325149c7c5662483f7 Mon Sep 17 00:00:00 2001 From: Bennett McAfee <105509324+bmcafee@users.noreply.github.com> Date: Wed, 8 Jul 2026 16:12:35 -0400 Subject: [PATCH 04/13] Add constrain.sign argument to mle and kalman --- R/metab.kalman.R | 50 ++++++++++++++++++++++--------- R/metab.mle.R | 73 ++++++++++++++++++++++++++++++++------------- man/metab.kalman.Rd | 14 ++++++++- man/metab.mle.Rd | 3 ++ 4 files changed, 105 insertions(+), 35 deletions(-) diff --git a/R/metab.kalman.R b/R/metab.kalman.R index 3ce46e9..322afe7 100644 --- a/R/metab.kalman.R +++ b/R/metab.kalman.R @@ -9,6 +9,7 @@ #'@param irr Vector of photosynthetically active radiation in \eqn{\mu mol\ m^{-2} s^{-1}}{micro mols / m^2 / s} #'@param wtr Vector of water temperatures in \eqn{^{\circ}C}{degrees C}. Used in scaling respiration with temperature #'@param n.boot Numeric of the number of bootstrap iterations. Set to 0 to bypass uncertainty quantification (default). Set to 2 or greater to quantify uncertainty via bootstrapping. +#'@param constrain.sign Logical. If TRUE, estimates of GPP and R coefficients with be constrained to positive and negative values, respectively. #'@param ... additional arguments; currently "datetime" is the only recognized argument passed through \code{...} #'@return #'A data.frame with columns corresponding to components of metabolism @@ -100,7 +101,7 @@ #' do.sat=do.sat, wtr=wtr[,2], #' k.gas=k.gas, do.obs=doobs[,2]) #'@export -metab.kalman <- function(do.obs, do.sat, k.gas, z.mix, irr, wtr, n.boot = 0, ...){ +metab.kalman <- function(do.obs, do.sat, k.gas, z.mix, irr, wtr, n.boot = 0, constrain.sign = TRUE, ...){ complete.inputs(do.obs=do.obs, do.sat=do.sat, k.gas=k.gas, z.mix=z.mix, irr=irr, wtr=wtr, error=TRUE) @@ -131,10 +132,19 @@ metab.kalman <- function(do.obs, do.sat, k.gas, z.mix, irr, wtr, n.boot = 0, ... } # Filter and fit - guesses <- c(log(1E-4),log(1E-4),log(5),log(5)) - fit <- optim(guesses, fn=KFnllDO, do.obs=do.obs, do.sat=do.sat, k.gas=(k.gas/freq), z.mix=z.mix, irr=irr, wtr=wtr) + if (constrain.sign == TRUE){ + guesses <- c(log(1E-4),log(1E-4),log(5),log(5)) + } else { + guesses <- c(1E-4,1E-4,log(5),log(5)) + } + fit <- optim(guesses, fn=KFnllDO, do.obs=do.obs, do.sat=do.sat, k.gas=(k.gas/freq), z.mix=z.mix, irr=irr, wtr=wtr, constrain.sign=constrain.sign) pars0 <- fit$par - pars <- c("gppCoeff"=exp(pars0[1]), "rCoeff"=-exp(pars0[2]), "Q"=exp(pars0[3]), "H"=exp(pars0[4])) + if (constrain.sign == TRUE){ + pars <- c("gppCoeff"=exp(pars0[1]), "rCoeff"=-exp(pars0[2]), "Q"=exp(pars0[3]), "H"=exp(pars0[4])) + } else { + pars <- c("gppCoeff"=pars0[1], "rCoeff"=pars0[2], "Q"=exp(pars0[3]), "H"=exp(pars0[4])) + } + # Smooth KFresults <- KFsmoothDO(pars, do.obs=do.obs, do.sat=do.sat, k.gas=(k.gas/freq), z.mix=z.mix, irr=irr, wtr=wtr) @@ -150,7 +160,7 @@ metab.kalman <- function(do.obs, do.sat, k.gas, z.mix, irr, wtr, n.boot = 0, ... } else if (n.boot >= 2){ - boot.results <- bootstrap.kalman(n.boot, Params = pars, KFresults, guesses, do.obs, do.sat, k.gas =(k.gas/freq), freq, z.mix, irr, wtr) + boot.results <- bootstrap.kalman(n.boot, Params = pars, KFresults, guesses, do.obs, do.sat, k.gas =(k.gas/freq), freq, z.mix, irr, wtr, constrain.sign) ci.GPP <- stats::quantile(boot.results$GPP, c(0.025, 0.975), na.rm=TRUE) # 95% confidence interval ci.R <- stats::quantile(boot.results$R, c(0.025, 0.975), na.rm=TRUE) # 95% confidence interval ci.NEP <- stats::quantile(boot.results$NEP, c(0.025, 0.975), na.rm=TRUE) # 95% confidence interval @@ -178,14 +188,19 @@ metab.kalman <- function(do.obs, do.sat, k.gas, z.mix, irr, wtr, n.boot = 0, ... # = Kalman filter/ nll = # ====================== # Main recursion written in C -KFnllDO <- function(Params, do.obs, do.sat, k.gas, z.mix, irr, wtr){ +KFnllDO <- function(Params, do.obs, do.sat, k.gas, z.mix, irr, wtr, constrain.sign){ # =========================== # = Unpack and set initials = # =========================== #!Pseudocode #1: Initial guesses for B, C, and Q t - c1 <- exp(Params[1]) #PAR coeff - c2 <- -exp(Params[2]) #log(Temp) coeff + if (constrain.sign == TRUE){ + c1 <- exp(Params[1]) #PAR coeff + c2 <- -exp(Params[2]) #log(Temp) coeff + } else { + c1 <- Params[1] #PAR coeff + c2 <- Params[2] #log(Temp) coeff + } Q <- exp(Params[3]) # Variance of the process error H <- exp(Params[4]) # Variance of observation error @@ -402,7 +417,7 @@ kalmanLoopR <- function(nlls, alpha, doobs, c1, c2, P, Q, H, beta, irr, wtr, kz, # ========================================= # = Function to bootstrap for uncertainty = # ========================================= -bootstrap.kalman <- function(n.boot, Params, KFresults, guesses, do.obs, do.sat, k.gas, freq, z.mix, irr, wtr, Hfac = NULL){ +bootstrap.kalman <- function(n.boot, Params, KFresults, guesses, do.obs, do.sat, k.gas, freq, z.mix, irr, wtr, constrain.sign, Hfac = NULL){ nobs <- length(do.obs) #d0 <- double(nobs-1) @@ -465,14 +480,21 @@ bootstrap.kalman <- function(n.boot, Params, KFresults, guesses, do.obs, do.sat, } # MLE - simFit <- optim(guesses, fn=KFnllDO, do.obs=do.sim, do.sat=do.sat, k.gas=k.gas, z.mix=z.mix, irr=irr, wtr=wtr) + simFit <- optim(guesses, fn=KFnllDO, do.obs=do.sim, do.sat=do.sat, k.gas=k.gas, z.mix=z.mix, irr=irr, wtr=wtr, constrain.sign=constrain.sign) simFitPar <- simFit$par - simFitGPP <- mean(exp(simFitPar[1]) * irr, na.rm = TRUE) * freq - simFitR <- mean(-exp(simFitPar[2]) * log(wtr), na.rm = TRUE) * freq + if (constrain.sign == TRUE){ + simFit.gppCoeff <- exp(simFitPar[1]) + simFit.rCoeff <- -exp(simFitPar[2]) + } else { + simFit.gppCoeff <- simFitPar[1] + simFit.rCoeff <- simFitPar[2] + } + simFitGPP <- mean(simFit.gppCoeff * irr, na.rm = TRUE) * freq + simFitR <- mean(simFit.rCoeff * log(wtr), na.rm = TRUE) * freq # Writing results to table - boot.results[j, "gppCoeff"] <- exp(simFitPar[1]) - boot.results[j, "rCoeff"] <- -exp(simFitPar[2]) + boot.results[j, "gppCoeff"] <- simFit.gppCoeff + boot.results[j, "rCoeff"] <- simFit.rCoeff boot.results[j, "Q"] <- exp(simFitPar[3]) boot.results[j, "H"] <- exp(simFitPar[4]) boot.results[j, "convergence"] <- simFit$convergence diff --git a/R/metab.mle.R b/R/metab.mle.R index 00c1632..e1dc6c7 100644 --- a/R/metab.mle.R +++ b/R/metab.mle.R @@ -10,6 +10,7 @@ #'@param error.type Option specifying if model should assume pure Process Error 'PE' or Observation Error 'OE'. Defaults to observation error 'OE'. #'@param n.boot Numeric of the number of bootstrap iterations. Set to 0 to bypass uncertainty quantification (default). Set to 2 or greater to quantify uncertainty via bootstrapping. #'@param ar1.resids Logical. If TRUE, the AR(1) structure of the residuals will be retained during randomization when bootstrapping. +#'@param constrain.sign Logical. If TRUE, estimates of GPP and R coefficients with be constrained to positive and negative values, respectively. #'@param ... additional arguments; currently "datetime" is the only recognized argument passed through \code{...} #'@return #'A data.frame with columns corresponding to components of metabolism @@ -109,7 +110,7 @@ #' #'metab.mle(doobs[,2], do.sat, k.gas, z.mix[,2], irr[,2], wtr[,3]) #'@export -metab.mle <- function(do.obs, do.sat, k.gas, z.mix, irr, wtr, error.type="OE", n.boot = 0, ar1.resids = TRUE, ...){ +metab.mle <- function(do.obs, do.sat, k.gas, z.mix, irr, wtr, error.type="OE", n.boot = 0, ar1.resids = TRUE, constrain.sign = TRUE, ...){ complete.inputs(do.obs=do.obs, do.sat=do.sat, k.gas=k.gas, z.mix=z.mix, irr=irr, wtr=wtr, error=TRUE) @@ -152,24 +153,38 @@ metab.mle <- function(do.obs, do.sat, k.gas, z.mix, irr, wtr, error.type="OE", n Q0 <- ((diff(range(do.obs,na.rm=TRUE)) - mean(do.obs,na.rm=TRUE))^2 / length(do.obs)) - guesses <- c(log(1E-4), log(1E-4), log(Q0)) + if (constrain.sign == TRUE){ + guesses <- c(log(1E-4), log(1E-4), log(Q0)) + } else { + guesses <- c(1E-4, 1E-4, log(Q0)) + } + #We have a different number of fitted parameters depending on error type of the model if(error.type=='OE'){ guesses <- c(guesses, do.obs[1]) - fit <- optim(guesses, fn=mleNllOE, do.obs=do.obs, do.sat=do.sat, k.gas=(k.gas/freq), z.mix=z.mix, irr=irr, wtr=wtr) + fit <- optim(guesses, fn=mleNllOE, do.obs=do.obs, do.sat=do.sat, k.gas=(k.gas/freq), z.mix=z.mix, irr=irr, wtr=wtr, constrain.sign = constrain.sign) pars0 <- fit$par - pars <- c("gppCoeff"=exp(pars0[1]), "rCoeff"=-exp(pars0[2]), "Q"=exp(pars0[3]), "nll"=fit$value, "doInit"=pars0[4]) + if (constrain.sign == TRUE){ + pars <- c("gppCoeff"=exp(pars0[1]), "rCoeff"=-exp(pars0[2]), "Q"=exp(pars0[3]), "nll"=fit$value, "doInit"=pars0[4]) + } else { + pars <- c("gppCoeff"=pars0[1], "rCoeff"=pars0[2], "Q"=exp(pars0[3]), "nll"=fit$value, "doInit"=pars0[4]) + } + }else if(error.type=='PE'){ guesses <- c(guesses) - fit <- optim(guesses, fn=mleNllPE, do.obs=do.obs, do.sat=do.sat, k.gas=(k.gas/freq), z.mix=z.mix, irr=irr, wtr=wtr) + fit <- optim(guesses, fn=mleNllPE, do.obs=do.obs, do.sat=do.sat, k.gas=(k.gas/freq), z.mix=z.mix, irr=irr, wtr=wtr, constrain.sign = constrain.sign) pars0 <- fit$par - pars <- c("gppCoeff"=exp(pars0[1]), "rCoeff"=-exp(pars0[2]), "Q"=exp(pars0[3]), "nll"=fit$value) + if (constrain.sign == TRUE){ + pars <- c("gppCoeff"=exp(pars0[1]), "rCoeff"=-exp(pars0[2]), "Q"=exp(pars0[3]), "nll"=fit$value, "doInit"=pars0[4]) + } else { + pars <- c("gppCoeff"=pars0[1], "rCoeff"=pars0[2], "Q"=exp(pars0[3]), "nll"=fit$value, "doInit"=pars0[4]) + } }else{ stop("error.type must be either 'OE' or 'PE', Observation Error or Process Error respectively.") @@ -187,7 +202,7 @@ metab.mle <- function(do.obs, do.sat, k.gas, z.mix, irr, wtr, error.type="OE", n if (n.boot == 0){ results <- list("params"=pars, "metab"=c("GPP"=GPP,"R"=R,"NEP"=GPP+R)) } else if (n.boot >= 2){ - boot.results <- bootstrap.mle(n.boot, do.obs, pars, guesses, do.sat, k.gas, freq, z.mix, irr, wtr, error.type, ar1.resids) + boot.results <- bootstrap.mle(n.boot, do.obs, pars, guesses, do.sat, k.gas, freq, z.mix, irr, wtr, error.type, ar1.resids, constrain.sign) ci.GPP <- stats::quantile(boot.results$GPP, c(0.025, 0.975), na.rm=TRUE) # 95% confidence interval ci.R <- stats::quantile(boot.results$R, c(0.025, 0.975), na.rm=TRUE) # 95% confidence interval ci.NEP <- stats::quantile(boot.results$NEP, c(0.025, 0.975), na.rm=TRUE) # 95% confidence interval @@ -227,9 +242,14 @@ mleLoopPE <- function(alpha, doobs, c1, c2, beta, irr, wtr, kz, dosat){ # ==================== # = mle NLL function = # ==================== -mleNllPE <- function(Params, do.obs, do.sat, k.gas, z.mix, irr, wtr){ - c1 <- Params[1] #PAR coeff - c2 <- Params[2] #log(Temp) coeff +mleNllPE <- function(Params, do.obs, do.sat, k.gas, z.mix, irr, wtr, constrain.sign){ + if (constrain.sign == TRUE){ + c1 <- exp(Params[1]) #PAR coeff + c2 <- -exp(Params[2]) #log(Temp) coeff + } else { + c1 <- Params[1] #PAR coeff + c2 <- Params[2] #log(Temp) coeff + } Q <- exp(Params[3]) # Variance of the process error # See KalmanDO_smooth.R comments for explanation of beta @@ -253,9 +273,14 @@ mleNllPE <- function(Params, do.obs, do.sat, k.gas, z.mix, irr, wtr){ # ==================== # = mle NLL function = # ==================== -mleNllOE <- function(Params, do.obs, do.sat, k.gas, z.mix, irr, wtr, error.type){ - c1 <- Params[1] #PAR coeff - c2 <- Params[2] #log(Temp) coeff +mleNllOE <- function(Params, do.obs, do.sat, k.gas, z.mix, irr, wtr, error.type, constrain.sign){ + if (constrain.sign == TRUE){ + c1 <- exp(Params[1]) #PAR coeff + c2 <- -exp(Params[2]) #log(Temp) coeff + } else { + c1 <- Params[1] #PAR coeff + c2 <- Params[2] #log(Temp) coeff + } Q <- exp(Params[3]) # Variance of the process error # See KalmanDO_smooth.R comments for explanation of beta @@ -274,7 +299,7 @@ mleNllOE <- function(Params, do.obs, do.sat, k.gas, z.mix, irr, wtr, error.type) # ========================================= # = Function to bootstrap for uncertainty = # ========================================= -bootstrap.mle <- function(n.boot, do.obs, pars, guesses, do.sat, k.gas, freq, z.mix, irr, wtr, error.type = "OE", ar1.resids = FALSE){ +bootstrap.mle <- function(n.boot, do.obs, pars, guesses, do.sat, k.gas, freq, z.mix, irr, wtr, error.type = "OE", ar1.resids = FALSE, constrain.sign = TRUE){ ## Calculate fitted values and residuals if (error.type == "OE"){ @@ -339,16 +364,24 @@ bootstrap.mle <- function(n.boot, do.obs, pars, guesses, do.sat, k.gas, freq, z. doSim <- alpha + simRes if (error.type == "OE"){ - simFit <- optim(guesses, fn = mleNllOE, do.obs = doSim, do.sat = do.sat, k.gas = k.gas / freq, z.mix = z.mix, irr = irr, wtr = wtr) + simFit <- optim(guesses, fn = mleNllOE, do.obs = doSim, do.sat = do.sat, k.gas = k.gas / freq, z.mix = z.mix, irr = irr, wtr = wtr, constrain.sign = constrain.sign) } else if (error.type == "PE"){ - simFit <- optim(guesses, fn = mleNllPE, do.obs = doSim, do.sat = do.sat, k.gas = k.gas / freq, z.mix = z.mix, irr = irr, wtr = wtr) + simFit <- optim(guesses, fn = mleNllPE, do.obs = doSim, do.sat = do.sat, k.gas = k.gas / freq, z.mix = z.mix, irr = irr, wtr = wtr, constrain.sign = constrain.sign) } + simFitPar <- simFit$par - simFitGPP <- mean(exp(simFitPar[1]) * irr, na.rm = TRUE) * freq - simFitR <- mean(-exp(simFitPar[2]) * log(wtr), na.rm = TRUE) * freq + if (constrain.sign == TRUE){ + simFit.gppCoeff <- exp(simFitPar[1]) + simFit.rCoeff <- -exp(simFitPar[2]) + } else { + simFit.gppCoeff <- simFitPar[1] + simFit.rCoeff <- simFitPar[2] + } + simFitGPP <- mean(simFit.gppCoeff * irr, na.rm = TRUE) * freq + simFitR <- mean(simFit.rCoeff * log(wtr), na.rm = TRUE) * freq - boot.results[i, "gppCoeff"] <- exp(simFitPar[1]) - boot.results[i, "rCoeff"] <- -exp(simFitPar[2]) + boot.results[i, "gppCoeff"] <- simFit.gppCoeff + boot.results[i, "rCoeff"] <- simFit.rCoeff boot.results[i, "Q"] <- exp(simFitPar[3]) boot.results[i, "doInit"] <- simFitPar[4] boot.results[i, "convergence"] <- simFit$convergence diff --git a/man/metab.kalman.Rd b/man/metab.kalman.Rd index c2ba52c..50a10db 100644 --- a/man/metab.kalman.Rd +++ b/man/metab.kalman.Rd @@ -4,7 +4,17 @@ \alias{metab.kalman} \title{Metabolism calculated from parameters estimated using a Kalman filter} \usage{ -metab.kalman(do.obs, do.sat, k.gas, z.mix, irr, wtr, n.boot = 0, ...) +metab.kalman( + do.obs, + do.sat, + k.gas, + z.mix, + irr, + wtr, + n.boot = 0, + constrain.sign = TRUE, + ... +) } \arguments{ \item{do.obs}{Vector of dissovled oxygen concentration observations, \eqn{mg O[2] L^{-1}}{mg O2 / L}} @@ -22,6 +32,8 @@ metab.kalman(do.obs, do.sat, k.gas, z.mix, irr, wtr, n.boot = 0, ...) \item{n.boot}{Numeric of the number of bootstrap iterations. Set to 0 to bypass uncertainty quantification (default). Set to 2 or greater to quantify uncertainty via bootstrapping.} +\item{constrain.sign}{Logical. If TRUE, estimates of GPP and R coefficients with be constrained to positive and negative values, respectively.} + \item{...}{additional arguments; currently "datetime" is the only recognized argument passed through \code{...}} } \value{ diff --git a/man/metab.mle.Rd b/man/metab.mle.Rd index c50bebb..e220b09 100644 --- a/man/metab.mle.Rd +++ b/man/metab.mle.Rd @@ -14,6 +14,7 @@ metab.mle( error.type = "OE", n.boot = 0, ar1.resids = TRUE, + constrain.sign = TRUE, ... ) } @@ -37,6 +38,8 @@ metab.mle( \item{ar1.resids}{Logical. If TRUE, the AR(1) structure of the residuals will be retained during randomization when bootstrapping.} +\item{constrain.sign}{Logical. If TRUE, estimates of GPP and R coefficients with be constrained to positive and negative values, respectively.} + \item{...}{additional arguments; currently "datetime" is the only recognized argument passed through \code{...}} } \value{ From 98a8ea47d14d7918ef6dd99c2ef025bdd52b915d Mon Sep 17 00:00:00 2001 From: Bennett McAfee <105509324+bmcafee@users.noreply.github.com> Date: Fri, 17 Jul 2026 10:31:09 -0400 Subject: [PATCH 05/13] Add warning for kalman bootstrapping --- R/metab.kalman.R | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/R/metab.kalman.R b/R/metab.kalman.R index 322afe7..221ce57 100644 --- a/R/metab.kalman.R +++ b/R/metab.kalman.R @@ -131,6 +131,15 @@ metab.kalman <- function(do.obs, do.sat, k.gas, z.mix, irr, wtr, n.boot = 0, con freq <- nobs } + # Additional warning about bootstrapping + if (n.boot == 0){ + # Do nothing + } else if (n.boot >= 2){ + warning("Bootstrapping with the Kalman filter is experimental. Please consider your confidence intervals carefully.") + } else { + stop("n.boot must be 0 (no uncertainty quantification) or greater than 2 (uncertainty quantification via bootstrapping).") + } + # Filter and fit if (constrain.sign == TRUE){ guesses <- c(log(1E-4),log(1E-4),log(5),log(5)) From 1305f3419fce9b0e2c797378ebbcec99456038c5 Mon Sep 17 00:00:00 2001 From: Bennett McAfee <105509324+bmcafee@users.noreply.github.com> Date: Fri, 17 Jul 2026 10:34:11 -0400 Subject: [PATCH 06/13] Change default constrain.sign to FALSE --- R/metab.kalman.R | 2 +- R/metab.mle.R | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/R/metab.kalman.R b/R/metab.kalman.R index 221ce57..1d4587e 100644 --- a/R/metab.kalman.R +++ b/R/metab.kalman.R @@ -101,7 +101,7 @@ #' do.sat=do.sat, wtr=wtr[,2], #' k.gas=k.gas, do.obs=doobs[,2]) #'@export -metab.kalman <- function(do.obs, do.sat, k.gas, z.mix, irr, wtr, n.boot = 0, constrain.sign = TRUE, ...){ +metab.kalman <- function(do.obs, do.sat, k.gas, z.mix, irr, wtr, n.boot = 0, constrain.sign = FALSE, ...){ complete.inputs(do.obs=do.obs, do.sat=do.sat, k.gas=k.gas, z.mix=z.mix, irr=irr, wtr=wtr, error=TRUE) diff --git a/R/metab.mle.R b/R/metab.mle.R index e1dc6c7..58b9fe7 100644 --- a/R/metab.mle.R +++ b/R/metab.mle.R @@ -299,7 +299,7 @@ mleNllOE <- function(Params, do.obs, do.sat, k.gas, z.mix, irr, wtr, error.type, # ========================================= # = Function to bootstrap for uncertainty = # ========================================= -bootstrap.mle <- function(n.boot, do.obs, pars, guesses, do.sat, k.gas, freq, z.mix, irr, wtr, error.type = "OE", ar1.resids = FALSE, constrain.sign = TRUE){ +bootstrap.mle <- function(n.boot, do.obs, pars, guesses, do.sat, k.gas, freq, z.mix, irr, wtr, error.type = "OE", ar1.resids = FALSE, constrain.sign = FALSE){ ## Calculate fitted values and residuals if (error.type == "OE"){ From caa4bdd59577b4138de14de15d20d34acb948824 Mon Sep 17 00:00:00 2001 From: Bennett McAfee <105509324+bmcafee@users.noreply.github.com> Date: Fri, 17 Jul 2026 10:38:06 -0400 Subject: [PATCH 07/13] Rebuild new default constrain.sign --- R/metab.mle.R | 2 +- man/metab.kalman.Rd | 2 +- man/metab.mle.Rd | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/R/metab.mle.R b/R/metab.mle.R index 58b9fe7..26bbcb0 100644 --- a/R/metab.mle.R +++ b/R/metab.mle.R @@ -110,7 +110,7 @@ #' #'metab.mle(doobs[,2], do.sat, k.gas, z.mix[,2], irr[,2], wtr[,3]) #'@export -metab.mle <- function(do.obs, do.sat, k.gas, z.mix, irr, wtr, error.type="OE", n.boot = 0, ar1.resids = TRUE, constrain.sign = TRUE, ...){ +metab.mle <- function(do.obs, do.sat, k.gas, z.mix, irr, wtr, error.type="OE", n.boot = 0, ar1.resids = TRUE, constrain.sign = FALSE, ...){ complete.inputs(do.obs=do.obs, do.sat=do.sat, k.gas=k.gas, z.mix=z.mix, irr=irr, wtr=wtr, error=TRUE) diff --git a/man/metab.kalman.Rd b/man/metab.kalman.Rd index 50a10db..2d67e5a 100644 --- a/man/metab.kalman.Rd +++ b/man/metab.kalman.Rd @@ -12,7 +12,7 @@ metab.kalman( irr, wtr, n.boot = 0, - constrain.sign = TRUE, + constrain.sign = FALSE, ... ) } diff --git a/man/metab.mle.Rd b/man/metab.mle.Rd index e220b09..3a3a751 100644 --- a/man/metab.mle.Rd +++ b/man/metab.mle.Rd @@ -14,7 +14,7 @@ metab.mle( error.type = "OE", n.boot = 0, ar1.resids = TRUE, - constrain.sign = TRUE, + constrain.sign = FALSE, ... ) } From a4cef0ca20cdbddeaa72bfb9d85a67cb6dced135 Mon Sep 17 00:00:00 2001 From: Bennett McAfee <105509324+bmcafee@users.noreply.github.com> Date: Wed, 5 Aug 2026 10:27:46 -0400 Subject: [PATCH 08/13] Update R/metab.kalman.R with error line break Co-authored-by: Jake Zwart --- R/metab.kalman.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/metab.kalman.R b/R/metab.kalman.R index 1d4587e..db0b209 100644 --- a/R/metab.kalman.R +++ b/R/metab.kalman.R @@ -187,7 +187,7 @@ metab.kalman <- function(do.obs, do.sat, k.gas, z.mix, irr, wtr, n.boot = 0, con stop("n.boot must be 0 (no uncertainty quantification) or greater than 2 (uncertainty quantification via bootstrapping).") } - + return(results) } From 492b572744b503db4a06ed717510ec667ffcf8cb Mon Sep 17 00:00:00 2001 From: Bennett McAfee <105509324+bmcafee@users.noreply.github.com> Date: Wed, 5 Aug 2026 10:28:30 -0400 Subject: [PATCH 09/13] Update R/metab.kalman.R bootstrap error message Co-authored-by: Jake Zwart --- R/metab.kalman.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/metab.kalman.R b/R/metab.kalman.R index db0b209..a1ba51c 100644 --- a/R/metab.kalman.R +++ b/R/metab.kalman.R @@ -137,7 +137,7 @@ metab.kalman <- function(do.obs, do.sat, k.gas, z.mix, irr, wtr, n.boot = 0, con } else if (n.boot >= 2){ warning("Bootstrapping with the Kalman filter is experimental. Please consider your confidence intervals carefully.") } else { - stop("n.boot must be 0 (no uncertainty quantification) or greater than 2 (uncertainty quantification via bootstrapping).") + stop("n.boot must be 0 (no uncertainty quantification) or 2 or greater (uncertainty quantification via bootstrapping).") } # Filter and fit From e4fe982278df1f1cca61fdb056eb53d7e552d8d7 Mon Sep 17 00:00:00 2001 From: Bennett McAfee <105509324+bmcafee@users.noreply.github.com> Date: Wed, 5 Aug 2026 10:29:12 -0400 Subject: [PATCH 10/13] Update R/metab.R documentation Co-authored-by: Jake Zwart --- R/metab.R | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/R/metab.R b/R/metab.R index aea9bc3..615c261 100644 --- a/R/metab.R +++ b/R/metab.R @@ -7,7 +7,8 @@ #' #' #'@usage -#'metab(data, method, verbose = TRUE, wtr.name="wtr", irr.name="irr", do.obs.name="do.obs", ...) +#'metab(data, method, verbose = TRUE, wtr.name="wtr", irr.name="irr", +#' do.obs.name="do.obs", ...) #' #'@param data #' a data.frame whose columns are From f081a574c396e0043c575053be894765e32cd41e Mon Sep 17 00:00:00 2001 From: Bennett McAfee <105509324+bmcafee@users.noreply.github.com> Date: Wed, 5 Aug 2026 15:29:45 -0400 Subject: [PATCH 11/13] AR1 fix + Documentation updates --- R/metab.R | 2 +- R/metab.kalman.R | 6 +++--- R/metab.mle.R | 8 ++++---- man/metab.Rd | 3 ++- man/metab.kalman.Rd | 6 +++--- man/metab.mle.Rd | 6 +++--- 6 files changed, 16 insertions(+), 15 deletions(-) diff --git a/R/metab.R b/R/metab.R index 615c261..64c8917 100644 --- a/R/metab.R +++ b/R/metab.R @@ -109,7 +109,7 @@ #' #'@export -metab <- function(data, method = NULL, verbose = TRUE, wtr.name="wtr", irr.name="irr", do.obs.name="do.obs", ...){ +metab <- function(data, method = NULL, verbose = interactive(), wtr.name="wtr", irr.name="irr", do.obs.name="do.obs", ...){ m.args <- list(...) diff --git a/R/metab.kalman.R b/R/metab.kalman.R index a1ba51c..2862603 100644 --- a/R/metab.kalman.R +++ b/R/metab.kalman.R @@ -8,7 +8,7 @@ #'@param z.mix Vector of mixed-layer depths in meters. To calculate, see \link[rLakeAnalyzer]{ts.meta.depths} #'@param irr Vector of photosynthetically active radiation in \eqn{\mu mol\ m^{-2} s^{-1}}{micro mols / m^2 / s} #'@param wtr Vector of water temperatures in \eqn{^{\circ}C}{degrees C}. Used in scaling respiration with temperature -#'@param n.boot Numeric of the number of bootstrap iterations. Set to 0 to bypass uncertainty quantification (default). Set to 2 or greater to quantify uncertainty via bootstrapping. +#'@param n.boot Numeric of the number of bootstrap iterations. Set to 0 to bypass uncertainty quantification (default). Set to 2 or greater to quantify uncertainty via bootstrapping. To ensure reproducible confidence intervals, use \link{set.seed}. #'@param constrain.sign Logical. If TRUE, estimates of GPP and R coefficients with be constrained to positive and negative values, respectively. #'@param ... additional arguments; currently "datetime" is the only recognized argument passed through \code{...} #'@return @@ -25,10 +25,10 @@ #'\item{GPP.cv}{coefficient of variation of the bootstrapped estimates of Gross Primary Production} #'\item{R.lci}{lower 95\% confidence interval of the estimate of Respiration, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} #'\item{R.uci}{upper 95\% confidence interval of the estimate of Respiration, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} -#'\item{R.cv}{coefficient of variation of the bootstrapped estimates of Respiration} +#'\item{R.cv}{coefficient of variation of the bootstrapped estimates of Respiration. As R is typically negative, the CV is also negative.} #'\item{NEP.lci}{lower 95\% confidence interval of the estimate of Net Ecosystem Production, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} #'\item{NEP.uci}{upper 95\% confidence interval of the estimate of Net Ecosystem Production, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} -#'\item{NEP.cv}{coefficient of variation of the bootstrapped estimates of Net Ecosystem Production} +#'\item{NEP.cv}{coefficient of variation of the bootstrapped estimates of Net Ecosystem Production. As the mean NEP can approach 0, the CV can become inflated.} #'} #' #' diff --git a/R/metab.mle.R b/R/metab.mle.R index 26bbcb0..c4611d5 100644 --- a/R/metab.mle.R +++ b/R/metab.mle.R @@ -8,7 +8,7 @@ #'@param irr Vector of photosynthetically active radiation in \eqn{\mu mol\ m^{-2} s^{-1}}{micro mols / m^2 / s} #'@param wtr Vector of water temperatures in \eqn{^{\circ}C}{degrees C}. Used in scaling respiration with temperature #'@param error.type Option specifying if model should assume pure Process Error 'PE' or Observation Error 'OE'. Defaults to observation error 'OE'. -#'@param n.boot Numeric of the number of bootstrap iterations. Set to 0 to bypass uncertainty quantification (default). Set to 2 or greater to quantify uncertainty via bootstrapping. +#'@param n.boot Numeric of the number of bootstrap iterations. Set to 0 to bypass uncertainty quantification (default). Set to 2 or greater to quantify uncertainty via bootstrapping. To ensure reproducible confidence intervals, use \link{set.seed}. #'@param ar1.resids Logical. If TRUE, the AR(1) structure of the residuals will be retained during randomization when bootstrapping. #'@param constrain.sign Logical. If TRUE, estimates of GPP and R coefficients with be constrained to positive and negative values, respectively. #'@param ... additional arguments; currently "datetime" is the only recognized argument passed through \code{...} @@ -26,10 +26,10 @@ #'\item{GPP.cv}{coefficient of variation of the bootstrapped estimates of Gross Primary Production} #'\item{R.lci}{lower 95\% confidence interval of the estimate of Respiration, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} #'\item{R.uci}{upper 95\% confidence interval of the estimate of Respiration, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} -#'\item{R.cv}{coefficient of variation of the bootstrapped estimates of Respiration} +#'\item{R.cv}{coefficient of variation of the bootstrapped estimates of Respiration. As R is typically negative, the CV is also negative.} #'\item{NEP.lci}{lower 95\% confidence interval of the estimate of Net Ecosystem Production, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} #'\item{NEP.uci}{upper 95\% confidence interval of the estimate of Net Ecosystem Production, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} -#'\item{NEP.cv}{coefficient of variation of the bootstrapped estimates of Net Ecosystem Production} +#'\item{NEP.cv}{coefficient of variation of the bootstrapped estimates of Net Ecosystem Production. As the mean NEP can approach 0, the CV can become inflated.} #'} #' The maximum likelihood estimates of model parameters can be accessed via \code{attributes(metab.mle(...))[["params"]]} #' @@ -331,7 +331,7 @@ bootstrap.mle <- function(n.boot, do.obs, pars, guesses, do.sat, k.gas, freq, z. n.obs <- length(do.obs) if (ar1.resids){ - 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) ar1.coeff <- as.numeric(ar1.lm$coefficients) ar1.sd <- sd(ar1.lm$residuals) } diff --git a/man/metab.Rd b/man/metab.Rd index 6a8c256..34cd853 100644 --- a/man/metab.Rd +++ b/man/metab.Rd @@ -4,7 +4,8 @@ \alias{metab} \title{Calculate metabolism} \usage{ -metab(data, method, verbose = TRUE, wtr.name="wtr", irr.name="irr", do.obs.name="do.obs", ...) +metab(data, method, verbose = TRUE, wtr.name="wtr", irr.name="irr", + do.obs.name="do.obs", ...) } \arguments{ \item{data}{a data.frame whose columns are diff --git a/man/metab.kalman.Rd b/man/metab.kalman.Rd index 2d67e5a..16be2e3 100644 --- a/man/metab.kalman.Rd +++ b/man/metab.kalman.Rd @@ -30,7 +30,7 @@ metab.kalman( \item{wtr}{Vector of water temperatures in \eqn{^{\circ}C}{degrees C}. Used in scaling respiration with temperature} -\item{n.boot}{Numeric of the number of bootstrap iterations. Set to 0 to bypass uncertainty quantification (default). Set to 2 or greater to quantify uncertainty via bootstrapping.} +\item{n.boot}{Numeric of the number of bootstrap iterations. Set to 0 to bypass uncertainty quantification (default). Set to 2 or greater to quantify uncertainty via bootstrapping. To ensure reproducible confidence intervals, use \link{set.seed}.} \item{constrain.sign}{Logical. If TRUE, estimates of GPP and R coefficients with be constrained to positive and negative values, respectively.} @@ -50,10 +50,10 @@ If n.boot is 2 or greater, the data frame will also contain the following column \item{GPP.cv}{coefficient of variation of the bootstrapped estimates of Gross Primary Production} \item{R.lci}{lower 95\% confidence interval of the estimate of Respiration, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} \item{R.uci}{upper 95\% confidence interval of the estimate of Respiration, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} -\item{R.cv}{coefficient of variation of the bootstrapped estimates of Respiration} +\item{R.cv}{coefficient of variation of the bootstrapped estimates of Respiration. As R is typically negative, the CV is also negative.} \item{NEP.lci}{lower 95\% confidence interval of the estimate of Net Ecosystem Production, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} \item{NEP.uci}{upper 95\% confidence interval of the estimate of Net Ecosystem Production, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} -\item{NEP.cv}{coefficient of variation of the bootstrapped estimates of Net Ecosystem Production} +\item{NEP.cv}{coefficient of variation of the bootstrapped estimates of Net Ecosystem Production. As the mean NEP can approach 0, the CV can become inflated.} } diff --git a/man/metab.mle.Rd b/man/metab.mle.Rd index 3a3a751..9839ddb 100644 --- a/man/metab.mle.Rd +++ b/man/metab.mle.Rd @@ -34,7 +34,7 @@ metab.mle( \item{error.type}{Option specifying if model should assume pure Process Error 'PE' or Observation Error 'OE'. Defaults to observation error 'OE'.} -\item{n.boot}{Numeric of the number of bootstrap iterations. Set to 0 to bypass uncertainty quantification (default). Set to 2 or greater to quantify uncertainty via bootstrapping.} +\item{n.boot}{Numeric of the number of bootstrap iterations. Set to 0 to bypass uncertainty quantification (default). Set to 2 or greater to quantify uncertainty via bootstrapping. To ensure reproducible confidence intervals, use \link{set.seed}.} \item{ar1.resids}{Logical. If TRUE, the AR(1) structure of the residuals will be retained during randomization when bootstrapping.} @@ -56,10 +56,10 @@ If n.boot is 2 or greater, the data frame will also contain the following column \item{GPP.cv}{coefficient of variation of the bootstrapped estimates of Gross Primary Production} \item{R.lci}{lower 95\% confidence interval of the estimate of Respiration, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} \item{R.uci}{upper 95\% confidence interval of the estimate of Respiration, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} -\item{R.cv}{coefficient of variation of the bootstrapped estimates of Respiration} +\item{R.cv}{coefficient of variation of the bootstrapped estimates of Respiration. As R is typically negative, the CV is also negative.} \item{NEP.lci}{lower 95\% confidence interval of the estimate of Net Ecosystem Production, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} \item{NEP.uci}{upper 95\% confidence interval of the estimate of Net Ecosystem Production, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} -\item{NEP.cv}{coefficient of variation of the bootstrapped estimates of Net Ecosystem Production} +\item{NEP.cv}{coefficient of variation of the bootstrapped estimates of Net Ecosystem Production. As the mean NEP can approach 0, the CV can become inflated.} } The maximum likelihood estimates of model parameters can be accessed via \code{attributes(metab.mle(...))[["params"]]} } From edf4d7fb35047388e522aec7583d905aaaf55f38 Mon Sep 17 00:00:00 2001 From: Bennett McAfee <105509324+bmcafee@users.noreply.github.com> Date: Wed, 5 Aug 2026 15:42:14 -0400 Subject: [PATCH 12/13] Revert metab method argument default --- R/metab.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/metab.R b/R/metab.R index 64c8917..44ccc0a 100644 --- a/R/metab.R +++ b/R/metab.R @@ -109,7 +109,7 @@ #' #'@export -metab <- function(data, method = NULL, verbose = interactive(), wtr.name="wtr", irr.name="irr", do.obs.name="do.obs", ...){ +metab <- function(data, method, verbose = interactive(), wtr.name="wtr", irr.name="irr", do.obs.name="do.obs", ...){ m.args <- list(...) From 87c94604f1d8269c831d192e8f144d057fdd99e6 Mon Sep 17 00:00:00 2001 From: Bennett McAfee <105509324+bmcafee@users.noreply.github.com> Date: Wed, 5 Aug 2026 15:49:14 -0400 Subject: [PATCH 13/13] Update metab documentation --- R/metab.R | 4 ++-- man/metab.Rd | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/R/metab.R b/R/metab.R index 44ccc0a..b79017d 100644 --- a/R/metab.R +++ b/R/metab.R @@ -7,7 +7,7 @@ #' #' #'@usage -#'metab(data, method, verbose = TRUE, wtr.name="wtr", irr.name="irr", +#'metab(data, method, verbose = interactive(), wtr.name="wtr", irr.name="irr", #' do.obs.name="do.obs", ...) #' #'@param data @@ -109,7 +109,7 @@ #' #'@export -metab <- function(data, method, verbose = interactive(), wtr.name="wtr", irr.name="irr", do.obs.name="do.obs", ...){ +metab <- function(data, method = NULL, verbose = interactive(), wtr.name="wtr", irr.name="irr", do.obs.name="do.obs", ...){ m.args <- list(...) diff --git a/man/metab.Rd b/man/metab.Rd index 34cd853..9ade0e4 100644 --- a/man/metab.Rd +++ b/man/metab.Rd @@ -4,7 +4,7 @@ \alias{metab} \title{Calculate metabolism} \usage{ -metab(data, method, verbose = TRUE, wtr.name="wtr", irr.name="irr", +metab(data, method, verbose = interactive(), wtr.name="wtr", irr.name="irr", do.obs.name="do.obs", ...) } \arguments{