# R/metab.bayesian.R In LakeMetabolizer: Tools for the Analysis of Ecosystem Metabolism

#### Documented in metab.bayesian

if(getRversion() >= "2.15.1")  utils::globalVariables("jags")

# ====================================
# = Function to write the jags model =
# ====================================
bayes.makeModel <- function(k.gas){

if(!requireNamespace("R2jags")){
stop("metab.bayesian requires R2jags.\ninstall.packages('R2jags')\n #Also install JAGS (http://mcmc-jags.sourceforge.net/)")
}

finite.1oK <- is.finite(1/k.gas)

choose.allK <- all(finite.1oK)
choose.noK <- all(!finite.1oK)
choose.bothK <- any(finite.1oK) & any(!finite.1oK)

choice.mod <- c("allK", "noK", "bothK")[c(choose.allK, choose.noK, choose.bothK)]

# Write the appropriate bayesian model into a temporary file
#modfile <- tempfile('jags.metab.bayes')
jags.dir = system.file('jags', package='LakeMetabolizer')
switch(choice.mod,
allK = {modfile = file.path(jags.dir, 'bayes.mod.allk.JAGS'); print("using allK model")},
noK = {modfile = file.path(jags.dir, 'bayes.mod.nok.JAGS'); print("using noK model")},
bothK = {modfile = file.path(jags.dir, 'bayes.mod.bothk.JAGS'); print("using bothK model")}
)

return(modfile)
}
#
#
# # ==================================
# # = Bayes model for all non-zero K =
# # ==================================
# bayes.mod.allK <- function(){
#
# 	# model process
# 	for(i in 2:N){
# 		Y[i] ~ dnorm(a[i], tauV) # observations (Y) are distributed with mean equivalent to true values, and precision tauV (1/tauV is variance of observation error)
# 		K[i-1] ~ dnorm(kP[i-1, 1], 1/kP[i-1, 2]) #distributin on K
#
# 		kz[i-1] <- K[i-1]/Zmix[i-1]
#
# 		a1[i] <- U[i-1,]%*%C + kz[i-1]*satO[i-1]
# 		aHat[i] <- a1[i]/kz[i-1] + -exp(-1*kz[i-1])*a1[i]/kz[i-1] + exp(-kz[i-1])*a[i-1]
#
# 		a[i] ~ dnorm(aHat[i], tauW) # true values have a mean equivalent to estimated values, but accompanied by process error (process precision is tauW)
# 	}
#
# 	# Starting values
# 	a <- a0
#
# 	#Priors on regression coefficients
# 	C ~ dnorm(cP[1,1], 1/cP[1,2])
# 	C ~ dnorm(cP[2,1], 1/cP[2,2])
#
# 	#Prior on errors
# 	tauV ~ dgamma(1.0E-3, 1.0E-3)
# 	tauW ~ dgamma(1.0E-3, 1.0E-3)
# 	sigmaV <- 1/sqrt(tauV)
# 	sigmaW <- 1/sqrt(tauW)
# }
#
# # ===============================
# # = Bayes model w/ K turned off =
# # ===============================
# bayes.mod.noK <- function(){
#
# 	# model process
# 	for(i in 2:N){
# 		Y[i] ~ dnorm(a[i], tauV) # observations (Y) are distributed with mean equivalent to true values, and precision tauV (1/tauV is variance of observation error)
# 		# NOT USED
# 		K[i-1] ~ dnorm(kP[i-1, 1], 1/kP[i-1, 2]) #distribution on K
#
# 		aHat[i] <- a[i-1] + U[i-1,]%*%C # the process
#
# 		a[i] ~ dnorm(aHat[i], tauW) # true values have a mean equivalent to estimated values, but accompanied by process error (process precision is tauW)
# 	}
#
# 	# Starting values
# 	a <- a0
#
# 	#Priors on regression coefficients
# 	C ~ dnorm(cP[1,1], 1/cP[1,2])
# 	C ~ dnorm(cP[2,1], 1/cP[2,2])
#
#
# 	#Prior on errors
# 	tauV ~ dgamma(1.0E-3, 1.0E-3)
# 	tauW ~ dgamma(1.0E-3, 1.0E-3)
# 	sigmaV <- 1/sqrt(tauV)
# 	sigmaW <- 1/sqrt(tauW)
# }
#
# # ===============================================
# # = Bayes model that can handle both 0 and !0 K =
# # ===============================================
# bayes.mod.bothK <- function(){
#
# 	# model process
# 	for(i in 2:N){
# 		Y[i] ~ dnorm(a[i], tauV) # observations (Y) are distributed with mean equivalent to true values, and precision tauV (1/tauV is variance of observation error)
# 		K[i-1] ~ dnorm(kP[i-1, 1], 1/kP[i-1, 2]) #distributin on K
#
# 		# jags cannot handle an expression that includes division by 0 (so can't do blah <- ifelse(kz==0, 1, 1/kz), b/c it'll do 1/kz even when kz==0)
# 		# so to avoid division by 0 when k.gas==0, have to make a "safe" kz that is 1 if kz is 0
# 		# when kzSafe is 1 (to avoid division by 0), that means that we just want aHat to be the bio process
# 		# so we have to cancel out all math that is done by a kzSafe==1 (b/c it is bogus) by multiplying by 0, and add bio process
# 		# but if kzSafe!=0, we should multiply the math involving kzSafe by 1 (to keep it), and multiply the added bio process by 0 (to remove it)
# 		# this is a pretty hacked solution, but I don't see another choice (the only control flow in jags is ifelse, no if(){}else{})
#
# 		kz[i-1] <- K[i-1]/Zmix[i-1]
# 		kzSafe[i-1] <- ifelse(kz[i-1]==0, 1, kz[i-1]) # if kz is 0, change to 1
# 		kzCancel[i-1] <- ifelse(kzSafe[i-1]==1, 0, 1) # if kzSafe is 1 (meaning kz is 0), kzCancel needs to be 0
# 		pCancel[i-1] <- ifelse(kzSafe[i-1]==1, 1, 0) # if kzSafe is just kz, then math involving kzSafe is not bogus, and need to cancel the added process
#
# 		a1[i] <- U[i-1,]%*%C + kz[i-1]*satO[i-1]
# 		aHat[i] <- (a1[i]/kzSafe[i-1] + -exp(-1*kzSafe[i-1])*a1[i]/kzSafe[i-1] + exp(-kzSafe[i-1])*a[i-1])*kzCancel[i-1] + (a[i-1] + a1[i])*pCancel[i-1]
#
# 		a[i] ~ dnorm(aHat[i], tauW) # true values have a mean equivalent to estimated values, but accompanied by process error (process precision is tauW)
# 	}
#
#
# 	# Starting values
# 	a <- a0
#
# 	#Priors on regression coefficients
# 	C ~ dnorm(cP[1,1], 1/cP[1,2])
# 	C ~ dnorm(cP[2,1], 1/cP[2,2])
#
# 	#Prior on errors
# 	tauV ~ dgamma(1.0E-3, 1.0E-3)
# 	tauW ~ dgamma(1.0E-3, 1.0E-3)
# 	sigmaV <- 1/sqrt(tauV)
# 	sigmaW <- 1/sqrt(tauW)
# }

