Continuous Variables: Descriptive Statistics, t-tests, ANOVA and Linear Regression

Welcome and how to use these notes

This document teaches the main concepts for analysing continuous outcome variables in health sciences, moving from the very basics (what is a variable? what is a mean?) to inferential tools (confidence intervals, t-tests, ANOVA) and finally modelling (simple and multiple linear regression).

All analyses use free, open data from the World Bank (World Development Indicators, WDI), downloaded directly in R with the WDI package. The unit of analysis is the country; the main outcome is life expectancy at birth, a classic population-health measure.

TipLearning objectives

By the end of these notes you should be able to:

1. Describe and visualise continuous data appropriately.
2. Explain the normal distribution, sampling distributions, standard errors and confidence intervals.
3. State and test hypotheses; interpret p-values, effect sizes and power.
4. Choose and run one-sample, independent and paired t-tests (and their non-parametric alternatives).
5. Run and interpret one-way and two-way ANOVA with post-hoc tests.
6. Measure association with correlation.
7. Fit, interpret, diagnose and validate simple and multiple linear regression models.
8. Report results in a clear, honest, publication-style way.

NoteCourse workflow

Every section follows the same pattern: Concept → Code → Interpretation → Your turn. Run the code chunk by chunk, change things, and break things. That is how statistics is learned.

1 Setup

1.1 Packages

# Run ONCE in your console (not when rendering)
install.packages(c("tidyverse", "WDI", "broom", "car", "lmtest",
                   "sandwich", "pwr", "knitr"))

library(tidyverse)   # data wrangling (dplyr, tidyr) + ggplot2
library(WDI)         # World Bank data
library(broom)       # tidy model output
library(car)         # Levene test, Type II/III ANOVA, VIF
library(lmtest)      # Breusch-Pagan test
library(sandwich)    # robust standard errors
library(pwr)         # power analysis
library(knitr)       # nice tables

theme_set(theme_minimal(base_size = 13))
set.seed(2381)       # reproducibility: same random numbers every run

1.2 Downloading the World Bank data

We pull ten indicators for all economies from 2000 to 2021. We cache the file locally so the notes still work offline after the first successful download.

Short name	WDI code	Description
life_exp	SP.DYN.LE00.IN	Life expectancy at birth (years)
gdp_pc	NY.GDP.PCAP.CD	GDP per capita (current US$)
health_exp_pc	SH.XPD.CHEX.PC.CD	Health expenditure per capita (current US$)
infant_mort	SP.DYN.IMRT.IN	Infant deaths per 1,000 live births
fertility	SP.DYN.TFRT.IN	Births per woman
urban_pct	SP.URB.TOTL.IN.ZS	Urban population (% of total)
water_pct	SH.H2O.BASW.ZS	People using at least basic drinking water (%)
sanitation_pct	SH.STA.BASS.ZS	People using at least basic sanitation (%)
physicians	SH.MED.PHYS.ZS	Physicians per 1,000 people
secondary_enrol	SE.SEC.ENRR	Gross secondary school enrolment (%)

indicators <- c(
  life_exp        = "SP.DYN.LE00.IN",
  gdp_pc          = "NY.GDP.PCAP.CD",
  health_exp_pc   = "SH.XPD.CHEX.PC.CD",
  infant_mort     = "SP.DYN.IMRT.IN",
  fertility       = "SP.DYN.TFRT.IN",
  urban_pct       = "SP.URB.TOTL.IN.ZS",
  water_pct       = "SH.H2O.BASW.ZS",
  sanitation_pct  = "SH.STA.BASS.ZS",
  physicians      = "SH.MED.PHYS.ZS",
  secondary_enrol = "SE.SEC.ENRR"
)

cache_file <- "wdi_hss2381_cache.csv"

if (file.exists(cache_file)) {
  wdi_raw <- read_csv(cache_file, show_col_types = FALSE)
} else {
  wdi_raw <- tryCatch(
    WDI(country = "all", indicator = indicators,
        start = 2000, end = 2021, extra = TRUE),
    error = function(e) NULL
  )
  if (is.null(wdi_raw)) {
    stop("Could not download from the World Bank. Check your internet connection ",
         "and re-run, or place a previously saved 'wdi_hss2381_cache.csv' in this folder.")
  }
  write_csv(wdi_raw, cache_file)
}

dim(wdi_raw)

[1] 5852   22

1.3 Cleaning and preparing the data

World Bank downloads include aggregates (“World”, “Europe & Central Asia”, “High income”…). These are not countries, so we remove them. We also create derived variables used later.

income_levels <- c("Low income", "Lower middle income",
                   "Upper middle income", "High income")

wdi <- wdi_raw |>
  filter(!is.na(region), region != "Aggregates",
         income %in% income_levels) |>
  mutate(
    income    = factor(income, levels = income_levels),   # ordered factor
    region    = factor(region),
    log_gdp   = log(gdp_pc),                              # log transforms (see Section on transformations)
    log_hexp  = log(health_exp_pc),
    log_infant = log(infant_mort)
  )

# Cross-section for 2019 (last pre-COVID year, good data coverage)
d19 <- wdi |>
  filter(year == 2019) |>
  mutate(
    urban_grp = factor(if_else(urban_pct >= median(urban_pct, na.rm = TRUE),
                               "More urban", "Less urban"),
                       levels = c("Less urban", "More urban"))
  )

nrow(d19)            # number of countries in 2019

[1] 213

glimpse(d19 |> select(country, region, income, life_exp, gdp_pc,
                      health_exp_pc, infant_mort, urban_pct, water_pct))

Rows: 213
Columns: 9
$ country       <chr> "Afghanistan", "Albania", "Algeria", "American Samoa", "…
$ region        <fct> South Asia, Europe & Central Asia, Middle East & North A…
$ income        <fct> Low income, Upper middle income, Upper middle income, Hi…
$ life_exp      <dbl> 62.94100, 79.46700, 75.68200, 72.75100, 84.09800, 63.051…
$ gdp_pc        <dbl> 496.6025, 6069.4390, 4468.4534, 12886.1360, 41257.8046, …
$ health_exp_pc <dbl> 74.06424, 366.34015, 242.07050, NA, 3021.07812, 64.44661…
$ infant_mort   <dbl> 56.6, 8.0, 20.7, NA, 2.9, 35.5, 8.0, 8.9, 10.7, NA, 3.2,…
$ urban_pct     <dbl> 25.14373, 56.76391, 72.48246, 83.05422, 88.53918, 66.663…
$ water_pct     <dbl> 69.82270, 94.75822, 91.99980, 99.50730, 100.00000, 64.15…

WarningImportant design note

We use one year (2019) for most analyses so that each country appears once. Using many years would repeat the same country many times, violating the independence assumption of t-tests, ANOVA and regression. Repeated measures are handled with the paired t-test (Section on paired t-tests).

2 Foundations: Data and variables

2.1 Types of variables

Type	Examples (from our data)	Typical summaries
Continuous (interval/ratio)	life expectancy, GDP per capita, urban %	mean, SD, median, IQR
Discrete (count)	number of countries, number of physicians	counts, rates
Categorical – nominal	region	frequencies, proportions
Categorical – ordinal	income group (Low < … < High)	frequencies, median category

A continuous variable can take any value in a range (e.g., 72.43 years). Our course focuses on continuous outcomes (the thing we measure: life expectancy) and continuous or categorical predictors/groups.

d19 |>
  summarise(across(c(life_exp, gdp_pc, urban_pct, region, income),
                   ~ class(.x)[1])) |>
  pivot_longer(everything(), names_to = "variable", values_to = "R class")

# A tibble: 5 × 2
  variable  `R class`
  <chr>     <chr>    
1 life_exp  numeric  
2 gdp_pc    numeric  
3 urban_pct numeric  
4 region    factor   
5 income    factor   

2.2 Outcome vs. explanatory variable

• Dependent / outcome / response (Y): life expectancy.
• Independent / explanatory / predictor (X): GDP per capita, income group, etc.

Choosing statistical methods starts from asking: What is my outcome? What are my groups or predictors? Are observations independent?

2.3 Missing data

d19 |>
  select(life_exp, gdp_pc, health_exp_pc, infant_mort, urban_pct,
         water_pct, sanitation_pct, physicians, secondary_enrol) |>
  summarise(across(everything(), ~ sum(is.na(.x)))) |>
  pivot_longer(everything(), names_to = "variable", values_to = "n_missing") |>
  mutate(pct_missing = round(100 * n_missing / nrow(d19), 1)) |>
  arrange(desc(n_missing))

# A tibble: 9 × 3
  variable        n_missing pct_missing
  <chr>               <int>       <dbl>
1 physicians             78        36.6
2 secondary_enrol        73        34.3
3 health_exp_pc          23        10.8
4 infant_mort            20         9.4
5 sanitation_pct         14         6.6
6 water_pct              11         5.2
7 gdp_pc                  6         2.8
8 life_exp                0         0  
9 urban_pct               0         0  

R drops rows with NA in most models (complete-case analysis). Always report how many observations were actually used. Missingness in country data is rarely random (poorer countries often report less), which can bias results.

3 Descriptive statistics

3.1 Measures of central tendency

• Mean \(\bar{x} = \frac{1}{n}\sum x_i\) — sensitive to outliers.
• Median — middle value; robust.
• Mode — most frequent value (for continuous data, the peak of the density).

x <- na.omit(d19$life_exp)

mean(x)

[1] 73.11814

median(x)

[1] 74.211

# Mode of a continuous variable = peak of the estimated density
dens <- density(x)
dens$x[which.max(dens$y)]

[1] 76.11789

3.2 Measures of spread

• Range, Interquartile range (IQR) = Q3 − Q1
• Variance \(s^2 = \frac{\sum (x_i-\bar{x})^2}{n-1}\) and standard deviation (SD) \(s=\sqrt{s^2}\)
• Coefficient of variation (CV) = SD / mean (unit-free comparison of spread)

range(x)

[1] 31.530 86.151

IQR(x)

[1] 11.4508

var(x)

[1] 62.61835

sd(x)

[1] 7.913175

sd(x) / mean(x)          # coefficient of variation

[1] 0.1082245

# Why n - 1?  Verify by hand:
sum((x - mean(x))^2) / (length(x) - 1)

[1] 62.61835

3.3 Quantiles, percentiles and the five-number summary

quantile(x, probs = c(0.05, 0.25, 0.5, 0.75, 0.95))

      5%      25%      50%      75%      95% 
60.12280 67.33700 74.21100 78.78780 83.33512 

fivenum(x)

[1] 31.5300 67.3370 74.2110 78.7878 86.1510

summary(x)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  31.53   67.34   74.21   73.12   78.79   86.15 

3.4 Shape: skewness and kurtosis

• Skewness > 0: long right tail (e.g., income, GDP); < 0: long left tail (e.g., life expectancy).
• Excess kurtosis: tail heaviness relative to the normal (0 for the normal).

skewness <- function(x) { x <- na.omit(x); m <- mean(x); mean((x - m)^3) / mean((x - m)^2)^1.5 }
kurtosis <- function(x) { x <- na.omit(x); m <- mean(x); mean((x - m)^4) / mean((x - m)^2)^2 - 3 }

