• Introduce probability model for count data.
• Interpretation of coefficients on log and non-log scale.
• Making predictions from a model in log and non-log counts.
• Running comparisons for nested models and interpreting the results.

• Generalised Linear Regression models for modelling counts as outcome variable.
• A count is the number of times something has happened or a frequency.
• Properties of counts
  • Non-negative integer (not -1, -4)
  • Does not allow decimals (1.4, -2.6)
• Events can’t happen 1.5 times (but 1.5 times on average)
• Examples?

• Number of crimes committed by someone / in a year.
• Number of times someone got sent to jail.
• Number of times someone had an affair.
• Number of words produced by a second grade writer.

In a Poisson model

\[ \begin{aligned} y_i &\sim \text{Poisson}(\lambda_i) \end{aligned} \]

we assume that \(y\) is distributed following a Poisson process.

A logarithmic link function maps the rate \(\lambda\) to a linear function of a predictor \(x\):

\[ \log(\lambda_i) = \beta_0 + \beta_1 \times x_i \]

This allows us to model rate changes relative to changes in one (more more) predictor variable(s).

Similar to binary logistic regression: probability is not a linear function of the predictor variables but log-odds are.

Number of each value of cigs: do any values stand out?

Distribution family is Poisson and link is log (default).

# Poisson model with intercept only
m_cigs_0 <- glm(cigs ~ 1,
              data = smoking_df,
              family = poisson(link = "log"))

# Return estimate for model parameters
coef(m_cigs_0)

(Intercept) 
   2.161769 

The intercept is the log number of cigarettes.

How does the number of cigs smoked vary as a function of age, education and income as predictors?

m_cigs <- glm(cigs ~ age + educ + lincome, data = smoking_df, family = poisson(link = "log"))

coef(m_cigs)

 (Intercept)          age         educ      lincome 
 0.797828756 -0.004992017 -0.045805445  0.219303695 

• For the slope estimates, notice the signs and related p-values
• Example interpretation: older individuals smoke less
• The average log number of cigarettes is lower.
• By how much does the average number of cigs smoked change for every unit change in age, educ, lincome?
• Estimates are the change in the logarithm of the mean.

What is the log mean (number of cigs) if age = 25, educ = 15, and lincome = 10?

betas <- coef(m_cigs) # assigning the output of the function to an object 'betas'
betas

 (Intercept)          age         educ      lincome 
 0.797828756 -0.004992017 -0.045805445  0.219303695 

Use the linear regression equation as usual; use indexing for the coefficients stored in betas.

What is the log mean (number of cigs) if age = 25, educ = 15, and lincome = 10?

# Saving the log average of the calculation to the object 'log_mean'
log_mean <- betas[1] + betas[2] * 25 + betas[3] * 15 + betas[4] * 10
log_mean

(Intercept) 
   2.178984 

Convert the log average number of cigs to the actual number of cigs.

For a hypothetical participant with the following characteristics, age = 18, educ = 4, lincome = 7, what is the probability that they smoke 3, 6, 10 cigarettes?

• Calculate the log average number of cigs using the regression equation.
• Take the exponential of that average.
• Used dpois to get the probability

Which number for x maximizes this probability value?

Alternative way to get these log means is to use predict or add_predictions

Let’s use three different values for age but, lets keep other variables constant.

# New dataframe with 3 values of age (21, 22, 23) but each have an educ of 15, lincome of 10
smoking_df2 <- tibble(age = c(21, 22, 23), educ = 15, lincome = 10)
smoking_df2

# A tibble: 3 × 3
    age  educ lincome
  <dbl> <dbl>   <dbl>
1    21    15      10
2    22    15      10
3    23    15      10

Predicted log of the mean values for these three different characteristics:

# Unless specified, these will be in the units of the model so logs)
library(modelr)
add_predictions(data = smoking_df2, model = m_cigs)

# A tibble: 3 × 4
    age  educ lincome  pred
  <dbl> <dbl>   <dbl> <dbl>
1    21    15      10  2.20
2    22    15      10  2.19
3    23    15      10  2.19

Predicted the mean values for these three different characteristics

# We can specify predictions to be the same unit as the data (number of cigs)
add_predictions(data = smoking_df2, model = m_cigs, type = 'response')

# A tibble: 3 × 4
    age  educ lincome  pred
  <dbl> <dbl>   <dbl> <dbl>
1    21    15      10  9.02
2    22    15      10  8.97
3    23    15      10  8.93

Implement the following two models with \(\text{cigs}_i \sim \text{Poisson}(\lambda_i)\):

\[ \mathcal{M_1}: \text{log}(\lambda_i) = \beta_0 + \beta_1 \times \log(\text{income}_i)\\ \mathcal{M_2}: \text{log}(\lambda_i) = \beta_0 + \beta_1 \times \log(\text{income}_i) + \beta_2 \times \text{education}_i + \beta_3 \times \text{white}_i \] We use white (binary), educ (as continuous variable), and the log of income (lincome).

Using anova(..., ...., test = 'Chisq') evaluate whether education and being white has a (combined) association with smoking (after controlling for income).

For binary logistic regression, we have seen that the deviance \(\mathcal{D}\) and the likelihood \(\mathcal{L}\) have the following correspondence:

\[ \mathcal{D} = -2 \times \log(\mathcal{L}) \]

This isn’t true for Poisson regression: The deviance of model m1 is

deviance(m1)

[1] 15740.23

but for the log likelihood we see that

-2 * logLik(m1) 

'log Lik.' 17210.8 (df=2)