Poisson regression

The data: Household size in the Philippines

The data fHH1.csv come from the 2015 Family Income and Expenditure Survey conducted by the Philippine Statistics Authority.

Goal: Understand the association between household size and various characteristics of the household

Response: - total: Number of people in the household other than the head

The data

hh_data <- read_csv("data/fHH1.csv")
hh_data |> slice(1:5) |> kable()

Why the least-squares model doesn’t work

The goal is to model \(\lambda\), the expected number of people in the household (other than the head), as a function of the predictors (covariates)

Poisson regression model

If \(Y_i \sim Poisson\) with \(\lambda = \lambda_i\) for the given values \(x_{i1}, \ldots, x_{ip}\), then

\[\log(\lambda_i) = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \dots + \beta_p x_{ip}\]

• Each observation can have a different value of \(\lambda\) based on its value of the predictors \(x_1, \ldots, x_p\)

• \(\lambda\) determines the mean and variance, so we don’t need to estimate a separate error term

Assumptions for Poisson regression

Model 1: Number in household vs. age

model1 <- glm(total ~ age, data = hh_data, family = poisson)

tidy(model1) |> 
  kable(digits = 4)

\[\log(\hat{\lambda}) = 1.5499 - 0.0047 ~ age\]

Answers

• Each additional year older the head of household is, the estimated average log of the number of people in the household is .0047 lower.

• Each additional year older the head of household is, the estimated average number of people in the household reduces by 0.5%.

• Each additional year older the head of household is the estimated average number of people in the household changes by a factor of .995.

• For every 10 years older the head of household is, the estimated average number of people in the household reduces by 5%.

Is the coefficient of age statistically significant?

\[H_0: \beta_1 = 0 \hspace{2mm} \text{ vs. } \hspace{2mm} H_a: \beta_1 \neq 0\]

What are plausible values for the coefficient of age?

95% confidence interval for the coefficient of age

\[\hat{\beta}_1 \pm Z^{*}\times SE(\hat{\beta}_1)\] \[-0.0047 \pm 1.96 \times 0.0009 = \mathbf{(-.0065, -0.0029)}\]

Model 2: Add a quadratic effect for age

We can determine whether to keep \(age^2\) in the model in two ways:

1️⃣ Use the p-value (or confidence interval) for the coefficient (since we are adding a single term to the model)

2️⃣ Conduct a drop-in-deviance test

Deviance

A deviance is a way to measure how the observed data deviates from the model predictions.

• It’s a measure unexplained variability in the response variable (similar to SSE in linear regression )

• Lower deviance means the model is a better fit to the data

Drop-in-Deviance Test

We can use a drop-in-deviance test to compare two models. To conduct the test

1️⃣ Compute the deviance for each model

2️⃣ Calculate the drop in deviance

\[\text{drop-in-deviance = Deviance(reduced model) - Deviance(larger model)}\]

Drop-in-deviance - Model1 and Model2

anova(model1, model2, test = "Chisq") |>
  kable(digits = 3)

Add location to the model?

Suppose we want to add location to the model, so we compare the following models:

Model A: \(\lambda_i = \beta_0 + \beta_1 ~ age_i + \beta_2 ~ age_i^2\)

Model B: \(\lambda_i = \beta_0 + \beta_1 ~ age_i + \beta_2 ~ age_i^2 + \beta_3 ~ Loc1_i + \beta_4 ~ Loc2_i + \beta_5 ~ Loc3_i + \beta_6 ~ Loc4_i\)

Question

Which of the following are reliable ways to determine if location should be added to the model?

Add location to the model?

model3 <- glm(total ~ age + age2 + location, data = hh_data, family = poisson)

anova(model2, model3, test = "Chisq") |>
  kable(digits = 3)

Model 3

Pearson residuals

We can calculate two types of residuals for Poisson regression: Pearson residuals and deviance residuals

Deviance residuals

The deviance residual indicates how much the observed data deviates from the fitted model

\[\text{deviance residual}_i = \text{sign}(y_i - \hat{\lambda}_i)\sqrt{2\Bigg[y_i\log\bigg(\frac{y_i}{\hat{\lambda}_i}\bigg) - (y_i - \hat{\lambda}_i)\Bigg]}\]

where

\[\text{sign}(y_i - \hat{\lambda}_i) = \begin{cases} 1 & \text{ if }(y_i - \hat{\lambda}_i) > 0 \\ -1 & \text{ if }(y_i - \hat{\lambda}_i) < 0 \\ 0 & \text{ if }(y_i - \hat{\lambda}_i) = 0 \end{cases}\]

Model 3: Residual plots

model3_aug_pearson <- augment(model3, type.residuals = "pearson") 
model3_aug_deviance <- augment(model3, type.residuals = "deviance")

p1 <- ggplot(data = model3_aug_pearson, aes(x = .fitted, y = .resid)) + 
  geom_point()  + 
  geom_smooth() + 
  labs(x = "Fitted values", 
       y = "Pearson residuals", 
       title = "Pearson residuals vs. fitted")