tibble(variable = c("life_exp", "gdp_pc", "log_gdp", "infant_mort", "log_infant"),
       skewness = c(skewness(d19$life_exp), skewness(d19$gdp_pc), skewness(d19$log_gdp),
                    skewness(d19$infant_mort), skewness(d19$log_infant)),
       kurtosis = c(kurtosis(d19$life_exp), kurtosis(d19$gdp_pc), kurtosis(d19$log_gdp),
                    kurtosis(d19$infant_mort), kurtosis(d19$log_infant)))

# A tibble: 5 × 3
  variable    skewness kurtosis
  <chr>          <dbl>    <dbl>
1 life_exp     -0.972     2.42 
2 gdp_pc        2.88     11.4  
3 log_gdp      -0.0580   -0.841
4 infant_mort   1.63      4.05 
5 log_infant   -0.182    -1.02 

Notice that logging GDP and infant mortality greatly reduces skewness.

3.5 A complete descriptive table (overall and by group)

desc_fun <- function(data, var) {
  data |>
    summarise(n = sum(!is.na({{ var }})),
              mean = mean({{ var }}, na.rm = TRUE),
              sd = sd({{ var }}, na.rm = TRUE),
              median = median({{ var }}, na.rm = TRUE),
              Q1 = quantile({{ var }}, .25, na.rm = TRUE),
              Q3 = quantile({{ var }}, .75, na.rm = TRUE),
              min = min({{ var }}, na.rm = TRUE),
              max = max({{ var }}, na.rm = TRUE))
}

d19 |> desc_fun(life_exp) |> kable(digits = 2, caption = "Life expectancy, all countries, 2019")

Life expectancy, all countries, 2019

n	mean	sd	median	Q1	Q3	min	max
213	73.12	7.91	74.21	67.34	78.79	31.53	86.15

d19 |> group_by(income) |> desc_fun(life_exp) |>
  kable(digits = 2, caption = "Life expectancy by World Bank income group, 2019")

Life expectancy by World Bank income group, 2019

income	n	mean	sd	median	Q1	Q3	min	max
Low income	25	61.92	7.60	62.18	60.18	65.80	31.53	72.93
Lower middle income	51	67.80	5.64	67.82	64.35	71.80	53.01	78.21
Upper middle income	54	73.29	4.15	73.79	71.26	76.66	62.30	80.30
High income	83	79.64	3.95	81.20	77.01	82.75	68.98	86.15

NoteWhich summary to report?

Roughly symmetric data → mean (SD). Skewed data or outliers → median (IQR). Always give n.

Your turn

1. Compute mean, SD, median and IQR of infant_mort. Why are mean and median so different?
2. Compute the CV of life_exp and gdp_pc. Which one varies more across countries (relatively)?

4 Data visualisation

Always plot before you test.

4.1 Histogram and density

ggplot(d19, aes(life_exp)) +
  geom_histogram(aes(y = after_stat(density)), bins = 20,
                 fill = "steelblue", colour = "white") +
  geom_density(linewidth = 1, colour = "firebrick") +
  geom_vline(xintercept = c(mean(d19$life_exp, na.rm = TRUE),
                            median(d19$life_exp, na.rm = TRUE)),
             linetype = c("solid", "dashed")) +
  labs(x = "Life expectancy at birth (years)", y = "Density",
       title = "Distribution of life expectancy, 2019",
       subtitle = "Solid line = mean, dashed = median (left-skewed: mean < median)")

4.2 Boxplots and violin plots by group

ggplot(d19, aes(income, life_exp, fill = income)) +
  geom_violin(alpha = .4, colour = NA) +
  geom_boxplot(width = .2, outlier.colour = "red") +
  geom_jitter(width = .1, alpha = .3, size = 1) +
  scale_fill_brewer(palette = "Blues") +
  guides(fill = "none") +
  labs(x = NULL, y = "Life expectancy (years)",
       title = "Life expectancy by income group, 2019") +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))

4.3 Scatterplot

p_raw <- ggplot(d19, aes(gdp_pc, life_exp)) +
  geom_point(alpha = .6) + geom_smooth(se = FALSE) +
  labs(x = "GDP per capita (US$)", y = "Life expectancy", title = "Raw scale: curved")

p_log <- ggplot(d19, aes(log_gdp, life_exp)) +
  geom_point(alpha = .6) + geom_smooth(method = "lm", se = TRUE) +
  labs(x = "log(GDP per capita)", y = "Life expectancy", title = "Log scale: roughly linear")

p_raw

p_log

4.4 Empirical cumulative distribution (ECDF) and Q-Q plot

ggplot(d19, aes(life_exp)) + stat_ecdf(geom = "step") +
  labs(x = "Life expectancy", y = "Proportion of countries \u2264 x",
       title = "Empirical CDF")

ggplot(d19, aes(sample = life_exp)) + stat_qq() + stat_qq_line(colour = "red") +
  labs(title = "Normal Q-Q plot of life expectancy",
       subtitle = "Points on the line = normal; curved ends = skew / heavy tails")

4.5 Correct reading of plots

Plot	Use
Histogram / density	shape, centre, spread, outliers
Boxplot / violin	compare groups
Scatterplot	relationship between two continuous variables
Q-Q plot	check normality
Residual plots	check regression assumptions (later)

5 The normal distribution and z-scores

The normal distribution \(N(\mu,\sigma^2)\) is bell-shaped and symmetric. The empirical rule: about 68%, 95%, 99.7% of values lie within 1, 2, 3 SD of the mean.

5.1 z-scores

\[z = \frac{x - \bar{x}}{s}\]

A z-score tells us how many SDs an observation is from the mean, allowing comparison across different scales.

m <- mean(d19$life_exp, na.rm = TRUE)
s <- sd(d19$life_exp,   na.rm = TRUE)

can <- d19 |> filter(country == "Canada") |> pull(life_exp)
z_can <- (can - m) / s
c(canada = can, z = z_can)

   canada         z 
82.158537  1.142448 

# Under normality, what percentile is Canada?
pnorm(z_can)

[1] 0.8733661

# Empirical percentile (no normality assumption)
ecdf(d19$life_exp)(can)

[1] 0.8779343

5.2 Probabilities with pnorm, qnorm, dnorm

# P(life expectancy < 65) under a normal model fitted to the data
pnorm(65, mean = m, sd = s)

[1] 0.1524688

# vs actual proportion
mean(d19$life_exp < 65, na.rm = TRUE)

[1] 0.1596244

# The value that cuts off the lowest 10%
qnorm(0.10, mean = m, sd = s)

[1] 62.977

# Verify the empirical rule
mean(abs(scale(d19$life_exp)) < 1, na.rm = TRUE)  # ~68%

[1] 0.6431925

mean(abs(scale(d19$life_exp)) < 2, na.rm = TRUE)  # ~95%

[1] 0.9812207

Life expectancy is somewhat left-skewed, so the normal model is only approximate — a good reminder to check, not assume.

6 Sampling distributions and the Central Limit Theorem

A statistic (e.g., the sample mean) varies from sample to sample. Its distribution is the sampling distribution. Its standard deviation is the standard error (SE):

\[SE(\bar{x}) = \frac{s}{\sqrt{n}}\]

6.1 Simulation: the Central Limit Theorem (CLT)

The CLT says: for large enough \(n\), the sampling distribution of \(\bar{x}\) is approximately normal even if the population is not. We treat our skewed infant_mort as the “population”.

pop <- na.omit(d19$infant_mort)       # a strongly right-skewed "population"

sizes <- c(2, 5, 30, 100)

# NOTE: the function argument is called `size` and the column `sample_size`.
# (Naming both `n` made tibble() use the new text column instead of the number.)
sim <- map_dfr(sizes, function(size) {
  tibble(sample_size = paste("n =", size),
         mean = replicate(3000, mean(sample(pop, size, replace = TRUE))))
}) |>
  mutate(sample_size = factor(sample_size, levels = paste("n =", sizes)))

ggplot(tibble(x = pop), aes(x)) +
  geom_histogram(bins = 30, fill = "grey60", colour = "white") +
  labs(title = "Population: infant mortality (skewed)", x = "Deaths per 1,000")

ggplot(sim, aes(mean)) +
  geom_histogram(bins = 30, fill = "steelblue", colour = "white") +
  facet_wrap(~ sample_size, scales = "free_y") +
  labs(title = "Sampling distribution of the mean becomes normal as n grows",
       x = "Sample mean")

# SE shrinks like 1/sqrt(n): simulated vs theoretical
sim |>
  group_by(sample_size) |>
  summarise(empirical_SE = sd(mean), .groups = "drop") |>
  mutate(theoretical_SE = sd(pop) / sqrt(sizes))

# A tibble: 4 × 3
  sample_size empirical_SE theoretical_SE
  <fct>              <dbl>          <dbl>
1 n = 2              13.9           13.8 
2 n = 5               8.69           8.75
3 n = 30              3.57           3.57
4 n = 100             1.95           1.96

7 Confidence intervals

A 95% confidence interval (CI) is produced by a procedure that captures the true parameter in 95% of repeated samples. It is not “a 95% probability the parameter is in this interval”.

\[\bar{x} \pm t_{1-\alpha/2,\;n-1}\times \frac{s}{\sqrt{n}}\]

7.1 By hand and with t.test()

x <- na.omit(d19$life_exp)
n <- length(x)
se <- sd(x) / sqrt(n)
tcrit <- qt(0.975, df = n - 1)

c(lower = mean(x) - tcrit * se, upper = mean(x) + tcrit * se)

   lower    upper 
72.04935 74.18694 

t.test(x)$conf.int          # identical

[1] 72.04935 74.18694
attr(,"conf.level")
[1] 0.95

7.2 Bootstrap confidence interval (no formula, few assumptions)

boot_means <- replicate(5000, mean(sample(x, replace = TRUE)))
quantile(boot_means, c(0.025, 0.975))

    2.5%    97.5% 
72.07644 74.15546 

# Bootstrap CI for a median (no simple formula exists)
boot_med <- replicate(5000, median(sample(x, replace = TRUE)))
quantile(boot_med, c(0.025, 0.975))

  2.5%  97.5% 
72.625 75.809 

7.3 Visualising what “95% confidence” means

true_mu <- mean(pop)
cis <- map_dfr(1:50, function(i) {
  smp <- sample(pop, 30)
  ci  <- t.test(smp)$conf.int
  tibble(i = i, mean = mean(smp), lo = ci[1], hi = ci[2],
         covers = lo <= true_mu & true_mu <= hi)
})

ggplot(cis, aes(i, mean, colour = covers)) +
  geom_pointrange(aes(ymin = lo, ymax = hi)) +
  geom_hline(yintercept = true_mu, linetype = "dashed") +
  coord_flip() +
  scale_colour_manual(values = c(`TRUE` = "steelblue", `FALSE` = "red")) +
  labs(title = "50 CIs from n = 30 samples", subtitle = "About 95% should cover the true mean",
       x = "Sample", y = "Mean infant mortality")

mean(cis$covers)

[1] 0.98

8 Logic of hypothesis testing

8.1 The framework

1. State H₀ (null: no difference/effect) and H₁ (alternative).
2. Choose α (usually 0.05) — the tolerated Type I error rate.
3. Compute a test statistic and p-value.
4. Decide and interpret in context, together with the effect size and CI.

	H₀ true	H₀ false
