compare_poisson_quasipoisson_and_nb.md

# suppose we have two populations for which we'd like to infer the incidence rate ratio --
# 
# we could use a poisson, quasipoisson, or negative binomial model 
# 
# we'd like to know if we get the same estimates for incidence rate ratios using each of these models 

# dependencies
library(MASS)
library(broom)


# simulate some population sizes
population1 <- 
  sample.int(n=1000, size = 1000, replace = TRUE)

population2 <- 
  sample.int(n=1000, size = 1000, replace = TRUE)

population3 <- 
  sample.int(n=1000, size = 1000, replace = TRUE)


# specify cases
population1_cases <- 
  round(sapply(population1, function(x) sample.int(size=1, n = x, replace = FALSE)))

# population 1 and 2 are designed to have rate ratio close to 1, so we expect to see a
# coefficient on groupb == 0
population2_cases <- 
  round(sapply(population2, function(x) sample.int(size=1, n = x, replace = FALSE)))

# population 3 is designed so we should expect to see a coefficient on groupc 10
population3_cases <- 
  round(sapply(population3, function(x) sample.int(size=1, n = x, replace = FALSE))*exp(10))

# specify dataframe 
df <- data.frame(
  popsize = c(population1, population2, population3),
  group = c(rep('a', 1000), rep('b', 1000), rep('c', 1000)),
  cases = c(population1_cases, population2_cases, population3_cases)
  )

# fit models
poisson_model <- glm(cases ~ group + offset(log(popsize)), family = poisson(link=log), data = df)

quasipoisson_model <- glm(cases ~ group + offset(log(popsize)), family = quasipoisson, data = df)

nb_model <- glm.nb(cases ~ group + offset(log(popsize)), data = df)

## Warning: glm.fit: algorithm did not converge

# print model summaries
summary(poisson_model)

## 
## Call:
## glm(formula = cases ~ group + offset(log(popsize)), family = poisson(link = log), 
##     data = df)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -4192.6    -13.2     -0.2     10.1   2836.2  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.715241   0.002021 -353.87   <2e-16 ***
## groupb       0.074672   0.002792   26.75   <2e-16 ***
## groupc      10.032063   0.002021 4963.34   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 1.4568e+10  on 2999  degrees of freedom
## Residual deviance: 2.2248e+09  on 2997  degrees of freedom
## AIC: 2224801088
## 
## Number of Fisher Scoring iterations: 5

summary(quasipoisson_model)

## 
## Call:
## glm(formula = cases ~ group + offset(log(popsize)), family = quasipoisson, 
##     data = df)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -4192.6    -13.2     -0.2     10.1   2836.2  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.71524    1.62440  -0.440    0.660    
## groupb       0.07467    2.24361   0.033    0.973    
## groupc      10.03206    1.62444   6.176 7.48e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 645913.1)
## 
##     Null deviance: 1.4568e+10  on 2999  degrees of freedom
## Residual deviance: 2.2248e+09  on 2997  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5

summary(nb_model)

## 
## Call:
## glm.nb(formula = cases ~ group + offset(log(popsize)), data = df, 
##     init.theta = 1.869506856, link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -4.3216  -0.8453  -0.0005   0.6011   1.0954  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.70898    0.02338 -30.324   <2e-16 ***
## groupb       0.05202    0.03303   1.575    0.115    
## groupc      10.01177    0.03289 304.430   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(1.8695) family taken to be 1)
## 
##     Null deviance: 3200447.7  on 2999  degrees of freedom
## Residual deviance:    3265.7  on 2997  degrees of freedom
## AIC: 56913
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  1.8695 
##           Std. Err.:  0.0457 
## 
##  2 x log-likelihood:  -56905.3750

# print model coefficient tables
tidy(poisson_model)

## # A tibble: 3 × 5
##   term        estimate std.error statistic   p.value
##   <chr>          <dbl>     <dbl>     <dbl>     <dbl>
## 1 (Intercept)  -0.715    0.00202    -354.  0        
## 2 groupb        0.0747   0.00279      26.7 1.29e-157
## 3 groupc       10.0      0.00202    4963.  0

tidy(quasipoisson_model)

## # A tibble: 3 × 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)  -0.715       1.62   -0.440  6.60e- 1
## 2 groupb        0.0747      2.24    0.0333 9.73e- 1
## 3 groupc       10.0         1.62    6.18   7.48e-10

tidy(nb_model)

## # A tibble: 3 × 5
##   term        estimate std.error statistic   p.value
##   <chr>          <dbl>     <dbl>     <dbl>     <dbl>
## 1 (Intercept)  -0.709     0.0234    -30.3  5.60e-202
## 2 groupb        0.0520    0.0330      1.58 1.15e-  1
## 3 groupc       10.0       0.0329    304.   0

This seemingly suggests that the poisson, quasipoisson, and glm.nb method all give similar results for larger coefficients, though the coefficients in the example that appear to be near zero slightly differ comparing the negative binomial model to the quasipoisson and poisson model coefficient estimates.

Perhaps another way to build intuition around the Poisson and Negative Binomial distributions and their differences is to use the fitdistr function from the {MASS} package. In the following examples, we can see that the point-estimates for the rates are estimated to be exactly the same, although the associated standard-error differ between the two models.

