This glossary collects the key concepts, methods and people in this lesson. It can be used as a reference while working through the material or as a review before assessments. Typing in the search box filters the entries.

Lesson 4 · Section 1

Analyzing Different Outcome Types Using GLMs

Every outcome has a type, and the type of the outcome decides which regression model to use.

Variable types

Two families of variables

Numeric

Continuous variables can take any value in a range, such as blood pressure or BMI.

Discrete variables take whole numbers only, such as the number of clinic visits.

Categorical

Nominal variables have groups with no order, such as usual place of care.

Ordinal variables have groups with a natural order, such as self-rated health.

Binary variables have exactly two groups, such as hypertension (yes or no).

Running example

Five outcomes from the PHAA survey

Review

Linear regression draws a straight line

Linear model for systolic blood pressure

\[ \color{#0B7B6B}{\text{SBP}} = \color{#6D28D9}{\beta_0} + \color{#1D4ED8}{\beta_1}\,\text{age} + \color{#C2410C}{\beta_2}\,\text{smoker} + \varepsilon \]

SBP outcome (mmHg) β₀ intercept β₁ change per year of age β₂ smokers minus non-smokers

0.45

mmHg higher per extra year of age

6.06

mmHg higher in smokers than non-smokers of the same age

Assumptions

Linear regression rests on four assumptions (LINE)

L: Linearity

The average outcome changes in a straight line with each continuous predictor.

I: Independence

Each person's outcome is unrelated to every other person's outcome.

N: Normality

The residuals follow a roughly bell-shaped (normal) distribution.

E: Equal variance

The residuals are spread about equally at every predicted value.

A fifth check looks for influential points that pull the line toward themselves.

Model checking

plot(lin) draws the four checks

Limits of the straight line

Other outcomes break these assumptions

Binary

The outcome is 0 or 1. A straight line can predict values below 0 or above 1, and the residuals cannot be normal.

Counts

Counts cannot be negative, pile up near zero, and spread out more as the average rises.

Categories

Ordinal and nominal groups have no numeric spacing, so a slope per one-unit step has no meaning.

Generalized linear models

A GLM keeps the line and changes two things

The GLM recipe

\[ \color{#C2410C}{g}\big(\color{#0B7B6B}{\text{average of } Y}\big) = \color{#6D28D9}{\beta_0 + \beta_1 X_1 + \beta_2 X_2} \]

g link function average of Y the mean or probability being modelled β₀ + β₁X₁ + … linear predictor

1. Distribution (family)

The distribution describes the kind of values the outcome takes, such as normal, binomial, Poisson or negative binomial.

2. Link function

The link function transforms the average, for example with the log or the logit (log-odds).

Model menu

One outcome type, one model

Outcome type	Model	R function	eβ is a…
Continuous	Linear regression	lm()	(not exponentiated) difference in means
Binary	Logistic regression	glm(family = binomial)	odds ratio
Ordinal	Ordinal logistic regression	MASS::polr()	odds ratio
Nominal	Multinomial logistic regression	nnet::multinom()	relative risk ratio
Count	Poisson regression	glm(family = poisson)	incidence rate ratio
Overdispersed count	Negative binomial regression	MASS::glm.nb()	incidence rate ratio
Time to event	Cox model (optional walkthrough)	survival::coxph()	hazard ratio

Review

Logistic regression and its checks

2.65

odds ratio for smokers (95% CI 1.04 to 6.30)

1.05

odds ratio per year of age (95% CI 1.01 to 1.08)

• The outcome has exactly two categories, and observations are independent.
• Each continuous predictor has a straight-line relationship with the log-odds.
• There are about 10 events per predictor (23 events here support about two predictors).

Carry forward

What to take into the next section

• The type of the outcome is identified first: continuous, binary, ordinal, nominal or count.
• Every model has assumptions, and each assumption has a check that is run after fitting.
• A GLM keeps the linear predictor and changes the distribution and the link, which sets the meaning of e^β.

Introduction and Overview

Regression models describe how an outcome changes with one or more predictors. Lesson 3 introduced two of them: linear regression for a continuous outcome and logistic regression for a binary (yes-or-no) outcome. Many outcomes in public health are neither continuous nor binary. Self-rated health is recorded in ordered categories, the place a person usually goes for care is a choice among unordered options, and the number of emergency visits is a count. This lesson presents one model for each of these outcome types. This first section sets out the types of variables, reviews linear regression and the assumptions on which it rests, and introduces generalized linear models (GLMs), the family to which every model in the lesson belongs.

Box 4.1 reviews four statistical ideas that are needed to read the results in this lesson: confidence intervals, p-values, probability and odds, and exponentiated coefficients. It ends with a retrieval question about an odds ratio that is reported later in this section.

Every model in this lesson is chosen according to the type of its outcome. The first step is therefore to identify the type of a variable, which the next part sets out.

Types of Variables

Every variable in a dataset belongs to one of two families. Numeric variables record an amount. They are continuous when any value in a range is possible (blood pressure, body mass index, age measured precisely) and discrete when only whole numbers are possible (the number of visits, cases or falls). Categorical variables place each person into a group. They are nominal when the groups have no natural order (usual place of care, type of cancer, province of residence) and ordinal when the groups have a natural order (self-rated health, pain rated none to severe, stage of disease). A categorical variable with exactly two groups is called binary.

Table 4.1 lists the five types with the family to which each belongs and an example of each from the PHAA data, the survey used throughout this lesson.

Table 4.1. The five types of variable, grouped into two families, with an example of each from the PHAA data.

The type of a variable is decided by what its values mean, and the way it is stored in a file can mislead. Self-rated health might be stored as the numbers 1 to 5, yet it is still ordinal, because the distance between "good" and "very good" is unknown. A postal code is made of characters and numbers, yet it is nominal. The cards below describe each type in more detail.

Worked Example 4.1 applies these definitions to the outcomes of four studies planned by a health unit.

As Worked Example 4.1 shows, the type of the outcome decides which model is suitable. The next part reviews the model for a continuous outcome, linear regression, and the assumptions on which it rests.

Linear Regression and Its Assumptions

Linear regression models the average of a continuous outcome as a straight-line function of the predictors. For systolic blood pressure (SBP) in the PHAA survey, with age and smoking as predictors, the model is written as in Equation 4.1.

When this model is fitted to the 795 people with complete data, the slope for age is 0.45 mmHg per year (95% CI 0.40 to 0.49) and the difference for smokers is 6.06 mmHg (95% CI 4.33 to 7.79). Activity 4.1 at the end of this section produces these numbers. The estimates can be trusted only if the assumptions below are reasonable. Each assumption is described with the way it is checked in R and with a simple response if it fails.

Figure 4.1 shows the four plots that plot(lin) draws for the blood pressure model, so that each of the checks described above can be seen on real output.

For the blood pressure model, the plots in Figure 4.1 support each of the assumptions that a plot can check. The next part turns to outcomes for which these assumptions fail.

Why Some Outcomes Need a Different Model

The same four plots look very different when a straight line is fitted to an outcome that is a count. Figure 4.2 shows the first two plots for a linear regression of urgent_visits (the number of urgent care visits in the past year) on age and smoking.

Each outcome type breaks the assumptions in its own way. A binary outcome takes only the values 0 and 1, so a straight line can predict impossible probabilities below 0 or above 1, and the residuals can never be normal. A count cannot be negative, often has many zeros, and becomes more variable as its average increases, which breaks the equal-variance assumption. An ordinal outcome has an order but no known spacing between its categories, so a slope "per category" has no clear meaning. A nominal outcome has no order at all, so any numbers attached to its categories are arbitrary labels.

Linear regression is therefore suited to a continuous outcome, and each of the other four types of outcome calls for a model that matches the values the outcome can take. The next part introduces the family of models that provides one for each type.

Generalized Linear Models

A generalized linear model (GLM) extends linear regression so that it can handle these outcomes. The extension was set out by Nelder and Wedderburn (1972) and is now the standard framework for regression in epidemiology. Every GLM has three parts. The distribution (also called the family) describes the kind of values the outcome can take. The linear predictor is the familiar sum β₀ + β₁X₁ + β₂X₂ + …, which combines the predictors exactly as in linear regression. The link function is a transformation of the average outcome that is set equal to the linear predictor.

Equation 4.2 shows how the link function connects the average outcome to the linear predictor.

The link function keeps predictions in a sensible range. A probability must lie between 0 and 1, and the logit link guarantees that it does. An expected count must be positive, and the log link guarantees that it is. The link also decides how coefficients are read. When the link is a log or a logit, the coefficient β is on a log scale, so it is exponentiated (e^β, written exp(coef(model)) in R) to turn it into a ratio: an odds ratio, a relative risk ratio or a rate ratio. A ratio of 1 means no association, a ratio above 1 means the outcome is more likely or more frequent, and a ratio below 1 means it is less likely or less frequent.

Table 4.2 lists each model used in this lesson with its type of outcome, its distribution and link, the reading of e^β, and the R function that fits it. The note beneath the table explains why its last row, the Cox model, is included.

Table 4.2. The regression models in this lesson, with the distribution, link, reading of e^β and R function of each.

Strictly, the ordinal and multinomial models are close relatives of the GLM family (they are often called multivariate GLMs), because they model several probabilities at once. They are fitted and read in the same way, and this course treats them as members of the family. The Cox model in the last row is not a GLM, because it leaves the shape of the baseline hazard unspecified. It is listed because a time-to-event outcome is one more outcome type, and it is covered only in the optional Survival Analysis in R walkthrough linked below.

Table 4.2 pairs each type of outcome with a model. The next part turns this pairing into a short guide for choosing a model.

Choosing a Model

The choice of model follows from the type of the outcome, and the predictors play no part in it. The guide in Figure 4.3 can be worked through for any new outcome. The first question is whether the outcome is numeric or categorical, and the second question narrows the choice to one model.

The last row of the numeric column in Figure 4.3 leads to the Cox model, which the note to Table 4.2 places outside the GLM family. Box 4.2 defines the terms used for time-to-event outcomes, and the card that follows it links to the optional Survival Analysis in R walkthrough.

With a model matched to each type of outcome, the rest of this section reviews logistic regression, the model for a binary outcome that Lesson 3 introduced. Sections 2, 3 and 4 then take up the ordinal, nominal and count models in turn.

Logistic Regression: A Brief Review

Logistic regression is the GLM for a binary outcome. It models the log of the odds of the event, and its exponentiated coefficients are odds ratios (Lesson 3 covers it in detail). In the PHAA data, the odds of hypertension are 2.65 times as high for smokers as for non-smokers of the same age (95% CI 1.04 to 6.30), and they are 1.05 times as high (about 5% higher) for each additional year of age (95% CI 1.01 to 1.08). Its assumptions differ from those of linear regression, because a binary outcome has no normal residuals and no constant spread to check.

Table 4.3 lists the assumptions of logistic regression with the way each is checked.

Table 4.3. The assumptions of logistic regression and how each is checked.

The rule of about 10 events per predictor in Table 4.3 matters in the PHAA data, because few participants have hypertension. Box 4.3 explains how the small number of events affects the confidence interval for smoking.

Once a model has been fitted and checked, its results are usually reported in a table. Box 4.4 shows how the output of the logistic model becomes a results table of the kind found in published papers.

This section has matched each type of outcome to a model and reviewed the two models from Lesson 3. The two R activities below put this review into practice. Activity 4.1 fits the linear model for systolic blood pressure and draws the plots in Figure 4.1, and Activity 4.2 fits the logistic model for hypertension and produces the adjusted odds ratios in Table 4.4. A knowledge check and a reflection close the section.

# install.packages(c("MASS", "nnet", "brant"))   # run once, if not yet installed
phaa <- read.csv("phaa_outcomes.csv")   # file must be in the working directory
phaa <- na.omit(phaa)                   # keep the 795 people with no missing values
phaa$smoker <- factor(phaa$smoker, levels = c("No", "Yes"))   # "No" = reference
nrow(phaa)                   # how many people are in the analysis?
summary(phaa$systolic_bp)    # a first look at the outcome

lin <- lm(systolic_bp ~ age + smoker, data = phaa)   # fit the linear model
summary(lin)                 # coefficients, standard errors and p-values
confint(lin)                 # 95% confidence intervals

par(mfrow = c(2, 2))         # show four plots on one screen
plot(lin)                    # the four standard checks of a linear model
par(mfrow = c(1, 1))         # back to one plot per screen

phaa <- read.csv("phaa_outcomes.csv")   # file must be in the working directory
phaa <- na.omit(phaa)                   # keep the 795 people with no missing values
phaa$smoker <- factor(phaa$smoker, levels = c("No", "Yes"))   # "No" = reference
phaa$hypertension <- factor(phaa$hypertension, levels = c("No", "Yes"))   # "Yes" is the event
table(phaa$hypertension)     # check: how many people had the event?

logit <- glm(hypertension ~ age + smoker, data = phaa, family = binomial)
summary(logit)
exp(cbind(OR = coef(logit), confint(logit)))   # odds ratios with 95% CIs

Lesson 4 · Section 2

Ordinal Logistic Regression

A model for outcomes whose categories have a natural order, such as self-rated health from poor to excellent.

Ordinal outcomes

Ordered categories with unknown spacing

Self-rated health runs from poor to excellent.

Pain is rated none, mild, moderate or severe.

Agreement runs from strongly disagree to strongly agree.

Cancer is staged from I to IV.

The order of the categories is known. The size of the step from one category to the next is unknown and may differ from step to step.

Two shortcuts

Why not linear or binary logistic regression?

Code as 1 to 5 and use lm()

This assumes every step is the same size and that the residuals are normal. Neither assumption fits five ordered categories.

Collapse to two groups and use glm()

This throws away information and depends on where the cut-off is placed.

Ordinal logistic regression

This model uses all of the categories and their order, and it makes no assumption about the size of each step.

The idea

Four splits, one odds ratio per predictor

Split 1PoorFairGoodVery goodExcellent

Split 2PoorFairGoodVery goodExcellent

Split 3PoorFairGoodVery goodExcellent

Split 4PoorFairGoodVery goodExcellent

below the splitat or above the split

Each split compares at or above with below. The model gives each split its own intercept and gives each predictor one odds ratio for all four splits.

Proportional odds

The model in symbols

Ordinal (proportional-odds) logistic regression

\[ \ln\frac{\color{#0B7B6B}{P(Y \ge j)}}{\color{#C2410C}{P(Y < j)}} = \color{#6D28D9}{\beta_{0j}} + \color{#1D4ED8}{\beta_1}\,\text{age} + \color{#1D4ED8}{\beta_2}\,\text{smoker} \]

P(Y ≥ j) at or above category j P(Y < j) below category j β_0j one intercept per split β₁, β₂ one slope per predictor, shared by all splits

R fits this model with polr() from the MASS package, and exp(coef(model)) gives an odds ratio for being in a higher category.

The PHAA example

Self-rated health by smoking

Reading the output

Odds ratios from polr()

0.50

odds ratio for smokers (95% CI 0.36 to 0.70)

0.71

odds ratio for 10 years of age (0.966 per year)

Smokers have about half the odds of reporting a higher health category than non-smokers of the same age, at every split of the scale.

Assumptions and checks

What to check in an ordinal model

Assumption	How to check it in R
The outcome has a genuine order	Reasoning about the categories
Observations are independent	Study design
Enough people in each category	table(outcome)
No strong multicollinearity	Correlations among the predictors
Proportional odds	brant(model) from the brant package

The Brant test

Does one odds ratio fit every split?

Test for	χ2	df	p-value	Assumption holds?
Omnibus (all predictors)	6.71	6	0.35	Yes
Age	4.30	3	0.23	Yes
Smoking	2.19	3	0.53	Yes

If a p-value were below 0.05, the options would be to combine sparse categories, to fit a multinomial model (Section 3), or to fit a partial proportional-odds model.

Carry forward

What to take into the next section

• Ordinal logistic regression uses the order of the categories by fitting all of the at-or-above splits at once.
• Each predictor has one odds ratio, read as the odds of being in a higher category.
• The shared odds ratio is the proportional-odds assumption, and brant() checks it.

Introduction and Overview

This section presents ordinal logistic regression, the model for an outcome with three or more ordered categories. The running example is self-rated health, a widely used single question in health surveys that asks people to rate their own health as poor, fair, good, very good or excellent. Self-rated health predicts later illness and death even after other measures of health are taken into account (Idler & Benyamini, 1997), which makes it a common outcome in public health research. In the PHAA data, self-rated health is a simulated variable added for this lesson. The section fits the model with age and smoking as predictors and checks its assumptions.

The section begins by setting out what makes an outcome ordinal and why the two common ways of simplifying such an outcome lose information.

What Makes an Outcome Ordinal?

An outcome is ordinal when its categories can be ranked but the distances between them are unknown. Self-rated health, pain severity, agreement scales, disease stage and education level are all ordinal. The defining feature is that the order carries information: a move from "good" to "very good" is a move in a known direction. What the order does not say is how large each move is. The step from "poor" to "fair" might reflect a much larger change in health than the step from "very good" to "excellent", and the data offer no way to measure the difference.

Two shortcuts are common in practice, and both lose something. Coding the categories as 1 to 5 and fitting a linear regression treats every step as equal and assumes normal residuals, which a five-category outcome cannot have. Collapsing the categories into two groups and fitting a logistic regression discards the detail within each group and makes the result depend on where the cut-off is placed. Ordinal logistic regression avoids both problems.

The three cards below set the two shortcuts beside ordinal logistic regression, so that what each approach assumes and what it loses can be compared.

Ordinal logistic regression keeps every category and uses their order. The next part shows how it does this by fitting a series of splits of the outcome.

How the Model Works: A Series of Splits

An outcome with five ordered categories can be split into "at or above" and "below" in four places. Each split is a yes-or-no comparison that a binary logistic regression could handle on its own. The ordinal model fits the four comparisons together.

Figure 4.4 shows the four splits for self-rated health.

The model gives each split its own intercept, because the proportion of people at or above the split differs from one split to the next: almost everyone is "fair or better", while only a small share is "excellent". Each predictor, however, receives a single coefficient that applies to every split. This is why the model is called the proportional-odds model, and it is the assumption that the Brant test checks later in this section (McCullagh, 1980).

Equation 4.3 writes the model for self-rated health with age and smoking as predictors.

Box 4.5 explains how R writes this model and what that means for the sign of each coefficient and for the intercepts.

Worked Example 4.2 calculates the crude odds ratio for smoking at each of the four splits in Figure 4.4, so that the single odds ratio estimated by the model can be set beside the split-specific values.

The next part reads the single odds ratio for smoking, and the rest of the output, from the model fitted to the PHAA data.

Interpreting the Results

Fitted to the 795 people in the PHAA data, the ordinal model gives an odds ratio for smoking of 0.50 (95% CI 0.36 to 0.70). Smokers have about half the odds of non-smokers of the same age of reporting a higher category of self-rated health, and the model applies this odds ratio at every split. The odds ratio for age is 0.966 per year (95% CI 0.957 to 0.975). Because one year is a small step, the odds ratio for ten years is easier to picture: 0.966¹⁰ = 0.71, so a person ten years older has about 29% lower odds of reporting a higher category. In R, the ten-year odds ratio is exp(10 * coef(ord)["age"]).

The output of polr() also lists four intercepts, labelled Poor|Fair, Fair|Good, Good|Very good and Very good|Excellent. These are the cut-points between neighbouring categories. They are needed to calculate predicted probabilities, but they are not usually reported or interpreted. The output gives t values but no p-values. The 95% confidence intervals from confint() answer the same question: when an interval for an odds ratio excludes 1, the association is statistically significant at the 5% level.

Figure 4.5 shows the distribution of self-rated health among non-smokers and smokers that lies behind the odds ratio for smoking.

These odds ratios can be trusted only if the assumptions of the model are reasonable, and the next part turns to those assumptions.

Assumptions and How to Check Them

The ordinal model shares the general assumptions of logistic regression and adds one of its own. Table 4.5 lists each assumption with the way it is checked.

Table 4.5. The assumptions of ordinal logistic regression and how each is checked.

The proportional-odds assumption in the last row of Table 4.5 is the one that the ordinal model adds. The accordion below explains how the Brant test checks it and what can be done if it fails.

This section has shown how the ordinal model is fitted, read and checked. The card below links to a narrated walkthrough of ordinal logistic regression in R, and Activity 4.3 fits the model of this section to the PHAA data, reproduces its odds ratios and runs the Brant test. A knowledge check and a reflection close the section.

library(MASS)    # polr() fits the ordinal logistic model
library(brant)   # brant() checks the proportional-odds assumption
phaa <- read.csv("phaa_outcomes.csv")   # file must be in the working directory
phaa <- na.omit(phaa)                   # keep the 795 people with no missing values
phaa$smoker <- factor(phaa$smoker, levels = c("No", "Yes"))   # "No" = reference
phaa$self_rated_health <- factor(phaa$self_rated_health,
  levels = c("Poor", "Fair", "Good", "Very good", "Excellent"), ordered = TRUE)
table(phaa$smoker, phaa$self_rated_health)   # check: people in each category, by smoking

ord <- polr(self_rated_health ~ age + smoker, data = phaa, Hess = TRUE)
summary(ord)

exp(cbind(OR = coef(ord), confint(ord)))   # odds ratios with 95% CIs

brant(ord)       # check the proportional-odds assumption

Lesson 4 · Section 3

Multinomial Logistic Regression

A model for outcomes with three or more categories that have no natural order, such as usual place of care.

Nominal outcomes

Categories with no order

People name one usual place of care.

Each cancer is diagnosed as one type.

Commuters have one main way of travelling to work.

Nicotine users choose a type of product.

No category is higher or lower than another, so each category is allowed its own relationship with the predictors.

The idea

Each category is compared with a reference

Walk-in clinicversusFamily doctor

Emergency departmentversusFamily doctor

No usual placeversusFamily doctor

With K categories there are K − 1 comparisons. Each comparison has its own intercept and its own slope for every predictor.

The PHAA example

Usual place of care by smoking

Reading the output

Relative risk ratios for smoking

Comparison (versus Family doctor)	RRR for smokers	95% CI
Walk-in clinic	1.72	1.08 to 2.74
Emergency department	3.01	1.57 to 5.75
No usual place	2.66	1.54 to 4.61

The relative risk ratio (RRR) compares the chance of one category against the reference for smokers and non-smokers of the same age.

Reading the output

Relative risk ratios for age (per year)

Comparison (versus Family doctor)	RRR per year	95% CI	RRR per 10 years
Walk-in clinic	0.977	0.96 to 0.99	0.79
Emergency department	0.983	0.96 to 1.00	0.84
No usual place	0.960	0.94 to 0.98	0.66

One predictor can have a clear association in one comparison and none in another. The ten-year values come from the unrounded coefficients.

Predicted probabilities

The model in plain numbers: two 30-year-olds

	Family doctor	Walk-in clinic	Emergency department	No usual place
Non-smoker, age 30	51%	28%	6%	15%
Smoker, age 30	32%	30%	12%	26%

Each row adds to 100%, apart from rounding. predict(mn, newdata, type = "probs") produces these values.

Assumptions and checks

What to check in a multinomial model

Assumption	How to check it
Categories have no order and do not overlap	Reasoning about the categories
Observations are independent	Study design
Enough people in every category	table(outcome, predictor)
No strong multicollinearity	Correlations among the predictors
Independence of irrelevant alternatives (IIA)	Reasoning: are any two categories close substitutes?

Choosing between the two

Ordinal or multinomial?

Ordered categories

Start with ordinal logistic regression. If the Brant test fails, a multinomial model is a reasonable fallback.

Unordered categories

Use multinomial logistic regression. The ordinal model does not apply.

For five categories and two predictors, the ordinal model estimates 2 slopes and a multinomial model estimates 8.

Carry forward

What to take into the next section

• A multinomial model compares each category with a reference category, giving K − 1 comparisons.
• Each relative risk ratio belongs to one comparison, and predicted probabilities are often the clearest summary.
• The model needs enough people in every category and assumes the categories are distinct alternatives.

Introduction and Overview

This section presents multinomial logistic regression, the model for an outcome with three or more categories that have no natural order. The running example is the usual place of care, a question asked in many Canadian health surveys because having a regular source of primary care is linked to better access to preventive services and continuity of care. In the PHAA data, usual place of care is a simulated variable added for this lesson, with four categories: family doctor, walk-in clinic, emergency department and no usual place. The section fits the model with age and smoking as predictors, interprets its relative risk ratios and predicted probabilities, and lists the assumptions to check.

The section begins with what makes an outcome nominal and why such an outcome needs a model of its own.

What Makes an Outcome Nominal?

An outcome is nominal when its categories are simply different from one another, with no ranking. A family doctor, a walk-in clinic and an emergency department are different kinds of places, and none of them is "more" than the others. Type of cancer, main mode of travel to work, type of nicotine product and reason for missing a vaccination are other examples. Because there is no order to use, a model for a nominal outcome must let each category have its own relationship with every predictor.

A smoker might be more likely than a non-smoker to use an emergency department and also more likely to have no usual place, and the two increases need not be the same size. An ordinal model, which gives each predictor a single odds ratio for "moving up the scale", cannot describe this pattern, because there is no scale to move along.

The next part shows how the multinomial model gives each category its own coefficients by comparing it with a reference category.

How the Model Works: Comparisons With a Reference Category

The multinomial model chooses one category as the reference category and compares every other category with it. In R, the reference is set with relevel(), and this lesson uses family doctor, the most common answer. Figure 4.6 sets out the comparisons that follow from this choice. With four categories, the model makes three comparisons:

In general, an outcome with K categories gives K − 1 comparisons. Each comparison is written like a small logistic regression, with its own intercept and its own slope for every predictor (Hosmer, Lemeshow & Sturdivant, 2013). The model fits all of the comparisons at the same time, so that the predicted probabilities of the four categories always add to 1.

Equation 4.4 writes the model for each comparison with the reference category, with age and smoking as predictors.

The choice of reference category changes how the results are presented. It does not change how well the model fits or the predicted probabilities. A different reference (for example, walk-in clinic) would give a different set of comparisons from the same fitted model. A good reference category is a large group that is a natural point of comparison for the research question.

Each comparison has its own coefficients, and the next part shows how their exponentiated values are read.

Interpreting the Results: Relative Risk Ratios

The exponentiated coefficients of a multinomial model are called relative risk ratios (RRRs). Each RRR describes one comparison. For smoking in the emergency department comparison, the RRR says how many times larger the ratio P(emergency department) ÷ P(family doctor) is for smokers than for non-smokers of the same age. An RRR above 1 means that the predictor is associated with a higher chance of the category relative to the reference, and an RRR below 1 means that it is associated with a lower chance of the category relative to the reference.

Worked Example 4.3 calculates a crude RRR for smoking from the counts in the PHAA data and compares it with the RRR from the model.

Table 4.6 gives the RRRs from the fitted model for smoking and for age in each of the three comparisons.

Table 4.6. Relative risk ratios for usual place of care from the multinomial model, with family doctor as the reference category.

Smoking is associated with all three comparisons: among people of the same age, smokers are more likely than non-smokers to name a walk-in clinic, an emergency department or no usual place, each compared with a family doctor. Age is clearly associated with two of the comparisons. Older people are less likely to use a walk-in clinic (RRR 0.79 per ten years) or to have no usual place (RRR 0.66 per ten years), each compared with having a family doctor. The interval for the emergency department comparison includes 1, so the data do not show a clear age difference for that comparison. The ten-year values are calculated from the unrounded coefficients, for example exp(10 * coef(mn)[, "age"]).

Box 4.6 explains what an RRR describes and what it leaves out.

The next part turns to these predicted probabilities.

Predicted Probabilities

The predict() function converts the fitted model into the probability of each category for any chosen combination of predictor values. Predicted probabilities are often the clearest way to present a multinomial model to readers who are not statisticians. Table 4.7 shows them for two 30-year-olds who differ only in smoking status.

Table 4.7. Predicted probabilities of each usual place of care for a non-smoker and a smoker aged 30.

For the 30-year-old smoker, the probability of naming an emergency department (0.12) is twice that for the non-smoker (0.06), and the probability of having no usual place rises from 0.15 to 0.26. These numbers can be stated directly in a report: in this population, about one in four 30-year-old smokers has no usual place of care.

Figure 4.7 shows the observed percentages in each category by smoking status, for all ages combined, which can be set beside the predicted probabilities in Table 4.7.

The relative risk ratios and predicted probabilities can be trusted only if the assumptions of the model hold, and the next part sets them out.

Assumptions and How to Check Them

Table 4.8 lists the assumptions of the multinomial model with the way each is checked. The first four resemble the assumptions of the ordinal model in Table 4.5 in Section 2, and the last, the independence of irrelevant alternatives (IIA), is new.

Table 4.8. The assumptions of multinomial logistic regression and how each is checked.

The accordion below explains the IIA assumption with an example and sets out why the multinomial model needs more data than the ordinal model.

With both models for categorical outcomes now described, the last part of this section compares them.

Ordinal or Multinomial?

When the categories of an outcome have no order, multinomial logistic regression is the appropriate choice. When the categories are ordered, ordinal logistic regression is the place to start, because it uses the order and summarizes each predictor with a single odds ratio. If the proportional-odds assumption fails, a multinomial model is a reasonable alternative, because it makes no assumption about order. The trade-off is that the multinomial model estimates many more coefficients and its results take longer to describe.

The card below links to a narrated walkthrough of multinomial logistic regression in R. Activity 4.4 then fits the model of this section to the PHAA data and reproduces the RRRs in Table 4.6 and the predicted probabilities in Table 4.7. A knowledge check and a reflection close the section.

library(nnet)    # multinom() fits the multinomial logistic model
phaa <- read.csv("phaa_outcomes.csv")   # file must be in the working directory
phaa <- na.omit(phaa)                   # keep the 795 people with no missing values
phaa$smoker <- factor(phaa$smoker, levels = c("No", "Yes"))   # "No" = reference
phaa$usual_care <- relevel(factor(phaa$usual_care), ref = "Family doctor")   # reference category
table(phaa$usual_care, phaa$smoker)   # check: people in each category, by smoking

mn <- multinom(usual_care ~ age + smoker, data = phaa, trace = FALSE)
summary(mn)

round(exp(coef(mn)), 2)      # relative risk ratios (RRR)
round(exp(confint(mn)), 2)   # 95% confidence intervals for the RRRs

new_people <- data.frame(age = 30, smoker = c("No", "Yes"))   # two 30-year-olds
round(predict(mn, newdata = new_people, type = "probs"), 2)  # predicted probabilities

Lesson 4 · Section 4

Poisson and Negative Binomial Regression Modelling

Two models for outcomes that count how many times something happened.

Count outcomes

What counts look like

• Counts are whole numbers: 0, 1, 2, 3, and so on.
• Counts can never be negative.
• Counts are usually skewed, with many zeros and a few large values.

464

people (58%) with no urgent care visits

0.80

average number of visits per person

The Poisson distribution

One number describes it: the mean

Poisson regression

A log link and incidence rate ratios

Poisson regression

\[ \ln\big(\color{#0B7B6B}{\text{expected visits}}\big) = \color{#6D28D9}{\beta_0} + \color{#1D4ED8}{\beta_1}\,\text{age} + \color{#1D4ED8}{\beta_2}\,\text{smoker} \qquad \color{#C2410C}{\text{IRR}} = e^{\beta} \]

expected visits the average count β₀ intercept β₁, β₂ slopes on the log scale IRR incidence rate ratio

1.44

IRR for smokers (95% CI 1.19 to 1.73)

1.13

IRR for 10 years of age (1.0126 per year)

Different follow-up times

The offset turns counts into rates

Poisson regression with an offset

\[ \ln\big(\color{#0B7B6B}{\text{expected count}}\big) = \color{#C2410C}{\ln(\text{follow-up time})} + \color{#6D28D9}{\beta_0} + \color{#1D4ED8}{\beta_1}\,\text{age} + \color{#1D4ED8}{\beta_2}\,\text{smoker} \]

ln(follow-up time) the offset, coefficient fixed at 1 expected count over each person's follow-up

glm(gp_visits ~ age + smoker + offset(log(fu_years)), family = poisson)

Use an offset whenever people (or places) are observed for different amounts of time or have different population sizes.

Checking the Poisson model

Is the variance close to the mean?

Check	What a good Poisson fit looks like	PHAA urgent care visits
Mean and variance of the counts	About equal	Mean 0.80, variance 1.76
Dispersion ratio	Close to 1	2.05
Zeros: observed versus expected	About equal	464 observed, 363 expected

All three checks point to overdispersion: the counts are more spread out than the Poisson model allows.

Why it matters

Overdispersion makes the Poisson model overconfident

Negative binomial regression

One extra number lets the variance exceed the mean

Variance of a negative binomial count

\[ \text{variance} = \color{#0B7B6B}{\mu} + \frac{\color{#0B7B6B}{\mu}^2}{\color{#C2410C}{\theta}} \]

μ the expected count θ theta: small values mean extra spread

Smoking	IRR	95% CI	p-value
Poisson	1.44	1.19 to 1.73	0.0001
Negative binomial (θ = 0.77)	1.44	1.08 to 1.91	0.012

Comparing the two models

Three checks favour the negative binomial model

Check	Poisson	Negative binomial
AIC (lower is better)	2155	1954
Dispersion ratio (close to 1 is good)	2.05	0.99
Zeros expected (464 observed)	363	463

A difference in AIC of 10 or more points is usually taken as strong support for the model with the lower value.

Assumptions and checks

Count models: what to check

Assumption	Poisson	Negative binomial	How to check it
Whole-number counts of 0 or more	Yes	Yes	table(outcome)
Independent observations	Yes	Yes	Study design
Predictors multiply the expected count	Yes	Yes	Log link (built in)
Offset if follow-up time differs	Yes	Yes	offset(log(time))
Variance equals the mean	Yes	No (variance can be larger)	Dispersion ratio, AIC, zeros

Lesson summary

Outcome type, model and main check

Outcome	Model	Main check
Continuous	Linear regression	Residual plots: plot(model)
Binary	Logistic regression	About 10 events per predictor
Ordinal	Ordinal logistic regression	Brant test: brant(model)
Nominal	Multinomial logistic regression	People in every category; distinct categories
Count	Poisson regression	Dispersion ratio close to 1
Overdispersed count	Negative binomial regression	Lower AIC; expected zeros

Introduction and Overview

This section presents the two standard models for counts, Poisson regression and negative binomial regression. The running example is the number of urgent care visits each PHAA participant made in the past 12 months, a simulated variable added for this lesson. The section describes the features of count data and the Poisson distribution, fits a Poisson regression and interprets its incidence rate ratios, checks the model for overdispersion, and then fits a negative binomial regression and compares the two models.

The section begins with the features that set counts apart from the other types of outcome in this lesson.

Count Outcomes

A count records how many times an event happened to a person, or in a place, over a period of time. Counts are whole numbers, they cannot be negative, and their distribution is usually skewed to the right, with many zeros and small values and a long tail of larger values. In the PHAA data, the 795 participants made between 0 and 12 urgent care visits in the past year: 464 (58%) made none, 174 made one and 84 made two, and only 4 made eight or more. The average was 0.80 visits per person.

A linear regression of these counts would break several assumptions at once. Its residuals would be skewed (not normal), their spread would grow with the predicted value (unequal variance), and the model could predict a negative number of visits. Figure 4.2 in Section 1, which showed a linear model fitted to urgent_visits, illustrates the first two problems.

Counts therefore need a distribution suited to whole numbers of zero or more, and the next part introduces the one on which Poisson regression is built.

The Poisson Distribution

The Poisson distribution describes the number of events in a fixed period when events happen independently of one another at a constant average rate. It is named after the French mathematician Siméon-Denis Poisson. The whole distribution is defined by one number, its mean (written μ, the Greek letter mu), and its most important property is that its variance equals its mean. A Poisson count with a mean of 3 has a variance of 3, and a Poisson count with a mean of 0.8 has a variance of 0.8.

Figure 4.8 shows how the shape of the distribution changes as its mean increases.

The mean-equals-variance property is what makes the Poisson distribution simple, and it is also what most often goes wrong when it is applied to real health data. The checks later in this section are built around it.

Poisson regression, the subject of the next part, uses this distribution to model how the expected count changes with the predictors.

Poisson Regression

Poisson regression is the GLM with a Poisson distribution and a log link. The log of the expected count is set equal to the linear predictor.

Equation 4.5 gives the model for urgent care visits with age and smoking as predictors.

Because the predictors act by multiplication, the exponentiated coefficient e^β is a ratio of expected counts, called the incidence rate ratio (IRR). An IRR of 1 means no association, an IRR above 1 means more events, and an IRR below 1 means fewer events. In the PHAA data, the Poisson model gives an IRR for smoking of 1.44 (95% CI 1.19 to 1.73): smokers make about 44% more urgent care visits than non-smokers of the same age. The IRR for age is 1.013 per year (95% CI 1.007 to 1.018). Calculated from the unrounded coefficient, exp(10 × 0.01252) = 1.13 per ten years, or about 13% more visits for each decade of age.

Worked Example 4.4 calculates a crude rate ratio for smoking from the average numbers of visits in the two groups and compares it with the IRR from the model.

The next part shows how the same model is adapted when the length of follow-up differs from person to person.

Different Follow-Up Times: The Offset

Every PHAA participant reported urgent care visits over the same 12 months, so their counts can be compared directly. Many studies instead follow people for different lengths of time, or count events in places with different population sizes. A person followed for ten years has more time to accumulate visits than a person followed for one year, and a city of a million people has more cases than a town of ten thousand. Poisson regression handles this with an offset: the log of each unit's follow-up time (or population size) is added to the model with its coefficient fixed at 1. The offset turns the model for counts into a model for rates, such as visits per person-year, which is why the exponentiated coefficients are called incidence rate ratios (Frome, 1983).

Worked code 4.1 fits a Poisson model with an offset to the PHAA follow-up file, in which the length of follow-up differs from person to person, and shows how its IRRs are then read as ratios of rates.

fu <- read.csv("phaa_followup.csv")      # GP visits over 0 to 10 years of follow-up
fu$smoker <- factor(fu$smoker, levels = c("No", "Yes"))
gp <- glm(gp_visits ~ age + smoker + offset(log(fu_years)), data = fu, family = poisson)
round(exp(cbind(IRR = coef(gp), confint(gp))), 3)

An offset accounts for differences in follow-up time. The next part turns to the main check of a Poisson model, which concerns its assumption that the variance equals the mean.

Checking the Poisson Model: Overdispersion

The Poisson model assumes that the variance of the counts equals their mean, after the predictors are taken into account. When the variance is larger than the mean, the data are overdispersed. Overdispersion is common in health data, because people differ in ways the predictors do not capture: some people use health care much more than others of the same age and smoking status. Three checks reveal it.

Table 4.9 sets out the three checks, with the result of each for the PHAA data.

Table 4.9. Three checks for overdispersion in a Poisson model, with the results for urgent care visits in the PHAA data.

All three checks point to overdispersion. The consequence is that the Poisson model is overconfident. Because it assumes a smaller variance than the data actually have, its standard errors are too small, its confidence intervals are too narrow and its p-values are too small. An analysis that ignored the problem could report an association as more certain than the data justify.

Figure 4.9 compares the observed numbers of people with each count of visits with the numbers expected by the Poisson model and by the negative binomial model introduced in the next part.

Box 4.7 describes a source of apparent overdispersion that should be ruled out before the Poisson model is replaced.

When the model includes the predictors and offset that it needs and the checks still show overdispersion, a model that allows the variance to exceed the mean is required. The next part introduces that model.

Negative Binomial Regression

Negative binomial regression keeps the log link and the IRR interpretation of Poisson regression and adds one extra parameter, θ (theta), that allows the variance to exceed the mean. The variance of a negative binomial count is μ + μ²/θ. A small θ means a lot of extra variation, and as θ becomes very large the second term shrinks toward zero and the model becomes a Poisson model (Hilbe, 2011). In R, the model is fitted with glm.nb() from the MASS package, and the output reports θ as Theta.

For the urgent care visits, θ is 0.77. The negative binomial model gives an IRR for smoking of 1.44 (95% CI 1.08 to 1.91), almost the same point estimate as the Poisson model. The confidence interval, however, is wider (it was 1.19 to 1.73 under the Poisson model), and the p-value rises from 0.0001 to 0.012. The IRR for age is 1.012 per year (95% CI 1.004 to 1.020). The associations remain, and the negative binomial model reports them with appropriate uncertainty.

Table 4.10 sets the results of the two models side by side.

Table 4.10. Results of the Poisson and negative binomial models for urgent care visits.

The comparison favours the negative binomial model on every check. The AIC (Akaike information criterion) measures how well a model fits the data, with a penalty for each parameter it estimates, and a lower value is better. A difference of about 10 points or more is usually taken as strong support for the model with the lower AIC (Burnham & Anderson, 2004), and here the difference is about 200 points. The dispersion ratio of the negative binomial model is 0.99, and it expects 463 zeros against the 464 observed.

The accordion below describes the models used when the counts contain many more zeros than even a negative binomial model expects.

The two count models share most of their assumptions, and the next part sets them out together.

Assumptions and How to Check Them

Table 4.11 lists the assumptions of the Poisson and negative binomial models, whether each model requires them and how each is checked. The two models differ only in the assumption about the variance.

Table 4.11. The assumptions of the Poisson and negative binomial models and how each is checked.

With the count models complete, the last part of this section draws the four sections of the lesson together.

Bringing the Lesson Together

The four sections of this lesson apply the same approach to five types of outcome. The type of the outcome is identified first, and it points to a model. The model is fitted with the same linear predictor (here, age and smoking), its exponentiated coefficients are read as the ratio that its link function implies, and its own assumptions are checked before the results are reported.

Table 4.12 summarises the five PHAA examples, with the model, the effect of smoking and the main check for each.

Table 4.12. The five types of outcome in the PHAA examples, with the model and R function, the effect of smoking (95% CI) and the main check for each.

The card below links to a narrated walkthrough of count models in R. Activity 4.5 fits the Poisson model for urgent care visits and runs the overdispersion checks in Table 4.9, and Activity 4.6 fits the negative binomial model and reproduces the comparison in Table 4.10. A knowledge check and a reflection close the section, and the final page of the lesson gathers its key takeaways.

library(MASS)    # glm.nb() fits the negative binomial model (next activity)
phaa <- read.csv("phaa_outcomes.csv")   # file must be in the working directory
phaa <- na.omit(phaa)                   # keep the 795 people with no missing values
phaa$smoker <- factor(phaa$smoker, levels = c("No", "Yes"))   # "No" = reference
table(phaa$urgent_visits)    # how many people had 0, 1, 2, ... visits
mean(phaa$urgent_visits)     # the average count
var(phaa$urgent_visits)      # the variance (spread) of the counts

pois <- glm(urgent_visits ~ age + smoker, data = phaa, family = poisson)
summary(pois)
exp(cbind(IRR = coef(pois), confint(pois)))   # incidence rate ratios with 95% CIs

# Check 1: the dispersion ratio is close to 1 when the Poisson model fits
sum(residuals(pois, type = "pearson")^2) / df.residual(pois)
# Check 2: compare the zeros we observed with the zeros the model expects
sum(phaa$urgent_visits == 0)
sum(dpois(0, fitted(pois)))

nb <- glm.nb(urgent_visits ~ age + smoker, data = phaa)
summary(nb)
exp(cbind(IRR = coef(nb), confint(nb)))       # incidence rate ratios with 95% CIs

AIC(pois, nb)                                               # lower AIC = better fit
sum(residuals(nb, type = "pearson")^2) / df.residual(nb)    # dispersion ratio for NB
sum(dnbinom(0, mu = fitted(nb), size = nb$theta))           # zeros expected by NB

Bringing It All Together

This lesson presented one regression model for each common type of outcome. Section 1 sorted variables into numeric (continuous or discrete) and categorical (binary, ordinal or nominal) types, reviewed the four assumptions of linear regression (linearity, independence, normality and equal variance of the residuals) and the plots that check them, and introduced the generalized linear model, which keeps the linear predictor and changes the distribution and the link function. It also reviewed logistic regression for binary outcomes and its requirement of enough events per predictor.

Section 2 fitted ordinal logistic regression to self-rated health. The model fits every at-or-above split of the outcome at once and gives each predictor one odds ratio for being in a higher category, an assumption that the Brant test checks. Section 3 fitted multinomial logistic regression to usual place of care. The model compares each category with a reference category, reports relative risk ratios for each comparison, and is often clearest when presented as predicted probabilities. Section 4 fitted Poisson regression to urgent care visits, found that the counts were overdispersed (a dispersion ratio of 2.05 and too few expected zeros), and fitted a negative binomial model that matched the data and gave wider confidence intervals that reflect the extra variation in the counts.

Across the lesson, the same two predictors were used with five different outcomes, and the same steps were followed each time: the outcome type was identified, the matching model was fitted, the exponentiated coefficients were read as the ratio that the link function implies, and the model's own assumptions were checked before the results were reported. The final assessment asks for these steps to be applied to new outcomes. The next lesson, Modelling Dependent Data, takes up the one assumption that every model in this lesson shares and that none of them can check from its own output: the independence of the observations.

The final assessment covers all four sections. All 15 questions must be answered correctly (100%), and the final reflection completed, to finish the lesson.

Final Knowledge Assessment