Reject H₀	Type I error (α, false positive)	Correct (power = 1 − β)
Do not reject H₀	Correct	Type II error (β, false negative)

ImportantWhat a p-value is (and is not)

The p-value is the probability of data at least as extreme as observed, assuming H₀ is true. It is not the probability that H₀ is true, not the size of the effect, and not a measure of importance. Statistical significance ≠ practical (clinical/public-health) significance.

8.2 Simulation: Type I error and power

# Type I error: both groups from the SAME population (H0 true)
p_null <- replicate(5000, t.test(rnorm(30), rnorm(30))$p.value)
mean(p_null < 0.05)                   # should be ~0.05

[1] 0.0536

# Power: true difference of 0.5 SD
p_alt <- replicate(5000, t.test(rnorm(30), rnorm(30, mean = 0.5))$p.value)
mean(p_alt < 0.05)                    # power with n = 30 per group

[1] 0.474

hist(p_null, main = "p-values when H0 is true are uniform", xlab = "p-value")

8.3 Multiple testing

If you run many tests, false positives accumulate: with 20 tests at α = 0.05, P(at least one false positive) ≈ 64%.

1 - (1 - 0.05)^20

[1] 0.6415141

# 20 tests on pure noise
pvals <- replicate(20, t.test(rnorm(30), rnorm(30))$p.value)
sum(pvals < 0.05)                              # unadjusted "discoveries"

[1] 2

sum(p.adjust(pvals, method = "bonferroni") < 0.05)

[1] 0

sum(p.adjust(pvals, method = "holm") < 0.05)

[1] 0

sum(p.adjust(pvals, method = "BH") < 0.05)     # false discovery rate control

[1] 0

8.4 Effect sizes

Measure	Use	Small / Medium / Large (Cohen)
Cohen’s d / Hedges’ g	difference between two means	0.2 / 0.5 / 0.8
η² (eta squared)	ANOVA: proportion of variance explained	0.01 / 0.06 / 0.14
r	correlation	0.1 / 0.3 / 0.5
R²	regression variance explained	0.02 / 0.13 / 0.26

9 Checking assumptions

Parametric tests (t-test, ANOVA, linear regression) share three core assumptions:

1. Independence of observations (determined by design, not testable by a formula).
2. Normality — of the outcome within groups (t-test/ANOVA) or of the residuals (regression). Matters most for small samples.
3. Homogeneity of variance (equal spread across groups / across fitted values).

9.1 Normality

# Formal test (H0: data are normal). With large n, trivial departures give tiny p-values!
shapiro.test(d19$life_exp)

    Shapiro-Wilk normality test

data:  d19$life_exp
W = 0.94039, p-value = 1.167e-07

shapiro.test(d19$log_infant)

    Shapiro-Wilk normality test

data:  d19$log_infant
W = 0.96447, p-value = 8.51e-05

# Normality WITHIN groups (what t-test/ANOVA require)
d19 |>
  filter(!is.na(life_exp)) |>
  group_by(income) |>
  summarise(shapiro_p = shapiro.test(life_exp)$p.value, n = n())

# A tibble: 4 × 3
  income              shapiro_p     n
  <fct>                   <dbl> <int>
1 Low income          0.0000248    25
2 Lower middle income 0.676        51
3 Upper middle income 0.0272       54
4 High income         0.000297     83

# Visual check within groups
ggplot(d19, aes(sample = life_exp)) + stat_qq() + stat_qq_line() +
  facet_wrap(~ income) + labs(title = "Q-Q plots by income group")

NoteRule of thumb

Rely on plots first; tests are secondary. With \(n > 30\) per group the t-test/ANOVA are fairly robust to non-normality (CLT). With strong skew or outliers and small n, use transformations or non-parametric tests.

9.2 Homogeneity of variance

# Levene's test (robust to non-normality): H0 = equal variances
leveneTest(life_exp ~ income, data = d19)

Levene's Test for Homogeneity of Variance (center = median)
       Df F value  Pr(>F)  
group   3  2.2581 0.08273 .
      209                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Bartlett's test (assumes normality; more sensitive)
bartlett.test(life_exp ~ income, data = d19)

    Bartlett test of homogeneity of variances

data:  life_exp by income
Bartlett's K-squared = 23.84, df = 3, p-value = 2.697e-05

# Rule of thumb: ratio of largest to smallest SD < 2
d19 |> group_by(income) |> summarise(sd = sd(life_exp, na.rm = TRUE)) |>
  summarise(ratio = max(sd) / min(sd))

# A tibble: 1 × 1
  ratio
  <dbl>
1  1.93

9.3 Outliers

# Tukey's rule: beyond 1.5 * IQR from the quartiles
d19 |>
  filter(!is.na(life_exp)) |>
  mutate(q1 = quantile(life_exp, .25), q3 = quantile(life_exp, .75),
         iqr = q3 - q1,
         outlier = life_exp < q1 - 1.5 * iqr | life_exp > q3 + 1.5 * iqr) |>
  filter(outlier) |>
  select(country, income, life_exp)

# A tibble: 1 × 3
  country                  income     life_exp
  <chr>                    <fct>         <dbl>
1 Central African Republic Low income     31.5

Never delete outliers just because they are inconvenient. Check for data errors; if genuine, keep them and use robust methods, or report results with and without them (sensitivity analysis).

9.4 Transformations

For right-skewed positive data (GDP, health spending, mortality), a log transformation often helps. Differences on the log scale correspond to ratios on the original scale.

par(mfrow = c(1, 2))
hist(d19$gdp_pc,  main = "GDP per capita", xlab = "US$", col = "grey80")
hist(d19$log_gdp, main = "log(GDP per capita)", xlab = "log US$", col = "steelblue")

par(mfrow = c(1, 1))

10 One-sample t-test

Question: Does the mean of a continuous variable differ from a known/benchmark value \(\mu_0\)?

\[t = \frac{\bar{x} - \mu_0}{s/\sqrt{n}}, \quad df = n-1\]

Example: Is mean country life expectancy in 2019 different from a (hypothetical) policy benchmark of 70 years?

• H₀: μ = 70 H₁: μ ≠ 70

t1 <- t.test(d19$life_exp, mu = 70)
t1

    One Sample t-test

data:  d19$life_exp
t = 5.7509, df = 212, p-value = 3.071e-08
alternative hypothesis: true mean is not equal to 70
95 percent confidence interval:
 72.04935 74.18694
sample estimates:
mean of x 
 73.11814 

# By hand
x <- na.omit(d19$life_exp)
t_stat <- (mean(x) - 70) / (sd(x) / sqrt(length(x)))
t_stat

[1] 5.750892

2 * pt(-abs(t_stat), df = length(x) - 1)    # two-sided p-value

[1] 3.071403e-08

# One-sided alternative: is the mean GREATER than 70?
t.test(d19$life_exp, mu = 70, alternative = "greater")

    One Sample t-test

data:  d19$life_exp
t = 5.7509, df = 212, p-value = 1.536e-08
alternative hypothesis: true mean is greater than 70
95 percent confidence interval:
 72.22239      Inf
sample estimates:
mean of x 
 73.11814 

# Effect size (Cohen's d for one sample)
(mean(x) - 70) / sd(x)

[1] 0.3940446

Interpretation template: “The mean life expectancy across 2019 countries was 73.1 years (SD = 7.9, n = 213), which differed significantly from the benchmark of 70 years, t(212) = 5.75, p = 3.07e-08, 95% CI [72, 74.2].”

NoteCIs and tests are two views of the same thing

If the 95% CI excludes \(\mu_0\), then p < 0.05 for a two-sided test at \(\mu_0\).

Your turn

Test whether mean urban_pct in 2019 differs from 50%. Check normality first.

11 Independent two-sample t-test

Question: Do two independent groups have different population means?

\[t = \frac{\bar{x}_1-\bar{x}_2}{SE_{diff}}\]

• Student’s t-test: assumes equal variances, \(SE = s_p\sqrt{1/n_1+1/n_2}\).
• Welch’s t-test (R default): does not assume equal variances; uses adjusted df. Preferred by default.

Example: Do high-income countries have different life expectancy than upper-middle-income countries?

two <- d19 |>
  filter(income %in% c("Upper middle income", "High income"), !is.na(life_exp)) |>
  mutate(income = droplevels(income))

two |> group_by(income) |> desc_fun(life_exp) |> kable(digits = 2)

income	n	mean	sd	median	Q1	Q3	min	max
Upper middle income	54	73.29	4.15	73.79	71.26	76.66	62.30	80.30
High income	83	79.64	3.95	81.20	77.01	82.75	68.98	86.15

ggplot(two, aes(income, life_exp, fill = income)) +
  geom_boxplot(alpha = .6) + geom_jitter(width = .1, alpha = .4) +
  guides(fill = "none") + labs(x = NULL, y = "Life expectancy")

11.1 Assumption checks, then the test

leveneTest(life_exp ~ income, data = two)

Levene's Test for Homogeneity of Variance (center = median)
       Df F value Pr(>F)
group   1  0.0296 0.8637
      135               

# Welch (default, unequal variances)
welch <- t.test(life_exp ~ income, data = two)
welch

    Welch Two Sample t-test

data:  life_exp by income
t = -8.9237, df = 109.36, p-value = 1.182e-14
alternative hypothesis: true difference in means between group Upper middle income and group High income is not equal to 0
95 percent confidence interval:
 -7.758169 -4.938351
sample estimates:
mean in group Upper middle income         mean in group High income 
                         73.29464                          79.64290 

# Student (equal variances)
t.test(life_exp ~ income, data = two, var.equal = TRUE)

    Two Sample t-test

data:  life_exp by income
t = -9.0175, df = 135, p-value = 1.669e-15
alternative hypothesis: true difference in means between group Upper middle income and group High income is not equal to 0
95 percent confidence interval:
 -7.740545 -4.955975
sample estimates:
mean in group Upper middle income         mean in group High income 
                         73.29464                          79.64290 

11.2 Effect size and tidy output

cohens_d <- function(x, y) {
  nx <- length(x); ny <- length(y)
  sp <- sqrt(((nx - 1) * var(x) + (ny - 1) * var(y)) / (nx + ny - 2))
  d  <- (mean(x) - mean(y)) / sp
  J  <- 1 - 3 / (4 * (nx + ny) - 9)          # small-sample correction
  c(cohens_d = d, hedges_g = d * J)
}

g1 <- two$life_exp[two$income == "High income"]
g2 <- two$life_exp[two$income == "Upper middle income"]
cohens_d(g1, g2)

cohens_d hedges_g 
1.576557 1.567782 

tidy(welch)    # estimate = difference in means, with CI, df, p

# A tibble: 1 × 10
  estimate estimate1 estimate2 statistic  p.value parameter conf.low conf.high
     <dbl>     <dbl>     <dbl>     <dbl>    <dbl>     <dbl>    <dbl>     <dbl>
1    -6.35      73.3      79.6     -8.92 1.18e-14      109.    -7.76     -4.94
# ℹ 2 more variables: method <chr>, alternative <chr>

