Code
library(dplyr)
library(ggplot2)
library(broom) # tidy model output
library(car) # vif(), crPlots(), influence tools
library(lmtest) # bptest()
library(knitr)
theme_set(theme_minimal(base_size = 12))Linear Regression in R: Fitting, Assumptions and Interpretation
A worked example with birth weight data
1 Purpose and scenario
This document walks through a complete linear regression workflow: exploring the data, fitting simple and multiple models, checking every assumption, and interpreting the output in plain language.
Scenario. A cohort of 500 singleton births from antenatal clinics. The research question is:
Which maternal and pregnancy factors are associated with infant birth weight, and by how much?
| Variable | Type | Description |
|---|---|---|
id |
identifier | Record number |
birth_weight_g |
numeric (outcome) | Birth weight in grams |
gestational_age_wks |
numeric | Gestational age at birth, weeks |
maternal_age_yrs |
numeric | Maternal age, years |
maternal_bmi |
numeric | Maternal BMI, kg/m² |
smoker |
binary | Smoked during pregnancy (no / yes) |
anc_visits |
count | Number of antenatal care visits |
residence |
binary | urban / rural |
infant_sex |
binary | female / male |
2 Setup
3 Load and inspect the data
Code
bw <- read.csv("birthweight.csv", stringsAsFactors = FALSE)
# Declare categorical variables as factors and set the reference level explicitly.
bw <- bw |>
mutate(
smoker = factor(smoker, levels = c("no", "yes")),
residence = factor(residence, levels = c("urban", "rural")),
infant_sex = factor(infant_sex, levels = c("female", "male"))
)
str(bw)'data.frame': 500 obs. of 9 variables:
$ id : int 1 2 3 4 5 6 7 8 9 10 ...
$ birth_weight_g : int 3272 3848 2555 3329 3511 3370 3324 3482 3284 3285 ...
$ gestational_age_wks: num 37.4 39.2 35.4 41.3 39.9 38.3 38.2 39.3 38.3 38.4 ...
$ maternal_age_yrs : int 38 31 21 21 22 30 34 16 30 32 ...
$ maternal_bmi : num 26 27.7 26.1 27.1 20 26.2 24.1 24.2 29.4 22.6 ...
$ smoker : Factor w/ 2 levels "no","yes": 1 1 1 1 2 1 1 1 1 1 ...
$ anc_visits : int 2 6 3 1 6 4 0 1 7 1 ...
$ residence : Factor w/ 2 levels "urban","rural": 1 1 2 2 1 2 2 2 1 2 ...
$ infant_sex : Factor w/ 2 levels "female","male": 2 1 1 2 1 2 2 2 1 1 ...Code
summary(bw) id birth_weight_g gestational_age_wks maternal_age_yrs
Min. : 1.0 Min. :1849 Min. :32.10 Min. :15.00
1st Qu.:125.8 1st Qu.:2939 1st Qu.:37.58 1st Qu.:23.75
Median :250.5 Median :3206 Median :38.90 Median :27.00
Mean :250.5 Mean :3189 Mean :38.82 Mean :27.09
3rd Qu.:375.2 3rd Qu.:3434 3rd Qu.:40.10 3rd Qu.:31.00
Max. :500.0 Max. :4293 Max. :42.00 Max. :44.00
maternal_bmi smoker anc_visits residence infant_sex
Min. :16.00 no :445 Min. : 0.000 urban:225 female:273
1st Qu.:21.60 yes: 55 1st Qu.: 3.000 rural:275 male :227
Median :24.20 Median : 4.000
Mean :24.15 Mean : 4.226
3rd Qu.:26.62 3rd Qu.: 6.000
Max. :37.20 Max. :10.000 Code
colSums(is.na(bw)) id birth_weight_g gestational_age_wks maternal_age_yrs
0 0 0 0
maternal_bmi smoker anc_visits residence
0 0 0 0
infant_sex
0 Always check types, ranges and missingness first. The first level of each factor becomes the reference category, so no, urban and female are the comparison groups.
4 Exploratory analysis
Regression assumptions are about the residuals, but looking at the raw data first tells you what to expect.
Code
ggplot(bw, aes(birth_weight_g)) +
geom_histogram(bins = 30, fill = "steelblue", colour = "white") +
labs(x = "Birth weight (g)", y = "Count")
Code
ggplot(bw, aes(gestational_age_wks, birth_weight_g, colour = smoker)) +
geom_point(alpha = 0.5) +
geom_smooth(method = "lm", se = FALSE) +
labs(x = "Gestational age (weeks)", y = "Birth weight (g)", colour = "Smoker")
Code
bw |>
select(birth_weight_g, gestational_age_wks, maternal_age_yrs,
maternal_bmi, anc_visits) |>
cor() |>
round(2) birth_weight_g gestational_age_wks maternal_age_yrs
birth_weight_g 1.00 0.60 0.06
gestational_age_wks 0.60 1.00 0.01
maternal_age_yrs 0.06 0.01 1.00
maternal_bmi 0.10 0.01 0.06
anc_visits 0.17 0.06 0.05
maternal_bmi anc_visits
birth_weight_g 0.10 0.17
gestational_age_wks 0.01 0.06
maternal_age_yrs 0.06 0.05
maternal_bmi 1.00 0.00
anc_visits 0.00 1.00Things to look for: an approximately linear trend with gestational age, roughly parallel lines for smokers and non-smokers, and no pair of predictors that is very strongly correlated with another (which would cause trouble later).
5 Simple linear regression
The model with one predictor:
\[ \text{BirthWeight}_i = \beta_0 + \beta_1 \,\text{GestAge}_i + \varepsilon_i, \qquad \varepsilon_i \sim N(0,\sigma^2) \]
Code
m_simple <- lm(birth_weight_g ~ gestational_age_wks, data = bw)
summary(m_simple)
Call:
lm(formula = birth_weight_g ~ gestational_age_wks, data = bw)
Residuals:
Min 1Q Median 3Q Max
-802.9 -223.8 3.3 223.7 916.5
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2076.898 314.716 -6.599 1.06e-10 ***
gestational_age_wks 135.657 8.099 16.750 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 325.3 on 498 degrees of freedom
Multiple R-squared: 0.3603, Adjusted R-squared: 0.3591
F-statistic: 280.5 on 1 and 498 DF, p-value: < 2.2e-16Code
confint(m_simple) 2.5 % 97.5 %
(Intercept) -2695.2320 -1458.5635
gestational_age_wks 119.7441 151.5694Interpretation.
- Slope (\(\hat\beta_1\) = 135.7): on average, each additional week of gestation is associated with a birth weight that is 136 g higher (95% CI 120 to 152).
- Intercept (\(\hat\beta_0\) = -2077): the predicted birth weight at 0 weeks of gestation. This is an extrapolation far outside the data and has no biological meaning. Centering the predictor fixes this (see below).
- \(R^2\) = 0.36: gestational age alone explains about 36% of the variation in birth weight.
- Residual standard error = 325 g: a typical observation falls about this far from the fitted line.
5.0.1 Centering to get a meaningful intercept
Code
m_centered <- lm(birth_weight_g ~ I(gestational_age_wks - 39), data = bw)
coef(m_centered) (Intercept) I(gestational_age_wks - 39)
3213.7154 135.6567 Now the intercept is the predicted birth weight at 39 weeks (3214 g). The slope is unchanged. Centering never changes the slope; it only changes what the intercept means.
6 Multiple linear regression
The research question involves several factors, and gestational age is a confounder for most of them, so we adjust for everything simultaneously:
\[ \begin{aligned} \text{BW}_i = \beta_0 &+ \beta_1\,\text{GA}_i + \beta_2\,\text{Age}_i + \beta_3\,\text{BMI}_i + \beta_4\,\text{Smoker}_i \\ &+ \beta_5\,\text{ANC}_i + \beta_6\,\text{Rural}_i + \beta_7\,\text{Male}_i + \varepsilon_i \end{aligned} \]
Code
m_multi <- lm(
birth_weight_g ~ gestational_age_wks + maternal_age_yrs + maternal_bmi +
smoker + anc_visits + residence + infant_sex,
data = bw
)
summary(m_multi)
Call:
lm(formula = birth_weight_g ~ gestational_age_wks + maternal_age_yrs +
maternal_bmi + smoker + anc_visits + residence + infant_sex,
data = bw)
Residuals:
Min 1Q Median 3Q Max
-772.03 -200.10 1.92 205.73 833.77
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2193.166 316.631 -6.927 1.35e-11 ***
gestational_age_wks 129.699 7.680 16.888 < 2e-16 ***
maternal_age_yrs 1.792 2.514 0.713 0.476286
maternal_bmi 10.650 3.837 2.776 0.005715 **
smokeryes -256.314 44.147 -5.806 1.15e-08 ***
anc_visits 24.344 6.721 3.622 0.000323 ***
residencerural -100.564 27.937 -3.600 0.000351 ***
infant_sexmale 49.318 27.645 1.784 0.075050 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 305.7 on 492 degrees of freedom
Multiple R-squared: 0.4417, Adjusted R-squared: 0.4338
F-statistic: 55.61 on 7 and 492 DF, p-value: < 2.2e-16A tidier table with confidence intervals:
Code
tidy(m_multi, conf.int = TRUE) |>
mutate(across(c(estimate, std.error, statistic, conf.low, conf.high), ~ round(.x, 1)),
p.value = format.pval(p.value, digits = 3, eps = 0.001)) |>
kable(col.names = c("Term", "Estimate", "SE", "t", "p", "95% CI low", "95% CI high"))| Term | Estimate | SE | t | p | 95% CI low | 95% CI high |
|---|---|---|---|---|---|---|
| (Intercept) | -2193.2 | 316.6 | -6.9 | < 0.001 | -2815.3 | -1571.1 |
| gestational_age_wks | 129.7 | 7.7 | 16.9 | < 0.001 | 114.6 | 144.8 |
| maternal_age_yrs | 1.8 | 2.5 | 0.7 | 0.47629 | -3.1 | 6.7 |
| maternal_bmi | 10.6 | 3.8 | 2.8 | 0.00571 | 3.1 | 18.2 |
| smokeryes | -256.3 | 44.1 | -5.8 | < 0.001 | -343.1 | -169.6 |
| anc_visits | 24.3 | 6.7 | 3.6 | < 0.001 | 11.1 | 37.6 |
| residencerural | -100.6 | 27.9 | -3.6 | < 0.001 | -155.5 | -45.7 |
| infant_sexmale | 49.3 | 27.6 | 1.8 | 0.07505 | -5.0 | 103.6 |
6.1 How to read the coefficients
Each coefficient is the expected change in birth weight for a one-unit change in that predictor, holding all other predictors in the model constant.
Code
b <- coef(m_multi)
ci <- confint(m_multi)| Term | Interpretation |
|---|---|
gestational_age_wks |
Each extra week of gestation: +130 g (95% CI 115 to 145), adjusted for the other variables. |
maternal_bmi |
Each extra BMI unit: +10.6 g (95% CI 3.1 to 18.2). |
smokeryes |
Smokers’ babies are on average -256 g lighter than non-smokers’ (95% CI -343 to -170). A categorical coefficient is a difference from the reference group. |
anc_visits |
Each additional antenatal visit: +24.3 g (95% CI 11.1 to 37.6). |
residencerural |
Rural births are -101 g relative to urban (95% CI -155 to -46). |
infant_sexmale |
Male infants: 49 g relative to female; the 95% CI (-5 to 104) is close to or includes zero, so the evidence is weak. |
maternal_age_yrs |
Each extra year of age: 1.8 g; the CI (-3.1 to 6.7) comfortably includes zero, so there is no clear association once the other variables are accounted for. |
(Intercept) |
Predicted weight for a female, urban, non-smoking infant with all numeric predictors equal to 0. Not meaningful on its own. |
Warning
Common interpretation errors.
- Saying “causes” when the design only supports “is associated with”.
- Forgetting “holding other variables constant”. A coefficient in a multiple model is not the same as the coefficient from a simple model.
- Treating a non-significant coefficient as proof of no effect. It means the data cannot distinguish the effect from zero; look at the confidence interval.
- Comparing coefficients across variables measured in different units. A BMI unit and a gestational week are not comparable. Standardise the predictors first if you need to compare them.
6.2 Overall model fit
Code
glance(m_multi) |>
select(r.squared, adj.r.squared, sigma, statistic, p.value, df, df.residual, AIC, BIC) |>
mutate(across(everything(), ~ round(.x, 3))) |>
kable()| r.squared | adj.r.squared | sigma | statistic | p.value | df | df.residual | AIC | BIC |
|---|---|---|---|---|---|---|---|---|
| 0.442 | 0.434 | 305.71 | 55.614 | 0 | 7 | 492 | 7151.512 | 7189.443 |
- \(R^2\) = 0.442: the model explains about 44% of the variability in birth weight. Adjusted \(R^2\) (0.434) penalises for the number of predictors and is the better figure for comparing models of different sizes.
- Residual SE = 306 g: this is the typical prediction error in the original units, often more useful than \(R^2\).
- F-test: tests whether at least one slope differs from zero. It is almost always significant when the individual predictors are, and tells you little beyond that.
6.3 Testing a group of terms
To test whether a set of predictors adds anything, compare nested models with an F-test. Here, do maternal age and infant sex jointly contribute?
Code
m_reduced <- update(m_multi, . ~ . - maternal_age_yrs - infant_sex)
anova(m_reduced, m_multi)Analysis of Variance Table
Model 1: birth_weight_g ~ gestational_age_wks + maternal_bmi + smoker +
anc_visits + residence
Model 2: birth_weight_g ~ gestational_age_wks + maternal_age_yrs + maternal_bmi +
smoker + anc_visits + residence + infant_sex
Res.Df RSS Df Sum of Sq F Pr(>F)
1 494 46332968
2 492 45981716 2 351252 1.8792 0.1538A large p-value means the two removed variables do not jointly improve the fit, and the simpler model is adequate. Whether to drop them is a scientific decision as well as a statistical one: if age and sex are known confounders, keep them regardless.
7 Checking the assumptions
The inference from a linear model (p-values, confidence intervals) relies on the LINE assumptions plus a few practical checks:
| Letter | Assumption | Main check |
|---|---|---|
| L | Linearity of the mean | Residuals vs fitted; component + residual plots |
| I | Independence of errors | Study design; residuals vs order/cluster |
| N | Normality of residuals | Q-Q plot, histogram |
| E | Equal variance (homoscedasticity) | Scale-location plot; Breusch-Pagan test |
| No severe multicollinearity | Variance inflation factors | |
| No unduly influential points | Cook’s distance, leverage |
First, collect the diagnostic quantities once:
Code
diag_df <- augment(m_multi) |>
mutate(obs = row_number())7.1 Linearity
Code
ggplot(diag_df, aes(.fitted, .resid)) +
geom_point(alpha = 0.5) +
geom_hline(yintercept = 0, linetype = "dashed") +
geom_smooth(method = "loess", se = FALSE, colour = "firebrick") +
labs(x = "Fitted values", y = "Residuals")
The red smoother should stay close to zero with no curve. A U-shape or arch means a predictor has a non-linear effect.
Code
crPlots(m_multi, terms = ~ gestational_age_wks + maternal_age_yrs +
maternal_bmi + anc_visits)
If the smoother departs clearly from the dashed line for a predictor, consider a polynomial or spline term for it.
7.2 Independence
This assumption cannot be fully tested from the residuals; it comes from the study design. Here every row is a different mother-infant pair, so independence is reasonable. It would be violated by:
- repeated measures on the same person,
- multiple infants per mother (twins, successive pregnancies),
- clustering such as women nested within health facilities or villages.
In those cases use mixed-effects models (lme4::lmer) or cluster-robust standard errors. A quick check for patterns related to record order:
Code
ggplot(diag_df, aes(obs, .resid)) +
geom_point(alpha = 0.5) +
geom_hline(yintercept = 0, linetype = "dashed") +
labs(x = "Record order", y = "Residuals")
7.3 Normality of residuals
Code
ggplot(diag_df, aes(sample = .std.resid)) +
stat_qq(alpha = 0.6) +
stat_qq_line(colour = "firebrick") +
labs(x = "Theoretical quantiles", y = "Standardised residuals")
Code
ggplot(diag_df, aes(.resid)) +
geom_histogram(aes(y = after_stat(density)), bins = 30,
fill = "grey70", colour = "white") +
stat_function(fun = dnorm,
args = list(mean = 0, sd = sd(diag_df$.resid)),
colour = "firebrick", linewidth = 1) +
labs(x = "Residuals", y = "Density")
Points should hug the line in the Q-Q plot. Mild departures in the tails are common and harmless.
Code
shapiro.test(residuals(m_multi))
Shapiro-Wilk normality test
data: residuals(m_multi)
W = 0.99675, p-value = 0.4161
Tip
Do not over-rely on formal normality tests. In large samples they flag trivial departures; in small samples they lack power. Use the plots as the primary judgement. Also note that normality matters mainly for small-sample inference and prediction intervals. With a reasonably large n, the coefficient estimates and their confidence intervals are robust to moderate non-normality (central limit theorem). Normality of the outcome itself is not required, only of the residuals.
7.4 Constant variance (homoscedasticity)
Code
ggplot(diag_df, aes(.fitted, sqrt(abs(.std.resid)))) +
geom_point(alpha = 0.5) +
geom_smooth(method = "loess", se = FALSE, colour = "firebrick") +
labs(x = "Fitted values", y = expression(sqrt("|Standardised residuals|")))
Code
bptest(m_multi)
studentized Breusch-Pagan test
data: m_multi
BP = 2.4104, df = 7, p-value = 0.9337The Breusch-Pagan test has the null hypothesis of constant variance. A small p-value suggests heteroscedasticity (a fan or funnel shape in the plot). Remedies: transform the outcome (log), use weighted least squares, or keep the model and report robust (heteroscedasticity-consistent) standard errors:
Code
library(sandwich)
coeftest(m_multi, vcov = vcovHC(m_multi, type = "HC3"))7.5 Multicollinearity
Code
vif(m_multi)gestational_age_wks maternal_age_yrs maternal_bmi smoker
1.017847 1.021874 1.005698 1.020795
anc_visits residence infant_sex
1.012900 1.033449 1.013540 A VIF measures how much the variance of a coefficient is inflated by correlation with the other predictors. Rough guide: values near 1 are ideal, above 5 deserves attention, above 10 is a serious problem. High VIFs inflate standard errors and make individual coefficients unstable, though they do not bias predictions from the model as a whole. Fixes include dropping or combining redundant predictors, or using ridge regression.
7.6 Influential observations and outliers
Code
cutoff <- 4 / nrow(bw)
ggplot(diag_df, aes(obs, .cooksd)) +
geom_col(width = 0.8, fill = "steelblue") +
geom_hline(yintercept = cutoff, linetype = "dashed", colour = "firebrick") +
geom_text(data = filter(diag_df, .cooksd > cutoff * 2.5),
aes(label = obs), vjust = -0.3, size = 3) +
labs(x = "Observation", y = "Cook's distance")
Code
diag_df |>
filter(.cooksd > cutoff) |>
arrange(desc(.cooksd)) |>
select(obs, birth_weight_g, gestational_age_wks, .fitted, .resid, .hat, .cooksd) |>
head(8) |>
mutate(across(where(is.numeric), ~ round(.x, 3))) |>
kable()| obs | birth_weight_g | gestational_age_wks | .fitted | .resid | .hat | .cooksd |
|---|---|---|---|---|---|---|
| 268 | 2107 | 36.5 | 2780.603 | -673.603 | 0.055 | 0.037 |
| 456 | 2080 | 35.9 | 2807.517 | -727.517 | 0.035 | 0.026 |
| 260 | 3315 | 36.2 | 2599.616 | 715.384 | 0.030 | 0.022 |
| 340 | 1849 | 33.6 | 2600.889 | -751.889 | 0.024 | 0.019 |
| 430 | 4163 | 42.0 | 3596.474 | 566.526 | 0.033 | 0.015 |
| 140 | 2962 | 34.4 | 2390.438 | 571.562 | 0.031 | 0.014 |
| 360 | 3248 | 37.1 | 2656.438 | 591.562 | 0.029 | 0.014 |
| 25 | 2668 | 40.0 | 3440.031 | -772.031 | 0.016 | 0.013 |
The line at \(4/n\) is a screening rule, not a verdict: some points will always exceed it. Investigate the flagged records. Are they data-entry errors? Biologically plausible extremes? Refit without them (a sensitivity analysis) and see if conclusions change. Do not delete observations merely because they are inconvenient.
Code
top_infl <- diag_df |> filter(.cooksd > cutoff * 2.5) |> pull(obs)
m_sens <- update(m_multi, data = bw[-top_infl, ])
cbind(full = coef(m_multi), without_influential = coef(m_sens)) |>
round(1) |>
kable()| full | without_influential | |
|---|---|---|
| (Intercept) | -2193.2 | -2162.2 |
| gestational_age_wks | 129.7 | 128.2 |
| maternal_age_yrs | 1.8 | 1.9 |
| maternal_bmi | 10.6 | 11.0 |
| smokeryes | -256.3 | -256.8 |
| anc_visits | 24.3 | 27.3 |
| residencerural | -100.6 | -91.2 |
| infant_sexmale | 49.3 | 52.9 |
If the coefficients barely move, the influential points are not driving your results.
7.7 Summary of the diagnostic checks
| Assumption | Verdict for this model |
|---|---|
| Linearity | No marked curvature expected in the plots. A quadratic term for BMI or gestational age may test as nominally significant (a chance finding in simulated data); judge by the size of the bend in the plot, not the p-value alone |
| Independence | Satisfied by design (one infant per mother, no clustering) |
| Normality | Q-Q plot close to the line |
| Constant variance | Scale-location plot flat; Breusch-Pagan p-value large |
| Multicollinearity | All VIFs close to 1 |
| Influence | No point shifts the coefficients materially |
8 What to do when assumptions fail
| Problem | Typical remedy |
|---|---|
| Non-linearity | Add polynomial or spline terms (poly(x, 2), splines::ns(x, 3)); transform the predictor |
| Heteroscedasticity | Log-transform the outcome; robust (HC3) SEs; weighted least squares |
| Non-normal residuals, skewed outcome | Transform the outcome; use a GLM (e.g. gamma family); bootstrap CIs |
| Correlated errors / clustering | Mixed-effects model; cluster-robust SEs; GEE |
| Multicollinearity | Drop or combine variables; centre before forming interactions; ridge |
| Influential points | Verify data; sensitivity analysis; robust regression (MASS::rlm) |
| Binary or count outcome | Linear regression is the wrong tool; use logistic or Poisson/negative binomial regression |
9 Effect modification (optional extension)
Is the effect of gestational age the same for smokers and non-smokers? Add an interaction:
Code
m_int <- lm(birth_weight_g ~ gestational_age_wks * smoker +
maternal_age_yrs + maternal_bmi + anc_visits +
residence + infant_sex, data = bw)
tidy(m_int, conf.int = TRUE) |>
filter(grepl("gestational_age_wks", term)) |>
mutate(across(where(is.numeric), ~ round(.x, 3))) |>
kable()| term | estimate | std.error | statistic | p.value | conf.low | conf.high |
|---|---|---|---|---|---|---|
| gestational_age_wks | 132.221 | 8.118 | 16.288 | 0.000 | 116.271 | 148.171 |
| gestational_age_wks:smokeryes | -24.557 | 25.592 | -0.960 | 0.338 | -74.839 | 25.726 |
Code
anova(m_multi, m_int)Analysis of Variance Table
Model 1: birth_weight_g ~ gestational_age_wks + maternal_age_yrs + maternal_bmi +
smoker + anc_visits + residence + infant_sex
Model 2: birth_weight_g ~ gestational_age_wks * smoker + maternal_age_yrs +
maternal_bmi + anc_visits + residence + infant_sex
Res.Df RSS Df Sum of Sq F Pr(>F)
1 492 45981716
2 491 45895651 1 86065 0.9207 0.3378The gestational_age_wks:smokeryes term is the difference in slopes between smokers and non-smokers. A large p-value means there is no evidence the gestational age effect differs by smoking status, so the simpler additive model is preferred. (The data were simulated without an interaction.)
10 Prediction
Code
new_mothers <- data.frame(
gestational_age_wks = c(39, 39, 36),
maternal_age_yrs = c(28, 28, 28),
maternal_bmi = c(24, 24, 24),
smoker = c("no", "yes", "no"),
anc_visits = c(4, 4, 4),
residence = c("urban", "urban", "urban"),
infant_sex = c("female", "female", "female")
)
preds <- cbind(
new_mothers[, c("gestational_age_wks", "smoker")],
round(predict(m_multi, new_mothers, interval = "confidence"), 0),
pred_low = round(predict(m_multi, new_mothers, interval = "prediction")[, "lwr"], 0),
pred_high = round(predict(m_multi, new_mothers, interval = "prediction")[, "upr"], 0)
)
kable(preds, col.names = c("GA (wks)", "Smoker", "Fitted", "Mean CI low",
"Mean CI high", "Individual PI low", "Individual PI high"))| GA (wks) | Smoker | Fitted | Mean CI low | Mean CI high | Individual PI low | Individual PI high |
|---|---|---|---|---|---|---|
| 39 | no | 3268 | 3219 | 3317 | 2666 | 3871 |
| 39 | yes | 3012 | 2925 | 3099 | 2405 | 3619 |
| 36 | no | 2879 | 2816 | 2942 | 2275 | 3483 |
Two kinds of interval:
- Confidence interval: uncertainty about the average birth weight of mothers with these characteristics. Narrow.
- Prediction interval: uncertainty about one new infant’s birth weight. Much wider, because it adds the residual scatter (\(\sigma\)).
Predict only within the range of the observed data. Extrapolating a linear model beyond it is risky.
11 Reporting the results
A model write-up should state what was done, the key estimates with confidence intervals, and the checks performed. Example:
We fitted a multiple linear regression of birth weight (g) on gestational age, maternal age, BMI, smoking during pregnancy, antenatal care visits, residence and infant sex (n = 500). Model assumptions (linearity, normality and constant variance of residuals, multicollinearity, influential observations) were assessed graphically and with formal tests, and no serious violations were found. After adjustment, each additional week of gestation was associated with 130 g higher birth weight (95% CI 115 to 145), and maternal smoking with 256 g lower birth weight (95% CI -343 to -170). The model explained 43% of the variance in birth weight (adjusted \(R^2\) = 0.43).
12 Exercises
- Refit the model with
log(birth_weight_g)as the outcome. How do the coefficients now need to be interpreted (hint:100 * (exp(beta) - 1)is the approximate percent change)? Do the diagnostics change? - Change the reference level of
residencetoruralwithrelevel(). What happens to the sign of the coefficient, and what stays the same? - Test whether the effect of
anc_visitsis non-linear by addingI(anc_visits^2). Useanova()to compare. - Standardise all numeric predictors with
scale()and refit. Which predictor has the largest standardised effect? - Add an interaction between
residenceandanc_visits. Does access to care matter more in rural areas? - Introduce an outlier by editing one birth weight to 6500 g. How do Cook’s distance and the coefficients react?
- Add
I((maternal_bmi - mean(maternal_bmi))^2)andI((gestational_age_wks - mean(gestational_age_wks))^2)to the model. Are the quadratic terms statistically significant, and is the curvature large enough to matter in practice (compare predictions with and without them)?