options(repos = c(CRAN = "https://cloud.r-project.org"))
needed <- c("dplyr","plm","lmtest","sandwich","ggplot2","scales",
"knitr","rmarkdown","gapminder")
to_install <- needed[!needed %in% installed.packages()[,"Package"]]
if (length(to_install)) install.packages(to_install)
library(dplyr); library(plm); library(ggplot2)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
##
## Attaching package: 'plm'
## The following objects are masked from 'package:dplyr':
##
## between, lag, lead
library(dplyr)
panel_raw <- tryCatch({
# --- try the R package first ---
if (!requireNamespace("gapminder", quietly = TRUE)) {
install.packages("gapminder", repos = "https://cloud.r-project.org")
}
library(gapminder)
gapminder::gapminder %>% as_tibble()
}, error = function(e) {
# --- fall back to the online CSV mirror of the same dataset ---
message("gapminder package unavailable -- pulling the same data online instead.")
read.csv(
"https://raw.githubusercontent.com/vincentarelbundock/Rdatasets/master/csv/gapminder/gapminder.csv",
stringsAsFactors = FALSE
) %>% select(-rownames) %>% as_tibble()
})
panel_raw <- panel_raw %>%
rename(gdp_pc = gdpPercap, life_exp = lifeExp, pop = pop) %>%
mutate(country = as.character(country), continent = as.character(continent))
write.csv(panel_raw, "data/panel_raw.csv", row.names = FALSE)
cat("Route A loaded:", nrow(panel_raw), "country-year observations,",
n_distinct(panel_raw$country), "countries, years",
min(panel_raw$year), "-", max(panel_raw$year), "\n")
## Route A loaded: 1704 country-year observations, 142 countries, years 1952 - 2007
library(dplyr)
library(plm)
panel_raw <- read.csv("data/panel_raw.csv", stringsAsFactors = FALSE)
panel <- panel_raw %>%
filter(!is.na(gdp_pc), !is.na(life_exp), gdp_pc > 0) %>%
arrange(country, year) %>%
mutate(
log_gdp_pc = log(gdp_pc),
log_gdp_pc_sq = log_gdp_pc^2, # tests diminishing returns
log_pop = log(pop),
period = ifelse(year < 1980, "1952-1977 (early)", "1982-2007 (recent)")
) %>%
group_by(country) %>%
mutate(
life_exp_growth = life_exp - lag(life_exp),
gdp_growth_pct = (gdp_pc - lag(gdp_pc)) / lag(gdp_pc) * 100
) %>%
ungroup()
waves_per_country <- panel %>% count(country) %>% pull(n)
n_waves <- max(waves_per_country)
balanced_countries <- panel %>% count(country) %>% filter(n == n_waves) %>% pull(country)
panel_balanced <- panel %>% filter(country %in% balanced_countries)
cat(sprintf(
"Cleaned panel: %d obs, %d countries (%d balanced, %d dropped for imbalance)\n",
nrow(panel), n_distinct(panel$country),
length(balanced_countries), n_distinct(panel$country) - length(balanced_countries)
))
## Cleaned panel: 1704 obs, 142 countries (142 balanced, 0 dropped for imbalance)
write.csv(panel, "data/panel_clean.csv", row.names = FALSE)
write.csv(panel_balanced, "data/panel_balanced.csv", row.names = FALSE)
panel_balanced$year_num <- panel_balanced$year
pdata <- pdata.frame(panel_balanced, index = c("country", "year"))
saveRDS(pdata, "data/pdata.rds")
library(plm)
library(lmtest)
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(sandwich)
library(dplyr)
panel <- read.csv("data/panel_clean.csv", stringsAsFactors = FALSE)
pdata <- readRDS("data/pdata.rds")
cluster_se <- function(model) coeftest(model, vcov = vcovHC(model, type = "HC1", cluster = "group"))
## ---- helper: pull estimate/SE/p-value from any coefficient matrix ----
## (works for lm summary, plm summary, and coeftest -- all have columns
## in the order: Estimate, Std.Error, test-stat, p-value)
get_coef <- function(mat, term) {
row <- mat[term, ]
list(est = unname(row[1]), se = unname(row[2]), p = unname(row[4]))
}
sig_stars <- function(p) {
if (p < 0.001) "***" else if (p < 0.01) "**" else if (p < 0.05) "*" else if (p < 0.1) "." else "(n.s.)"
}
sig_word <- function(p) if (p < 0.05) "statistically significant" else "not statistically significant"
# M1: Pooled OLS
m1 <- lm(life_exp ~ log_gdp_pc, data = panel)
# M2: Pooled OLS, quadratic (tests concavity)
m2 <- lm(life_exp ~ log_gdp_pc + log_gdp_pc_sq, data = panel)
# M3: Two-way Fixed Effects (country + year)
m3 <- plm(life_exp ~ log_gdp_pc + log_gdp_pc_sq + log_pop,
data = pdata, model = "within", effect = "twoways")
# M4: Random Effects (for Hausman comparison)
m4 <- plm(life_exp ~ log_gdp_pc + log_gdp_pc_sq + log_pop,
data = pdata, model = "random")
hausman <- phtest(m3, m4)
# M5: FE with continent interaction (regional heterogeneity)
m5 <- plm(life_exp ~ log_gdp_pc * continent + log_gdp_pc_sq + log_pop,
data = pdata, model = "within", effect = "twoways")
# M6: Has the gradient flattened? Early vs recent era
m6_early <- plm(life_exp ~ log_gdp_pc + log_gdp_pc_sq + log_pop,
data = pdata[as.numeric(as.character(pdata$year_num)) < 1980, ], model = "within")
m6_recent <- plm(life_exp ~ log_gdp_pc + log_gdp_pc_sq + log_pop,
data = pdata[as.numeric(as.character(pdata$year_num)) >= 1982, ], model = "within")
# ==============================================================================
# ---- Print everything, with automatic interpretation after each model ----
# ==============================================================================
## M1
cat("=== M1: Pooled OLS ===\n"); print(summary(m1))
## === M1: Pooled OLS ===
##
## Call:
## lm(formula = life_exp ~ log_gdp_pc, data = panel)
##
## Residuals:
## Min 1Q Median 3Q Max
## -32.778 -4.204 1.212 4.658 19.285
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -9.1009 1.2277 -7.413 1.93e-13 ***
## log_gdp_pc 8.4051 0.1488 56.500 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.62 on 1702 degrees of freedom
## Multiple R-squared: 0.6522, Adjusted R-squared: 0.652
## F-statistic: 3192 on 1 and 1702 DF, p-value: < 2.2e-16
c1 <- get_coef(summary(m1)$coefficients, "log_gdp_pc")
cat(sprintf(
"\n--- INTERPRETATION: M1 ---\nEach 1-log-point rise in GDP per capita is associated with %.2f more years of\nlife expectancy %s (p = %.4f, %s). Note: this is a naive pooled estimate with\nno control for confounding country characteristics -- see M3 for a cleaner estimate.\n",
c1$est, sig_word(c1$p), c1$p, sig_stars(c1$p)
))
##
## --- INTERPRETATION: M1 ---
## Each 1-log-point rise in GDP per capita is associated with 8.41 more years of
## life expectancy statistically significant (p = 0.0000, ***). Note: this is a naive pooled estimate with
## no control for confounding country characteristics -- see M3 for a cleaner estimate.
## M2
cat("\n=== M2: Pooled OLS, quadratic ===\n"); print(summary(m2))
##
## === M2: Pooled OLS, quadratic ===
##
## Call:
## lm(formula = life_exp ~ log_gdp_pc + log_gdp_pc_sq, data = panel)
##
## Residuals:
## Min 1Q Median 3Q Max
## -30.160 -4.150 1.277 4.476 20.961
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -36.4594 7.8704 -4.632 3.89e-06 ***
## log_gdp_pc 15.2089 1.9392 7.843 7.71e-15 ***
## log_gdp_pc_sq -0.4134 0.1175 -3.519 0.000445 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.594 on 1701 degrees of freedom
## Multiple R-squared: 0.6548, Adjusted R-squared: 0.6544
## F-statistic: 1613 on 2 and 1701 DF, p-value: < 2.2e-16
c2_lin <- get_coef(summary(m2)$coefficients, "log_gdp_pc")
c2_sq <- get_coef(summary(m2)$coefficients, "log_gdp_pc_sq")
cat(sprintf(
"\n--- INTERPRETATION: M2 ---\nThe squared term is %.3f (p = %.4f, %s), which is %s and %s.\n%s\n",
c2_sq$est, c2_sq$p, sig_stars(c2_sq$p), sig_word(c2_sq$p),
ifelse(c2_sq$est < 0, "negative", "positive"),
ifelse(c2_sq$est < 0 && c2_sq$p < 0.05,
"=> Confirms a concave (diminishing-returns) relationship: income buys progressively\n less additional life expectancy as income rises.",
"=> No clear evidence of diminishing returns in this specification.")
))
##
## --- INTERPRETATION: M2 ---
## The squared term is -0.413 (p = 0.0004, ***), which is statistically significant and negative.
## => Confirms a concave (diminishing-returns) relationship: income buys progressively
## less additional life expectancy as income rises.
## M3
cat("\n=== M3: Two-way Fixed Effects (clustered SE) ===\n")
##
## === M3: Two-way Fixed Effects (clustered SE) ===
fe3 <- cluster_se(m3); print(fe3)
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## log_gdp_pc 10.60874 4.04170 2.6248 0.008755 **
## log_gdp_pc_sq -0.45675 0.23116 -1.9759 0.048342 *
## log_pop 7.43313 1.20165 6.1858 7.893e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
c3_lin <- get_coef(fe3, "log_gdp_pc")
c3_sq <- get_coef(fe3, "log_gdp_pc_sq")
cat(sprintf(
"\n--- INTERPRETATION: M3 ---\nOnce country and year fixed effects absorb time-invariant confounders\n(institutions, geography) and common global shocks, a 1-log-point income rise is\nassociated with %.2f more years of life expectancy %s (p = %.4f, %s).\nThe quadratic term (%.3f, p = %.4f) is %s: diminishing returns %s after removing\nconfounding.\n",
c3_lin$est, sig_word(c3_lin$p), c3_lin$p, sig_stars(c3_lin$p),
c3_sq$est, c3_sq$p, sig_stars(c3_sq$p),
ifelse(c3_sq$est < 0 && c3_sq$p < 0.05, "survive", "do not clearly survive")
))
##
## --- INTERPRETATION: M3 ---
## Once country and year fixed effects absorb time-invariant confounders
## (institutions, geography) and common global shocks, a 1-log-point income rise is
## associated with 10.61 more years of life expectancy statistically significant (p = 0.0088, **).
## The quadratic term (-0.457, p = 0.0483) is *: diminishing returns survive after removing
## confounding.
## M4 + Hausman
cat("\n=== M4: Random Effects ===\n"); print(summary(m4))
##
## === M4: Random Effects ===
## Oneway (individual) effect Random Effect Model
## (Swamy-Arora's transformation)
##
## Call:
## plm(formula = life_exp ~ log_gdp_pc + log_gdp_pc_sq + log_pop,
## data = pdata, model = "random")
##
## Balanced Panel: n = 142, T = 12, N = 1704
##
## Effects:
## var std.dev share
## idiosyncratic 10.429 3.229 0.255
## individual 30.454 5.519 0.745
## theta: 0.8334
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -22.33649 -2.38604 0.13006 2.56770 12.81686
##
## Coefficients:
## Estimate Std. Error z-value Pr(>|z|)
## (Intercept) -152.49865 7.99466 -19.0751 < 2.2e-16 ***
## log_gdp_pc 16.11904 1.85089 8.7088 < 2.2e-16 ***
## log_gdp_pc_sq -0.55404 0.10827 -5.1171 3.103e-07 ***
## log_pop 7.49674 0.22319 33.5891 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 79804
## Residual Sum of Squares: 26530
## R-Squared: 0.66756
## Adj. R-Squared: 0.66697
## Chisq: 3413.66 on 3 DF, p-value: < 2.22e-16
c4_lin <- get_coef(summary(m4)$coefficients, "log_gdp_pc")
cat("\n=== Hausman test (FE vs RE) ===\n"); print(hausman)
##
## === Hausman test (FE vs RE) ===
##
## Hausman Test
##
## data: life_exp ~ log_gdp_pc + log_gdp_pc_sq + log_pop
## chisq = 29952, df = 3, p-value < 2.2e-16
## alternative hypothesis: one model is inconsistent
pct_overstate <- round((c4_lin$est - c3_lin$est) / c4_lin$est * 100, 1)
cat(sprintf(
"\n--- INTERPRETATION: M4 & Hausman ---\nHausman chi-sq = %.2f, p = %.4f.\n%s\nRandom Effects gives a coefficient of %.2f vs %.2f under Fixed Effects -- RE is\n%.1f%% %s than FE on this term.\n",
unname(hausman$statistic), hausman$p.value,
ifelse(hausman$p.value < 0.05,
"=> p < 0.05: reject the null that RE is consistent. Country effects ARE correlated\n with income, so Fixed Effects is the appropriate, unbiased specification.",
"=> p >= 0.05: fail to reject RE; Random Effects is not rejected as inconsistent here."),
c4_lin$est, c3_lin$est, abs(pct_overstate), ifelse(pct_overstate > 0, "higher", "lower")
))
##
## --- INTERPRETATION: M4 & Hausman ---
## Hausman chi-sq = 29952.40, p = 0.0000.
## => p < 0.05: reject the null that RE is consistent. Country effects ARE correlated
## with income, so Fixed Effects is the appropriate, unbiased specification.
## Random Effects gives a coefficient of 16.12 vs 10.61 under Fixed Effects -- RE is
## 34.2% higher than FE on this term.
## M5
cat("\n=== M5: FE with continent interaction ===\n")
##
## === M5: FE with continent interaction ===
fe5 <- cluster_se(m5); print(fe5)
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## log_gdp_pc 12.02592 3.72853 3.2254 0.001284 **
## log_gdp_pc_sq -0.69225 0.23992 -2.8853 0.003965 **
## log_pop 6.85415 1.19752 5.7236 1.251e-08 ***
## log_gdp_pc:continentAmericas 2.46433 1.73026 1.4243 0.154575
## log_gdp_pc:continentAsia 3.95055 1.48208 2.6655 0.007767 **
## log_gdp_pc:continentEurope 2.25863 1.54630 1.4607 0.144309
## log_gdp_pc:continentOceania -0.10825 1.80012 -0.0601 0.952054
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
interact_terms <- grep("log_gdp_pc:continent", rownames(fe5), value = TRUE)
cat("\n--- INTERPRETATION: M5 (regional heterogeneity) ---\n")
##
## --- INTERPRETATION: M5 (regional heterogeneity) ---
if (length(interact_terms) == 0) {
cat("No interaction terms were estimated (check that 'continent' has >1 level in this sample).\n")
} else {
for (term in interact_terms) {
region <- sub(".*continent", "", term)
ci <- get_coef(fe5, term)
cat(sprintf(
"%s: income-health slope differs from the reference region by %.2f (p = %.4f, %s) -- %s.\n",
region, ci$est, ci$p, sig_stars(ci$p), sig_word(ci$p)
))
}
cat("=> A significant interaction means the income-health gradient is not uniform\n across regions -- growth-based health strategies should be regionally tailored.\n")
}
## Americas: income-health slope differs from the reference region by 2.46 (p = 0.1546, (n.s.)) -- not statistically significant.
## Asia: income-health slope differs from the reference region by 3.95 (p = 0.0078, **) -- statistically significant.
## Europe: income-health slope differs from the reference region by 2.26 (p = 0.1443, (n.s.)) -- not statistically significant.
## Oceania: income-health slope differs from the reference region by -0.11 (p = 0.9521, (n.s.)) -- not statistically significant.
## => A significant interaction means the income-health gradient is not uniform
## across regions -- growth-based health strategies should be regionally tailored.
## M6a / M6b
cat("\n=== M6a: FE, early era (1952-1977) ===\n")
##
## === M6a: FE, early era (1952-1977) ===
fe6a <- cluster_se(m6_early); print(fe6a)
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## log_gdp_pc 12.33422 3.54778 3.4766 0.0005389 ***
## log_gdp_pc_sq -0.47913 0.20527 -2.3342 0.0198647 *
## log_pop 12.73931 1.06394 11.9737 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cat("\n=== M6b: FE, recent era (1982-2007) ===\n")
##
## === M6b: FE, recent era (1982-2007) ===
fe6b <- cluster_se(m6_recent); print(fe6b)
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## log_gdp_pc -0.44267 4.75331 -0.0931 0.9258
## log_gdp_pc_sq 0.23298 0.26354 0.8840 0.3770
## log_pop 7.83796 1.14283 6.8584 1.521e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
c6a <- get_coef(fe6a, "log_gdp_pc")
c6b <- get_coef(fe6b, "log_gdp_pc")
pct_shrink <- round((c6a$est - c6b$est) / c6a$est * 100, 1)
cat(sprintf(
"\n--- INTERPRETATION: M6 (era comparison) ---\nEarly era (1952-1977) coefficient: %.2f (p = %.4f, %s)\nRecent era (1982-2007) coefficient: %.2f (p = %.4f, %s)\n=> The income-health payoff has %s by %.1f%% between the two eras.\n%s\n",
c6a$est, c6a$p, sig_stars(c6a$p),
c6b$est, c6b$p, sig_stars(c6b$p),
ifelse(pct_shrink > 0, "shrunk", "grown"), abs(pct_shrink),
ifelse(pct_shrink > 0,
"Policy implication: as the world has gotten richer, growth alone now buys less\nhealth than it used to -- direct health-system investment becomes relatively more\nimportant at the margin.",
"Policy implication: the income-health link has strengthened over time in this sample.")
))
##
## --- INTERPRETATION: M6 (era comparison) ---
## Early era (1952-1977) coefficient: 12.33 (p = 0.0005, ***)
## Recent era (1982-2007) coefficient: -0.44 (p = 0.9258, (n.s.))
## => The income-health payoff has shrunk by 103.6% between the two eras.
## Policy implication: as the world has gotten richer, growth alone now buys less
## health than it used to -- direct health-system investment becomes relatively more
## important at the margin.
saveRDS(list(m1=m1, m2=m2, m3=m3, m4=m4, m5=m5,
m6_early=m6_early, m6_recent=m6_recent, hausman=hausman),
"output/models.rds")
cat("\n================ OVERALL SUMMARY ================\n")
##
## ================ OVERALL SUMMARY ================
cat(sprintf(
"Hausman p-value: %.4f -> %s\nFixed-effects income coefficient: %.2f %s\nDiminishing returns (quadratic term): %s\nGradient trend over time: %s (%.1f%% change)\n",
hausman$p.value, ifelse(hausman$p.value < 0.05, "use FIXED EFFECTS", "RE not rejected"),
c3_lin$est, sig_stars(c3_lin$p),
ifelse(c3_sq$est < 0 && c3_sq$p < 0.05, "confirmed", "not confirmed"),
ifelse(pct_shrink > 0, "weakening", "strengthening"), abs(pct_shrink)
))
## Hausman p-value: 0.0000 -> use FIXED EFFECTS
## Fixed-effects income coefficient: 10.61 **
## Diminishing returns (quadratic term): confirmed
## Gradient trend over time: weakening (103.6% change)
library(ggplot2)
library(dplyr)
library(plm) # needed so summary() dispatches correctly on saved plm objects
panel <- read.csv("data/panel_clean.csv", stringsAsFactors = FALSE)
models <- readRDS("output/models.rds")
un_theme <- theme_minimal(base_size = 13) +
theme(
plot.title = element_text(face = "bold", size = 16, margin = margin(b = 4)),
plot.subtitle = element_text(color = "grey35", size = 11, margin = margin(b = 12)),
plot.caption = element_text(color = "grey50", size = 8.5, hjust = 0),
panel.grid.minor = element_blank(),
legend.position = "bottom",
strip.text = element_text(face = "bold")
)
un_palette <- c("Africa" = "#E4572E", "Americas" = "#17BEBB", "Asia" = "#FFC914",
"Europe" = "#2E86AB", "Oceania" = "#A23B72")
# ==============================================================================
# Chart 1: Preston Curve
# ==============================================================================
latest <- panel %>% filter(year == max(year))
p1 <- ggplot(latest, aes(x = gdp_pc, y = life_exp, size = pop, color = continent)) +
geom_point(alpha = 0.75) +
geom_smooth(aes(x = gdp_pc, y = life_exp), method = "loess", se = FALSE,
color = "grey20", linewidth = 0.9, inherit.aes = FALSE) +
scale_x_log10(labels = scales::dollar_format()) +
scale_size_continuous(range = c(1, 16), guide = "none") +
scale_color_manual(values = un_palette, name = NULL) +
labs(title = "The Preston Curve: Income and Life Expectancy, 2007",
subtitle = "Each bubble is a country (size = population). Health gains from income growth diminish at higher incomes.",
x = "GDP per capita (current US$, log scale)", y = "Life expectancy at birth (years)",
caption = "Source: Gapminder Foundation (UN Population Division, UN Statistics Division & World Bank).") + un_theme
print(p1)
## `geom_smooth()` using formula = 'y ~ x'
# ggsave("output/01_preston_curve.png", p1, width = 9, height = 6, dpi = 200)
## ---- Automatic interpretation: Chart 1 ----
med_gdp <- median(latest$gdp_pc)
low_half <- latest %>% filter(gdp_pc <= med_gdp)
high_half <- latest %>% filter(gdp_pc > med_gdp)
slope_low <- coef(lm(life_exp ~ log(gdp_pc), data = low_half))["log(gdp_pc)"]
slope_high <- coef(lm(life_exp ~ log(gdp_pc), data = high_half))["log(gdp_pc)"]
pct_flatter <- round((1 - slope_high / slope_low) * 100, 1)
cat("\n--- INTERPRETATION: Chart 1 (Preston Curve) ---\n")
##
## --- INTERPRETATION: Chart 1 (Preston Curve) ---
cat(sprintf(
"Among the poorer half of countries (GDP/cap below $%s), each log-point of income is\nassociated with %.1f years more life expectancy.\n",
format(round(med_gdp), big.mark = ","), slope_low
))
## Among the poorer half of countries (GDP/cap below $6,124), each log-point of income is
## associated with 8.4 years more life expectancy.
cat(sprintf(
"Among the richer half, the same log-point income gain buys only %.1f years -\na slope %.1f%% flatter than the poorer half.\n",
slope_high, pct_flatter
))
## Among the richer half, the same log-point income gain buys only 5.6 years -
## a slope 33.4% flatter than the poorer half.
cat("=> Confirms the classic concave Preston Curve: income matters far more for health\n at low income levels than at high ones.\n")
## => Confirms the classic concave Preston Curve: income matters far more for health
## at low income levels than at high ones.
# ==============================================================================
# Chart 2: Regional trends
# ==============================================================================
p2 <- panel %>%
group_by(continent, year) %>%
summarise(life_exp = weighted.mean(life_exp, pop), .groups = "drop") %>%
ggplot(aes(x = year, y = life_exp, color = continent)) +
geom_line(linewidth = 1.1) + geom_point(size = 1.6) +
scale_color_manual(values = un_palette, name = NULL) +
labs(title = "Life Expectancy Gains by Region, 1952-2007",
subtitle = "Population-weighted regional averages", x = NULL,
y = "Life expectancy at birth (years)",
caption = "Source: Gapminder Foundation (UN & World Bank series).") + un_theme
print(p2)
# ggsave("output/02_regional_trends.png", p2, width = 9, height = 6, dpi = 200)
## ---- Automatic interpretation: Chart 2 ----
regional_gains <- panel %>%
group_by(continent, year) %>%
summarise(life_exp = weighted.mean(life_exp, pop), .groups = "drop") %>%
group_by(continent) %>%
summarise(
start_val = life_exp[which.min(year)],
end_val = life_exp[which.max(year)],
gain = end_val - start_val,
.groups = "drop"
) %>%
arrange(desc(gain))
cat("\n--- INTERPRETATION: Chart 2 (Regional Trends) ---\n")
##
## --- INTERPRETATION: Chart 2 (Regional Trends) ---
for (i in seq_len(nrow(regional_gains))) {
r <- regional_gains[i, ]
cat(sprintf("%-10s: %.1f -> %.1f years (+%.1f)\n", r$continent, r$start_val, r$end_val, r$gain))
}
## Asia : 42.9 -> 69.4 years (+26.5)
## Africa : 38.8 -> 54.6 years (+15.8)
## Americas : 60.2 -> 75.4 years (+15.1)
## Europe : 64.9 -> 77.9 years (+13.0)
## Oceania : 69.2 -> 81.1 years (+11.9)
cat(sprintf(
"=> %s gained the most life expectancy (+%.1f years); %s gained the least (+%.1f years).\n",
regional_gains$continent[1], regional_gains$gain[1],
regional_gains$continent[nrow(regional_gains)], regional_gains$gain[nrow(regional_gains)]
))
## => Asia gained the most life expectancy (+26.5 years); Oceania gained the least (+11.9 years).
# ==============================================================================
# Chart 3: Coefficient comparison
# ==============================================================================
coef_df <- data.frame(
model = c("Pooled OLS","Pooled OLS","Fixed Effects","Fixed Effects","Random Effects","Random Effects"),
term = rep(c("log(GDP per capita)", "log(GDP per capita)^2"), 3),
estimate = c(coef(models$m2)["log_gdp_pc"], coef(models$m2)["log_gdp_pc_sq"],
coef(models$m3)["log_gdp_pc"], coef(models$m3)["log_gdp_pc_sq"],
coef(models$m4)["log_gdp_pc"], coef(models$m4)["log_gdp_pc_sq"]),
se = c(summary(models$m2)$coefficients["log_gdp_pc",2], summary(models$m2)$coefficients["log_gdp_pc_sq",2],
summary(models$m3)$coefficients["log_gdp_pc",2], summary(models$m3)$coefficients["log_gdp_pc_sq",2],
summary(models$m4)$coefficients["log_gdp_pc",2], summary(models$m4)$coefficients["log_gdp_pc_sq",2])
) %>% mutate(model = factor(model, levels = c("Pooled OLS","Random Effects","Fixed Effects")))
p3 <- ggplot(coef_df, aes(x = estimate, y = model, color = model)) +
geom_vline(xintercept = 0, linetype = "dashed", color = "grey60") +
geom_pointrange(aes(xmin = estimate - 1.96*se, xmax = estimate + 1.96*se), linewidth = 0.9, size = 0.7) +
facet_wrap(~term, scales = "free_x") +
scale_color_manual(values = c("Pooled OLS"="#E4572E","Random Effects"="#FFC914","Fixed Effects"="#2E86AB")) +
labs(title = "Naive Pooled Estimates Overstate the Income-Health Link",
subtitle = "95% CI. Country & year fixed effects absorb confounders (institutions, geography, global shocks).",
x = "Coefficient estimate", y = NULL,
caption = "Fixed-effects standard errors clustered by country.") +
un_theme + theme(legend.position = "none", panel.spacing = unit(2.2, "lines"))
print(p3)
# ggsave("output/03_coefficient_comparison.png", p3, width = 9, height = 5, dpi = 200)
## ---- Automatic interpretation: Chart 3 ----
ols_est <- coef_df %>% filter(model == "Pooled OLS", term == "log(GDP per capita)") %>% pull(estimate)
fe_est <- coef_df %>% filter(model == "Fixed Effects", term == "log(GDP per capita)") %>% pull(estimate)
re_est <- coef_df %>% filter(model == "Random Effects", term == "log(GDP per capita)") %>% pull(estimate)
pct_overstate <- round((ols_est - fe_est) / ols_est * 100, 1)
quad_sign <- ifelse(coef(models$m3)["log_gdp_pc_sq"] < 0, "negative", "positive")
cat("\n--- INTERPRETATION: Chart 3 (Model Comparison) ---\n")
##
## --- INTERPRETATION: Chart 3 (Model Comparison) ---
cat(sprintf("Pooled OLS coefficient: %.2f\n", ols_est))
## Pooled OLS coefficient: 15.21
cat(sprintf("Random Effects coeff.: %.2f\n", re_est))
## Random Effects coeff.: 16.12
cat(sprintf("Fixed Effects coeff.: %.2f\n", fe_est))
## Fixed Effects coeff.: 10.61
cat(sprintf(
"=> Pooled OLS overstates the true income-health link by about %.1f%% relative to\n Fixed Effects, once time-invariant country traits (institutions, geography) are\n controlled for. The quadratic term is %s in every specification, confirming\n diminishing returns survive after removing confounding.\n",
pct_overstate, quad_sign
))
## => Pooled OLS overstates the true income-health link by about 30.2% relative to
## Fixed Effects, once time-invariant country traits (institutions, geography) are
## controlled for. The quadratic term is negative in every specification, confirming
## diminishing returns survive after removing confounding.
# ==============================================================================
# Chart 4: Has the gradient flattened over time?
# ==============================================================================
period_df <- data.frame(
period = c("Early era\n(1952-1977)", "Recent era\n(1982-2007)"),
estimate = c(coef(models$m6_early)["log_gdp_pc"], coef(models$m6_recent)["log_gdp_pc"]),
se = c(summary(models$m6_early)$coefficients["log_gdp_pc",2],
summary(models$m6_recent)$coefficients["log_gdp_pc",2])
)
p4 <- ggplot(period_df, aes(x = period, y = estimate, fill = period)) +
geom_col(width = 0.55) +
geom_errorbar(aes(ymin = estimate - 1.96*se, ymax = estimate + 1.96*se), width = 0.15) +
scale_fill_manual(values = c("#FFC914", "#2E86AB")) +
labs(title = "The Income-Health Payoff Has Shrunk Over Time",
subtitle = "Fixed-effects coefficient on log(GDP per capita), within-country variation only",
x = NULL, y = "Effect of a 1-log-point income rise\non life expectancy (years)",
caption = "Implication: economic growth alone delivers smaller health gains today than a generation ago -\ndirect investment in health systems matters more at the margin.") +
un_theme + theme(legend.position = "none")
print(p4)
# ggsave("output/04_gradient_over_time.png", p4, width = 7.5, height = 6, dpi = 200)
## ---- Automatic interpretation: Chart 4 ----
early_est <- period_df$estimate[1]
recent_est <- period_df$estimate[2]
pct_shrink <- round((early_est - recent_est) / early_est * 100, 1)
cat("\n--- INTERPRETATION: Chart 4 (Gradient Over Time) ---\n")
##
## --- INTERPRETATION: Chart 4 (Gradient Over Time) ---
cat(sprintf("Early era (1952-1977) coefficient: %.2f\n", early_est))
## Early era (1952-1977) coefficient: 12.33
cat(sprintf("Recent era (1982-2007) coefficient: %.2f\n", recent_est))
## Recent era (1982-2007) coefficient: -0.44
cat(sprintf(
"=> The income-health payoff has %s by about %.1f%% between the two eras.\n %s\n",
ifelse(pct_shrink > 0, "shrunk", "grown"), abs(pct_shrink),
ifelse(pct_shrink > 0,
"Policy implication: as economies develop, growth alone buys less health -\n direct investment in health systems becomes relatively more important.",
"Policy implication: the income-health link has strengthened over time.")
))
## => The income-health payoff has shrunk by about 103.6% between the two eras.
## Policy implication: as economies develop, growth alone buys less health -
## direct investment in health systems becomes relatively more important.