doi-usgs / egretci Goto Github PK

View Code? Open in Web Editor NEW

8.0 9.0 9.0 27.71 MB

A bootstrap method for estimating uncertainty of water quality trends. https://doi-usgs.github.io/EGRETci/

Home Page: http://doi-usgs.github.io/EGRETci/

License: Other

R 89.21% CSS 7.33% HTML 3.47%

r usgs egret bootstrap water-quality-trends rstats r-package

egretci's People

Contributors

Stargazers

Watchers

Forkers

ldecicco-usgs gpiorkowski msimoneau-mddelcc lindsayplatt wdwatkins thodson-usgs waternk jcmartinmu kranzfafka

egretci's Issues

Think about parallelizing wBT and flexFNci

The nBoot runs are independent, I think it would be OK. The screen messages would be either impossible or harder to produce, but that's probably a small price to pay. Adding the complexity is the real question.

Only one dashed line for CIs

There is only one dashed line for the 95% CIs in my plots from plotConcHistBoot(). Some times it looks like an upper 95% CI, other times I have no idea:

error in wBT runs- results must align with PA?

Hello,

Getting this reproducible error while running wBT and runPairs:
Error in if (z) "Reject Ho" else "Do Not Reject Ho" : missing value where TRUE/FALSE needed

Attached are example input text files used in model runs that produce error and my code as R script.

When using either input text file "as is", I get the error and change calculations are significantly off (Flux Change is actually NA):

"Susquehanna River at Danville, PA Nitrate-Nitrite

Water Year

Bootstrap process, for change from Water Year1986 to Water Year2019
data set runs from 1986 to Year 2019
Bootstrap block length in days 200
bootBreak is 3 confStop is 0.7

WRTDS estimated concentration change is 25256 mg/L
WRTDS estimated flux change is NA 10^6 kg/yr"

When I inputted "fake" results that complete the 1986 water year (in my case, I actually added fake sample results just for September or October), I no longer get the error.

Reproduced this same error for other input files and other sites.

Must the input sample result dates exactly coincide with the period of analysis (PA) in order for wBT to run without this error?

Thank you!

01540500_trend_git.zip
01540500dn - Original.txt
01540500dn.txt

Extremely high and low values in some iterations using iBoot

I am currently rerunning a model that looks at a 10 yr moving window with nBoot = 200, and I keep noticing odd iteration values. Here is some output from 1998-2007:

