Skip to content

Instantly share code, notes, and snippets.

Show Gist options
  • Save florianhartig/b3569979ea43e1f9b35d4fc1e248a200 to your computer and use it in GitHub Desktop.
Save florianhartig/b3569979ea43e1f9b35d4fc1e248a200 to your computer and use it in GitHub Desktop.
############ CREATE ZERO-INFLATED GLMM DATA #################
# This first part creates a dataset with beetles counts across an altitudinal gradient (several plots each observed several years), with a random intercept on year and zero-inflation.
altitude = rep(seq(0,1,len = 50), each = 20)
dataID = 1:1000
spatialCoordinate = rep(seq(0,30, len = 50), each = 20)
# random effects + zeroinflation
plot = rep(1:50, each = 20)
year = rep(1:20, times = 50)
yearRandom = rnorm(20, 0, 1)
plotRandom = rnorm(50, 0, 1)
overdispersion = rnorm(1000, sd = 0.5)
zeroinflation = rbinom(1000,1,0.6)
beetles <- rpois(1000, exp( 0 + 12*altitude - 12*altitude^2
# + overdispersion + plotRandom[plot]
+ yearRandom[year]) * zeroinflation )
data = data.frame(dataID, beetles, altitude, plot, year, spatialCoordinate)
plot(year, altitude, cex = beetles/50, pch =2, main = "Beetle counts across altitudinal gradient\n triangle is proportional to counts")
############## Analysis with JAGS ############################
library(R2jags)
modelData=as.list(data)
modelData = append(data, list(nobs=1000, nplots = 50, nyears = 20))
head(data)
# 1) Fit GLM only
modelstring="
model {
# Likelihood
for (i in 1:nobs) {
lambda[i] <- exp(intercept + alt * altitude[i] + alt2 * altitude[i] * altitude[i])
beetles[i]~dpois(lambda[i])
}
# Fixed effect priors
intercept ~ dnorm(0,0.0001)
alt ~ dnorm(0,0.0001)
alt2 ~ dnorm(0,0.0001)
# Posterior predictive simulations
for (i in 1:nobs) {
beetlesPred[i]~dpois(lambda[i])
}
Prediction <- sum(beetlesPred)
}
"
model=jags(model.file = textConnection(modelstring), data=modelData, n.iter=10000, parameters.to.save = c("intercept", "alt", "alt2", "beetlesPred", "lambda"), DIC = F)
library(DHARMa)
simulations = model$BUGSoutput$sims.list$beetlesPred
pred = apply(model$BUGSoutput$sims.list$lambda, 2, median)
dim(simulations)
sim = createDHARMa(simulatedResponse = t(simulations), observedResponse = data$beetles, fittedPredictedResponse = pred, integerResponse = T)
plotSimulatedResiduals(sim)
# 2) GLMM with random intercept on year, observation-level RE for overdispersion, and zero-inflation
modelstring="
model {
# Likelihood
for (i in 1:nobs) {
lambda[i] <- exp(intercept + alt * altitude[i] + alt2 * altitude[i] * altitude[i] + Ryear[year[i]] + RID[i] ) * Zero[i] + 0.00000001
beetles[i]~dpois(lambda[i])
}
# Fixed effect priors
intercept ~ dnorm(0,0.0001)
alt ~ dnorm(0,0.0001)
alt2 ~ dnorm(0,0.0001)
# Random effects
for (i in 1:nyears) {
Ryear[i]~dnorm(0,sigmaYear)
}
for (i in 1:nobs) {
RID[i]~dnorm(0,sigmaID)
}
# Variance priors
sigmaYear~dgamma(1,2)
sigmaID~dgamma(0.001,0.001)
# Zeroinflation
for (i in 1:nobs) {
Zero[i]~dbern(zeroMu + altZero * altitude[i])
}
zeroMu ~ dunif(0,1)
altZero ~ dnorm(0,0.0001)
# Posterior predictive simulations
for (i in 1:nobs) {
beetlesPred[i]~dpois(lambda[i])
}
Prediction <- sum(beetlesPred)
}
"
model=jags(model.file = textConnection(modelstring), data=modelData, n.iter=10000, parameters.to.save = c("intercept", "alt", "alt2", "beetlesPred", "Ryear", "sigmaYear", "lambda", "altZero", "zeroMu"), DIC = F)
library(DHARMa)
simulations = model$BUGSoutput$sims.list$beetlesPred
pred = apply(model$BUGSoutput$sims.list$lambda, 2, median)
dim(simulations)
sim = createDHARMa(simulatedResponse = t(simulations), observedResponse = data$beetles, fittedPredictedResponse = pred, integerResponse = T)
plotSimulatedResiduals(sim)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment