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..b79017d 100644 --- a/R/metab.R +++ b/R/metab.R @@ -7,7 +7,8 @@ #' #' #'@usage -#'metab(data, method, wtr.name="wtr", irr.name="irr", do.obs.name="do.obs", ...) +#'metab(data, method, verbose = interactive(), wtr.name="wtr", irr.name="irr", +#' do.obs.name="do.obs", ...) #' #'@param data #' a data.frame whose columns are @@ -31,6 +32,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 +61,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 +109,7 @@ #' #'@export -metab <- function(data, method = NULL, 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(...) @@ -163,6 +165,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 +178,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.kalman.R b/R/metab.kalman.R index 6eab269..2862603 100644 --- a/R/metab.kalman.R +++ b/R/metab.kalman.R @@ -8,6 +8,8 @@ #'@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. 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 #'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. 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. As the mean NEP can approach 0, the CV can become inflated.} +#'} #' #' #' Use \link{attributes} to access more model output: @@ -51,10 +65,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 +101,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, 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) - + 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 +130,64 @@ 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 } - + + # 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 2 or greater (uncertainty quantification via bootstrapping).") + } + # 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) + 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"=pars0[1], "rCoeff"=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 - 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, 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 + 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).") + } + + return(results) } @@ -138,80 +197,85 @@ metab.kalman <- function(do.obs, do.sat, k.gas, z.mix, irr, wtr, ...){ # = Kalman filter/ nll = # ====================== # 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 +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 + 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 + + # 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 +284,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 +418,104 @@ 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, constrain.sign, 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, constrain.sign=constrain.sign) + simFitPar <- simFit$par + 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"] <- 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 + 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/R/metab.mle.R b/R/metab.mle.R index 4df7de1..c4611d5 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,9 @@ #'@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. 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{...} #'@return #'A data.frame with columns corresponding to components of metabolism @@ -16,6 +19,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. 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. 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"]]} #' #'@details @@ -46,19 +61,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}. #' -#'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. +#'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}. #' -#'@importFrom stats dnorm optim +#'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 rnorm quantile #' #'@seealso #'\link{metab}, \link{metab.bookkeep}, \link{metab.ols}, \link{metab.kalman}, \link{metab.bayesian} @@ -87,8 +110,8 @@ #' #'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, 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) @@ -101,65 +124,101 @@ metab.mle <- function(do.obs, do.sat, k.gas, z.mix, irr, wtr, error.type="OE", . 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 - if(nobs!=freq){ # nobs and freq should agree, if they don't issue a warning - bad.date <- format.Date(datetime[1], format="%Y-%m-%d") - warning("number of observations on ", bad.date, " (", nobs, ") ", "does not equal estimated sampling frequency", " (", freq, ")", sep="") - } - }else{ # if datetime is *not* in the ... args - warning("datetime not found, inferring sampling frequency from # of observations") # issue a warning (note checks in addNAs) - # NOTE: because of the checks in addNA's, it is unlikely a user would receive this warning via metab() - # warning will only be seen through direct use of metab.bookkeep when datettime is not supplied - freq <- nobs - } - - chk.list <- list(do.obs, irr, do.sat, z.mix, k.gas, wtr) - if(!all(sapply(chk.list, is.numeric)) || !all(sapply(chk.list, is.vector))){ - stop('All metab.mle inputs must be numeric vectors.') - } - - if(!all(nobs==sapply(chk.list, length))){ - stop('All input data to metab.mle must be the same length') - } - - Q0 <- ((diff(range(do.obs,na.rm=TRUE)) - mean(do.obs,na.rm=TRUE))^2 / length(do.obs)) - - 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) - - pars0 <- fit$par - 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) - - pars0 <- fit$par - pars <- c("gppCoeff"=pars0[1], "rCoeff"=pars0[2], "Q"=exp(pars0[3]), "nll"=fit$value) - - }else{ - stop("error.type must be either 'OE' or 'PE', Observation Error or Process Error respectively.") - } - - # ==================================== - # = 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("params"=pars, "metab"=c("GPP"=GPP,"R"=R,"NEP"=GPP+R))) + 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 + if(nobs!=freq){ # nobs and freq should agree, if they don't issue a warning + bad.date <- format.Date(datetime[1], format="%Y-%m-%d") + warning("number of observations on ", bad.date, " (", nobs, ") ", "does not equal estimated sampling frequency", " (", freq, ")", sep="") + } + }else{ # if datetime is *not* in the ... args + warning("datetime not found, inferring sampling frequency from # of observations") # issue a warning (note checks in addNAs) + # NOTE: because of the checks in addNA's, it is unlikely a user would receive this warning via metab() + # warning will only be seen through direct use of metab.bookkeep when datettime is not supplied + freq <- nobs + } + + chk.list <- list(do.obs, irr, do.sat, z.mix, k.gas, wtr) + if(!all(sapply(chk.list, is.numeric)) || !all(sapply(chk.list, is.vector))){ + stop('All metab.mle inputs must be numeric vectors.') + } + + if(!all(nobs==sapply(chk.list, length))){ + stop('All input data to metab.mle must be the same length') + } + + Q0 <- ((diff(range(do.obs,na.rm=TRUE)) - mean(do.obs,na.rm=TRUE))^2 / length(do.obs)) + + 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, constrain.sign = constrain.sign) + + pars0 <- fit$par + 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, constrain.sign = constrain.sign) + + pars0 <- fit$par + 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.") + } + + # ==================================== + # = Use fits to calculate metabolism = + # ==================================== + GPP <- mean(pars[1]*irr, na.rm=TRUE) * freq + R <- mean(pars[2]*log(wtr), na.rm=TRUE) * freq + + # ================================= + # = 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, 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 + 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) } # ============================================ @@ -183,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 @@ -209,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 @@ -226,3 +295,101 @@ 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, constrain.sign = 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[2:n.obs] ~ resids[1:(n.obs - 1)] - 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, 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, constrain.sign = constrain.sign) + } + + simFitPar <- simFit$par + 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"] <- 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 + 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..9ade0e4 100644 --- a/man/metab.Rd +++ b/man/metab.Rd @@ -4,7 +4,8 @@ \alias{metab} \title{Calculate metabolism} \usage{ -metab(data, method, wtr.name="wtr", irr.name="irr", do.obs.name="do.obs", ...) +metab(data, method, verbose = interactive(), wtr.name="wtr", irr.name="irr", + do.obs.name="do.obs", ...) } \arguments{ \item{data}{a data.frame whose columns are @@ -28,6 +29,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 +104,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.kalman.Rd b/man/metab.kalman.Rd index b28d3fb..16be2e3 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, ...) +metab.kalman( + do.obs, + do.sat, + k.gas, + z.mix, + irr, + wtr, + n.boot = 0, + constrain.sign = FALSE, + ... +) } \arguments{ \item{do.obs}{Vector of dissovled oxygen concentration observations, \eqn{mg O[2] L^{-1}}{mg O2 / L}} @@ -20,6 +30,10 @@ 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. 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.} + \item{...}{additional arguments; currently "datetime" is the only recognized argument passed through \code{...}} } \value{ @@ -29,6 +43,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. 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. As the mean NEP can approach 0, the CV can become inflated.} +} Use \link{attributes} to access more model output: @@ -99,11 +125,14 @@ 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} } \author{ -Ryan Batt, Luke A. Winslow +Ryan Batt, Luke A. Winslow, Bennett McAfee } diff --git a/man/metab.mle.Rd b/man/metab.mle.Rd index 22fb799..9839ddb 100644 --- a/man/metab.mle.Rd +++ b/man/metab.mle.Rd @@ -4,10 +4,22 @@ \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, + constrain.sign = FALSE, + ... +) } \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 +34,12 @@ 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. 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.} + +\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{ @@ -31,6 +49,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. 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. 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"]]} } \description{ @@ -92,19 +122,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 }