iBoot, xConc and xFlux 127 -8.615097e+137 -2.975343e+136
iBoot, xConc and xFlux 128 -0.05028564 -1.459841
iBoot, xConc and xFlux 129 0.03666643 -0.1706541
iBoot, xConc and xFlux 130 0.05564319 0.03039397
iBoot, xConc and xFlux 131 0.04101899 -0.1198773
iBoot, xConc and xFlux 132 0.04211883 0.1540202
iBoot, xConc and xFlux 133 -0.09984916 -2.680429
iBoot, xConc and xFlux 134 -0.007287396 -0.3895768
iBoot, xConc and xFlux 135 0.002760677 -0.7174396
iBoot, xConc and xFlux 136 -0.001037794 -0.5118047
iBoot, xConc and xFlux 137 0.0370652 -0.2164727
iBoot, xConc and xFlux 138 -0.02083882 -0.5030731
iBoot, xConc and xFlux 139 0.02780439 -0.294982
iBoot, xConc and xFlux 140 0.003192806 -0.4554299
iBoot, xConc and xFlux 141 -0.05443917 -1.415746
iBoot, xConc and xFlux 142 0.03259536 -0.466467
iBoot, xConc and xFlux 143 0.01166099 -0.1782373
iBoot, xConc and xFlux 144 0.03090554 -0.2704868
iBoot, xConc and xFlux 145 -0.01288284 -1.037976
iBoot, xConc and xFlux 146 0.1241648 0.8563349
iBoot, xConc and xFlux 147 -0.03576329 -1.28158
iBoot, xConc and xFlux 148 2.897706e+30 1.57178e+32
iBoot, xConc and xFlux 149 0.08204259 0.3758117
iBoot, xConc and xFlux 150 -0.01351132 -0.4975678
iBoot, xConc and xFlux 151 0.01373099 -0.1204261
iBoot, xConc and xFlux 152 0.0691038 0.1389613
iBoot, xConc and xFlux 153 0.009624788 -0.4152171
iBoot, xConc and xFlux 154 -0.06310256 -1.229834
iBoot, xConc and xFlux 155 -0.03640907 -1.312824
iBoot, xConc and xFlux 156 0.02100833 -0.2485315
iBoot, xConc and xFlux 157 0.01610422 -0.3146633
iBoot, xConc and xFlux 158 -0.01941998 -1.072092
iBoot, xConc and xFlux 159 -0.06849002 -1.519826
iBoot, xConc and xFlux 160 0.06798843 0.2361849
iBoot, xConc and xFlux 161 0.01927036 -0.3341333
iBoot, xConc and xFlux 162 0.005404621 -0.6666091
iBoot, xConc and xFlux 163 0.04532531 -0.4035346
iBoot, xConc and xFlux 164 0.1073366 0.2855145
iBoot, xConc and xFlux 165 -0.02089126 -1.265881
iBoot, xConc and xFlux 166 -0.03983874 -1.183583
iBoot, xConc and xFlux 167 -1.307639e+79 -4.516112e+77
iBoot, xConc and xFlux 168 0.02232687 -0.2135413
iBoot, xConc and xFlux 169 -0.02690708 -1.169837
iBoot, xConc and xFlux 170 0.05313271 0.02626181
iBoot, xConc and xFlux 171 0.06698775 0.1982021
iBoot, xConc and xFlux 172 0.03341398 -0.7828989
iBoot, xConc and xFlux 173 0.07338573 0.1200608
iBoot, xConc and xFlux 174 0.07082601 0.122985
iBoot, xConc and xFlux 175 -7.434473e+101 -2.567598e+100
iBoot, xConc and xFlux 176 0.07073133 0.1410189
iBoot, xConc and xFlux 177 0.01842965 -0.3556271
iBoot, xConc and xFlux 178 -0.09302633 -2.252569
iBoot, xConc and xFlux 179 0.07934037 0.3939216
iBoot, xConc and xFlux 180 0.01684882 -0.2692689
iBoot, xConc and xFlux 181 0.02964044 -0.3813411
iBoot, xConc and xFlux 182 0.004074809 -0.3463523
iBoot, xConc and xFlux 183 -0.05169548 -1.534628
iBoot, xConc and xFlux 184 0.03390354 -0.1392435
iBoot, xConc and xFlux 185 0.06523865 0.07551456
iBoot, xConc and xFlux 186 0.0444889 -0.09759831
iBoot, xConc and xFlux 187 0.006336869 -0.4673323
iBoot, xConc and xFlux 188 -0.02044203 -1.04437
iBoot, xConc and xFlux 189 0.03901774 0.08431938
iBoot, xConc and xFlux 190 -0.04834053 -1.364399
iBoot, xConc and xFlux 191 -0.05268685 -1.252619
iBoot, xConc and xFlux 192 0.07842773 0.221899
iBoot, xConc and xFlux 193 -0.02316459 -1.329317
iBoot, xConc and xFlux 194 -0.09582452 -1.921
iBoot, xConc and xFlux 195 0.1049984 0.6393139
iBoot, xConc and xFlux 196 0.02325422 -0.2404297
iBoot, xConc and xFlux 197 -0.06165737 -1.684607
iBoot, xConc and xFlux 198 0.0514432 -0.08819077
iBoot, xConc and xFlux 199 -0.04280059 -1.591571
iBoot, xConc and xFlux 200 0.004167874 -0.5705411

It took 203 iterations to achieve 200 sucessful runs.

Maumee River at Waterville OH
TKN through 2007
Water Year

Change estimates are for 2007 minus 1998

Should we reject Ho that Flow Normalized Concentration Trend = 0 ? Do Not Reject Ho
best estimate of change in concentration is 0.0194 mg/L
Lower and Upper 90% CIs -1.19e-01 7.84e-02
also 95% CIs -1.46e+101 1.04e-01
and 50% CIs -2.75e-02 3.97e-02
approximate two-sided p-value for Conc 0.82
Likelihood that Flow Normalized Concentration is trending up = 0.59 is trending down = 0.41

Should we reject Ho that Flow Normalized Flux Trend = 0 ? Do Not Reject Ho
best estimate of change in flux is -0.409 10^6 kg/year
Lower and Upper 90% CIs -2.66e+00 2.80e-01
also 95% CIs -5.04e+99 4.97e-01
and 50% CIs -1.13e+00 -1.22e-01
approximate two-sided p-value for Flux 0.36
Likelihood that Flow Normalized Flux is trending up = 0.177 is trending down = 0.823

