Dat <- read.table('https://raw.githubusercontent.com/florianhartig/LearningBayes/master/data/LizardData.txt')
plot(Dat$Veg,Dat$Count)
7 Generalized linear mixed models
Note
In this chapter, we will discuss how to fit generalized linear mixed models (GLMMs) in a Bayesian framework, using a Poisson regression on lizard counts as a running example. You will learn
- How to fit a Bayesian Poisson GLM in JAGS and compare it to the frequentist
glm()fit - How to add posterior predictive simulations and check residuals with DHARMa
- How to extend the model with an overdispersion term
- How to add zero-inflation to a count model
7.1 Poisson regression
Standard frequentist GLM
fit <- glm(Count ~ Veg, data = Dat, family = "poisson")
summary(fit)
Call:
glm(formula = Count ~ Veg, family = "poisson", data = Dat)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.62494 0.01923 136.50 <2e-16 ***
Veg 0.26236 0.01736 15.11 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 3015.9 on 199 degrees of freedom
Residual deviance: 2785.7 on 198 degrees of freedom
AIC: 3433.1
Number of Fisher Scoring iterations: 6library(rjags)
library(DHARMa)
library(BayesianTools)
# Model specification
model = "
model{
# Likelihood
for(i in 1:n.dat){
# poisson model p(y|lambda)
y[i] ~ dpois(lambda[i])
# log link function
log(lambda[i]) <- mu[i]
# linear predictor on the log scale
mu[i] <- alpha + beta.Veg*Veg[i] + beta.Veg2*Veg2[i]
}
# Priors
alpha ~ dnorm(0,0.001)
beta.Veg ~ dnorm(0,0.001)
beta.Veg2 ~ dnorm(0,0.001)
}
"
###########################################################
# Setting up the JAGS run:
# Set up a list that contains all the necessary data
Model.Data <- list(y = Dat$Count, n.dat = nrow(Dat),
Veg = Dat$Veg, Veg2 = Dat$Veg^2)
# Specify a function to generate inital values for the parameters
inits.fn <- function() list(alpha = rnorm(1,0,1),
beta.Veg = rnorm(1,0,1),
beta.Veg2 = rnorm(1,0,1))
# Compile the model and run the MCMC for an adaptation (burn-in) phase
jagsModel <- jags.model(file= textConnection(model), data=Model.Data,
inits = inits.fn, n.chains = 3, n.adapt= 5000)Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 200
Unobserved stochastic nodes: 3
Total graph size: 1406
Initializing model# Specify parameters for which posterior samples are saved
para.names <- c('alpha','beta.Veg','beta.Veg2')
# Continue the MCMC runs with sampling
Samples <- coda.samples(jagsModel , variable.names = para.names, n.iter = 5000)
# Statistical summaries of the posterior distributions
summary(Samples)
Iterations = 5001:10000
Thinning interval = 1
Number of chains = 3
Sample size per chain = 5000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
alpha 3.2915 0.02315 0.0001890 0.0003286
beta.Veg 0.4526 0.02930 0.0002392 0.0003466
beta.Veg2 -0.8897 0.02881 0.0002352 0.0004405
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
alpha 3.2459 3.2759 3.2915 3.3070 3.3375
beta.Veg 0.3963 0.4328 0.4526 0.4723 0.5113
beta.Veg2 -0.9475 -0.9087 -0.8895 -0.8704 -0.8344# Plot MCMC samples
plot(Samples)
# Check convergence
gelman.diag(Samples)Potential scale reduction factors:
Point est. Upper C.I.
alpha 1 1
beta.Veg 1 1
beta.Veg2 1 1
Multivariate psrf
1# Correlation plot
correlationPlot(Samples)
7.2 Adding posterior predictions and residual checks
model ="
model{
# Likelihood
for(i in 1:n.dat){
# poisson model p(y|lambda)
y[i] ~ dpois(lambda[i])
# log link function
log(lambda[i]) <- mu[i]
# linear predictor on the log scale
mu[i] <- alpha + beta.Veg*Veg[i] + beta.Veg2*Veg2[i]
}
# Priors
alpha ~ dnorm(0,0.001)
beta.Veg ~ dnorm(0,0.001)
beta.Veg2 ~ dnorm(0,0.001)
# Model predictions
for(i in 1:n.pred){
y.pred[i] ~ dpois(lambda.pred[i])
log(lambda.pred[i]) <- mu.pred[i]
mu.pred[i] <- alpha + beta.Veg*Veg.pred[i] + beta.Veg2*Veg2.pred[i]
}
}
"
###########################################################
# Setting up the JAGS run:
# Set up a list that contains all the necessary data
# Note that for prediction (later used for model checking)
# we here use the original predictor variables.
Model.Data <- list(y = Dat$Count, n.dat = nrow(Dat),
Veg = Dat$Veg, Veg2 = Dat$Veg^2,
Veg.pred = Dat$Veg, Veg2.pred = Dat$Veg^2,
n.pred = nrow(Dat))
# Specify a function to generate inital values for the parameters
inits.fn <- function() list(alpha = rnorm(1,0,1),
beta.Veg = rnorm(1,0,1),
beta.Veg2 = rnorm(1,0,1))
# Compile the model and run the MCMC for an adaptation (burn-in) phase
jagsModel <- jags.model(file= textConnection(model), data=Model.Data,
inits = inits.fn, n.chains = 3, n.adapt= 5000)Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 200
Unobserved stochastic nodes: 203
Total graph size: 2007
Initializing model# Specify parameters for which posterior samples are saved
para.names <- c('alpha','beta.Veg','beta.Veg2')
# Continue the MCMC runs with sampling
Samples <- coda.samples(jagsModel , variable.names = para.names, n.iter = 5000)
# Statistical summaries of the posterior distributions
summary(Samples)
Iterations = 5001:10000
Thinning interval = 1
Number of chains = 3
Sample size per chain = 5000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
alpha 3.2909 0.02327 0.0001900 0.0003341
beta.Veg 0.4532 0.02926 0.0002389 0.0003448
beta.Veg2 -0.8889 0.02889 0.0002359 0.0004325
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
alpha 3.2449 3.2751 3.2910 3.3070 3.3364
beta.Veg 0.3962 0.4335 0.4530 0.4725 0.5119
beta.Veg2 -0.9462 -0.9079 -0.8888 -0.8695 -0.8321# Plot MCMC samples
plot(Samples)
# Check convergence
gelman.diag(Samples)Potential scale reduction factors:
Point est. Upper C.I.
alpha 1 1
beta.Veg 1 1
beta.Veg2 1 1
Multivariate psrf
1# Correlation plot
correlationPlot(Samples)
#############################################################
# Sample simulated posterior for lizard counts (y.pred)
Pred.Samples <- coda.samples(jagsModel,
variable.names = "y.pred",
n.iter = 5000)
# Transform mcmc.list object to a matrix
Pred.Mat <- as.matrix(Pred.Samples)
# Plot Model predictions against data
Pred.Q <- apply(Pred.Mat,2,quantile,prob=c(0.05,0.5,0.95))
plot(Dat$Veg, Dat$Count)
ord <- order(Dat$Veg)
lines(Dat$Veg[ord], Pred.Q['50%',ord],col='blue',lwd=2)
lines(Dat$Veg[ord], Pred.Q['5%',ord],col='blue')
lines(Dat$Veg[ord], Pred.Q['95%',ord],col='blue')
###########################################################
# Model checking with DHARMa
# Create model checking plots
res = createDHARMa(simulatedResponse = t(Pred.Mat),
observedResponse = Dat$Count,
fittedPredictedResponse = apply(Pred.Mat, 2, median),
integerResponse = T)
plot(res)
###########################################################7.3 Adding an overdispersion term
model = "
model{
# Likelihood
for(i in 1:n.dat){
# poisson model p(y|lambda)
y[i] ~ dpois(lambda[i])
# log link function
log(lambda[i]) <- mu[i] + eps[i]
# linear predictor on the log scale
mu[i] <- alpha + beta.Veg*Veg[i] + beta.Veg2*Veg2[i]
# overdispersion error terms
eps[i] ~ dnorm(0,tau.eps)
}
# Priors
alpha ~ dnorm(0,0.001)
beta.Veg ~ dnorm(0,0.001)
beta.Veg2 ~ dnorm(0,0.001)
tau.eps ~ dgamma(0.001,0.001)
# Model predictions
for(i in 1:n.pred){
y.pred[i] ~ dpois(lambda.pred[i])
log(lambda.pred[i]) <- mu.pred[i] + eps.pred[i]
mu.pred[i] <- alpha + beta.Veg*Veg.pred[i] + beta.Veg2*Veg2.pred[i]
eps.pred[i] ~ dnorm(0, tau.eps)
}
}
"
###########################################################
# Setting up the JAGS run:
# Set up a list that contains all the necessary data
# Note that for prediction (later used for model checking)
# we here use the original predictor variables.
Model.Data <- list(y = Dat$Count, n.dat = nrow(Dat),
Veg = Dat$Veg, Veg2 = Dat$Veg^2,
Veg.pred = Dat$Veg, Veg2.pred = Dat$Veg^2,
n.pred = nrow(Dat))
# Specify a function to generate inital values for the parameters
inits.fn <- function() list(alpha = rnorm(1,0,1),
beta.Veg = rnorm(1,0,1),
beta.Veg2 = rnorm(1,0,1),
tau.eps = 1
)
# Compile the model and run the MCMC for an adaptation (burn-in) phase
jagsModel <- jags.model(file= textConnection(model), data=Model.Data,
inits = inits.fn, n.chains = 3, n.adapt= 5000)Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 200
Unobserved stochastic nodes: 604
Total graph size: 3008
Initializing model# Specify parameters for which posterior samples are saved
para.names <- c('alpha','beta.Veg','beta.Veg2',
'tau.eps')
# Continue the MCMC runs with sampling
Samples <- coda.samples(jagsModel , variable.names = para.names, n.iter = 5000)
# Statistical summaries of the posterior distributions
summary(Samples)
Iterations = 5001:10000
Thinning interval = 1
Number of chains = 3
Sample size per chain = 5000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
alpha 2.9983 0.11799 0.0009634 0.009941
beta.Veg 0.6792 0.09252 0.0007554 0.004298
beta.Veg2 -1.0016 0.08074 0.0006592 0.004194
tau.eps 0.8481 0.14668 0.0011976 0.003949
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
alpha 2.7742 2.9141 3.0001 3.0800 3.2202
beta.Veg 0.5022 0.6162 0.6777 0.7399 0.8626
beta.Veg2 -1.1608 -1.0563 -1.0025 -0.9476 -0.8424
tau.eps 0.5957 0.7451 0.8353 0.9386 1.1673# Plot MCMC samples
plot(Samples)
# Check convergence
gelman.diag(Samples)Potential scale reduction factors:
Point est. Upper C.I.
alpha 1.06 1.18
beta.Veg 1.01 1.03
beta.Veg2 1.02 1.08
tau.eps 1.00 1.01
Multivariate psrf
1.05# Correlation plot
correlationPlot(Samples)
#############################################################
# Sample simulated posterior for lizard counts (y.pred)
Pred.Samples <- coda.samples(jagsModel,
variable.names = "y.pred",
n.iter = 5000)
# Transform mcmc.list object to a matrix
Pred.Mat <- as.matrix(Pred.Samples)
# Plot Model predictions against data
Pred.Q <- apply(Pred.Mat,2,quantile,prob=c(0.05,0.5,0.95))
plot(Dat$Veg, Dat$Count)
ord <- order(Dat$Veg)
lines(Dat$Veg[ord], Pred.Q['50%',ord],col='blue',lwd=2)
lines(Dat$Veg[ord], Pred.Q['5%',ord],col='blue')
lines(Dat$Veg[ord], Pred.Q['95%',ord],col='blue')
###########################################################
# Model checking with DHARMa
# Create model checking plots
res = createDHARMa(simulatedResponse = t(Pred.Mat),
observedResponse = Dat$Count,
fittedPredictedResponse = apply(Pred.Mat, 2, median),
integerResponse = T)
plot(res)
###########################################################7.4 Adding zero-inflation
model = "
model{
# Likelihood
for(i in 1:n.dat){
# poisson model p(y|lambda)
y[i] ~ dpois(lambda.eff[i])
# effective mean abundance
lambda.eff[i] <- lambda[i] * Inc[i]
# binary variable to indicate occupancy
Inc[i] ~ dbin(p.Inc,1)
# log link function
log(lambda[i]) <- mu[i] + eps[i]
# linear predictor on the log scale
mu[i] <- alpha + beta.Veg*Veg[i] + beta.Veg2*Veg2[i]
# overdispersion error terms
eps[i] ~ dnorm(0,tau.eps)
}
# Priors
alpha ~ dnorm(0,0.001)
beta.Veg ~ dnorm(0,0.001)
beta.Veg2 ~ dnorm(0,0.001)
tau.eps ~ dgamma(0.001,0.001)
p.Inc ~ dbeta(1,1)
# Model predictions
for(i in 1:n.pred){
y.pred[i] ~ dpois(lambda.eff.pred[i])
lambda.eff.pred[i] <- lambda.pred[i] * Inc.pred[i]
Inc.pred[i] ~ dbin(p.Inc, 1)
log(lambda.pred[i]) <- mu.pred[i] + eps.pred[i]
mu.pred[i] <- alpha + beta.Veg*Veg.pred[i] + beta.Veg2*Veg2.pred[i]
eps.pred[i] ~ dnorm(0, tau.eps)
}
}
"
###########################################################
# Setting up the JAGS run:
# Set up a list that contains all the necessary data
# Note that for prediction (later used for model checking)
# we here use the original predictor variables.
Model.Data <- list(y = Dat$Count, n.dat = nrow(Dat),
Veg = Dat$Veg, Veg2 = Dat$Veg^2,
Veg.pred = Dat$Veg, Veg2.pred = Dat$Veg^2,
n.pred = nrow(Dat))
# Specify a function to generate inital values for the parameters
inits.fn <- function() list(alpha = rnorm(1,0,1),
beta.Veg = rnorm(1,0,1),
beta.Veg2 = rnorm(1,0,1),
tau.eps = 1,
p.Inc = rbeta(1,1,1),
Inc = rep(1,nrow(Dat))
)
# Compile the model and run the MCMC for an adaptation (burn-in) phase
jagsModel <- jags.model(file= textConnection(model), data=Model.Data,
inits = inits.fn, n.chains = 3, n.adapt= 5000)Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 200
Unobserved stochastic nodes: 1005
Total graph size: 3810
Initializing model# Specify parameters for which posterior samples are saved
para.names <- c('alpha','beta.Veg','beta.Veg2',
'tau.eps','p.Inc')
# Continue the MCMC runs with sampling
Samples <- coda.samples(jagsModel , variable.names = para.names, n.iter = 5000)
# Statistical summaries of the posterior distributions
summary(Samples)
Iterations = 5001:10000
Thinning interval = 1
Number of chains = 3
Sample size per chain = 5000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
alpha 3.3352 0.04679 0.0003820 0.0016122
beta.Veg 0.2113 0.04998 0.0004081 0.0012269
beta.Veg2 -0.6883 0.04765 0.0003890 0.0014550
p.Inc 0.7315 0.03354 0.0002739 0.0004433
tau.eps 8.7415 1.93937 0.0158349 0.0557740
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
alpha 3.2437 3.3033 3.3357 3.3669 3.4261
beta.Veg 0.1168 0.1771 0.2104 0.2444 0.3109
beta.Veg2 -0.7829 -0.7200 -0.6877 -0.6559 -0.5962
p.Inc 0.6644 0.7090 0.7320 0.7545 0.7959
tau.eps 5.5991 7.3592 8.5155 9.8499 13.1798# Plot MCMC samples
plot(Samples)
# Check convergence
gelman.diag(Samples)Potential scale reduction factors:
Point est. Upper C.I.
alpha 1.00 1.01
beta.Veg 1.00 1.00
beta.Veg2 1.00 1.00
p.Inc 1.00 1.00
tau.eps 1.01 1.03
Multivariate psrf
1.01# Correlation plot
correlationPlot(Samples)
#############################################################
# Sample simulated posterior for lizard counts (y.pred)
Pred.Samples <- coda.samples(jagsModel,
variable.names = "y.pred",
n.iter = 5000)
# Transform mcmc.list object to a matrix
Pred.Mat <- as.matrix(Pred.Samples)
# Plot Model predictions against data
Pred.Q <- apply(Pred.Mat,2,quantile,prob=c(0.05,0.5,0.95))
plot(Dat$Veg, Dat$Count)
ord <- order(Dat$Veg)
lines(Dat$Veg[ord], Pred.Q['50%',ord],col='blue',lwd=2)
lines(Dat$Veg[ord], Pred.Q['5%',ord],col='blue')
lines(Dat$Veg[ord], Pred.Q['95%',ord],col='blue')
###########################################################
# Model checking with DHARMa
# Create model checking plots
res = createDHARMa(simulatedResponse = t(Pred.Mat),
observedResponse = Dat$Count,
fittedPredictedResponse = apply(Pred.Mat, 2, median),
integerResponse = T)
plot(res)