The V Lab
AudiencePublic health, epidemiology, and health science learners
Study timeAbout 120–180 minutes
PrerequisitesDescriptive statistics, confidence intervals, and basic R

About the tutorial data Every record is simulated with a fixed random seed and contains no real personal health information. The data deliberately include mild nonlinearity, heteroskedasticity, and missingness for diagnostic practice. Simulated relationships are teaching devices, not estimates for a real population.

How to use this tutorial

Work through the sequence question → data check → fit → diagnose → interpret → validate → report. Read each explanation before running its code, and attempt the exercises before opening the answers. Code is shown by default and can be folded with the page controls.

Learning objectives

By the end of this tutorial, you should be able to:

  • decide whether linear regression matches the research question, outcome, and data structure;
  • interpret intercepts and coefficients for numeric and categorical predictors in simple and multiple regression;
  • understand ordinary least squares through its geometric and loss-function intuition;
  • report coefficients, 95% confidence intervals, p-values, R2R^2, and adjusted R2R^2;
  • use meaningful units, interactions, and nonlinear terms to answer better-specified questions;
  • check residuals, normal Q–Q plots, heteroskedasticity, influence, and multicollinearity;
  • calculate HC3 heteroskedasticity-robust standard errors using base R;
  • distinguish a confidence interval for a mean response from a prediction interval for an individual;
  • assess training and test performance while avoiding overfitting and data leakage;
  • explain limitations arising from missing data, study design, and causal interpretation; and
  • write a clear, reproducible report without overstating the evidence.

1 From a Research Question to a Linear Model

1.1 When should linear regression be used?

Linear regression primarily describes how the conditional mean of a continuous outcome changes with one or more predictors. Examples include:

  • Does mean systolic blood pressure vary with age?
  • How is BMI associated with mean systolic blood pressure among people with the same modeled age, smoking status, and neighborhood?
  • After adjustment, how much does mean systolic blood pressure differ between current and non-current smokers?
  • Does the age–blood pressure association differ by smoking status?

A general model is:

Yi=β0+β1Xi1+⋯+βpXip+εi,E(εi∣Xi)=0. Y_i = \beta_0 + \beta_1X_{i1} + \cdots + \beta_pX_{ip} + \varepsilon_i, \qquad E(\varepsilon_i\mid X_i)=0.

“Linear” first means linear in the unknown parameters β\beta. Predictors may be squared, segmented, or transformed in prespecified ways; linear regression does not require every fitted relationship to be a straight line.

Identify the outcome before choosing a function Ordinary linear regression is usually not the first choice for a binary event, count, proportion, or survival time. Match the method to the outcome distribution, estimand, and sampling design rather than using it simply because lm() is convenient.

1.2 Prediction, description, and causation are different tasks

The same formula can serve different goals, but its evaluation criteria change:

Task Main question Primary emphasis
Describe an association How does the conditional mean vary in the observed sample? Coefficients, intervals, functional form, and transparent reporting
Predict How accurately can the model predict a new, unobserved person? Out-of-sample error, calibration, and target population
Estimate a causal effect How would the outcome change under an intervention on the exposure? Design, temporality, confounding control, and identification assumptions

A regression that “adjusts for several covariates” does not automatically estimate a causal effect. Unless a defensible design and causal assumptions support stronger language, all interpretations below concern conditional associations.

1.3 Ordinary least squares intuition

For participant ii, the model gives fitted value Ŷi\hat Y_i and residual ei=Yi−Ŷie_i=Y_i-\hat Y_i. Ordinary least squares (OLS) chooses coefficients that minimize the sum of squared residuals:

SSE=∑i=1n(Yi−Ŷi)2=∑i=1nei2. \mathrm{SSE}=\sum_{i=1}^{n}(Y_i-\hat Y_i)^2=\sum_{i=1}^{n}e_i^2.

Squaring prevents positive and negative residuals from cancelling and penalizes large deviations more strongly. OLS point estimates do not require the outcome itself to be normally distributed; normal residuals matter more directly for the exactness of classical small-sample t tests and intervals.

1.3.1 Matrix view (optional)

When the design matrix XX has full column rank, the OLS solution is:

𝛃̂=(XTX)−1XT𝐲. \hat{\boldsymbol\beta}=(X^TX)^{-1}X^T\mathbf y.

Perfect collinearity prevents a unique solution because XTXX^TX is then singular. Statistical software uses more numerically stable decompositions; fitting a model does not require manually inverting this matrix.

2 Meet the Simulated Data

2.1 Variables and unit of analysis

Each row represents one simulated participant. Systolic blood pressure in mmHg is the primary outcome; candidate predictors include age, BMI, current smoking status, sex, neighborhood, and weekly physical activity. Before any outcome exploration, 25% of the complete simulated data was randomly held out. From this point, ph_data contains development data only; the holdout is not summarized until the out-of-sample section.

str(ph_data)
## 'data.frame':    420 obs. of  11 variables:
##  $ participant_id            : chr  "P001" "P004" "P005" "P006" ...
##  $ age                       : num  23 31 38 67 44 28 33 30 48 53 ...
##  $ sex                       : Factor w/ 2 levels "Female","Male": 1 1 2 1 1 2 1 2 2 1 ...
##  $ neighborhood              : Factor w/ 4 levels "Central","North",..: 1 3 4 1 1 4 2 1 4 4 ...
##  $ smoking_status            : Factor w/ 2 levels "Not current",..: 1 1 1 1 1 2 1 1 1 1 ...
##  $ bmi                       : num  27.8 20.4 26.5 18.8 27.6 29 18.6 26.3 27.9 16.8 ...
##  $ physical_activity_min_week: num  113 69 116 105 32 361 NA 93 121 84 ...
##  $ systolic_bp               : num  117 109 120 151 111 ...
##  $ age_c10                   : num  -2.7 -1.9 -1.2 1.7 -0.6 -2.2 -1.7 -2 -0.2 0.3 ...
##  $ bmi_c5                    : num  0.56 -0.92 0.3 -1.24 0.52 0.8 -1.28 0.26 0.58 -1.64 ...
##  $ physical_activity_30      : num  3.77 2.3 3.87 3.5 1.07 ...

Record variable names, roles, and units in the analysis plan before examining results.

Variable Type Role or unit in this tutorial
systolic_bp Continuous Outcome, mmHg
age Continuous Years
bmi Continuous kg/m²
smoking_status Binary Current versus not current
sex Binary Male versus female
neighborhood Multicategory North, South, or West versus Central
physical_activity_min_week Continuous Minutes/week; contains missing values

2.2 Integrity checks before modeling

data_check <- data.frame(
  "Variable" = names(ph_data),
  "Class" = vapply(ph_data, function(x) class(x)[1], character(1)),
  "Missing" = vapply(ph_data, function(x) sum(is.na(x)), integer(1)),
  "Unique values" = vapply(
    ph_data,
    function(x) length(unique(x[!is.na(x)])),
    integer(1)
  ),
  check.names = FALSE
)
knitr::kable(data_check, caption = "Variable types, missingness, and unique values")
Variable types, missingness, and unique values
Variable Class Missing Unique values
participant_id participant_id character 0 420
age age numeric 0 67
sex sex factor 0 2
neighborhood neighborhood factor 0 4
smoking_status smoking_status factor 0 2
bmi bmi numeric 0 152
physical_activity_min_week physical_activity_min_week numeric 29 201
systolic_bp systolic_bp numeric 0 288
age_c10 age_c10 numeric 0 67
bmi_c5 bmi_c5 numeric 0 152
physical_activity_30 physical_activity_30 numeric 29 201

