Skip to content

Commit 51183f5

Browse files
dlorenz-usgsdlorenz-usgs
committed
new and revised vignettes, print bug fix
1 parent 4afdcc0 commit 51183f5

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

64 files changed

+1302
-78
lines changed

‎DESCRIPTION

+2-2
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,8 @@
11
Package: rloadest
22
Type: Package
33
Title: River Load Estimation
4-
Version: 0.4.2
5-
Date: 2015-07-20
4+
Version: 0.4.3
5+
Date: 2015-12-03
66
Author: Dave Lorenz, Rob Runkel, Laura De Cicco
77
Maintainer: Dave Lorenz <[email protected]>
88
Description: Collection of functions to make constituent load estimations

‎NAMESPACE

+3-1
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
# Generated by roxygen2 (4.1.0): do not edit by hand
1+
# Generated by roxygen2 (4.1.1): do not edit by hand
22

33
S3method(censoring,Surv)
44
S3method(coef,loadReg)
@@ -8,6 +8,7 @@ S3method(makepredictcall,center)
88
S3method(mean,character)
99
S3method(mean,factor)
1010
S3method(plot,loadReg)
11+
S3method(print,jackStats)
1112
S3method(print,loadReg)
1213
S3method(residuals,loadReg)
1314
S3method(rmse,loadReg)
@@ -16,6 +17,7 @@ export(AICc)
1617
export(c2load)
1718
export(center)
1819
export(dailyAg)
20+
export(jackStats)
1921
export(loadConvFactor)
2022
export(loadReg)
2123
export(loadReport)

‎R/center.R

+1-1
Original file line numberDiff line numberDiff line change
@@ -9,7 +9,7 @@
99
#' @note The centering value by default is computed by the method described
1010
#'in Cohn and others (1992).
1111
#' @return The centered value of \code{x}.
12-
#' @seealso \code{\link[USGSwsBase]{quadratic}}, \code{\link{scale}}
12+
#' @seealso \code{\link[smwrBase]{quadratic}}, \code{\link{scale}}
1313
#' @references Cohn, T.A., Caulder, D.L., Gilroy, E.J., Zynjuk, L.D., and
1414
#'Summers, R.M., 1992, The validity of a simple statistical model for
1515
#'estimating fluvial constituent loads---An empirical study involving nutrient

‎R/jackStats.R

