DOI-USGS / DOI-USGS/streamMetabolizer

metab_night calc of metabolism giving NA's

Open
#400 0 comments 0 reactions 0 assignees View on GitHub
Dominant language
Stan
Stars
47
Forks
25
PR merge metrics
No merged PRs in 30d

Description

Hi there! I ran a night time regression using metab_night in order to use these k600 values as priors for the metab() model by coercing my PAR to zero between the hours of 8pm and 5am (I'm using summer Arctic data, so PAR doesn't reach zero until August). However, I get several NA's for my k600 values after running the metab_night model--I am curious if this is an issue with my data/coding or if there is a way to mediate this?

I am also receiving an error message when running the metab() model that I have been unable to diagnose and ameliorate (note at end of code). Are you familiar with this issue? I would be very appreciative of any guidance you can provide regarding these questions.

```r
kup_in <- read.csv("2.5k_metabolizer_input.csv") #download Input file
kup_in$DateTime <- paste(kup_in$Date,kup_in$Time) #combine date/time cols
kup_in$DateTime <- lubridate::mdy_hms(kup_in$DateTime) #read date/time correctly
local.time <- kup_in[,14] #create vector name to change time format to solar
posix.local.time <- as.POSIXct(local.time, format="%Y-%m-%d %H:%M:%S", tz='America/Anchorage')

#CREATE METABOLIZER DATA COMPONENTS
solar.time <- streamMetabolizer::calc_solar_time(posix.local.time, longitude=-149.39)
kup_in$SolarTime <- streamMetabolizer::calc_solar_time(posix.local.time, longitude=-149.39)
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 sat DO

#GENERATE BAYESIAN INPUT DATA FILE
dat_bayes <- NULL
dat_bayes <- data.frame(solar.time,DO.obs,DO.sat,depth,tempC,light,Q)
colnames(dat_bayes)<-c("solar.time","DO.obs","DO.sat","depth","temp.water","light","discharge")

#new data.frame with the daily mean values of each parameter
daily <- dat_bayes %>%
group_by(date=as.Date(solar.time)) %>%
summarise_all(mean, na.rm=T)

#Prepare data for NTR, I have to force the light to go to 0 from 20-5
dat_night <- dat_bayes %>%
mutate(light=if_else(hour(solar.time) %in% c(0:5,20:24), 0, light) )

#Run night-time regression for all data
mm2 <- metab_night(specs = specs(mm_name("night")), dat_night)

#see output from this model
K600_nights_kup <- get_params(mm2, uncertainty='ci') %>%
select(date, K600.daily, K600.daily.lower, K600.daily.upper) %>%
left_join(daily) %>%
mutate(K_source="night") # 28 dates did not get assigned k600 values

#28 rows are missing from this data
ggplot(K600_nights_kup)+ geom_pointrange( aes(discharge, y=K600.daily, ymin=K600.daily.lower,
ymax=K600.daily.upper), size=1) + theme_bw()

#This is just for the priors so are quite uninformed
Q_K_NTR <- lm(K600.daily~poly(discharge,4), data=K600_nights_kup)
summary(Q_K_NTR)

#Calculate the bins of discharge, I do 4 given the shape of the K data,
#but I do some trial and error here to catch the shape of the data
bins_Q_kup <- data.frame(discharge = calc_bins(daily$discharge, 'number', n=4)$bounds)

#Now I add K values for the nodes of discharge
pred_K_NTR <- data.frame(K.600 = predict(Q_K_NTR, newdata=bins_Q_kup, se.fit = T), discharge=bins_Q_kup$discharge)

#comparison of the K nodes and the NTR data
ggplot()+
geom_point(data=K600_nights_kup, aes(discharge, K600.daily))+
geom_pointrange(data = pred_K_NTR, aes(x=discharge, y=K.600.fit,ymin=K.600.fit-K.600.se.fit, ymax=K.600.fit+K.600.se.fit ), color='red') + theme_bw()

# Set the model type.
bayes_name <- mm_name(type='bayes', pool_K600='binned',
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 <- revise(specs(bayes_name),
burnin_steps=1000, saved_steps=500,
day_start=3, day_end=27,
GPP_daily_mu = 1, GPP_daily_lower = 0,GPP_daily_sigma = 2,
ER_daily_mu = -6, ER_daily_upper = 0, ER_daily_sigma = 3,
K600_lnQ_nodes_centers= log(pred_K_NTR$discharge),
K600_lnQ_nodes_meanlog= log(pred_K_NTR$K.600.fit),
K600_daily_sigma_sigma=0.7, chains=1) #
bayes_specs

#run the model with the specs
bayes_fit_kup <- metab(bayes_specs, data=dat_bayes)

#Warning message below
Chain 3: Iteration: 1 / 1500 [ 0%] (Warmup)
[1] "Error in sampler$call_sampler(args_list[[i]]) : "
[2] " c++ exception (unknown reason)"
error occurred during calling the sampler; sampling not done
here are whatever error messages were returned
`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.Modeling failed
Warnings:
some chains had errors; consider specifying chains = 1 to debug
Stan model 'b_Kb_oipi_tr_plrckm' does not contain samples.
```
[2.5k_TXT_metabolizer_input.txt](https://github.com/USGS-R/streamMetabolizer/files/5495993/2.5k_TXT_metabolizer_input.txt)

Contributor guide

Open the contributing guide

Assessment

This issue has not been assessed yet.

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.