metab_night calc of metabolism giving NA's

[AbigailUVM]

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.

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