# R and Stan code for Bayesian Poisson regression
# Horseshoe crabs example: Section 3.3 of Agresti
# Regress Y = number of satellites on X = shell width
# Read in data
# You’ll need an internet connection for this to work
rm(list=ls()); set.seed(20210402);
filename <- "http://www.stat.ufl.edu/~aa/cat/data/Crabs.dat"
Crabs <- read.table(file=filename, header=T); rm(filename);
Crabs # The Data!
dim(Crabs); names(Crabs); # We'll only use sat & width for now
# Scatterplot of number of satellites versus width of shell
plot(jitter(sat) ~ jitter(width), data=Crabs, xlab="Width (cm)",
ylab="Number of satellites")
crabs.loess <- loess(sat ~ width, data=Crabs, family="symmetric")
range(Crabs$width)
x <- seq(21.0, 33.5, 0.1)
curve(predict(crabs.loess, data.frame(width=x)), add=T)
# Model fitting
library(rstan)
stan_model <- "
data{
int
int
vector [n] x;
}
parameters{
real alpha;
real beta;
}
transformed parameters{
vector[n] lambda;
lambda = exp(alpha + beta*x);
}
model{
# alpha ~ normal(0, 100);
# beta ~ normal(0, 100);
y ~ poisson(lambda) ;
}
”
n <- dim(Crabs)[1]; y <- Crabs$sat; x <- Crabs$width; data <- list(n=n, y=y, x=x) fit <- stan(model_code=stan_model, data=data, iter=2000, chains=5) sims <- extract(fit); names(sims); Results <- summary(fit)$summary; Results <- Results[1:2,] # Account for overdispersion by negative binomial regression ?!?!? # Do that in the next file