knitr::opts_chunk$set(echo = T, results = "hide")
For count response variables, the glm framework has two options. The Poisson family and the negative binomial family. In this vignette, we will consider both and learn when to use one or the other.
First lets get an idea of how the Poisson family looks. If lambda is near 1, the Poisson family is right skewed. As lambda increases, the Poisson family becomes more symmetrical.
library(stats) library(MASS) library(dplyr, warn.conflicts= FALSE) library(ggplot2) library(GlmSimulatoR) set.seed(1) #lambda = 1 poisson <- rpois(n = 10000, lambda = 1) poissonDF <- as.data.frame(x=poisson) ggplot(poissonDF, aes(x=poisson)) + geom_histogram(bins = 100) #lambda = 5 poisson <- rpois(n = 10000, lambda = 5) poissonDF <- as.data.frame(x=poisson) ggplot(poissonDF, aes(x=poisson)) + geom_histogram(bins = 100) #lambda = 10 poisson <- rpois(n = 10000, lambda = 10) poissonDF <- as.data.frame(x=poisson) ggplot(poissonDF, aes(x=poisson)) + geom_histogram(bins = 100)
rm(poisson, poissonDF)
Lets create data and train a model.
set.seed(1) simdata <- simulate_poisson(N = 10000, weights = c(.5, 1)) #Response looks similar to above histograms ggplot(simdata, aes(x=Y)) + geom_histogram(bins = 100) glmPoisson <- glm(Y ~ X1 + X2, data = simdata, family = poisson(link = "log")) summary(glmPoisson)
The estimated weights are close the the weights argument in simulate_poisson. The model estimates are accurate.
rm(simdata, glmPoisson)
One important characteristic of the Poisson family is $\mu$ equals $\sigma^2$. We see the sample statistics are nearly equal for different values of lambda.
set.seed(1) #lambda = 1 poisson <- rpois(n = 10000, lambda = 1) mean(poisson) var(poisson) #lambda = 5 poisson <- rpois(n = 10000, lambda = 5) mean(poisson) var(poisson) #lambda = 10 poisson <- rpois(n = 10000, lambda = 10) mean(poisson) var(poisson)
rm(poisson)
For modeling, the main difference between Poisson and the negative binomial is the extra parameter. The mean variance relationship is $\sigma^2$ equals $\mu + \mu^2/\theta$ for the negative binomial. Through $\theta$ many relationships are possible.
When the response variable is a count, but $\mu$ does not equal $\sigma^2$, the Poisson distribution is not applicable. Over dispersion can be detected by dividing the residual deviance by the degrees of freedom. If this quotient is much greater than one, the negative binomial distribution should be used. There is no hard cut off of "much larger than one", but a rule of thumb is 1.10 or greater is considered large.
When $\theta$ is large, $\mu + \mu^2/\theta$ is roughly equal to $\mu$ . Therefore the negative binomial will be similar to Poisson distribution.
set.seed(1) poisson <- rpois(n = 10000, lambda = 1) poissonDF <- as.data.frame(x=poisson) ggplot(poissonDF, aes(x=poisson)) + geom_histogram(bins = 100) negBin <- rnegbin(n = 10000, mu = 1, theta = 1000) negBinDF <- as.data.frame(x=negBin) ggplot(negBinDF, aes(x=negBin)) + geom_histogram(bins = 100)
rm(poisson, poissonDF, negBin, negBinDF)
If $\theta$ is small, $\mu + \mu^2/\theta$ does not approximately equal $\mu$. The negative binomial will look different than the Poisson distribution. Note the difference in scale on the x axis and y axis.
set.seed(1) poisson <- rpois(n = 10000, lambda = 1) poissonDF <- as.data.frame(x=poisson) ggplot(poissonDF, aes(x=poisson)) + geom_histogram(bins = 100) negBin <- rnegbin(n = 10000, mu = 1, theta = 1) negBinDF <- as.data.frame(x=negBin) ggplot(negBinDF, aes(x=negBin)) + geom_histogram(bins = 100)
rm(poisson, poissonDF, negBin, negBinDF)
Lets create data where the flexibility of the negative binomial is needed and compare the Poisson estimates to the negative binomial estimates.
set.seed(1) simdata <- simulate_negative_binomial(N = 10000, weights = c(.5, 1), ancillary = 5) #ancillary is theta. #Response looks like a negative binomial distribution. ggplot(simdata, aes(x=Y)) + geom_histogram(bins = 200) glmPoisson <- glm(Y ~ X1 + X2, data = simdata, family = poisson(link = "log")) glmNB <- glm.nb(Y ~ X1 + X2, data = simdata, link = "log") summary(glmPoisson) summary(glmNB)
rm(simdata, glmPoisson, glmNB)
The estimated slopes are very similar. However, the standard errors are not. The fact $\mu$ does not equal $\sigma^2$ has caused a major under estimation of standard errors for the Poisson GLM. For Poisson models, it is important to look at residual deviance divided by the degrees of freedom. When this quotient is larger than 1, the negative binomial should be considered. In the above the quotient was 6.33 (63327 / 9997) for the Poisson glm. For these data, the negative binomial model is more appropriate.
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.