+142
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,142 @@
1+
#' Jackknife Statistics
2+
#'
3+
#' Compute selected jackknife statistics for a rating-curve load-estimation model.
4+
#'
5+
#' @param fit an object of class "loadReg"---output from \code{loadReg}. Can also
6+
#'be an object of class "censReg."
7+
#' @param which a character string indicating the "load" or
8+
#'"concentration" model for an object of class "loadReg" or "censReg" for
9+
#'an object of class "censReg."
10+
#' @return An object of class "jackStats" containing these components:
11+
#'coef, the table of coefficient estimates, the jackknife bias and standard errors\cr
12+
#'coefficients, the jackknifed coefficients\cr
13+
#'pctcens, the percentage of left-censored values. \cr
14+
#'The PRESS statistic and individual jackknife differences are also returned
15+
#'when the percentage of censoring is 0.
16+
#' @note The \code{jackStats} function can only be used when the analysis is AMLE.
17+
#'
18+
#'Abdi and Williams (2010) describe the jackknife as refering to two related techniques: the first
19+
#'estimates the parameters, their bias and standard errors and the second evaluates the
20+
#'predictive performance of the model. The second technique is the PRESS statistic (Helsel
21+
#'and Hirsch, 2002), but can only be used on uncensored data; it is computed by \code{jackStats}
22+
#'when no data are censored. The first technique can be used to assess the coefficients of the
23+
#'regression---the bias should be small and the jackknife standard errors should not be much
24+
#'different from the standard errors reported for the regression. Efron and Tibshirani (1993)
25+
#'suggest that the bias is small if the relative bias (biuas divided by the jackknife standard
26+
#'error) is less than 0.25.
27+
#' @seealso \code{\link{loadReg}}
28+
#' @keywords utilities
29+
#' @references
30+
#' Abdi, H. and Williams, L.J., 2010, Jackknife, in encyclopedia of research design,
31+
#'Salkind, N.J., editor: Thousand Oaks, Calif., SAGE Publications, 1719 p.
32+
#'
33+
#'Efron, B. and Tibshirani, R.J., 1993, An introduction to the bootstrap: Boca Raton,
34+
#'Fla., Chapman and Hall/CRC, 436 p.
35+
#'
36+
#' Helsel, D.R. and Hirsch, R.M., 2002, Statistical methods in water resources:
37+
#'U.S. Geological Survey Techniques of Water-Resources Investigations, book 4,
38+
#'chap. A3, 522 p.
39+
#'Salkind,
40+
#' @examples
41+
#'# From application 1 in the vignettes
42+
#'data(app1.calib)
43+
#'app1.lr <- loadReg(Phosphorus ~ model(1), data = app1.calib,
44+
#' flow = "FLOW", dates = "DATES", conc.units="mg/L",
45+
#' station="Illinois River at Marseilles, Ill.")
46+
#'jackStats(app1.lr)
47+
#' @export
48+
jackStats <- function(fit, which="load") {
49+
## Compute some stats
50+
which <- match.arg(which, c("load", "concentration", "censReg"))
51+
if(which == "load") {
52+
# Verify AMLE
53+
if(fit$method != "AMLE") {
54+
stop("The analysis method must be AMLE")
55+
}
56+
# initial stuff
57+
NPAR <- fit$lfit$NPAR
58+
NOBS <- fit$lfit$NOBSC
59+
# Get the repsonse, X, and Yeff
60+
Y <- as.lcens(exp(fit$lfit$YLCAL), exp(fit$lfit$YD), fit$lfit$CENSFLAG)
61+
X <- fit$lfit$XLCAL
62+
Yeff <- fit$lfit$YLCAL
63+
## The code below computes an effective value for left-censored values
64+
## by computing the expected value from the prediction
65+
## provided if the method is ever published. No plans for now
66+
#Yeff[fit$lfit$CENSFLAG] <- fit$lfit$XLCAL[fit$lfit$CENSFLAG, ,drop=FALSE] %*%
67+
# fit$lfit$PARAML[seq(NPAR)] + fit$lfit$RESID[fit$lfit$CENSFLAG]
68+
#
69+
# Other Info
70+
dist <- "lognormal"
71+
parms <- fit$lfit$PARAML[seq(fit$lfit$NPAR)]
72+
parnames <- colnames(fit$lfit$XLCAL)
73+
} else if(which == "concentration") {
74+
# Verify AMLE
75+
if(fit$method != "AMLE") {
76+
stop("The analysis method must be AMLE")
77+
}
78+
# initial stuff
79+
NPAR <- fit$cfit$NPAR
80+
NOBS <- fit$cfit$NOBSC
81+
# Get the repsonse, X, and Yeff
82+
Y <- as.lcens(exp(fit$cfit$YLCAL), exp(fit$cfit$YD), fit$cfit$CENSFLAG)
83+
X <- fit$cfit$XLCAL
84+
Yeff <- fit$cfit$YLCAL
85+
## See comment above
86+
#Yeff[fit$cfit$CENSFLAG] <- fit$cfit$XLCAL[fit$cfit$CENSFLAG, ,drop=FALSE] %*%
87+
# fit$cfit$PARAML[seq(NPAR)] + fit$cfit$RESID[fit$cfit$CENSFLAG]
88+
#
89+
# Other Info
90+
dist <- "lognormal"
91+
parms <- fit$cfit$PARAML[seq(fit$cfit$NPAR)]
92+
parnames <- colnames(fit$cfit$XLCAL)
93+
} else { # must be censReg
94+
# Verify AMLE
95+
if(fit$method != "AMLE") {
96+
stop("The analysis method must be AMLE")
97+
}
98+
# initial stuff
99+
NPAR <- fit$NPAR
100+
NOBS <- fit$NOBSC
101+
dist <- fit$dist
102+
# Get the repsonse, X, and Yeff
103+
if(dist == "lognormal") {
104+
Y <- as.lcens(exp(fit$YLCAL), exp(fit$YD), fit$CENSFLAG)
105+
} else {
106+
Y <- as.lcens(fit$YLCAL, fit$YD, fit$CENSFLAG)
107+
}
108+
X <- fit$XLCAL
109+
Yeff <- fit$YLCAL
110+
## See comment above
111+
#Yeff[fit$CENSFLAG] <- fit$XLCAL[fit$CENSFLAG, ,drop=FALSE] %*%
112+
# fit$PARAML[seq(NPAR)] + fit$RESID[fit$CENSFLAG]
113+
#
114+
# Other Info
115+
parms <- fit$PARAML[seq(fit$NPAR)]
116+
parnames <- colnames(fit$XLCAL)
117+
}
118+
# do it
119+
# set up res and coeff storage
120+
pre <- numeric(NOBS)
121+
coeff <- matrix(0, nrow=NOBS, ncol=NPAR)
122+
for(i in seq(NOBS)) {
123+
tmp <- censReg_AMLE.fit(Y[-i], X[-i,], dist)
124+
coeff[i,] <- tmp$PARAML[-(NPAR + 1L)]
125+
pre[i] <- Yeff[i] - coeff[i,, drop=FALSE] %*% t(X[i,,drop=FALSE])
126+
}
127+
# Compute the jackknife bias and variance of coeffs
128+
out <- (NOBS - 1)*(rep(parms, each=NOBS) - coeff)
129+
bias <- -1/NOBS*colSums(out)
130+
var <- 1/(NOBS*(NOBS-1)) * (colSums(out^2) - NOBS * bias^2)
131+
coef <- cbind(est=parms, bias=bias, stderr=sqrt(var))
132+
rownames(coef) <- parnames
133+
retval <- list(coef=coef, coefficients=coeff,
134+
pctcens=pctCens(Y))
135+
# Include press if no censoring
136+
if(retval$pctcens == 0) {
137+
retval$press <- sum(pre^2)
138+
retval$pre <- pre
139+
}
140+
class(retval) <- "jackStats"
141+
return(retval)
142+
}