Interpretation: the estimated difference is 6.35 years in favour of high-income countries. Always interpret the size and CI of the difference, not just the p-value.

11.3 Skewed outcome: log-transform or non-parametric

# Infant mortality is skewed: compare on the log scale
t_log <- t.test(log_infant ~ income, data = d19 |>
                  filter(income %in% c("Upper middle income", "High income")) |>
                  mutate(income = droplevels(income)))
t_log

    Welch Two Sample t-test

data:  log_infant by income
t = 8.7504, df = 113.54, p-value = 2.276e-14
alternative hypothesis: true difference in means between group Upper middle income and group High income is not equal to 0
95 percent confidence interval:
 0.8283516 1.3131933
sample estimates:
mean in group Upper middle income         mean in group High income 
                         2.545743                          1.474970 

# Back-transform: difference of logs -> RATIO of geometric means
exp(t_log$estimate[1] - t_log$estimate[2])     # UMIC / HIC ratio

mean in group Upper middle income 
                         2.917632 

exp(t_log$conf.int)                             # 95% CI for the ratio

[1] 2.289542 3.718028
attr(,"conf.level")
[1] 0.95

11.4 Non-parametric alternative: Mann–Whitney U (Wilcoxon rank-sum)

Compares distributions/ranks and does not require normality. It tests whether one group tends to have larger values (a difference in medians only if shapes are similar).

wilcox.test(life_exp ~ income, data = two, conf.int = TRUE)

    Wilcoxon rank sum test with continuity correction

data:  life_exp by income
W = 603, p-value = 5.486e-13
alternative hypothesis: true location shift is not equal to 0
95 percent confidence interval:
 -7.654188 -5.107064
sample estimates:
difference in location 
             -6.265344 

Your turn

Compare fertility between “Low income” and “Lower middle income”. Check assumptions, run Welch’s t-test, compute Hedges’ g, and write one interpretation sentence.

12 Paired t-test

Question: Is the mean of within-pair differences zero? Use for before/after, matched, or repeated measures on the same unit.

\[t=\frac{\bar{d}}{s_d/\sqrt{n}},\quad df=n-1\]

Example: Did the same countries’ life expectancy change between 2000 and 2019?

paired <- wdi |>
  filter(year %in% c(2000, 2019)) |>
  select(iso3c, country, income, year, life_exp) |>
  pivot_wider(names_from = year, values_from = life_exp, names_prefix = "y") |>
  drop_na(y2000, y2019) |>
  mutate(diff = y2019 - y2000)

nrow(paired)

[1] 213

summary(paired$diff)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-15.079   3.285   4.753   5.467   6.857  18.893 

# Visualise the pairing (each line = one country)
paired |>
  pivot_longer(c(y2000, y2019), names_to = "year", values_to = "life_exp") |>
  ggplot(aes(year, life_exp, group = country)) +
  geom_line(alpha = .2) + geom_point(alpha = .3) +
  labs(x = NULL, y = "Life expectancy", title = "Each country: 2000 vs 2019")

12.1 Assumption and test

The assumption is on the differences, not on each time point.

ggplot(paired, aes(diff)) + geom_histogram(bins = 20, fill = "steelblue", colour = "white") +
  labs(x = "Change in life expectancy (2019 - 2000)")

qqnorm(paired$diff); qqline(paired$diff)

pt_res <- t.test(paired$y2019, paired$y2000, paired = TRUE)
pt_res

    Paired t-test

data:  paired$y2019 and paired$y2000
t = 20.998, df = 212, p-value < 2.2e-16
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
 4.953834 5.980280
sample estimates:
mean difference 
       5.467057 

# Equivalent: one-sample t-test on the differences
t.test(paired$diff, mu = 0)

    One Sample t-test

data:  paired$diff
t = 20.998, df = 212, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 4.953834 5.980280
sample estimates:
mean of x 
 5.467057 

# Effect size for paired data (Cohen's dz)
mean(paired$diff) / sd(paired$diff)

[1] 1.438772

12.2 Why pairing matters

# Wrong analysis (ignores pairing): usually loses power and gives wrong SE
t.test(paired$y2019, paired$y2000, paired = FALSE)$p.value

[1] 4.140333e-10

t.test(paired$y2019, paired$y2000, paired = TRUE)$p.value

[1] 1.093216e-53

# Correlation between the two time points shows why pairing helps
cor(paired$y2000, paired$y2019)

[1] 0.9243386

12.3 Non-parametric alternative: Wilcoxon signed-rank test

wilcox.test(paired$y2019, paired$y2000, paired = TRUE, conf.int = TRUE)

    Wilcoxon signed rank test with continuity correction

data:  paired$y2019 and paired$y2000
V = 22568, p-value < 2.2e-16
alternative hypothesis: true location shift is not equal to 0
95 percent confidence interval:
 4.691552 5.530227
sample estimates:
(pseudo)median 
      5.090226 

Your turn

Did mean infant_mort (log scale) decrease between 2000 and 2019? Back-transform to report a percentage reduction.

13 Power and sample size

Power = probability of detecting a true effect. It depends on effect size, sample size, α, and variability. Plan sample size before collecting data.

# Sample size per group to detect d = 0.5 with 80% power (two-sample t-test)
pwr.t.test(d = 0.5, sig.level = 0.05, power = 0.80, type = "two.sample")

     Two-sample t test power calculation 

              n = 63.76561
              d = 0.5
      sig.level = 0.05
          power = 0.8
    alternative = two.sided

NOTE: n is number in *each* group

# Power achieved with n = 25 per group for d = 0.5
pwr.t.test(n = 25, d = 0.5, sig.level = 0.05, type = "two.sample")$power

[1] 0.4101003

# Paired design: dz = 0.5
pwr.t.test(d = 0.5, sig.level = 0.05, power = 0.80, type = "paired")

     Paired t test power calculation 

              n = 33.36713
              d = 0.5
      sig.level = 0.05
          power = 0.8
    alternative = two.sided

NOTE: n is number of *pairs*

# Base R equivalent using raw units: detect a 3-year difference when SD = 6
power.t.test(delta = 3, sd = 6, sig.level = 0.05, power = 0.80)

     Two-sample t test power calculation 

              n = 63.76576
          delta = 3
             sd = 6
      sig.level = 0.05
          power = 0.8
    alternative = two.sided

NOTE: n is number in *each* group

# One-way ANOVA with 4 groups, f = 0.25 (medium effect)
pwr.anova.test(k = 4, f = 0.25, sig.level = 0.05, power = 0.80)

     Balanced one-way analysis of variance power calculation 

              k = 4
              n = 44.59927
              f = 0.25
      sig.level = 0.05
          power = 0.8

NOTE: n is number in each group

# Regression: R2 = 0.13 => f2 = R2/(1-R2); u = number of predictors
pwr.f2.test(u = 3, f2 = 0.13 / (1 - 0.13), sig.level = 0.05, power = 0.80)  # v = n - u - 1

     Multiple regression power calculation 

              u = 3
              v = 72.98495
             f2 = 0.1494253
      sig.level = 0.05
          power = 0.8

# Power curve
ns <- seq(10, 150, 5)
tibble(n = ns,
       `d = 0.2` = map_dbl(ns, ~ pwr.t.test(n = .x, d = .2)$power),
       `d = 0.5` = map_dbl(ns, ~ pwr.t.test(n = .x, d = .5)$power),
       `d = 0.8` = map_dbl(ns, ~ pwr.t.test(n = .x, d = .8)$power)) |>
  pivot_longer(-n, names_to = "effect", values_to = "power") |>
  ggplot(aes(n, power, colour = effect)) + geom_line(linewidth = 1) +
  geom_hline(yintercept = .8, linetype = "dashed") +
  labs(x = "n per group", y = "Power", title = "Power curves for two-sample t-test")

14 Correlation

Correlation measures the strength and direction of a linear (Pearson) or monotonic (Spearman, Kendall) association between two continuous variables. It ranges from −1 to +1 and is unit-free.

\[r=\frac{\sum (x_i-\bar{x})(y_i-\bar{y})}{(n-1)s_xs_y}\]

14.1 Why you must plot: Anscombe’s quartet

anscombe_long <- anscombe |>
  pivot_longer(everything(), names_to = c(".value", "set"), names_pattern = "(.)(.)")

anscombe_long |> group_by(set) |>
  summarise(mean_x = mean(x), mean_y = mean(y), r = cor(x, y))   # nearly identical!

# A tibble: 4 × 4
  set   mean_x mean_y     r
  <chr>  <dbl>  <dbl> <dbl>
1 1          9   7.50 0.816
2 2          9   7.50 0.816
3 3          9   7.5  0.816
4 4          9   7.50 0.817

ggplot(anscombe_long, aes(x, y)) + geom_point() + geom_smooth(method = "lm", se = FALSE) +
  facet_wrap(~ set) + labs(title = "Same r, completely different data")

14.2 Pearson, Spearman and Kendall

cor.test(d19$life_exp, d19$log_gdp)                       # Pearson (linear)

    Pearson's product-moment correlation

data:  d19$life_exp and d19$log_gdp
t = 21.45, df = 205, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.7843260 0.8694847
sample estimates:
      cor 
0.8317348 

cor.test(d19$life_exp, d19$gdp_pc, method = "spearman")   # rank-based, robust to monotone curvature

    Spearman's rank correlation rho

data:  d19$life_exp and d19$gdp_pc
S = 196954, p-value < 2.2e-16
alternative hypothesis: true rho is not equal to 0
sample estimates:
     rho 
0.866766 

cor.test(d19$life_exp, d19$log_gdp, method = "kendall")   # good for small n/ties

    Kendall's rank correlation tau

data:  d19$life_exp and d19$log_gdp
z = 14.412, p-value < 2.2e-16
alternative hypothesis: true tau is not equal to 0
sample estimates:
      tau 
0.6734675 

# r on raw vs log GDP: linearity matters for Pearson
c(raw = cor(d19$life_exp, d19$gdp_pc, use = "complete.obs"),
  log = cor(d19$life_exp, d19$log_gdp, use = "complete.obs"))

      raw       log 
0.6041549 0.8317348 

14.3 Correlation matrix

vars <- d19 |> select(life_exp, log_gdp, log_hexp, log_infant, fertility,
                      urban_pct, water_pct, sanitation_pct, secondary_enrol)
cm <- cor(vars, use = "pairwise.complete.obs")
round(cm, 2)

                life_exp log_gdp log_hexp log_infant fertility urban_pct
life_exp            1.00    0.83     0.82      -0.90     -0.82      0.60
log_gdp             0.83    1.00     0.96      -0.86     -0.76      0.66
log_hexp            0.82    0.96     1.00      -0.86     -0.77      0.67
log_infant         -0.90   -0.86    -0.86       1.00      0.78     -0.63
fertility          -0.82   -0.76    -0.77       0.78      1.00     -0.51
urban_pct           0.60    0.66     0.67      -0.63     -0.51      1.00
water_pct           0.80    0.77     0.76      -0.74     -0.84      0.60
sanitation_pct      0.84    0.76     0.78      -0.78     -0.83      0.56
secondary_enrol     0.70    0.74     0.76      -0.70     -0.70      0.52
                water_pct sanitation_pct secondary_enrol
