Quasi-likelihood models

Steve Simon

2021-09-20

With the generalized linear model, you have the option of choosing a model that works well with a variety of different distributions (normal, Poisson, Gamma, etc.). But you also have the option of working outside of these distributions. You can fit models that do not specify the distribtution, but which fit a particular transformation of the linear effect and particular mean-variance relationship.

suppressMessages(suppressWarnings(library(broom)))
suppressMessages(suppressWarnings(library(dplyr)))
suppressMessages(suppressWarnings(library(ggplot2)))
suppressMessages(suppressWarnings(library(magrittr)))

The quasi Poisson model

If you modify the variance function of Poisson regression, from

\(V(\mu)=\mu\)

to

\(V(\mu)=\phi \mu\)

and keep the link function the same, then you get a quasi Poisson model. There is no simple probability distribution that would produce this link and variance function, but it seems like a more reasonable approach, because there is more flexibility.

For technical reasons, the parameter \(\phi\) has to be greater than one, implying that the dispersion of the quasiPoisson model is greater than the dispersion of the Poissson model. The quasiPoisson model is often called an overdispersion model.

There are two ways to estimate \(\phi\). The first is very simple and direct.

\(\hat{\phi}= dispersion / df\)

The amount of overdispersion is equal to much much larger the dispersion is than the degrees of freedom.

A second approach estimates \(\phi\) iteratively along with the other terms in the regression model.

One important limitation of the quasi likelihood approach is that there is no natural way to estimate AIC (Akaike’s Information Criteria) an important statistic for comparing different models.

Start with Poisson regression

It is easiest to understand the quasiPoisson model by comparing it to the Poisson model.

Francis Huang has a [nice web page][hua1] showing two datasets and how to analyze them using the Poisson regression model as well as several more sophisticated models. The first data set is available for anyone to download.

fn <- "http://faculty.missouri.edu/huangf/data/poisson/articles.csv"
articles <- read.csv(fn)
head(articles)
##   fem      ment  phd mar kid5 art
## 1   0  7.999999 1.38   1    2   3
## 2   0  7.000000 4.29   0    0   0
## 3   0 47.000008 3.85   0    0   4
## 4   0 19.000002 3.59   1    1   1
## 5   0  0.000000 1.81   1    0   1
## 6   0  6.000000 3.59   1    1   1

The variables are + gender (fem), + marital status (mar), + number of children (kid5), + prestige of graduate program (phd), + the number of articles published by the individual’s mentor (ment) + the number of articles published by the scientist (art: the outcome).

To fit a Poisson regression model in R, use the glm function with the argument, family=“poisson”.

pois1 <- glm(art ~ fem + ment + phd + mar + kid5, family=poisson, data = articles)
summary(pois1)
## 
## Call:
## glm(formula = art ~ fem + ment + phd + mar + kid5, family = poisson, 
##     data = articles)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.5672  -1.5398  -0.3660   0.5722   5.4467  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.304617   0.102981   2.958   0.0031 ** 
## fem         -0.224594   0.054613  -4.112 3.92e-05 ***
## ment         0.025543   0.002006  12.733  < 2e-16 ***
## phd          0.012823   0.026397   0.486   0.6271    
## mar          0.155243   0.061374   2.529   0.0114 *  
## kid5        -0.184883   0.040127  -4.607 4.08e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 1817.4  on 914  degrees of freedom
## Residual deviance: 1634.4  on 909  degrees of freedom
## AIC: 3314.1
## 
## Number of Fisher Scoring iterations: 5

The interpretation of the coefficients should be done after exponentiating their values.

pois1 %>%
  tidy %>%
  mutate(lower_limit=estimate-1.96*std.error) %>%
  mutate(upper_limit=estimate+1.96*std.error) %>%
  select(term, estimate, lower_limit, upper_limit) %>%
  mutate(estimate=exp(estimate)) %>%
  mutate(lower_limit=exp(lower_limit)) %>%
  mutate(upper_limit=exp(upper_limit)) -> pois1_coefficients
pois1_coefficients
## # A tibble: 6 x 4
##   term        estimate lower_limit upper_limit
##   <chr>          <dbl>       <dbl>       <dbl>
## 1 (Intercept)    1.36        1.11        1.66 
## 2 fem            0.799       0.718       0.889
## 3 ment           1.03        1.02        1.03 
## 4 phd            1.01        0.962       1.07 
## 5 mar            1.17        1.04        1.32 
## 6 kid5           0.831       0.768       0.899

Fitting a quasi-Poisson model

quasi1 <- glm(art ~ fem + ment + phd + mar + kid5, family=quasipoisson, data = articles)
summary(quasi1)
## 
## Call:
## glm(formula = art ~ fem + ment + phd + mar + kid5, family = quasipoisson, 
##     data = articles)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.5672  -1.5398  -0.3660   0.5722   5.4467  
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.304617   0.139273   2.187 0.028983 *  
## fem         -0.224594   0.073860  -3.041 0.002427 ** 
## ment         0.025543   0.002713   9.415  < 2e-16 ***
## phd          0.012823   0.035700   0.359 0.719544    
## mar          0.155243   0.083003   1.870 0.061759 .  
## kid5        -0.184883   0.054268  -3.407 0.000686 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 1.829006)
## 
##     Null deviance: 1817.4  on 914  degrees of freedom
## Residual deviance: 1634.4  on 909  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5

Notice that there is no value listed for AIC.

quasi1 %>%
  tidy %>%
  mutate(lower_limit=estimate-1.96*std.error) %>%
  mutate(upper_limit=estimate+1.96*std.error) %>%
  select(term, estimate, lower_limit, upper_limit) %>%
  mutate(estimate=exp(estimate)) %>%
  mutate(lower_limit=exp(lower_limit)) %>%
  mutate(upper_limit=exp(upper_limit)) -> quasi1_coefficients
quasi1_coefficients
## # A tibble: 6 x 4
##   term        estimate lower_limit upper_limit
##   <chr>          <dbl>       <dbl>       <dbl>
## 1 (Intercept)    1.36        1.03        1.78 
## 2 fem            0.799       0.691       0.923
## 3 ment           1.03        1.02        1.03 
## 4 phd            1.01        0.944       1.09 
## 5 mar            1.17        0.993       1.37 
## 6 kid5           0.831       0.747       0.924

If you compare these confidence intervals to the earlier ones, they are wider. This reflects the overdispersion that is only accounted for properly by the quasi Poisson model.

References

David Lillis. Generalized Linear Models in R, Part 7: Checking for Overdispersion in Count Regression. The Analysis Factor blog. Available in html format