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 3 · Section 1

Comparing Two Variables and the Case for Regression

Simple tests compare two variables at a time, and regression compares many at once.

Running example

The PHAA survey and one question

800

adults in the PHAA survey

108

mean systolic BP (mmHg), SD 11

1

question: what goes with higher blood pressure?

The lesson asks which characteristics are associated with higher blood pressure, and by how much.

Two-variable tests

Parametric tests and their rank-based alternatives

Situation	Parametric test	Rank-based or exact alternative
Two numeric variables	Pearson correlation	Spearman’s ρ
Numeric outcome, two independent groups	Two-sample t-test	Mann–Whitney U (Wilcoxon rank-sum)
Numeric outcome, paired measurements	Paired t-test	Wilcoxon signed-rank
Numeric outcome, three or more groups	One-way ANOVA	Kruskal–Wallis
Two categorical variables	Pearson χ² test	Fisher’s exact test (small expected counts)

Rank-based tests suit skewed or ordinal outcomes; the reading lists the R function for every test.

Two-variable tests

Tests for paired, ordered and time-to-event data

Test	Used for
McNemar’s test	Paired yes-or-no outcomes, such as before-and-after measurements or matched pairs
Cochran–Armitage trend test	A steady rise or fall in a yes-or-no outcome across ordered categories
Log-rank test	Survival curves, or time to an event, in two or more groups

Every test in this section has a regression counterpart, listed in the reading.

Correlation

Age and blood pressure move together

t-test and ANOVA

Comparing group means

Smoking (t-test)

Smokers average 113.5 mmHg and non-smokers 107.2 mmHg, a difference of 6.3 mmHg (95% CI 4.2 to 8.4, p < 0.001).

Region (ANOVA)

The regional averages range from 106.1 to 109.1 mmHg, and the ANOVA p-value of 0.14 is consistent with chance.

Chi-square test

Elevated blood pressure by smoking

	Elevated BP: no	Elevated BP: yes	Share with elevated BP
Non-smokers	583	80	12.1%
Smokers	87	48	35.6%

Elevated blood pressure is defined here as a systolic reading of 120 mmHg or more (χ² = 44.2, p < 0.001).

The limit of two-variable tests

A crude association can be mostly confounding

The idea

Regression fits a line through the data

Simple linear regression