# ================================
# = Supply Data and run bayesFit =
# ================================
bayesFit <- function(data, params, mf, tend="median", ...){ #function that writes jags model, traces params, supplies data, etc

bf.args <- list(...)

jags.m <- R2jags::jags(data, NULL, parameters.to.save=params, mf)

tF <- function(x, tend){ # tendency function
switch(tend,
median=median(x), #median
mean=mean(x), #mean
mode1 = unique(x)[which.max(tabulate(match(x, unique(x))))], # mode --- most frequently observed value
mode2 = { # mode --- based on the highest peak of the posterior probability (from density plot)
xd <- density(x)
xd$x[which.max(xd$y)]
}
)
}
# medSim <- matrix(apply(jags2.m$BUGSoutput$sims.matrix, 2, median)[-(1)], nrow=115, ncol=3)
simOut <- jags.m$BUGSoutput$sims.matrix
ctSim <- apply(simOut, 2, tF, tend) # the central tendency metric
sdSim <- apply(simOut, 2, sd)

#Figure out the order of the sims.matrix columns ...
n.obs <- length(data$U[,1]) GPP <- mean(ctSim*data$U[,1], na.rm=TRUE) * n.obs # gpp coef * par, then sum
R <- mean(ctSim*data$U[,2], na.rm=TRUE) * n.obs # r coef * log(temp), then sum GPPsd <- sqrt(sum(sdSim^2*data$U[,1]^2))
Rsd <- sqrt(sum(sdSim^2*data$U[,2]^2)) NEPsd <- sqrt(GPPsd^2 + Rsd^2) return(list( "model" = jags.m, "params" = ctSim[1:2], "metab.sd" = matrix(c(GPPsd, Rsd, NEPsd), nrow=1, dimnames=list(NULL, c("GPPsd", "Rsd", "NEPsd"))), "metab" = matrix(c(GPP, R, GPP+R), nrow=1, dimnames=list(NULL, c("GPP", "R", "NEP"))) )) # need to clean up format, and maybe include a return of the sd's of the estimates } #'@title Metabolism model based on a bayesian parameter estimation framework #'@description This function runs the bayesian metabolism model on the supplied #'gas concentration and other supporting data. This allows for both estimates of #'metabolism along with uncertainty around the parameters. #'@param do.obs Vector of dissovled oxygen concentration observations, mg L^-1 #'@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} #'@param z.mix Vector of mixed-layer depths in meters. To calculate, see \link{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 priors Parameter priors supplied as a named numeric vector (example: c("gppMu"=0, "gppSig2"=1E5, "rMu"=0, "rSig2"=1E5, "kSig2"=NA)) #'@param ... additional arguments; currently "datetime" is the only recognized argument passed through \code{...} #' @return #' A list of length 4 with components: #' \item{model}{the jags model, including posterior draws (see \link[R2jags]{jags})} #' \item{params}{parameter estimates of interest from model (medians)} #' \item{metab.sd}{standard deviation of metabolism estimates} #' \item{metab}{daily metabolism estimates as a data.frame with columns corresponding to #' \describe{ #' \item{\code{GPP}}{numeric estimate of Gross Primary Production, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} #' \item{\code{R}}{numeric estimate of Respiration, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} #' \item{\code{NEP}}{numeric estimate of Net Ecosystem production, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} #' }} #'@references #'Holtgrieve, Gordon W., Daniel E. Schindler, Trevor a. Branch, and Z. Teresa A'mar. #'2010. \emph{Simultaneous Quantification of Aquatic Ecosystem Metabolism and Reaeration #'Using a Bayesian Statistical Model of Oxygen Dynamics}. #'Limnology and Oceanography 55 (3): 1047-1062. doi:10.4319/lo.2010.55.3.1047. #'http://www.aslo.org/lo/toc/vol_55/issue_3/1047.html. #'@seealso #'\link{metab.mle}, \link{metab.bookkeep}, \link{metab.kalman} #'@author Ryan Batt, Luke A. Winslow #' #'@importFrom stats sd density median #' #'@examples #'\dontrun{ #'library(rLakeAnalyzer) #' #'doobs = load.ts(system.file('extdata', #' 'sparkling.doobs', package="LakeMetabolizer")) #'wtr = load.ts(system.file('extdata', #' 'sparkling.wtr', package="LakeMetabolizer")) #'wnd = load.ts(system.file('extdata', #' 'sparkling.wnd', package="LakeMetabolizer")) #'irr = load.ts(system.file('extdata', #' 'sparkling.par', package="LakeMetabolizer")) #' #'#Subset a day #'mod.date = as.POSIXct('2009-07-08', 'GMT') #'doobs = doobs[trunc(doobs$datetime, 'day') == mod.date, ]
#'wtr = wtr[trunc(wtr$datetime, 'day') == mod.date, ] #'wnd = wnd[trunc(wnd$datetime, 'day') == mod.date, ]
#'irr = irr[trunc(irr$datetime, 'day') == mod.date, ] #' #'k600 = k.cole.base(wnd[,2]) #'k.gas = k600.2.kGAS.base(k600, wtr[,3], 'O2') #'do.sat = o2.at.sat(wtr[,1:2], altitude=300) #' #'metab.bayesian(irr=irr[,2], z.mix=rep(1, length(k.gas)), #' do.sat=do.sat[,2], wtr=wtr[,2], #' k.gas=k.gas, do.obs=doobs[,2]) #'} #'@export metab.bayesian <- function(do.obs, do.sat, k.gas, z.mix, irr, wtr, priors, ...){ complete.inputs(do.obs=do.obs, do.sat=do.sat, k.gas=k.gas, z.mix=z.mix, irr=irr, wtr=wtr, error=TRUE) if(any(z.mix <= 0)){ stop("z.mix must be greater than zero.") } if(any(wtr <= 0)){ stop("all wtr must be positive.") } mb.args <- list(...) nobs <- length(do.obs) # ========================================= # = Check for datetime and claculate freq = # ========================================= if("datetime"%in%names(mb.args)){ # check to see if datetime is in the ... args datetime <- mb.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
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
}