Also check impossible values, duplicate identifiers, incorrect units, coding changes, and eligibility criteria. An automatic summary() is useful but cannot replace a data dictionary and subject-matter knowledge.

summary(ph_data[c(
  "age", "bmi", "physical_activity_min_week", "systolic_bp",
  "smoking_status", "sex", "neighborhood"
)])
##       age            bmi       physical_activity_min_week  systolic_bp   
##  Min.   :18.0   Min.   :16.0   Min.   :  0                Min.   : 91.7  
##  1st Qu.:35.0   1st Qu.:22.2   1st Qu.: 47                1st Qu.:117.6  
##  Median :47.0   Median :24.9   Median : 94                Median :127.0  
##  Mean   :46.6   Mean   :25.0   Mean   :110                Mean   :128.0  
##  3rd Qu.:58.0   3rd Qu.:27.7   3rd Qu.:154                3rd Qu.:136.7  
##  Max.   :85.0   Max.   :38.5   Max.   :535                Max.   :170.7  
##                                NAs    :29                                
##      smoking_status     sex       neighborhood
##  Not current:353    Female:201   Central:136  
##  Current    : 67    Male  :219   North  : 85  
##                                  South  :110  
##                                  West   : 89  
##                                               
##                                               
##

2.3 Exploratory graphs: inspect shape before fitting a line

plot(
  ph_data$age, ph_data$systolic_bp,
  pch = 16, cex = 0.65,
  col = rgb(0, 114 / 255, 178 / 255, 0.38),
  xlab = "Age (years)", ylab = "Systolic blood pressure (mmHg)",
  main = "Check whether a linear relationship is plausible"
)
abline(
  lm(systolic_bp ~ age, data = ph_data),
  col = palette_ph["vermillion"], lwd = 2.5
)
lines(
  lowess(ph_data$age, ph_data$systolic_bp, f = 2 / 3),
  col = palette_ph["teal"], lwd = 2.5, lty = 2
)
legend(
  "topleft",
  legend = c("OLS line", "LOWESS smoother"),
  col = c(palette_ph["vermillion"], palette_ph["teal"]),
  lwd = 2.5, lty = c(1, 2), bty = "n"
)
Scatterplot of age and systolic blood pressure with a straight fitted line and a smooth curve.

Systolic blood pressure versus age. The straight line is a simple linear fit; the curve is a LOWESS smoother.

The smoother is a diagnostic for functional form, not the final model. Systematic separation between the smoother and straight line suggests considering scientifically meaningful nonlinear terms rather than blindly adding high-degree polynomials.

boxplot(
  systolic_bp ~ smoking_status,
  data = ph_data,
  names = c("Not current", "Current"),
  col = c("#D9EAF3", "#F4C6A6"),
  border = palette_ph["navy"],
  ylab = "Systolic blood pressure (mmHg)", xlab = "",
  main = "Distribution across a categorical predictor"
)
stripchart(
  systolic_bp ~ smoking_status,
  data = ph_data,
  vertical = TRUE, method = "jitter", pch = 16,
  col = rgb(23 / 255, 50 / 255, 77 / 255, 0.22),
  add = TRUE
)
Boxplots and points for systolic blood pressure among current and non-current smokers.

Systolic blood pressure by smoking status; jittered observations are overlaid on the boxplots.

The crude group difference may reflect differences in age, BMI, or neighborhood composition. Multiple regression can describe a conditional mean difference at the same modeled covariate values, but its credibility still depends on model specification and study design.

3 Simple Linear Regression

3.1 A model with one predictor

Begin by relating mean systolic blood pressure to age:

E(SBP∣age)=β0+β1age. E(\mathrm{SBP}\mid \mathrm{age})=\beta_0+\beta_1\mathrm{age}. 

simple_model <- lm(systolic_bp ~ age, data = ph_data)
summary(simple_model)
## 
## Call:
## lm(formula = systolic_bp ~ age, data = ph_data)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -35.45  -7.84   0.70   7.41  32.12 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 101.6771     1.7177    59.2   <2e-16 ***
## age           0.5660     0.0349    16.2   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.3 on 418 degrees of freedom
## Multiple R-squared:  0.386,  Adjusted R-squared:  0.384 
## F-statistic:  263 on 1 and 418 DF,  p-value: <2e-16
knitr::kable(
  coefficient_table(simple_model),
  digits = 3,
  caption = "Simple linear regression of systolic blood pressure on age"
)
Simple linear regression of systolic blood pressure on age
Term Estimate Standard error 95% CI lower 95% CI upper p-value
Intercept 101.677 1.718 98.301 105.054 0
Age (per 1 year) 0.566 0.035 0.497 0.635 0

The age coefficient is 0.57 mmHg per year (95% CI 0.5 to 0.63). It describes the mean blood-pressure difference between sampled people one year apart in age; it does not prove that changing age itself causes that difference.

The intercept extrapolates the mean to age zero, outside this adult sample. It mainly positions the line and usually has no substantive interpretation here. Centering age can make the intercept correspond to a meaningful reference age.

3.2 Observed values, fitted values, and residuals

simple_components <- data.frame(
  "Participant ID" = ph_data$participant_id[1:8],
  "Observed" = ph_data$systolic_bp[1:8],
  "Fitted" = fitted(simple_model)[1:8],
  "Residual" = residuals(simple_model)[1:8],
  check.names = FALSE
)
knitr::kable(
  simple_components,
  digits = 2,
  caption = "Observed values, fitted values, and residuals for the first 8 participants"
)
Observed values, fitted values, and residuals for the first 8 participants
Participant ID Observed Fitted Residual
1 P001 117 115 2.20
4 P004 109 119 -10.52
5 P005 120 123 -3.59
6 P006 151 140 11.70
7 P007 111 127 -15.18
8 P008 138 118 20.07
9 P009 111 120 -9.46
10 P010 122 119 3.34

In an OLS model with an intercept, residuals sum to zero up to numerical error. This does not mean the model predicts every person accurately.

c(
  residual_sum = sum(residuals(simple_model)),
  residual_mean = mean(residuals(simple_model)),
  residual_rmse = sqrt(mean(residuals(simple_model)^2))
)
##  residual_sum residual_mean residual_rmse 
##      3.55e-13      1.27e-15      1.13e+01

4 Multiple Linear Regression

4.1 Numeric and categorical predictors together

multiple_model <- lm(
  systolic_bp ~ age + bmi + smoking_status + sex + neighborhood,
  data = ph_data
)

knitr::kable(
  coefficient_table(multiple_model),
  digits = 3,
  caption = "Multiple linear regression of systolic blood pressure"
)
Multiple linear regression of systolic blood pressure
Term Estimate Standard error 95% CI lower 95% CI upper p-value
Intercept 84.046 3.500 77.165 90.926 0.000
Age (per 1 year) 0.557 0.033 0.491 0.622 0.000
BMI (per 1 unit) 0.640 0.126 0.393 0.886 0.000
Current smoking (reference: not current) 4.209 1.433 1.392 7.025 0.003
Male (reference: female) 2.083 1.058 0.003 4.162 0.050
North (reference: Central) 2.840 1.482 -0.074 5.754 0.056
South (reference: Central) -2.728 1.369 -5.419 -0.037 0.047
West (reference: Central) 2.262 1.460 -0.609 5.132 0.122

A numeric coefficient is the conditional mean difference for a one-unit increase holding the other modeled variables fixed. The age coefficient, for example, compares participants one year apart who have the same modeled BMI, smoking status, sex, and neighborhood.

