needed <- c("knitr", "jsonlite", "curl")
missing_packages <- needed[!vapply(needed, requireNamespace,
quietly = TRUE, FUN.VALUE = logical(1))]
if (length(missing_packages)) {
install.packages(missing_packages, repos = "https://cloud.r-project.org")
}
if (!all(vapply(needed, requireNamespace, quietly = TRUE,
FUN.VALUE = logical(1)))) {
stop("Install knitr, jsonlite and curl, then render again.")
}
knitr::opts_chunk$set(fig.width = 8, fig.height = 5, dpi = 120)
options(stringsAsFactors = FALSE, scipen = 5, timeout = 180)
set.seed(params$seed)
analysis_year <- as.integer(params$analysis_year)
baseline_year <- as.integer(params$baseline_year)
B <- as.integer(params$bootstrap_reps)
P <- as.integer(params$permutation_reps)
stopifnot(baseline_year < analysis_year, B >= 200, P >= 199)
dir.create("hss2381_cache", showWarnings = FALSE)
dir.create("hss2381_outputs", showWarnings = FALSE)
fmt <- function(x, digits = 2) formatC(x, digits = digits, format = "f")
pfmt <- function(p) format.pval(p, digits = 3, eps = 0.001)
tab <- function(x, caption = NULL, digits = 3) {
knitr::kable(x, caption = caption, digits = digits, row.names = FALSE)
}HSS 2381: Quantitative Methods in Health Sciences
Continuous variables from A to Z: a complete R workbook using World Bank open data
Welcome and scope
Health sciences students need more than a list of statistical tests. They need to understand what a variable means, how a dataset was made, what a graph shows, which comparison answers a question, and how much confidence a conclusion deserves. This workbook builds those skills in order, using one connected set of real, freely accessible World Bank indicators.
Every substantive dataset in the analyses comes from the World Bank Indicators API. There are no invented countries, fabricated health measurements or synthetic patient records. Bootstrap samples, permutations and teaching subsamples reuse the downloaded observations; they do not create new observed data. Numerical constants such as a 70-year reference mean or a two-year planning difference are explicitly hypothetical teaching choices, not World Bank estimates.
ImportantWhat these data can and cannot teach
The unit is a country or economy, not a patient. Country-level life expectancy is a population indicator, not an individual’s observed lifespan. Countries are not a simple random sample of people, and available economies are not necessarily a probability sample of all economies. Formal p-values and confidence intervals in this workbook are model-based teaching illustrations, conditional on assumptions about an underlying country-level process. They do not automatically justify population sampling claims, individual-level conclusions or causal effects. Geographic dependence, shared history and measurement uncertainty can invalidate simple independent-observation models.
0.1 Learning outcomes
By the end, you should be able to:
- Distinguish variables, observations, parameters, statistics, populations and samples.
- Import, audit, reshape and document real open data reproducibly.
- Describe centre, spread, shape, missingness and unusual observations.
- Use probability, normal and t distributions, standardisation and sampling distributions.
- Explain standard errors, confidence intervals, hypotheses, p-values, errors and power.
- Conduct and interpret one-sample, independent and paired t-tests.
- Conduct one-way ANOVA, Welch ANOVA, post-hoc comparisons and rank-based alternatives.
- Measure association with Pearson and Spearman correlations.
- Fit, interpret, diagnose and compare simple and multiple linear regression models.
- Distinguish statistical evidence, practical importance, prediction and causation.
- Report estimates, uncertainty, assumptions, sample sizes and limitations together.
- Complete a reproducible mini-project with tables, figures and a transparent analysis trail.
0.2 Suggested teaching sequence
| Block | Main topics | Suggested practical task |
|---|---|---|
| 1 | Research questions, measurement, R and provenance | Download and audit the indicators |
| 2 | Data preparation, missingness and descriptive statistics | Build a descriptive table |
| 3 | Visualisation, probability and distributions | Explain a histogram and a normal approximation |
| 4 | Sampling, standard errors and confidence intervals | Bootstrap a country-level mean |
| 5 | Hypothesis testing and one-sample t-tests | Compare with a stated reference |
| 6 | Independent and paired comparisons | Compare income groups and matched years |
| 7 | ANOVA and multiple comparisons | Compare four income groups |
| 8 | Non-parametric methods and correlation | Contrast rank and mean-based analyses |
| 9 | Simple linear regression | Model life expectancy and log expenditure |
| 10 | Multiple regression and diagnostics | Adjust for GDP and urbanisation |
| 11 | Extensions, sensitivity and prediction | Evaluate a held-out prediction model |
| 12 | Reporting and synthesis | Complete the capstone and explain limitations |
This sequence is a teaching suggestion, not a statement about the university’s timetable.
1 Getting started with R and Quarto
1.1 Software and rendering
Validation note: the document structure, YAML, chunk labels and code delimiter balance were checked during preparation. R and Quarto were not installed in the preparation environment, so the analyses have not been executed or render-tested here. The World Bank API also could not be downloaded from that environment; no data snapshot or numerical results are claimed. The code downloads the real data on your first successful render and retains a cache.
Install R, a current RStudio or another suitable editor, and Quarto. Save this file in a project folder. The first render needs internet access to download packages if absent and to obtain World Bank data. Later renders use a local RDS cache unless refresh_data: true is set in the YAML. The rendered HTML embeds its visual resources, but the QMD source itself is not an offline data snapshot.
Run the following once in the R console if your computer does not permit package installation during rendering:
install.packages(c("rmarkdown", "knitr", "jsonlite", "curl"),
repos = "https://cloud.r-project.org")Then use Render in RStudio, or run this in a terminal:
quarto render HSS2381_Continuous_Variables_World_Bank.qmd --to htmlThe Quarto R guide explains executable chunks and rendering. All executable chunks below are intended to run in document order in a fresh R session. A later chunk can depend on objects created earlier. Do not begin in the middle without running the setup and data preparation.
1.2 Reproducibility settings
A seed makes resampling reproducible for a fixed R version, dataset and code. It does not stop the World Bank from revising historical estimates. Keeping the cached data and session information is therefore important.
1.3 R from the very beginning
R stores values in named objects. <- assigns a value. A function performs an operation. # starts a comment. Names are case-sensitive: life_exp and Life_Exp are different objects. NA means missing, not zero.
course_code <- "HSS 2381" # Character string
reference_mean <- 70 # Hypothetical teaching reference
is_continuous_course <- TRUE # Logical value
reference_mean + 2[1] 72sqrt(reference_mean)[1] 8.3666c(analysis_year, baseline_year) # Combine values into a vector[1] 2023 2010class(reference_mean)[1] "numeric"A data frame is a rectangular table. $ selects a column, square brackets select rows and columns, and a logical condition identifies rows. mean(x, na.rm = TRUE) calculates a mean after excluding missing values; it does not explain why those values are missing. We will use base R for most analyses so the statistical operations remain visible and package requirements stay small.
2 Research questions, measurement and study design
2.1 Start with the question, not the test
Our connected questions are:
- What is the distribution of country-level life expectancy in a fixed year?
- How different are the means in low-income and high-income economies?
- How did life expectancy change between two years within the same economies?
- Do mean life expectancies differ across four income classifications?
- How is life expectancy associated with health expenditure, GDP and urbanisation?
- How well can a model predict life expectancy for held-out economies?
These are descriptive, associational or predictive questions. None becomes causal merely because a regression coefficient has a small p-value.
Before choosing a method, write down the outcome, predictor or grouping variable, unit, time, target quantity and design. For example: outcome = life expectancy in years; predictor = health expenditure per capita; unit = economy; time = 2023; target = change in expected life expectancy associated with a doubling of expenditure; design = ecological cross-sectional observational analysis.
2.2 Measurement scales
| Type | Meaning | Example in this workbook | Common mistake |
|---|---|---|---|
| Nominal | Categories without numerical order | World Bank region | Averaging category codes |
| Ordinal | Ordered categories | Income classification | Treating steps as equally spaced amounts |
| Continuous | Values on a numerical continuum | Life expectancy, GDP per capita | Assuming a rounded value is a count |
| Discrete count | Whole-number counts | Population | Treating a count as a proportion |
| Percentage | A proportion multiplied by 100 | Urban population percentage | Confusing percentage points with percentage change |
Income group is a categorical predictor even though the outcome is continuous. A continuous-variable course can therefore use group labels in t-tests and ANOVA. We do not turn life expectancy into a binary outcome here, because that discards information and changes the question.
Life expectancy and expenditure have meaningful units. A five-year difference is interpretable directly. A standardised difference is unitless and useful for some comparisons, but it does not replace the original units. GDP and expenditure in current US dollars are affected by prices and exchange rates. Cross-year financial comparisons require more care than a single-year comparison; this workbook uses paired years for life expectancy, not for nominal expenditure growth.
2.3 Design, bias and independence
A cross-sectional dataset has one row per economy in a selected year. A longitudinal dataset has repeated rows for the same economy. The paired t-test respects matching between years. A simple regression on all country-year rows would not automatically respect that dependence.
Selection bias can arise because countries with missing indicators differ from those with complete data. Confounding occurs when another variable helps explain an observed association. Measurement error can occur even in official indicators, which may include estimates or harmonised series. These problems cannot be diagnosed from a p-value alone.
3 Downloading and documenting real World Bank data
3.1 Indicators and sources
The World Bank Indicators API documentation describes the version 2 API and notes that authentication is not required. The country metadata documentation describes country identifiers, region and income fields. The aggregate documentation explains that the API also contains regional and income aggregates. We must remove those aggregate rows before analysing economies.
| R name | World Bank code | Indicator and unit |
|---|---|---|
life_exp |
SP.DYN.LE00.IN |
Life expectancy at birth, total, years |
health_pc |
SH.XPD.CHEX.PC.CD |
Current health expenditure per capita, current US dollars |
gdp_pc |
NY.GDP.PCAP.CD |
GDP per capita, current US dollars |
urban_pct |
SP.URB.TOTL.IN.ZS |
Urban population, percentage of total |
population |
SP.POP.TOTL |
Total population, people |
Indicator pages: life expectancy, health expenditure, GDP per capita, urban population, and total population. Read the metadata before interpreting any indicator. The World Bank republishes series from contributing organisations; the health expenditure page identifies WHO as its source. Preserve that attribution.
WarningA fixed year is not a latest-value dataset
The default analysis year is 2023, chosen as a teaching setting rather than a claim that it is the latest year. We request that exact year. Missing 2023 values stay missing. We never silently substitute a country’s most recent available year. Current metadata income classifications are retrieved at download time and are not historical 2023 classifications. All income-group examples use those current labels explicitly.
3.2 A transparent downloader with pagination, retries and caching
The next chunk downloads metadata and each indicator separately. It checks pagination rather than assuming the first page is complete. It records request URLs and retrieval time. If the API is unavailable and no valid cache exists, rendering stops with an informative error rather than fabricating data.
indicators <- c(life_exp = "SP.DYN.LE00.IN",
health_pc = "SH.XPD.CHEX.PC.CD",
gdp_pc = "NY.GDP.PCAP.CD",
urban_pct = "SP.URB.TOTL.IN.ZS",
population = "SP.POP.TOTL")
fetch_json <- function(url) {
last_error <- "Unknown download error"
for (attempt in seq_len(3L)) {
ans <- tryCatch({
h <- curl::new_handle(useragent = "HSS2381 educational R workbook",
connecttimeout = 30, timeout = 180)
response <- curl::curl_fetch_memory(url, handle = h)
if (response$status_code != 200L) {
stop("HTTP status ", response$status_code)
}
jsonlite::fromJSON(rawToChar(response$content), simplifyVector = FALSE)
}, error = function(e) e)
if (!inherits(ans, "error")) return(ans)
last_error <- conditionMessage(ans)
if (attempt < 3L) Sys.sleep(attempt)
}
stop("World Bank download failed: ", last_error,
"\nRequest: ", url,
"\nCheck internet/API availability or retain a previous valid cache.")
}
fetch_pages <- function(endpoint, query = "") {
records <- list()
urls <- character()
page <- 1L
repeat {
url <- paste0("https://api.worldbank.org/v2/", endpoint,
"?format=json&per_page=10000&page=", page, query)
obj <- fetch_json(url)
if (length(obj) != 2L || is.null(obj[[1]]$pages) ||
!is.list(obj[[2]])) {
stop("Unexpected or empty API response for ", endpoint)
}
records <- c(records, obj[[2]])
urls <- c(urls, url)
if (page >= as.integer(obj[[1]]$pages)) break
page <- page + 1L
}
list(records = records, urls = urls)
}
scalar <- function(x, default = NA_character_) {
if (is.null(x) || length(x) == 0L) default else as.character(x[[1]])
}
cache_file <- file.path("hss2381_cache",
paste0("wb_v1_", baseline_year, "_", analysis_year, ".rds"))
if (file.exists(cache_file) && !isTRUE(params$refresh_data)) {
bundle <- readRDS(cache_file)
} else {
meta_response <- fetch_pages("country")
metadata <- do.call(rbind, lapply(meta_response$records, function(z) {
data.frame(iso3c = scalar(z$id), country = scalar(z$name),
region = scalar(z$region$value),
region_id = scalar(z$region$id),
income_id = scalar(z$incomeLevel$id),
income = scalar(z$incomeLevel$value))
}))
# The metadata identifies aggregates with region_id == "NA".
metadata <- metadata[!is.na(metadata$region_id) &
metadata$region_id != "NA" &
metadata$region != "Aggregates", ]
requested_urls <- meta_response$urls
series <- vector("list", length(indicators))
names(series) <- names(indicators)
for (nm in names(indicators)) {
response <- fetch_pages(
paste0("country/all/indicator/", indicators[[nm]]),
paste0("&source=2&date=", baseline_year, ":", analysis_year))
series[[nm]] <- do.call(rbind, lapply(response$records, function(z) {
data.frame(iso3c = scalar(z$countryiso3code),
year = as.integer(scalar(z$date)),
value = as.numeric(scalar(z$value)))
}))
names(series[[nm]])[3] <- nm
requested_urls <- c(requested_urls, response$urls)
}
panel <- Reduce(function(a, b) merge(a, b, by = c("iso3c", "year"),
all = TRUE), series)
panel <- merge(panel, metadata, by = "iso3c", all = FALSE)
panel <- panel[order(panel$iso3c, panel$year), ]
rownames(panel) <- NULL
bundle <- list(panel = panel, metadata = metadata,
retrieved_utc = format(Sys.time(), tz = "UTC", usetz = TRUE),
indicators = indicators, requested_urls = requested_urls,
analysis_year = analysis_year, baseline_year = baseline_year)
saveRDS(bundle, cache_file)
}
stopifnot(is.list(bundle), is.data.frame(bundle$panel),
identical(bundle$indicators, indicators),
bundle$analysis_year == analysis_year,
bundle$baseline_year == baseline_year)
panel <- bundle$panel
stopifnot(!anyDuplicated(panel[c("iso3c", "year")]))
cat("Retrieved (UTC):", bundle$retrieved_utc, "\n")Retrieved (UTC): 2026-10-17 12:45:09 UTC cat("Country-year rows:", nrow(panel), "\n")Country-year rows: 3038 cat("Distinct economies:", length(unique(panel$iso3c)), "\n")Distinct economies: 217 The cache includes the metadata and all requested series. It is not a permanent guarantee of historical classification accuracy. Keep a copy with any submitted analysis. If you change indicator definitions in the source, use a new cache version or refresh deliberately.
3.3 Audit the imported values
str(panel)'data.frame': 3038 obs. of 12 variables:
$ iso3c : chr "ABW" "ABW" "ABW" "ABW" ...
$ year : int 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 ...
$ life_exp : num 74.3 74.6 74.8 75.1 75.3 ...
$ health_pc : num NA NA NA NA NA NA NA NA NA NA ...
$ gdp_pc : num 24093 25712 25120 25814 26130 ...
$ urban_pct : num 64.2 64 63.8 63.6 63.4 ...
$ population: num 101838 102591 104110 105675 106807 ...
$ country : chr "Aruba" "Aruba" "Aruba" "Aruba" ...
$ region : chr "Latin America & Caribbean " "Latin America & Caribbean " "Latin America & Caribbean " "Latin America & Caribbean " ...
$ region_id : chr "LCN" "LCN" "LCN" "LCN" ...
$ income_id : chr "HIC" "HIC" "HIC" "HIC" ...
$ income : chr "High income" "High income" "High income" "High income" ...head(panel[c("country", "iso3c", "year", names(indicators))])stopifnot(all(is.na(panel$life_exp) |
(panel$life_exp > 0 & panel$life_exp < 130)))
stopifnot(all(is.na(panel$urban_pct) |
(panel$urban_pct >= 0 & panel$urban_pct <= 100)))
stopifnot(all(is.na(panel$health_pc) | panel$health_pc >= 0))
stopifnot(all(is.na(panel$gdp_pc) | panel$gdp_pc >= 0))
stopifnot(all(is.na(panel$population) | panel$population > 0))These are broad plausibility checks, not a substitute for metadata review. A failed check should trigger investigation. Do not automatically delete a value merely because it is unusual.
4 Preparing analysis datasets
4.1 One economy, one selected year
cross <- panel[panel$year == analysis_year, ]
stopifnot(nrow(cross) > 0L, !anyDuplicated(cross$iso3c))
# Use stable metadata IDs rather than assuming English label spellings.
cross$income_group <- factor(cross$income_id,
levels = c("LIC", "LMC", "UMC", "HIC"),
labels = c("Low income", "Lower middle income",
"Upper middle income", "High income"))
cross$log2_health <- ifelse(!is.na(cross$health_pc) & cross$health_pc > 0,
log2(pmax(cross$health_pc, .Machine$double.xmin)), NA_real_)
cross$log2_gdp <- ifelse(!is.na(cross$gdp_pc) & cross$gdp_pc > 0,
log2(pmax(cross$gdp_pc, .Machine$double.xmin)), NA_real_)
cross$urban10 <- cross$urban_pct / 10
cross$high_income <- factor(ifelse(is.na(cross$income_group), NA,
ifelse(cross$income_id == "HIC", "High", "Other")),
levels = c("Other", "High"))
life <- cross[is.finite(cross$life_exp), ]
stopifnot(nrow(life) >= 30)
cat("Selected year:", analysis_year, "\n")Selected year: 2023 cat("Economies with observed life expectancy:", nrow(life), "\n")Economies with observed life expectancy: 217 tab(as.data.frame(table(cross$income_group, useNA = "ifany")),
"Current metadata income groups in the selected-year dataset")| Var1 | Freq |
|---|---|
| Low income | 25 |
| Lower middle income | 47 |
| Upper middle income | 59 |
| High income | 86 |
A base-2 logarithm makes doubling easy to interpret. If expenditure rises from 100 to 200, or from 1,000 to 2,000, its log base 2 increases by one. Zero or negative values cannot enter this log transformation. We mark them missing rather than adding an arbitrary constant. All transformations here are explicit analytical choices and are documented before modelling.
4.2 Missingness is part of the result
missing_table <- data.frame(
variable = names(indicators),
observed = vapply(cross[names(indicators)], function(x) sum(!is.na(x)), integer(1)),
missing = vapply(cross[names(indicators)], function(x) sum(is.na(x)), integer(1)))
missing_table$missing_percent <- 100 * missing_table$missing / nrow(cross)
tab(missing_table, "Indicator availability in the fixed analysis year")| variable | observed | missing | missing_percent |
|---|---|---|---|
| life_exp | 217 | 0 | 0.000 |
| health_pc | 191 | 26 | 11.982 |
| gdp_pc | 204 | 13 | 5.991 |
| urban_pct | 217 | 0 | 0.000 |
| population | 217 | 0 | 0.000 |
missing_by_income <- aggregate(
cbind(life_missing = as.numeric(is.na(cross$life_exp)),
expenditure_missing = as.numeric(is.na(cross$health_pc))) ~ income_group,
data = cross, FUN = mean)
tab(missing_by_income, "Proportion missing within current income groups")| income_group | life_missing | expenditure_missing |
|---|---|---|
| Low income | 0 | 0.040 |
| Lower middle income | 0 | 0.000 |
| Upper middle income | 0 | 0.051 |
| High income | 0 | 0.256 |
MCAR means missingness is unrelated to observed or unobserved values. MAR means missingness can depend on observed variables after conditioning appropriately. MNAR means it still depends on unobserved values. These are assumptions about a data-generating process, not labels that a simple missingness table proves.
Complete-case analysis is convenient but can change the population represented. Mean imputation artificially reduces variation and can distort relationships. This workbook does not impute missing values. Multiple imputation is an extension requiring a defensible imputation model, not a button to press automatically.
4.3 Build separate datasets for separate questions
independent <- droplevels(cross[
complete.cases(cross[c("life_exp", "income_group")]) &
cross$income_id %in% c("LIC", "HIC"), ])
# Direction is explicitly Low minus High throughout the independent comparison.
independent$income_group <- factor(independent$income_group,
levels = c("Low income", "High income"))
anova_data <- droplevels(cross[
complete.cases(cross[c("life_exp", "income_group")]), ])
reg_vars <- c("life_exp", "log2_health", "log2_gdp", "urban10")
reg <- cross[complete.cases(cross[reg_vars]), ]
reg <- reg[apply(reg[reg_vars], 1, function(z) all(is.finite(z))), ]
rownames(reg) <- NULL
stopifnot(all(table(independent$income_group) >= 3),
nlevels(anova_data$income_group) == 4,
all(table(anova_data$income_group) >= 3), nrow(reg) >= 30)
before <- panel[panel$year == baseline_year, c("iso3c", "country", "life_exp")]
after <- panel[panel$year == analysis_year, c("iso3c", "life_exp")]
names(before)[3] <- "before"
names(after)[2] <- "after"
paired <- merge(before, after, by = "iso3c", all = FALSE)
paired <- paired[complete.cases(paired[c("before", "after")]), ]
paired$change <- paired$after - paired$before
stopifnot(nrow(paired) >= 10, !anyDuplicated(paired$iso3c))
tab(data.frame(analysis = c("Life expectancy description", "Low vs high income",
"Four-group ANOVA", "Multiple regression", "Paired years"),
n_economies = c(nrow(life), nrow(independent), nrow(anova_data),
nrow(reg), nrow(paired))),
"Analysis-specific denominators")| analysis | n_economies |
|---|---|
| Life expectancy description | 217 |
| Low vs high income | 111 |
| Four-group ANOVA | 217 |
| Multiple regression | 187 |
| Paired years | 217 |
Matching is by identifier, never by row position. We do not assume the first country in one year is the first country in another. Different analyses legitimately use different denominators, but the report must say so.
5 Descriptive statistics: centre, spread and shape
5.1 Population quantities and sample summaries
A parameter describes a target population or model. A statistic is calculated from observed data. The arithmetic mean is
\[\bar{x}=\frac{1}{n}\sum_{i=1}^{n}x_i.\]
The sample variance and standard deviation are
\[s^2=\frac{\sum_i(x_i-\bar{x})^2}{n-1},\qquad s=\sqrt{s^2}.\]
The denominator \(n-1\) reflects estimation of the mean from the same observations in the usual sample-variance calculation. The SD has the original measurement unit. Variance has squared units. Neither is the same as a standard error.
The median is the middle ordered value. The interquartile range is \(Q_3-Q_1\) and describes the middle half of observations. A range uses only the minimum and maximum and is highly sensitive to extremes. A mean can be useful with skewed data, but it should not be presented without information about shape and spread.
summarise_numeric <- function(x) {
x <- x[is.finite(x)]
if (length(x) < 2) return(c(n = length(x), mean = NA, sd = NA,
median = NA, q1 = NA, q3 = NA, min = NA, max = NA))
c(n = length(x), mean = mean(x), sd = sd(x), median = median(x),
q1 = unname(quantile(x, .25)), q3 = unname(quantile(x, .75)),
min = min(x), max = max(x))
}
descriptives <- do.call(rbind, lapply(names(indicators), function(nm) {
data.frame(variable = nm, as.list(summarise_numeric(cross[[nm]])))
}))
tab(descriptives, "Descriptive statistics: each economy receives equal weight")| variable | n | mean | sd | median | q1 | q3 | min | max |
|---|---|---|---|---|---|---|---|---|
| life_exp | 217 | 73.791 | 7.123 | 74.552 | 68.484 | 79.594 | 54.462 | 86.372 |
| health_pc | 191 | 1433.837 | 2328.267 | 470.403 | 94.383 | 1611.698 | 17.440 | 13473.193 |
| gdp_pc | 204 | 22295.539 | 33559.351 | 7862.027 | 2945.927 | 30638.403 | 249.795 | 256799.911 |
| urban_pct | 217 | 62.610 | 23.243 | 63.663 | 43.250 | 81.165 | 14.678 | 100.000 |
| population | 217 | 37048618.129 | 142175383.814 | 6623183.000 | 826353.000 | 26418204.000 | 9816.000 | 1438069596.000 |
5.2 Verify the mean and SD by hand
x <- life$life_exp
n <- length(x)
manual_mean <- sum(x) / n
manual_variance <- sum((x - manual_mean)^2) / (n - 1)
manual_sd <- sqrt(manual_variance)
tab(data.frame(quantity = c("Mean", "Variance", "SD"),
manual = c(manual_mean, manual_variance, manual_sd),
built_in = c(mean(x), var(x), sd(x))))| quantity | manual | built_in |
|---|---|---|
| Mean | 73.791 | 73.791 |
| Variance | 50.735 | 50.735 |
| SD | 7.123 | 7.123 |
stopifnot(isTRUE(all.equal(manual_mean, mean(x))),
isTRUE(all.equal(manual_sd, sd(x))))5.3 Country averages versus population weighting
weighted_data <- cross[complete.cases(cross[c("life_exp", "population")]), ]
tab(data.frame(
summary = c("Unweighted mean across the same economies",
"Population-weighted mean of economy indicators"),
value = c(mean(weighted_data$life_exp),
weighted.mean(weighted_data$life_exp, weighted_data$population))))| summary | value |
|---|---|
| Unweighted mean across the same economies | 73.791 |
| Population-weighted mean of economy indicators | 73.310 |
The unweighted mean answers a question about an average economy. The population-weighted mean gives larger populations more influence. It is not automatically global life expectancy: life expectancy is derived from mortality schedules, and averaging national life expectancies does not reconstruct a global life table. Population weights are also not survey sampling weights. Do not use lm(weights = population) and claim that it solves the ecological inference problem.
5.4 Group summaries
group_summary <- do.call(rbind, lapply(levels(anova_data$income_group), function(g) {
z <- anova_data$life_exp[anova_data$income_group == g]
data.frame(group = g, as.list(summarise_numeric(z)))
}))
tab(group_summary, "Life expectancy by current metadata income classification")| group | n | mean | sd | median | q1 | q3 | min | max |
|---|---|---|---|---|---|---|---|---|
| Low income | 25 | 64.202 | 4.601 | 63.651 | 61.183 | 67.353 | 55.069 | 73.642 |
| Lower middle income | 47 | 67.910 | 4.896 | 67.649 | 65.053 | 71.712 | 54.462 | 77.817 |
| Upper middle income | 59 | 73.770 | 3.778 | 74.184 | 71.354 | 77.274 | 63.707 | 81.041 |
| High income | 86 | 79.807 | 4.177 | 81.134 | 77.367 | 82.941 | 62.109 | 86.372 |
Read the output: compare the mean and median within each group, the SD across groups, and the group sizes. A small group with a large SD has less precise mean estimation than a large group with a small SD, all else equal. Group differences remain country-level descriptions.
6 Visualisation before testing
6.1 Histograms and density curves
A histogram shows frequencies in intervals. Its appearance depends on bin boundaries and width. A density curve smooths the distribution and has total area one; its height is not a probability at an exact value.
oldpar <- par(mfrow = c(1, 2))
hist(life$life_exp, breaks = "FD", main = "Life expectancy",
xlab = "Years", ylab = "Number of economies")
hist(cross$health_pc[is.finite(cross$health_pc)], breaks = "FD",
main = "Health expenditure", xlab = "Current US dollars per capita",
ylab = "Number of economies")
par(oldpar)
z <- cross$log2_health[is.finite(cross$log2_health)]
hist(z, probability = TRUE, breaks = "FD", main = "Log expenditure",
xlab = "Log2 current health expenditure per capita", ylab = "Density")
lines(density(z), lwd = 2)
A log transformation can reduce right skew, but it changes the scale and therefore the interpretation. It does not guarantee normal residuals or make an observational analysis causal.
6.2 Boxplots and raw observations
boxplot(life_exp ~ income_group, data = anova_data, outline = FALSE,
xlab = "Current metadata income group", ylab = "Life expectancy (years)",
names = c("Low", "Lower middle", "Upper middle", "High"))
stripchart(life_exp ~ income_group, data = anova_data, vertical = TRUE,
method = "jitter", jitter = .12, add = TRUE, pch = 16, cex = .5)
The box spans the interquartile range. A conventional whisker rule flags observations beyond 1.5 IQR from a quartile. A flagged observation is not necessarily an error. The raw points reveal unequal sample sizes and patterns that a box alone can hide.
6.3 Scatterplots and a useful first look at association
plot(reg$log2_health, reg$life_exp, pch = 16,
xlab = "Log2 health expenditure per capita",
ylab = "Life expectancy (years)")
abline(lm(life_exp ~ log2_health, data = reg), lwd = 2)
Look for direction, curvature, clusters, changing spread and influential points. A single correlation coefficient cannot describe all of those features.
7 Standardisation and unusual observations
A z-score expresses distance from a mean in SD units:
\[z_i=\frac{x_i-\bar{x}}{s}.\]
A z-score of 1 means one SD above the mean, not one year above it. Standardising a variable does not make its distribution normal.
life$z_life <- as.numeric(scale(life$life_exp))
quartiles <- quantile(life$life_exp, c(.25, .75))
iqr <- diff(quartiles)
life$iqr_flag <- life$life_exp < quartiles[1] - 1.5 * iqr |
life$life_exp > quartiles[2] + 1.5 * iqr
tab(life[order(abs(life$z_life), decreasing = TRUE),
c("country", "life_exp", "z_life", "iqr_flag")][1:10, ],
"Largest absolute standardised deviations: investigate, do not automatically delete")| country | life_exp | z_life | iqr_flag |
|---|---|---|---|
| Nigeria | 54.462 | -2.714 | FALSE |
| Chad | 55.069 | -2.628 | FALSE |
| Lesotho | 57.375 | -2.305 | FALSE |
| Central African Republic | 57.408 | -2.300 | FALSE |
| South Sudan | 57.617 | -2.271 | FALSE |
| Somalia, Fed. Rep. | 58.816 | -2.102 | FALSE |
| Mali | 60.439 | -1.875 | FALSE |
| Guinea | 60.740 | -1.832 | FALSE |
| Benin | 60.774 | -1.827 | FALSE |
| Burkina Faso | 61.092 | -1.783 | FALSE |
An extreme outcome, a high-leverage predictor value and an influential regression observation are different concepts. We return to this distinction in the diagnostics section.
8 Probability and continuous distributions
8.1 Probability rules in a real dataset
An empirical probability can be a proportion of observed economies satisfying a condition. It describes the observed set unless a sampling or superpopulation interpretation is justified.
prob_over75 <- mean(life$life_exp > 75)
prob_65to80 <- mean(life$life_exp >= 65 & life$life_exp <= 80)
prob_high <- mean(anova_data$income_group == "High income")
prob_over75_given_high <- with(anova_data,
mean(life_exp[income_group == "High income"] > 75))
tab(data.frame(quantity = c("P(life expectancy > 75)",
"P(65 <= life expectancy <= 80)",
"P(current high-income classification)",
"P(life expectancy > 75 | current high income)"),
empirical_proportion = c(prob_over75, prob_65to80, prob_high,
prob_over75_given_high)))| quantity | empirical_proportion |
|---|---|
| P(life expectancy > 75) | 0.484 |
| P(65 <= life expectancy <= 80) | 0.631 |
| P(current high-income classification) | 0.396 |
| P(life expectancy > 75 | current high income) | 0.895 |
For events \(A\) and \(B\), \(P(A\mid B)=P(A\cap B)/P(B)\) when \(P(B)>0\). Independence means \(P(A\cap B)=P(A)P(B)\). Conditional and marginal probabilities are different. A difference between them is an association, not proof of a causal mechanism.
For a genuinely continuous random variable, probability is area under a density. The probability of one exact real-valued point is zero in the mathematical model, even though recorded data may be rounded.
8.2 The normal distribution
A normal model is described by mean \(\mu\) and SD \(\sigma\). It is symmetric and bell-shaped. In R, dnorm() gives density, pnorm() cumulative probability, and qnorm() a quantile. We use those functions with parameters estimated from the observed life expectancy distribution as an approximation, not as a claim that country life expectancy is exactly normal.
mu_hat <- mean(life$life_exp)
sigma_hat <- sd(life$life_exp)
tab(data.frame(
calculation = c("Normal approximation: P(X < 70)",
"Normal approximation: P(65 < X < 80)",
"Normal approximation: 90th percentile"),
result = c(pnorm(70, mu_hat, sigma_hat),
pnorm(80, mu_hat, sigma_hat) - pnorm(65, mu_hat, sigma_hat),
qnorm(.90, mu_hat, sigma_hat))))| calculation | result |
|---|---|
| Normal approximation: P(X < 70) | 0.297 |
| Normal approximation: P(65 < X < 80) | 0.700 |
| Normal approximation: 90th percentile | 82.919 |
hist(life$life_exp, probability = TRUE, breaks = "FD",
main = "Observed distribution and normal approximation",
xlab = "Life expectancy (years)", ylab = "Density")
curve(dnorm(x, mu_hat, sigma_hat), add = TRUE, lwd = 2)
8.3 The t and F distributions
The t distribution is symmetric with heavier tails than the standard normal, especially at low degrees of freedom. It appears when estimating a mean with an estimated SD under the usual assumptions. The F distribution is non-negative and appears in variance-ratio statistics such as ANOVA.
oldpar <- par(mfrow = c(1, 2))
curve(dnorm(x), from = -4, to = 4, ylab = "Density", xlab = "Statistic",
main = "Normal and t reference densities")
curve(dt(x, df = 5), add = TRUE, lty = 2)
curve(dt(x, df = 30), add = TRUE, lty = 3)
legend("topright", c("Normal", "t: df = 5", "t: df = 30"), lty = 1:3, bty = "n")
curve(df(x, df1 = 3, df2 = 100), from = 0, to = 5,
main = "F reference density", xlab = "F statistic", ylab = "Density")
par(oldpar)
The reference distributions are mathematical objects, not fabricated datasets. The numerical degrees of freedom here are teaching settings.
9 Sampling, standard errors and the central limit theorem
9.1 SD is not SE
The SD describes variation among observations. The standard error describes uncertainty in an estimator under a specified sampling or model framework. For an independent, identically distributed sample, the estimated SE of the mean is
\[SE(\bar{x})=\frac{s}{\sqrt{n}}.\]
Increasing \(n\) generally reduces this SE. It does not necessarily reduce the underlying SD, remove bias or make correlated observations independent.
x <- life$life_exp
tab(data.frame(n = length(x), mean = mean(x), sd = sd(x),
se_independent_model = sd(x) / sqrt(length(x))))| n | mean | sd | se_independent_model |
|---|---|---|---|
| 217 | 73.791 | 7.123 | 0.484 |
9.2 Empirical resampling illustration
The central limit theorem concerns a distribution of sample means under repeated sampling, not the shape of individual observations. The next exercise draws with replacement from the observed empirical distribution. It illustrates sampling behaviour conditional on that distribution, not a new World Bank dataset or a true survey design.
set.seed(params$seed + 1L)
sample_sizes <- c(5, 20, 60)
resampled_means <- lapply(sample_sizes, function(m) {
replicate(B, mean(sample(life$life_exp, size = m, replace = TRUE)))
})
oldpar <- par(mfrow = c(1, 3))
for (j in seq_along(sample_sizes)) {
hist(resampled_means[[j]], breaks = "FD",
main = paste("Resample n =", sample_sizes[j]),
xlab = "Mean life expectancy (years)", ylab = "Frequency")
}
par(oldpar)
Compare the centres and spreads. The means should concentrate as resample size increases. Do not conclude that a small sample is always adequate, or that the CLT repairs a biased selection mechanism. A simple finite-population random sample without replacement would require a different SE calculation, including a finite-population correction when appropriate.
10 Confidence intervals and bootstrap uncertainty
10.1 A t interval for a mean
A two-sided 95% t interval is
\[\bar{x}\pm t_{0.975,n-1}\frac{s}{\sqrt n}.\]
Under the model and repeated application of the procedure, 95% of such intervals cover the fixed target mean. A particular frequentist interval does not assign a 95% probability to a fixed parameter. It also does not contain 95% of individual country values.
x <- life$life_exp
se_mean <- sd(x) / sqrt(length(x))
critical_t <- qt(.975, df = length(x) - 1)
ci_mean <- mean(x) + c(-1, 1) * critical_t * se_mean
tab(data.frame(mean = mean(x), se = se_mean,
lower95 = ci_mean[1], upper95 = ci_mean[2], n = length(x)))| mean | se | lower95 | upper95 | n |
|---|---|---|---|---|
| 73.791 | 0.484 | 72.838 | 74.744 | 217 |
A narrower interval indicates greater precision under the model, not necessarily less bias. Changing from 95% to 99% makes the interval wider. The country-level independence assumption remains substantive.
10.2 Non-parametric bootstrap of the mean and median
set.seed(params$seed + 2L)
boot_stats <- replicate(B, {
z <- sample(life$life_exp, replace = TRUE)
c(mean = mean(z), median = median(z))
})
boot_intervals <- t(apply(boot_stats, 1, quantile, probs = c(.025, .975)))
tab(data.frame(statistic = rownames(boot_intervals),
estimate = c(mean(life$life_exp), median(life$life_exp)),
lower95 = boot_intervals[, 1], upper95 = boot_intervals[, 2]),
"Percentile bootstrap intervals conditional on the observed economy distribution")| statistic | estimate | lower95 | upper95 |
|---|---|---|---|
| mean | 73.791 | 72.841 | 74.726 |
| median | 74.552 | 72.975 | 76.142 |
The percentile bootstrap is easy to teach but is not universally accurate. It can perform poorly with small samples, heavy tails or boundary parameters. Resampling economies assumes the resampling units are exchangeable and approximately independent. It does not incorporate all uncertainty in World Bank source estimates or fix missing-data bias.
11 Hypothesis testing without shortcuts
11.1 The full reasoning chain
- State the question and target quantity before seeing the p-value.
- Specify a null hypothesis \(H_0\) and an alternative \(H_A\).
- Choose a method appropriate to the design and assumptions.
- Calculate an estimate, SE, test statistic and reference probability.
- Report effect size and uncertainty, not only a threshold decision.
- Interpret the result in the correct population and unit.
- Explain limitations and any sensitivity analyses.
A p-value is the probability, under the null model and its assumptions, of a result at least as incompatible with the null as the observed statistic. It is not the probability that the null is true, the probability that the result happened by chance, or the probability that a replication will be significant.
The significance level \(\alpha\) is a pre-specified decision threshold. A Type I error rejects a true null. A Type II error fails to reject a false null. Power is \(1-\beta\) for a specified alternative. Failure to reject does not prove equality. A tiny p-value does not guarantee an important effect, a valid design or an unbiased estimate.
11.2 Choosing a method
| Question | Design | Starting method | Alternative or sensitivity |
|---|---|---|---|
| Mean versus a fixed reference | One numerical sample | One-sample t-test | Bootstrap; signed-rank only for its own estimand |
| Two means | Independent groups | Welch t-test | Rank-sum or robust methods for a different target |
| Mean within-unit change | Matched pairs | Paired t-test | Signed-rank or sign test |
| Three or more means | Independent groups | Welch ANOVA | Classical ANOVA if appropriate; Kruskal-Wallis for ranks |
| Linear association | Paired numerical variables | Pearson correlation | Spearman correlation |
| Continuous outcome and predictors | Independent observational units | Linear regression | Transformations, robust SE, nonlinear terms |
| Three repeated measurements | Same units at each time | Repeated-measures model | Friedman test; mixed models as extension |
The table is a starting point, not an automatic decision tree. Independence comes from the design. Normality tests cannot verify it. Do not choose the method solely by whether Shapiro-Wilk returns p below .05.
12 One-sample t-test
12.1 Question, hypotheses and assumptions
Teaching question: is the mean country-level life expectancy different from a hypothetical reference of 70 years? The value 70 is an arbitrary teaching benchmark, not an official policy target or a World Bank world mean.
\[H_0:\mu=70,\qquad H_A:\mu\ne70.\]
The statistic is
\[t=\frac{\bar{x}-\mu_0}{s/\sqrt n},\qquad df=n-1.\]
The exact small-sample reference assumes independent normal observations. With sufficiently well-behaved independent data, large-sample approximations may be useful. Severe outliers and dependence still matter.
mu0 <- 70
one_test <- t.test(life$life_exp, mu = mu0, conf.level = .95)
one_test
One Sample t-test
data: life$life_exp
t = 7.8401, df = 216, p-value = 2.036e-13
alternative hypothesis: true mean is not equal to 70
95 percent confidence interval:
72.83790 74.74399
sample estimates:
mean of x
73.79094 manual_t <- (mean(life$life_exp) - mu0) /
(sd(life$life_exp) / sqrt(nrow(life)))
manual_p <- 2 * pt(-abs(manual_t), df = nrow(life) - 1)
tab(data.frame(manual_t = manual_t, r_t = unname(one_test$statistic),
manual_p = manual_p, r_p = one_test$p.value))| manual_t | r_t | manual_p | r_p |
|---|---|---|---|
| 7.84 | 7.84 | 0 | 0 |
stopifnot(isTRUE(all.equal(manual_t, unname(one_test$statistic))))12.2 Read and report the result
The estimate printed by t.test() is the mean, not the mean minus 70. The interval is for the mean. To obtain an interval for the difference from 70, subtract 70 from both limits.
one_difference <- mean(life$life_exp) - mu0
one_difference_ci <- one_test$conf.int - mu0
cat(sprintf(paste0("In %d, the unweighted mean among %d observed economies was ",
"%.2f years (95%% model-based CI %.2f to %.2f). ",
"The difference from the hypothetical 70-year reference was ",
"%.2f years; t(%d) = %.2f, p = %s.\n"),
analysis_year, nrow(life), mean(life$life_exp),
one_test$conf.int[1], one_test$conf.int[2], one_difference,
nrow(life) - 1, unname(one_test$statistic), pfmt(one_test$p.value)))In 2023, the unweighted mean among 217 observed economies was 73.79 years (95% model-based CI 72.84 to 74.74). The difference from the hypothetical 70-year reference was 3.79 years; t(216) = 7.84, p = <0.001.A two-sided test is used because either direction would be relevant. A one-sided test requires a justified direction chosen before inspecting results; it is not a way to halve an inconvenient p-value. See the R t-test documentation.
13 Independent-samples t-tests
13.1 Compare low-income and high-income economies
We compare life expectancy between two current metadata classifications. The contrast is Low income minus High income. The groups contain different economies. They are not matched pairs and were not randomised to income categories.
\[H_0:\mu_L-\mu_H=0,\qquad H_A:\mu_L-\mu_H\ne0.\]
Welch’s test estimates each group’s variance separately:
\[t=\frac{\bar{x}_L-\bar{x}_H}{\sqrt{s_L^2/n_L+s_H^2/n_H}}.\]
It uses a Satterthwaite degrees-of-freedom approximation. It does not require equal group variances. R’s t.test() uses Welch’s method by default for independent samples.
low <- independent$life_exp[independent$income_group == "Low income"]
high <- independent$life_exp[independent$income_group == "High income"]
welch_test <- t.test(low, high, var.equal = FALSE)
welch_test
Welch Two Sample t-test
data: low and high
t = -15.232, df = 36.292, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-17.68188 -13.52759
sample estimates:
mean of x mean of y
64.20236 79.80710 welch_se <- sqrt(var(low) / length(low) + var(high) / length(high))
welch_df <- (var(low) / length(low) + var(high) / length(high))^2 /
((var(low) / length(low))^2 / (length(low) - 1) +
(var(high) / length(high))^2 / (length(high) - 1))
tab(data.frame(contrast = "Low minus High", difference = mean(low) - mean(high),
se = welch_se, df = welch_df,
lower95 = welch_test$conf.int[1], upper95 = welch_test$conf.int[2],
p = welch_test$p.value))| contrast | difference | se | df | lower95 | upper95 | p |
|---|---|---|---|---|---|---|
| Low minus High | -15.605 | 1.024 | 36.292 | -17.682 | -13.528 | 0 |
13.2 The pooled-variance Student t-test
The pooled test assumes a common variance. Its pooled SD is
\[s_p=\sqrt{\frac{(n_L-1)s_L^2+(n_H-1)s_H^2}{n_L+n_H-2}}.\]
pooled_test <- t.test(low, high, var.equal = TRUE)
pooled_test
Two Sample t-test
data: low and high
t = -16.069, df = 109, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-17.52946 -13.68002
sample estimates:
mean of x mean of y
64.20236 79.80710 pooled_sd <- sqrt(((length(low) - 1) * var(low) +
(length(high) - 1) * var(high)) /
(length(low) + length(high) - 2))Use this comparison to understand the assumption, not to shop for the smaller p-value. A non-significant variance test does not establish equal variances. Welch’s method is often a sensible starting choice for unequal group sizes and variances, but it still requires a defensible independent-observation model.
13.3 Standardised effect sizes and a bootstrap interval
Cohen’s \(d=(\bar{x}_L-\bar{x}_H)/s_p\). Hedges’ \(g\) applies a small-sample correction. We use the common approximation \(J=1-3/(4df-1)\), with \(df=n_L+n_H-2\). The raw difference in years is the primary interpretable effect. The standardised effect depends on the chosen SD and is not a universal measure of health importance.
cohen_d <- (mean(low) - mean(high)) / pooled_sd
df_pooled <- length(low) + length(high) - 2
hedges_g <- (1 - 3 / (4 * df_pooled - 1)) * cohen_d
set.seed(params$seed + 3L)
boot_difference <- replicate(B,
mean(sample(low, replace = TRUE)) - mean(sample(high, replace = TRUE)))
boot_difference_ci <- quantile(boot_difference, c(.025, .975))
tab(data.frame(difference_years = mean(low) - mean(high),
cohen_d = cohen_d, hedges_g_approx = hedges_g,
bootstrap_lower95 = boot_difference_ci[1],
bootstrap_upper95 = boot_difference_ci[2]))| difference_years | cohen_d | hedges_g_approx | bootstrap_lower95 | bootstrap_upper95 |
|---|---|---|---|---|
| -15.605 | -3.651 | -3.626 | -17.46 | -13.592 |
Resampling separately within groups preserves group sizes. The bootstrap interval targets the raw mean difference, not \(g\). When variances differ substantially, explain the limitations of a pooled-SD effect size rather than presenting it as uniquely correct.
13.4 Report the comparison
cat(sprintf(paste0("Using current metadata income classifications, mean life ",
"expectancy was %.2f years in low-income economies (n = %d) ",
"and %.2f years in high-income economies (n = %d). ",
"The Low-minus-High difference was %.2f years ",
"(95%% model-based CI %.2f to %.2f), Welch t(%.1f) = %.2f, ",
"p = %s. This is an ecological association, not an income intervention effect.\n"),
mean(low), length(low), mean(high), length(high), mean(low) - mean(high),
welch_test$conf.int[1], welch_test$conf.int[2],
unname(welch_test$parameter), unname(welch_test$statistic),
pfmt(welch_test$p.value)))Using current metadata income classifications, mean life expectancy was 64.20 years in low-income economies (n = 25) and 79.81 years in high-income economies (n = 86). The Low-minus-High difference was -15.60 years (95% model-based CI -17.68 to -13.53), Welch t(36.3) = -15.23, p = <0.001. This is an ecological association, not an income intervention effect.14 Paired t-test: within-economy change
14.1 Why pairing changes the analysis
The same economy is observed in two years. The analysis uses differences \(d_i=\text{after}_i-\text{before}_i\):
\[t=\frac{\bar d-0}{s_d/\sqrt n},\qquad df=n-1.\]
The normality assumption concerns the differences, not a requirement that both marginal year distributions be normal. Independence is required between pairs, not within each pair. A paired t-test can be implemented as a one-sample t-test on differences.
paired_test <- t.test(paired$after, paired$before, paired = TRUE)
difference_test <- t.test(paired$change, mu = 0)
paired_test
Paired t-test
data: paired$after and paired$before
t = 15.547, df = 216, p-value < 2.2e-16
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
2.633461 3.398144
sample estimates:
mean difference
3.015802 stopifnot(isTRUE(all.equal(unname(paired_test$statistic),
unname(difference_test$statistic))))
tab(data.frame(n_pairs = nrow(paired), mean_change = mean(paired$change),
sd_change = sd(paired$change),
lower95 = paired_test$conf.int[1], upper95 = paired_test$conf.int[2],
p = paired_test$p.value))| n_pairs | mean_change | sd_change | lower95 | upper95 | p |
|---|---|---|---|---|---|
| 217 | 3.016 | 2.858 | 2.633 | 3.398 | 0 |
hist(paired$change, breaks = "FD", xlab = "Later minus earlier life expectancy (years)",
main = paste("Matched changes:", baseline_year, "to", analysis_year),
ylab = "Number of economies")
abline(v = 0, lty = 2)
14.2 Paired effect size and interpretation
paired_dz <- mean(paired$change) / sd(paired$change)
cat(sprintf(paste0("Among %d matched economies, mean life expectancy changed by ",
"%.2f years from %d to %d (95%% model-based CI %.2f to %.2f), ",
"paired t(%d) = %.2f, p = %s. Standardised mean change dz = %.2f.\n"),
nrow(paired), mean(paired$change), baseline_year, analysis_year,
paired_test$conf.int[1], paired_test$conf.int[2], nrow(paired) - 1,
unname(paired_test$statistic), pfmt(paired_test$p.value), paired_dz))Among 217 matched economies, mean life expectancy changed by 3.02 years from 2010 to 2023 (95% model-based CI 2.63 to 3.40), paired t(216) = 15.55, p = <0.001. Standardised mean change dz = 1.06.\(d_z\) uses the SD of changes. It is not interchangeable with an independent-groups Cohen’s \(d\). A temporal change is not automatically the effect of a programme. Many processes changed between the two years, and the same country’s estimates can be revised retrospectively.
15 One-way ANOVA: comparing more than two means
15.1 The question and the F statistic
For four current income groups, the null is
\[H_0:\mu_1=\mu_2=\mu_3=\mu_4.\]
The alternative is that at least one mean differs. It is not that all means differ. Classical one-way ANOVA partitions variation:
\[SS_T=SS_B+SS_W,\]
\[F=\frac{MS_B}{MS_W}=\frac{SS_B/(k-1)}{SS_W/(N-k)}.\]
\(SS_B\) describes between-group variation around the grand mean. \(SS_W\) describes within-group variation around group means. The classical F reference assumes independent observations, approximately normal errors within groups and a common error variance. With unequal group sizes, unequal variances deserve particular attention.
anova_fit <- aov(life_exp ~ income_group, data = anova_data)
summary(anova_fit) Df Sum Sq Mean Sq F value Pr(>F)
income_group 3 7037 2345.7 127.4 <2e-16 ***
Residuals 213 3922 18.4
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 115.2 Calculate the sums of squares directly
split_y <- split(anova_data$life_exp, anova_data$income_group)
grand_mean <- mean(anova_data$life_exp)
N <- nrow(anova_data)
k <- length(split_y)
ss_between <- sum(vapply(split_y, function(z) length(z) * (mean(z) - grand_mean)^2,
numeric(1)))
ss_within <- sum(vapply(split_y, function(z) sum((z - mean(z))^2), numeric(1)))
ss_total <- sum((anova_data$life_exp - grand_mean)^2)
f_manual <- (ss_between / (k - 1)) / (ss_within / (N - k))
tab(data.frame(source = c("Between", "Within", "Total"),
SS = c(ss_between, ss_within, ss_total),
df = c(k - 1, N - k, N - 1)))| source | SS | df |
|---|---|---|
| Between | 7037.01 | 3 |
| Within | 3921.81 | 213 |
| Total | 10958.82 | 216 |
stopifnot(isTRUE(all.equal(ss_total, ss_between + ss_within)))
tab(data.frame(F = f_manual, p = pf(f_manual, k - 1, N - k, lower.tail = FALSE)))| F | p |
|---|---|
| 127.397 | 0 |
15.3 ANOVA effect sizes
For one-way ANOVA, \(\eta^2=SS_B/SS_T\). A less upward-biased estimate is
\[\omega^2=\frac{SS_B-(k-1)MS_W}{SS_T+MS_W}.\]
ms_within <- ss_within / (N - k)
eta_squared <- ss_between / ss_total
omega_squared <- (ss_between - (k - 1) * ms_within) / (ss_total + ms_within)
tab(data.frame(eta_squared = eta_squared, omega_squared = omega_squared))| eta_squared | omega_squared |
|---|---|
| 0.642 | 0.636 |
A negative estimated omega squared can occur through sampling variation and is often interpreted as no explained variance; do not silently alter the calculation. These quantities describe variation associated with the grouping in the observed analysis, not a causal fraction of life expectancy determined by income.
16 ANOVA assumptions, Welch ANOVA and post-hoc comparisons
16.1 Examine residuals and spread
oldpar <- par(mfrow = c(1, 2))
plot(fitted(anova_fit), residuals(anova_fit), pch = 16,
xlab = "Fitted group mean", ylab = "Residual", main = "Residuals vs fitted")
abline(h = 0, lty = 2)
qqnorm(residuals(anova_fit), main = "ANOVA residual Q-Q plot")
qqline(residuals(anova_fit))
par(oldpar)
A Q-Q plot compares ordered residuals with theoretical normal quantiles. A roughly straight pattern is compatible with normality; departures can indicate skewness, heavy tails or outliers. It does not verify independence or correct model specification.
A median-centred Levene-style check, often called the Brown-Forsythe variance test, analyses absolute deviations from group medians. It is a diagnostic, not a gatekeeper that automatically selects a final test.
anova_data$abs_deviation <- with(anova_data,
abs(life_exp - ave(life_exp, income_group, FUN = median)))
summary(aov(abs_deviation ~ income_group, data = anova_data)) Df Sum Sq Mean Sq F value Pr(>F)
income_group 3 25.3 8.430 1.073 0.362
Residuals 213 1673.8 7.858 16.2 Welch ANOVA
welch_anova <- oneway.test(life_exp ~ income_group, data = anova_data,
var.equal = FALSE)
welch_anova
One-way analysis of means (not assuming equal variances)
data: life_exp and income_group
F = 115.17, num df = 3.000, denom df = 81.929, p-value < 2.2e-16Welch ANOVA compares means without requiring equal group variances. It is not the same as a rank test. See R’s one-way mean test documentation.
16.3 Tukey comparisons under the classical model
tukey <- TukeyHSD(anova_fit, conf.level = .95)
tukey Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = life_exp ~ income_group, data = anova_data)
$income_group
diff lwr upr
Lower middle income-Low income 3.707189 0.9567336 6.457643
Upper middle income-Low income 9.567413 6.9158571 12.218970
High income-Low income 15.604736 13.0800963 18.129377
Upper middle income-Lower middle income 5.860225 3.6878478 8.032602
High income-Lower middle income 11.897548 9.8820382 13.913058
High income-Upper middle income 6.037323 4.1590185 7.915627
p adj
Lower middle income-Low income 0.0032663
Upper middle income-Low income 0.0000000
High income-Low income 0.0000000
Upper middle income-Lower middle income 0.0000000
High income-Lower middle income 0.0000000
High income-Upper middle income 0.0000000Tukey’s method addresses the family of pairwise mean comparisons under the classical ANOVA framework. Its intervals are simultaneous within that family. It does not automatically become appropriate after Welch ANOVA, because the variance assumptions differ.
16.4 Pairwise Welch tests with Holm adjustment
pairwise_welch <- pairwise.t.test(anova_data$life_exp, anova_data$income_group,
pool.sd = FALSE, p.adjust.method = "holm")
pairwise_welch
Pairwise comparisons using t tests with non-pooled SD
data: anova_data$life_exp and anova_data$income_group
Low income Lower middle income Upper middle income
Lower middle income 0.0025 - -
Upper middle income 9.8e-11 3.3e-09 -
High income < 2e-16 < 2e-16 6.2e-15
P value adjustment method: holm # A separate table retains raw contrasts and their ordinary individual CIs.
comparisons <- combn(levels(anova_data$income_group), 2, simplify = FALSE)
pairwise_table <- do.call(rbind, lapply(comparisons, function(g) {
a <- anova_data$life_exp[anova_data$income_group == g[1]]
b <- anova_data$life_exp[anova_data$income_group == g[2]]
tt <- t.test(a, b)
data.frame(contrast = paste(g[1], "minus", g[2]),
difference = mean(a) - mean(b), lower95_individual = tt$conf.int[1],
upper95_individual = tt$conf.int[2], p_raw = tt$p.value)
}))
pairwise_table$p_holm <- p.adjust(pairwise_table$p_raw, method = "holm")
tab(pairwise_table, "Welch pairwise contrasts: CIs are individual, not multiplicity-adjusted")| contrast | difference | lower95_individual | upper95_individual | p_raw | p_holm |
|---|---|---|---|---|---|
| Low income minus Lower middle income | -3.707 | -6.045 | -1.370 | 0.002 | 0.002 |
| Low income minus Upper middle income | -9.567 | -11.679 | -7.456 | 0.000 | 0.000 |
| Low income minus High income | -15.605 | -17.682 | -13.528 | 0.000 | 0.000 |
| Lower middle income minus Upper middle income | -5.860 | -7.584 | -4.136 | 0.000 | 0.000 |
| Lower middle income minus High income | -11.898 | -13.577 | -10.218 | 0.000 | 0.000 |
| Upper middle income minus High income | -6.037 | -7.357 | -4.718 | 0.000 | 0.000 |
Do not describe the individual 95% intervals in this table as simultaneous intervals. Holm-adjusted p-values and unadjusted intervals answer different inferential questions. Games-Howell is another unequal-variance post-hoc option, but is not required for this dependency-light workbook.
16.5 Why multiplicity matters
With many tests, the chance of at least one false rejection can increase. Under independent true-null tests, that chance is \(1-(1-\alpha)^m\). Holm controls family-wise error for a defined family. Benjamini-Hochberg targets the false discovery rate under its relevant conditions. Neither method repairs biased data or a post-hoc choice of hypotheses.
m <- nrow(pairwise_table)
tab(data.frame(number_of_tests = m,
independent_null_family_error_at_05 = 1 - (1 - .05)^m))| number_of_tests | independent_null_family_error_at_05 |
|---|---|
| 6 | 0.265 |
tab(data.frame(contrast = pairwise_table$contrast,
raw = pairwise_table$p_raw,
holm = p.adjust(pairwise_table$p_raw, "holm"),
BH = p.adjust(pairwise_table$p_raw, "BH")))| contrast | raw | holm | BH |
|---|---|---|---|
| Low income minus Lower middle income | 0.002 | 0.002 | 0.002 |
| Low income minus Upper middle income | 0.000 | 0.000 | 0.000 |
| Low income minus High income | 0.000 | 0.000 | 0.000 |
| Lower middle income minus Upper middle income | 0.000 | 0.000 | 0.000 |
| Lower middle income minus High income | 0.000 | 0.000 | 0.000 |
| Upper middle income minus High income | 0.000 | 0.000 | 0.000 |
The R p-value adjustment documentation explains the available methods. Define the family of comparisons before choosing an adjustment.
17 Non-parametric methods
17.1 Rank-sum test for two independent groups
The Wilcoxon rank-sum or Mann-Whitney test uses ranks. It is not a general-purpose test of equality of medians. A location-shift or median interpretation requires additional assumptions about distribution shapes. A rank-based result can differ from a mean-based result because the estimand and assumptions differ.
rank_sum <- wilcox.test(low, high, exact = FALSE, conf.int = TRUE)
rank_sum
Wilcoxon rank sum test with continuity correction
data: low and high
W = 27, p-value = 1.418e-13
alternative hypothesis: true location shift is not equal to 0
95 percent confidence interval:
-17.91896 -14.07842
sample estimates:
difference in location
-15.8435 # Probability-of-superiority estimate with half weight for ties.
comparison_matrix <- outer(low, high, "-")
prob_low_exceeds_high <- mean(comparison_matrix > 0) +
.5 * mean(comparison_matrix == 0)
rank_biserial <- 2 * prob_low_exceeds_high - 1
tab(data.frame(P_low_exceeds_high_half_ties = prob_low_exceeds_high,
rank_biserial_low_minus_high = rank_biserial))| P_low_exceeds_high_half_ties | rank_biserial_low_minus_high |
|---|---|
| 0.013 | -0.975 |
The probability-of-superiority estimate compares a randomly selected observed low-income economy with a randomly selected observed high-income economy. It is not an individual patient’s probability of living longer. The confidence interval printed by wilcox.test() is for its location-shift parameter, not automatically the difference of sample medians. See the Wilcoxon documentation.
17.2 Signed-rank and sign tests for paired changes
The signed-rank test assumes an appropriate symmetry condition for a location interpretation. It is not simply a paired t-test with no assumptions. A sign test uses only whether each non-zero change is positive or negative and discards magnitude.
signed_rank <- wilcox.test(paired$change, mu = 0, exact = FALSE,
conf.int = TRUE)
signed_rank
Wilcoxon signed rank test with continuity correction
data: paired$change
V = 23203, p-value < 2.2e-16
alternative hypothesis: true location is not equal to 0
95 percent confidence interval:
2.366970 2.932022
sample estimates:
(pseudo)median
2.642955 nonzero_change <- paired$change[paired$change != 0]
sign_test <- binom.test(sum(nonzero_change > 0), length(nonzero_change), p = .5)
sign_test
Exact binomial test
data: sum(nonzero_change > 0) and length(nonzero_change)
number of successes = 207, number of trials = 217, p-value < 2.2e-16
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
0.9168867 0.9776833
sample estimates:
probability of success
0.9539171 The binomial calculation here supports a sign test for continuous paired changes; it is not a separate categorical-outcome curriculum. Zero differences are excluded explicitly. With a continuous change distribution and suitable assumptions, the sign test addresses whether the median change is zero.
17.3 Kruskal-Wallis for more than two groups
kw <- kruskal.test(life_exp ~ income_group, data = anova_data)
kw
Kruskal-Wallis rank sum test
data: life_exp by income_group
Kruskal-Wallis chi-squared = 141.68, df = 3, p-value < 2.2e-16pairwise.wilcox.test(anova_data$life_exp, anova_data$income_group,
p.adjust.method = "holm", exact = FALSE)
Pairwise comparisons using Wilcoxon rank sum test with continuity correction
data: anova_data$life_exp and anova_data$income_group
Low income Lower middle income Upper middle income
Lower middle income 0.0032 - -
Upper middle income 6.5e-10 1.4e-08 -
High income 5.7e-13 < 2e-16 2.1e-13
P value adjustment method: holm Kruskal-Wallis is an omnibus rank-based comparison. Under common-shape location-shift assumptions it can be interpreted as a location comparison. If spreads and shapes differ, a rejection is not automatically evidence that medians differ. Pairwise Wilcoxon comparisons with Holm adjustment are shown as a follow-up family; they are not Dunn’s test. See R’s Kruskal-Wallis documentation.
18 Correlation: association between continuous variables
18.1 Pearson correlation
Pearson’s correlation measures linear association:
\[r=\frac{\sum_i(x_i-\bar{x})(y_i-\bar{y})} {\sqrt{\sum_i(x_i-\bar{x})^2\sum_i(y_i-\bar{y})^2}}.\]
It lies between -1 and 1 and is unitless. A value near zero does not rule out a strong nonlinear relationship. Outliers can strongly affect it. The usual exact t test for zero correlation relies on independent paired observations and a bivariate normal model.
pearson <- cor.test(reg$log2_health, reg$life_exp, method = "pearson")
pearson
Pearson's product-moment correlation
data: reg$log2_health and reg$life_exp
t = 20.607, df = 185, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.7851204 0.8734768
sample estimates:
cor
0.8345893 18.2 Spearman correlation
Spearman’s correlation measures association between ranks and is useful for monotonic relationships. It does not require the same bivariate normal model as the usual Pearson test. It still requires a defensible independent-unit framework and can be affected by ties.
spearman <- cor.test(reg$log2_health, reg$life_exp,
method = "spearman", exact = FALSE)
spearman
Spearman's rank correlation rho
data: reg$log2_health and reg$life_exp
S = 167948, p-value < 2.2e-16
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
0.8458965 cor_matrix <- cor(reg[reg_vars], method = "pearson")
round(cor_matrix, 3) life_exp log2_health log2_gdp urban10
life_exp 1.000 0.835 0.858 0.563
log2_health 0.835 1.000 0.960 0.639
log2_gdp 0.858 0.960 1.000 0.632
urban10 0.563 0.639 0.632 1.000We use the same complete-case dataset for this matrix so every entry has the same denominator. Pairwise deletion would give different denominators and can produce a correlation matrix with problematic mathematical properties. See R’s correlation-test documentation.
Interpretation prompt: explain why correlation between GDP and expenditure can complicate a multiple regression that includes both. Do not use a fixed label such as ‘strong’ without considering the context and graph.
19 Simple linear regression
19.1 Model, coefficients and residuals
The simple model is
\[Y_i=\beta_0+\beta_1X_i+\epsilon_i,\]
where \(Y_i\) is life expectancy and \(X_i\) is log base 2 health expenditure. Ordinary least squares chooses coefficients that minimise the sum of squared residuals. A residual is \(e_i=y_i-\hat y_i\). It is an observed model discrepancy, not the unobservable true error itself.
simple_fit <- lm(life_exp ~ log2_health, data = reg)
summary(simple_fit)
Call:
lm(formula = life_exp ~ log2_health, data = reg)
Residuals:
Min 1Q Median 3Q Max
-16.9791 -1.7987 0.1428 2.7877 9.2032
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 51.0499 1.1141 45.82 <2e-16 ***
log2_health 2.5118 0.1219 20.61 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.959 on 185 degrees of freedom
Multiple R-squared: 0.6965, Adjusted R-squared: 0.6949
F-statistic: 424.6 on 1 and 185 DF, p-value: < 2.2e-16confint(simple_fit) 2.5 % 97.5 %
(Intercept) 48.851820 53.24795
log2_health 2.271315 2.75227The slope is the expected difference in country-level life expectancy associated with a doubling of health expenditure in this model. The intercept corresponds to expenditure of one current US dollar per capita because log2(1) = 0; it may be far outside the observed range and should not be given a substantive interpretation.
19.2 Derive the least-squares slope
x <- reg$log2_health
y <- reg$life_exp
b1 <- sum((x - mean(x)) * (y - mean(y))) / sum((x - mean(x))^2)
b0 <- mean(y) - b1 * mean(x)
tab(data.frame(term = c("Intercept", "Log2 expenditure"),
manual = c(b0, b1), fitted = unname(coef(simple_fit))))| term | manual | fitted |
|---|---|---|
| Intercept | 51.050 | 51.050 |
| Log2 expenditure | 2.512 | 2.512 |
stopifnot(isTRUE(all.equal(b1, unname(coef(simple_fit)[2]))))19.3 A readable coefficient table
coef_table <- function(model) {
sm <- coef(summary(model))
ci <- confint(model)
data.frame(term = rownames(sm), estimate = sm[, 1], se = sm[, 2],
lower95 = ci[, 1], upper95 = ci[, 2],
t = sm[, 3], p = sm[, 4], row.names = NULL)
}
simple_table <- coef_table(simple_fit)
tab(simple_table, "Simple linear regression: life expectancy in years")| term | estimate | se | lower95 | upper95 | t | p |
|---|---|---|---|---|---|---|
| (Intercept) | 51.050 | 1.114 | 48.852 | 53.248 | 45.820 | 0 |
| log2_health | 2.512 | 0.122 | 2.271 | 2.752 | 20.607 | 0 |
cat(sprintf(paste0("A doubling of health expenditure was associated with a ",
"%.2f-year difference in expected country-level life expectancy ",
"(95%% model-based CI %.2f to %.2f). This is not a causal estimate.\n"),
coef(simple_fit)[2], confint(simple_fit)[2, 1], confint(simple_fit)[2, 2]))A doubling of health expenditure was associated with a 2.51-year difference in expected country-level life expectancy (95% model-based CI 2.27 to 2.75). This is not a causal estimate.19.4 R-squared and residual standard error
\[R^2=1-\frac{\sum_i(y_i-\hat y_i)^2}{\sum_i(y_i-\bar y)^2}.\]
\(R^2\) describes the fraction of observed outcome variation accounted for by the fitted model relative to an intercept-only model. It is not the percentage of life expectancy caused by expenditure. A high \(R^2\) does not prove good causal identification or good prediction in new settings.
ss_residual <- sum(residuals(simple_fit)^2)
ss_outcome <- sum((reg$life_exp - mean(reg$life_exp))^2)
r2_manual <- 1 - ss_residual / ss_outcome
tab(data.frame(n = nobs(simple_fit), r_squared = r2_manual,
adjusted_r_squared = summary(simple_fit)$adj.r.squared,
residual_standard_error_years = sigma(simple_fit)))| n | r_squared | adjusted_r_squared | residual_standard_error_years |
|---|---|---|---|
| 187 | 0.697 | 0.695 | 3.959 |
# In simple regression with an intercept and the same rows, R-squared = r^2.
stopifnot(isTRUE(all.equal(r2_manual, cor(x, y)^2)))20 Multiple linear regression and confounding
20.1 Add GDP and urbanisation
Our working model is
\[E(Y_i\mid X_i)=\beta_0+\beta_1\log_2(H_i) +\beta_2\log_2(G_i)+\beta_3(U_i/10).\]
\(H_i\) is health expenditure per capita, \(G_i\) is GDP per capita, and \(U_i\) is urban population percentage. Dividing urban percentage by 10 makes its coefficient correspond to a 10-percentage-point difference. That is not a 10% relative increase.
multiple_fit <- lm(life_exp ~ log2_health + log2_gdp + urban10, data = reg)
summary(multiple_fit)
Call:
lm(formula = life_exp ~ log2_health + log2_gdp + urban10, data = reg)
Residuals:
Min 1Q Median 3Q Max
-14.473 -2.003 0.370 2.578 7.711
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 37.69450 2.70969 13.911 < 2e-16 ***
log2_health 0.39212 0.40942 0.958 0.339
log2_gdp 2.44799 0.46231 5.295 0.000000338 ***
urban10 0.08959 0.15714 0.570 0.569
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.697 on 183 degrees of freedom
Multiple R-squared: 0.7382, Adjusted R-squared: 0.7339
F-statistic: 172 on 3 and 183 DF, p-value: < 2.2e-16multiple_table <- coef_table(multiple_fit)
tab(multiple_table, "Multiple linear regression: adjusted country-level associations")| term | estimate | se | lower95 | upper95 | t | p |
|---|---|---|---|---|---|---|
| (Intercept) | 37.694 | 2.710 | 32.348 | 43.041 | 13.911 | 0.000 |
| log2_health | 0.392 | 0.409 | -0.416 | 1.200 | 0.958 | 0.339 |
| log2_gdp | 2.448 | 0.462 | 1.536 | 3.360 | 5.295 | 0.000 |
| urban10 | 0.090 | 0.157 | -0.220 | 0.400 | 0.570 | 0.569 |
The expenditure coefficient now compares economies with the same modelled GDP and urbanisation values. ‘Holding constant’ is a mathematical conditioning statement, not evidence that an experiment held those factors fixed. Adjustment can reduce some confounding only under substantive assumptions. It can also introduce bias if inappropriate variables, such as colliders or mediators for the intended causal question, are controlled.
20.2 Compare crude and adjusted coefficients on the same rows
crude_adjusted <- data.frame(
model = c("Expenditure only", "Expenditure + GDP + urbanisation"),
n = c(nobs(simple_fit), nobs(multiple_fit)),
expenditure_slope = c(coef(simple_fit)["log2_health"],
coef(multiple_fit)["log2_health"]))
tab(crude_adjusted)| model | n | expenditure_slope |
|---|---|---|
| Expenditure only | 187 | 2.512 |
| Expenditure + GDP + urbanisation | 187 | 0.392 |
Because both models use reg, a change in coefficient is not caused merely by changing the complete-case sample. It still does not prove which covariate is a confounder or establish causality. GDP and health expenditure can be highly related, which can make their separate conditional coefficients imprecise.
20.3 Nested model comparison
anova(simple_fit, multiple_fit)The partial F test asks whether the added GDP and urbanisation terms jointly improve fit under the classical linear-model assumptions. Models must be nested and fitted to the same observations. This is not a test of causal correctness. See R’s linear-model documentation.
21 Regression assumptions and diagnostics
21.1 What is actually assumed?
For the standard model, ask whether the conditional mean is specified appropriately, whether errors are independent, whether the error variance is constant for conventional SEs, and whether the design matrix has full rank. Normal errors support exact small-sample t and F inference. Predictors do not need to be normally distributed. The outcome does not need to have a normal marginal distribution.
A causal interpretation additionally needs strong identification assumptions, including an appropriate design and control of relevant confounding. Diagnostic plots cannot establish those assumptions.
oldpar <- par(mfrow = c(2, 2))
plot(multiple_fit, which = c(1, 2, 3, 5))
par(oldpar)
- Residuals versus fitted: curvature suggests a misspecified conditional mean; a fan suggests changing variance.
- Normal Q-Q: substantial tail departures can affect small-sample inference and intervals.
- Scale-location: helps inspect changing residual spread.
- Residuals versus leverage: combines unusual predictor combinations with residual behaviour.
21.2 Formal checks as supporting information
# Shapiro-Wilk is defined for 3 to 5000 non-missing observations.
e <- residuals(multiple_fit)
if (length(e) >= 3 && length(e) <= 5000) shapiro.test(e)
Shapiro-Wilk normality test
data: e
W = 0.9709, p-value = 0.0006128# Auxiliary n*R-squared check for variance related to these predictors.
# This is a teaching LM diagnostic with an asymptotic chi-square reference.
aux_data <- reg
aux_data$squared_residual <- residuals(multiple_fit)^2
aux_fit <- lm(squared_residual ~ log2_health + log2_gdp + urban10, data = aux_data)
lm_stat <- nrow(aux_data) * summary(aux_fit)$r.squared
tab(data.frame(variance_LM_statistic = lm_stat, df = 3,
asymptotic_p = pchisq(lm_stat, df = 3, lower.tail = FALSE)))| variance_LM_statistic | df | asymptotic_p |
|---|---|---|
| 9.775 | 3 | 0.021 |
A normality test can detect tiny deviations in a large sample and miss serious deviations in a small one. The auxiliary variance check is asymptotic and tests only variance patterns represented by its predictors. Neither result should replace graphical inspection and design reasoning.
21.3 Multicollinearity
A variance inflation factor for a numerical predictor is \(VIF_j=1/(1-R_j^2)\), where \(R_j^2\) comes from regressing that predictor on the others. This helper is restricted to the numerical main-effect predictors used here; it is not a generalised VIF implementation for multi-degree-of-freedom factors.
predictors <- c("log2_health", "log2_gdp", "urban10")
vifs <- vapply(predictors, function(v) {
others <- setdiff(predictors, v)
auxiliary <- lm(reformulate(others, response = v), data = reg)
1 / (1 - summary(auxiliary)$r.squared)
}, numeric(1))
tab(data.frame(predictor = names(vifs), VIF = unname(vifs)))| predictor | VIF |
|---|---|
| log2_health | 12.936 |
| log2_gdp | 12.751 |
| urban10 | 1.702 |
Large VIFs indicate difficulty separating conditional associations. They do not necessarily make prediction poor, and a fixed cutoff is not a universal law. Do not remove a substantively necessary confounder solely to reduce a VIF.
21.4 Outliers, leverage and influence
influence_table <- data.frame(country = reg$country, iso3c = reg$iso3c,
studentised_residual = rstudent(multiple_fit),
leverage = hatvalues(multiple_fit),
cooks_distance = cooks.distance(multiple_fit))
influence_table <- influence_table[order(influence_table$cooks_distance,
decreasing = TRUE), ]
tab(head(influence_table, 10), "Most influential observations by Cook's distance")| country | iso3c | studentised_residual | leverage | cooks_distance |
|---|---|---|---|---|
| Naoero | NRU | -4.181 | 0.044 | 0.185 |
| Nigeria | NGA | -3.750 | 0.022 | 0.073 |
| Tuvalu | TUV | -1.576 | 0.082 | 0.055 |
| Guyana | GUY | -1.776 | 0.059 | 0.049 |
| Lesotho | LSO | -2.099 | 0.041 | 0.046 |
| Sri Lanka | LKA | 2.138 | 0.030 | 0.035 |
| Lebanon | LBN | 2.043 | 0.032 | 0.034 |
| Chad | TCD | -2.568 | 0.020 | 0.032 |
| Bangladesh | BGD | 1.866 | 0.034 | 0.030 |
| Afghanistan | AFG | 1.266 | 0.062 | 0.026 |
An outlier has an unusual outcome conditional on the model. A high-leverage point has an unusual predictor combination. An influential point materially affects fitted coefficients or predictions. Cook’s distance combines aspects of residual size and leverage. Rules such as \(4/n\) are screening heuristics, not deletion criteria.
Robust uncertainty and sensitivity analysis
21.5 Extension: HC3 heteroskedasticity-robust standard errors
Heteroskedasticity-consistent covariance estimates allow non-constant error variance under suitable large-sample conditions. They do not fix dependence, confounding, incorrect functional form or measurement error. The next implementation calculates HC3 directly for a full-rank ordinary least-squares model. It uses a t reference with residual degrees of freedom as a conventional approximation, not an exact finite-sample guarantee.
hc3_table <- function(model) {
X <- model.matrix(model)
if (qr(X)$rank != ncol(X)) stop("HC3 helper requires a full-rank model.")
bread <- solve(crossprod(X))
h <- hatvalues(model)
if (any(h >= 1 - 1e-10)) stop("Leverage too close to one for this HC3 helper.")
weights_hc3 <- residuals(model)^2 / (1 - h)^2
meat <- crossprod(X, sweep(X, 1, weights_hc3, "*"))
covariance <- bread %*% meat %*% bread
se <- sqrt(diag(covariance))
b <- coef(model)
df <- df.residual(model)
critical <- qt(.975, df)
data.frame(term = names(b), estimate = unname(b), se_HC3 = unname(se),
lower95_approx = unname(b - critical * se),
upper95_approx = unname(b + critical * se),
p_approx = unname(2 * pt(-abs(b / se), df)))
}
robust_table <- hc3_table(multiple_fit)
tab(robust_table, "HC3 uncertainty: independent-unit large-sample sensitivity analysis")| term | estimate | se_HC3 | lower95_approx | upper95_approx | p_approx |
|---|---|---|---|---|---|
| (Intercept) | 37.694 | 3.227 | 31.327 | 44.062 | 0.000 |
| log2_health | 0.392 | 0.522 | -0.637 | 1.421 | 0.453 |
| log2_gdp | 2.448 | 0.598 | 1.268 | 3.628 | 0.000 |
| urban10 | 0.090 | 0.182 | -0.269 | 0.448 | 0.622 |
The coefficient estimates do not change because we changed the covariance estimator, not the least-squares fit. If economies are dependent, HC3 alone is insufficient. Cluster-robust or spatial methods require a defensible dependence structure and enough independent clusters.
21.6 Extension: leave-one-out coefficient sensitivity
loo_slopes <- vapply(seq_len(nrow(reg)), function(i) {
coef(lm(life_exp ~ log2_health + log2_gdp + urban10, data = reg[-i, ]))[
"log2_health"]
}, numeric(1))
tab(data.frame(full_sample_slope = coef(multiple_fit)["log2_health"],
minimum_leave_one_out = min(loo_slopes),
maximum_leave_one_out = max(loo_slopes)))| full_sample_slope | minimum_leave_one_out | maximum_leave_one_out |
|---|---|---|
| 0.392 | 0.249 | 0.617 |
This asks whether any single economy drives the expenditure coefficient. It does not justify removing the economy that produces an inconvenient answer. Report the full-sample analysis and describe sensitivity transparently.
Categorical predictors, interactions and ANCOVA
21.7 Extension: a categorical predictor is still a linear model
ANOVA and regression share the linear-model framework. A factor with four levels is represented by three indicator columns plus an intercept under treatment coding.
factor_fit <- lm(life_exp ~ income_group, data = anova_data)
tab(coef_table(factor_fit), "Factor regression: reference is the first income level")| term | estimate | se | lower95 | upper95 | t | p |
|---|---|---|---|---|---|---|
| (Intercept) | 64.202 | 0.858 | 62.511 | 65.894 | 74.811 | 0.000 |
| income_groupLower middle income | 3.707 | 1.062 | 1.613 | 5.801 | 3.490 | 0.001 |
| income_groupUpper middle income | 9.567 | 1.024 | 7.549 | 11.586 | 9.343 | 0.000 |
| income_groupHigh income | 15.605 | 0.975 | 13.683 | 17.527 | 16.005 | 0.000 |
head(model.matrix(factor_fit)) (Intercept) income_groupLower middle income income_groupUpper middle income
14 1 0 0
28 1 0 0
42 1 1 0
56 1 0 1
70 1 0 0
84 1 0 0
income_groupHigh income
14 1
28 0
42 0
56 0
70 1
84 1The intercept is the mean for the reference group. Each other coefficient is a difference from that reference. The overall F test from this model corresponds to the classical one-way ANOVA on the same rows.
21.8 Extension: two-factor ANOVA with an interaction
A second factor can be a World Bank region. The example uses two regions with the largest complete observed counts, chosen by a transparent coverage rule rather than a significance result. The second factor is current high-income versus other classified economies. This is an illustrative unbalanced observational factorial design, not a randomised experiment.
factorial_data <- cross[complete.cases(cross[c("life_exp", "region", "high_income")]), ]
region_counts <- sort(table(factorial_data$region), decreasing = TRUE)
chosen_regions <- names(region_counts)[seq_len(min(2L, length(region_counts)))]
factorial_data <- droplevels(factorial_data[factorial_data$region %in% chosen_regions, ])
factorial_data$region <- factor(factorial_data$region)
cell_counts <- table(factorial_data$region, factorial_data$high_income)
cell_counts
Other High
Europe & Central Asia 18 40
Sub-Saharan Africa 47 1factorial_fit <- NULL
if (nrow(cell_counts) == 2 && ncol(cell_counts) == 2 && all(cell_counts >= 3)) {
factorial_additive <- lm(life_exp ~ region + high_income, data = factorial_data)
factorial_fit <- lm(life_exp ~ region * high_income, data = factorial_data)
print(anova(factorial_additive, factorial_fit))
print(coef_table(factorial_fit))
} else {
cat("Insufficient coverage in one or more cells: factorial fit omitted, not fabricated.\n")
}Insufficient coverage in one or more cells: factorial fit omitted, not fabricated.The nested comparison tests the interaction jointly. In an unbalanced design, sequential sums of squares depend on term order; do not interpret a default ANOVA table as if all types of sums of squares were identical. A main effect in a model with an interaction is conditional on the other factor’s reference level.
21.9 Extension: ANCOVA and slope interactions
ANCOVA combines categorical groups with numerical covariates. Here we use current high-income classification and log expenditure. We first centre log expenditure so the group coefficient refers to a typical observed value rather than expenditure of one dollar.
interaction_data <- reg[!is.na(reg$high_income), ]
interaction_data$high_income <- droplevels(interaction_data$high_income)
interaction_centre <- mean(interaction_data$log2_health)
interaction_data$health_c <- interaction_data$log2_health - interaction_centre
ancova_fit <- NULL
interaction_fit <- NULL
if (nlevels(interaction_data$high_income) == 2 &&
all(table(interaction_data$high_income) >= 5)) {
ancova_fit <- lm(life_exp ~ health_c + high_income + log2_gdp + urban10,
data = interaction_data)
interaction_fit <- lm(life_exp ~ health_c * high_income + log2_gdp + urban10,
data = interaction_data)
print(coef_table(interaction_fit))
print(anova(ancova_fit, interaction_fit))
} else {
cat("Insufficient group coverage for the ANCOVA extension.\n")
} term estimate se lower95 upper95
1 (Intercept) 39.70063410 6.1610023 27.5440091 51.8572591
2 health_c 0.41310883 0.4311435 -0.4376050 1.2638226
3 high_incomeHigh -1.09551688 1.3073418 -3.6751075 1.4840737
4 log2_gdp 2.58992262 0.4859157 1.6311346 3.5487106
5 urban10 0.07802942 0.1580571 -0.2338421 0.3899009
6 health_c:high_incomeHigh 0.08128229 0.4670966 -0.8403726 1.0029372
t p
1 6.4438597 0.000000001022515
2 0.9581702 0.339255058705781
3 -0.8379728 0.403151146835412
4 5.3299832 0.000000289840230
5 0.4936786 0.622131695404800
6 0.1740160 0.862047345971205
Analysis of Variance Table
Model 1: life_exp ~ health_c + high_income + log2_gdp + urban10
Model 2: life_exp ~ health_c * high_income + log2_gdp + urban10
Res.Df RSS Df Sum of Sq F Pr(>F)
1 182 2488.2
2 181 2487.8 1 0.41621 0.0303 0.862health_c * high_income expands to the two main effects plus their interaction. The interaction coefficient is the difference in expenditure slopes between the groups. A significant interaction does not mean one subgroup has a significant slope and the other does not; the difference itself must be assessed. Limited overlap in expenditure across groups can make adjusted comparisons rely heavily on extrapolation.
Nonlinearity and transformations
21.10 Extension: a quadratic term
A model can be linear in its coefficients while nonlinear in a predictor. A quadratic term allows curvature:
\[E(Y)=\beta_0+\beta_1X_c+\beta_2X_c^2+\cdots.\]
reg$health_c <- reg$log2_health - mean(reg$log2_health)
linear_centred <- lm(life_exp ~ health_c + log2_gdp + urban10, data = reg)
quadratic_fit <- lm(life_exp ~ health_c + I(health_c^2) + log2_gdp + urban10,
data = reg)
anova(linear_centred, quadratic_fit)tab(coef_table(quadratic_fit))| term | estimate | se | lower95 | upper95 | t | p |
|---|---|---|---|---|---|---|
| (Intercept) | 41.142 | 5.962 | 29.379 | 52.905 | 6.901 | 0.000 |
| health_c | 0.385 | 0.410 | -0.425 | 1.195 | 0.938 | 0.349 |
| I(health_c^2) | -0.029 | 0.047 | -0.121 | 0.063 | -0.616 | 0.539 |
| log2_gdp | 2.465 | 0.464 | 1.550 | 3.380 | 5.313 | 0.000 |
| urban10 | 0.082 | 0.158 | -0.229 | 0.394 | 0.521 | 0.603 |
The slope now depends on \(X_c\): \(\beta_1+2\beta_2X_c\). Do not interpret \(\beta_1\) as one constant slope across the entire range. Polynomials can behave badly outside the observed range. A spline is another option, but flexibility must be justified and evaluated rather than added until p-values look attractive.
21.11 Extension: a log outcome changes the estimand
log_outcome_fit <- lm(log(life_exp) ~ log2_health + log2_gdp + urban10, data = reg)
log_outcome_table <- coef_table(log_outcome_fit)
tab(log_outcome_table)| term | estimate | se | lower95 | upper95 | t | p |
|---|---|---|---|---|---|---|
| (Intercept) | 3.794 | 0.039 | 3.717 | 3.872 | 96.822 | 0.000 |
| log2_health | 0.005 | 0.006 | -0.007 | 0.017 | 0.826 | 0.410 |
| log2_gdp | 0.035 | 0.007 | 0.021 | 0.048 | 5.171 | 0.000 |
| urban10 | 0.001 | 0.002 | -0.004 | 0.005 | 0.416 | 0.678 |
log_slope <- coef(log_outcome_fit)["log2_health"]
cat("Modelled multiplicative contrast in the conditional geometric mean for a doubling:",
fmt(exp(log_slope), 4), "\n")Modelled multiplicative contrast in the conditional geometric mean for a doubling: 1.0049 cat("Corresponding percentage contrast:", fmt(100 * (exp(log_slope) - 1)), "%\n")Corresponding percentage contrast: 0.49 %This model concerns the conditional mean of log life expectancy, or its exponentiated geometric-mean interpretation. Exponentiating a predicted log mean does not automatically give the arithmetic mean on the original scale. Back-transformation for arithmetic prediction requires attention to the residual distribution and possible smearing corrections. Life expectancy is not transformed here because it is inherently necessary; this is an explicit extension to teach scale changes.
22 Confidence intervals, prediction intervals and validation
22.1 Predict at observed covariate combinations
A confidence interval for a fitted mean concerns the modelled expected outcome at a predictor combination. A prediction interval concerns a new individual economy-level observation under the model and includes residual variation. It is not a patient-level prediction interval.
# Select actual rows near the 25th, 50th and 75th expenditure quantiles.
ord <- order(reg$log2_health)
positions <- pmax(1L, pmin(nrow(reg), round(c(.25, .50, .75) * nrow(reg))))
new_data <- reg[ord[positions], ]
mean_intervals <- predict(multiple_fit, newdata = new_data, interval = "confidence")
prediction_intervals <- predict(multiple_fit, newdata = new_data, interval = "prediction")
prediction_table <- data.frame(country_covariate_profile = new_data$country,
observed_life_exp = new_data$life_exp,
fitted = mean_intervals[, "fit"],
mean_lower95 = mean_intervals[, "lwr"],
mean_upper95 = mean_intervals[, "upr"],
prediction_lower95 = prediction_intervals[, "lwr"],
prediction_upper95 = prediction_intervals[, "upr"])
tab(prediction_table)| country_covariate_profile | observed_life_exp | fitted | mean_lower95 | mean_upper95 | prediction_lower95 | prediction_upper95 |
|---|---|---|---|---|---|---|
| Liberia | 62.163 | 64.419 | 62.891 | 65.947 | 56.967 | 71.871 |
| Botswana | 69.163 | 73.462 | 72.899 | 74.025 | 66.146 | 80.777 |
| Latvia | 75.427 | 77.921 | 77.234 | 78.609 | 70.595 | 85.248 |
These intervals use the conventional homoskedastic model, not the HC3 covariance extension. Robust coefficient SEs do not automatically produce valid heteroskedastic prediction intervals. The selected profiles are observed covariate combinations, not fabricated countries.
22.2 Extension: a held-out prediction exercise
Training performance is optimistic for evaluating prediction. We split the observed economies into training and test sets, fit only on the training set, and compare predictions with held-out outcomes. This is a teaching exercise; a random split may not reflect prediction to a new region or future year.
set.seed(params$seed + 4L)
train_id <- sample(seq_len(nrow(reg)), size = floor(.8 * nrow(reg)))
train <- reg[train_id, ]
test <- reg[-train_id, ]
train_fit <- lm(life_exp ~ log2_health + log2_gdp + urban10, data = train)
test_prediction <- predict(train_fit, newdata = test)
baseline_prediction <- rep(mean(train$life_exp), nrow(test))
rmse <- function(observed, predicted) sqrt(mean((observed - predicted)^2))
mae <- function(observed, predicted) mean(abs(observed - predicted))
tab(data.frame(model = c("Training mean baseline", "Multiple regression"),
test_n = nrow(test),
RMSE_years = c(rmse(test$life_exp, baseline_prediction),
rmse(test$life_exp, test_prediction)),
MAE_years = c(mae(test$life_exp, baseline_prediction),
mae(test$life_exp, test_prediction))))| model | test_n | RMSE_years | MAE_years |
|---|---|---|---|
| Training mean baseline | 38 | 6.826 | 5.210 |
| Multiple regression | 38 | 3.678 | 2.574 |
The baseline must use the training mean, not the test mean. Avoid using test outcomes to select predictors, tune transformations or choose a final model. One split is unstable. Repeated cross-validation can assess that variability, while geographic or temporal validation may better match the intended use.
Repeated measurements beyond two years
22.3 Extension: three matched years and Friedman test
The Friedman test compares repeated measurements using within-unit ranks. It is not an ordinary independent-groups Kruskal-Wallis test. We use three actual years and retain only economies observed at all three times.
middle_year <- as.integer(floor((baseline_year + analysis_year) / 2))
repeat_years <- unique(c(baseline_year, middle_year, analysis_year))
repeat_wide <- NULL
repeat_long <- NULL
if (length(repeat_years) == 3) {
pieces <- lapply(repeat_years, function(yr) {
z <- panel[panel$year == yr, c("iso3c", "life_exp")]
names(z)[2] <- paste0("y", yr)
z
})
repeat_wide <- Reduce(function(a, b) merge(a, b, by = "iso3c"), pieces)
repeat_wide <- repeat_wide[complete.cases(repeat_wide), ]
repeat_long <- do.call(rbind, lapply(repeat_years, function(yr) {
data.frame(iso3c = repeat_wide$iso3c, year = factor(yr, levels = repeat_years),
life_exp = repeat_wide[[paste0("y", yr)]])
}))
repeat_long$iso3c <- factor(repeat_long$iso3c)
if (nrow(repeat_wide) >= 5) {
print(friedman.test(life_exp ~ year | iso3c, data = repeat_long))
} else cat("Too few complete three-year profiles.\n")
} else cat("The chosen years do not provide three distinct time points.\n")
Friedman rank sum test
data: life_exp and year and iso3c
Friedman chi-squared = 339.31, df = 2, p-value < 2.2e-1622.4 Extension: classical repeated-measures ANOVA
if (!is.null(repeat_long) && nrow(repeat_wide) >= 5) {
repeated_fit <- aov(life_exp ~ year + Error(iso3c / year), data = repeat_long)
print(summary(repeated_fit))
}
Error: iso3c
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 216 40717 188.5
Error: iso3c:year
Df Sum Sq Mean Sq F value Pr(>F)
year 2 987.3 493.7 162.1 <2e-16 ***
Residuals 432 1315.6 3.0
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1This classical analysis is shown to explain the repeated-measures error structure. With more than two time points, the usual univariate F test involves a sphericity assumption. The code does not estimate Greenhouse-Geisser or Huynh-Feldt corrections, so do not treat its uncorrected p-value as automatically valid. Time intervals are not necessarily equal. Mixed-effects or generalised least-squares models with a justified covariance structure are more flexible extensions, especially for incomplete panels. A country’s repeated rows are never independent replicates merely because the dataset is long.
Permutation inference and its assumptions
22.5 Extension: a studentised two-group permutation exercise
Permutation tests rearrange labels under an exchangeability assumption. Income labels were not randomised, so this is not randomisation-based causal inference. A raw mean-difference permutation test can be problematic for a weak equal-means null with unequal variances. We therefore use a Welch-style studentised statistic, while still making the exchangeability limitation explicit.
welch_stat <- function(a, b) (mean(a) - mean(b)) /
sqrt(var(a) / length(a) + var(b) / length(b))
observed_stat <- welch_stat(low, high)
combined <- c(low, high)
n_low <- length(low)
set.seed(params$seed + 5L)
permuted_stat <- replicate(P, {
idx <- sample.int(length(combined), n_low, replace = FALSE)
welch_stat(combined[idx], combined[-idx])
})
permutation_p <- (1 + sum(abs(permuted_stat) >= abs(observed_stat))) / (P + 1)
tab(data.frame(observed_studentised_statistic = observed_stat,
Monte_Carlo_permutations = P, two_sided_p = permutation_p))| observed_studentised_statistic | Monte_Carlo_permutations | two_sided_p |
|---|---|---|
| -15.232 | 4999 | 0 |
The plus-one correction avoids reporting a zero Monte Carlo p-value. Studentisation can improve large-sample behaviour for unequal variances, but it does not turn observational group labels into randomised assignments or remove country dependence.
23 Power, sample size and practical importance
23.1 Planning a study is not interpreting an observed p-value
Power depends on the target difference, variability, sample size, significance level and design. A scientifically important difference should be chosen from context before data collection. The next example uses a hypothetical two-year difference and the observed pooled SD only as a planning illustration. It is not a recommendation for a clinical study or a claim that economies can be recruited at will.
planning_difference <- 2
power_plan <- power.t.test(delta = planning_difference, sd = pooled_sd,
sig.level = .05, power = .80,
type = "two.sample", alternative = "two.sided")
power_plan
Two-sample t test power calculation
n = 72.65903
delta = 2
sd = 4.273942
sig.level = 0.05
power = 0.8
alternative = two.sided
NOTE: n is number in *each* groupcat("Rounded-up planned observations per group:", ceiling(power_plan$n), "\n")Rounded-up planned observations per group: 73 n is per group for this two-sample calculation. It assumes the corresponding equal-variance t-test framework and equal group sizes. It is not a Welch-specific or clustered design calculation. See R’s power calculation documentation.
23.2 Paired planning uses the SD of differences
paired_plan <- power.t.test(delta = planning_difference, sd = sd(paired$change),
sig.level = .05, power = .80, type = "paired")
paired_plan
Paired t test power calculation
n = 18.03738
delta = 2
sd = 2.857546
sig.level = 0.05
power = 0.8
alternative = two.sided
NOTE: n is number of *pairs*, sd is std.dev. of *differences* within pairscat("Rounded-up planned pairs:", ceiling(paired_plan$n), "\n")Rounded-up planned pairs: 19 Do not substitute the SD of one year’s measurements for the SD of changes. Pairing can improve precision when measurements within a unit are positively related, but its benefit depends on that correlation and the design.
23.3 Extension: precision-based planning
For a rough normal-approximation interval half-width \(h\), \(n\approx(z_{.975}s/h)^2\). The example below uses a hypothetical one-year half-width and the observed country SD. It ignores finite-population and dependence corrections and is not an exact t-based design.
hypothetical_half_width <- 1
rough_n <- ceiling((qnorm(.975) * sd(life$life_exp) / hypothetical_half_width)^2)
tab(data.frame(target_half_width_years = hypothetical_half_width,
approximate_independent_observations = rough_n))| target_half_width_years | approximate_independent_observations |
|---|---|
| 1 | 195 |
Avoid ‘observed power’ calculated from the observed effect as a substitute for interpreting a confidence interval. A wide interval already communicates imprecision. A non-significant result is not proof that a practically important effect is absent.
23.4 Extension: equivalence is not non-significance
To argue that an effect is small enough to be practically negligible, specify an equivalence margin in advance and use an equivalence procedure. As a teaching illustration, a two-one-sided-tests approach at alpha .05 can be related to whether a 90% CI lies entirely within the pre-specified equivalence bounds. Here the hypothetical margin is plus or minus two years for the paired mean change. This margin is a teaching choice, not an established health standard.
equivalence_margin <- 2
paired_ci90 <- t.test(paired$change, conf.level = .90)$conf.int
within_margin <- paired_ci90[1] > -equivalence_margin &&
paired_ci90[2] < equivalence_margin
tab(data.frame(lower90 = paired_ci90[1], upper90 = paired_ci90[2],
lower_equivalence_bound = -equivalence_margin,
upper_equivalence_bound = equivalence_margin,
CI_entirely_within_bounds = within_margin))| lower90 | upper90 | lower_equivalence_bound | upper_equivalence_bound | CI_entirely_within_bounds |
|---|---|---|---|---|
| 2.695 | 3.336 | -2 | 2 | FALSE |
This is not the same as failing to reject a zero mean change. The paired t assumptions and the limitations of country-level inference still apply.
24 Building a complete analysis report
24.1 What every results section should contain
A complete report states the data source and retrieval date, analysis year, unit, exclusions, missingness, analysis-specific sample sizes, transformations, method, assumptions, estimates, intervals, p-values when relevant, sensitivity analyses and limitations. A reader should be able to identify exactly what was compared without reading the R code.
Do not write ‘health expenditure increases life expectancy’ from these models. Write ‘higher country-level expenditure was associated with higher or lower expected life expectancy under the fitted model’, using the actual coefficient direction. Avoid claiming that a country’s income category caused its observed outcome.
24.2 A compact inferential summary
inferential_summary <- data.frame(
analysis = c("One-sample mean versus hypothetical 70 years",
"Independent mean difference: Low minus High",
"Paired mean change: later minus earlier"),
n_units = c(nrow(life), nrow(independent), nrow(paired)),
estimate = c(mean(life$life_exp) - 70, mean(low) - mean(high), mean(paired$change)),
lower95 = c(one_test$conf.int[1] - 70, welch_test$conf.int[1],
paired_test$conf.int[1]),
upper95 = c(one_test$conf.int[2] - 70, welch_test$conf.int[2],
paired_test$conf.int[2]),
p = c(one_test$p.value, welch_test$p.value, paired_test$p.value))
tab(inferential_summary, "Different questions and denominators: all estimates are in years")| analysis | n_units | estimate | lower95 | upper95 | p |
|---|---|---|---|---|---|
| One-sample mean versus hypothetical 70 years | 217 | 3.791 | 2.838 | 4.744 | 0 |
| Independent mean difference: Low minus High | 111 | -15.605 | -17.682 | -13.528 | 0 |
| Paired mean change: later minus earlier | 217 | 3.016 | 2.633 | 3.398 | 0 |
These are separate teaching questions, not one pre-registered confirmatory family. If a real project treats them as a single family of confirmatory tests, specify and apply the relevant multiplicity strategy.
24.3 A coefficient plot with confidence intervals
plot_table <- multiple_table[multiple_table$term != "(Intercept)", ]
ypos <- seq_len(nrow(plot_table))
oldpar <- par(mar = c(5, 9, 3, 1))
plot(plot_table$estimate, ypos,
xlim = range(c(plot_table$lower95, plot_table$upper95, 0)),
ylim = c(.5, nrow(plot_table) + .5), yaxt = "n", pch = 16,
xlab = "Adjusted coefficient (life expectancy years)", ylab = "",
main = "Adjusted country-level associations")
segments(plot_table$lower95, ypos, plot_table$upper95, ypos)
axis(2, at = ypos, labels = c("Expenditure: doubling", "GDP: doubling",
"Urban share: +10 points"), las = 1)
abline(v = 0, lty = 2)
par(oldpar)
The three coefficients use different predictor units, so their absolute sizes are not a direct variable-importance ranking. Standardisation changes units but does not solve confounding or collinearity.
24.4 Export tables, data and models
write.csv(descriptives, "hss2381_outputs/descriptive_statistics.csv", row.names = FALSE)
write.csv(group_summary, "hss2381_outputs/income_group_descriptives.csv", row.names = FALSE)
write.csv(missing_table, "hss2381_outputs/missingness.csv", row.names = FALSE)
write.csv(inferential_summary, "hss2381_outputs/t_test_summary.csv", row.names = FALSE)
write.csv(pairwise_table, "hss2381_outputs/pairwise_welch_comparisons.csv", row.names = FALSE)
write.csv(multiple_table, "hss2381_outputs/regression_coefficients.csv", row.names = FALSE)
write.csv(robust_table, "hss2381_outputs/regression_HC3.csv", row.names = FALSE)
write.csv(cross, "hss2381_outputs/analysis_year_data.csv", row.names = FALSE)
write.csv(paired, "hss2381_outputs/paired_year_data.csv", row.names = FALSE)
write.csv(influence_table, "hss2381_outputs/influence_diagnostics.csv", row.names = FALSE)
saveRDS(list(simple = simple_fit, multiple = multiple_fit, quadratic = quadratic_fit,
anova = anova_fit), "hss2381_outputs/fitted_models.rds")
writeLines(c(paste("Retrieved UTC:", bundle$retrieved_utc),
paste("Analysis year:", analysis_year),
paste("Baseline year:", baseline_year),
"Income groups: current metadata at retrieval, not historical labels.",
"Source: World Bank Indicators API, source=2 (WDI).",
"Original contributing organisations are identified on indicator pages.",
"Request URLs:", bundle$requested_urls),
"hss2381_outputs/data_provenance.txt")
writeLines(capture.output(sessionInfo()), "hss2381_outputs/session_info.txt")
cat("Outputs saved in hss2381_outputs; raw download bundle retained in hss2381_cache.\n")Outputs saved in hss2381_outputs; raw download bundle retained in hss2381_cache.The export filenames are fixed and are overwritten on a new render. Archive a project folder or rename outputs if you need to preserve multiple versions. A self-contained HTML is convenient for reading; reproducibility still requires the source, cached data and software information.
25 Guided practical exercises
25.1 Lab 1: understand the dataset
Tasks
- Identify the observational unit and explain why it is not a patient.
- List each indicator’s unit and source page.
- Count economies with observed life expectancy and expenditure in the selected year.
- Explain why regional aggregate rows were removed.
- Explain why current income classifications cannot be described as historical classifications.
Starter analysis
lab1 <- data.frame(
year = analysis_year,
economies = nrow(cross),
life_observed = sum(!is.na(cross$life_exp)),
expenditure_observed = sum(!is.na(cross$health_pc)),
both_observed = sum(complete.cases(cross[c("life_exp", "health_pc")]))
)
tab(lab1)| year | economies | life_observed | expenditure_observed | both_observed |
|---|---|---|---|---|
| 2023 | 217 | 217 | 191 | 191 |
Interpretation check: the denominator for both_observed can be smaller than either single-indicator denominator. No observation is created to fill the gap.
25.2 Lab 2: describe and visualise a continuous variable
Tasks
- Summarise expenditure with mean, median, SD and IQR.
- Draw raw and log-scale histograms.
- Explain why the mean and median can differ.
- Identify one limitation of a boxplot without raw points.
- Explain why log transformation does not prove normality.
health_observed <- cross$health_pc[is.finite(cross$health_pc) & cross$health_pc > 0]
tab(data.frame(mean = mean(health_observed), median = median(health_observed),
sd = sd(health_observed), IQR = IQR(health_observed)))| mean | median | sd | IQR |
|---|---|---|---|
| 1433.837 | 470.403 | 2328.267 | 1517.316 |
Interpretation check: the SD is in dollars per capita; the IQR is also in dollars per capita. Variance would be in squared dollars per capita. Log-scale summaries have different units and cannot be relabelled as dollar summaries.
25.3 Lab 3: inference for one mean
Tasks
- Use a hypothetical 75-year reference instead of 70.
- State the null and alternative before reading the result.
- Calculate the difference and its CI.
- Explain what changes and what stays the same when the reference changes.
- Explain why the reference is not an estimated World Bank world mean.
lab_reference <- 75
lab_one <- t.test(life$life_exp, mu = lab_reference)
tab(data.frame(difference = mean(life$life_exp) - lab_reference,
lower95 = lab_one$conf.int[1] - lab_reference,
upper95 = lab_one$conf.int[2] - lab_reference,
p = lab_one$p.value))| difference | lower95 | upper95 | p |
|---|---|---|---|
| -1.209 | -2.162 | -0.256 | 0.013 |
Interpretation check: the observed mean and its SE do not change when only the reference changes. The difference, test statistic and p-value do change.
25.4 Lab 4: independent versus paired comparisons
Tasks
- Explain why low-income versus high-income economies are independent groups in this design.
- Explain why earlier versus later values within an economy are paired.
- Identify the direction of each reported contrast.
- Explain why matching by row number is unsafe.
- Explain why a paired temporal change is not a programme effect.
lab_independent <- t.test(low, high)
lab_paired <- t.test(paired$after, paired$before, paired = TRUE)
tab(data.frame(design = c("Independent: Low minus High", "Paired: later minus earlier"),
estimate = c(mean(low) - mean(high), mean(paired$change)),
se = c(lab_independent$stderr, lab_paired$stderr)))| design | estimate | se |
|---|---|---|
| Independent: Low minus High | -15.605 | 1.024 |
| Paired: later minus earlier | 3.016 | 0.194 |
Interpretation check: these analyses answer different questions on different sets of units. Comparing their p-values does not tell you which scientific relationship is more important.
25.5 Lab 5: ANOVA and follow-up comparisons
Tasks
- State the four-group omnibus null.
- Compare classical and Welch ANOVA assumptions.
- Explain why an omnibus rejection does not identify every differing pair.
- Interpret one pairwise mean difference and its direction.
- Distinguish Holm-adjusted p-values from individual confidence intervals.
summary(anova_fit) Df Sum Sq Mean Sq F value Pr(>F)
income_group 3 7037 2345.7 127.4 <2e-16 ***
Residuals 213 3922 18.4
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1welch_anova
One-way analysis of means (not assuming equal variances)
data: life_exp and income_group
F = 115.17, num df = 3.000, denom df = 81.929, p-value < 2.2e-16head(pairwise_table)Interpretation check: an omnibus test asks whether all means can be equal under the model. A pairwise contrast asks about a particular difference. Neither proves an income intervention would change life expectancy.
25.6 Lab 6: correlation and regression
Tasks
- Interpret the expenditure slope per doubling.
- Explain why a log2 predictor is not interpreted per one-dollar increase.
- Compare crude and adjusted slopes on the same rows.
- Inspect the residual plot and identify one concern, if present.
- Explain why an adjusted coefficient is not automatically causal.
lab_slope <- coef(multiple_fit)["log2_health"]
# A quadrupling is two doublings in the additive log2 model.
tab(data.frame(expenditure_contrast = c("Doubling", "Quadrupling"),
modelled_life_expectancy_difference = c(lab_slope, 2 * lab_slope)))| expenditure_contrast | modelled_life_expectancy_difference |
|---|---|
| Doubling | 0.392 |
| Quadrupling | 0.784 |
Interpretation check: this calculation assumes the model’s constant log-scale slope and unchanged other predictors. It does not guarantee that such a change is feasible or causal.
26 Capstone: a complete country-level health analysis
26.1 Research brief
Prepare a short report answering: How is country-level life expectancy associated with health expenditure, and how does the association change after accounting for GDP and urbanisation? Add one group comparison and one paired-year analysis to demonstrate the broader methods.
Your report should include:
- A clear ecological research question, target year and unit.
- A data dictionary, retrieval date and current-classification caveat.
- A missingness table and analysis-specific denominators.
- A descriptive table and at least three informative figures.
- A Welch two-group comparison with raw difference and CI.
- A paired comparison with matching explained.
- A four-group mean comparison with an appropriate follow-up strategy.
- Pearson and Spearman correlations with their different interpretations.
- Simple and multiple regression on the same rows.
- Residual, influence and collinearity diagnostics.
- At least one sensitivity analysis, clearly labelled exploratory.
- A conclusion separating association, practical importance and causation.
26.2 Suggested assessment rubric
| Component | Suggested share | Evidence of understanding |
|---|---|---|
| Question, design and provenance | 15% | Correct unit, year, source and target quantity |
| Data preparation and missingness | 15% | Reproducible exclusions and honest denominators |
| Description and visualisation | 15% | Appropriate summaries, units and readable figures |
| Inference and effect interpretation | 20% | Correct contrasts, assumptions, estimates and uncertainty |
| Regression and diagnostics | 20% | Correct coefficient units and substantive diagnostics |
| Reporting and limitations | 15% | No causal overclaiming; reproducible outputs |
This rubric is a proposed teaching tool, not an official grading policy.
26.3 Reusable analysis template
The following function is an additional reusable example. It fits a continuous-outcome model after complete-case selection, checks numerical variables, and returns the data, model, coefficient table and fit statistics. It does not automatically select predictors or establish that the model is scientifically appropriate.
analyse_continuous <- function(data, outcome, predictors) {
required <- c(outcome, predictors)
if (!all(required %in% names(data))) stop("One or more variables are absent.")
d <- data[complete.cases(data[required]), , drop = FALSE]
numeric_vars <- required[vapply(d[required], is.numeric, logical(1))]
if (length(numeric_vars)) {
keep <- apply(d[numeric_vars], 1, function(z) all(is.finite(z)))
d <- d[keep, , drop = FALSE]
}
if (!is.numeric(d[[outcome]])) stop("Outcome must be numerical.")
d <- droplevels(d)
model <- lm(reformulate(predictors, response = outcome), data = d)
if (df.residual(model) < 5) stop("Too few residual degrees of freedom.")
if (anyNA(coef(model))) stop("Rank-deficient model: inspect predictors.")
list(data = d, model = model, coefficients = coef_table(model),
fit = data.frame(n = nobs(model), R2 = summary(model)$r.squared,
adjusted_R2 = summary(model)$adj.r.squared,
residual_SD = sigma(model)))
}
# Change these names only after checking the new question and assumptions.
DV <- "life_exp"
IV <- c("log2_health", "log2_gdp", "urban10")
capstone <- analyse_continuous(cross, DV, IV)
tab(capstone$coefficients)| term | estimate | se | lower95 | upper95 | t | p |
|---|---|---|---|---|---|---|
| (Intercept) | 37.694 | 2.710 | 32.348 | 43.041 | 13.911 | 0.000 |
| log2_health | 0.392 | 0.409 | -0.416 | 1.200 | 0.958 | 0.339 |
| log2_gdp | 2.448 | 0.462 | 1.536 | 3.360 | 5.295 | 0.000 |
| urban10 | 0.090 | 0.157 | -0.220 | 0.400 | 0.570 | 0.569 |
tab(capstone$fit)| n | R2 | adjusted_R2 | residual_SD |
|---|---|---|---|
| 187 | 0.738 | 0.734 | 3.697 |
A reusable function improves consistency, but it does not make every combination of variables meaningful. Recheck units, transformations, missingness, functional form, overlap and design each time.
27 Worked conceptual answers and common errors
27.1 Short-answer checks
1. Why is a country’s life expectancy not a patient’s lifespan? Life expectancy is a population-level indicator derived from mortality conditions. The dataset contains economies, not individual follow-up records. A country-level relationship cannot be assigned directly to individuals.
2. Does a normal-looking histogram prove a t-test is valid? No. It says nothing about independence, selection bias, the target population or the correctness of matching. For paired tests, the relevant distribution is that of differences.
3. Is a standard error a smaller version of the SD? They measure different things. SD describes observational variation; SE describes uncertainty in an estimator under a model or design. The formula \(s/\sqrt n\) depends on independence and related assumptions.
4. Does p = .03 mean a 3% probability that the null is true? No. It is a tail probability under the null model and its assumptions. It is not a posterior probability for the null.
5. Does p > .05 prove two means are equal? No. The data may be imprecise. Examine the interval and consider whether it includes differences that matter. Equivalence requires a specified margin and a suitable procedure.
6. Why prefer Welch as a starting independent-groups comparison? It estimates variances separately and does not impose a common variance. It still requires appropriate units, independence and a defensible inferential framework.
7. Why is an ANOVA rejection not a complete answer? It says at least one mean differs under the model. It does not identify all pairs, quantify every contrast or establish causation.
8. Is Kruskal-Wallis always a test of medians? No. A median or location interpretation needs additional common-shape assumptions. Distributional differences can involve spread and shape.
9. Must regression predictors be normal? No. The standard normal-error assumption concerns the conditional errors, not the marginal predictor distributions.
10. What does the expenditure coefficient mean here? In the additive log2 model, it is a difference in expected country-level life expectancy associated with a doubling of expenditure, conditional on included predictors. It is not an effect per dollar and not automatically causal.
11. Does robust SE solve all assumption violations? No. HC3 addresses a covariance issue under independent-unit large-sample conditions. It does not fix dependence, confounding, omitted nonlinearities or selection bias.
12. Why can a prediction interval be wider than a confidence interval? Predicting a new outcome includes residual variation in addition to uncertainty about the fitted mean. Both intervals remain conditional on model assumptions.
13. Why not pool all years and run an ordinary t-test? Repeated measurements within an economy are dependent. Ignoring that structure can misstate uncertainty and changes the scientific question.
14. Why retain current income classifications as a caveat? The metadata retrieved today does not automatically reconstruct a country’s classification in the selected historical year. Historical classification analysis needs a separate dated source.
15. Does population weighting create a global life table? No. A weighted average of national life expectancy indicators does not reconstruct mortality schedules for a global life table.
27.2 Common errors and corrections
| Error | Better practice |
|---|---|
| Replace missing expenditure with zero | Keep missingness explicit and justify any imputation |
| Use latest available values for a fixed-year comparison | Request the exact year and report unavailable values |
| Include regional aggregates as if they were countries | Filter using country metadata before analysis |
| Match paired data by row position | Join by stable economy identifier |
| Choose one-sided testing after seeing the result | Pre-specify direction and rationale |
| Delete every boxplot outlier | Investigate source, context and influence |
| Use a normality p-value as the entire method decision | Consider design, estimand, graphs and robustness |
| Report only whether p is below .05 | Report estimate, interval, unit, n and limitations |
| Interpret log2 expenditure per dollar | Interpret one unit as a doubling |
| Call a 10-point urban-share difference a 10% increase | Say 10 percentage points |
| Interpret rank tests as universal median tests | State the rank/location estimand and assumptions |
| Compare nested models on different rows | Use a common complete-case dataset |
| Call HC3 a correction for country clustering | Model dependence explicitly if needed |
| Rank predictor importance by raw coefficient magnitude | Account for units, collinearity and the question |
| Claim adjusted regression proves causality | Describe conditional association and identification limits |
28 Formula and R command reference
| Concept | Formula or idea | Main R command |
|---|---|---|
| Mean | \(\sum x_i/n\) | mean(x) |
| Median | Middle ordered value | median(x) |
| Sample variance | \(\sum(x_i-\bar x)^2/(n-1)\) | var(x) |
| Standard deviation | Square root of variance | sd(x) |
| Interquartile range | \(Q_3-Q_1\) | IQR(x) |
| Standardisation | \((x-\bar x)/s\) | scale(x) |
| Mean SE | \(s/\sqrt n\) under independence | sd(x)/sqrt(length(x)) |
| Normal cumulative probability | Area to the left | pnorm() |
| Normal quantile | Inverse cumulative probability | qnorm() |
| t critical value | Quantile with degrees of freedom | qt() |
| One-sample t | Mean versus reference | t.test(x, mu = ...) |
| Welch t | Independent means, separate variances | t.test(x, y) |
| Pooled t | Independent means, common variance | t.test(x, y, var.equal = TRUE) |
| Paired t | One-sample test of matched differences | t.test(after, before, paired = TRUE) |
| Classical ANOVA | Between/within mean squares | aov(y ~ group) |
| Welch ANOVA | Unequal-variance mean comparison | oneway.test(y ~ group) |
| Tukey comparisons | Classical pairwise family | TukeyHSD() |
| Multiplicity | Defined family adjustment | p.adjust() |
| Rank-sum | Independent rank comparison | wilcox.test(x, y) |
| Signed-rank | Paired rank comparison | wilcox.test(after, before, paired = TRUE) |
| Kruskal-Wallis | Multi-group rank comparison | kruskal.test(y ~ group) |
| Friedman | Repeated rank comparison | friedman.test(y ~ time \| id) |
| Pearson correlation | Linear association | cor.test(x, y) |
| Spearman correlation | Rank/monotonic association | cor.test(x, y, method = "spearman") |
| Linear regression | Conditional mean model | lm(y ~ x1 + x2) |
| Coefficient interval | Parameter uncertainty | confint(model) |
| Prediction | Mean or new-observation interval | predict(model, interval = ...) |
| Influence | Change in fitted model | cooks.distance(model) |
| Power planning | Difference, SD, alpha and design | power.t.test() |
29 Troubleshooting and reproducibility checklist
29.1 If rendering fails
- Package unavailable: install the named package in the same R library used by Quarto. Restart R and render from the beginning.
- API timeout or HTTP error: check network access and the official API. A valid local cache can support later renders. The document never substitutes invented values.
- Unexpected API response: inspect the request and indicator metadata. A series may have changed, been retired or returned an error payload.
- Too few complete observations: inspect the selected year and missingness table. Do not silently change the year or relax the question.
- A group has no observations: inspect current metadata and group coverage. Do not create observations to make a test run.
- Object not found: run all earlier chunks in order or render in a fresh session.
- Singular model: inspect duplicated predictors, empty factor levels and perfect collinearity.
- Different results on another date: compare cached data, retrieval dates, parameter settings and package versions.
29.2 Before submitting a report
Sources and further reading
This workbook’s explanations and code are original teaching material. The links below document the course scope, data and software methods; no course slides or proprietary teaching notes are reproduced.
- University of Ottawa: Health Sciences course catalogue, including HSS 2381.
- World Bank: About the Indicators API.
- World Bank: Country API queries and metadata.
- World Bank: Aggregate API queries.
- World Bank: Life expectancy at birth, total.
- World Bank: Current health expenditure per capita.
- World Bank: GDP per capita, current US dollars.
- World Bank: Urban population percentage.
- World Bank: Total population.
- R documentation: t-tests.
- R documentation: one-way mean tests, including Welch ANOVA.
- R documentation: Wilcoxon tests.
- R documentation: Kruskal-Wallis test.
- R documentation: correlation tests.
- R documentation: linear models.
- R documentation: p-value adjustment.
- R documentation: t-test power calculations.
- Quarto: using R and rendering executable documents.
- OpenStax: Introductory Statistics 2e, a free companion textbook.
Indicator pages identify their original contributing organisations and applicable licences. Retain attribution when redistributing downloaded data. Software documentation can describe newer versions than the installed R; the session information below records the actual local environment.
Session information
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
loaded via a namespace (and not attached):
[1] digest_0.6.39 fastmap_1.2.0 xfun_0.55 knitr_1.51
[5] htmltools_0.5.9 rmarkdown_2.30 cli_3.6.5 vctrs_0.6.5
[9] compiler_4.5.2 rstudioapi_0.17.1 tools_4.5.2 curl_7.0.0
[13] evaluate_1.0.5 yaml_2.3.12 otel_0.2.0 rlang_1.3.0
[17] jsonlite_2.0.0 htmlwidgets_1.6.4