Skip to content

Commit

Permalink
Merge pull request #63 from ldecicco-USGS/master
Browse files Browse the repository at this point in the history
Randomized residuals.
  • Loading branch information
ldecicco-USGS committed Sep 9, 2015
2 parents f97ee76 + bd7a88a commit b13a80c
Show file tree
Hide file tree
Showing 31 changed files with 516 additions and 195 deletions.
7 changes: 4 additions & 3 deletions DESCRIPTION
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
Package: EGRET
Type: Package
Title: Exploration and Graphics for RivEr Trends (EGRET)
Version: 2.3.1
Version: 2.4.1
Authors@R: c( person("Robert", "Hirsch", role = c("aut"),
email = "rhirsch@usgs.gov"),
person("Laura", "DeCicco", role = c("aut","cre"),
Expand All @@ -10,7 +10,7 @@ Description: Statistics and graphics for streamflow history,
water quality trends, and the statistical modeling algorithm: Weighted
Regressions on Time, Discharge, and Season (WRTDS).
License: CC0
Date: 2015-07-23
Date: 2015-09-08
Depends:
R (>= 3.0)
Imports:
Expand All @@ -22,7 +22,8 @@ Imports:
utils,
graphics,
stats,
grDevices
grDevices,
truncnorm
Suggests:
xtable,
knitr,
Expand Down
2 changes: 2 additions & 0 deletions NAMESPACE
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,7 @@ export(is.egret)
export(logPretty1)
export(logPretty3)
export(makeAnnualSeries)
export(makeAugmentedSample)
export(mergeReport)
export(modelEstimation)
export(monthInfo)
Expand Down Expand Up @@ -115,3 +116,4 @@ import(survival)
import(utils)
importFrom(fields,interp.surface)
importFrom(lubridate,decimal_date)
importFrom(truncnorm,rtruncnorm)
4 changes: 4 additions & 0 deletions NEWS
Original file line number Diff line number Diff line change
@@ -1,3 +1,7 @@
EGRET 2.3.1
===========
* Removed unnecessary disclaimer.

EGRET 2.2.0
===========
* Minor bug fixes
Expand Down
27 changes: 19 additions & 8 deletions R/boxConcThree.R
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@
#' @param cex numerical value giving the amount by which plotting symbols should be magnified
#' @param tinyPlot logical variable, if TRUE plot is designed to be plotted small as part of a multi-plot figure, default is FALSE.
#' @param customPar logical defaults to FALSE. If TRUE, par() should be set by user before calling this function
#' @param rResid logical option to plot randomized residuals.
#' @param \dots arbitrary graphical parameters that will be passed to genericEGRETDotPlot function (see ?par for options)
#' @keywords graphics water-quality statistics
#' @seealso \code{\link[graphics]{boxplot}}
Expand All @@ -31,7 +32,7 @@
#' boxConcThree(eList)
boxConcThree<-function (eList, tinyPlot=FALSE,
printTitle = TRUE, moreTitle = "WRTDS",customPar=FALSE,
font.main=2,cex=0.8,cex.main = 1.1, cex.axis = 1.1,...){
font.main=2,cex=0.8,cex.main = 1.1, cex.axis = 1.1,rResid=FALSE,...){

localINFO <- getInfo(eList)
localSample <- getSample(eList)
Expand All @@ -56,7 +57,6 @@ boxConcThree<-function (eList, tinyPlot=FALSE,
index2 <- rep(2, nS)
index3 <- rep(3, nD)
index <- c(index1, index2,index3)
concV <- c(localSample$ConcAve,localSample$ConcHat,localDaily$ConcDay)

plotTitle <- if (printTitle) {
paste(localINFO$shortName, ",", localINFO$paramShortName,
Expand All @@ -65,10 +65,6 @@ boxConcThree<-function (eList, tinyPlot=FALSE,
""
}

yMax<-max(concV,na.rm=TRUE)
yTicks<-yPretty(yMax)
yTop<-yTicks[length(yTicks)]

if (tinyPlot) {
yLab <- paste("Conc. (",localINFO$param.units,")",sep="")
if (!customPar) par(mar=c(4,5,1,0.1),tcl=0.5,cex.lab=cex.axis)
Expand All @@ -83,19 +79,34 @@ boxConcThree<-function (eList, tinyPlot=FALSE,
name3 <- "All day\nestimates"
groupNames <- c(name1,name2,name3)

if(!rResid){
concV <- c(localSample$ConcAve,localSample$ConcHat,localDaily$ConcDay)

} else {
if(!("rObserved" %in% names(localSample))){
eList <- makeAugmentedSample(eList)
localSample <- eList$Sample
}
concV <- c(localSample$rObserved,localSample$ConcHat,localDaily$ConcDay)
}

yMax<-max(concV,na.rm=TRUE)
yTicks<-yPretty(yMax)
yTop<-yTicks[length(yTicks)]

boxplot(concV ~ index,varwidth=TRUE,
names=groupNames,xlab="",ylab=yLab,
ylim=c(0,yTop),axes=FALSE,
main=plotTitle,font.main=font.main,cex=cex,
cex.main=cex.main,
las=1,yaxs="i",
...)
...)

axis(1,tcl=0.5,at=c(1,2,3),labels=groupNames,cex.axis=cex.axis*0.5454)
axis(2,tcl=0.5,las=1,at=yTicks,cex.axis=cex.axis)
axis(3,tcl=0.5,at=c(1,2,3),labels=FALSE)
axis(4,tcl=0.5,at=yTicks,labels=FALSE)
box()
if (!tinyPlot) mtext(title2,side=3,line=-1.5)

invisible(eList)
}
26 changes: 19 additions & 7 deletions R/boxResidMonth.R
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,7 @@
#' @param font.main font to be used for plot main titles
#' @param customPar logical defaults to FALSE. If TRUE, par() should be set by user before calling this function
#' @param las numeric in {0,1,2,3}; the style of axis labels
#' @param rResid logical option to plot randomized residuals.
#' @param \dots arbitrary graphical parameters that will be passed to genericEGRETDotPlot function (see ?par for options)
#' @keywords graphics water-quality statistics
#' @seealso \code{\link[graphics]{boxplot}}
Expand All @@ -35,7 +36,7 @@
#' boxResidMonth(eList)
boxResidMonth<-function(eList, stdResid = FALSE, las=1,
printTitle = TRUE, cex=0.8, cex.axis=1.1, cex.main=1.1,
font.main=2, tinyPlot=FALSE, customPar=FALSE,...) {
font.main=2, tinyPlot=FALSE, customPar=FALSE,rResid=FALSE,...) {

localINFO <- getInfo(eList)
localSample <- getSample(eList)
Expand Down Expand Up @@ -63,18 +64,28 @@ boxResidMonth<-function(eList, stdResid = FALSE, las=1,
}

plotTitle<-if(printTitle) paste(localINFO$shortName,"\n",localINFO$paramShortName,"\nBoxplots of residuals by month") else ""
resid<-log(localSample$ConcAve) - localSample$yHat
resid<-if(stdResid) resid/localSample$SE else resid

if(!rResid){
resid<-log(localSample$ConcAve) - localSample$yHat
} else {
if(!("rResid" %in% names(localSample))){
eList <- makeAugmentedSample(eList)
localSample <- eList$Sample
}
resid <- localSample$rResid
}

if(stdResid) {
resid<- resid/localSample$SE
}

singleMonthList <- sapply(c(1:12),function(x){monthInfo[[x]]@monthAbbrev})

namesListFactor <- factor(singleMonthList, levels=singleMonthList)
monthList <- as.character(apply(localSample, 1, function(x) monthInfo[[as.numeric(x[["Month"]])]]@monthAbbrev))
monthList <- factor(monthList, namesListFactor)

tempDF <- data.frame(month=monthList, resid=resid)

boxplot(tempDF$resid ~ tempDF$month,

boxplot(resid ~ monthList,
varwidth=TRUE,
xlab="Month",ylab=yLab,
main=plotTitle,
Expand All @@ -86,4 +97,5 @@ boxResidMonth<-function(eList, stdResid = FALSE, las=1,
...)
abline(h=0)
if (!tinyPlot) mtext(title2,side=3,line=-1.5)
invisible(eList)
}
8 changes: 4 additions & 4 deletions R/censoredSegments.R
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@
#' @export
#' @examples
#' x <- c(1,2,3,4,5,6)
#' y <- c(1,3,4,3.3,4.4,2)
#' y <- c(1,3,4,3.3,4.4,7)
#' xlim <- c(min(x)*.75,max(x)*1.25)
#' ylim <- c(0,1.25*max(y))
#' xlab <- "Date"
Expand All @@ -28,9 +28,9 @@
#' plotTitle="Test"
#' )
#' yBottom <- 0
#' yLow <- c(NA,3,4,3.3,4.4,2)
#' yHigh <- c(1,3,4,3.3,4.4,2)
#' Uncen <- c(0,1,1,1,1,1)
#' yLow <- c(NA,3,4,3.3,4,7)
#' yHigh <- c(1,3,4,3.3,5,NA)
#' Uncen <- c(0,1,1,1,0,0)
#' censoredSegments(yBottom=yBottom,yLow=yLow,yHigh=yHigh,x=x,Uncen=Uncen)
censoredSegments <- function(yBottom,yLow,yHigh,x,Uncen,col="black",lwd=1){
yLowVal<-ifelse(is.na(yLow),yBottom,yLow) #yLow would be NA if "simple" censored....so giving it a value here
Expand Down
21 changes: 12 additions & 9 deletions R/fluxBiasEight.R
Original file line number Diff line number Diff line change
Expand Up @@ -21,12 +21,14 @@
#' @param cex.axis magnification to be used for axis annotation relative to the current setting of cex
#' @param col color of points on plot, see ?par 'Color Specification'
#' @param lwd number line width
#' @param rResid logical option to plot censored residuals as segments, or randomized points.
#' @param \dots arbitrary graphical parameters that will be passed to genericEGRETDotPlot function (see ?par for options)
#' @keywords graphics water-quality statistics
#' @export
#' @examples
#' eList <- Choptank_eList
#' fluxBiasMulti(eList)
#' fluxBiasMulti(eList, rResid=TRUE)
#' # Water year:
#' \dontrun{
#' pdf("fluxBiasMulti.pdf", height=9, width=8)
Expand All @@ -40,7 +42,7 @@
#' }
fluxBiasMulti<-function (eList, qUnit = 2, fluxUnit = 3, moreTitle = "WRTDS",
cex = 0.7, cex.axis = 1.1,cex.main=1.1,
col="black", lwd=1,...){
col="black", lwd=1,rResid=FALSE,...){

localINFO <- getInfo(eList)
localSample <- getSample(eList)
Expand Down Expand Up @@ -69,29 +71,29 @@ fluxBiasMulti<-function (eList, qUnit = 2, fluxUnit = 3, moreTitle = "WRTDS",
title2<-if(paLong==12) "" else setSeasonLabelByUser(paStartInput=paStart,paLongInput=paLong)

par(oma = c(0, 10, 4, 10),mfrow=c(4,2))
plotResidPred(eList,
eList <- plotResidPred(eList,
stdResid = FALSE, tinyPlot=TRUE, printTitle = FALSE,cex=cex,
cex.axis = cex.axis, col=col,lwd=lwd,...)
cex.axis = cex.axis, col=col,lwd=lwd, rResid=rResid,...)
plotResidQ(eList,
qUnit, tinyPlot = TRUE, printTitle = FALSE,cex=cex,
cex.axis = cex.axis, col=col,lwd=lwd,...)
cex.axis = cex.axis, col=col,lwd=lwd, rResid=rResid,...)
plotResidTime(eList,
printTitle = FALSE, tinyPlot=TRUE,cex=cex,
cex.axis = cex.axis, col=col,lwd=lwd,...)
cex.axis = cex.axis, col=col,lwd=lwd, rResid=rResid,...)
boxResidMonth(eList,
printTitle = FALSE, tinyPlot=TRUE,cex=cex,
cex.axis = cex.axis,lwd=lwd,...)
cex.axis = cex.axis,lwd=lwd, rResid=rResid,...)
boxConcThree(eList,
localINFO = localINFO, printTitle=FALSE, tinyPlot=TRUE,cex=cex,
cex.axis = cex.axis, lwd=lwd,...)
cex.axis = cex.axis, lwd=lwd, rResid=rResid,...)
plotConcPred(eList, printTitle=FALSE,
tinyPlot=TRUE,cex=cex,
cex.axis = cex.axis, col=col,lwd=lwd,...)
cex.axis = cex.axis, col=col,lwd=lwd, rResid=rResid,...)
boxQTwice(eList, printTitle = FALSE, qUnit = qUnit,tinyPlot=TRUE,cex=cex,
cex.axis = cex.axis, lwd=lwd,...)
plotFluxPred(eList,
fluxUnit, tinyPlot = TRUE, printTitle = FALSE,cex=cex,
cex.axis = cex.axis, col=col,lwd=lwd,...)
cex.axis = cex.axis, col=col,lwd=lwd, rResid=rResid,...)
fluxBias <- fluxBiasStat(localSample)
fB <- as.numeric(fluxBias[3])
fB <- format(fB, digits = 3)
Expand All @@ -105,5 +107,6 @@ fluxBiasMulti<-function (eList, qUnit = 2, fluxUnit = 3, moreTitle = "WRTDS",
}

par(mfcol = c(1, 1), oma = c(0, 0, 0, 0))
invisible(eList)

}
40 changes: 40 additions & 0 deletions R/makeRandomResiduals.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,40 @@
#' Create Randomized Residuals and Observations
#'
#' This function is used to add two columns to the Sample data frame: rResid and rObserved.
#' rResid is the randomized residual value computed in log concentration units, and rObserved
#' is the randomized 'observed' value of concentration in concentration units.
#'
#' @param eList named list with at least the Daily dataframe
#' @keywords water-quality statistics
#' @return eList named list with modified Daily data frame.
#' @export
#' @importFrom truncnorm rtruncnorm
#' @examples
#' eList <- Choptank_eList
#' eList <- makeAugmentedSample(eList)
makeAugmentedSample <- function(eList){

localSample <- eList$Sample
numSamples<-length(localSample$Uncen)
a <- ifelse(localSample$Uncen==0&!is.na(localSample$ConcLow),log(localSample$ConcLow)-localSample$yHat,-Inf)
b <- ifelse(localSample$Uncen==1,+Inf,log(localSample$ConcHigh) - localSample$yHat)
mean <- ifelse(localSample$Uncen==1,log(localSample$ConcHigh) - localSample$yHat,0)
sd <- ifelse(localSample$Uncen==1,0,localSample$SE)
localSample$rResid<-truncnorm::rtruncnorm(numSamples,a,b,mean,sd)
localSample$rObserved <- exp(localSample$rResid + localSample$yHat)

eList <- as.egret(eList$INFO, eList$Daily, localSample, eList$surfaces)
return(eList)
}

# makeRandomResiduals <- function(eList){
# localSample <- eList$Sample
# numSamples<-length(localSample$Uncen)
# a <- ifelse(localSample$Uncen==0&!is.na(localSample$ConcLow),log(localSample$ConcLow)-localSample$yHat,-Inf)
# b <- ifelse(localSample$Uncen==1,+Inf,log(localSample$ConcHigh) - localSample$yHat)
# mean <- ifelse(localSample$Uncen==1,log(localSample$ConcHigh) - localSample$yHat,0)
# sd <- ifelse(localSample$Uncen==1,0,localSample$SE)
# localSample$rResid<-truncnorm::rtruncnorm(numSamples,a,b,mean,sd)
# eList <- as.egret(eList$INFO, eList$Daily, localSample, eList$surfaces)
# return(eList)
# }
8 changes: 5 additions & 3 deletions R/multiPlotDataOverview.R
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@
#' @param cex.main magnification to be used for main titles relative to the current setting of cex
#' @param logScaleConc logical if TRUE y in concentration graphs plotted in log axis. Default is TRUE.
#' @param logScaleQ logical if TRUE y in streamflow graphs plotted in log axis. Default is TRUE.
#' @param rResid logical option to plot randomized residuals.
#' @keywords graphics water-quality statistics
#' @seealso \code{\link{plotConcQ}}, \code{\link{boxConcMonth}}, \code{\link{plotConcTime}}, \code{\link{boxQTwice}}
#' @export
Expand All @@ -24,7 +25,7 @@
#' eList <- setPA(eList, paStart=6,paLong=3)
#' multiPlotDataOverview(eList, qUnit=1)
multiPlotDataOverview<-function (eList, qUnit = 2,cex.main=1.2,
logScaleConc=TRUE, logScaleQ=TRUE){
logScaleConc=TRUE, logScaleQ=TRUE,rResid=FALSE){

localINFO <- getInfo(eList)

Expand All @@ -39,9 +40,10 @@ multiPlotDataOverview<-function (eList, qUnit = 2,cex.main=1.2,
title2<-if(paLong==12) "" else setSeasonLabelByUser(paStartInput=paStart,paLongInput=paLong)

par(mfcol=c(2,2),oma=c(0,2.4,4.5,2.4),tcl=0.5)
plotConcQ(eList, qUnit = qUnit, tinyPlot = TRUE, printTitle = FALSE,rmSciX=TRUE,logScale=logScaleConc)
plotConcQ(eList, qUnit = qUnit, tinyPlot = TRUE, printTitle = FALSE,
rmSciX=TRUE,logScale=logScaleConc,rResid=rResid)
boxConcMonth(eList, printTitle = FALSE, tinyPlot=TRUE,logScale=logScaleConc)
plotConcTime(eList, printTitle = FALSE, tinyPlot = TRUE,logScale=logScaleConc)
plotConcTime(eList, printTitle = FALSE, tinyPlot = TRUE,logScale=logScaleConc,rResid=rResid)
boxQTwice(eList, printTitle = FALSE, qUnit = qUnit, tinyPlot=TRUE,logScale=logScaleQ)
title<-paste(localINFO$shortName,"\n",localINFO$paramShortName)

Expand Down

0 comments on commit b13a80c

Please sign in to comment.