• Another type of generalised linear model
• Type of count model (one of three we will be looking at)

• A regression model 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

• Number of times someone got sent to jail

• Number of times someone had an affair

• Sometimes / often people treat counts as a continuous variable

• Sometimes that’s fine (large quantities) sometimes it is not (small quantities, i.e. when events are rare)

In a Poison model

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

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

When the rate changes as a function of changes in predictor variables, we require a logarithmic link function to map the rate \(\lambda\) to a linear function of a set of predictors:

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

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

Number of each value of cigs: does any value stand out?

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

m_cigs_0 <- glm(cigs ~ 1,
              data = smoking_df,
              family = poisson(link = "log"))

coef(m_cigs_0)

(Intercept) 
   2.161769 

Hence, 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.

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

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

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

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 

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?

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})\\ \mathcal{M_2}: \text{log}(\lambda_i) = \beta_0 + \beta_1 \times \log(\text{income}) + \beta_2 \times \text{education} + \beta_3 \times \text{white} \] 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).

What is the deviance of each of the two models?

• Difference in the deviance which can be obtained using deviance() (the lower the deviance the better fit the model).
• Is their difference larger than 2 times the difference in the number of parameters between the two models?
• A difference greater than 4 usually means that one model is better than the other.
• Null hypothesis test using a \(\chi^2\) distribution which returns the p-value of pchisq()
• The null hypothesis here is that models are identical: if p is below 0.05, the null hypothesis is unlikely.

• Name of the probability model for count data.
• Making predictions from a model in log and non-log values.
• Interpretation of coefficients.
• Running and interpreting nested model comparison