life_exp             0.80           0.84            0.70
log_gdp              0.77           0.76            0.74
log_hexp             0.76           0.78            0.76
log_infant          -0.74          -0.78           -0.70
fertility           -0.84          -0.83           -0.70
urban_pct            0.60           0.56            0.52
water_pct            1.00           0.90            0.74
sanitation_pct       0.90           1.00            0.74
secondary_enrol      0.74           0.74            1.00

as.data.frame(as.table(cm)) |>
  ggplot(aes(Var1, Var2, fill = Freq)) +
  geom_tile() + geom_text(aes(label = round(Freq, 2)), size = 3) +
  scale_fill_gradient2(low = "firebrick", mid = "white", high = "steelblue",
                       midpoint = 0, limits = c(-1, 1)) +
  labs(x = NULL, y = NULL, fill = "r") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

WarningCorrelation ≠ causation (and ecological fallacy)

• A strong correlation between GDP and life expectancy does not prove wealth causes longevity (confounding, reverse causation, shared causes such as education, institutions).
• Ecological fallacy: associations between country averages need not hold for individuals within countries.
• Correlation is sensitive to outliers and restricted ranges, and captures only monotonic/linear patterns.

15 One-way ANOVA

Question: Do the means of three or more independent groups differ?

Why not many t-tests? Six pairwise tests among 4 groups inflate the family-wise Type I error. ANOVA tests everything with one omnibus test:

• H₀: μ₁ = μ₂ = … = μₖ
• H₁: at least one mean differs

Idea: partition total variability.

\[SS_{total}=SS_{between}+SS_{within},\qquad F=\frac{MS_{between}}{MS_{within}}=\frac{SS_B/(k-1)}{SS_W/(N-k)}\]

If H₀ is true, \(F\approx 1\). Large F → group means differ more than expected by chance.

Example: Does mean life expectancy differ across the four World Bank income groups?

15.1 Step 1: Explore

anova_df <- d19 |> filter(!is.na(life_exp), !is.na(income))

anova_df |> group_by(income) |> desc_fun(life_exp) |> kable(digits = 2)

income	n	mean	sd	median	Q1	Q3	min	max
Low income	25	61.92	7.60	62.18	60.18	65.80	31.53	72.93
Lower middle income	51	67.80	5.64	67.82	64.35	71.80	53.01	78.21
Upper middle income	54	73.29	4.15	73.79	71.26	76.66	62.30	80.30
High income	83	79.64	3.95	81.20	77.01	82.75	68.98	86.15

ggplot(anova_df, aes(income, life_exp, fill = income)) +
  geom_boxplot(alpha = .6, outlier.shape = NA) +
  geom_jitter(width = .15, alpha = .4) +
  stat_summary(fun = mean, geom = "point", shape = 23, size = 3, fill = "white") +
  scale_fill_brewer(palette = "Blues") + guides(fill = "none") +
  labs(x = NULL, y = "Life expectancy (years)",
       title = "Life expectancy by income group", subtitle = "Diamonds = group means")

15.2 Step 2: Fit the model

fit_aov <- aov(life_exp ~ income, data = anova_df)
summary(fit_aov)

             Df Sum Sq Mean Sq F value Pr(>F)    
income        3   8111  2703.5   109.4 <2e-16 ***
Residuals   209   5165    24.7                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Same model as a linear model (ANOVA IS regression with dummy variables)
fit_lm <- lm(life_exp ~ income, data = anova_df)
anova(fit_lm)

Analysis of Variance Table

Response: life_exp
           Df Sum Sq Mean Sq F value    Pr(>F)    
income      3 8110.6 2703.52  109.41 < 2.2e-16 ***
Residuals 209 5164.5   24.71                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

tidy(fit_aov)

# A tibble: 2 × 6
  term         df sumsq meansq statistic   p.value
  <chr>     <dbl> <dbl>  <dbl>     <dbl>     <dbl>
1 income        3 8111. 2704.       109.  1.29e-42
2 Residuals   209 5165.   24.7       NA  NA       

15.3 Step 3: Check assumptions

par(mfrow = c(1, 2))
plot(fit_aov, which = 1)    # residuals vs fitted: equal spread?
plot(fit_aov, which = 2)    # Q-Q of residuals

par(mfrow = c(1, 1))

shapiro.test(residuals(fit_aov))

    Shapiro-Wilk normality test

data:  residuals(fit_aov)
W = 0.92341, p-value = 4.482e-09

leveneTest(life_exp ~ income, data = anova_df)

Levene's Test for Homogeneity of Variance (center = median)
       Df F value  Pr(>F)  
group   3  2.2581 0.08273 .
      209                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

15.4 Step 4: Effect size

ss <- summary(fit_aov)[[1]][["Sum Sq"]]
df_aov <- summary(fit_aov)[[1]][["Df"]]
ms_within <- ss[2] / df_aov[2]

eta_sq   <- ss[1] / sum(ss)                                   # biased upward
omega_sq <- (ss[1] - df_aov[1] * ms_within) / (sum(ss) + ms_within) # less biased
c(eta_squared = eta_sq, omega_squared = omega_sq)

  eta_squared omega_squared 
    0.6109620     0.6042529 

15.5 Step 5: Post-hoc comparisons (only if the omnibus F is significant)

# Tukey HSD: all pairwise comparisons, controls family-wise error
tuk <- TukeyHSD(fit_aov)
tuk

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = life_exp ~ income, data = anova_df)

$income
                                             diff       lwr       upr    p adj
Lower middle income-Low income           5.876697  2.733563  9.019831 1.48e-05
Upper middle income-Low income          11.371204  8.256922 14.485485 0.00e+00
High income-Low income                  17.719464 14.782396 20.656532 0.00e+00
Upper middle income-Lower middle income  5.494506  2.980746  8.008267 3.00e-07
High income-Lower middle income         11.842767  9.552218 14.133315 0.00e+00
High income-Upper middle income          6.348260  4.097467  8.599053 0.00e+00

plot(tuk, las = 1)   # intervals not crossing 0 = significant

# Alternative: pairwise t-tests with Holm correction
pairwise.t.test(anova_df$life_exp, anova_df$income, p.adjust.method = "holm")

    Pairwise comparisons using t tests with pooled SD 

data:  anova_df$life_exp and anova_df$income 

                    Low income Lower middle income Upper middle income
Lower middle income 2.5e-06    -                   -                  
Upper middle income < 2e-16    9.9e-08             -                  
High income         < 2e-16    < 2e-16             1.7e-11            

P value adjustment method: holm 

15.6 Planned contrasts (hypothesis-driven alternative to post-hoc)

If you specify questions in advance, use orthogonal contrasts: more power than testing everything.

# Level order: Low, Lower-middle, Upper-middle, High
contr <- cbind(
  "High vs rest"     = c(-1, -1, -1,  3),
  "Low vs LM & UM"   = c( 2, -1, -1,  0),
  "LM vs UM"         = c( 0,  1, -1,  0)
)
crossprod(contr)          # off-diagonals = 0 -> orthogonal

               High vs rest Low vs LM & UM LM vs UM
High vs rest             12              0        0
Low vs LM & UM            0              6        0
LM vs UM                  0              0        2

anova_c <- anova_df
contrasts(anova_c$income) <- contr
fit_c <- aov(life_exp ~ income, data = anova_c)
summary(fit_c, split = list(income = list("High vs rest" = 1,
                                          "Low vs LM & UM" = 2,
                                          "LM vs UM" = 3)))

                          Df Sum Sq Mean Sq F value   Pr(>F)    
income                     3   8111    2704  109.41  < 2e-16 ***
  income: High vs rest     1   5790    5790  234.29  < 2e-16 ***
  income: Low vs LM & UM   1   1529    1529   61.88 1.93e-13 ***
  income: LM vs UM         1    792     792   32.04 4.93e-08 ***
Residuals                209   5165      25                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

15.7 When assumptions fail

# (a) Welch's ANOVA: unequal variances
oneway.test(life_exp ~ income, data = anova_df, var.equal = FALSE)

    One-way analysis of means (not assuming equal variances)

data:  life_exp and income
F = 91.971, num df = 3.000, denom df = 77.047, p-value < 2.2e-16

# Games-Howell-like follow-up: pairwise Welch t-tests with Holm correction
pairwise.t.test(anova_df$life_exp, anova_df$income, pool.sd = FALSE, p.adjust.method = "holm")

    Pairwise comparisons using t tests with non-pooled SD 

data:  anova_df$life_exp and anova_df$income 

                    Low income Lower middle income Upper middle income
Lower middle income 0.0015     -                   -                  
Upper middle income 2.2e-07    3.4e-07             -                  
High income         2.9e-11    < 2e-16             5.9e-14            

P value adjustment method: holm 

# (b) Kruskal-Wallis: non-parametric, based on ranks
kruskal.test(life_exp ~ income, data = anova_df)

    Kruskal-Wallis rank sum test

data:  life_exp by income
Kruskal-Wallis chi-squared = 138.1, df = 3, p-value < 2.2e-16

pairwise.wilcox.test(anova_df$life_exp, anova_df$income, p.adjust.method = "holm", exact = FALSE)

    Pairwise comparisons using Wilcoxon rank sum test with continuity correction 

data:  anova_df$life_exp and anova_df$income 

                    Low income Lower middle income Upper middle income
Lower middle income 0.00029    -                   -                  
Upper middle income 4.0e-10    8.6e-07             -                  
High income         3.5e-13    < 2e-16             2.2e-12            

P value adjustment method: holm 

15.8 Reporting

“A one-way ANOVA showed that mean life expectancy differed significantly across income groups, F(3, 209) = 109.4, p <0.001, ω² = 0.6. Tukey post-hoc tests indicated that every pair of income groups differed…” (verify against the Tukey table above).

Your turn

Run a one-way ANOVA comparing log_infant across regions (region). Which assumption looks worst? Run Kruskal–Wallis as a check.

16 Two-way ANOVA

Question: Do two categorical factors (and their interaction) affect a continuous outcome?

• Main effect of A, main effect of B
• Interaction A×B: the effect of A depends on the level of B.

Example: Life expectancy by income group and urbanisation level (more vs less urban than the median country).

tw <- d19 |> filter(!is.na(life_exp), !is.na(urban_grp))

table(tw$income, tw$urban_grp)           # cell sizes: UNBALANCED -> use Type II/III SS

                     
                      Less urban More urban
  Low income                  23          2
  Lower middle income         40         11
  Upper middle income         24         30
  High income                 19         64

tw |> group_by(income, urban_grp) |>
  summarise(n = n(), mean = mean(life_exp), sd = sd(life_exp), .groups = "drop") |>
  kable(digits = 2)

income	urban_grp	n	mean	sd
Low income	Less urban	23	61.05	7.28
Low income	More urban	2	71.95	1.38
Lower middle income	Less urban	40	67.00	5.39
Lower middle income	More urban	11	70.70	5.83
Upper middle income	Less urban	24	72.99	4.07
Upper middle income	More urban	30	73.54	4.26
High income	Less urban	19	77.81	4.12
High income	More urban	64	80.19	3.76