p2 <-  ggplot(data = model3_aug_deviance, aes(x = .fitted, y = .resid)) + 
  geom_point()  + 
  geom_smooth() + 
  labs(x = "Fitted values", 
       y = "Deviance residuals", 
       title = "Deviance residuals vs. fitted")

p1 + p2

Goodness-of-fit

• Goal: Use the (residual) deviance to assess how much the predicted values differ from the observed values. Recall \((\text{deviance}) = \sum_{i=1}^{n}(\text{deviance residual})_i^2\)

• If the model sufficiently fits the data, then

\[\text{deviance} \sim \chi^2_{df}\]

where \(df\) is the model’s residual degrees of freedom

Model 3: Goodness-of-fit calculations

model3$deviance

[1] 2187.8

model3$df.residual

[1] 1493

pchisq(model3$deviance, model3$df.residual, lower.tail = FALSE)

[1] 3.153732e-29

The probability of observing a deviance greater than 2187.8 is \(\approx 0\), so there is significant evidence of lack-of-fit.

Lack-of-fit

There are a few potential reasons for lack-of-fit:

• Missing important interactions or higher-order terms

• Missing important variables (perhaps this means a more comprehensive data set is required)

• There could be extreme observations causing the deviance to be larger than expected (assess based on the residual plots)

• There could be a problem with the Poisson model

  • May need more flexibility in the model to handle overdispersion

Overdispersion

Overdispersion: There is more variability in the response than what is implied by the Poisson model

hh_data |>
  summarise(mean = mean(total), var = var(total)) |>
  kable(digits = 3)

hh_data |>
  group_by(location) |>
  summarise(mean = mean(total), var = var(total)) |>
  kable(digits = 3)

Why overdispersion matters

If there is overdispersion, then there is more variation in the response than what’s implied by a Poisson model. This means

❌ The standard errors of the model coefficients are artificially small

❌ The p-values are artificially small

❌ This could lead to models that are more complex than what is needed

Dispersion parameter

The dispersion parameter is represented by \(\phi\)

\[\hat{\phi} =\frac{\text{deviance}}{\text{residual df}} = \frac{\sum_{i=1}^{n}(\text{Pearson residuals})^2}{n - p}\]

where \(p\) is the number of terms in the model (including the intercept)

• If there is no overdispersion \(\hat{\phi} = 1\)

• If there is overdispersion \(\hat{\phi} > 1\)

Accounting for dispersion in the model

We inflate the standard errors of the coefficient by multiplying the variance by \(\hat{\phi}\)

Model 3: Quasi-Poisson model

model3_q <- glm(total ~ age + age2 + location, data = hh_data, 
                family = quasipoisson) #<<

Poisson vs. Q-Poisson

Q-Poisson: Inference for coefficients

Negative binomial regression model

Another approach to handle overdispersion is to use a negative binomial regression model

• This has more flexibility than the quasi-Poisson model, because there is a new parameter in addition to \(\lambda\)

Negative binomial regression model

• Main idea: Generate a \(\lambda\) for each observation (household) and generate a count using the Poisson random variable with parameter \(\lambda\)
  • Makes the counts more dispersed than with a single parameter
• Think of it as a Poisson model such that \(\lambda\) is also random

Negative binomial regression in R

library(MASS)
model3_nb <- glm.nb(total ~ age + age2 + location, data = hh_data)
tidy(model3_nb) |> 
  kable(digits = 4)

Data: Airbnbs in NYC

The data set NYCairbnb-sample.csv contains information about a random sample of 1000 Airbnbs in New York City. It is a subset of the data on 40628 Airbnbs scraped by Awad et al. (2017).

Variables

• number_of_reviews: Number of reviews for the unit on Airbnb (proxy for number of rentals)
• price: price per night in US dollars
• room_type: Entire home/apartment, private room, or shared room
• days: Number of days the unit has been listed (date when info scraped - date when unit first listed on Airbnb)

Data: Airbnbs in NYC

airbnb <- read_csv("data/NYCairbnb-sample.csv") |>
  dplyr::select(id, number_of_reviews, days, room_type, price)

Goal: Use the price and room type of Airbnbs to describe variation in the number of reviews (a proxy for number of rentals).

EDA

p1 <- ggplot(data = airbnb, aes(x = number_of_reviews)) + 
  geom_histogram() + 
   labs(x = "Number of reviews",
    title = "Distribution of number of reviews")

p2 <- airbnb |>
  filter(price <= 2000) |>
  group_by(price) |>
  summarise(log_mean = log(mean(number_of_reviews))) |>
  ggplot(aes(x = price, y = log_mean)) + 
  geom_point(alpha= 0.7) + 
  geom_smooth() + 
  labs(x = "Price in  US dollars",
    y = "Log(mean # reviews)", 
    title = "Log mean #  of reviews vs. price", 
    subtitle = "Airbnbs $2000 or less")

p3 <- airbnb |>
  filter(price <= 500) |>
  group_by(price) |>
  summarise(log_mean = log(mean(number_of_reviews))) |>
  ggplot(aes(x = price, y = log_mean)) + 
  geom_point(alpha= 0.7) + 
  geom_smooth() + 
  labs(x = "Price in  US dollars",
    y = "Log(mean # reviews)", 
    title = "Log mean # of reviews vs. price", 
    subtitle = "Airbnbs $500 or less")