Do you have any idea why some of the iterations are so off? Thank you!

Add set.seed argument for reproducibility

Add uncertainty to CI bands

dif <- flexBoot$z11[1:500] - flexBoot$z00[1:500]
summary(dif)
sortDif <- sort(dif)
trueCI <- quantile(sortDif, probs=c(0.05,0.95), type = 6)
trueWidth <- trueCI[2] - trueCI[1]
nReps <- 500
bootWidth <- rep(NA,nReps)
for (iRep in 1:nReps){
difB <- sample(dif,500,replace = TRUE)
difB <- sort(difB)
CI <- quantile(difB, probs=c(0.05, 0.95), type = 6)
bootWidth[iRep] <- CI[2] - CI[1]
}
print(sd(bootWidth))
summary(bootWidth)
cat("\n", trueWidth)
pct <- 100 *sd(bootWidth) / trueWidth
cat("\nCV of widths in percent", pct)

Please verify xFlux computation

The doc for runPairsBoot.R states:
"xConc, xFlux are vectors of length iboot, of the change in flow normalized concentration of flux computed by each bootstrap replicate (mg/L)"

and xConc, xFlux are computed on line 161-162 of runPairsBoot.R as:
xConc_here <- (2 * regDeltaConc) - (c22 - c11)
xFlux_here <- (2 * regDeltaFlux) - (f22 - f11)

Where regDeltaConc is the change in flow normalized flux computed by runPairs.R:
(c22 - c11)

My question is, if both xConc and regDeltaConc represent the change in flow normalized concentration, why is one: (2 * regDeltaConc) - (c22 - c11)
and the other: (c22 - c11).

Should they both be c22-c11 (and the same for flux)?

plotHistogramTrend() does not plot all the data

Since plotHistogramTrend() has xaxs and yaxs set to "i" it is common for the function not to include all the data (horizontally as well as vertically). I think it would work better if they were set to the default, "r", even though the x axis results a little higher.

more tests....

water year vs calendar year.

I am using the EGRETci package to evaluate water quality trends for data in the Mississippi, Minnesota, and St. Croix Rivers. We investigate trends on a calendar year basis.  I was just looking for some clarification to make sure I am understanding the bootstrap results correctly.

Because I am investigating trends on a calendar year basis, I made sure for my eList that paStart=1 and paLong=12.  For this site the first water quality sample was 2/8/1985 and the last sample was 12/19/2014. I then wanted to see the trend for calendar year 2005 to calendar year 2014 (last 10 years).  Here are the first few lines of my bootstrap output text file:

Minnesota River at Jordan  Chloride

Calendar Year

  Bootstrap process, for change from Water Year 2005 to Water Year 2014
                   data set runs from WaterYear 1985 to Water Year 2015
  Bootstrap block length in days 200
  bootBreak is 39  confStop is 0.7


I am confused by how the program summarizes data by water year.  The program knows I want to evaluate by calendar year, but the summary still says I am looking at water year 2005 to 2014. It then says I have samples from water year 1985 to water year 2015, because I had samples after 10/1/2014.  So does the output mean I am looking at change from calendar year 2005 to 2014 but says “water year” anyways? Or did it change my request to look at water year 2005 to water year 2014 after all?  I just want to make sure my 2014 October –December samples are being used and from this summary I’m not sure they are.

I took a look at this question and I think I know what is going on.  My question to you is ask you to run look at your original WRTDS results for calendar year for the period of interest to you.  

That means you want to do this:

eList <- setPA(paStart=1, paLong=12)
tableChange(eList, yearPoints=c(2005,2014))

from that output take note of the change values for concentration in mg/L and for flux in 10^6 kg/yr.

Now use EGRETci and do this:

caseSetUp <- trendSetUp(eList)

then it will ask you for first water year, tell it 2005,  then last water year, tell it 2014.  