ggplot(tw, aes(income, life_exp, colour = urban_grp, group = urban_grp)) +
  stat_summary(fun = mean, geom = "line", linewidth = 1) +
  stat_summary(fun = mean, geom = "point", size = 3) +
  stat_summary(fun.data = mean_se, geom = "errorbar", width = .1) +
  labs(x = NULL, y = "Mean life expectancy (\u00b11 SE)", colour = NULL,
       title = "Interaction plot",
       subtitle = "Non-parallel lines suggest an interaction") +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))

16.1 Fit the model

With unbalanced data the order of terms matters in summary(aov()) (Type I SS). Use Type II (no interaction) or Type III (with interaction, needs sum-to-zero contrasts) via car::Anova().

old <- options(contrasts = c("contr.sum", "contr.poly"))    # required for Type III

fit_tw <- lm(life_exp ~ income * urban_grp, data = tw)
Anova(fit_tw, type = 3)

Anova Table (Type III tests)

Response: life_exp
                 Sum Sq  Df    F value    Pr(>F)    
(Intercept)      412241   1 17825.0106 < 2.2e-16 ***
income             2818   3    40.6228 < 2.2e-16 ***
urban_grp           383   1    16.5420 6.785e-05 ***
income:urban_grp    193   3     2.7807   0.04212 *  
Residuals          4741 205                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

options(old)                                                # restore defaults

16.2 Interpretation guide

1. Look at the interaction first.
2. If the interaction is significant, interpret main effects with caution; describe the simple effects (effect of A within each level of B).
3. If it is not significant, drop it and interpret the main effects from an additive model.

fit_add <- lm(life_exp ~ income + urban_grp, data = tw)
Anova(fit_add, type = 2)

Anova Table (Type II tests)

Response: life_exp
          Sum Sq  Df F value    Pr(>F)    
income    4768.8   3  67.013 < 2.2e-16 ***
urban_grp  230.5   1   9.718  0.002083 ** 
Residuals 4934.0 208                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Does the interaction improve the model? (nested model F-test)
anova(fit_add, fit_tw)

Analysis of Variance Table

Model 1: life_exp ~ income + urban_grp
Model 2: life_exp ~ income * urban_grp
  Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
1    208 4934.0                              
2    205 4741.1  3    192.93 2.7807 0.04212 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Post hoc for the additive model
TukeyHSD(aov(life_exp ~ income + urban_grp, data = tw))

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = life_exp ~ income + urban_grp, data = tw)

$income
                                             diff       lwr       upr   p adj
Lower middle income-Low income           5.876697  2.797017  8.956378 9.4e-06
Upper middle income-Low income          11.371204  8.319793 14.422614 0.0e+00
High income-Low income                  17.719464 14.841689 20.597239 0.0e+00
Upper middle income-Lower middle income  5.494506  3.031494  7.957519 2.0e-07
High income-Lower middle income         11.842767  9.598460 14.087073 0.0e+00
High income-Upper middle income          6.348260  4.142906  8.553614 0.0e+00

$urban_grp
                          diff       lwr    upr     p adj
More urban-Less urban 1.768085 0.4522695 3.0839 0.0086914

16.3 Diagnostics

par(mfrow = c(1, 2)); plot(fit_tw, which = 1:2); par(mfrow = c(1, 1))

leveneTest(life_exp ~ income * urban_grp, data = tw)

Levene's Test for Homogeneity of Variance (center = median)
       Df F value Pr(>F)
group   7  1.2982 0.2527
      205               

Your turn

Add a different second factor (e.g., a median split of health_exp_pc) and repeat. Is there an interaction?

17 Simple linear regression

Question: How does the mean of a continuous outcome Y change with a continuous predictor X?

\[Y_i=\beta_0+\beta_1X_i+\varepsilon_i,\qquad \varepsilon_i\sim N(0,\sigma^2)\]

• \(\beta_0\) intercept (mean Y when X = 0)
• \(\beta_1\) slope (change in mean Y per 1-unit increase in X)
• Estimated by ordinary least squares (OLS): minimise \(\sum (y_i-\hat{y}_i)^2\).

Example: Life expectancy vs log GDP per capita.

17.1 Fit and interpret

slr_df <- d19 |> drop_na(life_exp, log_gdp)

m1 <- lm(life_exp ~ log_gdp, data = slr_df)
summary(m1)

Call:
lm(formula = life_exp ~ log_gdp, data = slr_df)

Residuals:
     Min       1Q   Median       3Q      Max 
-28.8677  -2.3296   0.4004   2.9568   7.5833 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  32.9010     1.9027   17.29   <2e-16 ***
log_gdp       4.5021     0.2099   21.45   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.397 on 205 degrees of freedom
Multiple R-squared:  0.6918,    Adjusted R-squared:  0.6903 
F-statistic: 460.1 on 1 and 205 DF,  p-value: < 2.2e-16

tidy(m1, conf.int = TRUE)

# A tibble: 2 × 7
  term        estimate std.error statistic  p.value conf.low conf.high
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>    <dbl>     <dbl>
1 (Intercept)    32.9      1.90       17.3 6.49e-42    29.1      36.7 
2 log_gdp         4.50     0.210      21.5 2.71e-54     4.09      4.92

glance(m1) |> select(r.squared, adj.r.squared, sigma, statistic, p.value, df, df.residual, AIC)

# A tibble: 1 × 8
  r.squared adj.r.squared sigma statistic  p.value    df df.residual   AIC
      <dbl>         <dbl> <dbl>     <dbl>    <dbl> <dbl>       <int> <dbl>
1     0.692         0.690  4.40      460. 2.71e-54     1         205 1205.

How to read the output

Piece	Meaning
Estimate (log_gdp)	change in mean life expectancy per 1 unit of log GDP
Std. Error, t, p	test of H₀: β₁ = 0
R²	proportion of variability in Y explained by X
Residual standard error (σ̂)	typical prediction error in Y units
F-statistic	overall model test (equals t² in simple regression)

Interpreting a log-transformed predictor: a 1% increase in GDP per capita is associated with about \(\beta_1\times\ln(1.01)\approx\beta_1/100\) years higher life expectancy.

b <- coef(m1)
b["log_gdp"] / 100                          # years per 1% increase in GDP

   log_gdp 
0.04502098 

b["log_gdp"] * log(2)                       # years per DOUBLING of GDP

 log_gdp 
3.120616 

17.2 OLS by hand

x <- slr_df$log_gdp; y <- slr_df$life_exp
b1 <- cor(x, y) * sd(y) / sd(x)             # slope = r * sy / sx
b0 <- mean(y) - b1 * mean(x)                # line passes through (xbar, ybar)
c(b0 = b0, b1 = b1)

       b0        b1 
32.901039  4.502098 

# R-squared = squared correlation in simple regression
cor(x, y)^2

[1] 0.6917828

summary(m1)$r.squared

[1] 0.6917828

# Matrix solution
X <- cbind(1, x)
solve(t(X) %*% X) %*% t(X) %*% y

       [,1]
  32.901039
x  4.502098

17.3 Visualising the fit

ggplot(slr_df, aes(log_gdp, life_exp)) +
  geom_point(aes(colour = income), alpha = .7) +
  geom_smooth(method = "lm", colour = "black") +
  labs(x = "log(GDP per capita)", y = "Life expectancy (years)", colour = "Income group",
       title = "Simple linear regression: life expectancy ~ log GDP")

17.4 Prediction: confidence vs prediction intervals

• Confidence interval: uncertainty about the mean of Y at a given X.
• Prediction interval: uncertainty about a single new observation (always wider).

new <- tibble(gdp = c(1000, 10000, 50000), log_gdp = log(gdp))
bind_cols(new,
          predict(m1, new, interval = "confidence") |> as_tibble() |> rename(ci_lwr = lwr, ci_upr = upr),
          predict(m1, new, interval = "prediction") |> as_tibble() |> select(pi_lwr = lwr, pi_upr = upr))

# A tibble: 3 × 7
    gdp log_gdp   fit ci_lwr ci_upr pi_lwr pi_upr
  <dbl>   <dbl> <dbl>  <dbl>  <dbl>  <dbl>  <dbl>
1  1000    6.91  64.0   63.0   65.0   55.3   72.7
2 10000    9.21  74.4   73.8   75.0   65.7   83.1
3 50000   10.8   81.6   80.6   82.6   72.9   90.3

WarningAvoid extrapolation

Only predict within the observed range of X. Predictions outside it assume the linear pattern continues, which is often false.

17.5 Why transform? Compare raw vs log GDP

m_raw <- lm(life_exp ~ gdp_pc, data = slr_df)
c(R2_raw = summary(m_raw)$r.squared, R2_log = summary(m1)$r.squared)

   R2_raw    R2_log 
0.3650032 0.6917828 

par(mfrow = c(1, 2))
plot(m_raw, which = 1, main = "Raw GDP: curved residuals")
plot(m1,    which = 1, main = "log GDP: better")

par(mfrow = c(1, 1))

17.6 Regression assumptions (LINE)

Letter	Assumption	How to check
L	Linearity of mean Y vs X	scatter, residuals vs fitted (no pattern)
I	Independence of errors	study design (clustering, time, space)
N	Normality of residuals	Q-Q plot, histogram (matters for small n)
E	Equal variance (homoscedasticity)	residuals vs fitted (no funnel), Breusch–Pagan

par(mfrow = c(2, 2)); plot(m1); par(mfrow = c(1, 1))

aug <- augment(m1)

ggplot(aug, aes(.fitted, .resid)) + geom_point(alpha = .6) +
  geom_hline(yintercept = 0, linetype = "dashed") + geom_smooth(se = FALSE) +
  labs(title = "Residuals vs fitted", x = "Fitted", y = "Residual")

ggplot(aug, aes(.resid)) + geom_histogram(bins = 20, fill = "steelblue", colour = "white") +
  labs(title = "Residual distribution")

shapiro.test(aug$.resid)

    Shapiro-Wilk normality test

data:  aug$.resid
W = 0.89095, p-value = 4.072e-11

bptest(m1)                              # H0: constant variance

    studentized Breusch-Pagan test

data:  m1
BP = 7.2563, df = 1, p-value = 0.007065

17.7 Influential observations

infl <- augment(m1) |>
  mutate(country = slr_df$country) |>
  filter(.cooksd > 4 / nrow(slr_df) | abs(.std.resid) > 3) |>
  select(country, life_exp, .fitted, .std.resid, .hat, .cooksd) |>
  arrange(desc(.cooksd))
infl

# A tibble: 7 × 6
  country                  life_exp .fitted .std.resid    .hat .cooksd
  <chr>                       <dbl>   <dbl>      <dbl>   <dbl>   <dbl>
1 Central African Republic     31.5    60.4      -6.64 0.0232   0.524 
2 Nigeria                      53.0    69.2      -3.70 0.00660  0.0454
3 Chad                         53.0    63.5      -2.41 0.0154   0.0453
4 Greenland                    71.4    81.9      -2.40 0.0134   0.0391
5 Lesotho                      55.2    64.4      -2.09 0.0136   0.0300
6 Tajikistan                   71.0    63.4       1.74 0.0156   0.0240
7 Burundi                      62.2    57.5       1.09 0.0326   0.0200

• Outlier: large residual (unusual Y given X).
• Leverage: unusual X (hat value > 2(p+1)/n).
• Influence (Cook’s D): removing the point would change the fit materially. Rule of thumb: D > 4/n deserves a look.

