doseresNMA:

# goodness of fit plot with some statistics (I took it from BUGSnet package)
require(R2jags)
require(dplyr)

fitplot<- function(m){
  jagssamples <- as.mcmc(m) # I NEED R2jags and coda packages to run that
  samples = do.call(rbind, jagssamples)%>% data.frame()
  dev <- samples %>% select(., starts_with("resdev"))
  totresdev <- samples$totresdev %>% mean()
  pmdev <- colMeans(dev)
  rhat <- samples %>% select(., starts_with("rhat"))
  r <- samples %>% select(., starts_with("r."))
  n <- samples %>% select(., starts_with("n"))
  rtilde <- rhat %>%
    colMeans() %>%
    matrix(nrow=1, ncol=ncol(rhat)) %>%
    data.frame() %>%
    slice(rep(1:n(), each = nrow(rhat)))
  
  pmdev_fitted <- 2*(r*log(r/rtilde)+(n-r)*log((n-r)/(n-rtilde)))[1,]
  leverage = pmdev-pmdev_fitted
  DIC = sum(leverage,na.rm = T) + totresdev
  sign = sign(colMeans(r)-colMeans(rhat))
  w = sign*sqrt(as.numeric(pmdev))
  pD = sum(leverage,na.rm = T)
  eq = function(x,c){c-x^2}
  x=seq(-3, 3, 0.001)
  c1=eq(x, c=rep(1, 6001))
  
  plot(w, leverage, xlab=expression('w'[ik]), ylab=expression('leverage'[ik]),
       ylim=c(0, max(c1+3, na.rm=TRUE)*1.15), xlim=c(-3,3),las=1)
  points(x, ifelse(c1 < 0, NA, c1),   lty=1, col="firebrick3",    type="l")
  points(x, ifelse(c1 < -1, NA, c1)+1, lty=2, col="chartreuse4",   type="l")
  points(x, ifelse(c1 < -2, NA, c1)+2, lty=3, col="mediumpurple3", type="l")
  points(x, ifelse(c1 < -3, NA, c1)+3, lty=4, col="deepskyblue3",  type="l")
  text(2, 4.3, paste("pD=", round(pD, 2)), cex = 0.8)
  text(2, 3.9, paste("Dres=", round(totresdev,2)), cex = 0.8)
  text(2, 3.5, paste("DIC=", round(DIC,2)), cex = 0.8)
}