# How to choose best multiple regression model (Poisson/quasipoisson/negative binomial)? - R

I'm creating some multiple regression models on some national statistic data checking whether there is a divide in infant deaths between the north and sound of the UK.

I have created models for males, females and all the data (males+females). This data contains deaths, population, divide(north/south) and gender per year from 1965-2016.

As an example using the female data, I at first created a Poisson model, this was the output:

glm(formula = deaths ~ Divide + year + offset(log(population)),
family = poisson(link = "log"), data = nsfemaleMerge)

Deviance Residuals:
Min       1Q   Median       3Q      Max
-6.5420  -2.1335   0.4779   2.3226   5.7547

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 64.5738825  0.3792783  170.25   <2e-16 ***
DivideSouth -0.1726417  0.0054365  -31.76   <2e-16 ***
year        -0.0348758  0.0001914 -182.23   <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: 38817.15  on 103  degrees of freedom
Residual deviance:   890.77  on 101  degrees of freedom
AIC: 1817.8

Number of Fisher Scoring iterations: 4


As you can see there is overdispersion present. I then did a quasipoisson and negative Binomial mode, quasipossion resulted in:

glm(formula = deaths ~ Divide + year + offset(log(population)),
family = quasipoisson(link = "log"), data = nsfemaleMerge)

Deviance Residuals:
Min       1Q   Median       3Q      Max
-6.5420  -2.1335   0.4779   2.3226   5.7547

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 64.5738825  1.1243161   57.43   <2e-16 ***
DivideSouth -0.1726417  0.0161158  -10.71   <2e-16 ***
year        -0.0348758  0.0005673  -61.47   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for quasipoisson family taken to be 8.787408)

Null deviance: 38817.15  on 103  degrees of freedom
Residual deviance:   890.77  on 101  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 4


and negative binomial resulted in:

Call:
glm.nb(formula = deaths ~ Divide + year + offset(log(population)),
data = nsfemaleMerge, init.theta = 142.5507992, link = log)

Deviance Residuals:
Min       1Q   Median       3Q      Max
-2.3861  -0.6426   0.1220   0.7028   2.2366

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 62.6865274  1.1806850   53.09   <2e-16 ***
DivideSouth -0.1810315  0.0177523  -10.20   <2e-16 ***
year        -0.0339236  0.0005934  -57.17   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(142.5508) family taken to be 1.023361)

Null deviance: 3655.61  on 103  degrees of freedom
Residual deviance:  104.97  on 101  degrees of freedom
AIC: 1262.2

Number of Fisher Scoring iterations: 1


I then compared the parameter estimates and standard errors:

The standard errors for quasipoisson and NB seem extremely close whereas for poisson seem to be underestimated.

I also did an ANOVA:

Analysis of Deviance Table

Model 1: deaths ~ Divide + year + offset(log(population))
Model 2: deaths ~ Divide + year + offset(log(population))
Model 3: deaths ~ Divide + year + offset(log(population))
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1       101     890.77
2       101     890.77  0     0.00
3       101     104.97  0   785.81


My question is how do I choose which model is best to use? Also how would I represent my model or visualise it say for a research paper?

I'm still getting into stats so any help would be appreciated, thanks!