#libraries
#install.packages(list("readxl", "MASS", "car", "boot", "actuar", "EnvStats", "gridExtra", "dobson", "knitr", "nnet", "questionr"))
library(caret)

## Loading required package: lattice

## Loading required package: ggplot2

library(readxl)
library(MASS)
library(car)

## Warning: package 'car' was built under R version 3.5.1

## Loading required package: carData

library(boot)

## Warning: package 'boot' was built under R version 3.5.1

## 
## Attaching package: 'boot'

## The following object is masked from 'package:car':
## 
##     logit

## The following object is masked from 'package:lattice':
## 
##     melanoma

library(actuar)

## Warning: package 'actuar' was built under R version 3.5.1

## 
## Attaching package: 'actuar'

## The following object is masked from 'package:grDevices':
## 
##     cm

library(dobson) #data sets from text

## 
## Attaching package: 'dobson'

## The following objects are masked from 'package:boot':
## 
##     aids, dogs, melanoma, remission, survival

## The following object is masked from 'package:MASS':
## 
##     housing

## The following object is masked from 'package:lattice':
## 
##     melanoma

library(gridExtra)

## Warning: package 'gridExtra' was built under R version 3.5.1

library(knitr)

## Warning: package 'knitr' was built under R version 3.5.1

library(nnet)

## Warning: package 'nnet' was built under R version 3.5.1

library(questionr)

## Warning: package 'questionr' was built under R version 3.5.1

library(purrr)

## Warning: package 'purrr' was built under R version 3.5.1

## 
## Attaching package: 'purrr'

## The following object is masked from 'package:car':
## 
##     some

## The following object is masked from 'package:caret':
## 
##     lift

library(tidyverse)

## -- Attaching packages ------------------------------------------------------------------------------------------------------------ tidyverse 1.2.1 --

## v tibble  1.4.2     v dplyr   0.7.6
## v tidyr   0.8.1     v stringr 1.3.1
## v readr   1.1.1     v forcats 0.3.0

## -- Conflicts --------------------------------------------------------------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::combine() masks gridExtra::combine()
## x dplyr::filter()  masks stats::filter()
## x dplyr::lag()     masks stats::lag()
## x purrr::lift()    masks caret::lift()
## x dplyr::recode()  masks car::recode()
## x dplyr::select()  masks MASS::select()
## x purrr::some()    masks car::some()

#paremeters
n = 1000

#data
distributions <- data_frame(
  exponential = rexp(n, rate = 1),
  gamma = rgamma(n, shape = 5),
  weibull = rweibull(n, shape = 1),
  pareto = rpareto(n, shape = 20, scale = 2),
  lognormal = exp(rnorm(n)),
  beta = rbeta(n, shape1 = 1, shape2 = 2)
  )

#functions
dist_plot <- function(dist_name){
  
  p1 <-   distributions %>% 
          ggplot(aes(get(dist_name), fill = "area")) + 
          geom_density() + 
          ggtitle(paste0("Empirical PDF: ", dist_name)) + 
          guides(fill=FALSE) + 
          xlab("range of X")
  
  p2 <-   distributions %>% 
          ggplot(aes(get(dist_name))) + 
          stat_ecdf() + 
          ggtitle(paste0("Empirical CDF: ", dist_name))
  
  grid.arrange(p1, p2, nrow = 1)
}

feature_plot <- function(x_input, y_input){
  transparentTheme(trans = .9)
  caret::featurePlot(x = x_input,
                     y = y_input,
            plot = "pairs",
            scales = list(x = list(relation="free"), 
                          y = list(relation="free")), 
            adjust = 1.5, 
            pch = "+",
            auto.key = list(columns = 2))
}

Background

Commonly-used Response Distributions

Exponential Disribution

Memoryless
Poisson interarrival times
Moment Generating Function = lamda/(lamda - t)

dist_plot("exponential")

Gamma Distribution

Sum of Exponentials
Poisson arrival times, totals for exponential claims
MGF = (1/(1 - theta*t)^alpha) for t < 1/theta

dist_plot("gamma")

Weibull Distribution

Power Transformation of Expontential: X^(1/tao)
MGF is too complicated

dist_plot("weibull")

Pareto Distribution

Only finite moments up to alpha, and no MGF beyond
Thick-tailed in the probability of extreme values drops slowly
Hazard function is alpha(x + theta)
Mixture of Exponential Distributions

dist_plot("pareto")

Lognormal Distribution

Exponential Transformation of Normal: e^z
Right-Skewed

dist_plot("lognormal")

Beta Distribution

Transformation of X1 = Gamma(alpha, 1), X2 = Gamma (beta, 1): X1/(X1 + X2)

dist_plot("beta")

Example 1: Continuous Response with ANOVA

Data

The data is the percentages of total calories obtained from complex carbohydrates, for twenty male insulin-dependent diabetics who had been on a high-carbohydrate diet for six months. The question we are trying to answer is whether or not age plays a role in the percent of total calories obtained from complex carbs.

data("carbohydrate")
head(carbohydrate)

## # A tibble: 6 x 4
##   carbohydrate   age weight protein
##          <dbl> <dbl>  <dbl>   <dbl>
## 1           33    33    100      14
## 2           40    47     92      15
## 3           37    49    135      18
## 4           27    35    144      12
## 5           30    46    140      15
## 6           43    52    101      15

We can test models with and without age included to see if it plays a significant role.

m0 = glm(carbohydrate ~ age + weight + protein, 
         family = gaussian, data = carbohydrate)
m0

## 
## Call:  glm(formula = carbohydrate ~ age + weight + protein, family = gaussian, 
##     data = carbohydrate)
## 
## Coefficients:
## (Intercept)          age       weight      protein  
##     36.9601      -0.1137      -0.2280       1.9577  
## 
## Degrees of Freedom: 19 Total (i.e. Null);  16 Residual
## Null Deviance:       1093 
## Residual Deviance: 567.7     AIC: 133.7

