Nothing
#' Make Monthly Prediction Intervals
#'
#' Month statistics using WRTDSKalman bootstrapping approach. The input to this
#' function is the dailyBootOut matrix which contains nReplicate sets of
#' daily flux values for the period of interest.
#' The results are in the form of quantiles of concentration and of flux
#' for each of these months.
#'
#' @param dailyBootOut data frame returned from \code{\link{genDailyBoot}}
#' @param fluxUnit number representing entry in pre-defined fluxUnit class array.
#' \code{\link[EGRET]{printFluxUnitCheatSheet}}
#' @param eList named list with at least the Daily, Sample, and INFO
#' dataframes. Created from the EGRET package, after running \code{\link[EGRET]{modelEstimation}}.
#' @importFrom stats aggregate
#' @return a list of 2 data frames, one for average concentration, in mg/L
#' and one for flux (unit depends on fluxUnit argument)
#' In each data frame the first column is monthSeq that corresponds to the months
#' in the "MonthSeq" column in the eList$Daily data frame. The remaining columns are
#' quantiles of the flux or concentration (depending on the data frame).
#' @export
#' @examples
#' eList <- EGRET::Choptank_eList
#' # This example is only based on 4 iterations
#' # Actual prediction intervals should be calculated on
#' # a much larger number of iterations (several hundred).
#' dailyBoot <- Choptank_dailyBootOut
#' monthPcts <- makeMonthPI(dailyBoot, eList)
#' head(monthPcts[["flux"]])
#' head(monthPcts[["conc"]])
#'
makeMonthPI <- function(dailyBootOut, eList, fluxUnit = 3){
nDaily <- length(dailyBootOut[,1])
nIter <- length(dailyBootOut[1,])
firstMonth <- eList$Daily$MonthSeq[1]
lastMonth <- eList$Daily$MonthSeq[nDaily]
nMonths <- lastMonth - firstMonth + 1
monthsConc <- matrix(data = NA, nrow = nMonths, ncol = nIter)
attr(monthsConc, "firstMonth") <- firstMonth
monthsFlux <- monthsConc
attr(monthsFlux, "firstMonth") <- firstMonth
if (is.numeric(fluxUnit)){
fluxUnit <- EGRET::fluxConst[shortCode=fluxUnit][[1]]
} else if (is.character(fluxUnit)){
fluxUnit <- EGRET::fluxConst[fluxUnit][[1]]
}
for(iter in 1:nIter){
monthsFlux[,iter] <- aggregate(dailyBootOut[,iter], by = list(eList$Daily$MonthSeq), mean)$x * fluxUnit@unitFactor
aveFlow <- aggregate(eList$Daily$Q, by = list(eList$Daily$MonthSeq), mean)$x
monthsConc[,iter] <- monthsFlux[,iter] / (86.4 * aveFlow)
}
DecYear <- aggregate(eList$Daily$DecYear, by = list(eList$Daily$MonthSeq), mean)$x
# now do the percentiles by month
probs <- c(0.01, 0.025, 0.05, 0.1, 0.25, 0.5, 0.75, 0.90, 0.95, 0.975, 0.99)
monthFluxPcts <- matrix(data = NA, nrow = nMonths, ncol = 11)
monthConcPcts <- matrix(data = NA, nrow = nMonths, ncol = 11)
for(iMonth in 1:nMonths){
monthFluxPcts[iMonth,] <- quantile(monthsFlux[iMonth,], probs = probs, type = 6, na.rm = TRUE)
monthConcPcts[iMonth,] <- quantile(monthsConc[iMonth,], probs = probs, type = 6, na.rm = TRUE)
}
lastMonth <- firstMonth + length(monthFluxPcts[,1]) - 1
monthSeqVec <- seq(firstMonth, lastMonth)
month_flux_results <- data.frame(monthSeq = monthSeqVec,
DecYear = DecYear,
p1 = monthFluxPcts[,1],
p2.5 = monthFluxPcts[,2],
p5 = monthFluxPcts[,3],
p10 = monthFluxPcts[,4],
p25 = monthFluxPcts[,5],
p50 = monthFluxPcts[,6],
p75 = monthFluxPcts[,7],
p90 = monthFluxPcts[,8],
p95 = monthFluxPcts[,9],
p97.5 = monthFluxPcts[,10],
p99 = monthFluxPcts[,11])
month_conc_results <- data.frame(monthSeq = monthSeqVec,
DecYear = DecYear,
p1 = monthConcPcts[,1],
p2.5 = monthConcPcts[,2],
p5 = monthConcPcts[,3],
p10 = monthConcPcts[,4],
p25 = monthConcPcts[,5],
p50 = monthConcPcts[,6],
p75 = monthConcPcts[,7],
p90 = monthConcPcts[,8],
p95 = monthConcPcts[,9],
p97.5 = monthConcPcts[,10],
p99 = monthConcPcts[,11])
return(list(flux = month_flux_results,
conc = month_conc_results))
}
#' monthSeqToDec
#'
#' Convert a sequence of month integers into their decimal years.
#'
#' @export
#' @param monthSeq integer vector of months. Month 1 is considered Jan. 1850.
#' @examples
#'
#' months <- 1558:1600
#' monthSeqToDec(months)
monthSeqToDec <- function(monthSeq){
# function that turns a monthSeq value into a decYear value
dec <- floor((monthSeq - 1)/12) +
1850 + (1/24) +
((monthSeq-1)%%12/12)
return(dec)
}
#' Make Annual Prediction Intervals
#'
#' This function takes the output from \code{\link{genDailyBoot}} and
#' calculates the quantiles for an annual (based on paStart/paLong) aggregation.
#' This means that the function can be used for seasons.
#'
#' @param dailyBootOut data frame returned from \code{\link{genDailyBoot}}
#' @param fluxUnit number representing entry in pre-defined fluxUnit class array.
#' \code{\link[EGRET]{printFluxUnitCheatSheet}}
#' @param eList named list with at least the Daily, Sample, and INFO
#' dataframes. Created from the EGRET package, after running \code{\link[EGRET]{modelEstimation}}.
#' @param paLong numeric integer specifying the length of the period of analysis, in months, 1<=paLong<=12, default is 12
#' @param paStart numeric integer specifying the starting month for the period of analysis, 1<=paStart<=12, default is 10
#' @importFrom stats aggregate
#'
#' @return a list of 2 data frames, one for average concentration, in mg/L
#' and one for flux (unit depends on fluxUnit argument)
#' In each data frame the first column is DecYear. The remaining columns are
#' quantiles of the flux or concentration (depending on the data frame).
#' @export
#' @examples
#' eList <- EGRET::Choptank_eList
#' # This example is only based on 4 iterations
#' # Actual prediction intervals should be calculated on
#' # a much larger number of iterations (several hundred).
#' dailyBoot <- Choptank_dailyBootOut
#' annualPcts <- makeAnnualPI(dailyBoot, eList)
#' head(annualPcts[["flux"]])
#' head(annualPcts[["conc"]])
#'
makeAnnualPI <- function(dailyBootOut, eList,
paLong = 12, paStart = 10, fluxUnit = 3){
nIter <- length(dailyBootOut[1,])
localDaily <- eList$Daily
Annual_results <- data.frame()
if (is.numeric(fluxUnit)){
fluxUnit <- EGRET::fluxConst[shortCode=fluxUnit][[1]]
} else if (is.character(fluxUnit)){
fluxUnit <- EGRET::fluxConst[fluxUnit][[1]]
}
Annual_results <- data.frame()
for(i in 1:nIter){
localDaily$GenFlux <- dailyBootOut[,i]
localDaily$GenConc <- localDaily$GenFlux /(localDaily$Q *86.40)
AnnualResults_i <- EGRET::setupYears(localDaily,
paLong = paLong, paStart = paStart)
AnnualResults_i <- AnnualResults_i[,c("DecYear", "GenFlux", "GenConc")]
Annual_results <- rbind(Annual_results, AnnualResults_i)
}
probs <- c(0.01, 0.025, 0.05, 0.1, 0.25, 0.5, 0.75, 0.90, 0.95, 0.975, 0.99)
years <- unique(Annual_results$DecYear)
yearFluxPcts <- matrix(data = NA, nrow = length(years), ncol = length(probs))
yearConcPcts <- matrix(data = NA, nrow = length(years), ncol = length(probs))
for(year in years){
i <- which(years %in% year)
flux_conc_year <- Annual_results[Annual_results$DecYear == year,]
yearFluxPcts[i,] <- quantile(flux_conc_year$GenFlux, probs = probs, type = 6, na.rm = TRUE) * fluxUnit@unitFactor
yearConcPcts[i,] <- quantile(flux_conc_year$GenConc, probs = probs, type = 6, na.rm = TRUE)
}
year_flux_results <- data.frame(DecYear = years,
p1 = yearFluxPcts[,1],
p2.5 = yearFluxPcts[,2],
p5 = yearFluxPcts[,3],
p10 = yearFluxPcts[,4],
p25 = yearFluxPcts[,5],
p50 = yearFluxPcts[,6],
p75 = yearFluxPcts[,7],
p90 = yearFluxPcts[,8],
p95 = yearFluxPcts[,9],
p97.5 = yearFluxPcts[,10],
p99 = yearFluxPcts[,11])
year_conc_results <- data.frame(DecYear = years,
p1 = yearConcPcts[,1],
p2.5 = yearConcPcts[,2],
p5 = yearConcPcts[,3],
p10 = yearConcPcts[,4],
p25 = yearConcPcts[,5],
p50 = yearConcPcts[,6],
p75 = yearConcPcts[,7],
p90 = yearConcPcts[,8],
p95 = yearConcPcts[,9],
p97.5 = yearConcPcts[,10],
p99 = yearConcPcts[,11])
return(list(flux = year_flux_results,
conc = year_conc_results))
}
#' Make Daily Prediction Intervals
#'
#' This function takes the output from \code{\link{genDailyBoot}} and
#' calculates the quantiles for a daily aggregation.
#'
#' @param dailyBootOut data frame returned from \code{\link{genDailyBoot}}
#' @param fluxUnit number representing entry in pre-defined fluxUnit class array.
#' \code{\link[EGRET]{printFluxUnitCheatSheet}}
#' @param eList named list with at least the Daily, Sample, and INFO
#' dataframes. Created from the EGRET package, after running \code{\link[EGRET]{modelEstimation}}.
#' @importFrom stats aggregate
#'
#' @return a list of 2 data frames, one for average concentration, in mg/L
#' and one for flux (unit depends on fluxUnit argument)
#' In each data frame the first column is Date. The remaining columns are
#' quantiles of the flux or concentration (depending on the data frame).
#' @export
#' @examples
#' eList <- EGRET::Choptank_eList
#' # This example is only based on 4 iterations
#' # Actual prediction intervals should be calculated on
#' # a much larger number of iterations (several hundred).
#' dailyBoot <- Choptank_dailyBootOut
#' \donttest{
#' dailyPcts <- makeDailyPI(dailyBoot, eList)
#' head(dailyPcts[["flux"]])
#' head(dailyPcts[["conc"]])
#' }
makeDailyPI <- function(dailyBootOut, eList, fluxUnit = 3){
nDaily <- length(dailyBootOut[,1])
nIter <- length(dailyBootOut[1,])
if (is.numeric(fluxUnit)){
fluxUnit <- EGRET::fluxConst[shortCode=fluxUnit][[1]]
} else if (is.character(fluxUnit)){
fluxUnit <- EGRET::fluxConst[fluxUnit][[1]]
}
# now do the percentiles by month
probs <- c(0.01, 0.025, 0.05, 0.1, 0.25, 0.5, 0.75, 0.90, 0.95, 0.975, 0.99)
dailyFluxPcts <- matrix(data = NA, nrow = nDaily, ncol = 11)
dailyConcPcts <- matrix(data = NA, nrow = nDaily, ncol = 11)
for(iDay in 1:nDaily){
dailyFluxPcts[iDay,] <- quantile(dailyBootOut[iDay,], probs = probs, type = 6, na.rm = TRUE)
dailyConcPcts[iDay,] <- quantile(dailyBootOut[iDay,]/ (86.4 * eList$Daily$Q[iDay]), probs = probs, type = 6, na.rm = TRUE)
}
daily_flux_results <- data.frame(Date = eList$Daily$Date,
DecYear = eList$Daily$DecYear,
p1 = dailyFluxPcts[,1],
p2.5 = dailyFluxPcts[,2],
p5 = dailyFluxPcts[,3],
p10 = dailyFluxPcts[,4],
p25 = dailyFluxPcts[,5],
p50 = dailyFluxPcts[,6],
p75 = dailyFluxPcts[,7],
p90 = dailyFluxPcts[,8],
p95 = dailyFluxPcts[,9],
p97.5 = dailyFluxPcts[,10],
p99 = dailyFluxPcts[,11])
daily_conc_results <- data.frame(Date = eList$Daily$Date,
DecYear = eList$Daily$DecYear,
p1 = dailyConcPcts[,1],
p2.5 = dailyConcPcts[,2],
p5 = dailyConcPcts[,3],
p10 = dailyConcPcts[,4],
p25 = dailyConcPcts[,5],
p50 = dailyConcPcts[,6],
p75 = dailyConcPcts[,7],
p90 = dailyConcPcts[,8],
p95 = dailyConcPcts[,9],
p97.5 = dailyConcPcts[,10],
p99 = dailyConcPcts[,11])
return(list(flux = daily_flux_results,
conc = daily_conc_results))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.