\[ \color{#0B7B6B}{Y} = \color{#6D28D9}{\beta_0} + \color{#C2410C}{\beta_1}\color{#1D4ED8}{X} + \varepsilon \]

Y outcome β₀ intercept β₁ slope per unit of X X predictor

Least squares chooses the line with the smallest sum of squared residuals.

One framework

The t-test is a regression with one binary predictor

Analysis	Result
t.test(systolic_bp ~ smoker, var.equal = TRUE)	Difference 6.30 mmHg, t = 6.01, p = 2.9 × 10−9
lm(systolic_bp ~ smoker)	Slope 6.30 mmHg, t = 6.01, p = 2.9 × 10−9

A regression with region as a predictor reproduces the ANOVA in the same way.

Multivariable regression

Each coefficient is adjusted for the others

Multivariable linear regression

\[ \begin{aligned} \color{#0B7B6B}{Y} = {} & \color{#6D28D9}{\beta_0} + \color{#C2410C}{\beta_1}X_1 \\ & + \color{#C2410C}{\beta_2}X_2 + \dots + \varepsilon \end{aligned} \]

Adjusted estimate

Each slope compares people who have the same values of the other predictors in the model.

Its limit

Regression adjusts only for confounders that were measured and included in the model.

Carry forward

What to take into the next section

• The types of the two variables decide which two-variable test is used, and rank-based tests suit skewed or ordinal data.
• A two-variable comparison can be distorted by confounding.
• Regression contains these tests and adjusts each coefficient for the other predictors.

Introduction and Overview

This section begins with the familiar tests that compare two variables at a time and then introduces regression as the method that compares many variables at once. All of the examples use phaa_survey_clean.csv, the cleaned PHAA survey of 800 adults from Lesson 2, with systolic blood pressure as the main outcome.

Choosing a Two-Variable Test

A two-variable test asks whether two variables are associated. The choice of test depends on the type of each variable, which Lesson 2 described as numeric (continuous or discrete) or categorical (binary, nominal or ordinal). The table below lists the four tests used most often in this course.

Each test has its own assumptions. The correlation coefficient measures a straight-line relationship. The t-test and ANOVA assume roughly normal outcomes within each group (or reasonably large groups), and the classic versions assume similar spread in each group; R's default t.test() uses the Welch version, which does not require equal spread, and var.equal = TRUE requests the classic version used in this section. The chi-square test needs an expected count of about 5 or more in each cell of the table. When these assumptions are doubtful, a rank-based or exact version of the test is used instead, as the box below the next example shows.

What Two-Variable Tests Cannot Do

A two-variable test compares people who differ in the predictor, but those people may also differ in many other ways. When a third variable is related to both the predictor and the outcome, it can create, hide or exaggerate an association. Such a variable is called a confounder.

Among the 768 participants with blood pressure, age and physical activity recorded, each additional 100 minutes of activity is associated with a systolic blood pressure 2.81 mmHg lower (95% CI 1.86 to 3.76 mmHg lower). Activity and age are negatively correlated (r = −0.27), and blood pressure rises by about 0.43 mmHg per year of age. When age is added to the model, the difference per 100 minutes falls to 0.91 mmHg (95% CI 0.06 to 1.76 mmHg lower, p = 0.037), a reduction of about two thirds. Most of the crude association reflects the fact that the more active participants are younger.

Regression: One Framework for These Tests

Linear regression describes how the average of a continuous outcome changes with one or more predictors. With one predictor, it fits a straight line.

The fitted line for age is SBP = 87.9 + 0.45 × age. The intercept, 87.9 mmHg, is the predicted value at age 0, which is outside the data and has no practical meaning here; it fixes the height of the line. The slope, 0.45 mmHg per year, is the quantity of interest. The residual for each person is the observed blood pressure minus the value on the line, and ordinary least squares (OLS) chooses the intercept and slope that make the sum of the squared residuals as small as possible.

The Simple Tests Are Special Cases of Regression

A regression with a single binary predictor reproduces the t-test. In lm(systolic_bp ~ smoker), the intercept (107.16 mmHg) is the mean for non-smokers, the reference group, and the slope for smokerYes (6.30 mmHg) is the difference between the two means. Its t statistic (6.01) and p-value (2.9 × 10⁻⁹) match those of t.test(..., var.equal = TRUE) exactly. A regression with a categorical predictor of three or more groups reproduces the one-way ANOVA, and the F-test printed at the bottom of its summary() is the ANOVA F-test. With one numeric predictor, the slope equals the correlation multiplied by the ratio of the standard deviations of the outcome and the predictor, and the R² of the regression equals the squared correlation (0.54² = 0.29 for age).

Regression with More Than One Predictor

A multivariable regression includes several predictors. Each coefficient then describes the change in the average outcome for a one-unit change in its predictor among people with the same values of the other predictors. In the physical activity example, the coefficient of −0.0091 mmHg per minute in the model with age compares people of the same age. This property, adjustment for the other variables in the model, is the main reason regression is used in epidemiology.

phaa <- read.csv("phaa_survey_clean.csv")   # file must be in the working directory
phaa$smoker <- factor(phaa$smoker, levels = c("No", "Yes"))              # "No" = reference
phaa$gender <- factor(phaa$gender, levels = c("Man", "Woman", "Non-binary"))   # "Man" = reference
cor.test(phaa$age, phaa$systolic_bp)          # correlation: two numeric variables

t.test(systolic_bp ~ smoker, data = phaa, var.equal = TRUE)   # t-test: two group means

round(tapply(phaa$systolic_bp, phaa$region, mean, na.rm = TRUE), 1)   # mean by region
summary(aov(systolic_bp ~ region, data = phaa))                       # ANOVA: three or more means

phaa$elevated_bp <- factor(ifelse(phaa$systolic_bp >= 120, "Yes", "No"),
                           levels = c("No", "Yes"))   # "Yes" (120 mmHg or more) is the event
tab <- table(phaa$smoker, phaa$elevated_bp)
tab                                            # counts
round(prop.table(tab, margin = 1), 3)          # share with elevated BP in each smoking group
chisq.test(tab)                                # chi-square test: two categorical variables

phaa <- read.csv("phaa_survey_clean.csv")   # file must be in the working directory
phaa$smoker <- factor(phaa$smoker, levels = c("No", "Yes"))              # "No" = reference
phaa$gender <- factor(phaa$gender, levels = c("Man", "Woman", "Non-binary"))   # "Man" = reference
fit_smk <- lm(systolic_bp ~ smoker, data = phaa)   # regression with one binary predictor
summary(fit_smk)$coefficients

act <- na.omit(phaa[, c("systolic_bp", "phys_act_min", "age")])   # same people in both models
crude <- lm(systolic_bp ~ phys_act_min, data = act)            # physical activity alone
adjusted <- lm(systolic_bp ~ phys_act_min + age, data = act)   # adding age
round(summary(crude)$coefficients, 4)
round(summary(adjusted)$coefficients, 4)

100 * coef(crude)[["phys_act_min"]]      # difference per 100 minutes, crude
100 * coef(adjusted)[["phys_act_min"]]   # difference per 100 minutes, adjusted for age

Lesson 3 · Section 2

Linear Regression and Model Checking

Fitting a linear model, reading its output, and checking that its assumptions hold.

The model

Five predictors of systolic blood pressure

lm(systolic_bp ~ age + gender + smoker + bmi + dep_score, data = cc)

Who is included

The model uses the 795 people with no missing values for the outcome and the five predictors.

Categorical predictors

Gender enters as two indicator variables (woman and non-binary), each compared with men.

Reading the coefficients

Adjusted differences in systolic blood pressure

Predictor	Estimate (mmHg)	95% CI	p-value
Age (per year)	0.43	0.38 to 0.47	< 0.001
Woman vs man	−0.34	−1.58 to 0.90	0.59
Non-binary vs man	−0.61	−4.03 to 2.81	0.73
Smoker vs non-smoker	6.15	4.49 to 7.80	< 0.001
BMI (per unit)	0.43	0.23 to 0.63	< 0.001
Depression score (per point)	0.38	0.30 to 0.46	< 0.001

t-tests and confidence intervals

Is each coefficient different from zero?

t statistic and 95% confidence interval

\[ t = \frac{\color{#C2410C}{\hat\beta}}{\color{#1D4ED8}{SE(\hat\beta)}} \] \[ \color{#C2410C}{\hat\beta} \pm 1.96 \times \color{#1D4ED8}{SE(\hat\beta)} \]

β̂ estimated coefficient SE its standard error

For smoking, t = 6.15 ÷ 0.84 = 7.3 and the 95% CI is 4.5 to 7.8 mmHg.

Goodness of fit

How much does the model explain?

0.41

R² (adjusted R² also 0.41)

8.7

residual standard error (mmHg)

92.7

F statistic on 6 and 788 df, p < 0.001

A low R² is common in health research and does not by itself mean that a model is wrong.

Assumptions

The LINE assumptions and two further checks

Assumption	Main check
Linearity	Residuals vs fitted plot; a squared term for age
Independence	The study design (one row per person)
Normal residuals	Q-Q plot; Shapiro-Wilk test
Equal variance	Scale-location plot; Breusch-Pagan test
No undue influence	Residuals vs leverage plot; Cook's distance
No strong collinearity	Variance inflation factors (VIF)

Diagnostic plots

Four plots, all adequate

Formal tests

Tests support the plots

Check	R code	Result
Linearity in age	anova(model, model_sq)	F = 0.43, p = 0.51
Equal variance	bptest(model)	BP = 4.68, p = 0.58
Normal residuals	shapiro.test(residuals(model))	W = 0.998, p = 0.62

The test asks whether there is a departure, and the plot shows how large it is.

Influence and collinearity

No influential points and no collinearity

Cook's distance

The largest value is 0.061, far below 0.5, although 35 rows pass the liberal 4/n screening line.

Variance inflation factors

All values are below 1.12, far below the warning thresholds of 5 and 10.

When a check fails

Common problems and first responses

What the check shows	First response
A curve in residuals vs fitted	Add a squared term (or a spline) for the predictor
A funnel of growing spread	Log-transform a positive, skewed outcome, or use sandwich (HC3) standard errors
Skewed residuals in a small sample	Transform the outcome
A highly influential point	Check for an error; report the model with and without it
Clustered or repeated observations	Use a mixed model or GEE (Lesson 5)

Transformations

A log transformation can remove a funnel

Carry forward

What to take into the next section

• Each coefficient is an adjusted difference in the average outcome, reported with its confidence interval.
• R², the residual standard error and the F-test describe how well the model fits.
• The diagnostic plots come first, and each failed check has a standard response.

Introduction and Overview

This section fits a multivariable linear regression for systolic blood pressure, reads its output, checks its assumptions, and describes what to do when an assumption fails. The model has five predictors (age, gender, smoking, BMI and the depression score) and uses the 795 PHAA participants with complete data for these variables.

Fitting and Reading the Model

The model is written in R as lm(systolic_bp ~ age + gender + smoker + bmi + dep_score, data = cc), where cc holds the complete cases. The table interprets each coefficient. Every estimate is adjusted for the other predictors in the model.

The smoking estimate in this model (6.15 mmHg) is close to the unadjusted difference of 6.30 mmHg from Section 1, but Section 4 shows that this similarity hides two changes in opposite directions: adjustment for the confounders lowers the estimate, and adjustment for BMI, a mediator, raises it.

Testing the Coefficients and the Whole Model

Each row of the summary() output gives a t-test of whether that coefficient differs from zero, computed as the estimate divided by its standard error. The 95% confidence interval, from confint(), is approximately the estimate plus or minus 1.96 standard errors. The bottom of the output describes the whole model: R² is the share of the variation in the outcome explained by the predictors (0.414), the adjusted R² applies a penalty for the number of predictors (0.409), the residual standard error is the typical size of a residual (8.70 mmHg), and the F-test asks whether the predictors together explain more variation than chance (F = 92.7 on 6 and 788 degrees of freedom, p < 0.001).

Assumptions of Linear Regression

The t-tests, confidence intervals and F-test above are trustworthy only when the model's assumptions hold. The four main assumptions are summarised by the word LINE, and two further checks are usually added.

Reading the Four Diagnostic Panels

Calling plot() on a fitted lm object in R produces four panels, each of which addresses one or more of the assumptions above. The table describes what each panel shows and how a problem appears in it.

The panels are read together. A funnel in the first panel that reappears as a rising line in the third is strong evidence of non-constant variance, whereas a single stray point that appears in the second panel and again in the fourth points to influence, and the distributional assumption may still hold for the remaining observations. R labels the three most extreme observations in each panel by row name, which makes it straightforward to inspect those rows directly.

What Adequate and Problem Panels Look Like

The four panels share a vocabulary. A fitted value is the model's prediction for an observation, and a residual is the observed outcome minus the fitted value. A standardised residual is a residual divided by its estimated standard deviation, so that values beyond about ±2 are unusual and values beyond ±3 are rare. The red line in three of the panels is a smoother that follows the average pattern of the points, and the numbers printed beside a few points are row names, which identify the observations R considers most extreme in that panel.

The first figure shows the four panels R produces for the blood pressure model in this section. All four look as an adequate model should. The figures that follow place a simulated adequate panel beside simulated problem panels so that each pattern can be recognised; the R box at the end of this subsection contains the R code that generates the simulated data, so the pictures can be reproduced.

Residuals vs Fitted

Each point is one observation, placed at its predicted value on the horizontal axis and its residual on the vertical axis. In an adequate model (panel A) the red smoother stays close to the dashed zero line and the vertical spread is similar everywhere. A curved smoother (panel B) means the predictions are systematically too low at both ends and too high in the middle, which is the signature of a non-linear relationship; the response is to model the curve, as described under remedies below. A band that widens from left to right (panel C) means the errors are larger for some fitted values than for others, which is non-constant variance.

Q-Q Residuals

R sorts the standardised residuals from smallest to largest and plots each one against the value expected for its rank if the errors were exactly normal; these expected values are the theoretical quantiles. Normal errors fall along the dashed diagonal. Small wiggles and a few points drifting at the extreme ends occur in any real dataset, as in panel A. A systematic bend at one end (panel B) indicates skewed errors, here a long tail of large positive residuals. An S-shape with both ends pulling away from the line (panel C) indicates heavy tails, meaning more extreme residuals in both directions than a normal distribution would produce.

Scale-Location

The vertical axis is the square root of the absolute standardised residual, which measures the size of each residual regardless of its sign. A flat red line (panel A) means residuals are of similar size across the range of fitted values. A rising red line (panel B) means residuals grow as the predictions grow, the signature of non-constant variance. This panel often shows the problem more clearly than the Residuals vs Fitted panel, because positive and negative residuals are combined.

Residuals vs Leverage

The horizontal axis is leverage, which measures how unusual an observation's combination of predictor values is, and the vertical axis is the standardised residual, which measures how poorly the model predicts it. The dashed grey curves are contours of Cook's distance at 0.5 and 1. Points near the right-hand corners combine high leverage with a large residual and therefore have the greatest influence on the coefficients. In panel A every point lies well inside the contours. In panel B, row 41 has a leverage of 0.25 and a standardised residual of −3.5, which places it beyond the contour of 1 (its Cook's distance is 1.39). Removing that single row moves the slope for x1 from 0.36 to 0.46.

suppressPackageStartupMessages({library(car); library(lmtest); library(sandwich)})
set.seed(410)
n  <- 300
x1 <- runif(n, 0, 10)
x2 <- rnorm(n)
ok     <- data.frame(x1, x2, y = 2 + 0.5 * x1 + 0.3 * x2 + rnorm(n))
curved <- data.frame(x1, x2, y = 2 + 0.5 * x1 + 0.12 * (x1 - 5)^2 + 0.3 * x2 + rnorm(n))
fan    <- data.frame(x1, x2, y = 2 + 0.5 * x1 + 0.3 * x2 + rnorm(n) * (0.2 + 0.35 * x1))
skew   <- data.frame(x1, x2, y = 2 + 0.5 * x1 + 0.3 * x2 + (rexp(n, 1) - 1) * 1.2)
heavy  <- data.frame(x1, x2, y = 2 + 0.5 * x1 + 0.3 * x2 + rt(n, df = 3) * 0.8)
model_ok     <- lm(y ~ x1 + x2, data = ok)
model_curved <- lm(y ~ x1 + x2, data = curved)
model_fan    <- lm(y ~ x1 + x2, data = fan)
model_skew   <- lm(y ~ x1 + x2, data = skew)
model_heavy  <- lm(y ~ x1 + x2, data = heavy)
set.seed(7)
ni  <- 40
xi  <- runif(ni, 0, 10)
xx2 <- rnorm(ni)
e   <- rnorm(ni)
inf <- data.frame(x1 = c(xi, 16), x2 = c(xx2, 0))
inf$y <- c(2 + 0.5 * xi + 0.3 * xx2 + e, 5)
model_inf <- lm(y ~ x1 + x2, data = inf)

par(mfrow = c(1, 3))
plot(model_ok, which = 1); plot(model_curved, which = 1); plot(model_fan, which = 1)
plot(model_ok, which = 2); plot(model_skew, which = 2); plot(model_heavy, which = 2)
par(mfrow = c(1, 2))
plot(model_ok, which = 3); plot(model_fan, which = 3)
plot(lm(y ~ x1 + x2, data = inf[1:40, ]), which = 5); plot(model_inf, which = 5)
par(mfrow = c(1, 1))

Formal Checks

The formal tests for this model all agree with the plots: the squared age term does not improve the fit (F = 0.43, p = 0.51), the Breusch-Pagan test finds no evidence of unequal variance (BP = 4.68 on 6 df, p = 0.58), and the Shapiro-Wilk test finds no evidence of non-normal residuals (W = 0.998, p = 0.62). The accordion below explains each test and its limits.

When an Assumption Fails

A failed check calls for a change to the analysis that matches the cause of the problem. The blood pressure model passed every check, so the examples below use simulated data to show what the problems and their remedies look like.

Choosing a Remedy in Practice

A failed check is a signal to change the analysis in a way that matches the cause of the problem. The table links the pattern seen in the diagnostics to its usual cause, a first response, and the R code that carries the response out.

Worked Example: A Log Transformation for a Skewed, Funnel-Shaped Outcome

Costs, lengths of stay and many biomarker concentrations are positive, right-skewed, and more variable when they are larger. Modelling the logarithm of such an outcome often removes the funnel. The simulated example below behaves in this way.

R Worked code: the outcome as measured and on the log scale

library(car)   # for ncvTest()
set.seed(411)
x    <- runif(300, 0, 10)
cost <- data.frame(x, y = exp(1 + 0.15 * x + rnorm(300, sd = 0.45)))
model_cost    <- lm(y ~ x, data = cost)        # outcome as measured
model_logcost <- lm(log(y) ~ x, data = cost)   # log-transformed outcome
ncvTest(model_cost)
ncvTest(model_logcost)
coef(model_logcost)
exp(coef(model_logcost)["x"])   # multiplicative change per unit of x

Console output

Non-constant Variance Score Test Variance formula: ~ fitted.values Chisquare = 58.8403, Df = 1, p = 1.71e-14 Non-constant Variance Score Test Variance formula: ~ fitted.values Chisquare = 0.3549518, Df = 1, p = 0.55132 (Intercept) x 0.9638985 0.1585067 x 1.17176

The variance test moves from p = 1.7 × 10⁻¹⁴ to p = 0.55 after the transformation. The coefficient on the log scale, 0.159, is interpreted after exponentiation, the step that reverses the logarithm: exp(0.159) = 1.17, so each one-unit increase in x multiplies the geometric mean of the outcome by 1.17, a 17% increase. The general rule is that a coefficient β for a log-transformed outcome corresponds to a percentage change of 100 × (e^β − 1), where e is the mathematical constant, about 2.718, on which natural logarithms are based. Predictions from such a model, transformed back with exp(), estimate the geometric mean of the outcome, which for a skewed outcome lies close to the median and below the arithmetic mean. A transformation requires an outcome that is always positive; for outcomes with zeros or negative values, sandwich standard errors are the simpler choice.

phaa <- read.csv("phaa_survey_clean.csv")   # file must be in the working directory
phaa$smoker <- factor(phaa$smoker, levels = c("No", "Yes"))              # "No" = reference
phaa$gender <- factor(phaa$gender, levels = c("Man", "Woman", "Non-binary"))   # "Man" = reference
cc <- na.omit(phaa[, c("systolic_bp", "age", "gender", "smoker", "bmi", "dep_score")])
nrow(cc)                                       # people with complete data for the model
model <- lm(systolic_bp ~ age + gender + smoker + bmi + dep_score, data = cc)
summary(model)

round(confint(model), 2)                     # 95% confidence intervals

# install.packages(c("lmtest", "car", "sandwich"))   # run once, if not yet installed
library(lmtest); library(car); library(sandwich)
phaa <- read.csv("phaa_survey_clean.csv")   # file must be in the working directory
phaa$smoker <- factor(phaa$smoker, levels = c("No", "Yes"))              # "No" = reference
phaa$gender <- factor(phaa$gender, levels = c("Man", "Woman", "Non-binary"))   # "Man" = reference
cc <- na.omit(phaa[, c("systolic_bp", "age", "gender", "smoker", "bmi", "dep_score")])
nrow(cc)                                       # people with complete data for the model
model <- lm(systolic_bp ~ age + gender + smoker + bmi + dep_score, data = cc)
par(mfrow = c(2, 2)); plot(model); par(mfrow = c(1, 1))   # the four diagnostic plots

model_sq <- lm(systolic_bp ~ age + I(age^2) + gender + smoker + bmi + dep_score, data = cc)
anova(model, model_sq)                         # linearity: does a squared age term help?
bptest(model)                                  # equal spread (Breusch-Pagan test)
shapiro.test(residuals(model))                 # normality of the residuals

cooks <- cooks.distance(model)
sum(cooks > 4 / nrow(cc))                      # rows above the 4/n screening line
round(max(cooks), 3)                           # the largest Cook's distance
vif(model)                                     # collinearity (values near 1 are good)

# A check that does not rely on equal spread: sandwich (HC3) standard errors
round(coeftest(model, vcov = vcovHC(model, type = "HC3"))[, 1:2], 3)

Lesson 3 · Section 3

Logistic Regression and Model Checking

A model for yes-or-no outcomes, read as odds ratios and checked for fit and discrimination.

The outcome

Elevated blood pressure: yes or no

128

of 798 people with SBP of 120 mmHg or more

16%

have elevated blood pressure

23

have the hypertension variable, too few to model

Why not linear regression?

A straight line breaks the rules for probabilities

• A straight line can predict probabilities below 0 or above 1.
• The residuals of a 0 or 1 outcome cannot be normally distributed.
• The spread of the residuals depends on the predicted probability.

Probability, odds and log odds

From a probability to an odds ratio

	Probability	Odds = p ÷ (1 − p)
Smokers	0.356	0.55
Non-smokers	0.121	0.14

The crude odds ratio is 0.55 ÷ 0.14 ≈ 4.0, and the log odds, ln[p ÷ (1 − p)], is the scale of the model.

The logistic model

A straight line on the log-odds scale is an S-curve for probability

Fitting the model

From a log-odds coefficient to an odds ratio

glm(elevated_bp ~ age + gender + smoker + bmi + dep_score, family = binomial, data = cc)

Odds ratio from a coefficient

\[ \text{OR} = e^{\color{#C2410C}{\beta}} \] \[ e^{1.754} = 5.78 \]

The coefficient for smoking is 1.754 on the log-odds scale, which is an odds ratio of 5.78.

Odds ratios

Smoking and ten years of age have the largest odds ratios

Predicted probabilities

Two men aged 50, differing only in smoking

Non-smoker

The predicted probability of elevated blood pressure is 9.7%.

Smoker

The predicted probability is 38.4%, about 3.9 times as high, while the odds ratio is 5.8.

Predictions are made with predict(model, newdata, type = "response").

Model checking 1

Events, linearity, collinearity and influence

Check	Result
Events per coefficient (about 10 or more)	126 events ÷ 6 coefficients ≈ 21
Linearity on the log-odds scale	Squared age term: likelihood ratio test p = 0.43
Collinearity	All VIFs below 1.14
Influence	Largest Cook's distance 0.065

Model checking 2

Overall fit, calibration and discrimination

Carry forward

What to take into the next section

• Logistic regression models the log odds, and exp(β) is an adjusted odds ratio.
• Predicted probabilities help to explain what an odds ratio means.
• The checks cover events per coefficient, linearity, calibration and discrimination.

Introduction and Overview

This section models a binary outcome, elevated blood pressure (a systolic reading of 120 mmHg or more), with logistic regression. It uses the same five predictors as Section 2 and the 795 PHAA participants with complete data, of whom 126 have elevated blood pressure.

Why Linear Regression Does Not Suit a Binary Outcome

A binary outcome takes only two values, coded 1 for the event (elevated blood pressure) and 0 otherwise. A linear regression fitted to such an outcome can predict values below 0 or above 1, its residuals cannot be normally distributed, and their spread depends on the predicted probability. Logistic regression avoids all three problems.

Probability, Odds and Log Odds

For an outcome with probability p, the odds are p ÷ (1 − p). A probability of 0.5 corresponds to odds of 1, a probability of 0.2 to odds of 0.25, and a probability of 0.8 to odds of 4. The odds ratio (OR) divides the odds in one group by the odds in another. Among smokers, 48 of 135 have elevated blood pressure (odds 48 ÷ 87 = 0.552), and among non-smokers 80 of 663 do (odds 80 ÷ 583 = 0.137), so the crude odds ratio is 0.552 ÷ 0.137 = 4.02. The log odds, or logit, is the natural logarithm of the odds; it ranges from −∞ to +∞, which allows it to be modelled with a straight line.

Fitting the Model and Reading the Odds Ratios

The model is fitted with glm(elevated_bp ~ age + gender + smoker + bmi + dep_score, family = binomial, data = cc), where elevated_bp is a factor with levels "No" and "Yes", so that R models the probability of "Yes". The coefficients are on the log-odds scale, and exp(cbind(OR = coef(logit), confint(logit))) turns them into odds ratios with 95% confidence intervals.

The adjusted odds ratio for smoking (5.78) is larger than the crude odds ratio (4.02 in all 798 people with a reading, and 4.03 in the 795 people with complete data whom the model uses). This change should not be read as evidence of confounding. Smokers and non-smokers in the model have almost the same average age (45.5 and 45.1 years), yet adding age alone raises the odds ratio from 4.03 to 5.28, because an odds ratio changes when strong predictors of the outcome are added to the model, even when they are not confounders. Part of the change also reflects adjustment for BMI, which Section 4 treats as a mediator. The accordion below explains how each type of predictor is interpreted.

Checking a Logistic Model

Logistic regression has its own checks in place of the normality and equal-variance checks of the linear model. The first four parallel the linear model; the last three describe how well the model works as a whole.

The deviance measures how far a model is from fitting the data perfectly, and the likelihood ratio test compares the deviance of two nested models; here the deviance falls from 695.1 for the model with no predictors to 495.6 for the full model. McFadden's pseudo-R², 1 − 495.6 ÷ 695.1 = 0.29, summarises the improvement on a scale from 0 to 1, but it is not the share of variance explained and is usually much lower than R² from a linear model. The Hosmer-Lemeshow test, like other formal tests, becomes very sensitive in large samples, so the calibration plot is read alongside it.

phaa <- read.csv("phaa_survey_clean.csv")   # file must be in the working directory
phaa$smoker <- factor(phaa$smoker, levels = c("No", "Yes"))              # "No" = reference
phaa$gender <- factor(phaa$gender, levels = c("Man", "Woman", "Non-binary"))   # "Man" = reference
phaa$elevated_bp <- factor(ifelse(phaa$systolic_bp >= 120, "Yes", "No"),
                           levels = c("No", "Yes"))   # "Yes" (120 mmHg or more) is the event
table(phaa$hypertension)                       # the original hypertension variable
table(phaa$elevated_bp)                        # the outcome used in this section

cc <- na.omit(phaa[, c("elevated_bp", "age", "gender", "smoker", "bmi", "dep_score")])
logit <- glm(elevated_bp ~ age + gender + smoker + bmi + dep_score,
             family = binomial, data = cc)
summary(logit)

round(exp(cbind(OR = coef(logit), confint(logit))), 3)   # odds ratios with 95% CIs
round(exp(10 * coef(logit)[["age"]]), 2)                 # odds ratio for 10 years of age

two_people <- data.frame(age = 50, gender = "Man", smoker = c("No", "Yes"),
                         bmi = 21, dep_score = 21)   # same age, BMI and depression score
round(predict(logit, newdata = two_people, type = "response"), 3)   # predicted probabilities

# install.packages(c("ResourceSelection", "pROC"))   # run once, if not yet installed
library(ResourceSelection); library(pROC)
phaa <- read.csv("phaa_survey_clean.csv")   # file must be in the working directory
phaa$smoker <- factor(phaa$smoker, levels = c("No", "Yes"))              # "No" = reference
phaa$gender <- factor(phaa$gender, levels = c("Man", "Woman", "Non-binary"))   # "Man" = reference
phaa$elevated_bp <- factor(ifelse(phaa$systolic_bp >= 120, "Yes", "No"),
                           levels = c("No", "Yes"))   # "Yes" (120 mmHg or more) is the event
cc <- na.omit(phaa[, c("elevated_bp", "age", "gender", "smoker", "bmi", "dep_score")])
logit <- glm(elevated_bp ~ age + gender + smoker + bmi + dep_score,
             family = binomial, data = cc)
null_model <- glm(elevated_bp ~ 1, family = binomial, data = cc)
anova(null_model, logit, test = "Chisq")      # does the model beat no predictors at all?

logit_sq <- update(logit, . ~ . + I(age^2))
anova(logit, logit_sq, test = "Chisq")        # linearity on the log-odds scale

hoslem.test(as.numeric(cc$elevated_bp == "Yes"), fitted(logit), g = 10)   # calibration
round(1 - logit$deviance / logit$null.deviance, 3)                       # McFadden pseudo R-squared

roc_model <- roc(cc$elevated_bp, fitted(logit), levels = c("No", "Yes"), direction = "<")
auc(roc_model)                                 # discrimination: area under the ROC curve
plot(roc_model, legacy.axes = TRUE)           # x-axis shows 1 - specificity

Lesson 3 · Section 4

Building and Comparing Models

Choosing the variables for a model, and comparing one model with another.

Two goals

Explanation or prediction?

Explanation

The aim is to estimate the effect of an exposure, so the variables are chosen from subject-matter knowledge and confounders are kept.

Prediction

The aim is accurate forecasts for new people, so any variable that improves accuracy on new data can be kept.

Causal diagram

Confounders go in, mediators stay out

Overfitting

How many predictors can the data support?

Linear regression

A common guide is at least 10 to 20 observations for each estimated coefficient.

Logistic regression

A common guide is at least 10 events for each estimated coefficient.

A model with too many terms fits the noise in its sample and predicts poorly for new people.

Same data

Compare models on the same people

12

candidate variables, including the outcome

748

people with complete data for all of them

1

data set used for every comparison

Nested models

Do the extra terms improve the fit?

Model	Predictors	Adjusted R²
Smaller	age, gender, smoker, depression	0.382
Larger	smaller model + BMI + physical activity	0.396

anova(m_small, m_large) gives F = 9.54 on 2 and 740 df, p < 0.001.

AIC and BIC

Lower is better, and only differences matter

Model	AIC	BIC
Smaller (age, gender, smoker, depression)	5389.7	5422.1
Larger (+ BMI, physical activity)	5374.7	5416.3
Difference	15.0	5.8

Selection strategies

Backward elimination by AIC

What it kept

From 11 candidate predictors, step() kept age, smoking, depression, BMI and physical activity.

What it cannot do

It cannot tell a confounder from a mediator, and it dropped gender, a confounder in the causal diagram.

Worked example

The research question determines the model

Interactions

Does the smoking effect depend on age?

lm(systolic_bp ~ age * smoker + gender + dep_score) compared with the model without age:smoker gives F = 1.08, p = 0.30.

Rule 1

An interaction term is always fitted together with both of its main effects.

Rule 2

Interactions are chosen in advance from subject-matter knowledge, which limits false positives.

Lesson summary

From two-variable tests to model building

Section	Main tool	Main check
1. Two-variable tests	Correlation, t-test, ANOVA, chi-square and their rank-based versions	Is the comparison confounded?
2. Linear regression	lm(), adjusted differences in means	Four diagnostic plots, VIF, Cook's distance
3. Logistic regression	glm(family = binomial), odds ratios	Events per coefficient, calibration, AUC
4. Model building	Causal diagram, F-test, LRT, AIC, BIC	The same people in every model compared

Introduction and Overview

This section describes how to decide which variables enter a regression model and how to compare models. It uses the PHAA survey with systolic blood pressure as the outcome and smoking as the exposure of interest, and every comparison uses the 748 participants with complete data for twelve candidate variables.

Two Goals for a Model

A model built to explain estimates the effect of an exposure on an outcome; its variables are chosen to remove confounding, and its coefficients are interpreted. A model built to predict forecasts the outcome for new people; its variables are chosen for accuracy, and its individual coefficients matter less. The cards below summarise the two goals and the principle of parsimony, which applies to both.

Starting From a Causal Diagram

For an explanatory question, the variables are chosen before the results are seen, from subject-matter knowledge expressed as a causal diagram (also called a directed acyclic graph, or DAG). Each arrow is an assumed causal effect. A confounder has arrows to both the exposure and the outcome and is included in the model. A mediator lies on a path from the exposure to the outcome and is left out when the total effect of the exposure is of interest, because adjusting for it removes the part of the effect that travels through it.

The steps below describe the full model-building process, from specifying the largest model worth considering to presenting the result.

How Many Predictors Can the Data Support?

Every estimated coefficient uses up information. A common guide asks for at least 10 to 20 observations per coefficient in a linear model and at least 10 events per coefficient in a logistic model. A categorical predictor with k categories uses k − 1 coefficients. When a model has too many terms for its data, it starts to fit the random noise of its sample, which is called overfitting: its R² looks excellent, but its predictions for new people are poor. The simulator below shows the pattern.

Comparing Models

Models are compared only when they are fitted to the same observations. Because different predictors have different missing values, the safest approach is to create one complete-case data set for all candidate variables first, as the R activity does with na.omit(). The tabs below describe the main criteria.

Variable Selection Strategies

When many candidate variables are available, selection strategies apply a criterion such as AIC or a p-value threshold automatically.

In R, step(m_full, direction = "backward") performs backward elimination by AIC. Starting from all 11 candidate predictors (age, gender, smoking, depression, BMI, physical activity, anxiety, social support, discrimination, education and region), it keeps age, smoking, depression, BMI and physical activity. AIC-based elimination keeps any term whose removal would raise AIC, which corresponds roughly to keeping terms with a p-value below about 0.16, so physical activity (p = 0.14 in the full model) is retained. Gender is dropped, although the causal diagram treats it as a confounder.

Worked Example: The Smoking Estimate Across Models

Model	Smoking estimate (mmHg)	95% CI	AIC	What it estimates
Smoking only	6.15	4.03 to 8.28	5714.9	The crude difference, affected by any confounding.
Confounders from the causal diagram (age, gender, depression)	5.37	3.66 to 7.09	5389.7	The total effect of smoking, if the diagram is correct.
Confounders + BMI and physical activity	6.01	4.28 to 7.75	5374.7	The effect of smoking other than through BMI and activity.
Backward elimination by AIC	6.02	4.29 to 7.75	5371.7	A mixture chosen by fit, which drops gender and keeps the mediators.

The models with the mediators have the lower AIC, mainly because BMI predicts blood pressure well. In these data, smokers have a lower average BMI (19.3 against 21.1 among these 748 people), so part of the effect of smoking on blood pressure is counteracted through BMI, and adjusting for BMI removes that counteracting path and raises the smoking estimate. For the question "how much higher is blood pressure because of smoking?", the confounder model, with an estimate of 5.37 mmHg, is the appropriate one. Fit statistics then help to choose among models that answer the same question, for example a straight-line term for age against a curved one.

Interactions

An interaction (effect modification) is present when the effect of one predictor depends on the value of another. In R, age * smoker adds age, smoking and their product age:smoker. Comparing lm(systolic_bp ~ age * smoker + gender + dep_score) with the model without the product gives F = 1.08 on 1 and 741 df (p = 0.30), so there is no evidence that the effect of smoking changes with age, and the simpler model is kept. An interaction term is always fitted with both of its main effects, and the interactions to test are best chosen in advance; the accordion below lists common strategies, of which the theory-driven ones are preferred.

Reporting the Final Model

A report of a regression model states the goal of the analysis and how the variables were chosen (with the causal diagram, if one was used), the number of people analysed out of the number available, each estimate on its own scale (a difference in means for a linear model, an odds ratio for a logistic model) with its 95% confidence interval, the checks that were carried out and their results, and any sensitivity analyses, such as the same model with an alternative set of variables.

phaa <- read.csv("phaa_survey_clean.csv")   # file must be in the working directory
phaa$smoker <- factor(phaa$smoker, levels = c("No", "Yes"))              # "No" = reference
phaa$gender <- factor(phaa$gender, levels = c("Man", "Woman", "Non-binary"))   # "Man" = reference
vars <- c("systolic_bp", "age", "gender", "smoker", "dep_score", "bmi", "phys_act_min",
          "anx_score", "social_support_score", "discrimination_score", "education", "region")
build <- na.omit(phaa[, vars])                 # the same people for every model compared
nrow(build)

m_small <- lm(systolic_bp ~ age + gender + smoker + dep_score, data = build)
m_large <- lm(systolic_bp ~ age + gender + smoker + dep_score + bmi + phys_act_min, data = build)
anova(m_small, m_large)                       # nested F-test: do the added terms help?
round(c(adjR2_small = summary(m_small)$adj.r.squared,
        adjR2_large = summary(m_large)$adj.r.squared), 3)
AIC(m_small, m_large)                         # lower AIC is better
BIC(m_small, m_large)                         # BIC penalises extra terms more

m_full <- lm(systolic_bp ~ ., data = build)  # every candidate predictor
m_back <- step(m_full, direction = "backward", trace = 0)   # backward elimination by AIC
formula(m_back)

m_int <- lm(systolic_bp ~ age * smoker + gender + dep_score, data = build)   # age:smoker added
anova(m_small, m_int)                         # does the smoking effect differ by age?

Bringing It All Together

This lesson moved from comparing two variables to building and comparing regression models. Section 1 matched each pair of variable types to a two-variable test, listed the rank-based and exact alternatives used when an outcome is skewed or ordinal or when counts are small, applied the correlation coefficient, the t-test, ANOVA and the chi-square test to the PHAA survey, and showed with physical activity and blood pressure that a two-variable comparison can be distorted by confounding: the difference of 2.8 mmHg per 100 minutes of activity fell to 0.9 mmHg once age was taken into account. It then introduced the regression line, least squares and the multivariable model, in which each coefficient is adjusted for the other predictors, and showed that the t-test is a regression with one binary predictor.

Section 2 fitted a linear regression for systolic blood pressure with five predictors, in which smokers averaged 6.1 mmHg higher than comparable non-smokers and the model explained 41% of the variation. The model passed each check: the four diagnostic plots, a squared age term, the Breusch-Pagan and Shapiro-Wilk tests, Cook's distance and variance inflation factors. The section also described the standard responses when a check fails, such as a squared term, a log transformation or sandwich standard errors. Section 3 modelled elevated blood pressure with logistic regression, in which smokers had 5.8 times the odds of comparable non-smokers, explained odds, log odds and predicted probabilities, and checked the model for events per coefficient, linearity on the log-odds scale, calibration (Hosmer-Lemeshow p = 0.59) and discrimination (AUC = 0.86). Section 4 showed that the choice of variables starts from the goal of the analysis and, for an explanatory question, from a causal diagram; that models are compared on the same people with F-tests, likelihood ratio tests, AIC and BIC; and that the best-fitting model, which included two mediators of smoking, was not the right model for estimating the total effect of smoking.

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