# loadReg2
# We initially extended the current loadReg class to implement
# loadModelInterface, but that implementation was risky - getOriginalData could
# work correctly one time and not the next, for example.
#### loadReg2 class ####
#' A load model class specific to loadReg objects produced by the USGS rloadest
#' package.
#'
#' This class is closely related to loadModel except that the model slot is more
#' strictly defined, and the associated methods are a little different.
#'
#' @include loadModel.R
#' @include loadReg.R
#' @rdname loadReg2-class
#' @name loadReg2-class
#' @importFrom methods setClass
#' @exportClass loadReg2
#' @family load.model.classes
setClass(
"loadReg2",
# S3 objects can't be listed as types here, so we'll enforce loadReg types in validity() below.
# slots=c(
# #fit="loadReg"),
contains="loadModel",
# Validity of parent objects is checked first, automatically
validity=function(object) {
errorstrs <- character()
if(!is(object@fit, "loadReg")) {
errorstrs <- c(errorstrs, "object@fit should be a loadReg object")
}
md <- object@metadata
if(md@load.rate.units != paste0(md@load.units, " d^-1")) {
errorstrs <- c(errorstrs, "for loadReg2 models, object@metadata@load.rate.units must equal object@metadata@load.rate / day")
}
# Return
if(length(errorstrs) == 0) {
TRUE
} else {
errorstrs
}
}
)
#' Create a fitted loadReg2 object.
#'
#' Generates a new model of class loadReg2 (\code{\link{loadReg2-class}}).
#' loadReg2s are wrappers for loadReg objects produced by the USGS
#' \pkg{rloadest} package. \code{loadReg2}s can implement the
#' \code{\link{loadModelInterface}} more reliably than is possible for a
#' \code{loadReg} object. Some metadata information is set within the loadReg()
#' call; entries that can be edited may be given as arguments to the loadReg2()
#' call (e.g., site.id, lat, lon, basin.area).
#'
#' @importFrom methods is new
#' @param load.reg An unevaluated call to \code{\link[rloadest]{loadReg}}. This
#' call will be parsed and evaluated within \code{loadReg2} to create a fully
#' functional load model for use within \pkg{loadflex}.
#' @param pred.format character. Should predictions be made for 'flux' (load
#' rate) or 'conc' (concentration)? \pkg{rloadest}, and therefore loadReg2,
#' uses different models for flux and concentration, though fitted to the same
#' data and with the same model structure except for whether the left-hand
#' side of the model formula is load rate or concentration. The model specific
#' to \code{pred.format} will be used to generate predictions.
#' @param store One or more character strings specifying which information to
#' write within the model. Options are 'data': the original fitting data;
#' 'fitting.function': a fitting function that can produce a new loadComp
#' object from new data
#' @inheritParams metadata
#' @param ... Other arguments passed to this model.
#' @return A fitted loadReg2 model.
#'
#' @importFrom rloadest loadReg
#' @export
#' @family load.model.inits
#' @examples
#' \dontrun{
#' library(rloadest)
#' no3_lr <- suppressWarnings(
#' loadReg2(loadReg(NO3 ~ model(9), data=get(data(lamprey_nitrate)),
#' flow="DISCHARGE", dates="DATE", time.step="instantaneous",
#' flow.units="cfs", conc.units="mg/L", load.units="kg",
#' station='Lamprey River, NH')))
#' }
loadReg2 <- function(
load.reg, pred.format=c("flux","conc"), store=c("data","fitting.function"),
consti.name="", load.rate="",
site.id="", lat=as.numeric(NA), lon=as.numeric(NA), basin.area=as.numeric(NA),
flow.site.name="", flow.site.id="", flow.lat=as.numeric(NA), flow.lon=as.numeric(NA), flow.basin.area=as.numeric(NA),
basin.area.units="km^2", custom=NULL,
...) {
# Validate arguments
pred.format <- match.arg.loadflex(pred.format)
store <- match.arg.loadflex(store, choices=c("data","fitting.function"), several.ok=TRUE)
# Require that load.reg is a call to loadReg (and not an already-generated
# loadReg) to ensure that we can extract the right information from it.
loadReg_call <- substitute(load.reg)
if(is(loadReg_call, "name")) {
# The refitting function [only] will pass the call as a symbol (name)
# containing language that is itself a call to loadReg. A clever user could
# probably trick this test into thinking they're inputting a refit...that's
# on them to avoid.
loadReg_call <- load.reg
is_refit <- TRUE
} else {
# If loadReg_call isn't a symbol, then we're fitting this model for the
# first time (at least, it's not an official refitting)
is_refit <- FALSE
}
# By now, loadReg_call should always be a call to loadReg(). Check that it is.
call_ok <- FALSE
if(is(loadReg_call, "call")) {
if(loadReg_call[[1]] == "loadReg") {
call_ok <- TRUE
# one option for the future is to replace the function name with the
# namespace::function version, as in "loadReg_call[[1]] <-
# quote(rloadest::loadReg)". someday this may be sufficient to evaluate an
# rloadest model, but we're not there yet. For one thing, loadReg calls
# may include additional loadReg functions such as center(). For another,
# rloadest:::loadestTadj appears to call dectime without declaring imports
# smwrBase, which means you need to have rloadest actually loaded so that
# its dependencies (including smwrBase) are also loaded.
}
}
if(!call_ok) {
stop("load.reg must be a call of the form loadReg(...)")
}
# Require that the rloadest namespace has been loaded by the user. Can't do
# the package loading here because we would get the WARNING "'library' or
# 'require' call to 'rloadest' in package code." during R CMD check. Can't
# import the relevant rloadest functions into this function because they're
# not openly named anywhere in this function code, which means R CMD check
# would worry that we're importing things without using them. So instead we'll
# require the user to do the loading for us.
checkRloadestStatus()
# Fit or refit the model
if(is_refit) {
# Refitting occurs in the context (environment) of refit.envir, which is an
# argument that we've left off the list of documented arguments (it's in
# ...) because normal users shouldn't be aware of it or using it explicitly.
if(!("refit.envir" %in% names(list(...)))) {
stop("Expecting a refit.envir when load.reg is a symbol")
}
refit_envir <- list(...)[["refit.envir"]]
# Now refit in the refit_envir. enclos works because refit_envir is actually a list.
load.reg <- eval(loadReg_call, envir=refit_envir, enclos=parent.frame())
} else {
# if first time, fit in parent environment as a straight loadReg call would
load.reg <- eval(loadReg_call, envir=parent.frame())
}
# Get the metadata, which we'll need to do the following data checks
meta <- getMetadata(load.reg)
# Fill the metadata fields not fixed by the rloadest loadReg object
meta <- updateMetadata(
meta,
consti.name=consti.name, load.rate=load.rate,
site.id=site.id, lat=lat, lon=lon, basin.area=basin.area,
flow.site.name=flow.site.name, flow.site.id=flow.site.id, flow.lat=flow.lat, flow.lon=flow.lon, flow.basin.area=flow.basin.area,
basin.area.units=basin.area.units, custom=custom)
# Parse the loadReg call, attaching argument names as needed. This will be
# useful both for the refitting function and for saving the original data.
if("fitting.function" %in% store | "data" %in% store) {
loadReg_call <- match.call(definition=loadReg, call=loadReg_call)
}
# Pull the data into local memory for error checking. Look for easy problems -
# namely, 0's or negative numbers in the discharge or concentration
# measurements, which will have been logged without noting problems in
# loadReg().
load.model.data <- eval(loadReg_call[["data"]], envir=parent.frame())
if(any(getCol(metadata=meta, data=load.model.data, "flow", attach.units=FALSE) <= 0, na.rm=TRUE)) {
warning("Some discharge values are <= 0; this can cause flawed loadReg model fits without further warning")
}
if(any(getCol(metadata=meta, data=load.model.data, "conc", attach.units=FALSE) <= 0, na.rm=TRUE)) {
warning("Some concentration values are <= 0; this can cause flawed loadReg model fit without further warning")
}
## REFITTING FUNCTION ##
if("fitting.function" %in% store) {
# Create an environment (or identify the existing one, if is_refit) where
# each argument to loadReg() is saved as an object.
if(!is_refit) {
# Start creating the loadReg piece of the refitting function. We will
# eventually keep only the reference to loadReg and the formula, so
# everything else is simply a placeholder.
loadReg_recall <- loadReg_call
# In the refitting function, the loadReg() call should refer to
# training.data instead of the original fitting data.
loadReg_recall[["data"]] <- as.symbol("training.data")
# All other arguments in loadReg_recall should be the same as in the
# original call. But if those arguments contain any symbols, then the
# value of those symbols might change, and the behavior of
# refitting_function could be quite unpredictable. To avoid this, we will
# copy the value of each argument into a new, safe environment, and we
# will change the references in loadReg_recall to refer to those saved
# copies instead of the originals.
# Here's the new, safe environment.
refit_envir <- new.env(parent=parent.frame(n=2))
# Change all arguments except "formula" and "data" to point into refit_envir
save_these <- !(names(loadReg_call) %in% c("", "formula", "data"))
for(argname in names(loadReg_call)[save_these]) {
savename <- paste0("saved.", argname)
# Copy the argument value into refit_envir under a new name (savename)
assign(savename, eval(loadReg_call[[argname]], envir=parent.frame()), envir=refit_envir)
# Change the calling function to refer to this saved copy
loadReg_recall[[argname]] <- as.symbol(savename)
}
# For all that talk about environments, it's actually easier for
# refit_envir to be a list; this will permit us to use eval() with an
# enclos argument that actually gets used instead of ignored.
refit_envir <- as.list(refit_envir)
} else {
# If the refit.envir argument was supplied, then (1) we have already
# loaded refit_envir above, and (2) we will assume loadReg_call already
# contains references to objects within refit_envir - i.e., is already the
# call we want to use for refitting.
loadReg_recall <- loadReg_call
}
# The full refitting function wraps loadReg_recall within a loadReg2 call.
# refitting_function is a closure that will remember the current values of
# loadReg_recall, pred.format, store, and refit_envir exactly as they
# are now.
load.model.fitting.function <- function(training.data, store=c()) {
loadReg2(loadReg_recall, pred.format=pred.format, store=store, refit.envir=refit_envir)
}
}
# Create and return the model
load.model <- new(
"loadReg2",
fit=load.reg,
pred.format=pred.format,
metadata=meta,
data=if("data" %in% store) load.model.data else NULL,
fitting.function=if("fitting.function" %in% store) load.model.fitting.function else NULL)
}
#' Make flux or concentration predictions from a loadReg2 model.
#'
#' Makes instantaneous predictions (at the temporal resolution of
#' \code{newdata}) from a fitted \code{\link{loadReg2}} model. See
#' \code{\link{predictSolute}} for details.
#'
#' @inheritParams predictSolute
#' @param load.model A loadReg2 object.
#' @param newdata \code{data.frame}, optional. Predictor data. Column names
#' should match those given in the \code{loadReg2} metadata. If \code{newdata}
#' is not supplied, the original fitting data will be used.
#' @return A vector of data.frame of predictions, as for the generic
#' \code{\link{predictSolute}}.
#' @importFrom rloadest predConc predLoad
#' @export
#' @family predictSolute
predictSolute.loadReg2 <- function(
load.model, flux.or.conc=c("flux","conc"), newdata,
interval=c("none","confidence","prediction"), level=0.95,
lin.or.log=c("linear","log"), se.fit=FALSE, se.pred=FALSE,
date=FALSE, attach.units=FALSE,
agg.by=c("unit", "day", "month", "water year", "calendar year", "total",
"mean water year", "mean calendar year", "[custom]"), ...) {
# Validate arguments
flux.or.conc <- match.arg.loadflex(flux.or.conc)
interval <- match.arg.loadflex(interval)
lin.or.log <- match.arg.loadflex(lin.or.log)
agg.by <- match.arg(agg.by)
# Validate rloadest status
checkRloadestStatus()
# Extract the metadata in object into our standard metadata format, and
# confirm that load.units are as we first passed them to loadReg
metadata <- getMetadata(load.model)
loadReg2.load.units <- sapply(c('load.units','load.rate.units'), getInfo, metadata=getMetadata(load.model))
rloadest.load.units <- sapply(c('load.units','load.rate.units'), getInfo, metadata=getMetadata(getFittedModel(load.model)))
if(!isTRUE(all.equal(loadReg2.load.units, rloadest.load.units))) {
stop("loadReg2 metadata on load.units and load.rate.units must not be modified from the rloadest (loadReg) values. ",
"rloadest always predicts load rate in load.units per day")
}
if(length(list(...)) > 0) {
args <- list(...)
argnames <- names(args)
accepted_argnames <- c("seopt", "print") # "load.units","conf.int"
if(is.null(argnames) || any(!(argnames %in% accepted_argnames))) {
stop("unrecognized args: ", paste0(setdiff(argnames, accepted_argnames), collapse=", "),
". args in ... should be in this list: ", paste0(accepted_argnames, collapse=", "))
}
}
# If newdata==NULL, we should return the predictions from the fitting data.
if(missing(newdata)) {
newdata <- getFittingData(load.model)
}
# If time.step != "instantaneous", input data are aggregated to time.step by
# rloadest::predLoad()/predConc() before being used for prediction. This is
# logical - predictor data should be at the same level of aggregation as the
# data with which the model was calibrated. But predict.loadReg() is supposed
# to act like other predict() methods, and so we'll emit a message if there's
# any aggregation, to make it clearer that something non-standard is
# happening.
if(load.model@fit$time.step != "instantaneous") {
if(load.model@fit$time.step == "day" & class(getCol(metadata, newdata, "date"))[1] == "Date") {
# do nothing in this case because there won't be aggregation
} else {
# but if it's not the day-day case, tell the user what's going to happen
message("Predictor data will be aggregated to time.step='",load.model@fit$time.step,
"' before applying model. To override, set time.step='instantaneous'")
}
# NB: Whenever time.step == "instantaneous", loadReg/loadConc requires that
# dates be in POSIXct format. We'll let those functions deal with any user
# errors on this point.
}
#deal with mean water/calendar year options
#need to filter data down to compelete years, and then send 'total' option to rloadest
#data freq, leap years, whatever is in newdata ob
if(agg.by %in% c("mean water year", "mean calendar year")) {
newdata <- filterCompleteYears(data = newdata, time.step = load.model@fit$time.step,
agg.by = agg.by)
agg.by = "total" #going to rloadest
if(nrow(newdata) > 176000) {
stop(paste(strwrap("Sorry, rloadest can't handle more than 176,000 data points at a time.
Please change the agg.by argument to a shorter period, or else
supply a shorter period of data for newdata"), collapse = "\n"))
}
}
# Get around data size constraints of predLoad and predConc by splitting the
# prediction into smaller chunks as necessary. chunk.size is the maximum
# number of values per chunk; rloadest can handle at most 176000 values by
# default. nchunks is the number of chunks required.
chunk.size <- 175872 #multiple of 96, so breaks iv data on day boundaries
nchunks <- ceiling(nrow(newdata) / chunk.size)
datachunks <- lapply(1:nchunks, function(i) {
newdata[((i-1)*chunk.size + 1):min(i*chunk.size, nrow(newdata)),]
})
# Now do the prediction for each chunk and then reassemble the full set of
# predictions
preds_lin_raw <- lapply(datachunks, function(datachunk) {
switch(
flux.or.conc,
"flux"={
# ... may contain seopt, print
predLoad(fit=load.model@fit, newdata=datachunk, by=agg.by, allow.incomplete=FALSE, conf.int=level, ...)
},
"conc"={
predLoad_args <- intersect(c('seopt','print'), names(list(...))) # load.units and conf.int have already been rejected above
if(length(predLoad_args) > 0) warning(paste("these args are ignored for flux.or.conc='conc':", paste(predLoad_args, collapse=', ')))
#will fail here if agg.by != 'unit' or 'day'
predConc(fit=load.model@fit, newdata=datachunk, by=agg.by, allow.incomplete=FALSE, conf.int=level)
}
)
})
preds_lin_raw <- do.call(rbind, preds_lin_raw)
preds_col <- .sentenceCase(flux.or.conc)
# Generate predictions in log space in case these were requested. Uses the
# assumption, as in predLoad and predConc, that the Ferguson (1986) method is
# adequate for re-retransformation. (predLoad and predConc use these to
# compute prediction intervals)
preds_log_raw <- preds_lin_raw[,1,drop=FALSE]
preds_log_raw$meanlog <- preds_log_raw$sdlog <- NA
KDays <- !is.na(preds_lin_raw[["SEP"]])
preds_log_raw[KDays,c("meanlog","sdlog")] <-
linToLog(meanlin=preds_lin_raw[[preds_col]][KDays], sdlin=preds_lin_raw$SEP[KDays])
# Format predictions; add intervals if requested
if(interval == "none") {
preds <- switch(
lin.or.log,
'linear' = preds_lin_raw[[preds_col]],
'log' = preds_log_raw[["meanlog"]]
)
} else if(interval == "confidence") {
stop("confidence intervals not implemented for loadReg2 models")
} else if(interval=="prediction") {
# From the predLoad/predConc documentation: "The term confidence interval is
# used here as in the original documentation for LOADEST, but the values
# that are reported are the prediction intervals, computed from the SEP."
pred.int.num <- as.character(round(level*100, 0))
preds <- switch(
lin.or.log,
'linear' =
data.frame(
fit = preds_lin_raw[[preds_col]],
lwr = preds_lin_raw[[paste0("L",pred.int.num)]],
upr = preds_lin_raw[[paste0("U",pred.int.num)]]),
'log' =
# preds_lin CI was computed in log space, so a straight log transformation
# is fine here. Use the bias re-corrected mean, though
data.frame(
fit = preds_log_raw$meanlog,
lwr = log(preds_lin_raw[[paste0("L",pred.int.num)]]),
upr = log(preds_lin_raw[[paste0("U",pred.int.num)]]))
)
} # other possibilities prevented by match.arg.loadflex(interval) above
# If se.fit==TRUE or se.pred==TRUE, format the output to approximately
# parallel the output of other predict functions: a list of vectors with
# SE-related information
if(se.fit | se.pred) {
if(!is.data.frame(preds)) {
preds <- data.frame(fit=preds)
}
if(se.fit) {
preds$se.fit <- switch(
lin.or.log,
'linear' = preds_lin_raw[["Std.Err"]],
'log' = NA # wouldn't be compatible with the meanlog computed from the SEP
)
}
if(se.pred) {
preds$se.pred <- switch(
lin.or.log,
'linear' = preds_lin_raw[["SEP"]],
'log' = preds_log_raw[["sdlog"]]
)
}
}
# Attach units if requested. Applies to fit, lwr, upr, se.fit, and se.pred.
# formatPreds with same from and to formats doesn't do anything if
# attach.units==FALSE, so don't bother unless attach.units==TRUE
if(attach.units) {
if(lin.or.log == 'log') {
warning("units requested but unavailable for lin.or.log='log'")
} else {
if(is.data.frame(preds)) {
preds <- as.data.frame(lapply(preds, function(pr) {
formatPreds(
pr, from.format=flux.or.conc, to.format=flux.or.conc,
newdata=NULL, metadata=metadata, attach.units=attach.units)
}))
} else {
preds <- formatPreds(
preds, from.format=flux.or.conc, to.format=flux.or.conc,
newdata=NULL, metadata=metadata, attach.units=attach.units)
}
}
}
# Add dates if requested
if(date) {
if(!is.data.frame(preds)) {
preds <- data.frame(fit=preds)
}
# prepend the date column
preds <- data.frame(date=getCol(metadata, newdata, "date"), preds)
#output column names changes depending on period
date_col_name <- intersect(names(preds_lin_raw), c("Period", "Date"))
preds <- data.frame(date=preds_lin_raw[date_col_name], preds)
}
# Bring out units if hidden in a data.frame
if(attach.units && is.data.frame(preds)) {
preds <- u(preds)
}
return(preds)
}
#' Produce a set of predictions that reflect the coefficient uncertainty and
#' possibly also natural variation.
#'
#' This function resamples the coefficients from their joint distribution, then
#' makes predictions whose individual errors are sampled from a time series with
#' the same first-order autocorrelation as the original series of errors.
#'
#' @importFrom stats arima.sim
#' @inheritParams simulateSolute
#' @param load.model A loadReg2 object.
#' @param newdata \code{data.frame}, optional. Predictor data. Column names
#' should match those given in the \code{loadReg2} metadata. If
#' \code{newdata} is not supplied, the original fitting data will be used.
#' @return A vector of predictions that are distributed according to the
#' uncertainty of the coefficients and the estimated natural variability +
#' measurement error.
#' @export
#' @family simulateSolute
simulateSolute.loadReg2 <- function(load.model, flux.or.conc=c("flux","conc"), newdata,
method=c("parametric", "non-parametric"), from.interval=c("confidence", "prediction"), rho, ...) {
# Validate arguments
flux.or.conc <- match.arg.loadflex(flux.or.conc)
from.interval <- match.arg.loadflex(from.interval, c("confidence","prediction"))
method <- match.arg.loadflex(method, c("parametric", "non-parametric"))
if(missing(newdata)) newdata <- getFittingData(load.model) # need to do this for predictSolute call in case non-parametric method overwrites load.model@data
checkRloadestStatus()
# Generate bootstrap model
if(method=="parametric") {
# Resample the coefficients and put them into the load.model@fit
load.model@fit <- resampleCoefficients.loadReg(load.model@fit, flux.or.conc)
} else if(method=="non-parametric") {
# Refit the model to new (resampled) data - this probably takes longer but has nice statistical properties
df <- getFittingData(load.model)
load.model <- getFittingFunction(load.model)(
training.data=df[sample.int(nrow(df), replace=TRUE), ],
store=c())
}
# Make predictions for the new data
fitting.preds <- predictSolute(load.model, flux.or.conc=flux.or.conc, newdata=newdata, interval="none")
## Resample Errors ##
# EXPERIMENTAL: To keep the autocorrelation that's necessary for loadComp to
# make sense (even though it's assumed not to exist in the calibration data),
# either (1) estimate autocorrelation from the calibration data, or (2) accept
# an estimate of autocorrelation from the user. Then use that autocorrelation
# and the residual standard deviation to add noise to the predictions with an
# autocorrelation structure based on that of the residuals (or supplied by the
# caller)
if(from.interval == "prediction") {
if(missing(rho)) {
# NOT FINISHED: The abs.or.rel and use.log arguments here are specific to most but not all models - make more robust
rho.out <- estimateRho(load.model, flux.or.conc=flux.or.conc, abs.or.rel.resids="absolute", use.log=TRUE,
newdata=newdata, irregular.timesteps.ok=TRUE, plot.acf=FALSE)
rho <- rho.out$rho # implies that the time.step for newata is the one we want
}
# Simulate the residuals so that they have the autocorrelation structure
# that was estimated from the original residuals or supplied by the caller.
# Not sure whether rho takes care of matching the distribution of resampled
# residuals to the original; it seems not to.
noise <- as.numeric(arima.sim(model=list(order=c(1,0,0), ar=rho), n=nrow(newdata)))
# Normalize the noise to have the same sd as the original residuals. Is that reasonable?
error <- "print.function"
error("need to calculate s.hat")
#noise <- noise * s.hat/sd(noise)
# Add the noise back into the predictions, first converting back to linear space
fitting.preds <- exp(log(fitting.preds) + noise)
}
return(fitting.preds)
}
#' Extract model summary statistics from a loadReg2 model
#'
#' Produce a 1-row data.frame of model metrics. The relevant metrics for
#' loadReg2 models are largely the same as those reported by the `rloadest`
#' package, though reported in this streamlined data.frame format for bulk
#' reporting. `summarizeModel.loadReg2` should rarely be accessed directly;
#' instead, call `summarizeModel()` on a [loadReg2] object.
#'
#' @md
#' @inheritParams summarizeModel
#' @return Returns a 1-row data frame with the following columns:
#'
#' * `site.id` - the unique identifier of the site, as in [metadata()]
#'
#' * `constituent` - the unique identifier of the constituent, as in
#' [metadata()]
#'
#' * `eqn` - the regression equation, possibly in the form `const ~ model(x)`
#' where `x` is the `Number` of a pre-defined equation in [rloadest::Models]
#'
#' * `RMSE` - the square root of the mean squared error. Errors will be
#' computed from either fluxes or concentrations, as determined by the value
#' of `pred.format` that was passed to [loadReg2()] when this model was
#' created
#'
#' * `r.squared` - the r-squared value, generalized for censored data,
#' describing the amount of observed variation explained by the model
#'
#' * `p.value` - the p-value for the overall model fit
#'
#' * `cor.resid` - the serial correlation of the model residuals
#'
#' * `PPCC` - the probability plot correlation coefficient measuring the
#' normality of the residuals
#'
#' * `Intercept`, `lnQ`, `lnQ2`, `DECTIME`, `DECTIME2`, `sin.DECTIME`,
#' `cos.DECTIME`, etc. - the fitted value of the intercept and other terms
#' included in this model (list differs by model equation)
#'
#' * `Intercept.SE`, `lnQ.SE`, `lnQ2.SE`, `DECTIME.SE`, `DECTIME2.SE`,
#' `sin.DECTIME.SE`, `cos.DECTIME.SE`, etc. - the standard error of the fitted
#' estimates of these terms
#'
#' * `Intercept.p.value`, `lnQ.p.value`, `lnQ2.p.value`, `DECTIME.p.value`,
#' `DECTIME2.p.value`, `sin.DECTIME.p.value`, `cos.DECTIME.p.value` - the
#' p-values for each of these model terms
#'
#' @importFrom dplyr select everything
#' @export
#' @family summarizeModel
#' @examples
#' \dontrun{
#' library(rloadest)
#' no3_lr <- suppressWarnings(
#' loadReg2(loadReg(NO3 ~ model(9), data=get(data(lamprey_nitrate)),
#' flow="DISCHARGE", dates="DATE", time.step="instantaneous",
#' flow.units="cfs", conc.units="mg/L", load.units="kg",
#' station='Lamprey River, NH')))
#' summarizeModel(no3_lr)
#' }
summarizeModel.loadReg2 <- function(load.model, ...) {
# collect and combine the default summary and the loadReg summary
out.info <- NextMethod()
out.loadReg <- summarizeModel(getFittedModel(load.model), flux.or.conc=load.model@pred.format, ...)
out <- bind_cols(out.info, out.loadReg)
# return
return(out)
}
#remove incomplete water/calendar years
filterCompleteYears <- function(data, time.step, agg.by) {
if(agg.by == "mean water year") {
data$year <- smwrBase::waterYear(data$DATE, numeric = TRUE)
} else {
data$year <- lubridate::year(data$DATE)
}
daySteps <- ifelse(time.step == "instantaneous", yes = 96, no = 1)
data_summary <- dplyr::group_by(data, year) %>%
dplyr::summarise(n = n()) %>%
dplyr::mutate(leap = lubridate::leap_year(year),
completeThresh = daySteps * 365 + leap * daySteps)
completeThresh <- '.dplyrvar'
year <- '.dplyrvar'
leap <- '.dplyrvar'
completeYears <- filter(data_summary, n >= completeThresh)
returnData <- filter(data, year %in% completeYears$year) %>% select(-year)
if(nrow(returnData) == 0) {
stop("No complete periods. Please select a different aggregation period")
}
return(returnData)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.