Sorry about that, please find here:
suppressPackageStartupMessages({
library(tidyverse)
library(streamMetabolizer)
library(ggpubr)
library(lubridate)
library(rstan)
library(StanHeaders)
library(chron)
library(Rcpp)
library(RcppArmadillo)
library(reprex)
library(datapasta)
})
search()
packageVersion("tidyverse")
# SET DIRECTORY AND ADJUST DATE/TIME; ADDED UPDATED SOLAR TIME EDITS
setwd("C:\\Users\\Administrator\\Desktop")
getwd()
kup_in <- read.csv("reprex_data.csv") #download Input file
names(kup_in)[1:2] <- c("Date","Time") #Had to manually set colnames bc Windows messed up formatting
#Converts excel date/times to R date/times
options(chron.origin = c(month=1, day=1, year=1900))
TC <- chron((kup_in$UStc-2), format = c(dates="Day/Month/Year", times = "h:m:s"))
#Chron datetimes, change format
chron.time <- chron::chron(TC)
time.format <- "%Y-%m-%d %H:%M:%S"
text.time <- format(chron.time, time.format) #as.POSIXct works poorly with chron
posix.time.localtz <- as.POSIXct(text.time, format=time.format, tz='America/Anchorage')
#Convert to POSIXct time and force timzone to America/Anchorage
solar.time <- streamMetabolizer::calc_solar_time(posix.time.localtz, longitude=-149.39)
#CREATE METABOLIZER DATA COMPONENTS
kup_in$SolarTime <- streamMetabolizer::calc_solar_time(posix.time.localtz, longitude=-149.39)
new.solar.time <- kup_in$SolarTime
BP <- kup_in$BP
DO.obs <- kup_in$USDO
light <- kup_in$Light
Q <- kup_in$Q
depth <- kup_in$Depth
tempC <- kup_in$USt
kup_in$tempK <- paste((kup_in$USt)+273.15) #create col for temp in Kelvin
tempK <- kup_in$tempK #create vector for kelvin temp
tempK <- as.numeric(tempK) #change class to numeric
kup_in$DO.sat <- exp(-139.34411+
((1.575701*10^5)/tempK)-
((6.642308*10^7)/tempK^2)+
((1.243800*10^10)/tempK^3)-
((8.621949*10^11)/tempK^4)) #create col in data set for DOsat
DO.sat <- exp(-139.34411+
((1.575701*10^5)/tempK)-
((6.642308*10^7)/tempK^2)+
((1.243800*10^10)/tempK^3)-
((8.621949*10^11)/tempK^4)) #Use Benson&Krause eq. to calc satDO
est.k600 <- kup_in$k600
#GENERATE BAYESIAN INPUT DATA FILE
dat_bayes <- NULL
dat_bayes <- data.frame(new.solar.time,DO.obs,DO.sat,depth,tempC,light,Q,est.k600)
colnames(dat_bayes)<- c("solar.time","DO.obs","DO.sat","depth","temp.water","light","discharge","k600")
str(dat_bayes)
# Set the model type.
bayes_name <- mm_name(type='bayes', pool_K600='none',
err_obs_iid=TRUE, err_proc_acor=FALSE, err_proc_iid=TRUE,
ode_method = 'trapezoid', deficit_src='DO_mod', engine='stan')
#### Set the model specs. For final run pump up the burnin steps to several thousand
bayes_specs_2 <- revise(specs(bayes_name),
burnin_steps=1000, saved_steps=500,
day_start=3, day_end=27,
GPP_daily_mu = 2, GPP_daily_lower = 0,GPP_daily_sigma = 2,
ER_daily_mu = -4, ER_daily_upper = 0, ER_daily_sigma = 3,
K600_daily_meanlog = log(12),
K600_daily_sdlog = 0.01,
chains=1)
bayes_specs_2
#remove k600 and discharge from dat_bayes
dat_bayes <- dat_bayes[,c(-7:-8)]
#run the model with the specs
bayes_fit_kup_2 <- metab(bayes_specs_2, data=dat_bayes)
data.frame(
stringsAsFactors = FALSE,
Date = c("6/28/2018","6/28/2018",
"6/28/2018","6/28/2018","6/28/2018"),
Time = c("12:00:00 AM","12:05:00 AM",
"12:10:00 AM","12:15:00 AM","12:20:00 AM"),
UStc = c(43279,43279.00347,43279.00694,
43279.01042,43279.01389),
USt = c(9.13, 9.1, 9.07, 9.03, 9),
Light = c(55.37360911,55.93355151,
53.89058287,56.48904704,58.9929764),
Depth = c(0.37036501,0.370455948,
0.370455948,0.370000078,0.370000078)
suppressPackageStartupMessages({
library(tidyverse)
library(streamMetabolizer)
library(ggpubr)
library(lubridate)
library(rstan)
library(StanHeaders)
library(chron)
library(Rcpp)
library(RcppArmadillo)
library(reprex)
})
search()
#> [1] ".GlobalEnv" "package:reprex"
#> [3] "package:RcppArmadillo" "package:Rcpp"
#> [5] "package:chron" "package:rstan"
#> [7] "package:StanHeaders" "package:lubridate"
#> [9] "package:ggpubr" "package:streamMetabolizer"
#> [11] "package:forcats" "package:stringr"
#> [13] "package:dplyr" "package:purrr"
#> [15] "package:readr" "package:tidyr"
#> [17] "package:tibble" "package:ggplot2"
#> [19] "package:tidyverse" "package:stats"
#> [21] "package:graphics" "package:grDevices"
#> [23] "package:utils" "package:datasets"
#> [25] "package:methods" "Autoloads"
#> [27] "tools:callr" "package:base"
packageVersion("tidyverse")
#> [1] '1.3.0'
# SET DIRECTORY AND ADJUST DATE/TIME; ADDED UPDATED SOLAR TIME EDITS
setwd("C:\\Users\\Administrator\\Desktop")
getwd()
#> [1] "C:/Users/Administrator/Desktop"
kup_in <- read.csv("reprex_data.csv") #download Input file
names(kup_in)[1:2] <- c("Date","Time") #Had to manually set colnames bc Windows messed up formatting
#Converts excel date/times to R date/times
options(chron.origin = c(month=1, day=1, year=1900))
TC <- chron((kup_in$UStc-2), format = c(dates="Day/Month/Year", times = "h:m:s"))
#Chron datetimes, change format
chron.time <- chron::chron(TC)
time.format <- "%Y-%m-%d %H:%M:%S"
text.time <- format(chron.time, time.format) #as.POSIXct works poorly with chron
posix.time.localtz <- as.POSIXct(text.time, format=time.format, tz='America/Anchorage')
#Convert to POSIXct time and force timzone to America/Anchorage
solar.time <- streamMetabolizer::calc_solar_time(posix.time.localtz, longitude=-149.39)
#CREATE METABOLIZER DATA COMPONENTS
kup_in$SolarTime <- streamMetabolizer::calc_solar_time(posix.time.localtz, longitude=-149.39)
new.solar.time <- kup_in$SolarTime
BP <- kup_in$BP
DO.obs <- kup_in$USDO
light <- kup_in$Light
Q <- kup_in$Q
depth <- kup_in$Depth
tempC <- kup_in$USt
kup_in$tempK <- paste((kup_in$USt)+273.15) #create col for temp in Kelvin
tempK <- kup_in$tempK #create vector for kelvin temp
tempK <- as.numeric(tempK) #change class to numeric
kup_in$DO.sat <- exp(-139.34411+
((1.575701*10^5)/tempK)-
((6.642308*10^7)/tempK^2)+
((1.243800*10^10)/tempK^3)-
((8.621949*10^11)/tempK^4)) #create col in data set for DOsat
DO.sat <- exp(-139.34411+
((1.575701*10^5)/tempK)-
((6.642308*10^7)/tempK^2)+
((1.243800*10^10)/tempK^3)-
((8.621949*10^11)/tempK^4)) #Use Benson&Krause eq. to calc satDO
est.k600 <- kup_in$k600
#GENERATE BAYESIAN INPUT DATA FILE
dat_bayes <- NULL
dat_bayes <- data.frame(new.solar.time,DO.obs,DO.sat,depth,tempC,light,Q,est.k600)
colnames(dat_bayes)<- c("solar.time","DO.obs","DO.sat","depth","temp.water","light","discharge","k600")
str(dat_bayes)
#> 'data.frame': 287 obs. of 8 variables:
#> $ solar.time: POSIXct, format: "2018-06-27 22:04:04" "2018-06-27 22:09:03" ...
#> $ DO.obs : num 10.6 10.6 10.6 10.6 10.6 ...
#> $ DO.sat : num 11.5 11.5 11.5 11.6 11.6 ...
#> $ depth : num 0.37 0.37 0.37 0.37 0.37 ...
#> $ temp.water: num 9.13 9.1 9.07 9.03 9 8.97 8.94 8.91 8.87 8.84 ...
#> $ light : num 55.4 55.9 53.9 56.5 59 ...
#> $ discharge : num 4.86 4.87 4.87 4.84 4.84 ...
#> $ k600 : num 0.726 0.726 0.726 0.725 0.727 ...
# Set the model type.
bayes_name <- mm_name(type='bayes', pool_K600='none',
err_obs_iid=TRUE, err_proc_acor=FALSE, err_proc_iid=TRUE,
ode_method = 'trapezoid', deficit_src='DO_mod', engine='stan')
#### Set the model specs. For final run pump up the burnin steps to several thousand
bayes_specs_2 <- revise(specs(bayes_name),
burnin_steps=1000, saved_steps=500,
day_start=3, day_end=27,
GPP_daily_mu = 2, GPP_daily_lower = 0,GPP_daily_sigma = 2,
ER_daily_mu = -4, ER_daily_upper = 0, ER_daily_sigma = 3,
K600_daily_meanlog = log(12),
K600_daily_sdlog = 0.01,
chains=1)
bayes_specs_2
#> Model specifications:
#> model_name b_np_oipi_tr_plrckm.stan
#> engine stan
#> split_dates FALSE
#> keep_mcmcs TRUE
#> keep_mcmc_data TRUE
#> day_start 3
#> day_end 27
#> day_tests full_day, even_timesteps, complete_data, pos_disch...
#> required_timestep NA
#> GPP_daily_mu 2
#> GPP_daily_lower 0
#> GPP_daily_sigma 2
#> ER_daily_mu -4
#> ER_daily_upper 0
#> ER_daily_sigma 3
#> K600_daily_meanlog 2.484906649788
#> K600_daily_sdlog 0.01
#> err_obs_iid_sigma_scale 0.03
#> err_proc_iid_sigma_scale 5
#> params_in GPP_daily_mu, GPP_daily_lower, GPP_daily_sigma, ER...
#> params_out GPP, ER, DO_R2, GPP_daily, ER_daily, K600_daily, e...
#> n_chains 4
#> n_cores 4
#> burnin_steps 1000
#> saved_steps 500
#> thin_steps 1
#> verbose FALSE
#> chains 1
#remove k600 and discharge from dat_bayes
dat_bayes <- dat_bayes[,c(-7:-8)]
#run the model with the specs
bayes_fit_kup_2 <- metab(bayes_specs_2, data=dat_bayes)
#> Warning: `select_()` is deprecated as of dplyr 0.7.0.
#> Please use `select()` instead.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_warnings()` to see where this warning was generated.
#> Warning: `select_vars()` is deprecated as of dplyr 0.8.4.
#> Please use `tidyselect::vars_select()` instead.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_warnings()` to see where this warning was generated.
#> Warning: `mutate_()` is deprecated as of dplyr 0.7.0.
#> Please use `mutate()` instead.
#> See vignette('programming') for more help
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_warnings()` to see where this warning was generated.
#> Warning: `data_frame()` is deprecated as of tibble 1.1.0.
#> Please use `tibble()` instead.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_warnings()` to see where this warning was generated.
#> Error: All columns in a tibble must be vectors.
#> x Column `solar.time` is NULL.
#> Timing stopped at: 1.44 0 1.43