demoMCMC: Fit a Poisson GLM with MCMC

In ASMbook: Functions for the Book "Applied Statistical Modeling for Ecologists"

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
)

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?

Value

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

Author(s)

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)


# 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)

ASMbook documentation built on Sept. 11, 2024, 5:38 p.m.