ASMbook (version 1.0.2)

demoMCMC: Fit a Poisson GLM with MCMC

Description

This is a demo function that fits a Poisson GLM with one continuous covariate to some data (y, x) using a random-walk Metropolis Markov chain Monte Carlo algorithm.

Usage

demoMCMC(
  y,
  x,
  true.vals = c(2.5, 0.14),
  inits = c(0, 0),
  prior.sd.alpha = 100,
  prior.sd.beta = 100,
  tuning.params = c(0.1, 0.1),
  niter = 10000,
  nburn = 1000,
  quiet = TRUE,
  show.plots = TRUE
)

Value

A list containing input settings, acceptance probabilities and MCMC samples.

Arguments

y

A vector of counts, e.g., y in the Swiss bee-eater example

x

A vector of a continuous explanatory variable, e.g. year x in the bee-eaters

true.vals

True intercept and slope if known (i.e., when run on simulated data)

inits

Initial values in the MCMC algorithm for alpha, beta

prior.sd.alpha

SD of normal prior for alpha

prior.sd.beta

SD of normal prior for beta

tuning.params

SD of the Gaussian proposal distributions for alpha, beta

niter

Total chain length (before burnin)

nburn

Burn-in length

quiet

Logical, suppress console output

show.plots

Logical, should diagnostic plots be shown?

Author

Marc Kéry

Examples

# Load the real data used in the publication by Mueller (Vogelwarte, 2021)

# Counts of known pairs in the country 1990-2020
y <- c(0,2,7,5,1,2,4,8,10,11,16,11,10,13,19,31,
      20,26,19,21,34,35,43,53,66,61,72,120,102,159,199)
year <- 1990:2020     # Define year range
x <- (year-1989)      # Scaled, but not centered year as a covariate
x <- x-16             # Now it's centered

oldpar <- par(no.readonly = TRUE)
par(mfrow = c(1, 2), mar = c(5,5,5,2), cex.lab = 1.5, cex.axis = 1.5, cex.main = 1.5)
plot(table(y), xlab = 'Count (y)', ylab = 'Frequency', frame = FALSE,
     type = 'h', lend = 'butt', lwd = 5, col = 'gray20', main = 'Frequency distribution of counts')
plot(year, y, xlab = 'Year (x)', ylab = 'Count (y)', frame = FALSE, cex = 1.5,
     pch = 16,  col = 'gray20', main = 'Relationship y ~ x')
fm <- glm(y ~ x, family = 'poisson')        # Add Poisson GLM line of best fit
lines(year, predict(fm, type = 'response'), lwd = 3, col = 'red', lty = 3)

# \donttest{
# Execute the function with default function args
# In a real test you should run more iterations
par(mfrow = c(1,1))
str(tmp <- demoMCMC(niter=100, nburn=50))

# Use data created above
par(mfrow = c(1,1))
str(tmp <- demoMCMC(y = y, x = x, niter=100, nburn=50))

# }
par(oldpar)