Your turn

Regress life_exp on urban_pct. Interpret slope, R² and the diagnostics. Is the relationship linear?

18 Multiple linear regression

Question: What is the association between X₁ and Y holding other variables constant?

\[Y=\beta_0+\beta_1X_1+\beta_2X_2+\dots+\beta_pX_p+\varepsilon\]

Each \(\beta_j\) = expected change in Y per 1-unit increase in \(X_j\) with all other predictors held fixed (an adjusted association). Uses: adjusting for confounders, improving prediction, comparing relative importance.

18.1 Build a common analysis dataset

All candidate models must use the same rows to be comparable.

reg_df <- d19 |>
  drop_na(life_exp, log_gdp, log_hexp, urban_pct, water_pct, income)
nrow(reg_df)

[1] 181

18.2 Fit models step by step

m1 <- lm(life_exp ~ log_gdp,                                        data = reg_df)
m2 <- lm(life_exp ~ log_gdp + log_hexp,                             data = reg_df)
m3 <- lm(life_exp ~ log_gdp + log_hexp + urban_pct + water_pct,     data = reg_df)
m4 <- lm(life_exp ~ log_gdp + log_hexp + urban_pct + water_pct + income, data = reg_df)

summary(m3)

Call:
lm(formula = life_exp ~ log_gdp + log_hexp + urban_pct + water_pct, 
    data = reg_df)

Residuals:
     Min       1Q   Median       3Q      Max 
-22.7976  -1.7480   0.6424   2.4665   8.0720 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 28.04324    3.27175   8.571 5.08e-15 ***
log_gdp      2.90807    0.79399   3.663  0.00033 ***
log_hexp     0.13077    0.67910   0.193  0.84752    
urban_pct    0.01163    0.01794   0.648  0.51772    
water_pct    0.19809    0.03134   6.322 2.06e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.938 on 176 degrees of freedom
Multiple R-squared:  0.7685,    Adjusted R-squared:  0.7632 
F-statistic:   146 on 4 and 176 DF,  p-value: < 2.2e-16

tidy(m3, conf.int = TRUE) |> kable(digits = 3)

term	estimate	std.error	statistic	p.value	conf.low	conf.high
(Intercept)	28.043	3.272	8.571	0.000	21.586	34.500
log_gdp	2.908	0.794	3.663	0.000	1.341	4.475
log_hexp	0.131	0.679	0.193	0.848	-1.209	1.471
urban_pct	0.012	0.018	0.648	0.518	-0.024	0.047
water_pct	0.198	0.031	6.322	0.000	0.136	0.260

Notice how the log_gdp coefficient changes as predictors are added — evidence of confounding/shared variance (GDP and health spending are correlated).

map_dfr(list(m1 = m1, m2 = m2, m3 = m3, m4 = m4), tidy, .id = "model") |>
  filter(term == "log_gdp") |> select(model, estimate, std.error, p.value) |> kable(digits = 3)

model	estimate	std.error	p.value
m1	4.754	0.226	0
m2	4.031	0.854	0
m3	2.908	0.794	0
m4	3.499	0.908	0

18.3 Adjusted R² and model comparison

\(R^2\) always increases when adding predictors; adjusted R² penalises useless ones. AIC/BIC balance fit and complexity (lower = better).

models <- list(`1: GDP` = m1, `2: + health exp` = m2, `3: + urban, water` = m3, `4: + income` = m4)

map_dfr(models, glance, .id = "model") |>
  select(model, r.squared, adj.r.squared, sigma, AIC, BIC, df) |>
  kable(digits = 3)

model	r.squared	adj.r.squared	sigma	AIC	BIC	df
1: GDP	0.712	0.710	4.356	1050.347	1059.942	1
2: + health exp	0.713	0.710	4.359	1051.564	1064.358	2
3: + urban, water	0.768	0.763	3.938	1016.743	1035.934	4
4: + income	0.777	0.768	3.897	1015.827	1044.614	7

# Nested model F-tests: does the added block improve fit?
anova(m1, m2, m3, m4)

Analysis of Variance Table

Model 1: life_exp ~ log_gdp
Model 2: life_exp ~ log_gdp + log_hexp
Model 3: life_exp ~ log_gdp + log_hexp + urban_pct + water_pct
Model 4: life_exp ~ log_gdp + log_hexp + urban_pct + water_pct + income
  Res.Df    RSS Df Sum of Sq       F    Pr(>F)    
1    179 3396.5                                   
2    178 3381.8  1     14.66  0.9657   0.32713    
3    176 2729.0  2    652.83 21.4984 4.585e-09 ***
4    173 2626.7  3    102.30  2.2458   0.08473 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

18.4 Categorical predictors (dummy coding)

R creates k − 1 dummy variables for a factor with k levels; the first level is the reference category.

levels(reg_df$income)                 # "Low income" is the reference

[1] "Low income"          "Lower middle income" "Upper middle income"
[4] "High income"        

head(model.matrix(m4), 3)             # see the dummy variables

  (Intercept)  log_gdp log_hexp urban_pct water_pct incomeLower middle income
1           1 6.207790 4.304933  25.14373  69.82270                         0
2           1 8.711021 5.903562  56.76391  94.75822                         0
3           1 8.404798 5.489229  72.48246  91.99980                         0
  incomeUpper middle income incomeHigh income
1                         0                 0
2                         1                 0
3                         1                 0

tidy(m4, conf.int = TRUE) |> filter(str_detect(term, "income")) |> kable(digits = 2)

term	estimate	std.error	statistic	p.value	conf.low	conf.high
incomeLower middle income	-3.09	1.37	-2.26	0.03	-5.80	-0.39
incomeUpper middle income	-4.45	1.86	-2.40	0.02	-8.11	-0.79
incomeHigh income	-4.53	2.53	-1.79	0.08	-9.53	0.47

# Change the reference group
m4b <- lm(life_exp ~ log_gdp + log_hexp + urban_pct + water_pct +
            relevel(income, ref = "High income"), data = reg_df)

Each income coefficient = adjusted mean difference in life expectancy versus the reference group.

18.5 Interaction (effect modification)

Does the GDP–life expectancy slope differ by income group?

m5 <- lm(life_exp ~ log_gdp * income + log_hexp + urban_pct + water_pct, data = reg_df)
anova(m4, m5)                         # test the interaction block

Analysis of Variance Table

Model 1: life_exp ~ log_gdp + log_hexp + urban_pct + water_pct + income
Model 2: life_exp ~ log_gdp * income + log_hexp + urban_pct + water_pct
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1    173 2626.7                           
2    170 2606.7  3    19.982 0.4344 0.7287

ggplot(reg_df, aes(log_gdp, life_exp, colour = income)) +
  geom_point(alpha = .5) + geom_smooth(method = "lm", se = FALSE) +
  labs(title = "Different slopes by income group?", x = "log GDP", y = "Life expectancy")

18.6 Centering and standardising

Centering makes the intercept meaningful (value at the average X). Standardising (z-scores) puts coefficients on a common SD scale to compare predictor importance.

reg_c <- reg_df |> mutate(across(c(log_gdp, log_hexp, urban_pct, water_pct), ~ as.numeric(scale(.x)),
                                 .names = "z_{.col}"))

m3_std <- lm(life_exp ~ z_log_gdp + z_log_hexp + z_urban_pct + z_water_pct, data = reg_c)
tidy(m3_std, conf.int = TRUE) |> kable(digits = 3)   # slope = years per 1 SD of predictor

term	estimate	std.error	statistic	p.value	conf.low	conf.high
(Intercept)	72.614	0.293	248.093	0.000	72.036	73.191
z_log_gdp	4.176	1.140	3.663	0.000	1.926	6.427
z_log_hexp	0.215	1.119	0.193	0.848	-1.993	2.424
z_urban_pct	0.263	0.405	0.648	0.518	-0.537	1.063
z_water_pct	2.944	0.466	6.322	0.000	2.025	3.864

18.7 Multicollinearity

When predictors are highly correlated, coefficients become unstable with inflated SEs. Use the Variance Inflation Factor: VIF > 5 (or 10) is a concern.

vif(m3)

  log_gdp  log_hexp urban_pct water_pct 
15.093756 14.535341  1.908117  2.518328 

cor(reg_df |> select(log_gdp, log_hexp, urban_pct, water_pct))

            log_gdp  log_hexp urban_pct water_pct
log_gdp   1.0000000 0.9642504 0.6793559 0.7690261
log_hexp  0.9642504 1.0000000 0.6752707 0.7602783
urban_pct 0.6793559 0.6752707 1.0000000 0.5867053
water_pct 0.7690261 0.7602783 0.5867053 1.0000000

18.8 Non-linearity: polynomial terms

m_quad <- lm(life_exp ~ log_gdp + I(log_gdp^2) + log_hexp + urban_pct + water_pct, data = reg_df)
anova(m3, m_quad)       # does curvature help?

Analysis of Variance Table

Model 1: life_exp ~ log_gdp + log_hexp + urban_pct + water_pct
Model 2: life_exp ~ log_gdp + I(log_gdp^2) + log_hexp + urban_pct + water_pct
  Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
1    176 2729.0                              
2    175 2644.8  1    84.186 5.5704 0.01937 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

18.9 Diagnostics for the multiple regression

par(mfrow = c(2, 2)); plot(m3); par(mfrow = c(1, 1))

# Partial regression (added-variable) plots: each predictor adjusted for the others
avPlots(m3)

# Formal checks
bptest(m3)                       # heteroscedasticity

    studentized Breusch-Pagan test

data:  m3
BP = 27.959, df = 4, p-value = 1.272e-05

shapiro.test(residuals(m3))      # normality of residuals

    Shapiro-Wilk normality test

data:  residuals(m3)
W = 0.90837, p-value = 3.554e-09

# Influence
infl_m3 <- augment(m3) |>
  mutate(country = reg_df$country) |>
  filter(.cooksd > 4 / nrow(reg_df)) |>
  select(country, life_exp, .fitted, .std.resid, .cooksd) |> arrange(desc(.cooksd))
infl_m3

# A tibble: 10 × 5
   country                  life_exp .fitted .std.resid .cooksd
   <chr>                       <dbl>   <dbl>      <dbl>   <dbl>
 1 Central African Republic     31.5    54.3      -6.10  0.814 
 2 Nigeria                      53.0    67.5      -3.76  0.111 
 3 Lesotho                      55.2    64.4      -2.43  0.105 
 4 Tuvalu                       66.1    74.1      -2.09  0.0614
 5 Congo, Dem. Rep.             60.2    54.1       1.64  0.0573
 6 Ethiopia                     65.7    57.7       2.11  0.0528
 7 Eswatini                     59.8    67.9      -2.07  0.0321
 8 Papua New Guinea             65.3    61.1       1.12  0.0285
 9 Uganda                       66.4    59.2       1.86  0.0266
10 Madagascar                   63.5    57.0       1.69  0.0257

# Sensitivity: refit without the influential points
m3_noinf <- update(m3, data = reg_df |> filter(!country %in% infl_m3$country))
bind_rows(all = tidy(m3), without_influential = tidy(m3_noinf), .id = "fit") |>
  select(fit, term, estimate, p.value) |>
  pivot_wider(names_from = fit, values_from = c(estimate, p.value)) |> kable(digits = 3)