# ======================================
# = # Check for priors supplied in ... =
# ======================================
if(!missing(priors)){
t.priors <- priors
p.name.logic <- all(c(c("gppMu", "gppSig2", "rMu", "rSig2", "kSig2"))%in%names(t.priors))
p.class.logic <- is.integer(t.priors) | is.numeric(t.priors)
if(p.name.logic & p.class.logic){
priors <- t.priors
}else{
stop("supplied priors was not a (properly) named numeric/integer vector")
}
}else{
priors <- c("gppMu"=0, "gppSig2"=1E5, "rMu"=0, "rSig2"=1E5, "kSig2"=NA)
}

if(!all(c(is.numeric(do.obs), is.numeric(do.sat), is.numeric(k.gas), is.numeric(z.mix), is.numeric(irr), is.numeric(wtr)))){
stop('All inputs to metab.bayes must be numeric vectors')
}

requireNamespace("R2jags")

# Define model and write to file
# Model choice depends on k values (all 0, all non-0, mixture)
modfile <- bayes.makeModel(k.gas=(k.gas/freq))

# ===========================================
# = Define objects to be used in jags model =
# ===========================================
#Supply elements of U (PAR, log(temp))
U <- matrix(NA, nrow=length(irr), ncol=2)
U[,1] <- irr # PAR Values
U[,2] <- log(wtr) # log(temp) values