“Holding fixed” describes a conditional comparison within the model. It does not guarantee that perfectly matched people exist in the sample, nor does it make the comparison causal.

4.2 Categorical predictors and reference levels

R uses the first factor level as the reference by default. Set it deliberately before fitting rather than changing it in response to whichever comparison looks most significant.

list(
  smoking_status = levels(ph_data$smoking_status),
  sex = levels(ph_data$sex),
  neighborhood = levels(ph_data$neighborhood)
)
## $smoking_status
## [1] "Not current" "Current"    
## 
## $sex
## [1] "Female" "Male"  
## 
## $neighborhood
## [1] "Central" "North"   "South"   "West"

Thus, smoking_statusCurrent compares current with non-current smoking. The three neighborhood coefficients compare North, South, and West with Central. An unordered factor with kk levels generally produces k−1k-1 coefficients.

To use North as the reference when that comparison matches the question:

ph_reference_demo <- ph_data
ph_reference_demo$neighborhood <- relevel(
  ph_reference_demo$neighborhood,
  ref = "North"
)
reference_demo_model <- lm(
  systolic_bp ~ age + bmi + smoking_status + sex + neighborhood,
  data = ph_reference_demo
)
coef(reference_demo_model)[grep("neighborhood", names(coef(reference_demo_model)))]
## neighborhoodCentral   neighborhoodSouth    neighborhoodWest 
##              -2.840              -5.568              -0.578

Changing the reference level changes how selected coefficients are expressed, but not participant-level fitted values, residuals, or overall fit.

4.3 Rescaling and centering in meaningful units

Effects per one year or one BMI unit can be too granular. This tutorial defines:

  • age_c10 = (age - 50) / 10: per 10 years, with zero at age 50;
  • bmi_c5 = (bmi - 25) / 5: per 5 BMI units, with zero at BMI 25.
scaled_model <- lm(
  systolic_bp ~ age_c10 + bmi_c5 + smoking_status + sex + neighborhood,
  data = ph_data
)

knitr::kable(
  coefficient_table(scaled_model),
  digits = 3,
  caption = "Multiple model expressed in meaningful units"
)
Multiple model expressed in meaningful units
Term Estimate Standard error 95% CI lower 95% CI upper p-value
Intercept 127.86 1.085 125.730 129.997 0.000
Age (per 10 years; centered at 50) 5.57 0.335 4.907 6.224 0.000
BMI (per 5 units; centered at 25) 3.20 0.628 1.963 4.432 0.000
Current smoking (reference: not current) 4.21 1.433 1.392 7.025 0.003
Male (reference: female) 2.08 1.058 0.003 4.162 0.050
North (reference: Central) 2.84 1.482 -0.074 5.754 0.056
South (reference: Central) -2.73 1.369 -5.419 -0.037 0.047
West (reference: Central) 2.26 1.460 -0.609 5.132 0.122

After adjustment for the other modeled variables, a 10-year age difference corresponds to a mean systolic blood-pressure difference of 5.57 mmHg. A 5-unit BMI difference corresponds to 3.2 mmHg. Centering leaves model fit unchanged but makes the intercept correspond to a 50-year-old, BMI-25, non-current-smoking female participant in Central.

4.4 Confidence intervals for coefficients

A point estimate is the value that best fits this sample; its confidence interval expresses repeated-sampling uncertainty. Interval width depends on sample size, residual variability, the predictor distribution, and collinearity.

scaled_ci <- confint(scaled_model)
keep_coef <- rownames(scaled_ci) != "(Intercept)"
coef_values <- coef(scaled_model)[keep_coef]
ci_values <- scaled_ci[keep_coef, , drop = FALSE]
plot(
  coef_values, seq_along(coef_values),
  xlim = range(ci_values), yaxt = "n",
  pch = 19, col = palette_ph["blue"],
  xlab = "Mean systolic blood-pressure difference (mmHg)", ylab = "",
  main = "Coefficients and 95% confidence intervals"
)
segments(
  ci_values[, 1], seq_along(coef_values),
  ci_values[, 2], seq_along(coef_values),
  col = palette_ph["navy"], lwd = 2
)
abline(v = 0, lty = 2, col = "grey45")
axis(
  2, at = seq_along(coef_values),
  labels = label_terms(names(coef_values)), las = 1, cex.axis = 0.72
)
A coefficient plot with point estimates and horizontal confidence intervals.

Multiple-regression coefficients and 95% confidence intervals; the intercept is omitted for readability.

Do not equate “the interval includes zero” with “there is no association,” or “the interval excludes zero” with “the association matters for public health.” Interpret the magnitude, compatible range, units, and context together.

5 Model Fit: R² and Adjusted R²

5.1 What do these measures answer?

R2=1−∑i(Yi−Ŷi)2∑i(Yi−Y‾)2. R^2=1-\frac{\sum_i(Y_i-\hat Y_i)^2}{\sum_i(Y_i-\bar Y)^2}.

R2R^2 is the proportion of outcome variability explained by the model in the current sample. Adding a predictor cannot lower training-sample R2R^2, even when that predictor has no practical value. Adjusted R2R^2 penalizes complexity and can decline, but it still does not replace out-of-sample evaluation.

fit_statistics <- function(model, model_name) {
  s <- summary(model)
  data.frame(
    "Model" = model_name,
    "Sample size" = nobs(model),
    "Parameters" = length(coef(model)),
    "R-squared" = s$r.squared,
    "Adjusted R-squared" = s$adj.r.squared,
    "Residual standard error" = s$sigma,
    check.names = FALSE
  )
}

fit_comparison <- rbind(
  fit_statistics(simple_model, "Age only"),
  fit_statistics(scaled_model, "Age + BMI + smoking + sex + neighborhood")
)
knitr::kable(fit_comparison, digits = 3, caption = "In-sample fit measures")
In-sample fit measures
Model Sample size Parameters R-squared Adjusted R-squared Residual standard error
Age only 420 2 0.386 0.384 11.3
Age + BMI + smoking + sex + neighborhood 420 8 0.461 0.452 10.6

A low R2R^2 does not invalidate a well-estimated association, and a high R2R^2 does not prove that a model is unbiased, transportable, or causal.

6 Interaction: Does an Association Differ Across Groups?

6.1 Fit an age × smoking-status interaction

A no-interaction model assumes the same age slope in both smoking groups. An interaction permits a different slope among current smokers:

interaction_model <- lm(
  systolic_bp ~ age_c10 * smoking_status +
    bmi_c5 + sex + neighborhood,
  data = ph_data
)

knitr::kable(
  coefficient_table(interaction_model),
  digits = 3,
  caption = "Model with an age-by-smoking-status interaction"
)
Model with an age-by-smoking-status interaction
Term Estimate Standard error 95% CI lower 95% CI upper p-value
Intercept 127.89 1.082 125.760 130.012 0.000
Age (per 10 years; centered at 50) 5.30 0.359 4.599 6.010 0.000
Current smoking (reference: not current) 4.30 1.429 1.491 7.107 0.003
BMI (per 5 units; centered at 25) 3.21 0.626 1.983 4.444 0.000
Male (reference: female) 2.07 1.054 0.000 4.145 0.050
North (reference: Central) 2.61 1.482 -0.301 5.524 0.079
South (reference: Central) -2.84 1.365 -5.525 -0.157 0.038
West (reference: Central) 2.05 1.459 -0.823 4.914 0.162
Age × current smoking 1.93 0.976 0.010 3.847 0.049

