#' Rating-curve generalized additive model
#'
#' Fits generalized additive models to log-log relationships using `mgcv::gam`,
#' but incorporates aspects
#' specific to hydrologic rating curves. Specifically, tracks the geometric mean of
#' flow and concentration, adjusts for transformation bias, and calculates relevant
#' fit statistics such as Nash-Sutcliffe efficiency (NSE)
#'
#' Dots get passed to mgcv::gam() call.
#' All logarithms are base e (natural log)
#' Arguments ending in "var" are names of objects that will be assigned in the function.
#' logflowvar and logcvar will be shifted to have mean zero. POSSIBLY CHANGE THIS?
#'
#' @param formula a GAM formula. See `help("gam", package = "mgcv")` for details.
#' @param data an object of class rcData, generated using `makeModelData`, or a data frame that makeModelData can use.
#' @param timeout the amount of time in seconds to try mgcv::gam before stopping with a timeout error
#' @param ... further arguments passed to `mgcv::gam`
#' @importFrom mgcv gam
#' @importFrom R.utils withTimeout
#' @export
rcgam <- function(formula, data, timeout = 3, ...) {
cl <- match.call()
if (!is(data, "rcData"))
data = makeModelData(data)
formula = as.formula(formula)
out = withTimeout(mgcv::gam(formula = formula, data = data, ...), timeout = timeout,
onTimeout = "error")
out$call = cl
al = attributes(data)
conc = al$transform$cinvert(data$c)
conc.pred = al$transform$cinvert(out$fitted.values)
scoef = mean(conc) / mean(conc.pred)
if (scoef > 1.5) message(paste0("smearing coefficient is high: ", scoef))
fitted.retrans = conc.pred * scoef
resid.retrans = conc - fitted.retrans
NSE = 1 - sum(resid.retrans ^ 2) / sum((conc - mean(conc))^2)
newbits = list(data = data,
stats = c(al$stats),
yname = "conc",
smearCoef = scoef,
transform = al$transform,
units = al$units,
fitted.retrans = fitted.retrans,
resid.retrans = resid.retrans, NSE = NSE)
out = c(out, newbits)
structure(out, class = c("rcgam", "gam", "glm", "lm"))
}
#' Predict method for rcgam model objects
#'
#' @param object an object of type rcgam, with which to make predictions
#' @param smear Unbias the predictions using the smearing coefficient? Defaults to TRUE
#' @param retrans Retransform the predictions into the original units? Defaults to TRUE
#' @param type What to predict. If not "response" (the default), then the following are coerced to be FALSE:
#' retrans, restrict, smear. See ?predict.gam.
#' @param restrict Restrict the range of predictions to those observed in the calibration data? Defaults to FALSE.
#' @param ... Other arguments passed to mgcv::predict.gam
#'
#'
#' @return a list with element "fit" and potentially others,
#' depending on optional arguments passed. See help("predict.gam", package = "mgcv") for more info
#' The difference is that here the returned value will always be a list, possibly of length 1.
#' @export
predict.rcgam <- function(object, smear = TRUE, retransform = TRUE,
type = "response", restrict = FALSE,
what = c("concentration", "load"), ...) {
arglist = list(...)
what = match.arg(what)
if ("newdata" %in% names(arglist) && !is(arglist$newdata, "rcData"))
arglist$newdata = makePredData(arglist$newdata, object = object)
predfun <- get("predict.gam", asNamespace("mgcv"))
out <- do.call("predfun", args = c(list(object = object, type = type), arglist))
if (type != "response") {
if (smear || retransform || restrict)
warning("The following arguments are not allowed to be TRUE when type != 'response': smear, retransform, restrict")
if (what == "load")
warning("type != response -- coercing what to be 'concentration'.")
return(out)
}
if (!is(out, "list"))
out = list(fit = out)
if (restrict) {
cmax <- max(object$model$c)
cmin <- min(object$model$c)
out$fit[out$fit > cmax] <- cmax
out$fit[out$fit < cmin] <- cmin
}
if (retransform){
out$fit = object$transform$cinvert(out$fit)
if (smear)
out$fit = out$fit * object$smearCoef
if (what == "load") {
newdata <- arglist[["newdata"]]
if (is.null(newdata))
stop("newdata argument must be supplied for load prediction.")
flow <- attr(newdata, "transform")$qinvert(newdata$q)
assert_that(is(flow, "numeric"))
outload <- loadTS(flow = flow, conc = out$fit, datetime = newdata$Date,
load.units = "kg/day",
flow.units = attr(newdata, "units")["qunits"],
conc.units = attr(newdata, "units")["cunits"])
out$fit <- outload$load
}
} else if (smear)
warning("smearing not applicable with non-retransformed predictions. Ignoring this argument.")
out
}
#' Convert an rcgam to a gam
#' Useful for applying gam-specific methods that may not work with e.g. predict.rcgam
#' @export
as.gam <- function(object) {
stopifnot(is(object, "rcgam"))
toRemove <- c("data", "stats", "yname", "smearCoef", "transform", "units",
"fitted.retrans", "resid.retrans", "NSE")
out <- within.list(object, rm(list = toRemove))
nulls <- vapply(out, is.null, logical(1))
out <- structure(out[!nulls], class = c("gam", "glm", "lm"))
out
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.