DOI-USGS / DOI-USGS/streamMetabolizer
Setting priors for K
- Dominant language
- Stan
- Stars
- 47
- Forks
- 25
- PR merge metrics
- No merged PRs in 30d
Description
I have previously measured reaeration rates that I would like to use as priors, but my problem is I don't know how to set that in the code. I have a k rate of 66/day. In issue #367, Bob Hall said to use the following specs:
K600_daily_meanlog_meanlog = 0
K600_daily_meanlog_sdlog =1
I don't have a background in Bayesian models, so I don't quite understand what these specs mean. In issue #367 , I believe the person had very low k600, so I interpret the specs given by Bob would limit the modeled k600 between 0 and 1. Is that what's going on?
Furthermore, what is the best way to convert k to k600? I assume that my measured reaeration rate isn't the standardized gas exchange rate and must be converted somehow.
tldr:
I have measured reaeration rates that I would like to use as priors for K600 but don't know how to specify that properly in the code.
Thanks!
```
library(streamMetabolizer)
setwd("C:/Users/slb0035/Google Drive/School/Working/METABOLISM/streamMetabolizer/INPUT/")
data<-read.csv("sb4_spring2018.csv", stringsAsFactors = FALSE)
data$solar.time<- as.POSIXct(data$solar.time, tz='America/New_York')
lubridate::tz(data$solar.time)
data$solar.time <- streamMetabolizer::calc_solar_time(data$solar.time, longitude=-85) #FBMI longitude
calc_light(data$solar.time, 32, -85)
library(unitted)
data$light <- calc_light(u(data$solar.time), u(32, 'degN'), u(-85, 'degE'), u(3004, 'umol m^-2 s^-1'), attach.units=TRUE)
bayes_name <- mm_name(type='bayes', pool_K600="normal", err_obs_iid=TRUE, err_proc_iid=TRUE)
bayes_specs <- specs(bayes_name, K600_daily_meanlog_meanlog=0, K600_daily_meanlog_sdlog=66)
mm <- metab(bayes_specs, data=data)
output<-as.data.frame(get_params(mm))
write.csv(output, 'sb4_spring2018-output_0-6_kmany.csv') #whatever filename is
pdf('plot_DO_preds_sb4_spring2018output_0-2_kmany.pdf')
plot_DO_preds(mm)
dev.off()
pdf('plot_metab_preds_sb4_spring2018output_0-2_kmany.pdf')
plot_metab_preds(mm)
dev.off()
mcmc <- get_mcmc(mm)
pdf('mcmc_sb4_spring2018output_0-2_kmany.pdf')
rstan::traceplot(mcmc, pars='K600_daily', nrow=3)
dev.off()
```
```
Session info ---------------------------------------------------------------------------------------------
setting value
version R version 3.5.1 (2018-07-02)
system x86_64, mingw32
ui RStudio (1.1.456)
language (EN)
collate English_United States.1252
tz America/Chicago
date 2018-10-01
Packages -------------------------------------------------------------------------------------------------
package * version date source
assertthat 0.2.0 2017-04-11 CRAN (R 3.5.1)
base * 3.5.1 2018-07-02 local
bindr 0.1.1 2018-03-13 CRAN (R 3.5.1)
bindrcpp * 0.2.2 2018-03-29 CRAN (R 3.5.1)
bitops 1.0-6 2013-08-17 CRAN (R 3.5.0)
codetools 0.2-15 2016-10-05 CRAN (R 3.5.1)
colorspace 1.3-2 2016-12-14 CRAN (R 3.5.1)
compiler 3.5.1 2018-07-02 local
crayon 1.3.4 2017-09-16 CRAN (R 3.5.1)
curl 3.2 2018-03-28 CRAN (R 3.5.1)
datasets * 3.5.1 2018-07-02 local
deSolve 1.21 2018-05-09 CRAN (R 3.5.0)
devtools 1.13.6 2018-06-27 CRAN (R 3.5.1)
digest 0.6.17 2018-09-12 CRAN (R 3.5.1)
dplyr 0.7.6 2018-06-29 CRAN (R 3.5.1)
ggplot2 * 3.0.0 2018-07-03 CRAN (R 3.5.1)
glue 1.3.0 2018-07-17 CRAN (R 3.5.1)
graphics * 3.5.1 2018-07-02 local
grDevices * 3.5.1 2018-07-02 local
grid 3.5.1 2018-07-02 local
gridExtra 2.3 2017-09-09 CRAN (R 3.5.1)
gtable 0.2.0 2016-02-26 CRAN (R 3.5.1)
httr 1.3.1 2017-08-20 CRAN (R 3.5.1)
inline 0.3.15 2018-05-18 CRAN (R 3.5.1)
jsonlite 1.5 2017-06-01 CRAN (R 3.5.1)
labeling 0.3 2014-08-23 CRAN (R 3.5.0)
LakeMetabolizer 1.5.0 2016-06-23 CRAN (R 3.5.1)
lazyeval 0.2.1 2017-10-29 CRAN (R 3.5.1)
lubridate 1.7.4 2018-04-11 CRAN (R 3.5.1)
magrittr 1.5 2014-11-22 CRAN (R 3.5.1)
memoise 1.1.0 2017-04-21 CRAN (R 3.5.1)
methods * 3.5.1 2018-07-02 local
munsell 0.5.0 2018-06-12 CRAN (R 3.5.1)
parallel 3.5.1 2018-07-02 local
pillar 1.3.0 2018-07-14 CRAN (R 3.5.1)
pkgconfig 2.0.2 2018-08-16 CRAN (R 3.5.1)
plyr 1.8.4 2016-06-08 CRAN (R 3.5.1)
purrr 0.2.5 2018-05-29 CRAN (R 3.5.1)
R6 2.2.2 2017-06-17 CRAN (R 3.5.1)
Rcpp 0.12.18 2018-07-23 CRAN (R 3.5.1)
RCurl 1.95-4.11 2018-07-15 CRAN (R 3.5.1)
reshape2 1.4.3 2017-12-11 CRAN (R 3.5.1)
rLakeAnalyzer 1.11.4 2018-03-14 CRAN (R 3.5.1)
rlang 0.2.2 2018-08-16 CRAN (R 3.5.1)
rstan * 2.17.3 2018-01-20 CRAN (R 3.5.1)
rstudioapi 0.7 2017-09-07 CRAN (R 3.5.1)
scales 1.0.0 2018-08-09 CRAN (R 3.5.1)
StanHeaders * 2.17.2 2018-01-20 CRAN (R 3.5.1)
stats * 3.5.1 2018-07-02 local
stats4 3.5.1 2018-07-02 local
streamMetabolizer * 0.10.9 2018-05-09 local
stringi 1.1.7 2018-03-12 CRAN (R 3.5.0)
stringr 1.3.1 2018-05-10 CRAN (R 3.5.1)
tibble 1.4.2 2018-01-22 CRAN (R 3.5.1)
tidyr 0.8.1 2018-05-18 CRAN (R 3.5.1)
tidyselect 0.2.4 2018-02-26 CRAN (R 3.5.1)
tools 3.5.1 2018-07-02 local
unitted * 0.2.9 2018-05-09 local
utils * 3.5.1 2018-07-02 local
withr 2.1.2 2018-03-15 CRAN (R 3.5.1)
XML 3.98-1.16 2018-08-19 CRAN (R 3.5.1)
yaml 2.2.0 2018-07-25 CRAN (R 3.5.1)
```
Contributor guide
Research direction
Start with the mm_name(...), specs(...), and metab(...) calls in the report, then compare the K600_daily_meanlog settings with issue #367. Run the supplied example and inspect get_params(mm) and get_mcmc(mm) to verify the prior behaves as intended; done means the project documents how a measured 66/day reaeration rate should be expressed as a K600 prior.
Written by the indexing model from the issue text.
Assessment
- Tech stack
- r
- Domain
- data
- Issue type
- Documentation
- Difficulty
- 5/5
- Estimated time
- Over a week
- Activity status
- Stale
- Clarity
- Needs clarification
- Newbie friendliness
- 25/100