We notice that the significance level for age is suspiciously high at 0.31389. We would expect to see the residual deviance to be less than the degrees of freedom, but 567.66 is much greater than 16.

summary(m0)

## 
## Call:
## glm(formula = carbohydrate ~ age + weight + protein, family = gaussian, 
##     data = carbohydrate)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -10.3424   -4.8203    0.9897    3.8553    7.9087  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept) 36.96006   13.07128   2.828  0.01213 * 
## age         -0.11368    0.10933  -1.040  0.31389   
## weight      -0.22802    0.08329  -2.738  0.01460 * 
## protein      1.95771    0.63489   3.084  0.00712 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 35.47893)
## 
##     Null deviance: 1092.80  on 19  degrees of freedom
## Residual deviance:  567.66  on 16  degrees of freedom
## AIC: 133.67
## 
## Number of Fisher Scoring iterations: 2

Another way of seeing this is by looking at the decrease in deviance by adding age to the model. As seen in the ANOVA table below, this decrease is small, which again suggests that it does not need to be included.

anova(m0)

## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: carbohydrate
## 
## Terms added sequentially (first to last)
## 
## 
##         Df Deviance Resid. Df Resid. Dev
## NULL                       19    1092.80
## age      1     3.82        18    1088.98
## weight   1   183.98        17     905.00
## protein  1   337.34        16     567.66

We fit a new model without age.

m1 = glm(carbohydrate ~ weight + protein, 
         family = gaussian, data = carbohydrate)
summary(m1)

## 
## Call:
## glm(formula = carbohydrate ~ weight + protein, family = gaussian, 
##     data = carbohydrate)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -10.6812   -3.9135    0.9464    4.0880    9.7948  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept) 33.13032   12.57155   2.635  0.01736 * 
## weight      -0.22165    0.08326  -2.662  0.01642 * 
## protein      1.82429    0.62327   2.927  0.00941 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 35.64835)
## 
##     Null deviance: 1092.80  on 19  degrees of freedom
## Residual deviance:  606.02  on 17  degrees of freedom
## AIC: 132.98
## 
## Number of Fisher Scoring iterations: 2

anova(m1)

## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: carbohydrate
## 
## Terms added sequentially (first to last)
## 
## 
##         Df Deviance Resid. Df Resid. Dev
## NULL                       19    1092.80
## weight   1   181.38        18     911.42
## protein  1   305.40        17     606.02

The AIC for the model without age is also lower.

summary <- data_frame(
  Model = c("carbohydrates ~ age + weight + protein", "carbohydrates ~ weight + protein"),
  AIC = c(AIC(m0), AIC(m1))
    )

kable(summary)

Model	AIC
carbohydrates ~ age + weight + protein	133.6734
carbohydrates ~ weight + protein	132.9812

We can compute the AIC manually using the AIC formula -2log_likelihood + 2p where p is the number of parameters being estimated. In this case, that is 4 + 1 including the estimate for sigma^2.

logLik(m1)

## 'log Lik.' -62.49061 (df=4)

-2*logLik(m1) + 2*4

## 'log Lik.' 132.9812 (df=4)

Finally, we can use an ANOVA test to compare the model with age included to the model without it.

H0: age should be included H1: age should not be included

Reject Region: p < 0.05