# Priors (kP) for k.gas
kP <- matrix(NA, nrow=length(k.gas), ncol=2)
# variances for K
kP[,1] <- k.gas/freq # means for K, adjusted to sampling frequency (velocity denominator in time step, not day)
if(is.na(priors["kSig2"])){
k0.logic <- !is.finite(1/kP[,1]) # test for when k is 0
kP[,2] <- sum(kP[,1])/sum(!k0.logic)*0.1 # k variance = mean of the non-zero K, times 0.1
kP[k0.logic,2] <- 1E-9
}else{
kP[,2] <- priors["kSig2"]
}

# Priors for regression coefficients (cP)
cP <- matrix(NA, nrow=2, ncol=2)
cP[1,1] <- priors["gppMu"] # prior mean of GPP coefficient (C[1,1]*PAR=GPP)
cP[1,2] <- priors["gppSig2"] # prior variance of GPP coefficient
cP[2,1] <- priors["rMu"] # prior mean of R coefficient (C[2,1]*log(Temp)=R)
cP[2,2] <- priors["rSig2"] # prior variance of R coefficient

# Put in final format supplied to jags
data <- list(Y=do.obs, N=length(do.obs), U=U, kP=kP, cP=cP, satO=do.sat, a0=do.obs, Zmix=z.mix)
params <- c("C", "sigmaV", "sigmaW")

output <- bayesFit(data, params, mf=modfile)
return(output)

}


## Try the LakeMetabolizer package in your browser

Any scripts or data that you put into this service are public.

LakeMetabolizer documentation built on Nov. 16, 2022, 1:09 a.m.
{"ft_lang_label":"__label__en","ft_lang_prob":0.51526403,"math_prob":0.8711155,"size":13878,"snap":"2023-14-2023-23","text_gpt3_token_len":4656,"char_repetition_ratio":0.12584691,"word_repetition_ratio":0.20283763,"special_character_ratio":0.39061823,"punctuation_ratio":0.19123784,"nsfw_num_words":0,"has_unicode_error":false,"math_prob_llama3":0.9935035,"pos_list":[0],"im_url_duplicate_count":[null],"WARC_HEADER":"{\"WARC-Type\":\"response\",\"WARC-Date\":\"2023-03-31T09:37:37Z\",\"WARC-Record-ID\":\"<urn:uuid:e778ce90-3a3c-4fe4-bb12-c25d3c23b7d9>\",\"Content-Length\":\"97895\",\"Content-Type\":\"application/http; msgtype=response\",\"WARC-Warcinfo-ID\":\"<urn:uuid:fdb93174-af32-4d93-b00f-030deefbac49>\",\"WARC-Concurrent-To\":\"<urn:uuid:f477d6ac-3914-47c9-81a5-55c46e1e0f0b>\",\"WARC-IP-Address\":\"51.81.83.12\",\"WARC-Target-URI\":\"https://rdrr.io/cran/LakeMetabolizer/src/R/metab.bayesian.R\",\"WARC-Payload-Digest\":\"sha1:4XMQKYFHYZZRPMYU6SJICI7J3HFUL3FW\",\"WARC-Block-Digest\":\"sha1:EZ7OB6OOVA3KRBXXRCH4O4RAFDQ4FQLY\",\"WARC-Identified-Payload-Type\":\"text/html\",\"warc_filename\":\"/cc_download/warc_2023/CC-MAIN-2023-14/CC-MAIN-2023-14_segments_1679296949598.87_warc_CC-MAIN-20230331082653-20230331112653-00263.warc.gz\"}"}