p1  / (p2 + p3) 

EDA

Considerations for modeling

We would like to fit the Poisson regression model

\[\log(\lambda) = \beta_0 + \beta_1 ~ price + \beta_2 ~ room\_type1 + \beta_3 ~ room\_type2\]

Offset

• Sometimes counts are not directly comparable because the observations differ based on some characteristic directly related to the counts, i.e. the sampling effort.

• An offset can be used to adjust for differences in sampling effort.

Adding an offset to the Airbnb model

We will add the offset \(\log(days)\) to the model. This accounts for the fact that we would expect Airbnbs that have been listed longer to have more reviews.

\[\log(\lambda_i) = \beta_0 + \beta_1 ~ price_i + \beta_2 ~ room\_type1_i + \beta_3 ~ room\_type2_i + \log(days_i)\]

Note: The response variable for the model is still \(\log(\lambda_i)\), the log mean number of reviews

Detail on the offset

We want to adjust for the number of days, so we are interested in \(\frac{reviews}{days}\).

Airbnb model in R

airbnb_model <- glm(number_of_reviews ~ price + room_type, 
                    data = airbnb, family = poisson, 
                    offset = log(days)) #<<

Interpretations

Question:

• Interpret the coefficient of price.
• Interpret the coefficient of room_typePrivate room

Goodness-of-fit

\[\begin{aligned}&H_0: \text{ The model is a good fit for the data}\\ &H_a: \text{ There is significant lack of fit}\end{aligned}\]

pchisq(airbnb_model$deviance, airbnb_model$df.residual, lower.tail = F)

[1] 0

Data: Weekend drinking

The data weekend-drinks.csv contains information from a survey of 77 students in a introductory statistics course on a dry campus.

Variables

• drinks: Number of drinks they had in the past weekend
• off_campus: 1 - lives off campus, 0 otherwise
• first_year: 1 - student is a first-year, 0 otherwise
• sex: f - student identifies as female, m - student identifies as male

EDA: Response variable

Observed vs. expected response

Two types of zeros

There are two types of zeros

• Those who happen to have a zero in the data set (people who drink but happened to not drink last weekend)
• Those who will always report a value of zero (non-drinkers)
  • These are called true zeros

Zero-inflated Poisson model

Zero-inflated Poisson (ZIP) model has two parts

Details of the ZIP model

• The ZIP model is a special case of a latent variable model
  • A type of mixture model where observations for one or more groups occur together but the group membership unknown
• Zero-inflated models are a common type of mixture model; they apply beyond Poisson regression

Tidy output

Use the tidy function from the poissonreg package for tidy model output.

library(poissonreg)

tidy(drinks_zip, type = "count") %>% kable(digits = 3)

Interpreting the model coefficients

Questions

• Interpret the intercept.
• Interpret the coefficients of off_campus and sexm.

Estimated proportion zeros

Questions:

Based on the model…

• What is the probability a first-year student is a non-drinker?
• What is the probability a upperclass student (sophomore, junior, senior) is a non-drinker?

These are just a few of the many models…

• Use the Vuong Test to compare the fit of the ZIP model to a regular Poisson model
  • Why can’t we use a drop-in-deviance test?
• We’ve just discussed the ZIP model here, but there are other model applications beyond the standard Poisson regression model (e.g., hurdle models, Zero-inflated Negative Binomial models, etc. )

Estimating coefficients in Poisson model

• Least squares estimation would not work because the normality and equal variance assumptions don’t hold for Poisson regression

• Maximum likelihood estimation is used to estimate the coefficients of Poisson regression.

• The likelihood is the product of the probabilities for the \(n\) independent observations in the data

Likelihood for regular Poisson model

Let’s go back to example about household size in the Philippines. We will focus on the model using the main effect of age to understand variability in mean household size.

Likelihood for regular Poisson model

We will use the log likelihood to make finding the MLE easier

Likelihood for regular Poisson model

Given the age of the head of the household, we fit the model

\[\log(\lambda_i) = \beta_0 + \beta_1~age_i\]

Then we replace each \(\lambda_i\) in \(\log(L)\) with \(e^{\beta_0 + \beta_1~age_i}\).

Probabilities under ZIP model

There are three different types of observations in the data:

• Observed zero and will always be 0 (true zeros)
• Observed 0 but will not always be 0
• Observed non-zero count and will not always be 0

Probabilities under ZIP model

True zeros

\[P(0 | \text{true zero})= \alpha\]

Probabilities under ZIP model

Putting this all together. Let \(y_i\) be an observed response then

\[P(Y_i = y_i | x_i) = \begin{cases} \alpha + (1 - \alpha)e^{-\lambda_i} && \text{ if } y_i = 0 \\ (1 - \alpha)\frac{e^{-\lambda_i}\lambda_i^{y_i}}{y_i!} && \text{ if } y_i > 0 \end{cases}\]