then it will ask for nBoot, bootBreak and blockLength,  You can give it 5 and 5 and 200.  (It won't really matter here, I just want you to be able to get the result quickly - when you really run it you will want bigger values for the first two of these).

Then. 

eBoot <- wBT(eList,caseSetUp)

when you look at the results, near the top you will have two lines that say something like this:  

WRTDS estimated concentration change is xxx  mg/L
 WRTDS estimated flux change is          yyy  10^6 kg/yr

The xxx should match the concentration change and the yyy the flux change you got in tableChange.  

I tried it with an example of my own and it worked fine.  Assuming the numbers work out for you then the interpretation is this.  This is really a trend test on change from calendar year 2005 to calendar year 2014.

The upshot is that we need to change what our automatic output says, but the computations are doing exactly what you want.  When we started designing the software we made it work just for Water Year.  Then we changed to making it more flexible, but we never went back and changed the message that prints out.

Let us know if this experiment works.

Laura: can you look at a way of changing the code to reflect what is really being tested?  We also needs to see that it works for any arbitrary choice of paStart and paLong so it might not even be a whole year.  The paLong and paStart in eList$INFO are driving the test, but the stuff we print doesn't really reflect that.  Probably need to say what the PA actually is and defines the starting and ending  of the period being evaluated.

Error-plotHistogramTrend for Flux

Error in hist.default(reps, breaks = xSeq, yaxs = "i", xaxs = "i", tcl = 0.5, :
some 'x' not counted; maybe 'breaks' do not span range of 'x

Here is the list:
rep Concentration | value nPos post_p Lower Upper | 1 -0.0309 0 0.25 -0.0309 -0.0309 | 2 -0.0168 0 0.167 -0.0309 -0.0168 | 3 -0.03 0 0.125 -0.0309 -0.0168 | 4 -0.0181 0 0.1 -0.0309 -0.0168 | 5 -0.0169 0 0.0833 -0.0309 -0.0168 | 6 -0.0101 0 0.0714 -0.0309 -0.0101 | 7 -0.0636 0 0.0625 -0.0636 -0.0101 | 8 -0.0528 0 0.0556 -0.0636 -0.0101 | 9 -0.0171 0 0.05 -0.0636 -0.0101 | 10 -0.0307 0 0.0455 -0.0636 -0.0101 | 11 -0.013 0 0.0417 -0.0636 -0.0101 | 12 -0.017 0 0.0385 -0.0636 -0.0101 | 13 -0.0452 0 0.0357 -0.0636 -0.0101 | 14 -0.0345 0 0.0333 -0.0636 -0.0101 | 15 -0.0337 0 0.0312 -0.0636 -0.0101 | 16 -0.0265 0 0.0294 -0.0636 -0.0101 | 17 -0.0541 0 0.0278 -0.0636 -0.0101 | 18 -0.0287 0 0.0263 -0.0636 -0.0101 | 19 -0.0486 0 0.025 -0.0636 -0.0101 | 20 -0.0262 0 0.0238 -0.0632 -0.0103 | 21 -0.0859 0 0.0227 -0.0836 -0.0104 | 22 -0.0391 0 0.0217 -0.0825 -0.0105 | 23 -0.00791 0 0.0208 -0.0814 -0.00835 | 24 -0.042 0 0.02 -0.0803 -0.00846 | 25 -0.047 0 0.0192 -0.0792 -0.00857 26 -0.0191 0 0.0185 -0.0781 -0.00868 | 27 -0.0403 0 0.0179 -0.077 -0.00879 | 28 -0.0671 0 0.0172 -0.0774 -0.0089 | 29 -0.021 0 0.0167 -0.0765 -0.00901 | 30 -0.0259 0 0.0161 -0.0756 -0.00912 31 -0.00347 0 0.0156 -0.0746 -0.00613 | 32 -0.0465 0 0.0152 -0.0737 -0.00635 | 33 -0.0117 0 0.0147 -0.0727 -0.00658 | 34 -0.0362 0 0.0143 -0.0718 -0.0068 | 35 -0.00482 0 0.0139 -0.0709 -0.00455 | 36 -0.0547 0 0.0135 -0.0699 -0.00462 | 37 -0.0518 0 0.0132 -0.069 -0.00469 | 38 -0.0171 0 0.0128 -0.0681 -0.00476 | 39 -0.0405 0 0.0125 -0.0671 -0.00482 | Flux
value nPos post_p Lower Upper
0.283 1 0.75 0.283 0.283
0.6621 2 0.833 0.283 0.6621
0.3014 3 0.875 0.283 0.6621
0.7496 4 0.9 0.283 0.7496
0.5048 5 0.917 0.283 0.7496
0.7315 6 0.929 0.283 0.7496
-0.4792 6 0.812 -0.4792 0.7496
-0.6414 6 0.722 -0.6414 0.7496
0.6196 7 0.75 -0.6414 0.7496
0.5275 8 0.773 -0.6414 0.7496
0.8981 9 0.792 -0.6414 0.8981
0.7489 10 0.808 -0.6414 0.8981
0.4748 11 0.821 -0.6414 0.8981
0.2384 12 0.833 -0.6414 0.8981
0.2557 13 0.844 -0.6414 0.8981
0.3773 14 0.853 -0.6414 0.8981
-0.3951 14 0.806 -0.6414 0.8981
0.6072 15 0.816 -0.6414 0.8981
-0.5431 15 0.775 -0.6414 0.8981
0.4925 16 0.786 -0.6364 0.8907
-0.7203 16 0.75 -0.7124 0.8833
-0.0548 16 0.717 -0.7085 0.8759
0.9121 17 0.729 -0.7045 0.9093
0.878 18 0.74 -0.7006 0.9086
| -0.02709 18 0.712 -0.6966 0.9079
0.5373 19 0.722 -0.6927 0.9072
0.4758 20 0.732 -0.6887 0.9065
0.0856 21 0.741 -0.6848 0.9058
0.2997 22 0.75 -0.6808 0.9051
| 0.007642 23 0.758 -0.6769 0.9044
0.6906 24 0.766 -0.6729 0.9037
-0.6041 24 0.742 -0.669 0.903
0.06099 25 0.75 -0.665 0.9023
0.8749 26 0.757 -0.6611 0.9016
0.4456 27 0.764 -0.6571 0.9009
-0.8404 27 0.743 -0.7383 0.9002
-0.2126 27 0.724 -0.7323 0.8995
0.3157 28 0.731 -0.7263 0.8988
0.3362 29 0.738 -0.7203 0.8981

NOTICE: upcoming default branch name change

The master branch of this repository will soon be renamed from master to main, as part of a coordinated change across the USGS-R and USGS-VIZLAB organizations. This is part of a broader effort across the git community to use more inclusive language. For instance, git, GitHub, and GitLab have all changed or are in the process of changing their default branch name.

We will make this change early in the week of February 6, 2022. The purpose of this issue is to give notification of the change and provide information on how to make it go smoothly.

When this change is made, forks and clones of the repository will need to be updated as well. Here are instructions for making the necessary changes.

On GitHub, change the default branch name from master to main (<your repository> -> Settings -> Branches).
Change the default branch on any forks you have, as well as local clones. (See details below)
If you have collaborators on this repository, let them know that they will need to change their forked/local repos. Point them to this issue to facilitate the process!
Search within your repository for "master" so that you can change references (e.g. URLs pointing to specific commit points of files) to point to "main" instead.
When you are done, feel free to close this issue!

Changing default branches on forks and local clones

First, update your local git configuration so that new repositories created locally will have the correct default branch: git config --global init.defaultBranch main.

Now, for any forks, do the following:

Go to <your fork> -> Settings -> Branches and edit the default branch from master to main.
Update the settings for your local clone of this fork to match this change.

git branch -m master main
git fetch origin
git branch -u origin/main main
git remote set-head origin -a

Add run.parallel option for runPairsBoot, runGroupsBoot?

Consider doing it like we did with WRTDS:
https://github.com/DOI-USGS/EGRET/blob/main/R/runSurvReg.R#L67

BlankTime() and cibands

After a call to blankTime to filter a time span with missing data, it would seem that a call to ciBands (after running the required commands prior to calling ciBands) doesn't work? The error message is shown below, it is one I haven't seen before and I think it may be related to the fact that I've applied blankTime higher up in the script:

...
CIAnnualResults <- ciBands(eList, repAnnual, probs)
#Error in quantile.default(manyAnnualResults[iYear, 1, 1:nBoot], prob = probs,  : 
#  missing values and NaN's not allowed if 'na.rm' is FALSE

Improve documentation of runPairBoot

Specifically, define the stuff coming out of bootOut

flexFNci could probably use same logic as wBT

For calculating the right slice of "surface"s

wbt() function documentation rendered incorrectly

The function wBT() returns a list object (eBoot). The vectors "pFlux" and "pConc" are the reported change in flow-normalized values expressed in percent change. However, the documentation that describes what this functions returns does not properly render the percent symbol.

https://github.com/USGS-R/EGRETci/blob/d1c4845e53a8599faef68015596ce20dc5cd316d/R/wBTCode.R#L226

library(EGRETci)
?wBT

Improve error message

If you ask for wBT for years past the model results, you get a cryptic message.

wBT and ciCalculations take a REALLY long time

It takes about 10 minutes for wBT to run and I haven't successfully gotten ciCalculations to run. Any ideas?