# R code for Poisson regression and negative binomial 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(20210401);
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 # Loglinear model fit by maximum likelihood (Poisson GLM with log link) fit.log <- glm(sat ~ width, data=Crabs, family=poisson(link="log")) summary(fit.log) curve(predict(fit.log, data.frame(width=x), type="response"), add=T, col="red") # Account for overdispersion by negative binomial regression library(MASS) fit.negbin <- glm.nb(sat ~ width, data=Crabs) summary(fit.negbin) # The theta parameter in glm.nb is the m parameter of Reich and Ghosh m.hat <- fit.negbin$theta; m.hat; # Estimate that Var = lambda + lambda^2 / 0.905