‎R/plot.loadReg.R

+1-1
Original file line numberDiff line numberDiff line change
@@ -22,7 +22,7 @@
2222
#' @param span the span to use for the loess smooth. Set to 0 to suppress.
2323
#' @param \dots further arguments passed to or from other methods.
2424
#' @return The object \code{x} is returned invisibly.
25-
#' @note This plotting function uses the core routines in the \code{USGSwsGraphs}
25+
#' @note This plotting function uses the core routines in the \code{smwrGraphs}
2626
#'package. It requires a graphics page that is set up from the functions
2727
#'in that package (\code{setpage} or \code{setPDF}) if \code{set.up} is
2828
#'\code{FALSE}. The graphs that are produced by this function are based

‎R/predConc.R

+6-4
Original file line numberDiff line numberDiff line change
@@ -190,9 +190,10 @@ predConc <- function(fit, newdata, by="day",
190190
names(retval)[6L:7L] <- Names
191191
} else if(by == "day") {
192192
## One or more obs per day
193+
## KDate are the dates of the days
193194
## KDays are the indexes to days
194195
## Kin are the indexes to the good model inputs
195-
## Kdy are the indexes to the days
196+
## Kdy are the indexes to the days of model inputs
196197
## KinAll are the indexes to all days
197198
##
198199
## Preserve flow for later
@@ -210,11 +211,12 @@ predConc <- function(fit, newdata, by="day",
210211
Kin <- Kin[is.finite(rowSums(model.inp))]
211212
KDate <- as.Date(as.POSIXlt(newdata[[dates]]))
212213
Kdy <- as.integer(KDate)
213-
KDate <- unique(KDate)
214-
Kdy <- Kdy - Kdy[1L] + 1L # make relative to first day (Index)
215-
if(length(unique(table(Kdy))) > 1L && !allow.incomplete) {
214+
Ktbl <- table(Kdy)
215+
if(length(unique(Ktbl)) > 1L && !allow.incomplete) {
216216
warning("Variable observations per day, either set the allow.incomplete argument to TRUE or use the resampleUVdata function to construct a uniform series")
217217
}
218+
Kdy <- rep(seq(1, length(Ktbl)), Ktbl) # Create the index to days in input
219+
KDate <- unique(KDate)
218220
KinAll <- unique(Kdy)
219221
## Make it daily flow, Flow0 indicates a partial 0 flow
220222
Flow0 <- tapply(Flow, Kdy, min)

‎R/predLoad.R

+20-18
Original file line numberDiff line numberDiff line change
@@ -215,9 +215,10 @@ predLoad <- function(fit, newdata, load.units=fit$load.units, by="total",
215215
names(retval)[6L:7L] <- Names
216216
} else if(by == "day") {
217217
## One or more obs per day
218+
## KDate are the dates of the days
218219
## KDays are the indexes to days
219220
## Kin are the indexes to the good model inputs
220-
## Kdy are the indexes to the days
221+
## Kdy are the indexes to the days of model inputs
221222
## KinAll are the indexes to all days
222223
##
223224
## Preserve flow for later
@@ -235,11 +236,12 @@ predLoad <- function(fit, newdata, load.units=fit$load.units, by="total",
235236
Kin <- Kin[is.finite(rowSums(model.inp))]
236237
KDate <- as.Date(as.POSIXlt(newdata[[dates]]))
237238
Kdy <- as.integer(KDate)
238-
KDate <- unique(KDate)
239-
Kdy <- Kdy - Kdy[1L] + 1L # make relative to first day (Index)
240-
if(length(unique(table(Kdy))) > 1L && !allow.incomplete) {
239+
Ktbl <- table(Kdy)
240+
if(length(unique(Ktbl)) > 1L && !allow.incomplete) {
241241
warning("Variable observations per day, either set the allow.incomplete argument to TRUE or use the resampleUVdata function to construct a uniform series")
242242
}
243+
Kdy <- rep(seq(1, length(Ktbl)), Ktbl) # Create the index to days in input
244+
KDate <- unique(KDate)
243245
KinAll <- unique(Kdy)
244246
## Make it daily flow, Flow0 indicates a partial 0 flow
245247
Flow0 <- tapply(Flow, Kdy, min)
@@ -392,17 +394,17 @@ predLoad <- function(fit, newdata, load.units=fit$load.units, by="total",
392394
if(any(drop)) { # Any missing values?
393395
OK <- FALSE
394396
retby <- data.frame(Period="period", Ndays=NDays,
395-
Flux=NA_real_, Std.Err=NA_real_,
396-
SEP=NA_real_, L95=NA_real_, U95=NA_real_)
397+
Flux=NA_real_, Std.Err=NA_real_,
398+
SEP=NA_real_, L95=NA_real_, U95=NA_real_)
397399
} else {
398400
## Check for complete periods
399401
Period <- as.character(newdata[[by]][period[1L]])
400402
if(!is.null(targN <- checkDays[[Period]])) {
401403
if(targN != NDays) {
402404
OK <- FALSE
403405
retby <- data.frame(Period="period", Ndays=NDays,
404-
Flux=NA_real_, Std.Err=NA_real_,
405-
SEP=NA_real_, L95=NA_real_, U95=NA_real_)
406+
Flux=NA_real_, Std.Err=NA_real_,
407+
SEP=NA_real_, L95=NA_real_, U95=NA_real_)
406408

407409
}
408410
}
@@ -440,16 +442,16 @@ predLoad <- function(fit, newdata, load.units=fit$load.units, by="total",
440442
if(out.data$IERR > 0) {
441443
warning(" *** Error (",out.data$IERR,") occurred in processing data. ***\n",sep="")
442444
retby <- data.frame(Period="period", Ndays=NA_integer_,
443-
Flux=NA_real_, Std.Err=NA_real_,
444-
SEP=NA_real_, L95=NA_real_, U95=NA_real_)
445+
Flux=NA_real_, Std.Err=NA_real_,
446+
SEP=NA_real_, L95=NA_real_, U95=NA_real_)
445447
} else {
446448
## OK, we've had a successful run, pack the data into a data frame
447449
retby <- data.frame(Period="period", Ndays=NDays,
448-
Flux=out.data$load * lcf,
449-
Std.Err=sqrt(out.data$loadvar) * lcf,
450-
SEP=out.data$loadsep * lcf,
451-
L95=out.data$loadlow * lcf,
452-
U95=out.data$loadup * lcf)
450+
Flux=out.data$load * lcf,
451+
Std.Err=sqrt(out.data$loadvar) * lcf,
452+
SEP=out.data$loadsep * lcf,
453+
L95=out.data$loadlow * lcf,
454+
U95=out.data$loadup * lcf)
453455
}
454456
}
455457
return(retby)
@@ -476,14 +478,14 @@ predLoad <- function(fit, newdata, load.units=fit$load.units, by="total",
476478
"calibration data set streamflow. Load estimates require extrapolation.\n\n",
477479
sep="")
478480
} else if(Qsum[2L, 1L] < Qsum[1L, 1L]) {
479-
cat("WARNING: The minimum estimation data set steamflow exceeds the minimum\n",
481+
cat("WARNING: The minimum estimation data set steamflow exceeds the minimum\n",
480482
"calibration data set streamflow. Load estimates require extrapolation.\n\n",
481483
sep="")
482-
} else {
484+
} else {
483485
cat("The maximum estimation data set streamflow does not exceed the maximum\n",
484486
"calibration data set streamflow. No extrapolation is required.\n\n",
485487
sep="")
486-
}
488+
}
487489
if(!(by %in% c("day", "unit"))) {
488490
cat("\n-------------------------------------------------------------\n",
489491
"Constituent Output File Part IIb: Estimation (Load Estimates)\n",

‎R/print.jackStats.R

+31
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,31 @@
1+
#' Print Results
2+
#'
3+
#' Print the results of a jackknife analysis of a censored regression.
4+
#'
5+
#' @param x an object of class "jackStats"---output from \code{jackStats}.
6+
#' @param digits the number of significant digits to print.
7+
#' @param \dots further arguments passed to or from other methods.
8+
#' @return The object \code{x} is returned invisibly.
9+
#' @note The printed output includes the original original estimate,
10+
#'the jackknife bias and standard error and the relative bias for
11+
#'each parameter in the regression model.
12+
#'
13+
#' @seealso \code{\link{loadReg}}
14+
#' @keywords utilities
15+
#' @export
16+
#' @method print jackStats
17+
print.jackStats <- function(x, digits=4, ...) {
18+
##
19+
if(x$pctcens == 0) {
20+
cat("jackknife estimates:\n\n",
21+
"PRESS: ", signif(x$press, digits), "\nCoefficients:\n", sep="")
22+
} else {
23+
cat("jackknife estimates:\n\nCoefficients:\n", sep="")
24+
}
25+
coef <- x$coef
26+
# compute the relative bias
27+
coef <- cbind(coef, abs(coef[,2L]/coef[,3L]))
28+
colnames(coef) <- c("Estimate", "Bias", "Std. Error", "Rel. Bias")
29+
print(coef, digits=digits)
30+
invisible(x)
31+
}

‎R/print.loadReg.R

+21-3
Original file line numberDiff line numberDiff line change
@@ -76,8 +76,15 @@ print.loadReg <- function(x, digits=4, brief=TRUE, load.only=brief, ...) {
7676
" Constituent Output File Part Ia: Calibration (Load Regression)\n",
7777
"----------------------------------------------------------------------\n\n",
7878
sep="")
79-
# CENSFLAG is logical for AMLE, integer for MLE, this protects
80-
Ncen <- sum(x$lfit$CENSFLAG != 0)
79+
# CENSFLAG is logical for AMLE, integer for MLE, this protects,
80+
# but does not work for other than left-censored values.
81+
if(x$method == "AMLE") {
82+
Ncen <- sum(x$lfit$CENSFLAG != 0)
83+
} else {
84+
Centbl <- table(x$lfit$survreg$y[,3L])
85+
Ncen <- x$lfit$NOBSC - Centbl["1"]
86+
names(Centbl) <- c("0"="right", "1"="none", "2"="left", "3"="interval")[names(Centbl)]
87+
}
8188
cat(" Number of Observations: ", x$lfit$NOBSC, "\n",
8289
"Number of Uncensored Observations: ",
8390
x$lfit$NOBSC - Ncen, "\n",
@@ -88,6 +95,11 @@ print.loadReg <- function(x, digits=4, brief=TRUE, load.only=brief, ...) {
8895
" Period of record: ",
8996
as.character(x$PoR[1L]), " to ", as.character(x$PoR[2L]),
9097
"\n\n", sep="")
98+
if(x$method == "MLE") {
99+
cat("Censoring of data:\n")
100+
print(Centbl)
101+
cat("\n\n")
102+
}
91103
if(!brief || nrow(x$model.eval) > 1) { # Capture best model selection
92104
cat("Model Evaluation Criteria Based on ", x$method, " Results\n",
93105
"-----------------------------------------------\n\n", sep="")
@@ -286,7 +298,13 @@ print.loadReg <- function(x, digits=4, brief=TRUE, load.only=brief, ...) {
286298
" Constituent Output File Part Ia: Calibration (Concentration Regression)\n",
287299
"--------------------------------------------------------------------------\n\n",
288300
sep="")
289-
Ncen <- sum(x$cfit$CENSFLAG)
301+
if(x$method == "AMLE") {
302+
Ncen <- sum(x$cfit$CENSFLAG != 0)
303+
} else {
304+
Centbl <- table(x$cfit$survreg$y[,3L])
305+
Ncen <- x$cfit$NOBSC - Centbl["1"]
306+
names(Centbl) <- c("0"="right", "1"="none", "2"="left", "3"="interval")[names(Centbl)]
307+
}
290308
if(!brief) { # Capture best model selection
291309
cat("Model # ", x$model.no, " selected\n\n", sep="")
292310
}

0 commit comments

Comments
 (0)