In this model:

  • age_c10 is the 10-year age slope among non-current smokers;
  • smoking_statusCurrent is the group difference when age_c10 = 0, or at age 50;
  • the interaction coefficient is the difference between the two age slopes; and
  • the age slope among current smokers is the age main effect plus the interaction.

6.2 Obtain group-specific slopes as linear combinations

linear_combination_ci <- function(model, weights, level = 0.95) {
  beta <- coef(model)
  V <- vcov(model)
  L <- setNames(rep(0, length(beta)), names(beta))
  L[names(weights)] <- weights
  estimate <- sum(L * beta)
  se <- sqrt(drop(t(L) %*% V %*% L))
  critical <- qt(1 - (1 - level) / 2, df.residual(model))
  c(
    estimate = estimate,
    standard_error = se,
    lower = estimate - critical * se,
    upper = estimate + critical * se
  )
}

slope_not_current <- linear_combination_ci(
  interaction_model,
  c(age_c10 = 1)
)
slope_current <- linear_combination_ci(
  interaction_model,
  c(age_c10 = 1, "age_c10:smoking_statusCurrent" = 1)
)
interaction_slopes <- rbind(
  "Not current" = slope_not_current,
  "Current" = slope_current
)
knitr::kable(
  interaction_slopes,
  digits = 2,
  caption = "Ten-year systolic blood-pressure slopes by smoking status (mmHg)"
)
Ten-year systolic blood-pressure slopes by smoking status (mmHg)
estimate standard_error lower upper
Not current 5.30 0.36 4.60 6.01
Current 7.23 0.91 5.45 9.02 
age_grid <- seq(min(ph_data$age), max(ph_data$age), length.out = 120)
interaction_grid <- rbind(
  data.frame(
    age = age_grid,
    age_c10 = (age_grid - 50) / 10,
    bmi_c5 = 0,
    smoking_status = factor("Not current", levels = levels(ph_data$smoking_status)),
    sex = factor("Female", levels = levels(ph_data$sex)),
    neighborhood = factor("Central", levels = levels(ph_data$neighborhood))
  ),
  data.frame(
    age = age_grid,
    age_c10 = (age_grid - 50) / 10,
    bmi_c5 = 0,
    smoking_status = factor("Current", levels = levels(ph_data$smoking_status)),
    sex = factor("Female", levels = levels(ph_data$sex)),
    neighborhood = factor("Central", levels = levels(ph_data$neighborhood))
  )
)
interaction_grid$prediction <- predict(interaction_model, newdata = interaction_grid)

plot(
  systolic_bp ~ age, data = ph_data,
  pch = 16, cex = 0.45, col = "grey80",
  xlab = "Age (years)", ylab = "Predicted mean systolic BP (mmHg)",
  main = "Interpret an interaction graphically"
)
for (group_index in seq_along(levels(ph_data$smoking_status))) {
  rows <- interaction_grid$smoking_status ==
    levels(ph_data$smoking_status)[group_index]
  lines(
    interaction_grid$age[rows], interaction_grid$prediction[rows],
    col = c(palette_ph["blue"], palette_ph["vermillion"])[group_index],
    lwd = 3
  )
}
legend(
  "topleft", legend = c("Not current", "Current"),
  col = c(palette_ph["blue"], palette_ph["vermillion"]),
  lwd = 3, bty = "n"
)
Two fitted age and predicted blood-pressure lines, one for each smoking group.

Model-predicted age–blood-pressure relationships by smoking status, with BMI fixed at 25, sex at female, and neighborhood at Central.

The hierarchy principle When a model contains X×ZX\times Z, retain the XX and ZZ main effects in most settings, even if one has a large p-value. Let a prespecified scientific question motivate an interaction, and explain it through stratified predictions or marginal effects over a meaningful covariate range.

7 Nonlinear Relationships

7.1 Represent curvature with a centered quadratic term

If the age slope is not constant, add the square of centered age:

linear_age_model <- lm(
  systolic_bp ~ age_c10 + bmi_c5 + smoking_status + sex + neighborhood,
  data = ph_data
)
quadratic_model <- lm(
  systolic_bp ~ age_c10 + I(age_c10^2) +
    bmi_c5 + smoking_status + sex + neighborhood,
  data = ph_data
)

knitr::kable(
  coefficient_table(quadratic_model),
  digits = 3,
  caption = "Model with a centered quadratic age term"
)
Model with a centered quadratic age term
Term Estimate Standard error 95% CI lower 95% CI upper p-value
Intercept 125.556 1.154 123.286 127.825 0.000
Age (per 10 years; centered at 50) 6.152 0.347 5.471 6.834 0.000
Age squared 0.889 0.180 0.536 1.243 0.000
BMI (per 5 units; centered at 25) 3.225 0.611 2.023 4.426 0.000
Current smoking (reference: not current) 4.566 1.396 1.823 7.310 0.001
Male (reference: female) 2.630 1.035 0.595 4.665 0.011
North (reference: Central) 3.198 1.444 0.360 6.036 0.027
South (reference: Central) -3.307 1.337 -5.935 -0.679 0.014
West (reference: Central) 1.940 1.422 -0.855 4.735 0.173

With a quadratic term present, age_c10 is the local slope near age 50 rather than one slope shared across all ages. Do not interpret the squared term by itself; graph predictions using both age terms.

knitr::kable(
  as.data.frame(anova(linear_age_model, quadratic_model)),
  digits = 3,
  caption = "Comparison of nested linear-age and quadratic-age models"
)
Comparison of nested linear-age and quadratic-age models
Res.Df RSS Df Sum of Sq F Pr(>F)
412 46714 NA NA NA NA
411 44094 1 2620 24.4 0

This F test compares two nested models; it cannot establish that a quadratic curve is the true mechanism. Judge functional form using graphs, prediction goals, prior knowledge, and external validation as well.

quadratic_grid <- data.frame(
  age_c10 = (age_grid - 50) / 10,
  bmi_c5 = 0,
  smoking_status = factor("Not current", levels = levels(ph_data$smoking_status)),
  sex = factor("Female", levels = levels(ph_data$sex)),
  neighborhood = factor("Central", levels = levels(ph_data$neighborhood))
)
quadratic_prediction <- predict(
  quadratic_model,
  newdata = quadratic_grid,
  interval = "confidence"
)
plot(
  age_grid, quadratic_prediction[, "fit"], type = "n",
  ylim = range(quadratic_prediction),
  xlab = "Age (years)", ylab = "Predicted mean systolic BP (mmHg)",
  main = "Graph the curve rather than guessing from one coefficient"
)
polygon(
  c(age_grid, rev(age_grid)),
  c(quadratic_prediction[, "lwr"], rev(quadratic_prediction[, "upr"])),
  col = rgb(86 / 255, 180 / 255, 233 / 255, 0.28), border = NA
)
lines(
  age_grid, quadratic_prediction[, "fit"],
  col = palette_ph["blue"], lwd = 3
)
A curved age and predicted mean blood-pressure line surrounded by a shaded confidence band.

Adjusted quadratic relationship between age and predicted mean systolic blood pressure with a 95% confidence band.

High-degree polynomials can oscillate sharply near data boundaries, making extrapolation especially dangerous. More complex relationships may call for piecewise lines, splines, or subject-matter transformations, with complexity matched to sample size and purpose.

8 Modeling Assumptions and Diagnostics

8.1 Assumption checklist