# dependencies
library(MASS)
library(broom)
library(ggplot2)
library(dplyr)
library(magrittr)

# simulate some population sizes
N <- 100
population1 <- population2 <- population3 <- rep(N, N)

# specify cases
population1_cases <- 
  round(sapply(population1, function(x) sample.int(size=1, n = x, replace = FALSE)))

population2_cases <- 
  round(sapply(population2, function(x) sample.int(size=1, n = x, replace = FALSE))/2)

population3_cases <- 
  round(sapply(population3, function(x) sample.int(size=1, n = x, replace = FALSE))*2)

# specify dataframe 
df <- data.frame(
  popsize = c(population1, population2, population3),
  group = c(rep('a', N), rep('b', N), rep('c', N)),
  cases = c(population1_cases, population2_cases, population3_cases)
  )

# check the distribution of cases
ggplot(df, aes(x = cases, fill = group)) + 
  geom_histogram(alpha = .7) + 
  facet_grid(group~.)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# use MASS::fitdistr to estimate negative binomial and poisson fits to each group
fitdist_df <- 
  df %>% nest_by(group) %>% 
    rowwise() %>% 
    mutate(
      fitdist_dnbinom = list(fitdistr(data$cases, "Negative Binomial")),
      dnbinom = list(dnbinom(1:(2*N), # this estimates the density at x=1:(2*N) for the fit negative binomial model
          mu = fitdist_dnbinom[[1]][[2]], 
          size = fitdist_dnbinom[[1]][[1]])),
      fitdist_poisson = list(fitdistr(data$cases, "Poisson")),
       dpois = list(dpois(x = 1:(2*N), # this estimates the density at x=1:(2*N) for the fit poisson model
           lambda = fitdist_poisson[[1]][[1]]))
        )

# put the density of the fit poisson and negative binomial next to a column 'x' with values 1:(2*N) 
fitdist_df %<>% rowwise() %>% mutate(
  dnbinom = list(data.frame(x = 1:(2*N), dnbinom = dnbinom)),
  dpois = list(data.frame(x = 1:(2*N), dpois = dpois)),
)

# pull out the fits for tabling separately
rate_estimates <- fitdist_df %>% select(group, fitdist_dnbinom, fitdist_poisson)

# extract the rate estimates and their standard errors
rate_estimates %<>% rowwise() %>% mutate(
  dnbinom_mu = fitdist_dnbinom[[1]][[2]],
  dpoisson_lambda = fitdist_poisson[[1]][[1]],
  dnbinom_mu_std_error = fitdist_dnbinom[[2]][[2]],
  dpoisson_lambda_std_error = fitdist_poisson[[2]][[1]],
)

# pivot longer so that we have everything in unique rows
rate_estimates %<>% pivot_longer(
  cols = c('dnbinom_mu', 'dpoisson_lambda', 'dnbinom_mu_std_error', 'dpoisson_lambda_std_error'),
  values_to = 'rate_estimate')

# separate the variable names 
rate_estimates %<>% separate(name, into = c('family', 'term', 'type'))

# distinguish estimates from their associated standard errors
rate_estimates$type <- ifelse(is.na(rate_estimates$type), 'estimate', 'std')

# pivot wider to put the standard errors in a column next to the corresponding estimates
rate_estimates %<>% pivot_wider(
  id_cols = c('group', 'family', 'term'),
  names_from = 'type',
  values_from = 'rate_estimate')

# remove the fit objects
fitdist_df %<>% select(-fitdist_poisson, -fitdist_dnbinom, -data)

# remove the x column from one of the nested density data frames
fitdist_df %<>% mutate(
  dnbinom = list(dnbinom %>% select(-x))) 

# unnest so we have the density at each x for the fit dnbinom and dpois
fitdist_df %<>% tidyr::unnest(cols = c(dnbinom, dpois))

# visualize distributions and their point estimates for the rates
ggplot(df, aes(x = cases, fill = group)) + 
  geom_histogram(alpha = .7) + 
  geom_area(data = fitdist_df, mapping = aes(x = x, y = dnbinom * 1000, group = group), 
    alpha = .5) + 
  geom_area(data = fitdist_df, mapping = aes(x = x, y = dpois * 1000, group = group), 
    alpha = .5) + 
  geom_vline(
    data = rate_estimates, 
    mapping = aes(xintercept = estimate, color = group, linetype = family)) + 
  facet_grid(group~.) + 
  labs(linetype = "Rate Estimate by Family") + 
  ylab("Density") + 
  xlab("Observed Events") + 
  ggtitle("Comparing Poisson and Negative Binomial Fits to Observed Counts")

print(rate_estimates)

## # A tibble: 6 × 5
##   group family   term   estimate   std
##   <chr> <chr>    <chr>     <dbl> <dbl>
## 1 a     dnbinom  mu         49.7 3.47 
## 2 a     dpoisson lambda     49.7 0.705
## 3 b     dnbinom  mu         24.6 1.65 
## 4 b     dpoisson lambda     24.6 0.496
## 5 c     dnbinom  mu        104.  6.28 
## 6 c     dpoisson lambda    104.  1.02