term	estimate_all	estimate_without_influential	p.value_all	p.value_without_influential
(Intercept)	28.043	34.740	0.000	0.000
log_gdp	2.908	1.299	0.000	0.058
log_hexp	0.131	1.515	0.848	0.009
urban_pct	0.012	0.018	0.518	0.215
water_pct	0.198	0.188	0.000	0.000

18.10 Robust standard errors (when variance is not constant)

coeftest(m3)                                       # classical SEs

t test of coefficients:

             Estimate Std. Error t value  Pr(>|t|)    
(Intercept) 28.043241   3.271753  8.5713 5.084e-15 ***
log_gdp      2.908075   0.793987  3.6626   0.00033 ***
log_hexp     0.130770   0.679101  0.1926   0.84752    
urban_pct    0.011631   0.017944  0.6482   0.51772    
water_pct    0.198090   0.031335  6.3217 2.063e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

coeftest(m3, vcov = vcovHC(m3, type = "HC3"))      # heteroscedasticity-robust SEs

t test of coefficients:

             Estimate Std. Error t value  Pr(>|t|)    
(Intercept) 28.043241   6.250260  4.4867 1.303e-05 ***
log_gdp      2.908075   1.119213  2.5983  0.010162 *  
log_hexp     0.130770   1.033157  0.1266  0.899423    
urban_pct    0.011631   0.016907  0.6880  0.492389    
water_pct    0.198090   0.059562  3.3258  0.001073 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

18.11 Variable selection and parsimony

full <- lm(life_exp ~ log_gdp + log_hexp + urban_pct + water_pct + income, data = reg_df)
step_model <- step(full, direction = "both", trace = 0)   # AIC-based stepwise
formula(step_model)

life_exp ~ log_gdp + water_pct + income

WarningCaution with automated selection

Stepwise methods produce over-optimistic p-values and unstable models. Prefer theory-driven variable choice (draw a causal diagram: what is the exposure? what are confounders?). Use selection mainly for prediction, and validate (below).

18.12 Model validation: out-of-sample prediction

A model always fits its own data best. Evaluate on unseen data using k-fold cross-validation.

set.seed(2381)
k <- 5
folds <- sample(rep(1:k, length.out = nrow(reg_df)))

cv_rmse <- function(formula) {
  map_dbl(1:k, function(i) {
    train <- reg_df[folds != i, ]
    test  <- reg_df[folds == i, ]
    fit   <- lm(formula, data = train)
    sqrt(mean((test$life_exp - predict(fit, newdata = test))^2))
  }) |> mean()
}

tibble(model = c("1: GDP", "2: + health exp", "3: + urban, water", "4: + income"),
       CV_RMSE = c(cv_rmse(life_exp ~ log_gdp),
                   cv_rmse(life_exp ~ log_gdp + log_hexp),
                   cv_rmse(life_exp ~ log_gdp + log_hexp + urban_pct + water_pct),
                   cv_rmse(life_exp ~ log_gdp + log_hexp + urban_pct + water_pct + income))) |>
  kable(digits = 2)

model	CV_RMSE
1: GDP	4.23
2: + health exp	4.25
3: + urban, water	4.07
4: + income	4.13

18.13 Interpreting and reporting a final model

final <- m3
tidy(final, conf.int = TRUE) |>
  mutate(across(where(is.numeric), ~ round(.x, 3))) |>
  kable(caption = "Final model: life expectancy (years), n = 2019 countries with complete data")

Final model: life expectancy (years), n = 2019 countries with complete data

term	estimate	std.error	statistic	p.value	conf.low	conf.high
(Intercept)	28.043	3.272	8.571	0.000	21.586	34.500
log_gdp	2.908	0.794	3.663	0.000	1.341	4.475
log_hexp	0.131	0.679	0.193	0.848	-1.209	1.471
urban_pct	0.012	0.018	0.648	0.518	-0.024	0.047
water_pct	0.198	0.031	6.322	0.000	0.136	0.260

c(R2 = summary(final)$r.squared, adjR2 = summary(final)$adj.r.squared,
  n = nobs(final), sigma = sigma(final))

         R2       adjR2           n       sigma 
  0.7684521   0.7631897 181.0000000   3.9377130 

Example write-up: “In a multiple linear regression of 2019 life expectancy (n = 181 countries), each doubling of GDP per capita was associated with 2.02 additional years of life expectancy (95% CI 0.93 to 3.1), adjusting for health expenditure, urbanisation and access to basic water. The model explained 76% of the variance (adjusted R²). Results are ecological associations and do not establish causation.”

Your turn

1. Use log_infant as the outcome and fit a model with log_gdp, water_pct and secondary_enrol. Interpret a coefficient on the ratio scale: \(e^{\beta}\).
2. Add fertility. Check VIF. What happens to the other coefficients?

19 Putting it all together

19.1 Which method do I use?

flowchart TD
  A[Continuous outcome] --> B{Question}
  B -->|Compare mean to a value| C[One-sample t-test]
  B -->|Compare 2 groups| D{Independent or paired?}
  B -->|Compare 3+ groups| E{1 factor or 2?}
  B -->|Relationship with continuous X| F{Adjust for others?}
  D -->|Independent| G[Welch t-test / Mann-Whitney]
  D -->|Paired| H[Paired t-test / Wilcoxon signed-rank]
  E -->|1| I[One-way ANOVA + Tukey / Kruskal-Wallis]
  E -->|2| J[Two-way ANOVA with interaction]
  F -->|No| K[Correlation / Simple regression]
  F -->|Yes| L[Multiple linear regression]

19.2 Summary table

Method	Outcome	Predictor / groups	Key assumptions	R function	Effect size
One-sample t	continuous	none (vs μ₀)	normal (or large n)	t.test(x, mu=)	d
Welch t	continuous	2 independent groups	normal within groups	t.test(y ~ g)	d / g
Paired t	continuous	2 paired measures	differences normal	t.test(x, y, paired=TRUE)	dz
Mann–Whitney / Wilcoxon	continuous/ordinal	2 groups	independent, similar shapes	wilcox.test()	r, rank-biserial
One-way ANOVA	continuous	1 factor, k ≥ 3	normal, equal var., independent	aov(), oneway.test()	η², ω²
Kruskal–Wallis	continuous/ordinal	1 factor, k ≥ 3	independent	kruskal.test()	ε²
Two-way ANOVA	continuous	2 factors	as above	lm() + car::Anova()	partial η²
Pearson correlation	2 continuous	—	linear, bivariate normal-ish	cor.test()	r
Spearman	2 continuous/ordinal	—	monotonic	cor.test(method="spearman")	ρ
Simple linear regression	continuous	1 continuous	LINE	lm(y ~ x)	R², β
Multiple regression	continuous	≥ 2 predictors	LINE, no severe collinearity	lm(y ~ x1 + x2)	adj. R², β, partial R²

19.3 A-to-Z analysis checklist

1. Question & design — outcome, exposure/groups, independence, sample size/power.
2. Data — import, clean, check ranges, handle missing data, define variable types.
3. Describe — summaries (mean/SD or median/IQR), plots.
4. Check assumptions — independence, normality, equal variance, linearity, outliers.
5. Choose and run the test/model — match to the question (see flowchart).
6. Diagnose — residuals, influential points, sensitivity analyses.
7. Quantify — effect size with 95% CI, not just p.
8. Interpret in context — practical significance, limitations, confounding, ecological fallacy.
9. Report reproducibly — code, versions, exact n, assumptions checked, software.

Capstone exercises

Use the 2019 World Bank data (or change year == 2019 to another year and compare).

1. Descriptives: Produce a Table 1 of life_exp, infant_mort, fertility, urban_pct by income group (mean (SD) or median (IQR) as appropriate).
2. t-tests: (a) one-sample: is mean fertility different from 2.1 (replacement level)? (b) independent: water_pct in High vs Low-middle income. (c) paired: fertility 2000 vs 2019.
3. ANOVA: Test whether sanitation_pct differs among regions. Check assumptions, report effect size, run suitable post-hoc tests.
4. Correlation: Which variable has the strongest correlation with life_exp? Which has the strongest with infant_mort? Justify Pearson vs Spearman.
5. Regression: Build a model to explain life_exp using at most five predictors. Justify each, check diagnostics, validate with cross-validation, and write a 150-word results paragraph including limitations.
6. Critical thinking: Give two reasons why the observed association between health expenditure and life expectancy could be non-causal.

Appendix: Reproducibility and references

sessionInfo()

R version 4.5.2 (2025-10-31)
Platform: aarch64-apple-darwin20
Running under: macOS Tahoe 26.6.2

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: America/Toronto
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] knitr_1.51      pwr_1.3-0       sandwich_3.1-1  lmtest_0.9-40  
 [5] zoo_1.8-15      car_3.1-3       carData_3.0-5   broom_1.0.11   
 [9] WDI_2.7.9       lubridate_1.9.4 forcats_1.0.1   stringr_1.6.0  
[13] dplyr_1.1.4     purrr_1.2.0     readr_2.1.6     tidyr_1.3.2    
[17] tibble_3.3.0    ggplot2_4.0.1   tidyverse_2.0.0

loaded via a namespace (and not attached):
 [1] gtable_0.3.6       xfun_0.55          htmlwidgets_1.6.4  lattice_0.22-7    
 [5] tzdb_0.5.0         vctrs_0.6.5        tools_4.5.2        generics_0.1.4    
 [9] parallel_4.5.2     pkgconfig_2.0.3    Matrix_1.7-4       RColorBrewer_1.1-3
[13] S7_0.2.1           lifecycle_1.0.4    compiler_4.5.2     farver_2.1.2      
[17] htmltools_0.5.9    yaml_2.3.12        Formula_1.2-5      pillar_1.11.1     
[21] crayon_1.5.3       abind_1.4-8        nlme_3.1-168       tidyselect_1.2.1  
[25] digest_0.6.39      stringi_1.8.7      labeling_0.4.3     splines_4.5.2     
[29] fastmap_1.2.0      grid_4.5.2         archive_1.1.12.1   cli_3.6.5         
[33] magrittr_2.0.4     dichromat_2.0-0.1  utf8_1.2.6         withr_3.0.2       
[37] scales_1.4.0       backports_1.5.0    bit64_4.6.0-1      timechange_0.3.0  
[41] rmarkdown_2.30     bit_4.6.0          otel_0.2.0         hms_1.1.4         
[45] evaluate_1.0.5     mgcv_1.9-3         rlang_1.3.0        glue_1.8.0        
[49] rstudioapi_0.17.1  vroom_1.6.7        jsonlite_2.0.0     R6_2.6.1          

Data source: World Bank, World Development Indicators, https://data.worldbank.org, licensed CC BY 4.0. Accessed via the WDI R package.

Suggested reading:

• Field, A. Discovering Statistics Using R.
• Dalgaard, P. Introductory Statistics with R.
• Vittinghoff, E. et al. Regression Methods in Biostatistics.
• Wickham, H. & Grolemund, G. R for Data Science (https://r4ds.hadley.nz).
• Quarto documentation: https://quarto.org.