Classical linear-regression inference commonly considers:

  1. A reasonable conditional-mean specification: omitted nonlinearities, interactions, or structure do not remain systematically in the residuals.
  2. Independent observations: repeated measures and clustering within families, schools, neighborhoods, or hospitals require specific treatment.
  3. Homoskedasticity: conditional error variability is approximately constant.
  4. Approximately normal residuals: this matters chiefly for exact classical small-sample intervals and tests; predictors need not be normal.
  5. No perfect collinearity: no design-matrix column is an exact combination of the others.
  6. Credible measurement and sampling: residual plots cannot automatically reveal serious measurement error, selection bias, or inappropriate weights.
  7. A defined domain of use: do not extrapolate to populations or numeric ranges unsupported by the sample.

Diagnostics are not a one-time pass/fail examination. They help identify disagreement between the model and data, assess sensitivity, and improve reporting.

8.2 A working model for complete diagnostics

This working model allows both age curvature and an age-by-smoking interaction:

working_model <- lm(
  systolic_bp ~ age_c10 * smoking_status + I(age_c10^2) +
    bmi_c5 + sex + neighborhood,
  data = ph_data
)
summary(working_model)
## 
## Call:
## lm(formula = systolic_bp ~ age_c10 * smoking_status + I(age_c10^2) + 
##     bmi_c5 + sex + neighborhood, data = ph_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -30.272  -6.761   0.021   6.692  26.530 
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    125.623      1.152  109.02  < 2e-16 ***
## age_c10                          5.916      0.372   15.92  < 2e-16 ***
## smoking_statusCurrent            4.636      1.393    3.33  0.00095 ***
## I(age_c10^2)                     0.871      0.180    4.84  1.8e-06 ***
## bmi_c5                           3.238      0.610    5.31  1.8e-07 ***
## sexMale                          2.610      1.033    2.53  0.01187 *  
## neighborhoodNorth                2.994      1.445    2.07  0.03889 *  
## neighborhoodSouth               -3.392      1.334   -2.54  0.01140 *  
## neighborhoodWest                 1.762      1.422    1.24  0.21618    
## age_c10:smoking_statusCurrent    1.657      0.952    1.74  0.08253 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.3 on 410 degrees of freedom
## Multiple R-squared:  0.495,  Adjusted R-squared:  0.484 
## F-statistic: 44.6 on 9 and 410 DF,  p-value: <2e-16

8.3 Residual, Q–Q, scale–location, and influence plots

old_par <- par(mfrow = c(2, 2), mar = c(4, 4, 2, 1))
plot(working_model, which = 1:4, caption = rep("", 4))
Residual-versus-fitted, normal Q-Q, scale-location, and residual-leverage plots for a linear model.

Four standard diagnostics for the working model. Look for systematic shape, tail departures, a funnel pattern, and influential observations. 

par(old_par)
  • Residuals vs Fitted: an unstructured cloud around zero is desirable; curvature suggests a functional-form problem.
  • Normal Q–Q: emphasize systematic and tail departures; minor departures may be unimportant in a large sample.
  • Scale–Location: an upward trend or funnel suggests changing conditional variance.
  • Cook’s distance: identifies observations with substantial influence on the fit; it is not an automatic deletion rule.

8.4 Quantitative screening for heteroskedasticity

First divide fitted values into five groups and compare residual standard deviations:

fitted_groups <- cut(
  fitted(working_model),
  breaks = quantile(fitted(working_model), probs = seq(0, 1, 0.2)),
  include.lowest = TRUE
)
variance_by_fit <- data.frame(
  "Fitted-value group" = levels(fitted_groups),
  "Sample size" = as.vector(table(fitted_groups)),
  "Residual standard deviation" = as.vector(
    tapply(residuals(working_model), fitted_groups, sd)
  ),
  check.names = FALSE
)
knitr::kable(
  variance_by_fit,
  digits = 2,
  caption = "Residual variability across fitted-value quintiles"
)
Residual variability across fitted-value quintiles
Fitted-value group Sample size Residual standard deviation
[107,119] 84 10.65
(119,124] 84 9.55
(124,130] 84 9.55
(130,137] 84 9.63
(137,166] 84 11.73

The following one-predictor Breusch–Pagan-style screen regresses squared residuals on fitted values. Its statistic is compared approximately with a chi-squared distribution with one degree of freedom.

hetero_auxiliary <- lm(I(residuals(working_model)^2) ~ fitted(working_model))
hetero_statistic <- nobs(working_model) * summary(hetero_auxiliary)$r.squared
hetero_p_value <- pchisq(hetero_statistic, df = 1, lower.tail = FALSE)
c(statistic = hetero_statistic, df = 1, p_value = hetero_p_value)
## statistic        df   p_value 
##     1.319     1.000     0.251

This simplified screen is not a final verdict. A significant result may arise from a misspecified mean; a non-significant result does not prove constant variance. Combine it with residual plots, knowledge of the data-generating process, and robust sensitivity analyses.

8.5 Leverage, studentized residuals, and Cook’s distance

influence_summary <- data.frame(
  "Participant ID" = ph_data$participant_id,
  "Leverage" = hatvalues(working_model),
  "Studentized residual" = rstudent(working_model),
  "Cook's distance" = cooks.distance(working_model),
  check.names = FALSE
)
top_influence <- order(
  influence_summary[["Cook's distance"]],
  decreasing = TRUE
)[1:8]
knitr::kable(
  influence_summary[top_influence, ],
  digits = 3,
  caption = "Eight observations with the largest Cook's distances"
)
Eight observations with the largest Cook’s distances
Participant ID Leverage Studentized residual Cook’s distance
486 P486 0.108 -1.91 0.044
399 P399 0.035 -3.01 0.032
235 P235 0.026 -2.91 0.022
543 P543 0.037 -2.35 0.021
71 P071 0.034 -2.40 0.020
375 P375 0.135 -1.10 0.019
367 P367 0.056 1.76 0.018
213 P213 0.035 -2.23 0.018

Common heuristics include leverage above 2p/n2p/n, Cook’s distance above 4/n4/n, and an absolute studentized residual above 2 or 3. These are not universal deletion thresholds.

p_working <- ncol(model.matrix(working_model))
n_working <- nobs(working_model)
c(
  leverage_2p_over_n = 2 * p_working / n_working,
  cooks_4_over_n = 4 / n_working,
  count_abs_studentized_over_2 = sum(abs(rstudent(working_model)) > 2)
)
##           leverage_2p_over_n               cooks_4_over_n 
##                      0.04762                      0.00952 
## count_abs_studentized_over_2 
##                     18.00000

For a high-influence observation, verify the source record, units, and eligibility, then compare analyses with and without it. Exclusion requires a defensible data-quality or target-population reason; making a p-value smaller is not one.

8.6 Multicollinearity

Collinearity does not necessarily bias predictions, but it can inflate coefficient standard errors and make coefficients sensitive to small data changes. The function below calculates a VIF for each design-matrix column:

vif_columns <- function(model) {
  X <- model.matrix(model)
  X <- X[, colnames(X) != "(Intercept)", drop = FALSE]
  values <- vapply(seq_len(ncol(X)), function(j) {
    outcome_column <- X[, j]
    other_columns <- X[, -j, drop = FALSE]
    r_squared_j <- summary(lm(outcome_column ~ other_columns))$r.squared
    if (1 - r_squared_j < .Machine$double.eps^0.5) {
      Inf
    } else {
      1 / (1 - r_squared_j)
    }
  }, numeric(1))
  data.frame(
    "Design-matrix column" = label_terms(colnames(X)),
    "VIF" = values,
    row.names = NULL,
    check.names = FALSE
  )
}

