#' @title FitModel
#'
#' @description This function fits a generalized additive model
#' to a loads dataset
#'
#' @usage FitModel(x, parms = list(flow = FALSE, seasonal = FALSE,
#' RFlimb = FALSE, MA = c(MA1day = FALSE, MA2days = FALSE,
#' MAweek = FALSE, MAmonth = FALSE, MA6months = FALSE,
#' MA12months = FALSE), trend = FALSE, correlation = FALSE))
#'
#' @param x The loads dataset, CQ object resulting from the \code{\link{CreateData}}
#' function.
#' @param parms A list comprising terms to fit in a model:
#' flow: linear (\code{"linear"}) and quadratic (\code{"quadratic"}) terms for flow.
#' Default: \code{FALSE}
#' seasonal: logical term indicating whether seasonal terms are fitted. Default:\code{FALSE}
#' RFlimb: logical term indicating whether a rising/falling limb term is fitted.
#' Default: \code{FALSE},
#' MA: vector indicating whether one or more moving average terms are fitted.
#' Default:\code{c(MA1day = FALSE, MA2days = FALSE, MAweek = FALSE, MAmonth = FALSE,
#' MA6months = FALSE, MA12months = FALSE)}
#' trend: logical term indicating whether a long term trend is fitted. Default: \code{FALSE}
#' correlation: AR1 correlation structure incorporated. Default: \code{FALSE}
#'
#' @details This is a wrapper function for fitting a generalised additive model using the
#' \code{link{mgcv}} package. Terms are fitted in the model by "switching" terms on and off.
#' Apart from the \code{flow} term, all other terms are logical. If a term is intended to be
#' included in a model, then that term can be set to \code{TRUE}.
#'
#' To include additional terms to those specified, the user will need to run the \code{gam}
#' function from the \code{\link{mgcv}} package that this method is based upon. The following is
#' an example where an additional variable, "dam" is included into the model.
#' \code{mod <- gam(log(Conc) ~ log(pQ) + I(log(pQ)^2) + per(days, k=2) + s(MA1day) + s(trend) + dam,
#' family = gaussian, data = loaddata$CQ, control = gam.control(keepData = TRUE)) }
#' where \code{loaddata$CQ} contains an additional column of data labelled, "dam".
#'
#' @return an object of class \code{fitmodel} with elements as described in
#' \code{gamObject} from the \code{mgcv} package. Standard \code{summary}
#' and \code{anova} functions apply.
#'
#' @references Kuhnert, P.M., Henderson, B.L., Lewis, S.E.,
#' Bainbridge, Z.T., Wilkinson, S.N. and Brodie, J.E. (2012)
#' Quantifying total suspended sediment export from the Burdekin
#' River catchment using the loads regression estimator tool,
#' Water Resources Research, 48, W04533,doi:10.1029/2011WR011080.
#'
#' @importFrom "methods" "Quote"
#' @import "mgcv"
#'
#' @seealso ReadInData, CreateData, gam, diagnostic
#'
#' @keywords gam, loads
#'
#' @export
FitModel <- function(x, parms = list(flow = FALSE, seasonal = FALSE,
RFlimb = FALSE, MA = c(MA1day = FALSE, MA2days = FALSE, MAweek = FALSE, MAmonth = FALSE, MA6months = FALSE, MA12months = FALSE),
trend = FALSE, correlation = FALSE)){
# Might need to put a clause in here to take into account models with only
# linear terms. If this happens, use glm and gls to fit the models with and
# without autocorrelation. This means changing the plot, summary and diagnostic functions
if(parms$RFlimb == TRUE & length(table(x$limb)) != 3){
parms$RFlimb <- FALSE
warning("RFlimb set to FALSE because not all limb categories represented. \n")
}
xterms <- mod <- list()
xknots <- NULL
######################################################################################
# Flow terms
######################################################################################
if(parms$flow == "quadratic")
xterms$flow <- Quote(LQflow(pQ, quad = TRUE))
else if(parms$flow == "linear")
xterms$flow <- Quote(LQflow(pQ, quad = FALSE))
else{
warning("Model specified is not listed. Defaulting to linear term for flow only.\n")
xterms$flow <- Quote(LQflow(pQ, quad = FALSE))
}
######################################################################################
# Seasonal term
######################################################################################
# seasonal
if(parms$seasonal){
xterms$seasonal <- Quote(s(month, bs = "cc"))
xknots <- Quote(list(month = seq(1, 13, length = 10)))
}
######################################################################################
# Rising falling limb term
######################################################################################
# Rising/Falling Limb
if(parms$RFlimb)
xterms$limb <- Quote(limb)
######################################################################################
# MA terms
######################################################################################
xterms$MA <- vector(mode = "list", length = 6)
names(xterms$MA) <- c("MA1day", "MA2days", "MAweek", "MAmonth", "MA6months", "MA12months")
# MA terms
if(parms$MA["MA1day"])
xterms$MA$MA1day <- Quote(s(MA1day, bs = "cr"))
if(parms$MA["MA2days"])
xterms$MA$MA2days <- Quote(s(MA2days, bs = "cr"))
if(parms$MA["MAweek"])
xterms$MA$MAweek <- Quote(s(MAweek, bs = "cr"))
if(parms$MA["MAmonth"])
xterms$MA$MAmonth <- Quote(s(MAmonth, bs = "cr"))
if(parms$MA["MA6months"])
xterms$MA$MA6months <- Quote(s(MA6months, bs = "cr"))
if(parms$MA["MA12months"])
xterms$MA$MA12months <- Quote(s(MA12months, bs = "cr"))
######################################################################################
# Trend term
######################################################################################
# Trend
if(parms$trend)
xterms$trend <- Quote(s(trendY, bs = "cr"))
######################################################################################
# Fitting model (with and without correlation)
######################################################################################
if(parms$correlation){
cat("Fitting a Generalized Additive Model with an AR1 correlation structure. \nThis can take time for longer time series ...\n")
fitCommand <- Quote({
mod <- try(gamm(as.numeric(log(Conc)) ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10, family = gaussian, data = x,
knots = k, correlation = corAR1(form = ~ time.index), keep.data = TRUE), silent = TRUE)
})
}
else{
cat("Fitting a Generalized Additive Model ...\n")
fitCommand <- Quote({
mod <- try(gam(log(Conc) ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10, family = gaussian, data = x,
knots = k, control = gam.control(keepData = TRUE)), silent = TRUE)
})
}
thisFitCommand <- do.call("substitute", list(fitCommand,
list(x1 = xterms$flow, x2 = xterms$seasonal, x3 = xterms$limb,
x4 = xterms$MA[[1]], x5 = xterms$MA[[2]], x6 = xterms$MA[[3]],
x7 = xterms$MA[[4]], x8 = xterms$MA[[5]], x9 = xterms$MA[[6]],
x10 = xterms$trend, k = xknots)))
eval(thisFitCommand)
oldClass(mod) <- c("fitmodel", class(mod))
mod
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.