| Title: | River Load Estimation |
|---|---|
| Description: | Collection of functions to make constituent load estimations based on the LOADEST program. |
| Authors: | Dave Lorenz [cre], Rob Runkel [aut], Laura De Cicco [aut], David Watkins [ctb], Joseph Stachelek [ctb] |
| Maintainer: | Dave Lorenz <[email protected]> |
| License: | CC0 |
| Version: | 0.4.5 |
| Built: | 2024-10-27 04:20:37 UTC |
| Source: | https://github.com/ldecicco-USGS/rloadest |
| Package: | rloadest |
| Type: | Package |
| License: | CC0 |
| LazyLoad: | yes |
This package is intended to replicate and extend the LOADEST program for estimating constituent loads in streams and rivers. Some subtle differences between the output for LOADEST and rloadest include:
The least absolute deviation (LAD) method is not supported in rloadest.
LOADEST uses centered time when computing the sine and cosine terms in model numbers 4, 6, 7, 8, and 9, but the functions in rloadest use the actual decimal time so that the seasonality can more easily be assessed by the user.
The order of the terms in the predefined models is different between LOADEST and the rloadest functions.
The printed output of the model descriptions from rloadest matches the format the most users of R would recognize from other linear model output rather then the printed output from LOADEST.
Furthermore, the model building capability in the rloadest functions make easier to explore other forms of rating-curve models than LOADEST.
Compute Akaike's An Information Criterion with Correction (AICc) for for finite sample sizes.
AICc(object)AICc(object)
object |
the output from |
A numeric value corresponding to the AICc of object.
The penalty that AIC applies for adding explanatory variables
is biased low when the number of samples is small. As a result, models with
small smaple sizes can be overfitted. AICc can be used to identify more
parsimonious models.
Hurvitch, C.M. and Tsai, C.L., 1989, Regression and time series model selection in small samples: Biometrika, v. 76, no. 2, p. 297–307.
# From application 1 in the vignettes data(app1.calib) app1.lr <- loadReg(Phosphorus ~ model(1), data = app1.calib, flow = "FLOW", dates = "DATES", conc.units="mg/L", station="Illinois River at Marseilles, Ill.") AICc(app1.lr)# From application 1 in the vignettes data(app1.calib) app1.lr <- loadReg(Phosphorus ~ model(1), data = app1.calib, flow = "FLOW", dates = "DATES", conc.units="mg/L", station="Illinois River at Marseilles, Ill.") AICc(app1.lr)
Illinois River at Marseilles, Illinois (Helsel & Hirsch, 2002)
app1.calibapp1.calib
Data frame with 96 rows and 4 columns
| Name | Type | Description |
| DATES | Date | Date of daily value |
| TIMES | character | Time corresponding to noon of daily value |
| FLOW | numeric | Daily mean streamflow |
| Phosphorus | numeric | Daily mean phosphorus concentration (assumed) |
Example calibration dataset from LOADEST
Runkel, R.G., Crawford, C.G., and Cohn, T.A., 2004, Load Estimator (LOADEST): a FORTRAN program for estimating constituent loads in streams and rivers: U.S. Geological Survey Techniques and Methods Book 4, Chapter A5, 69 p.
## Not run: data(app1.calib) # Plot concentration vs. flow with(app1.calib, plot(FLOW, Phosphorus, log="xy")) ## End(Not run)## Not run: data(app1.calib) # Plot concentration vs. flow with(app1.calib, plot(FLOW, Phosphorus, log="xy")) ## End(Not run)
St.Joseph River near Newville, IN (Station # 04178000)
app2.calibapp2.calib
Data frame with 32 rows and 4 columns
| Name | Type | Description |
| DATES | Date | Date of daily value |
| TIMES | integer | Time that sample was actually taken |
| FLOW | numeric | Daily mean streamflow |
| Atrazine | numeric | Daily mean atrazine concentration (assumed) |
Example calibration dataset from LOADEST
Runkel, R.G., Crawford, C.G., and Cohn, T.A., 2004, Load Estimator (LOADEST): a FORTRAN program for estimating constituent loads in streams and rivers: U.S. Geological Survey Techniques and Methods Book 4, Chapter A5, 69 p.
## Not run: data(app2.calib) # Plot concentration vs. flow with(app2.calib, plot(FLOW, Atrazine, log="xy")) ## End(Not run)## Not run: data(app2.calib) # Plot concentration vs. flow with(app2.calib, plot(FLOW, Atrazine, log="xy")) ## End(Not run)
St.Joseph River near Newville, IN (Station # 04178000)
app2.estapp2.est
Data frame with 730 rows and 3 columns
| Name | Type | Description |
| DATES | Date | Date of daily value |
| TIMES | integer | Time corresponding to noon of daily value |
| FLOW | numeric | Daily mean streamflow |
Example estimation dataset from LOADEST
Runkel, R.G., Crawford, C.G., and Cohn, T.A., 2004, Load Estimator (LOADEST): a FORTRAN program for estimating constituent loads in streams and rivers: U.S. Geological Survey Techniques and Methods Book 4, Chapter A5, 69 p.
## Not run: data(app2.est) summary(app2.est) ## End(Not run)## Not run: data(app2.est) summary(app2.est) ## End(Not run)
White River at Hazleton, Ind. (Station Number 03374100)
app4.calibapp4.calib
Data frame with 45 rows and 9 columns
| Name | Type | Description |
| DATES | Date | Date of daily value |
| TIMES | character | Time that sample was actually taken |
| FLOW | numeric | Daily mean streamflow |
| Buty.rmk | character | Remark code for butylate concentration |
| Buty | numeric | Daily mean butylate concentration (assumed) |
| Atra | numeric | Daily mean atrazine concentration (assumed) |
| Alach.rmk | character | Remark code for alachlor concentration |
| Alach | numeric | Daily mean alachlor concentration (assumed) |
| SuspSed | numeric | Daily mean suspended sediment concentration (assumed) |
Obtained from Charlie Crawford, 5 July 2001
Runkel, R.G., Crawford, C.G., and Cohn, T.A., 2004, Load Estimator (LOADEST): a FORTRAN program for estimating constituent loads in streams and rivers: U.S. Geological Survey Techniques and Methods Book 4, Chapter A5, 69 p.
## Not run: data(app4.calib) # Plot atrazine concentration vs. flow with(app4.calib, plot(FLOW, Atra, log="xy")) ## End(Not run)## Not run: data(app4.calib) # Plot atrazine concentration vs. flow with(app4.calib, plot(FLOW, Atra, log="xy")) ## End(Not run)
White River at Hazleton, Ind. (Station Number 03374100)
app4.estapp4.est
Data frame with 730 rows and 3 columns
| Name | Type | Description |
| DATES | Date | Date of daily value |
| TIMES | character | Time corresponding to noon of daily value |
| FLOW | numeric | Daily mean streamflow |
Obtained from Charlie Crawford, 5 July 2001
Runkel, R.G., Crawford, C.G., and Cohn, T.A., 2004, Load Estimator (LOADEST): a FORTRAN program for estimating constituent loads in streams and rivers: U.S. Geological Survey Techniques and Methods Book 4, Chapter A5, 69 p.
## Not run: data(app4.est) summary(app4.est) ## End(Not run)## Not run: data(app4.est) summary(app4.est) ## End(Not run)
Little Arkansas River near Halstead, Kansas (Station Number 07143672)
app5.1999app5.1999
Data frame with 13 rows and 3 columns
| Name | Type | Description |
| DATES | Date | Date of daily value |
| TIMES | character | Time that sample was actually taken |
| Alkalinity | numeric | Daily mean alkalinity (assumed) |
Retrieved from NWISweb on July 26, 2013 from URL: https://nwis.waterdata.usgs.gov/ks/nwis/qwdata
data(app5.1999) head(app5.1999)data(app5.1999) head(app5.1999)
Little Arkansas River near Halstead, Kansas (Station Number 07143672)
app5.calibapp5.calib
Data frame with 103 rows and 5 columns
| Name | Type | Description |
| DATES | Date | Date of daily value |
| TIMES | character | Time that sample was actually taken |
| FLOW | numeric | Daily mean streamflow |
| SC | numeric | Daily mean specific conductance |
| Alkalinity | numeric | Daily mean alkalinity (assumed) |
Example calibration dataset from LOADEST
Runkel, R.G., Crawford, C.G., and Cohn, T.A., 2004, Load Estimator (LOADEST): a FORTRAN program for estimating constituent loads in streams and rivers: U.S. Geological Survey Techniques and Methods Book 4, Chapter A5, 69 p.
## Not run: data(app5.calib) # Plot concentration vs. specific conductance with(app5.calib, plot(SC, Alkalinity, log="xy")) ## End(Not run)## Not run: data(app5.calib) # Plot concentration vs. specific conductance with(app5.calib, plot(SC, Alkalinity, log="xy")) ## End(Not run)
Little Arkansas River near Halstead, Kansas (Station Number 07143672)
app5.estapp5.est
Data frame with 358 rows and 4 columns
| Name | Type | Description |
| DATES | Date | Date of daily value |
| TIMES | character | ime corresponding to noon of daily value |
| FLOW | numeric | Daily mean streamflow |
| SC | numeric | Daily mean specific conductance |
Example estimation dataset from LOADEST
Runkel, R.G., Crawford, C.G., and Cohn, T.A., 2004, Load Estimator (LOADEST): a FORTRAN program for estimating constituent loads in streams and rivers: U.S. Geological Survey Techniques and Methods Book 4, Chapter A5, 69 p.
## Not run: data(app5.est) summary(app5.est) ## End(Not run)## Not run: data(app5.est) summary(app5.est) ## End(Not run)
Convert concentration and flow to load (flux).
c2load(conc, flow, flow.units = "cfs", conc.units = "", load.units = "kg", ignore.censoring = TRUE)c2load(conc, flow, flow.units = "cfs", conc.units = "", load.units = "kg", ignore.censoring = TRUE)
conc |
the concentration data missing values are permitted and result in missing values in the output. |
flow |
the flow data missing values are permitted and result in missing values in the output. |
flow.units |
character string describing the flow unit. |
conc.units |
character string describing the concentration unit. |
load.units |
character string describing the load unit. |
ignore.censoring |
logical, see Value. |
If ignore.censoring is TRUE, the default,
then return a vector of numeric values with censored values replaced
by 1/2 the detection limit. Otherwise, return a vector that retains
the censoring—if conc is numeric, then uncensored; if
conc is of class "qw," then the returned data would be of
class "lcens" or "mcens."
will need some.
# These calls return the conversion factors c2load(1, 1, conc.units="mg/L") c2load(1, 1, conc.units="mg/L", load.units="tons")# These calls return the conversion factors c2load(1, 1, conc.units="mg/L") c2load(1, 1, conc.units="mg/L", load.units="tons")
Gets the type of censoring ("none," "left," "multiple") for an object.
## S3 method for class 'Surv' censoring(x)## S3 method for class 'Surv' censoring(x)
x |
the object to get the type of censoring. For an object of class "Surv," the type must be "interval." |
A character string "none," "left," or "multiple" describing the type of censoring present in the object.
This function is mostly used within other functions to determine the 'best' technique to use for analysis.
## Not run: library(survival) censoring(Surv(2.3, 2.3, type="interval2")) ## End(Not run)## Not run: library(survival) censoring(Surv(2.3, 2.3, type="interval2")) ## End(Not run)
Computes centered values. Used primarily in a linear regression formula.
center(x, center = NULL)center(x, center = NULL)
x |
a numeric vector for which to compute the centered values. Missing values are permitted and result in corresponding missing values in the output. |
center |
an optional value to use for the center of |
The centered value of x.
The centering value by default is computed by the method described in Cohn and others (1992).
Cohn, T.A., Caulder, D.L., Gilroy, E.J., Zynjuk, L.D., and Summers, R.M., 1992, The validity of a simple statistical model for estimating fluvial constituent loads—An empirical study involving nutrient loads entering Chesapeake Bay: Water Resources research, v. 28, no. 5, p. 937–942.
# trivial ceneterd values from 1 to 10 center(seq(10))# trivial ceneterd values from 1 to 10 center(seq(10))
Extract the model coefficients from a load regression.
## S3 method for class 'loadReg' coef(object, summary = FALSE, which = "load", ...)## S3 method for class 'loadReg' coef(object, summary = FALSE, which = "load", ...)
object |
the output from |
summary |
include standard errors and other information? |
which |
string indicating which coefficients to return; "load" returns the load model and "concentration" returns the concentration model coefficients. |
... |
further arguments passed to or from other methods. |
Either a names vector of the coefficients, if summary is
FALSE or a matrix of the coefficients, their standard errors,
z-scores, and attained p-values, if summary is TRUE.
The attained p-values are computed from the log-likelihood test for AMLE regression and from a Wald chi-square test for MLE regression.
# From application 1 in the vignettes data(app1.calib) app1.lr <- loadReg(Phosphorus ~ model(1), data = app1.calib, flow = "FLOW", dates = "DATES", conc.units="mg/L", station="Illinois River at Marseilles, Ill.") # Extract the coefficients coef(app1.lr)# From application 1 in the vignettes data(app1.calib) app1.lr <- loadReg(Phosphorus ~ model(1), data = app1.calib, flow = "FLOW", dates = "DATES", conc.units="mg/L", station="Illinois River at Marseilles, Ill.") # Extract the coefficients coef(app1.lr)
Compute daily mean water-quality values for data frames containing water-quality data.
dailyAg(x, dates = "sample_dt", times = "sample_tm")dailyAg(x, dates = "sample_dt", times = "sample_tm")
x |
a data frame containing water-quality data. |
dates |
the name of the sample date column. |
times |
the name of the sample time column. |
A data frame like x but with the means for each
column by day and the sample time set to "1200."
Extract the fitted values of a load regression.
## S3 method for class 'loadReg' fitted(object, suppress.na.action = FALSE, which = "load", ...)## S3 method for class 'loadReg' fitted(object, suppress.na.action = FALSE, which = "load", ...)
object |
an object of class "loadReg"—output from |
suppress.na.action |
logical, suppress the effects of the
|
which |
a character string indicating the type of fitted values. Must be either "load" or "concentration." |
... |
further arguments passed to or from other methods. |
The fitted values from the regression. Note that these are not back- transformed but are in natrual log units.
# From application 1 in the vignettes data(app1.calib) app1.lr <- loadReg(Phosphorus ~ model(1), data = app1.calib, flow = "FLOW", dates = "DATES", conc.units="mg/L", station="Illinois River at Marseilles, Ill.") # Extract the fitted values fitted(app1.lr)# From application 1 in the vignettes data(app1.calib) app1.lr <- loadReg(Phosphorus ~ model(1), data = app1.calib, flow = "FLOW", dates = "DATES", conc.units="mg/L", station="Illinois River at Marseilles, Ill.") # Extract the fitted values fitted(app1.lr)
Compute selected jackknife statistics for a rating-curve load-estimation model.
jackStats(fit, which = "load")jackStats(fit, which = "load")
fit |
an object of class "loadReg"—output from |
which |
a character string indicating the "load" or "concentration" model for an object of class "loadReg" or "censReg" for an object of class "censReg." |
An object of class "jackStats" containing these components:
coef, the table of coefficient estimates, the jackknife bias and standard errors
coefficients, the jackknifed coefficients
pctcens, the percentage of left-censored values.
The PRESS statistic and individual jackknife differences are also returned
when the percentage of censoring is 0.
The jackStats function can only be used when the analysis is AMLE.
Abdi and Williams (2010) describe the jackknife as refering to two related techniques: the first
estimates the parameters, their bias and standard errors and the second evaluates the
predictive performance of the model. The second technique is the PRESS statistic (Helsel
and Hirsch, 2002), but can only be used on uncensored data; it is computed by jackStats
when no data are censored. The first technique can be used to assess the coefficients of the
regression—the bias should be small and the jackknife standard errors should not be much
different from the standard errors reported for the regression. Efron and Tibshirani (1993)
suggest that the bias is small if the relative bias (biuas divided by the jackknife standard
error) is less than 0.25.
Abdi, H. and Williams, L.J., 2010, Jackknife, in encyclopedia of research design, Salkind, N.J., editor: Thousand Oaks, Calif., SAGE Publications, 1719 p.
Efron, B. and Tibshirani, R.J., 1993, An introduction to the bootstrap: Boca Raton, Fla., Chapman and Hall/CRC, 436 p.
Helsel, D.R. and Hirsch, R.M., 2002, Statistical methods in water resources: U.S. Geological Survey Techniques of Water-Resources Investigations, book 4, chap. A3, 522 p. Salkind,
# From application 1 in the vignettes data(app1.calib) app1.lr <- loadReg(Phosphorus ~ model(1), data = app1.calib, flow = "FLOW", dates = "DATES", conc.units="mg/L", station="Illinois River at Marseilles, Ill.") jackStats(app1.lr)# From application 1 in the vignettes data(app1.calib) app1.lr <- loadReg(Phosphorus ~ model(1), data = app1.calib, flow = "FLOW", dates = "DATES", conc.units="mg/L", station="Illinois River at Marseilles, Ill.") jackStats(app1.lr)
Computes the conversion factor to compute the flux units (load units) from the concentration and flow units
loadConvFactor(flow.units, conc.units, load.units)loadConvFactor(flow.units, conc.units, load.units)
flow.units |
character string describing the flow unit. The only valid values are "cubic meter per second," "cms," "cubic feet per second," or "cfs." |
conc.units |
character string describing the concentration unit. The valid units are "mg/l," "mg/L," "ug/l," "ug/L," "ng/l," "ng/L," "milligrams per liter," "micrograms per liter," "nanograms per liter," "col/100mL," "col/dL," or "colonies per 100mL." |
load.units |
character string describing the load unit. The valid values are "pounds," "tons," "mg," "milligrams," "grams," "g," "kilograms," "kg," "metric tons," "Mg," or "million colonies." |
The conversion factor.
loadConvFactor("cubic meter per second","milligrams per liter","pounds")loadConvFactor("cubic meter per second","milligrams per liter","pounds")
Internal support function for rloadest that computes the adjustment to flow.
loadestQadj(x, round = options("round.flow"))loadestQadj(x, round = options("round.flow"))
x |
the calibration flow data. |
round |
either a numeric value indicating the number
of decimal places, or a list containing a value indicating
the number of decimal places. If |
The centering value for flow.
Internal support function for rloadest that computes the adjustment to time.
loadestTadj(x, round = options("round.time"))loadestTadj(x, round = options("round.time"))
x |
the calibration date data. |
round |
either a numeric value indicating the number
of decimal places, or a list containing a value indicating
the number of decimal places. If |
A vector of length 2 containing the base (reference) year and the centering time correction value. The user value for centered time would be sum.
Build a rating-curve (regression) model for river load estimation.
loadReg(formula, data, subset, na.action, flow, dates, flow.units = "cfs", conc.units = "", load.units = "kg", time.step = "day", station = "")loadReg(formula, data, subset, na.action, flow, dates, flow.units = "cfs", conc.units = "", load.units = "kg", time.step = "day", station = "")
formula |
a formula describing the regression model. See Details. |
data |
the data to search for the variables in |
subset |
an expression to select a subset of the data. |
na.action |
what to do with missing values. |
flow |
character string indicating the name of the flow column. |
dates |
character string indicating the name of the date column. |
flow.units |
character string describing the flow units. See Details. |
conc.units |
character string describing the concentration unit. See Details. |
load.units |
character string describing the load unit. See Details. |
time.step |
character string describing the time step of the calibration data. Must be one of "instantaneous," "2 hours," "3 hours," "4 hours," "6 hours," "12 hours," or "day." The default is "day." |
station |
character string description of the station. |
The left-hand side of the formula may be any numeric variable, just as with
lm or a variable of class "lcens," "mcens," or "qw." Also permitted are
variables constructed using Surv of type "right,", "interval," or
"interval2" (for left-censored data, use as.lcens.
For un- or left-censored data, AMLE is used unless weights are specified in
the model, then MLE is used, through a call to survreg. For any other
censored data, MLE is used.
Typically, loadReg expects the response variable to have units of
concentration, mass per volume. For these models, See loadConvFactor
for details about valid values for flow.units, conc.units and
load.units. For some applications, like bed load estimation, the response
variable can have units of flux, mass per time. For these models, conc.units
can be expressed as any valid load.units per day. The rate must be expressed
in terms of "/d," "/dy," or "/day."
An object of class "loadReg."
Runkel, R.L., Crawford, C.G., and Cohn, T.A., 2004, Load estimator (LOADEST): a FORTRAN program for estimating constituent loads in streams and rivers: U.S. Geological Survey Techniques and Methods book 4, chap. A5, 69 p.
censReg, as.lcens, as.mcens,
Surv
# From application 1 in the vignettes data(app1.calib) app1.lr <- loadReg(Phosphorus ~ model(1), data = app1.calib, flow = "FLOW", dates = "DATES", conc.units="mg/L", station="Illinois River at Marseilles, Ill.") print(app1.lr)# From application 1 in the vignettes data(app1.calib) app1.lr <- loadReg(Phosphorus ~ model(1), data = app1.calib, flow = "FLOW", dates = "DATES", conc.units="mg/L", station="Illinois River at Marseilles, Ill.") print(app1.lr)
Create a 2-page pdf file report of a rating-curve load model. The report contains the text output and 6 diagnostic plots.
loadReport(x, file)loadReport(x, file)
x |
the load model. |
file |
the output file base name; the .pdf suffix
is appended to make the actual file name. if missing, then the
name of |
The actual file name is returned invisibly.
Compute some summary statistics for a rating-curve load-estimation model.
loadStats(fit, which = "load")loadStats(fit, which = "load")
fit |
an object of class "loadReg"—output from |
which |
a character string indicating the "load" or "concentration" model. |
A list containing outSum, selected summary statistics;
of the observed and estimated values and outBias, the bias
statistics.
# From application 1 in the vignettes data(app1.calib) app1.lr <- loadReg(Phosphorus ~ model(1), data = app1.calib, flow = "FLOW", dates = "DATES", conc.units="mg/L", station="Illinois River at Marseilles, Ill.") loadStats(app1.lr)# From application 1 in the vignettes data(app1.calib) app1.lr <- loadReg(Phosphorus ~ model(1), data = app1.calib, flow = "FLOW", dates = "DATES", conc.units="mg/L", station="Illinois River at Marseilles, Ill.") loadStats(app1.lr)
Computes the factor to convert from between load units.
loadUnitConv(from, to)loadUnitConv(from, to)
from |
character string describing the load units to convert from. |
to |
character string describing the load units to convert to. |
The conversion factor.
loadUnitConv("kilograms", "tons")loadUnitConv("kilograms", "tons")
Compute the log-likelihood statistics for a load regression.
## S3 method for class 'loadReg' logLik(object, ...)## S3 method for class 'loadReg' logLik(object, ...)
object |
the output from |
... |
further arguments passed to or from other methods. |
An object of class "logLik" containing the log-likelihood and the attributes "df" (degrees of freedom) and "nobs" (number of observations).
# From application 1 in the vignettes data(app1.calib) app1.lr <- loadReg(Phosphorus ~ model(1), data = app1.calib, flow = "FLOW", dates = "DATES", conc.units="mg/L", station="Illinois River at Marseilles, Ill.") logLik(app1.lr)# From application 1 in the vignettes data(app1.calib) app1.lr <- loadReg(Phosphorus ~ model(1), data = app1.calib, flow = "FLOW", dates = "DATES", conc.units="mg/L", station="Illinois River at Marseilles, Ill.") logLik(app1.lr)
A utility to help model.frame.default create the right matrices
when predicting from models with center term. Used only internally.
## S3 method for class 'center' makepredictcall(var, call)## S3 method for class 'center' makepredictcall(var, call)
var |
a variable. |
call |
the term in the formula, as a call. |
A replacement for call for the prediction variable.
Method functions to compute the "mean" for factors and characters. These
funcitons are intended primarily as support functions when aggregating
unit values data in predLoad.
## S3 method for class 'factor' mean(x, ...) ## S3 method for class 'character' mean(x, ...)## S3 method for class 'factor' mean(x, ...) ## S3 method for class 'character' mean(x, ...)
x |
an object of class "factor" or "character." |
... |
further arguments passed to or from other methods. |
If all values are identical, then the unique value, otherwise
the missing value (NA).
Support function for building a predefined rating curve load model.
model(model.no, data, flow, time)model(model.no, data, flow, time)
model.no |
the model number. |
data |
the dataset. |
flow |
character string indicating the name of the flow column. |
time |
character string indicating the name of the date/time column. |
Row number to select from data to build a predefined model.
Models A dataset that describes the equivalent formula for each of the predefined model number.
ModelsModels
Data frame with 9 rows and 2 columns
| Name | Type | Description |
| Number | integer | The predefined model number |
| Formula | factor | The equivalent formula for the predefined model number |
data(Models) print(Models)data(Models) print(Models)
Compute the Nash-Sutcliffe efficiency rating for model estiamtes.
nashSutcliffe(obs, est, na.rm = TRUE)nashSutcliffe(obs, est, na.rm = TRUE)
obs |
a vector of the observed values. |
est |
a vector of the model estiamted values. Each value must
pair with each value in |
na.rm |
remove missing values from |
A single numeric value representing the Nash-Sutcliffe efficiency rating for the observed and estiamted data.
Example data representing atrazine
Plot rating-curve load model diagnostics.
## S3 method for class 'loadReg' plot(x, which = "All", set.up = TRUE, span = 1, ...)## S3 method for class 'loadReg' plot(x, which = "All", set.up = TRUE, span = 1, ...)
x |
an object of class "loadReg"—output from |
which |
either "All" or any of a sequence from 1 to 7 indicating which plot, see Details. |
set.up |
set up the graphics page? |
span |
the span to use for the loess smooth. Set to 0 to suppress. |
... |
further arguments passed to or from other methods. |
Seven graphs can be produced by this function. If which is "All," then all plots are produced.
The argument which can also be the name of an explanatory variable so that a partial residual
plot is created for a single variable. Or which can be any of a sequence of numbers from 1 thorugh 7.
Numeric values for which:
Observed vs. fitted.
Fitted vs. Residual
S-L plot
A correlogram if dates are available in the model or in the data set
Q-normal
Tukey boxplots for oberved and estimated
Partial residual plots for each explanatory variable
The object x is returned invisibly.
This plotting function uses the core routines in the smwrGraphs
package. It requires a graphics page that is set up from the functions
in that package (setpage or setPDF) if set.up is
FALSE. The graphs that are produced by this function are based
on the publication guidelines of the USGS.
# From application 1 in the vignettes ## Not run: data(app1.calib) app1.lr <- loadReg(Phosphorus ~ model(1), data = app1.calib, flow = "FLOW", dates = "DATES", conc.units="mg/L", station="Illinois River at Marseilles, Ill.") # Produce the full suite of diagnostic plots plot(app1.lr) ## End(Not run)# From application 1 in the vignettes ## Not run: data(app1.calib) app1.lr <- loadReg(Phosphorus ~ model(1), data = app1.calib, flow = "FLOW", dates = "DATES", conc.units="mg/L", station="Illinois River at Marseilles, Ill.") # Produce the full suite of diagnostic plots plot(app1.lr) ## End(Not run)
Estimate concentrations from a rating-curve model from loadReg
for a new data frame.
predConc(fit, newdata, by = "day", allow.incomplete = FALSE, conf.int = 0.95)predConc(fit, newdata, by = "day", allow.incomplete = FALSE, conf.int = 0.95)
fit |
the output from |
newdata |
a data frame of the prediction variables. MIssing values
are not permitted in any column in |
by |
the time frame for estimates. See |
allow.incomplete |
compute loads for periods withing
missing values or incomplete record? See |
conf.int |
the confidence interval to compute for concentrations.
See |
The time frame specified by by must be either "unit" or
"day."
If allow.incomplete is TRUE, then concentrations will be
computed based on all nonmissing values, otherwise missing values
NAs will be returned. For this application, missing values
includes NAs and incomplete days. For prediction by "day" when
there are variable number of unit values per day, allow.incomplete
must be set to TRUE.
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.
A data frame containing the concnetration estimates.
# From application 1 in the vignettes data(app1.calib) app1.lr <- loadReg(Phosphorus ~ model(1), data = app1.calib, flow = "FLOW", dates = "DATES", conc.units="mg/L", station="Illinois River at Marseilles, Ill.") predConc(app1.lr, app1.calib)# From application 1 in the vignettes data(app1.calib) app1.lr <- loadReg(Phosphorus ~ model(1), data = app1.calib, flow = "FLOW", dates = "DATES", conc.units="mg/L", station="Illinois River at Marseilles, Ill.") predConc(app1.lr, app1.calib)
Estimate loads from a rating-curve model from loadReg
for a new data frame, aggregating the loads by specified time
periods.
predLoad(fit, newdata, load.units = fit$load.units, by = "total", seopt = "exact", allow.incomplete = FALSE, conf.int = 0.95, print = FALSE)predLoad(fit, newdata, load.units = fit$load.units, by = "total", seopt = "exact", allow.incomplete = FALSE, conf.int = 0.95, print = FALSE)
fit |
the output from |
newdata |
a data frame of the prediction variables. Missing values
are not permitted in any column in |
load.units |
a character string indicating the units of the
predicted loads/fluxes. By default, uses the value specified in
|
by |
the time frame for estimates. See |
seopt |
a character string indicating how to comute the standard error of the aggregated load estimates, must be either "exact" or "approximate." Only the first letter is necessary. |
allow.incomplete |
compute loads for periods withing
missing values or incomplete record? See |
conf.int |
the confidence interval to compute for loads
computed by "day" or "unit." The confidence interval is fixed at
0.95 for any other value for |
print |
print a report summary of the load estimate? |
The time frame specified by by can be "unit,"
"day," "month," "water year," "calendar year," "total," or
the name of a column in newdata that can be used to
group the data.
If allow.incomplete is TRUE, then loads will be
computed based on all nonmissing values, otherwise missing values
NAs will be returned. For this application, missing values
includes NAs and gaps in the record, except for by
set to "total" or user defined groups where missing values only
includes NAs. For prediction by "day" when there are variable
number of unit values per day, allow.incomplete must be
set to TRUE.
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.
A data frame containing the load estimates.
# From application 1 in the vignettes data(app1.calib) app1.lr <- loadReg(Phosphorus ~ model(1), data = app1.calib, flow = "FLOW", dates = "DATES", conc.units="mg/L", station="Illinois River at Marseilles, Ill.") predLoad(app1.lr, app1.calib)# From application 1 in the vignettes data(app1.calib) app1.lr <- loadReg(Phosphorus ~ model(1), data = app1.calib, flow = "FLOW", dates = "DATES", conc.units="mg/L", station="Illinois River at Marseilles, Ill.") predLoad(app1.lr, app1.calib)
Print the results of a jackknife analysis of a censored regression.
## S3 method for class 'jackStats' print(x, digits = 4, ...)## S3 method for class 'jackStats' print(x, digits = 4, ...)
x |
an object of class "jackStats"—output from |
digits |
the number of significant digits to print. |
... |
further arguments passed to or from other methods. |
The object x is returned invisibly.
The printed output includes the original original estimate, the jackknife bias and standard error and the relative bias for each parameter in the regression model.
Print the results of an load rating-curve regression.
## S3 method for class 'loadReg' print(x, digits = 4, brief = TRUE, load.only = brief, ...)## S3 method for class 'loadReg' print(x, digits = 4, brief = TRUE, load.only = brief, ...)
x |
an object of class "loadReg"—output from |
digits |
the number of significant digits to print. |
brief |
print the brief output? See Note. |
load.only |
print only the load model and not concentration model results. |
... |
further arguments passed to or from other methods. |
The object x is returned invisibly.
The printed output replicates the output described in Runkel (2004) and
includes a short section summarizing the data, the load model and coefficients,
regression statistics, and comparison of observed and estimated loads. If
load.only is set to FALSE, then similar output is generated for the
concetration model. If brief is FALSE, then additional descriptions
of selected sections of the output are produced.
If the estimation method is "MLE," then the estimated loads used in the comparison
to observed loads are approximate becuase they are estimated using MLE, rather than
AMLE, wihch is used for predLoad and predConc. The bias is very small
when the residual varaince is less than 0.5, but can be large when the residual
variance is greater than 1.
Runkel, R.L., Crawford, C.G., and Cohn, T.A., 2004, Load estimator (LOADEST): A FORTRAN program for estimating constituent loads in streams and rivers: U.S. Geological Survey Techniques and Methods book 4, chap. A5, 69 p.
Unit-value data can be recorded at any arbitrary intervals. For some
applications, such as load estimates, a uniform series is required. The
resampleUVdata function resamples the orginal unit-value data to a consistent
time interval.
resampleUVdata(UVdata, time.step = 15, start.date = "", end.date = "", max.diff = "2 hours", missing.days = "exclude")resampleUVdata(UVdata, time.step = 15, start.date = "", end.date = "", max.diff = "2 hours", missing.days = "exclude")
UVdata |
the dataset containing the unit-values data. Must have one column that represents the time of the observation that is class "POSIXt." Missing values are not permitted in that column. |
time.step |
the time step of the new data in minutes; must divide an hour exactly evenly. The default value is 15 minutes. |
start.date |
a character string indicating the first day of the output dataset. The
default value ("") indicates use the first day in |
end.date |
a character string indicating the last day of the output dataset. The
default value ("") indicates use the last day in |
max.diff |
a character string indicating the maximum difference in time to
sucessfully resample the unit-value data. The default is "2 hours" see
|
missing.days |
a characer string indicating what action should be taken for days
not present in |
A data frame like UVdata but having a uniform time step.
Extract the residuals from the load or concetration regression model. The residuals will be the same unless the log of flow is not an explanatory variable.
## S3 method for class 'loadReg' residuals(object, type = "working", suppress.na.action = FALSE, model = c("load", "concentration"), ...)## S3 method for class 'loadReg' residuals(object, type = "working", suppress.na.action = FALSE, model = c("load", "concentration"), ...)
object |
an object of class "loadReg"—output from |
type |
The type of residuals, see Details. |
suppress.na.action |
logical, suppress the effects of the
|
model |
the type of model, must be either "load" or "concentration." |
... |
not used, required for other methods. |
The value for type can be any one of the following:
| Value | Description |
| "working" | Residuals with censored residuals replaced by their expected values |
| "response" | Residuals from the linear predictor |
| "influence" | An estimate of Cook's D values based on "working" residuals |
| "leverage" | The hat diagonals |
| "S-L" | The square-root of the absolute value of the residuals with censored residuals replaced by their expected value |
Also, any other value of type for residuals.survreg
can be used to obtain those residuals. Note that "working" and "response"
are defined in the table above, in keeping with older versions of
loadReg.
The residuals from the regression as specified by type.
The residuals from the load regression are the same as those from the concentration regression, so there is no option to distinguish among those models.
Compute the root-mean-squared error (RMSE) of the difference between observed values and the fitted values for the load or concentration model. The RMSEs will be the same unless the log of flow is not an explanatory variable.
## S3 method for class 'loadReg' rmse(x, model = c("load", "concentration"), ...)## S3 method for class 'loadReg' rmse(x, model = c("load", "concentration"), ...)
x |
the output from |
model |
the type of model, must be either "load" or "concentration." |
... |
not used, required for other methods. |
The estimated root-mean-squared error, also know as the residual standard error.
# From application 1 in the vignettes data(app1.calib) app1.lr <- loadReg(Phosphorus ~ model(1), data = app1.calib, flow = "FLOW", dates = "DATES", conc.units="mg/L", station="Illinois River at Marseilles, Ill.") rmse(app1.lr)# From application 1 in the vignettes data(app1.calib) app1.lr <- loadReg(Phosphorus ~ model(1), data = app1.calib, flow = "FLOW", dates = "DATES", conc.units="mg/L", station="Illinois River at Marseilles, Ill.") rmse(app1.lr)
Support function for building a segmented rating curve load model. Required
in the formula in segLoadReg to define the segmented model.
seg(x, N)seg(x, N)
x |
the data to segment. Missing values are permitted and result corresponing in missing values in output. |
N |
the number of breaks in the segmented model. |
The data in x with attributes to build the segmented model.
Build a segmented rating-curve (regression) model for river load estimation.
segLoadReg(formula, data, subset, na.action, flow, dates, flow.units = "cfs", conc.units = "", load.units = "kg", time.step = "day", station = "", trace = TRUE)segLoadReg(formula, data, subset, na.action, flow, dates, flow.units = "cfs", conc.units = "", load.units = "kg", time.step = "day", station = "", trace = TRUE)
formula |
a formula describing the regression model. See Details. |
data |
the data to search for the variables in |
subset |
an expression to select a subset of the data. |
na.action |
what to do with missing values. |
flow |
character string indicating the name of the flow column. |
dates |
character string indicating the name of the date column. |
flow.units |
character string describing the flow unit. |
conc.units |
character string describing the concentration unit. |
load.units |
character string describing the load unit. |
time.step |
character string describing the time step of the calibration data. |
station |
character string description of the station. |
trace |
if logical, then if |
The left-hand side of the formula can be any numeric variable, just
as with lm or a variable of class "lcens," "mcens," "qw," or "Surv."
The response variable must be a column in data; it cannot be constructed using
as.lcens, as.mcens, or Surv. The initial segmented model is
based on uncensored data—simple substitution of 1/2 the reporting limit is used
for left-censored values and multiply censored values cause the analysis to fail.
The first term of right-hand side must be defined by the seg function
with the number of segments. The first term may be followed by any number of
additional terms. The final model will place the segmeted term in the last postition
and seg will be replaced by the proper call to segment.
An object of class "loadReg."
will need some.
# From application 1 in the vignettes data(app1.calib) app1.lr <- loadReg(Phosphorus ~ model(1), data = app1.calib, flow = "FLOW", dates = "DATES", conc.units="mg/L", station="Illinois River at Marseilles, Ill.") print(app1.lr)# From application 1 in the vignettes data(app1.calib) app1.lr <- loadReg(Phosphorus ~ model(1), data = app1.calib, flow = "FLOW", dates = "DATES", conc.units="mg/L", station="Illinois River at Marseilles, Ill.") print(app1.lr)
Computes the basis for a segmented model. Used primarily in a linear regression or load model.
segment(x, psi)segment(x, psi)
x |
the data to segment. Missing values are permitted and result corresponding in missing values in output. |
psi |
a numeric vector of the breakpoints. |
A matrix contining a column named X of the data in x and
paired U and P columns for each of the breaks in psi that form the
basis for that segment.
Select the "best" predefined rating-curve (regression) model for river load estimation.
selBestModel(constituent, data, subset, na.action, flow, dates, flow.units = "cfs", conc.units = "", load.units = "kg", time.step = "day", station = "", criterion = c("AIC", "SPCC", "AICc"))selBestModel(constituent, data, subset, na.action, flow, dates, flow.units = "cfs", conc.units = "", load.units = "kg", time.step = "day", station = "", criterion = c("AIC", "SPCC", "AICc"))
constituent |
a character string giving the name of the response variable for which loads are to be computed. |
data |
the data to search for the variables in |
subset |
an expression to select a subset of the data. |
na.action |
what to do with missing values. |
flow |
character string indicating the name of the flow column. |
dates |
character string indicating the name of the date column. |
flow.units |
character string describing the flow unit. |
conc.units |
character string describing the concentration unit. |
load.units |
character string describing the load unit. |
time.step |
character string describing the time step of the calibration data. |
station |
character string description of the station. |
criterion |
the criterion to use for subset selection, must be one of "AIC," "SPCC," or "AICc." |
An object of class "loadReg."
will need some.
# Use the data from application 1 in the vignettes data(app1.calib) app1.lr <- selBestModel("Phosphorus", data = app1.calib, flow = "FLOW", dates = "DATES", conc.units="mg/L", station="Illinois River at Marseilles, Ill.") # Extract the fitted values print(app1.lr)# Use the data from application 1 in the vignettes data(app1.calib) app1.lr <- selBestModel("Phosphorus", data = app1.calib, flow = "FLOW", dates = "DATES", conc.units="mg/L", station="Illinois River at Marseilles, Ill.") # Extract the fitted values print(app1.lr)
Select the "best" subset of a user-defined rating-curve (regression) model for rver load estimation.
selBestSubset(formula, min.formula = ~1, data, subset, na.action, flow, dates, flow.units = "cfs", conc.units = "", load.units = "kg", time.step = "day", station = "", criterion = c("AIC", "SPCC", "AICc"))selBestSubset(formula, min.formula = ~1, data, subset, na.action, flow, dates, flow.units = "cfs", conc.units = "", load.units = "kg", time.step = "day", station = "", criterion = c("AIC", "SPCC", "AICc"))
formula |
a formula describing the regression model. See |
min.formula |
a formula containing the minimum variables to use in the final model. The default is to only force the intercept term, whihc will normally be acceptable. In some rare cases, the log of flow may be needed. |
data |
the data to search for the variables in |
subset |
an expression to select a subset of the data. |
na.action |
what to do with missing values. |
flow |
character string indicating the name of the flow column. |
dates |
character string indicating the name of the date column. |
flow.units |
character string describing the flow unit. |
conc.units |
character string describing the concentration unit. |
load.units |
character string describing the load unit. |
time.step |
character string describing the time step of the calibration data. |
station |
character string description of the station. |
criterion |
the criterion to use for subset selection, must be one of "AIC," "SPCC," or "AICc." |
An object of class "loadReg."
The printed output of the model inlcudes the anova component from
step. That table summarizes the step wise selection and the criterion used
for each step. The statistics possibly represent a smaller sample size than used
for the final model because the step function requires a data set with no
missing values. If missing values are found a warning is printed.
Internal support function to extract the model matrix for one of the 9 predefined models in LOADEST.
setXLDat(data, flow, dates, Qadj, Tadj, model.no)setXLDat(data, flow, dates, Qadj, Tadj, model.no)
data |
the dataset containing dates and flow columns. |
flow |
character string indicating the name of the flow column. |
dates |
character string indicating the name of the date column. |
Qadj |
the centering value for flow. |
Tadj |
the centering value for decimal time. |
model.no |
the model number, must be in the range 1-9. |
The model matrix corresponding to the selected model number.
Computes the variance inflation factor (Helsel and Hirsch, 2002) for each variable in a load regression.
## S3 method for class 'loadReg' vif(model, ...)## S3 method for class 'loadReg' vif(model, ...)
model |
an object of class "loadReg"—output from |
... |
further arguments passed to or from other methods. |
further arguments passed to or from other methods.
# From application 1 in the vignettes data(app1.calib) app1.lr <- loadReg(Phosphorus ~ model(1), data = app1.calib, flow = "FLOW", dates = "DATES", conc.units="mg/L", station="Illinois River at Marseilles, Ill.") vif(app1.lr, app1.calib)# From application 1 in the vignettes data(app1.calib) app1.lr <- loadReg(Phosphorus ~ model(1), data = app1.calib, flow = "FLOW", dates = "DATES", conc.units="mg/L", station="Illinois River at Marseilles, Ill.") vif(app1.lr, app1.calib)