knitr::kable(
  vif_columns(working_model),
  digits = 2,
  caption = "Column-level VIFs for the working model design matrix"
)
Column-level VIFs for the working model design matrix
Design-matrix column VIF
Age (per 10 years; centered at 50) 1.35
Current smoking (reference: not current) 1.02
Age squared 1.17
BMI (per 5 units; centered at 25) 1.03
Male (reference: female) 1.05
North (reference: Central) 1.33
South (reference: Central) 1.35
West (reference: Central) 1.33
Age × current smoking 1.17

Treat VIF cutoffs such as 5 or 10 as heuristics. Column-level VIFs need particular care with multilevel factors, interactions, and polynomials. Centering can reduce unnecessary numerical collinearity, but it cannot repair two distinct variables that nearly measure the same construct.

9 HC3 Heteroskedasticity-Robust Standard Errors

9.1 Why use robust standard errors?

If the conditional-mean model is reasonable but error variance is not constant, OLS coefficients can still describe conditional-mean associations while classical standard errors may be unreliable. A heteroskedasticity-consistent covariance matrix replaces the common-variance assumption with observation-specific residual information.

HC3 strongly adjusts high-leverage observations. Its “meat” matrix uses:

Ω̂ii=ei2(1−hii)2,V̂HC3=(XTX)−1XTΩ̂X(XTX)−1. \widehat\Omega_{ii}=\frac{e_i^2}{(1-h_{ii})^2}, \qquad \widehat V_{HC3}=(X^TX)^{-1}X^T\widehat\Omega X(X^TX)^{-1}.

9.2 Implement HC3 using base R only

hc3_table <- function(model, level = 0.95) {
  if (anyNA(coef(model))) {
    stop("The model has non-estimable coefficients; resolve perfect collinearity first.")
  }
  X <- model.matrix(model)
  residual <- residuals(model)
  leverage <- hatvalues(model)
  omega_hc3 <- residual^2 / (1 - leverage)^2

  bread <- solve(crossprod(X))
  meat <- crossprod(X, X * omega_hc3)
  covariance_hc3 <- bread %*% meat %*% bread
  standard_error <- sqrt(diag(covariance_hc3))

  estimate <- coef(model)
  statistic <- estimate / standard_error
  degrees_freedom <- df.residual(model)
  critical <- qt(1 - (1 - level) / 2, degrees_freedom)

  data.frame(
    term = names(estimate),
    estimate = unname(estimate),
    robust_se_hc3 = standard_error,
    lower = estimate - critical * standard_error,
    upper = estimate + critical * standard_error,
    p_value = 2 * pt(
      abs(statistic),
      df = degrees_freedom,
      lower.tail = FALSE
    ),
    row.names = NULL
  )
}
working_hc3 <- hc3_table(working_model)
classic_se <- summary(working_model)$coefficients[, "Std. Error"]
hc3_comparison <- data.frame(
  "Term" = label_terms(working_hc3$term),
  "Estimate" = working_hc3$estimate,
  "Classical SE" = classic_se[working_hc3$term],
  "HC3 SE" = working_hc3$robust_se_hc3,
  "HC3 lower" = working_hc3$lower,
  "HC3 upper" = working_hc3$upper,
  "HC3 p-value" = working_hc3$p_value,
  row.names = NULL,
  check.names = FALSE
)
knitr::kable(
  hc3_comparison,
  digits = 3,
  caption = "Classical and HC3 robust standard errors"
)
Classical and HC3 robust standard errors
Term Estimate Classical SE HC3 SE HC3 lower HC3 upper HC3 p-value
Intercept 125.623 1.152 1.150 123.363 127.883 0.000
Age (per 10 years; centered at 50) 5.916 0.372 0.384 5.160 6.672 0.000
Current smoking (reference: not current) 4.636 1.393 1.549 1.592 7.680 0.003
Age squared 0.871 0.180 0.188 0.500 1.241 0.000
BMI (per 5 units; centered at 25) 3.238 0.610 0.602 2.054 4.422 0.000
Male (reference: female) 2.610 1.033 1.051 0.544 4.676 0.013
North (reference: Central) 2.994 1.445 1.528 -0.010 5.998 0.051
South (reference: Central) -3.392 1.334 1.373 -6.091 -0.693 0.014
West (reference: Central) 1.762 1.422 1.377 -0.946 4.469 0.202
Age × current smoking 1.657 0.952 0.990 -0.290 3.603 0.095

Robust standard errors do not repair every problem HC3 changes uncertainty estimates, not OLS coefficients. It cannot fix a misspecified conditional mean, serious measurement error, uncontrolled confounding, selection bias, clustered dependence, or unsupported extrapolation. Clustered data require standard errors or models matched to that dependence structure.

10 Confidence Intervals and Prediction Intervals

10.1 Two intervals answer different questions

  • Confidence interval: How uncertain is the population mean systolic blood pressure at a specified covariate combination?
  • Prediction interval: Where might the outcome of one new person with that covariate combination fall?

A prediction interval also includes individual residual variation around the conditional mean, so it is usually wider.

new_people <- data.frame(
  age = c(35, 50, 70),
  bmi = c(23, 25, 30),
  smoking_status = factor(
    c("Not current", "Current", "Not current"),
    levels = levels(ph_data$smoking_status)
  ),
  sex = factor(
    c("Female", "Male", "Female"),
    levels = levels(ph_data$sex)
  ),
  neighborhood = factor(
    c("Central", "North", "West"),
    levels = levels(ph_data$neighborhood)
  )
)
new_people$age_c10 <- (new_people$age - 50) / 10
new_people$bmi_c5 <- (new_people$bmi - 25) / 5

mean_intervals <- predict(
  working_model, newdata = new_people,
  interval = "confidence", level = 0.95
)
individual_intervals <- predict(
  working_model, newdata = new_people,
  interval = "prediction", level = 0.95
)
interval_table <- cbind(
  new_people[c("age", "bmi", "smoking_status", "sex", "neighborhood")],
  "Mean fitted value" = mean_intervals[, "fit"],
  "Mean CI lower" = mean_intervals[, "lwr"],
  "Mean CI upper" = mean_intervals[, "upr"],
  "Individual PI lower" = individual_intervals[, "lwr"],
  "Individual PI upper" = individual_intervals[, "upr"]
)
knitr::kable(
  interval_table,
  digits = 1,
  caption = "Confidence intervals for means and prediction intervals for individuals"
)
Confidence intervals for means and prediction intervals for individuals
age bmi smoking_status sex neighborhood Mean fitted value Mean CI lower Mean CI upper Individual PI lower Individual PI upper
35 23 Not current Female Central 117 115 120 97 138
50 25 Current Male North 136 132 139 115 156
70 30 Not current Female West 146 142 149 125 166

predict.lm(..., interval = "prediction") uses the classical homoskedastic model. Under substantial heteroskedasticity, individual variability may change across predictor values, so classical prediction intervals may be miscalibrated. Robust coefficient standard errors alone do not automatically yield valid individual prediction intervals.

