In FPRgroup/FPEMcountry: fpemlocal

knitr::opts_chunk$set(echo = TRUE)
knitr::opts_knit$set(root.dir = rprojroot::find_rstudio_root_file())
# knitr::opts_knit$get("root.dir")  # alternative to the previous line
# the default autosave location will depend on this being setup
options(warn=-1)

#make more dynamic later
jagsout <- readRDS("output/runs/68.rds")$run$jagsout
bugs <- jagsout$BUGSoutput
# ptm <- proc.time()
# #### to be replaced with current run
# proc.time() - ptm
# bugs$summary
# bugs_15000 <- res$mod$BUGSoutput
mcmc.array <- bugs$sims.array
library(coda)
invlogit <- function(x){
  exp(x)/(1+exp(x))
}
# "mod.ct", 
# "unmet.ct", 
# "trad.ct", 
# " mu.jn", 
# "logitratio.yunmet.hat.j",
# "logitRmax.c",
# "logitomega.c",
# "logitpmax.c",
# "logitRomega.c",
# "setlevel.c",
# "RT.c",
# "unmet.intercept.c"
par <- bugs$sims.list$unmet.ct %>% invlogit()
unmet <- c()
for (i in 1:dim(par)[3]) {
unmet[i] <- effectiveSize(par[,,i]) %>% round()
}
par <- bugs$sims.list$mod.ct %>% invlogit()
mod <- c()
for (i in 1:dim(par)[3]) {
mod[i] <- effectiveSize(par[,,i]) %>% round()
}
par <- bugs$sims.list$trad.ct %>% invlogit()
trad <- c()
for (i in 1:dim(par)[3]) {
trad[i] <- effectiveSize(par[,,i]) %>% round
}
temp <- data.frame(unmet = unmet, mod = mod, trad = trad)
write.table(temp, "output/neff_logit.txt")



# checking some parameters
PlotTrace <- function(#Traceplot for one parameter
  ### Trace plot for one parameter and add loess smoother for each chain
  parname, mcmc.array,##<< needs to be 3-dimensional array!
  n.chains= NULL, n.sim= NULL, main = NULL){
  if (is.null(main)) main <- parname
  if (is.null(n.sim)) n.sim <- dim(mcmc.array)[1]
  if (is.null(n.chains)) n.chains <- dim(mcmc.array)[2]
  plot(c(mcmc.array[,1,parname]), type = "l", ylab = parname,  main = main,
       ylim = c(min(mcmc.array[,,parname]),max(mcmc.array[,,parname])))
  for (chain in 1:n.chains){
    lines(c(mcmc.array[,chain,parname]), type = "l", col = chain)
  }
  for (chain in 1:n.chains){
    curve(predict(loess(c(mcmc.array[,chain,parname])~seq(1,n.sim)),x), lty = 2, lwd = 3, add = TRUE, type = "l", col = chain)
  }
}

max(bugs$summary[,"Rhat"])
bugs$summary[bugs$summary[,"Rhat"]>1.1,]
write.table(bugs$summary[bugs$summary[,"Rhat"]>1.1,], "rhats.txt")



temp <- MCMCvis::MCMCsummary(jagsout)
temp$Rhat %>% max
max(jagsout$BUGSoutput$summary[,"Rhat"])

MCMCvis::MCMCtrace(jagsout,
          params = "all",
                 ISB = FALSE,
                 pdf = TRUE,
                 # prior = pr,
                 post_zm = TRUE,
          open_pdf = FALSE,
          Rhat = TRUE,
          filename = "output/mcmctrace.pdf")
          #wd = getwd())



# first_year <- 1975
# last_year <- 2030
# I <- 8 #observations
# pdf("traceplots.pdf")
# # PlotTrace("logitRomega.c", mcmc.array)
# PlotTrace("mu.in[4,1]", mcmc.array)
# for (t in 1:length(first_year:last_year)) PlotTrace(paste0("mod.ct[1,",t,"]"), mcmc.array)
# for (t in 1:length(first_year:last_year)) PlotTrace(paste0("trad.ct[1,",t,"]"), mcmc.array)
# for (t in 1:length(first_year:last_year)) PlotTrace(paste0("unmet.ct[1,",t,"]"), mcmc.array)
# for (i in 1:I) PlotTrace(paste0("logitratio.yunmet.hat.i[",i,"]"), mcmc.array)
# dev.off()