Merge pull request #9 from ldecicco-USGS/master

Updates to historgrams

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: EGRETci
 Type: Package
 Title: Exploration and Graphics for RivEr Trends (EGRET) Confidence Intervals
-Version: 1.0.1
+Version: 1.0.2
 Authors@R: c( person("Robert", "Hirsch", role = c("aut"),
     email = "[email protected]"),
     person("Laura", "DeCicco", role = c("aut","cre"),

R/plotCIs.R

@@ -319,13 +319,14 @@ ciBands <- function(eList, repAnnualResults, probs=c(0.05,0.95)){
 #' @param eList named list with at least the Daily, Sample, and INFO dataframes. Created from the EGRET package, after running \code{\link[EGRET]{modelEstimation}}.
 #' @param eBoot named list. Returned from \code{\link{wBT}}.
 #' @param caseSetUp data frame. Returned from \code{\link{trendSetUp}}.
-#' @param xSeq vector defaults to seq(-100,100,10). It is recommended to try the default
-#' first. The first argument in the seq function needs to be lower than the minimum value, the second argument 
-#' needs to be higher than the highest value, both should probably be multiples of 10 or 20, 
-#' and the third argument should probably be 5 or 10.  Finally, it is good to have the first and second arguments straddle zero. 
 #' @param flux logical if TRUE, plots flux results, if FALSE plots concentration
+#' @param xMin minimum bin value, it is good to have the xMin and xMax arguments straddle zero. 
+#' @param xMax maximum bin value
+#' @param xStep step size, should probably be multiples of 10 or 20
 #' @param printTitle logical if TRUE, includes title
 #' @param cex.main numeric title font size
+#' @param cex.axis numeric axis font size
+#' @param cex.lab numeric label font size
 #' @param col.fill character fill color
 #' @param \dots base R graphical parameters that can be passed to the hist function
 #' @export
@@ -348,8 +349,9 @@ ciBands <- function(eList, repAnnualResults, probs=c(0.05,0.95)){
 #' plotHistogramTrend(eList, eBoot, caseSetUp, 
 #'                    flux=TRUE, xSeq = seq(-20,60,5))
 #' }
-plotHistogramTrend <- function (eList, eBoot, caseSetUp, xSeq=seq(-100,100,10), 
-                                flux=TRUE, printTitle=TRUE, cex.main=1.1, col.fill="grey",...){
+plotHistogramTrend <- function (eList, eBoot, caseSetUp, 
+                                flux = TRUE, xMin = NA, xMax = NA, xStep = NA,
+                                printTitle=TRUE, cex.main=1.1, cex.axis = 1.1, cex.lab = 1.1, col.fill="grey",...){
 
   periodName <- setSeasonLabel(data.frame(PeriodStart = eList$INFO$paStart, 
                                           PeriodLong = eList$INFO$paLong))
@@ -365,17 +367,25 @@ plotHistogramTrend <- function (eList, eBoot, caseSetUp, xSeq=seq(-100,100,10),
     titleWord <- "Concentration"
   }
 
-  titleToPrint <- ifelse(printTitle, 
-                         paste("Histogram of trend in", 
-                               eList$INFO$paramShortName, "\nFlow Normalized", titleWord, 
-                               caseSetUp$year1, "to", caseSetUp$year2, "\n", eList$INFO$shortName, 
-                               periodName), "")
-  hist(reps, breaks = xSeq, yaxs = "i", xaxs = "i", tcl = 0.5, 
+  titleToPrint <- ifelse(printTitle, paste("Trend magnitude in", 
+                                           eList$INFO$paramShortName, "\nFlow Normalized", titleWord, 
+                                           caseSetUp$year1, "to", caseSetUp$year2, "\n", eList$INFO$shortName, 
+                                           periodName), "")
+  minReps <- min(reps,na.rm = TRUE)
+  maxReps <- max(reps,na.rm = TRUE)
+  xMin <- if(is.na(xMin)) min(-10,minReps) else xMin
+  xMax <- if(is.na(xMax)) max(10,maxReps) else xMax
+  xStep <- if(is.na(xStep)) (xMax-xMin) / 10 else xStep
+  xSeq <- seq(xMin,xMax,xStep)
+  hist(reps, breaks = xSeq, yaxs = "i", xaxs = "i", axes = FALSE, ylab = "",
        main = titleToPrint, freq = FALSE, xlab = xlabel, col = col.fill, 
-       cex.main = cex.main, ...)
+       cex.main = cex.main, cex.lab = cex.lab, ...)
   abline(v = change, lwd = 3, lty = 2)
   abline(v = 0, lwd = 3)
   box()
+  axis(1, tcl = 0.5, labels = TRUE, cex.axis = cex.axis)
+  axis(2, tcl = 0.5, labels = TRUE, las = 1, cex.axis = cex.axis)
+  title(ylab = "Density", line = 4.5, cex.lab = cex.lab)
   axis(3, tcl = 0.5, labels = FALSE)
   axis(4, tcl = 0.5, labels = FALSE)
 }

inst/doc/EGRETci.R

@@ -47,21 +47,37 @@ eList <- Choptank_eList # Example data from EGRET package
 eBoot <- Choptank_eBoot
 caseSetUp <- Choptank_caseSetUp
 
-#Concentration:
+#Concentration an initial run:
 plotHistogramTrend(eList, eBoot, caseSetUp,  
-                   flux=FALSE, xSeq = seq(-20,60,5))
+                   flux=FALSE)
 
-#Flux
+#Flux an initial run:
 plotHistogramTrend(eList, eBoot, caseSetUp,
-                   flux=TRUE, xSeq = seq(-20,60,5))
+                   flux=TRUE)
 
 
 ## ---- histExampleCombo, fig.width=7, fig.height=4---------
 par(mfrow=c(1,2))
-plotHistogramTrend(eList, eBoot, caseSetUp, flux=FALSE,
-                   printTitle=FALSE, ylim=c(0,0.07), xSeq=seq(-10,70,10))
-plotHistogramTrend(eList, eBoot, caseSetUp, flux=TRUE,
-                   printTitle=FALSE, ylim=c(0,0.07), xSeq=seq(-10,70,10))
+#Concentration, presentation version:
+plotHistogramTrend(eList, eBoot, caseSetUp,  
+                   flux=FALSE, xMin = -5, xMax = 65, xStep = 5)
+
+#Flux, presentation version:
+plotHistogramTrend(eList, eBoot, caseSetUp,
+                   flux=TRUE, xMin = -5, xMax = 55, xStep = 5)
+# To return to figures printing in 1 row, 1 columns:
+par(mfrow=c(1,1))                   
+
+## ---- histExampleCombo2, fig.width=7, fig.height=10-------
+par(mfrow=c(2,1))
+#Concentration, presentation version:
+plotHistogramTrend(eList, eBoot, caseSetUp,  
+                   flux=FALSE, xMin = -5, xMax = 65, xStep = 5)
+
+#Flux, presentation version:
+plotHistogramTrend(eList, eBoot, caseSetUp,
+                  flux=TRUE, xMin = -5, xMax = 65, xStep = 5,
+                  printTitle = FALSE)
 # To return to figures printing in 1 row, 1 columns:
 par(mfrow=c(1,1))                   
 

inst/doc/EGRETci.Rmd

@@ -186,11 +186,15 @@ This will result in the creation of the text file called `outputText.txt`, and a
 
 ## Histograms
 
-The function `plotHistogramTrend` plots a histogram of all of the trend magnitudes (expressed as percentage change over the period) from the full set of replicates created by running `wBT`.  The function will plot either Flow-Normalized Concentration trends or Flow-Normalized Flux trends.  The trend magnitudes for all the replicates are stored in the eBoot list.  The first example workflow can be used to plot the two histograms: one for trends in Flow Normalized Concentration and one for Flow Normalized Flux.  
+The function `plotHistogramTrend` plots a histogram of all of the trend magnitudes (expressed in percentage change over the selected period) from the full set of replicates created by `wBT`.  These magnitudes are stored in the eBoot list.  These histograms serve the purpose of providing a graphical impression of the central tendency and the uncertainty about the size of a trend.  They can be used to help answer questions like: "How sure are we that the trend is positive?" or "How sure are we that the decrease is at least a 20% decrease over the period of interest?" The `plotHistogramTrend` function can produce histograms either for Flow Normalized Concentration or for Flow Normalized Flux.    
 
-A comment about the final argument in these functions: The default is `xSeq = seq(-100,100,10)`.  This means that the bins for the histogram will run from -100 to +100 in steps of 10 (the units are percentage change).  This will almost certainly be too broad a set of bins, but the first run, with the xSeq left at its devault values (meaning that you can just leave it out of the argument list) will give you an idea of the true range that you need to consider.  Then, run it again with a much tighter range that still includes all of the replicate trends.  For example, suppose we see that the replicates seen in the first plot run tell us this: the lowest value is between -30 and -20 and the highest value is between 40 and 50, then when we give the command again we should say `xSeq = seq(-30,50,10)`.  If it turns out that all the values are much closer to zero we might use something like `xSeq = seq(-15,5,5)` so the bins are a width of 5 rather than 10.  It is also a good idea to have the full width of the histogram include zero.  For example, in the case where the values were all in the range of +20 to +70, you might use something like `xSeq = seq(0,70,10)`.  Note that if the desired output histogram needs to be smoother, then the original `wBT` run needs to set `bootBreak` to a larger value such as 100 or 200.  The larger `bootBreak` is, the smoother the histogram will be.    
+When creating a histogram it is important to carefully determine what "bins" to use.  In this plot it is a good idea to make sure that the plot covers a range of both negative and positive values even though it might be the case that the entire range of the trend magnitudes is positive (or negative).  The function is set up so that it defaults to setting up 10 bins and the bins are set to cover the full range of the magnitudes and at least cover magnitudes as low as -10% and as high as +10%.  The bins are set up based on a minimum value (xMin), a maximum value (xMax), and a step size (xStep).  The defaults used in an initial run can simply leave out these three arguments and will produce a plot.  This first plot would not be suitable for publication or display, but is intended to guide the analyst to set these three arguments to produce a suitable plot.  
 
-Here is an example workflow where the analyst has already determined what a good xSeq argument might be.
+Ulimately we want a plot with these characteristics:  The value of xStep should be an integer such as 5, 10, or 20 so the width of the bins is 5%, 10% or 20%.  xStep should be small enough so that the histogram looks relatively smooth.  If it is too small the histogram bars will tend to be irregular.  xStep should not be so large that there is poor resolution of the shape of the distribution.  If a very smooth histogram is desired then the analyst needs to have run the `wBT` function with a high number of replicates (say `nBoot = 200`). The value of xStep should be less than the minimum value observed in the set of bootstrap replicates.  xMin should be some negative integer multiple of step size.  If no values were negative, then setting it at -1 times xStep would be a good choice.  xMax should be larger than the maximum value in the set of bootstrap replicates and should be some positive integer multiple of xStep.  If there are no positive values then 1 times xStep would be a good choice.
+
+Note that the plots always show a vertical solid line at 0% trend, simply as a point of reference for the "no trend" result.  They also show a dashed vertical line at the value of the ordinary WRTDS estimate of the trend magnitude in percent.  This vertical line fall near the median value for all the bootstrap replicates.
+
+This initial run will give the analyst the basic information needed to set up a good set of values for xMin, xMax and xStep.
 
 ```{r, fig.height=6, fig.width=6}
 library(EGRET)
@@ -201,30 +205,51 @@ eList <- Choptank_eList # Example data from EGRET package
 eBoot <- Choptank_eBoot
 caseSetUp <- Choptank_caseSetUp
 
-#Concentration:
+#Concentration an initial run:
 plotHistogramTrend(eList, eBoot, caseSetUp,  
-                   flux=FALSE, xSeq = seq(-20,60,5))
+                   flux=FALSE)
 
-#Flux
+#Flux an initial run:
 plotHistogramTrend(eList, eBoot, caseSetUp,
-                   flux=TRUE, xSeq = seq(-20,60,5))
+                   flux=TRUE)
 
 ```
 
-Alternatively, the two plots can be shown side-by-side using a workflow like this.
+Having seen these initial versions of the plot one can run the functions again with a set of arguments that will produce a plot that is suitable for presentation or publication. 
+
+Alternatively, the two plots can be shown side-by-side using a workflow. Use the `par` function to set up both functions to plot side-by-side:
 
-Use the `par` function to set up both functions to plot side-by-side:
 
 ```{r , histExampleCombo, fig.width=7, fig.height=4}
 par(mfrow=c(1,2))
-plotHistogramTrend(eList, eBoot, caseSetUp, flux=FALSE,
-                   printTitle=FALSE, ylim=c(0,0.07), xSeq=seq(-10,70,10))
-plotHistogramTrend(eList, eBoot, caseSetUp, flux=TRUE,
-                   printTitle=FALSE, ylim=c(0,0.07), xSeq=seq(-10,70,10))
+#Concentration, presentation version:
+plotHistogramTrend(eList, eBoot, caseSetUp,  
+                   flux=FALSE, xMin = -5, xMax = 65, xStep = 5)
+
+#Flux, presentation version:
+plotHistogramTrend(eList, eBoot, caseSetUp,
+                   flux=TRUE, xMin = -5, xMax = 55, xStep = 5)
 # To return to figures printing in 1 row, 1 columns:
 par(mfrow=c(1,1))                   
 ```
 
+Or:
+```{r , histExampleCombo2, fig.width=7, fig.height=10}
+par(mfrow=c(2,1))
+#Concentration, presentation version:
+plotHistogramTrend(eList, eBoot, caseSetUp,  
+                   flux=FALSE, xMin = -5, xMax = 65, xStep = 5)
+
+#Flux, presentation version:
+plotHistogramTrend(eList, eBoot, caseSetUp,
+                  flux=TRUE, xMin = -5, xMax = 65, xStep = 5,
+                  printTitle = FALSE)
+# To return to figures printing in 1 row, 1 columns:
+par(mfrow=c(1,1))                   
+```
+
+
+
 
 # Confidence Bands