simple_grid <- data.frame(
  age = seq(min(ph_data$age), max(ph_data$age), length.out = 160)
)
simple_confidence <- predict(simple_model, simple_grid, interval = "confidence")
simple_prediction <- predict(simple_model, simple_grid, interval = "prediction")
plot(
  ph_data$age, ph_data$systolic_bp,
  pch = 16, cex = 0.42, col = "grey75",
  xlab = "Age (years)", ylab = "Systolic blood pressure (mmHg)",
  main = "A confidence band is not a prediction band"
)
polygon(
  c(simple_grid$age, rev(simple_grid$age)),
  c(simple_prediction[, "lwr"], rev(simple_prediction[, "upr"])),
  col = rgb(230 / 255, 159 / 255, 0, 0.16), border = NA
)
polygon(
  c(simple_grid$age, rev(simple_grid$age)),
  c(simple_confidence[, "lwr"], rev(simple_confidence[, "upr"])),
  col = rgb(0, 114 / 255, 178 / 255, 0.28), border = NA
)
lines(
  simple_grid$age, simple_confidence[, "fit"],
  col = palette_ph["navy"], lwd = 2.5
)
legend(
  "topleft",
  legend = c("Fitted mean", "Mean 95% CI", "Individual 95% PI"),
  col = c(palette_ph["navy"], palette_ph["blue"], palette_ph["orange"]),
  lwd = c(2.5, 8, 8), bty = "n"
)
Age and blood-pressure points with a narrow mean confidence band and a wider individual prediction band.

In simple regression, the 95% confidence band for the mean response and the 95% prediction band for a new individual.

11 Out-of-Sample Predictive Performance

11.1 Why separate training and test data?

Adding variables cannot increase the residual sum of squares in training data, yet performance on new data may deteriorate. When prediction is a goal, set aside test data before modeling, and perform all variable selection, transformations, and tuning within the training data.

train_data <- ph_data

prediction_formula <- systolic_bp ~
  age_c10 * smoking_status + I(age_c10^2) +
  bmi_c5 + sex + neighborhood

train_model <- lm(prediction_formula, data = train_data)
test_prediction <- predict(train_model, newdata = test_data)

This is the first point at which the untouched holdout outcome is opened and summarized.

split_composition <- data.frame(
  "Dataset" = c("Development/training", "Untouched holdout test"),
  "Sample size" = c(nrow(train_data), nrow(test_data)),
  "Mean systolic BP" = c(
    mean(train_data$systolic_bp),
    mean(test_data$systolic_bp)
  ),
  check.names = FALSE
)
knitr::kable(
  split_composition,
  digits = 1,
  caption = "Prespecified training data and the newly opened holdout test data"
)
Prespecified training data and the newly opened holdout test data
Dataset Sample size Mean systolic BP
Development/training 420 128
Untouched holdout test 140 131 
prediction_metrics <- function(observed, predicted) {
  c(
    RMSE = sqrt(mean((observed - predicted)^2)),
    MAE = mean(abs(observed - predicted)),
    R2_test = 1 - sum((observed - predicted)^2) /
      sum((observed - mean(observed))^2)
  )
}

model_metrics <- prediction_metrics(test_data$systolic_bp, test_prediction)
baseline_prediction <- rep(mean(train_data$systolic_bp), nrow(test_data))
baseline_metrics <- prediction_metrics(test_data$systolic_bp, baseline_prediction)
performance_table <- rbind(
  "Regression model" = model_metrics,
  "Training-mean baseline" = baseline_metrics
)
knitr::kable(
  performance_table,
  digits = 2,
  caption = "Predictive performance in the untouched holdout test data"
)
Predictive performance in the untouched holdout test data
RMSE MAE R2_test
Regression model 10.7 8.73 0.51
Training-mean baseline 15.7 11.87 -0.04

RMSE is more sensitive to large errors, whereas MAE averages absolute errors. Test-set R2R^2 can be negative, meaning the model performs worse than a simple mean-based reference. A single random split has sampling uncertainty; a formal prediction study generally needs cross-validation and external validation.

Prevent data leakage If the analyst uses the full dataset to select variables, handle unusual observations, or choose transformations before splitting, the test data have already influenced training. Estimate imputation, standardization, and feature-selection rules within each training fold and then apply those rules to validation data.

12 Missing Data and Causal Interpretation

12.1 Complete-case analysis changes who is analyzed

By default, lm() drops a row when any variable in the formula is missing. Adding physical activity therefore reduces the analysis sample:

activity_model <- lm(
  systolic_bp ~ age_c10 * smoking_status + I(age_c10^2) +
    bmi_c5 + sex + neighborhood + physical_activity_30,
  data = ph_data,
  na.action = na.exclude
)

missing_model_comparison <- rbind(
  fit_statistics(working_model, "Without physical activity"),
  fit_statistics(activity_model, "Physical activity added; complete cases")
)
knitr::kable(
  missing_model_comparison,
  digits = 3,
  caption = "Change in the analysis sample after adding a variable with missingness"
)
Change in the analysis sample after adding a variable with missingness
Model Sample size Parameters R-squared Adjusted R-squared Residual standard error
Without physical activity 420 10 0.495 0.484 10.3
Physical activity added; complete cases 391 11 0.497 0.483 10.4

Compare coefficients from the two models cautiously: a difference may reflect adjustment for physical activity, a change in the analyzed sample, or both. Report variable-level missingness, the model sample size, assumptions about missingness, and the handling method.

Complete-case analysis is easiest to justify under strong conditions such as missing completely at random. When missingness relates to observed information, multiple imputation may be appropriate; dependence on unobserved values also calls for sensitivity analysis. This tutorial does not use single mean imputation because it understates variability and distorts relationships.

12.2 More adjustment is not always better

Choose variables from the estimand and causal structure:

  • Confounders may need control.
  • Adjusting for a mediator that occurs after exposure changes the effect being estimated.
  • Adjusting for a collider affected by causes of both exposure and outcome can introduce bias.
  • Automated stepwise selection yields unstable models, exaggerates significance, and ignores subject-matter knowledge.
  • Ordinary lm() independence-based inference is inappropriate for complex samples, clusters, or repeated measures unless the design is handled explicitly.

Causal language requires additional support “The association remained significant after controlling for age and BMI” does not establish causation. A causal interpretation requires clear temporality, a well-defined intervention, confounding assumptions, comparability, justified missingness and selection mechanisms, and an estimator appropriate to the study design.

13 Reporting Linear Regression

13.1 Minimum reporting checklist

A reproducible report should state:

  1. the target population, analysis sample, outcome, primary exposure, and estimand;
  2. numeric-variable units, centering, transformations, and categorical reference groups;
  3. the rationale for covariates and whether interactions or nonlinear terms were prespecified;
  4. the analysis sample size, key-variable missingness, and missing-data handling;
  5. coefficients and 95% CIs, with p-values when useful but never as substitutes for effect interpretation;
  6. R2R^2, adjusted R2R^2, or an out-of-sample measure aligned with the objective;
  7. checks of residuals, heteroskedasticity, influence, collinearity, and independence;
  8. whether standard errors are classical or robust, and why;
  9. limitations involving design, measurement, selection, residual confounding, and transportability; and
  10. software, code, random seeds, and data provenance sufficient for reproduction.

13.2 A four-sentence results example

Among 420 simulated participants, we used OLS to describe conditional associations of systolic blood pressure with age, BMI, smoking status, sex, and neighborhood, allowing a prespecified quadratic age term and age-by-smoking interaction. Among non-current smokers at the same modeled values of other covariates, a 10-year age difference near age 50 corresponded to 5.9 mmHg in mean systolic blood pressure (HC3 95% CI 5.2 to 6.7). The in-sample adjusted R2R^2 was 0.48; residual checks suggested non-constant variance, so we report HC3 robust standard errors. These simulated cross-sectional-style data do not support a causal effect of aging, and the result should not be extrapolated beyond the observed covariate ranges.

A full report should also interpret the interaction explicitly, noting that the displayed age coefficient belongs to the non-current-smoking reference group and graphing predictions for both groups.

