## ----global_options, include=FALSE--------------------------------------------
knitr::opts_chunk$set(fig.width=8.5, fig.height=5.5, fig.align='center', warning=FALSE, message=FALSE)

## ----echo = F, message = F----------------------------------------------------
library(DHARMa)
set.seed(123)

## ----eval = F-----------------------------------------------------------------
# library(rjags)
# library(BayesianTools)
# 
# set.seed(123)
# 
# dat <- DHARMa::createData(200, overdispersion = 0.2)
# 
# Data = as.list(dat)
# Data$nobs = nrow(dat)
# Data$nGroups = length(unique(dat$group))
# 
# modelCode = "model{
# 
#   for(i in 1:nobs){
#     observedResponse[i] ~ dpois(lambda[i])  # poisson error distribution
#     lambda[i] <- exp(eta[i]) # inverse link function
#     eta[i] <- intercept + env*Environment1[i]  # linear predictor
#   }
# 
#   intercept ~ dnorm(0,0.0001)
#   env ~ dnorm(0,0.0001)
# 
#   # Posterior predictive simulations
#   for (i in 1:nobs) {
#     observedResponseSim[i]~dpois(lambda[i])
#   }
# 
# }"
# 
# jagsModel <- jags.model(file= textConnection(modelCode), data=Data, n.chains = 3)
# para.names <- c("intercept","env", "lambda", "observedResponseSim")
# Samples <- coda.samples(jagsModel, variable.names = para.names, n.iter = 5000)
# 
# x = BayesianTools::getSample(Samples)
# 
# colnames(x) # problem: all the variables are in one array - this is better in STAN, where this is a list - have to extract the right columns by hand
# posteriorPredDistr = x[,3:202] # this is the uncertainty of the mean prediction (lambda)
# posteriorPredSim = x[,203:402] # these are the simulations
# 
# sim = createDHARMa(simulatedResponse = t(posteriorPredSim), observedResponse = dat$observedResponse, fittedPredictedResponse = apply(posteriorPredDistr, 2, median), integerResponse = T)
# plot(sim)

## ----eval=F-------------------------------------------------------------------
#   # Posterior predictive simulations
#   for (i in 1:nobs) {
#     observedResponseSim[i]~dpois(lambda[i])
#   }

## ----eval = F-----------------------------------------------------------------
#   for(i in 1:nobs){
#     observedResponse[i] ~ dpois(lambda[i])  # poisson error distribution
#     lambda[i] <- exp(eta[i]) # inverse link function
#     eta[i] <- intercept + env*Environment1[i] + RE[group[i]] # linear predictor
#   }
# 
#   for(j in 1:nGroups){
#    RE[j] ~ dnorm(0,tauRE)
#   }

## ----eval=F-------------------------------------------------------------------
#   for(j in 1:nGroups){
#    RESim[j] ~ dnorm(0,tauRE)
#   }
# 
#   for (i in 1:nobs) {
#     observedResponseSim[i] ~ dpois(lambdaSim[i])
#     lambdaSim[i] <- exp(etaSim[i])
#     etaSim[i] <- intercept + env*Environment1[i] + RESim[group[i]]
#   }