anova(m0, m1, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: carbohydrate ~ age + weight + protein
## Model 2: carbohydrate ~ weight + protein
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1        16     567.66                     
## 2        17     606.02 -1  -38.359   0.2984

We see that the p-value 0.2984 > 0.05 indicates that we are unable to reject the null hypothesis and say that age should be included.

glm.diag.plots(m1)

We can check for multi-colinearity by looking at the variance inflation factors (VIF). Because all of these values are close to 1, age, weight, and protein are not significantly correlated.

vif(m0)

##      age   weight  protein 
## 1.044506 1.027723 1.065696

Example 2: Continuous Response with ANCOVA

ANCOVA is the same as ANOVA with the addition of categerorical variables. For simplicity, let’s create a new category and add it to the data. Imagine that there is a variable exercise, indicating whether or not the male patient tested participated in at least 60 minutes of exercise during each week. This could be related to carbohydrate in that patients who exercised ate more carbs in order to stay energized.

carbohydrate_w_exercise <- carbohydrate %>% 
  mutate(exercise = ifelse(carbohydrate > 40, yes = "Active", no = "Inactive" ))

m2 <- glm(carbohydrate ~ weight + protein + exercise, 
         family = gaussian, 
         data = carbohydrate_w_exercise)

summary(m2)

## 
## Call:
## glm(formula = carbohydrate ~ weight + protein + exercise, family = gaussian, 
##     data = carbohydrate_w_exercise)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -8.135  -1.600   0.222   1.982   5.356  
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       39.74001    8.10837   4.901 0.000160 ***
## weight            -0.11689    0.05686  -2.056 0.056499 .  
## protein            1.09331    0.42197   2.591 0.019699 *  
## exerciseInactive -10.12876    1.98870  -5.093 0.000109 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 14.44968)
## 
##     Null deviance: 1092.80  on 19  degrees of freedom
## Residual deviance:  231.19  on 16  degrees of freedom
## AIC: 115.71
## 
## Number of Fisher Scoring iterations: 2

The p-value on the new exercise field of 0.000109 indicates that this should be included in the model. In other words, we can say that exercise plays a role in determining the response variable. And this should be the case, given that we designed it this way !

anova(m2, test = "F")

## Analysis of Deviance Table
## 
## Model: gaussian, link: identity
## 
## Response: carbohydrate
## 
## Terms added sequentially (first to last)
## 
## 
##          Df Deviance Resid. Df Resid. Dev      F    Pr(>F)    
## NULL                        19    1092.80                     
## weight    1   181.38        18     911.42 12.552 0.0027066 ** 
## protein   1   305.40        17     606.02 21.135 0.0002974 ***
## exercise  1   374.83        16     231.19 25.940 0.0001085 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Example 3: Binary Response on Grouped Data

Instead of having individual records, we have groups of records with sample sizes equal to the value of n below.

#data.  This will need to be moved to google sheets or some other online data storage location
beetles_raw <- read_csv("//FILE-NA1-02/USERDATA2$/sam82554/Desktop/MAS-I/R/TIA Data/beetle_mortality.csv")

## Parsed with column specification:
## cols(
##   dose = col_double(),
##   number = col_integer(),
##   killed = col_integer()
## )

head(beetles_raw)

## # A tibble: 6 x 3
##    dose number killed
##   <dbl>  <int>  <int>
## 1  1.69     59      6
## 2  1.72     60     13
## 3  1.76     62     18
## 4  1.78     56     28
## 5  1.81     63     52
## 6  1.84     59     53

Because we have different sample sizes at each dosage level, we need to predict the number that die and survive at each dosage level. This allows for the groups with more beetles in total to have both more killed as well as more survivors. An alternative way of modeling this would be to look at the percentage of survival as killed/total, but I’ll leave it as the author has it set up originally.

beetles <- beetles_raw %>% 
  mutate(alive = number - killed)

y <- beetles %>% 
  select(killed, alive) %>% 
  as.matrix()

x <- beetles %>% 
  select(dose) %>% 
  unlist()

m0 <- glm(y ~ x,
          family = binomial(link = "logit"))

m0

## 
## Call:  glm(formula = y ~ x, family = binomial(link = "logit"))
## 
## Coefficients:
## (Intercept)            x  
##      -60.72        34.27  
## 
## Degrees of Freedom: 7 Total (i.e. Null);  6 Residual
## Null Deviance:       284.2 
## Residual Deviance: 11.23     AIC: 41.43

We see below that while the coefficients are significant, the fit could be better. There is still variance that is not being explained as we can see from the residual deviance of 11.232 as compared to 6 degrees of freedom. Ideally, this would be less than the degrees of freedom.

summary(m0)

## 
## Call:
## glm(formula = y ~ x, family = binomial(link = "logit"))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.5941  -0.3944   0.8329   1.2592   1.5940  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -60.717      5.181  -11.72   <2e-16 ***
## x             34.270      2.912   11.77   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 284.202  on 7  degrees of freedom
## Residual deviance:  11.232  on 6  degrees of freedom
## AIC: 41.43
## 
## Number of Fisher Scoring iterations: 4

When switch to a probit link function we see that the standard error decreases, which is good, but the residual deviance still high.

m1 <- glm(y ~ x,
          family = binomial(link = "probit"))
summary(m1)

## 
## Call:
## glm(formula = y ~ x, family = binomial(link = "probit"))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.5714  -0.4703   0.7501   1.0632   1.3449  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -34.935      2.648  -13.19   <2e-16 ***
## x             19.728      1.487   13.27   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 284.20  on 7  degrees of freedom
## Residual deviance:  10.12  on 6  degrees of freedom
## AIC: 40.318
## 
## Number of Fisher Scoring iterations: 4

The complimentary log log function (AKA the extreme value tolerance distribution) is a better fit because the reisual deviance is lower than the mean degrees of freedom, the AIC is lower,and the p-values are still low.

m2 <- glm(y ~ x,
          family = binomial(link = "cloglog"))
summary(m2)

## 
## Call:
## glm(formula = y ~ x, family = binomial(link = "cloglog"))
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.80329  -0.55135   0.03089   0.38315   1.28883  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -39.572      3.240  -12.21   <2e-16 ***
## x             22.041      1.799   12.25   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 284.2024  on 7  degrees of freedom
## Residual deviance:   3.4464  on 6  degrees of freedom
## AIC: 33.644
## 
## Number of Fisher Scoring iterations: 4

We can compute the Deviance by taking the sum of the squared deviance residuals:

sum(residuals(m2, type = "deviance")^2)

## [1] 3.446439

as well as the Pearson Chi-Square statistic:

sum(residuals(m2, type = "pearson")^2)

## [1] 3.294694

The graph below shows how cloglog fit is closer on average to the true number of beetles killed. This means that the residual is smaller.

beetles %>%
  mutate(`logit fitted` = number*(m0$fitted.values),
         `complimentary log log fitted` = number*(m2$fitted.values)) %>% #multiplye percentage of killed times sample size to get number killed
  gather(key = obs_type, #name of column which will have rows of the old column names
         value = killed, #name of the value column which will have been stacked
         - c(dose, number, alive)) %>% #columns which you DON'T want to stack
  ggplot(aes(dose, killed, shape = obs_type, color = obs_type)) + 
  geom_point() + 
  ggtitle("Actual Verses Fitted on Grouped Binary Response")

## Warning: attributes are not identical across measure variables;
## they will be dropped

We can compute the pseudo R^2 value by first computing the null model, that is, the model with only an intercept term.

m_null = glm(y ~ 1, family = binomial(link="logit"))
1 - logLik(m2)/logLik(m_null)

## 'log Lik.' 0.904496 (df=2)

Let’s look again at the diagnostic plots.

glmdiag <- glm.diag(m2)   #Creating diagnostic information
glm.diag.plots(m2, glmdiag)

A final note from the author: The model could be computed on individual observations instead of grouped data. For either case, the coefficients of the model will come out the same; however the diagnostic statistics for the model will be different. These differences relate to the models varying ability to predict group rates of events as compared to its ability to predict individual outcomes.

Example 4: Nominal Logistic Regression

We look at a data set of car preferences to determine the importance of power steering and air conditioning by sex. The response is the level of importance of power steering as no/little importance, important, or very important.

car_preferences <- read_csv("//FILE-NA1-02/USERDATA2$/sam82554/Desktop/MAS-I/R/TIA Data/car_preferences.csv") %>% 
  modify_if(is.character, as.factor)

## Parsed with column specification:
## cols(
##   sex = col_character(),
##   age = col_character(),
##   age_num = col_integer(),
##   response = col_character(),
##   freq = col_integer()
## )

glimpse(car_preferences)

## Observations: 18
## Variables: 5
## $ sex      <fct> women, women, women, women, women, women, women, wome...
## $ age      <fct> 18-23, 18-23, 18-23, 24-40, 24-40, 24-40, > 40, > 40,...
## $ age_num  <int> 0, 0, 0, 1, 1, 1, 2, 2, 2, 0, 0, 0, 1, 1, 1, 2, 2, 2
## $ response <fct> no/little, important, very important, no/little, impo...
## $ freq     <int> 26, 12, 7, 9, 21, 15, 5, 14, 41, 40, 17, 8, 17, 15, 1...

summary(car_preferences)

##     sex       age       age_num            response      freq      
##  men  :9   > 40 :6   Min.   :0   important     :6   Min.   : 5.00  
##  women:9   18-23:6   1st Qu.:0   no/little     :6   1st Qu.: 9.75  
##            24-40:6   Median :1   very important:6   Median :15.00  
##                      Mean   :1                      Mean   :16.67  
##                      3rd Qu.:2                      3rd Qu.:17.75  
##                      Max.   :2                      Max.   :41.00

car_preferences %>% 
  map(~unique(.x))

## $sex
## [1] women men  
## Levels: men women
## 
## $age
## [1] 18-23 24-40 > 40 
## Levels: > 40 18-23 24-40
## 
## $age_num
## [1] 0 1 2
## 
## $response
## [1] no/little      important      very important
## Levels: important no/little very important
## 
## $freq
##  [1] 26 12  7  9 21 15  5 14 41 40 17  8 18

We set the reference levels

car_preferences <- within(car_preferences, sex <- relevel(sex, ref = "women"))
car_preferences <- within(car_preferences, response <- relevel(response, ref = "no/little"))
car_preferences <- within(car_preferences, age <- relevel(age, ref = "18-23"))

The summary below is difficult to interpret

m0 <- multinom(response ~ age + sex, weights = freq, data =  car_preferences)

## # weights:  15 (8 variable)
## initial  value 329.583687 
## iter  10 value 290.490920
## final  value 290.351098 
## converged

summary(m0)

## Call:
## multinom(formula = response ~ age + sex, data = car_preferences, 
##     weights = freq)
## 
## Coefficients:
##                (Intercept)  age> 40 age24-40     sexmen
## important       -0.5907992 1.587709 1.128268 -0.3881301
## very important  -1.0390726 2.916757 1.478104 -0.8130202
## 
## Std. Errors:
##                (Intercept)   age> 40  age24-40    sexmen
## important        0.2839756 0.4028997 0.3416449 0.3005115
## very important   0.3305014 0.4229276 0.4009256 0.3210382
## 
## Residual Deviance: 580.7022 
## AIC: 596.7022

So we can instead look at the odds ratios. What this shows us is that the odds of importance is going up with age. As the age increases from 18 - 20 to 24-40, the odds ratio increases by 3.09. As this increases to > 40, the odds ratio increases to 4.9. The relative importance for power steering was less for men than women as the odds ratios for very important/sexmen and important/sexmen are less than 1.

odds.ratio(m0)

##                                  OR    2.5 %  97.5 %         p    
## important/(Intercept)       0.55388  0.31747  0.9664 0.0374836 *  
## important/age> 40           4.89253  2.22118 10.7766 8.124e-05 ***
## important/age24-40          3.09030  1.58195  6.0368 0.0009584 ***
## important/sexmen            0.67832  0.37639  1.2225 0.1965078    
## very important/(Intercept)  0.35378  0.18510  0.6762 0.0016670 ** 
## very important/age> 40     18.48125  8.06742 42.3378 5.327e-12 ***
## very important/age24-40     4.38463  1.99832  9.6206 0.0002272 ***
## very important/sexmen       0.44352  0.23640  0.8321 0.0113262 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

We can compare the fitted probabilities with the expeted probabilities to compute a pearson statistic.

m_null <- multinom(response ~ 1, weights = freq, data =  car_preferences)

## # weights:  6 (2 variable)
## initial  value 329.583687 
## final  value 329.272024 
## converged

summary(m_null)

## Call:
## multinom(formula = response ~ 1, data = car_preferences, weights = freq)
## 
## Coefficients:
##                (Intercept)
## important      -0.11066559
## very important -0.03883986
## 
## Std. Errors:
##                (Intercept)
## important        0.1419933
## very important   0.1393729
## 
## Residual Deviance: 658.544 
## AIC: 662.544

The likelihood ratio chi-square statistic is large, indicating that the fitted m0 is not explaining a lot of the variance in the data.

2*(logLik(m0) - logLik(m_null))

## 'log Lik.' 77.84185 (df=8)

Example 5: Cumulative Odds Model

To simplify, we encode Not Important, Important, Very Important, as 1, 2, or 3. Instead of modeling the probability of being in each of the categories, the proportional odds model looks at the cumulative probabilities of case 1, case 1 OR case 2, or case 1 OR case 2 OR case 3.

We notice that the AIC is improved over the nominal model. The deviance has increased slightly due to the the decrease in the number of parameters.

car_preferences %>% 
  mutate(response_case = case_when(
    response %% "no/little " == 0 ~ 1,
    response %% "important" == 0 ~ 2,
    response %% "very important" == 0 ~ 3,
    TRUE ~ as.double(response)
  ))

## Warning in Ops.factor(response, "no/little "): '%%' not meaningful for
## factors

## Warning in Ops.factor(response, "important"): '%%' not meaningful for
## factors

## Warning in Ops.factor(response, "very important"): '%%' not meaningful for
## factors

## # A tibble: 18 x 6
##    sex   age   age_num response        freq response_case
##    <fct> <fct>   <int> <fct>          <int>         <dbl>
##  1 women 18-23       0 no/little         26             1
##  2 women 18-23       0 important         12             2
##  3 women 18-23       0 very important     7             3
##  4 women 24-40       1 no/little          9             1
##  5 women 24-40       1 important         21             2
##  6 women 24-40       1 very important    15             3
##  7 women > 40        2 no/little          5             1
##  8 women > 40        2 important         14             2
##  9 women > 40        2 very important    41             3
## 10 men   18-23       0 no/little         40             1
## 11 men   18-23       0 important         17             2
## 12 men   18-23       0 very important     8             3
## 13 men   24-40       1 no/little         17             1
## 14 men   24-40       1 important         15             2
## 15 men   24-40       1 very important    12             3
## 16 men   > 40        2 no/little          8             1
## 17 men   > 40        2 important         15             2
## 18 men   > 40        2 very important    18             3

m_pro_odds <- polr(response ~ age + sex, 
                   weights = freq, 
                   data = car_preferences)

summary(m_pro_odds)

## 
## Re-fitting to get Hessian

## Call:
## polr(formula = response ~ age + sex, data = car_preferences, 
##     weights = freq)
## 
## Coefficients:
##            Value Std. Error t value
## age> 40   2.2325     0.2915   7.659
## age24-40  1.1471     0.2776   4.132
## sexmen   -0.5762     0.2262  -2.548
## 
## Intercepts:
##                          Value   Std. Error t value
## no/little|important       0.0435  0.2323     0.1874
## important|very important  1.6550  0.2556     6.4744
## 
## Residual Deviance: 581.2956 
## AIC: 591.2956

Example 6: Count Data

Data is deaths as related collorary failure for smoking and non-smoking patients. We are looking for a relationship between smoking and deaths.

smoking_death <- read_csv("//FILE-NA1-02/USERDATA2$/sam82554/Desktop/MAS-I/R/TIA Data/smoking_death.csv") %>% 
  modify_if(is.character, as.factor)

## Parsed with column specification:
## cols(
##   age = col_character(),
##   age_num = col_integer(),
##   smoking = col_character(),
##   deaths = col_integer(),
##   `person-years` = col_integer()
## )

glimpse(smoking_death)

## Observations: 10
## Variables: 5
## $ age            <fct> 35 to 44, 45 to 54, 55 to 64, 65 to 74, 75 to 8...
## $ age_num        <int> 1, 2, 3, 4, 5, 1, 2, 3, 4, 5
## $ smoking        <fct> smoker, smoker, smoker, smoker, smoker, non-smo...
## $ deaths         <int> 32, 104, 206, 186, 102, 2, 12, 28, 28, 31
## $ `person-years` <int> 52407, 43248, 28612, 12663, 5317, 18790, 10673,...

We can look at the box plot and see tha average number of deaths for smoking is higher than that without. There is little correlation with age. We also see that non-smokers are younger on average than smokers.

featurePlot(x = smoking_death %>% select( -smoking, - age, - age_num),
            y = smoking_death$smoking,
            plot = "box",
            scales = list(y = list(relation="free"),
                          x = list(rot = 90)),
            layout = c(4,1 ), 
            auto.key = list(columns = 2)
)

smoking_death %>% 
  mutate(percent_coronary_failure =deaths/`person-years`) %>% 
  ggplot(aes(age, percent_coronary_failure)) + 
  geom_boxplot() + 
  ggtitle("Percent of coronary by age increases quantratically")

We can include a quadratic term in the model to take this into account. We also will include an interaction with age and smoking. We will use a Poisson family due to the fact that we are dealing with count data. We also include an offset term of person-years, which forces the coefficient to be 1.

m0 <- glm(deaths ~ age_num + I(age_num^2) + smoking + age_num*smoking + offset(log(`person-years`)), 
          family = "poisson", 
          data = smoking_death)

The summary below shows a great fit. The AIC is low, the residual deviance is less than the degrees of freedom, all of the variables are highly-significant.

summary(m0)

## 
## Call:
## glm(formula = deaths ~ age_num + I(age_num^2) + smoking + age_num * 
##     smoking + offset(log(`person-years`)), family = "poisson", 
##     data = smoking_death)
## 
## Deviance Residuals: 
##        1         2         3         4         5         6         7  
##  0.43820  -0.27329  -0.15265   0.23393  -0.05700  -0.83049   0.13404  
##        8         9        10  
##  0.64107  -0.41058  -0.01275  
## 
## Coefficients:
##                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           -10.79176    0.45008 -23.978  < 2e-16 ***
## age_num                 2.37648    0.20795  11.428  < 2e-16 ***
## I(age_num^2)           -0.19768    0.02737  -7.223 5.08e-13 ***
## smokingsmoker           1.44097    0.37220   3.872 0.000108 ***
## age_num:smokingsmoker  -0.30755    0.09704  -3.169 0.001528 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 935.0673  on 9  degrees of freedom
## Residual deviance:   1.6354  on 5  degrees of freedom
## AIC: 66.703
## 
## Number of Fisher Scoring iterations: 4

We can inspect the actual verses fitted values.

p1 <- 
smoking_death %>% 
  mutate(fitted = fitted.values(m0)) %>% 
  filter(smoking == "smoker") %>% 
  rename(actual = deaths) %>% 
  select(age, actual, fitted) %>% 
  gather(obs_type, deaths, - age) %>% 
  ggplot(aes(age, deaths, colour = obs_type, shape = obs_type)) + 
  geom_point() + 
  ggtitle("Actual Vs. Fitted for Smokers")

p2 <- 
  smoking_death %>% 
  mutate(fitted = fitted.values(m0)) %>% 
  filter(smoking == "non-smoker") %>% 
  rename(actual = deaths) %>% 
  select(age, actual, fitted) %>% 
  gather(obs_type, deaths, - age) %>% 
  ggplot(aes(age, deaths, colour = obs_type, shape = obs_type)) + 
  geom_point() + 
  ggtitle("Actual Vs. Fitted for Non-Smokers")

grid.arrange(p1, p2, ncol = 2)

The exponentials of the coefficients shows the relativities. The signs of the coefficients on smoker below indicates that the likelihood of dying due to collorary failure for smokers is higher than for non-smokers, but this impact decreases as age increases because the interaction of smoker and age is less than 1. The other coefficients can be interpreted similarly.

m0 %>% coef() %>% exp()

##           (Intercept)               age_num          I(age_num^2) 
##          2.056824e-05          1.076692e+01          8.206353e-01 
##         smokingsmoker age_num:smokingsmoker 
##          4.224800e+00          7.352475e-01

Example 7: Continency Tables with Poisson Regression

The data for has information relating aspirin usage to ulcers. The ulcer field is the type of ulcer, the casecontrol field indicates if the patient was tested in the case group, those with known ulcers, and control group individuals who were similar to the case group but not known to have an ulcer, and aspirin indicating whether or not a patient uses aspirin regularly. The frequency column indicates the number of patients in each of these groups. Just like in the previous example, this means that we are dealing with “compressed” or summarized data instead of line-items.

aspirin_ulcers <- read_csv("//FILE-NA1-02/USERDATA2$/sam82554/Desktop/MAS-I/R/TIA Data/aspirin_ulcers.csv") %>% 
  modify_if(is.character, as.factor)

## Parsed with column specification:
## cols(
##   ulcer = col_character(),
##   casecontrol = col_character(),
##   aspirin = col_character(),
##   frequency = col_integer()
## )

head(aspirin_ulcers)

## # A tibble: 6 x 4
##   ulcer    casecontrol aspirin  frequency
##   <fct>    <fct>       <fct>        <int>
## 1 gastric  control     non-user        62
## 2 gastric  control     user             6
## 3 gastric  case        non-user        39
## 4 gastric  case        user            25
## 5 duodenal control     non-user        53
## 6 duodenal control     user             8

The goal is to be able to predict the number of patients in each of these categories based on the 3 binary variables. This means that we have a total of 2^3 = 8 possible combinations. Which of these should we include? To start, we look only at the interaction with ulcer and casecontrol.

We see that the model is poor given the high deviance and p-values.

m0 <- glm(frequency ~ ulcer + casecontrol + ulcer*casecontrol, 
          family = "poisson", 
          data = aspirin_ulcers)
summary(m0)

## 
## Call:
## glm(formula = frequency ~ ulcer + casecontrol + ulcer * casecontrol, 
##     family = "poisson", data = aspirin_ulcers)
## 
## Deviance Residuals: 
##      1       2       3       4       5       6       7       8  
##  4.301  -5.932   1.196  -1.287   3.684  -4.857   3.480  -4.547  
## 
## Coefficients:
##                                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                      3.349904   0.132453  25.291   <2e-16 ***
## ulcergastric                     0.115832   0.182123   0.636    0.525    
## casecontrolcontrol               0.067823   0.184221   0.368    0.713    
## ulcergastric:casecontrolcontrol -0.007198   0.253512  -0.028    0.977    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 127.75  on 7  degrees of freedom
## Residual deviance: 126.71  on 4  degrees of freedom
## AIC: 174.32
## 
## Number of Fisher Scoring iterations: 5

Just adding aspirin improves the fit.

m1 <- glm(frequency ~ aspirin + ulcer + casecontrol + ulcer*casecontrol, 
          family = "poisson", 
          data = aspirin_ulcers)
summary(m1)

## 
## Call:
## glm(formula = frequency ~ aspirin + ulcer + casecontrol + ulcer * 
##     casecontrol, family = "poisson", data = aspirin_ulcers)
## 
## Deviance Residuals: 
##       1        2        3        4        5        6        7        8  
##  0.8952  -2.1191  -1.8828   3.2603   0.4872  -1.0836   0.3954  -0.8691  
## 
## Coefficients:
##                                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                      3.834796   0.135904  28.217   <2e-16 ***
## aspirinuser                     -1.463058   0.161872  -9.038   <2e-16 ***
## ulcergastric                     0.115832   0.182123   0.636    0.525    
## casecontrolcontrol               0.067823   0.184221   0.368    0.713    
## ulcergastric:casecontrolcontrol -0.007198   0.253512  -0.028    0.977    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 127.749  on 7  degrees of freedom
## Residual deviance:  21.789  on 3  degrees of freedom
## AIC: 71.404
## 
## Number of Fisher Scoring iterations: 5

m2 <- glm(frequency ~ aspirin + ulcer + casecontrol + ulcer*casecontrol + aspirin*casecontrol, 
          family = "poisson", 
          data = aspirin_ulcers)
summary(m2)

## 
## Call:
## glm(formula = frequency ~ aspirin + ulcer + casecontrol + ulcer * 
##     casecontrol + aspirin * casecontrol, family = "poisson", 
##     data = aspirin_ulcers)
## 
## Deviance Residuals: 
##       1        2        3        4        5        6        7        8  
##  0.1766  -0.5251  -1.1381   1.6950  -0.1879   0.5191   1.1388  -2.1123  
## 
## Coefficients:
##                                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                      3.724598   0.143676  25.924  < 2e-16 ***
## aspirinuser                     -0.980829   0.204122  -4.805 1.55e-06 ***
## ulcergastric                     0.115832   0.182123   0.636  0.52477    
## casecontrolcontrol               0.271396   0.194885   1.393  0.16374    
## ulcergastric:casecontrolcontrol -0.007198   0.253511  -0.028  0.97735    
## aspirinuser:casecontrolcontrol  -1.125046   0.348984  -3.224  0.00127 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 127.749  on 7  degrees of freedom
## Residual deviance:  10.538  on 2  degrees of freedom
## AIC: 62.153
## 
## Number of Fisher Scoring iterations: 4

As we add additional interactions, the number of parameters increases, which means that the degrees of freedom, 8 - p, where p is the number of coefficients, decreases. If we wanted to have perfect predictions on the training data, we would just include 8 parameters in order to have 1 prediction per row in the data.

m3 <- glm(frequency ~ aspirin + ulcer + casecontrol + ulcer*casecontrol + aspirin*casecontrol + aspirin*ulcer, 
          family = "poisson", 
          data = aspirin_ulcers)
summary(m3)

## 
## Call:
## glm(formula = frequency ~ aspirin + ulcer + casecontrol + ulcer * 
##     casecontrol + aspirin * casecontrol + aspirin * ulcer, family = "poisson", 
##     data = aspirin_ulcers)
## 
## Deviance Residuals: 
##       1        2        3        4        5        6        7        8  
##  0.4486  -1.2085  -0.5393   0.7280  -0.4661   1.4673   0.5073  -1.0829  
## 
## Coefficients:
##                                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                      3.81846    0.14515  26.307  < 2e-16 ***
## aspirinuser                     -1.37910    0.29514  -4.673 2.97e-06 ***
## ulcergastric                    -0.06977    0.20415  -0.342  0.73254    
## casecontrolcontrol               0.21517    0.19172   1.122  0.26174    
## ulcergastric:casecontrolcontrol  0.10574    0.26147   0.404  0.68590    
## aspirinuser:casecontrolcontrol  -1.14288    0.35207  -3.246  0.00117 ** 
## aspirinuser:ulcergastric         0.70005    0.34603   2.023  0.04306 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 127.749  on 7  degrees of freedom
## Residual deviance:   6.283  on 1  degrees of freedom
## AIC: 59.898
## 
## Number of Fisher Scoring iterations: 4

When we add in the last interaction term, we create what is known as a saturated model. This has deviance equal to zero (1.7764e-14 below), because deviance is the difference in log likelihood between the saturated model and the model being tested.

m4 <- glm(frequency ~ aspirin + ulcer + casecontrol + ulcer*casecontrol + aspirin*casecontrol + aspirin*casecontrol*ulcer, 
          family = "poisson", 
          data = aspirin_ulcers)
summary(m4)

## 
## Call:
## glm(formula = frequency ~ aspirin + ulcer + casecontrol + ulcer * 
##     casecontrol + aspirin * casecontrol + aspirin * casecontrol * 
##     ulcer, family = "poisson", data = aspirin_ulcers)
## 
## Deviance Residuals: 
## [1]  0  0  0  0  0  0  0  0
## 
## Coefficients:
##                                             Estimate Std. Error z value
## (Intercept)                                  3.89182    0.14286  27.243
## aspirinuser                                 -1.81238    0.38132  -4.753
## ulcergastric                                -0.22826    0.21459  -1.064
## casecontrolcontrol                           0.07847    0.19818   0.396
## ulcergastric:casecontrolcontrol              0.38510    0.28469   1.353
## aspirinuser:casecontrolcontrol              -0.07847    0.53784  -0.146
## aspirinuser:ulcergastric                     1.36769    0.45940   2.977
## aspirinuser:ulcergastric:casecontrolcontrol -1.81222    0.73329  -2.471
##                                             Pr(>|z|)    
## (Intercept)                                  < 2e-16 ***
## aspirinuser                                 2.01e-06 ***
## ulcergastric                                 0.28747    
## casecontrolcontrol                           0.69214    
## ulcergastric:casecontrolcontrol              0.17614    
## aspirinuser:casecontrolcontrol               0.88400    
## aspirinuser:ulcergastric                     0.00291 ** 
## aspirinuser:ulcergastric:casecontrolcontrol  0.01346 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 1.2775e+02  on 7  degrees of freedom
## Residual deviance: 1.7764e-14  on 0  degrees of freedom
## AIC: 55.615
## 
## Number of Fisher Scoring iterations: 3

The AIC continues to decrease as we add more terms. But this does not mean that the model would improve at making predictions on the test set, but only on the training set.

grid.arrange(
  data_frame(
    number_of_parameters = c(4, 5, 6, 7, 8),
    AIC = list(m0, m1, m2, m3, m4) %>% map_dbl(AIC)
    ) %>% 
    ggplot(aes(number_of_parameters, AIC, label = round(AIC, 1))) + 
    geom_point() + 
    geom_line() + 
    geom_text(nudge_y = 2, nudge_x = 0.1)
  ,
  data_frame(
    number_of_parameters = c(4, 5, 6, 7, 8),
    Deviance = list(m0, m1, m2, m3, m4) %>% map_dbl(deviance)
    ) %>% 
    ggplot(aes(number_of_parameters, Deviance, label = round(Deviance, 1))) + 
    geom_point() + 
    geom_line() + 
    geom_text(nudge_y = 2, nudge_x = 0.1),
  ncol = 2
  )

The third model seems to be the best trade-off in terms of decreasing the error while not adding too many parameters, frequency ~ aspirin + ulcer + casecontrol + ulcer*casecontrol + aspirin*casecontrol, where the above shows the decrease in AIC from 71.4 to 62.2.

This is an example of when we need to rely on statistics for model selection as we do not see a great fit from just looking at the actual verses fitted values.

aspirin_ulcers %>% 
  mutate(fitted = fitted.values(m2)) %>% 
  rename(actual = frequency) %>% 
  gather(obs_type, frequency, - ulcer, - casecontrol, - aspirin) %>% 
  group_by(obs_type) %>% 
  mutate(index = row_number()) %>% 
  ungroup() %>% 
  ggplot(aes(index, frequency, colour = obs_type, shape = obs_type)) + 
  geom_point() + 
  ggtitle("Actual Vs. Fitted still shows an imperfect match")

## Warning: attributes are not identical across measure variables;
## they will be dropped

Example 7: Overdispersion in Count Data

For a Poisson distribution, one of the nice properties is that the mean is equal to the variance. In Poisson regression, in order to fit a model, we need the empirical mean to be close to the empirical variance. Overdispersion is when this is not the case, and special treatment is needed.

The data set contains information for third party claims by geography for 176 different locations. Each location has a number of claims, population, population density, and number of accidents. We apply log transforms in order to aid modeling.

third_party_claims <- read_csv("//FILE-NA1-02/USERDATA2$/sam82554/Desktop/MAS-I/R/TIA Data/third_party_claims.csv") %>% 
  modify_if(is.character, as.factor) %>% 
  mutate(
    log_claims = log(claims),
    log_accidents = log(accidents),
    log_population = log(population)
  )

## Parsed with column specification:
## cols(
##   lga = col_character(),
##   sd = col_integer(),
##   claims = col_integer(),
##   accidents = col_integer(),
##   ki = col_integer(),
##   population = col_integer(),
##   pop_density = col_double()
## )

head(third_party_claims)

## # A tibble: 6 x 10
##   lga          sd claims accidents    ki population pop_density log_claims
##   <fct>     <int>  <int>     <int> <int>      <int>       <dbl>      <dbl>
## 1 ASHFIELD      1   1103      2304   920     124850      0.499        7.01
## 2 AUBURN        1   1939      2660  1465     143500      0.148        7.57
## 3 BANKSTOWN     1   4339      7381  3864     470700      0.205        8.38
## 4 BAULKHAM~     1   1491      3217  1554     311300      0.0259       7.31
## 5 BLACKTOWN     1   3801      6655  4175     584900      0.0812       8.24
## 6 BOTANY        1    387      2013   854     106350      0.178        5.96
## # ... with 2 more variables: log_accidents <dbl>, log_population <dbl>

The graph below shows a linear relation ship between the log transforms of accidents and claims.

third_party_claims %>% 
  ggplot(aes(log_claims, log_accidents)) + 
  geom_point()

We try a Poisson model, but the residual deviance is large.

m0 <- glm(claims ~ log_accidents, 
          family = "poisson",
          offset = log_population,
          data = third_party_claims) 

summary(m0)

## 
## Call:
## glm(formula = claims ~ log_accidents, family = "poisson", data = third_party_claims, 
##     offset = log_population)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -38.957   -3.551    0.116    3.842   45.965  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -7.093809   0.026992 -262.81   <2e-16 ***
## log_accidents  0.259103   0.003376   76.75   <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: 22393  on 175  degrees of freedom
## Residual deviance: 15837  on 174  degrees of freedom
## AIC: 17066
## 
## Number of Fisher Scoring iterations: 4

From Baysian theory, the negative binomial is related to the Poisson-Gamma mixture, which has a large variance than a poisson model itself. When we switch this response distribution to the negative binomial, we see a better fit in that the AIC decreases significantly. Notice that we know have a dispersion parameter.

m1 <- glm.nb(claims ~ log_accidents + offset(log_population),
             data = third_party_claims)
summary(m1)

## 
## Call:
## glm.nb(formula = claims ~ log_accidents + offset(log_population), 
##     data = third_party_claims, init.theta = 5.830937458, link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.5448  -0.8172  -0.1964   0.4260   3.7295  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -6.95443    0.15837  -43.91   <2e-16 ***
## log_accidents  0.25389    0.02472   10.27   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(5.8309) family taken to be 1)
## 
##     Null deviance: 298.16  on 175  degrees of freedom
## Residual deviance: 192.33  on 174  degrees of freedom
## AIC: 2041.3
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  5.831 
##           Std. Err.:  0.671 
## 
##  2 x log-likelihood:  -2035.255

We look at the rate of claims against the log of the accidents and see that the negative binomial fit (in red) does a slightly better job of capturing the higher rates.

plot(claims/population ~ log_accidents, data = third_party_claims, pch=16, cex=.8, las=1, cex.axis=1.1, cex.lab=1.1)
curve(exp(-6.95443 +0.25389*x), add=TRUE, lwd=4, col = "red")
curve(exp(-7.09381 + 0.25910*x), add=TRUE, lwd=3, col = "blue")

An alternative way of dealing with overdispersion is through quais-likelihood. This results in the same coefficient estimates, but they have larger variances. So you get the same estimates for the Poisson mean, but with larger predicted variances.

m2 <- glm(claims ~ log_accidents, 
          family = quasi(link="log",variance="mu"),
          offset = log_population,
          data = third_party_claims)
summary(m2)

## 
## Call:
## glm(formula = claims ~ log_accidents, family = quasi(link = "log", 
##     variance = "mu"), data = third_party_claims, offset = log_population)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -38.957   -3.551    0.116    3.842   45.965  
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -7.09381    0.27223 -26.058  < 2e-16 ***
## log_accidents  0.25910    0.03405   7.609 1.66e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasi family taken to be 101.7172)
## 
##     Null deviance: 22393  on 175  degrees of freedom
## Residual deviance: 15837  on 174  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4

This gives us the same model as the Poisson model, only the variances are now larger. The means of the response stays the same, but the variances increases according to the dispersion parameter of 101.71 as seen below. Also note that the AIC cannot be computed as we do not have a real likelihood value.

Where does R come up with the dispersion parameter? From taking Pearson Chi-Square statistic and dividing by the degrees of freedom. This is just an estimate of the variance.

sum(residuals(m2, type = "pearson")^2)/174

## [1] 101.7168

Sources:

Regression on count data: http://data.princeton.edu/wws509/notes/c4a.pdf

Linear Models Part I

Samuel Castillo

August 18, 2018