14 Exercises and Answers

14.1 Exercise 1: Identify a suitable outcome

Which question is most directly suited to ordinary linear regression?

  1. Whether a person is hospitalized in the next 30 days.
  2. The number of hospitalizations in the next 30 days.
  3. Systolic blood pressure in mmHg at discharge.
  4. Time from enrollment to first hospitalization.
Show answer Answer: 3. Systolic blood pressure is continuous, so linear regression can describe its conditional mean. Binary hospitalization, counts, and event times generally need methods matched to their distributions and censoring structures.

14.2 Exercise 2: Interpret a numeric coefficient

How should the bmi_c5 coefficient in scaled_model be interpreted? Can it be called a causal effect of BMI?

Show answer It is the mean systolic blood-pressure difference for a 5-unit BMI difference among participants at the same modeled age, smoking status, sex, and neighborhood. This observational regression alone does not establish causation because confounding, selection, measurement, and temporality have not been resolved.

14.3 Exercise 3: Interpret a categorical reference group

What comparison does neighborhoodSouth represent in multiple_model? Would model predictions change if South became the reference?

Show answer It compares South with the first factor level, Central, holding the other modeled variables fixed. Changing the reference to South re-expresses the intercept and neighborhood coefficients, but fitted values, residuals, R2R^2, and the overall neighborhood information remain unchanged within the same model space.

14.4 Exercise 4: Fit a model

Fit systolic blood pressure as a function of age, BMI, and smoking status, and report 95% CIs. The answer uses only base R and the helper defined above.

Show answer and code 
exercise_model <- lm(
  systolic_bp ~ age_c10 + bmi_c5 + smoking_status,
  data = ph_data
)
knitr::kable(
  coefficient_table(exercise_model),
  digits = 3,
  caption = "Exercise model coefficients and 95% confidence intervals"
)
Exercise model coefficients and 95% confidence intervals
Term Estimate Standard error 95% CI lower 95% CI upper p-value
Intercept 129.24 0.595 128.07 130.41 0.000
Age (per 10 years; centered at 50) 5.39 0.340 4.72 6.06 0.000
BMI (per 5 units; centered at 25) 3.13 0.637 1.87 4.38 0.000
Current smoking (reference: not current) 4.17 1.458 1.30 7.03 0.004
Age and BMI are expressed per 10 years and 5 units, respectively. The smoking coefficient uses non-current smoking as the reference.

14.5 Exercise 5: Read an interaction

Suppose the age main effect is 5.0 mmHg per 10 years and the age-by-current-smoking interaction is 1.8 mmHg per 10 years. What is the age slope among current smokers? At what age does the smoking_statusCurrent main effect compare groups?

Show answer The age slope among current smokers is 5.0+1.8=6.85.0+1.8=6.8 mmHg per 10 years. Because age is centered at 50, the smoking main effect compares groups at age 50 and the same specified values of other covariates.

14.6 Exercise 6: Diagnose a residual funnel

What would you do if the Residuals vs Fitted plot showed residual spread increasing with fitted values?

Show answer First verify outcome units, unusual records, and functional form. Then examine grouped residual variances and ask whether the changing variability is scientifically plausible. If the conditional-mean model is reasonable, report an HC3-style robust sensitivity analysis. Individual prediction additionally requires modeling or validating the changing variance. Robust standard errors alone do not fix omitted nonlinearity, clustering, or bias.

14.7 Exercise 7: Confidence or prediction interval?

A health department wants uncertainty for the mean systolic blood pressure of 60-year-old, BMI-25, non-current-smoking women. A clinician wants a plausible range for the next individual with those characteristics. Which interval serves each question?

Show answer The first question uses a confidence interval for the conditional mean; the second uses a prediction interval for a new individual. The prediction interval includes individual residual variability and is therefore usually much wider.

14.8 Exercise 8: What does a high training R² mean?

A model with many variables and interactions has training R2=0.90R^2=0.90. Is that enough to claim excellent predictive performance?

Show answer No. Training R2R^2 cannot decrease as model complexity increases and may reflect overfitting. Evaluate RMSE, MAE, calibration, and the intended population in data that played no role in variable selection or tuning, with external validation when possible.

15 Quick Reference

15.1 Common base R code

Purpose Code pattern
Fit a linear model lm(y ~ x1 + x2, data = d)
View a full summary summary(model)
Extract coefficients coef(model)
Classical covariance matrix vcov(model)
95% coefficient intervals confint(model)
Fitted values and residuals fitted(model); residuals(model)
Standard diagnostic plots plot(model)
Leverage and Cook’s distance hatvalues(model); cooks.distance(model)
Studentized residuals rstudent(model)
Confidence interval for a mean predict(model, newdata, interval = "confidence")
Prediction interval for an individual predict(model, newdata, interval = "prediction")
Compare nested models anova(smaller, larger)
Inspect the design matrix model.matrix(model)

15.2 Interpretation sequence

When reading model output, proceed in this order:

Question and estimand → analysis sample → units and reference groups → functional form → coefficients and intervals → interaction or nonlinearity → diagnostics → out-of-sample performance → limitations and domain of use

15.3 Core formulas

Concept Formula or meaning
Conditional mean E(Y∣X)=XβE(Y\mid X)=X\beta
Residual ei=Yi−Ŷie_i=Y_i-\hat Y_i
OLS objective Minimize ∑iei2\sum_i e_i^2
OLS solution β̂=(XTX)−1XTY\hat\beta=(X^TX)^{-1}X^TY
R2R^2 1−SSE/SST1-\mathrm{SSE}/\mathrm{SST}
RMSE n−1∑i(Yi−Ŷi)2\sqrt{n^{-1}\sum_i(Y_i-\hat Y_i)^2}
HC3 weight ei2/(1−hii)2e_i^2/(1-h_{ii})^2
Age slope among current smokers Age main effect + age × smoking interaction

Final model check

Before publishing results, confirm that:

  • the outcome is a continuous variable suitable for conditional-mean modeling;
  • units, coding, eligibility criteria, and reference groups have been verified;
  • functional form and interactions are supported by the question, graphs, and subject-matter knowledge;
  • every reported coefficient can be explained in a complete sentence;
  • the interval type matches the question;
  • residuals, heteroskedasticity, influence, collinearity, and dependence have been considered;
  • missingness and changes in the analysis sample are reported;
  • in-sample fit has not been mislabeled as out-of-sample predictive ability;
  • an adjusted association is not called causal without sufficient support; and
  • code, seeds, and software information are adequate for reproduction.
sessionInfo()
## R version 4.6.1 (2026-06-24)
## Platform: aarch64-apple-darwin23
## Running under: macOS Tahoe 26.5.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
## 
## locale:
## [1] C.UTF-8/C.UTF-8/C.UTF-8/C/C.UTF-8/C.UTF-8
## 
## time zone: America/Edmonton
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## loaded via a namespace (and not attached):
##  [1] digest_0.6.39   R6_2.6.1        fastmap_1.2.0   xfun_0.60
##  [5] cachem_1.1.0    knitr_1.51      htmltools_0.5.9 rmarkdown_2.31
##  [9] lifecycle_1.0.5 cli_3.6.6       sass_0.4.10     jquerylib_0.1.4
## [13] compiler_4.6.1  tools_4.6.1     evaluate_1.0.5  bslib_0.12.0
## [17] yaml_2.3.12     rlang_1.3.0     jsonlite_2.0.0

All data and results in this tutorial are simulated for statistical education. Real public health decisions require attention to study design, data quality, subject-matter evidence, and appropriate professional review.