The full analytical workflow for this project, including data preparation, exploratory analysis, panel model estimation, diagnostics, and visualization, is documented in a companion R Markdown appendix. The appendix contains all R code used to generate the results presented here and was rendered to ensure reproducibility. To maintain clarity in the presentation, detailed code and supplemental figures are provided separately and are available at the following RPubs link.
library(tidyverse)
library(readr)
library(jsonlite)
library(plotly)
library(plm)
library(lme4)
library(nlme)
library(mgcv)
library(splines)
library(broom)
library(modelsummary)
library(performance)
df <- read_csv("life-expectancy.csv")
metadata <- fromJSON("life-expectancy.metadata.json")
df_panel <- df %>%
mutate(Code = toupper(Code)) %>%
filter(!is.na(Code)) %>%
select(Entity, Code, Year, `Period life expectancy at birth`) %>%
rename(life_expectancy = `Period life expectancy at birth`) %>%
group_by(Code, Entity) %>%
arrange(Year, .by_group = TRUE) %>%
tidyr::complete(Year = seq(min(Year), max(Year), by = 1)) %>%
tidyr::fill(Entity, .direction = "downup") %>%
tidyr::fill(life_expectancy, .direction = "down") %>%
ungroup() %>%
mutate(time_since_1950 = Year - 1950,
decade = factor(floor(Year / 10) * 10),
Code = factor(Code),
Entity = factor(Entity))
global_trend <- df_panel %>% group_by(Year) %>% summarise(mean_life = mean(life_expectancy, na.rm = TRUE))
country_summary <- df_panel %>% group_by(Code) %>% summarise(mean_life = mean(life_expectancy), sd_life = sd(life_expectancy), first_year = min(Year), last_year = max(Year))
ggplot(global_trend, aes(x = Year, y = mean_life)) +
geom_line(color = "#9C3E59", linewidth = 1.01) +
labs(title = "Global Life Expectancy Trend", x = "Year", y = "Average Life Expectancy") +
theme_minimal(base_family = "Cormorant Garamond") +
theme(panel.background = element_rect(fill = NA, color = NA),
plot.background = element_rect(fill = NA, color = NA),
panel.grid.major = element_line(color = "#422E12"),
panel.grid.minor = element_blank(),
text = element_text(color = "#E0D7C5"),
plot.title = element_text(color = "#E0D7C5", size = 20, face = "bold"),
axis.text = element_text(color = "#E0D7C5", size = 14, face = "bold"),
axis.title = element_text(color = "#E0D7C5",size = 14),)
selected_codes <- df_panel %>% group_by(Code) %>% summarise(max_year = max(Year)) %>% arrange(desc(max_year)) %>% slice_head(n = 12) %>% pull(Code)
plot_df <- df_panel %>% filter(Code %in% selected_codes)
plot_ly(plot_df, x = ~Year, y = ~life_expectancy, color = ~Code, type = "scatter", mode = "lines",
colors = c("#6C5280", "#9C3E59", "#BD6262", "#E8DB74", "#74852A", "#DB8F48")) %>%
layout(title = list(text = "<b>Life Expectancy Paths for Selected Countries</b>", font = list(size = 22, color = "#e5e7eb")),
paper_bgcolor = "rgba(0,0,0,0)",
plot_bgcolor = "rgba(0,0,0,0)",
font = list(color = "#e5e7eb", family = "Cormorant Garamond"),
xaxis = list(title = "Year", gridcolor = "#422E12"),
yaxis = list(title = "Life Expectancy", gridcolor = "#422E12"),
hoverlabel = list(bgcolor = "#E0D7C5", bordercolor = "#8A623A",
font = list(family = "Cormorant Garamond, serif", size = 16, color = "#1f2933")),
margin = list(l = 70, r = 50, t = 60, b = 70))
df_panel %>%
mutate(period = case_when(Year < 1975 ~ "1950-1974", Year < 2000 ~ "1975-1999", TRUE ~ "2000-2023")) %>%
ggplot(aes(x = period, y = life_expectancy, fill = period)) +
geom_violin(color = "#0b1e2d", alpha = 0.8) +
geom_boxplot(width = 0.2, outlier.shape = NA, color = "#0b1e2d", alpha = 0.6) +
scale_fill_manual(values = c("#7f5539", "#b08968", "#e0afa0")) +
labs(title = "Distribution Shifts Over Time") +
theme_minimal(base_family = "Cormorant Garamond") +
theme(panel.background = element_rect(fill = NA, color = NA),
plot.background = element_rect(fill = NA, color = NA),
panel.grid.major = element_line(color = "#422E12"),
panel.grid.minor = element_blank(),
text = element_text(color = "#1f2933"),
plot.title = element_text(color = "#1f2933",size = 20, face = "bold"),
axis.text = element_text(color = "#1f2933",size = 14, face = "bold"),
axis.title = element_text(color = "#1f2933",size = 14, face = "bold"),
legend.position = "none")
country_volatility <- df_panel %>% group_by(Code) %>% summarise(volatility = sd(life_expectancy), trend = coef(lm(life_expectancy ~ Year))[2])
country_volatility %>%
ggplot(aes(x = trend, y = volatility, color = trend)) +
geom_point(alpha = 0.78) +
geom_smooth(method = "lm", color = "#8A623A", se = FALSE, linewidth = 1) +
scale_color_gradientn(colors = c("#4B3F52", "#9C3E59", "#BD6262", "#E8DB74", "#74852A")) +
labs(title = "Trend vs Variability Across Countries", x = "Slope of Trend", y = "SD of Life Expectancy") +
theme_minimal(base_family = "Cormorant Garamond") +
theme(panel.background = element_rect(fill = NA, color = NA),
plot.background = element_rect(fill = NA, color = NA),
panel.grid.major = element_line(color = "#422E12"),
panel.grid.minor = element_blank(),
text = element_text(color = "#1f2933"),
plot.title = element_text(color = "#1f2933",size = 20, face = "bold"),
axis.text = element_text(color = "#1f2933",size = 14, face = "bold"),
axis.title = element_text(color = "#1f2933",size = 14, face = "bold"),
legend.position = "none")
panel_df <- pdata.frame(df_panel, index = c("Code", "Year"))
m1_pool <- plm(life_expectancy ~ time_since_1950, data = panel_df, model = "pooling")
m2_pool_quad <- plm(life_expectancy ~ poly(time_since_1950, 2, raw = TRUE), data = panel_df, model = "pooling")
m3_fe <- plm(life_expectancy ~ time_since_1950, data = panel_df, model = "within")
m4_fe_poly <- plm(life_expectancy ~ poly(time_since_1950, 2, raw = TRUE), data = panel_df, model = "within")
m5_re <- plm(life_expectancy ~ time_since_1950, data = panel_df, model = "random")
m6_mixed_ri <- lmer(life_expectancy ~ time_since_1950 + (1 | Code), data = df_panel, REML = TRUE)
m7_mixed_rs <- lmer(life_expectancy ~ time_since_1950 + (time_since_1950 | Code), data = df_panel, REML = TRUE)
m8_spline_fe <- lm(life_expectancy ~ ns(time_since_1950, df = 6) + Code, data = df_panel)
m9_gam <- gam(life_expectancy ~ s(time_since_1950, k = 12) + s(Code, bs = "re"), data = df_panel, method = "REML")
model_list <- list(m1_pool = m1_pool, m2_pool_quad = m2_pool_quad, m3_fe = m3_fe, m4_fe_poly = m4_fe_poly, m5_re = m5_re, m6_mixed_ri = m6_mixed_ri, m7_mixed_rs = m7_mixed_rs, m8_spline_fe = m8_spline_fe, m9_gam = m9_gam)
pool_diag <- tibble(
fitted = fitted(m1_pool),
resid = residuals(m1_pool)
)
pool_diag %>%
ggplot(aes(x = fitted, y = resid, fill = resid)) +
geom_point(
shape = 21,
color = "#1f2933",
stroke = 0.1,
alpha = 0.55
) +
scale_fill_gradientn(
colors = c("#4B3F52", "#9C3E59", "#BD6262", "#E8DB74", "#74852A")
) +
geom_hline(yintercept = 0, color = "#7f5539", linewidth = 1.1) +
labs(title = "Pooled OLS Residuals vs Fitted") +
theme_minimal(base_family = "Cormorant Garamond") +
theme(
panel.background = element_rect(fill = NA, color = NA),
plot.background = element_rect(fill = NA, color = NA),
panel.grid.major = element_line(color = "#422E12"),
panel.grid.minor = element_blank(),
text = element_text(color = "#1f2933"),
plot.title = element_text(color = "#1f2933",size = 20, face = "bold"),
axis.text = element_text(color = "#1f2933",size = 14, face = "bold"),
axis.title = element_text(color = "#1f2933",size = 14, face = "bold"),
legend.position = "none"
)
fe_diag <- tibble(
fitted = fitted(m3_fe),
resid = residuals(m3_fe)
)
fe_diag %>%
ggplot(aes(x = fitted, y = resid, fill = fitted)) +
geom_point(
shape = 21,
color = "#1f2933", # border color
stroke = 0.1,
alpha = 0.55
) +
scale_fill_gradientn(
colors = c("#52396E", "#9C3E59", "#BD6262", "#E8DB74", "#74852A")
) +
geom_hline(yintercept = 0, color = "#e0afa0", linewidth = 1.1) +
labs(title = "Fixed Effects Residuals vs Fitted") +
theme_minimal(base_family = "Cormorant Garamond") +
theme(
panel.background = element_rect(fill = NA, color = NA),
plot.background = element_rect(fill = NA, color = NA),
panel.grid.major = element_line(color = "#422E12"),
panel.grid.minor = element_blank(),
text = element_text(color = "#1f2933"),
plot.title = element_text(color = "#1f2933",size = 20, face = "bold"),
axis.text = element_text(color = "#1f2933",size = 14, face = "bold"),
axis.title = element_text(color = "#1f2933",size = 14, face = "bold"),
legend.position = "none"
)
mixed_diag <- tibble(resid = resid(m6_mixed_ri))
ggplot(mixed_diag, aes(sample = resid)) +
stat_qq(
shape = 21,
fill = "#BD6262", # interior color
color = "#1f2933", # border color
stroke = 0.1,
alpha = 0.65
) +
stat_qq_line(color = "#74852A", linewidth = 1.05) +
labs(title = "Random Intercept QQ Plot") +
theme_minimal(base_family = "Cormorant Garamond") +
theme(
panel.background = element_rect(fill = NA, color = NA),
plot.background = element_rect(fill = NA, color = NA),
panel.grid.major = element_line(color = "#422E12"),
panel.grid.minor = element_blank(),
text = element_text(color = "#1f2933"),
plot.title = element_text(color = "#1f2933",size = 20, face = "bold"),
axis.text = element_text(color = "#1f2933",size = 14, face = "bold"),
axis.title = element_text(color = "#1f2933",size = 20, face = "bold")
)
re_effects <- ranef(m7_mixed_rs)
re_df <- tibble(
Code = rownames(re_effects[[1]]),
intercept = re_effects[[1]][, 1],
slope = re_effects[[1]][, 2]
)
re_df %>%
ggplot(aes(x = intercept, y = slope, fill = slope, size = abs(intercept))) +
geom_point(
shape = 21,
color = "#1f2933", # border
stroke = 0.6,
alpha = 0.78
) +
scale_fill_gradientn(
colors = c("#4B3F52", "#9C3E59", "#BD6262", "#E8DB74", "#74852A")
) +
geom_hline(yintercept = 0, color = "#8A623A", linewidth = 0.9, linetype = "dashed") +
geom_vline(xintercept = 0, color = "#8A623A", linewidth = 0.9, linetype = "dashed") +
labs(title = "Random Intercepts and Slopes") +
theme_minimal(base_family = "Cormorant Garamond") +
theme(
panel.background = element_rect(fill = NA, color = NA),
plot.background = element_rect(fill = NA, color = NA),
panel.grid.major = element_line(color = "#422E12"),
panel.grid.minor = element_blank(),
text = element_text(color = "#1f2933"),
plot.title = element_text(color = "#1f2933",size = 20, face = "bold"),
axis.text = element_text(color = "#1f2933",size = 14, face = "bold"),
axis.title = element_text(color = "#1f2933",size = 14, face = "bold"),
legend.position = "none"
)
comparison_df <- tibble(model = names(model_list),
AIC = map_dbl(model_list, AIC),
BIC = map_dbl(model_list, BIC))
comparison_df
## # A tibble: 9 × 3
## model AIC BIC
## <chr> <dbl> <dbl>
## 1 m1_pool 175544. 175568.
## 2 m2_pool_quad 170015. 170048.
## 3 m3_fe 150684. 150700.
## 4 m4_fe_poly 135992. 136016.
## 5 m5_re 150973. 150997.
## 6 m6_mixed_ri 152167. 152199.
## 7 m7_mixed_rs 127310. 127358.
## 8 m8_spline_fe 131403. 133377.
## 9 m9_gam 129653. 131661.
modelsummary(model_list[c("m1_pool", "m2_pool_quad", "m3_fe", "m4_fe_poly", "m5_re", "m8_spline_fe")], statistic = "{statistic}")
| m1_pool | m2_pool_quad | m3_fe | m4_fe_poly | m5_re | m8_spline_fe | |
|---|---|---|---|---|---|---|
| (Intercept) | 55.839 | 52.245 | 53.203 | 37.728 | ||
| 683.214 | 611.860 | 94.976 | 47.221 | |||
| time_since_1950 | 0.182 | 0.257 | 0.257 | |||
| 123.233 | 221.685 | 221.455 | ||||
| poly(time_since_1950, 2, raw = TRUE)1 | 0.270 | 0.325 | ||||
| 157.255 | 337.592 | |||||
| poly(time_since_1950, 2, raw = TRUE)2 | 0.001 | 0.001 | ||||
| 79.151 | 143.106 | |||||
| ns(time_since_1950, df = 6)1 | 18.929 | |||||
| 31.263 | ||||||
| ns(time_since_1950, df = 6)2 | 33.585 | |||||
| 52.844 | ||||||
| ns(time_since_1950, df = 6)3 | 33.124 | |||||
| 52.874 | ||||||
| ns(time_since_1950, df = 6)4 | 49.603 | |||||
| 111.291 | ||||||
| ns(time_since_1950, df = 6)5 | 13.959 | |||||
| 10.922 | ||||||
| ns(time_since_1950, df = 6)6 | 61.178 | |||||
| 202.882 | ||||||
| CodeAFG | -24.480 | |||||
| -33.807 | ||||||
| CodeAGO | -24.418 | |||||
| -34.771 | ||||||
| CodeAIA | -0.427 | |||||
| -0.590 | ||||||
| CodeALB | -1.373 | |||||
| -1.896 | ||||||
| CodeAND | 7.845 | |||||
| 10.833 | ||||||
| CodeARE | -2.350 | |||||
| -3.246 | ||||||
| CodeARG | 0.367 | |||||
| 0.585 | ||||||
| CodeARM | -3.838 | |||||
| -5.301 | ||||||
| CodeASM | -1.202 | |||||
| -1.660 | ||||||
| CodeATG | 1.189 | |||||
| 1.642 | ||||||
| CodeAUS | 10.028 | |||||
| 15.797 | ||||||
| CodeAUT | 2.886 | |||||
| 4.621 | ||||||
| CodeAZE | -8.830 | |||||
| -12.194 | ||||||
| CodeBDI | -22.177 | |||||
| -30.626 | ||||||
| CodeBEL | 7.368 | |||||
| 12.090 | ||||||
| CodeBEN | -20.383 | |||||
| -28.150 | ||||||
| CodeBES | 3.191 | |||||
| 4.407 | ||||||
| CodeBFA | -23.133 | |||||
| -31.947 | ||||||
| CodeBGD | -17.119 | |||||
| -27.245 | ||||||
| CodeBGR | 0.575 | |||||
| 0.887 | ||||||
| CodeBHR | -2.424 | |||||
| -3.347 | ||||||
| CodeBHS | -1.237 | |||||
| -1.709 | ||||||
| CodeBIH | -2.065 | |||||
| -2.852 | ||||||
| CodeBLM | 2.380 | |||||
| 3.287 | ||||||
| CodeBLR | -4.178 | |||||
| -6.452 | ||||||
| CodeBLZ | -4.561 | |||||
| -6.299 | ||||||
| CodeBMU | 3.111 | |||||
| 4.296 | ||||||
| CodeBOL | -15.265 | |||||
| -23.571 | ||||||
| CodeBRA | -9.043 | |||||
| -13.963 | ||||||
| CodeBRB | 0.628 | |||||
| 0.868 | ||||||
| CodeBRN | -1.941 | |||||
| -2.681 | ||||||
| CodeBTN | -16.952 | |||||
| -23.411 | ||||||
| CodeBWA | -12.817 | |||||
| -17.701 | ||||||
| CodeCAF | -25.062 | |||||
| -34.611 | ||||||
| CodeCAN | 8.593 | |||||
| 14.192 | ||||||
| CodeCHE | 8.820 | |||||
| 14.036 | ||||||
| CodeCHL | -4.533 | |||||
| -7.000 | ||||||
| CodeCHN | -8.028 | |||||
| -11.725 | ||||||
| CodeCIV | -21.393 | |||||
| -29.544 | ||||||
| CodeCMR | -20.073 | |||||
| -29.250 | ||||||
| CodeCOD | -21.779 | |||||
| -30.077 | ||||||
| CodeCOG | -14.948 | |||||
| -20.643 | ||||||
| CodeCOK | -5.592 | |||||
| -7.723 | ||||||
| CodeCOL | -7.450 | |||||
| -11.504 | ||||||
| CodeCOM | -17.370 | |||||
| -23.988 | ||||||
| CodeCPV | -8.045 | |||||
| -11.110 | ||||||
| CodeCRI | -2.465 | |||||
| -3.928 | ||||||
| CodeCUB | -1.855 | |||||
| -2.869 | ||||||
| CodeCUW | 0.686 | |||||
| 0.947 | ||||||
| CodeCYM | -0.715 | |||||
| -0.988 | ||||||
| CodeCYP | 1.756 | |||||
| 2.731 | ||||||
| CodeCZE | 1.926 | |||||
| 3.068 | ||||||
| CodeDEU | 4.619 | |||||
| 7.359 | ||||||
| CodeDJI | -16.988 | |||||
| -23.460 | ||||||
| CodeDMA | -2.908 | |||||
| -4.016 | ||||||
| CodeDNK | 9.913 | |||||
| 16.811 | ||||||
| CodeDOM | -9.376 | |||||
| -13.695 | ||||||
| CodeDZA | -12.901 | |||||
| -19.135 | ||||||
| CodeECU | -4.861 | |||||
| -6.714 | ||||||
| CodeEGY | -13.407 | |||||
| -19.716 | ||||||
| CodeERI | -21.150 | |||||
| -29.208 | ||||||
| CodeESH | -17.837 | |||||
| -24.633 | ||||||
| CodeESP | 1.895 | |||||
| 2.995 | ||||||
| CodeEST | 2.076 | |||||
| 3.220 | ||||||
| CodeETH | -22.631 | |||||
| -31.254 | ||||||
| CodeFIN | 7.118 | |||||
| 12.147 | ||||||
| CodeFJI | -7.160 | |||||
| -9.887 | ||||||
| CodeFLK | -8.070 | |||||
| -11.144 | ||||||
| CodeFRA | 7.510 | |||||
| 12.512 | ||||||
| CodeFRO | 5.293 | |||||
| 7.310 | ||||||
| CodeFSM | -8.930 | |||||
| -12.333 | ||||||
| CodeGAB | -12.891 | |||||
| -17.802 | ||||||
| CodeGBR | 9.596 | |||||
| 16.458 | ||||||
| CodeGEO | -2.602 | |||||
| -3.594 | ||||||
| CodeGGY | 6.849 | |||||
| 9.459 | ||||||
| CodeGHA | -17.630 | |||||
| -26.256 | ||||||
| CodeGIB | 5.689 | |||||
| 7.857 | ||||||
| CodeGIN | -23.547 | |||||
| -32.519 | ||||||
| CodeGLP | 0.826 | |||||
| 1.140 | ||||||
| CodeGMB | -21.437 | |||||
| -29.604 | ||||||
| CodeGNB | -23.569 | |||||
| -32.549 | ||||||
| CodeGNQ | -20.602 | |||||
| -28.451 | ||||||
| CodeGRC | 3.023 | |||||
| 4.806 | ||||||
| CodeGRD | -1.351 | |||||
| -1.865 | ||||||
| CodeGRL | -7.519 | |||||
| -10.385 | ||||||
| CodeGTM | -15.580 | |||||
| -24.058 | ||||||
| CodeGUF | -2.909 | |||||
| -4.017 | ||||||
| CodeGUM | 0.476 | |||||
| 0.657 | ||||||
| CodeGUY | -9.203 | |||||
| -13.964 | ||||||
| CodeHKG | 5.197 | |||||
| 7.177 | ||||||
| CodeHND | -11.478 | |||||
| -17.128 | ||||||
| CodeHRV | 1.279 | |||||
| 1.767 | ||||||
| CodeHTI | -17.787 | |||||
| -24.564 | ||||||
| CodeHUN | -0.207 | |||||
| -0.319 | ||||||
| CodeIDN | -12.999 | |||||
| -19.116 | ||||||
| CodeIMN | 1.299 | |||||
| 1.794 | ||||||
| CodeIND | -15.293 | |||||
| -24.205 | ||||||
| CodeIRL | 6.269 | |||||
| 9.665 | ||||||
| CodeIRN | -9.770 | |||||
| -13.493 | ||||||
| CodeIRQ | -9.554 | |||||
| -13.194 | ||||||
| CodeISL | 7.956 | |||||
| 13.080 | ||||||
| CodeISR | 5.452 | |||||
| 7.529 | ||||||
| CodeITA | 4.066 | |||||
| 6.499 | ||||||
| CodeJAM | -2.670 | |||||
| -4.225 | ||||||
| CodeJEY | 5.013 | |||||
| 6.923 | ||||||
| CodeJOR | -4.719 | |||||
| -6.516 | ||||||
| CodeJPN | 2.397 | |||||
| 3.857 | ||||||
| CodeKAZ | -9.086 | |||||
| -14.579 | ||||||
| CodeKEN | -16.875 | |||||
| -24.816 | ||||||
| CodeKGZ | -7.943 | |||||
| -10.970 | ||||||
| CodeKHM | -18.696 | |||||
| -26.238 | ||||||
| CodeKIR | -11.756 | |||||
| -16.236 | ||||||
| CodeKNA | -5.183 | |||||
| -7.158 | ||||||
| CodeKOR | -6.727 | |||||
| -10.258 | ||||||
| CodeKWT | -7.523 | |||||
| -11.472 | ||||||
| CodeLAO | -17.539 | |||||
| -24.222 | ||||||
| CodeLBN | -2.125 | |||||
| -2.935 | ||||||
| CodeLBR | -23.034 | |||||
| -31.810 | ||||||
| CodeLBY | -8.803 | |||||
| -12.157 | ||||||
| CodeLCA | -3.250 | |||||
| -4.488 | ||||||
| CodeLIE | 6.142 | |||||
| 8.482 | ||||||
| CodeLKA | -5.472 | |||||
| -8.437 | ||||||
| CodeLSO | -18.199 | |||||
| -25.133 | ||||||
| CodeLTU | 0.613 | |||||
| 0.946 | ||||||
| CodeLUX | 5.250 | |||||
| 8.094 | ||||||
| CodeLVA | 2.107 | |||||
| 3.272 | ||||||
| CodeMAC | 3.705 | |||||
| 5.117 | ||||||
| CodeMAF | 3.482 | |||||
| 4.809 | ||||||
| CodeMAR | -11.832 | |||||
| -16.341 | ||||||
| CodeMCO | 7.902 | |||||
| 10.912 | ||||||
| CodeMDA | -5.101 | |||||
| -7.045 | ||||||
| CodeMDG | -17.489 | |||||
| -24.153 | ||||||
| CodeMDV | -11.214 | |||||
| -15.486 | ||||||
| CodeMEX | -8.103 | |||||
| -12.636 | ||||||
| CodeMHL | -9.432 | |||||
| -13.025 | ||||||
| CodeMKD | -1.812 | |||||
| -2.502 | ||||||
| CodeMLI | -26.386 | |||||
| -36.439 | ||||||
| CodeMLT | 4.442 | |||||
| 6.135 | ||||||
| CodeMMR | -16.273 | |||||
| -23.984 | ||||||
| CodeMNE | 1.182 | |||||
| 1.633 | ||||||
| CodeMNG | -13.094 | |||||
| -18.083 | ||||||
| CodeMNP | -0.449 | |||||
| -0.620 | ||||||
| CodeMOZ | -23.200 | |||||
| -32.039 | ||||||
| CodeMRT | -14.674 | |||||
| -20.265 | ||||||
| CodeMSR | -2.202 | |||||
| -3.040 | ||||||
| CodeMTQ | 2.188 | |||||
| 3.022 | ||||||
| CodeMUS | -6.632 | |||||
| -9.816 | ||||||
| CodeMWI | -24.591 | |||||
| -33.961 | ||||||
| CodeMYS | -2.639 | |||||
| -3.644 | ||||||
| CodeMYT | -2.058 | |||||
| -2.842 | ||||||
| CodeNAM | -15.009 | |||||
| -20.728 | ||||||
| CodeNCL | -4.156 | |||||
| -5.740 | ||||||
| CodeNER | -24.421 | |||||
| -36.371 | ||||||
| CodeNGA | -25.291 | |||||
| -34.927 | ||||||
| CodeNIC | -13.415 | |||||
| -20.019 | ||||||
| CodeNIU | -5.046 | |||||
| -6.968 | ||||||
| CodeNLD | 8.723 | |||||
| 14.217 | ||||||
| CodeNOR | 12.111 | |||||
| 19.800 | ||||||
| CodeNPL | -16.920 | |||||
| -23.367 | ||||||
| CodeNRU | -11.054 | |||||
| -15.266 | ||||||
| CodeNZL | 5.532 | |||||
| 7.690 | ||||||
| CodeOMN | -10.008 | |||||
| -13.821 | ||||||
| CodeOWID_KOS | -7.634 | |||||
| -10.543 | ||||||
| CodeOWID_USS | -11.478 | |||||
| -8.920 | ||||||
| CodeOWID_WRL | -4.710 | |||||
| -8.002 | ||||||
| CodePAK | -16.595 | |||||
| -24.715 | ||||||
| CodePAN | -3.246 | |||||
| -4.741 | ||||||
| CodePER | -8.505 | |||||
| -12.111 | ||||||
| CodePHL | -7.215 | |||||
| -10.330 | ||||||
| CodePLW | -8.297 | |||||
| -11.458 | ||||||
| CodePNG | -14.267 | |||||
| -19.961 | ||||||
| CodePOL | 0.613 | |||||
| 0.893 | ||||||
| CodePRI | 3.085 | |||||
| 4.260 | ||||||
| CodePRK | -9.383 | |||||
| -14.309 | ||||||
| CodePRT | 2.063 | |||||
| 2.937 | ||||||
| CodePRY | -7.595 | |||||
| -11.727 | ||||||
| CodePSE | -7.146 | |||||
| -9.868 | ||||||
| CodePYF | 0.502 | |||||
| 0.693 | ||||||
| CodeQAT | -0.016 | |||||
| -0.022 | ||||||
| CodeREU | 0.039 | |||||
| 0.054 | ||||||
| CodeROU | -2.220 | |||||
| -3.228 | ||||||
| CodeRUS | -6.045 | |||||
| -9.388 | ||||||
| CodeRWA | -19.198 | |||||
| -26.512 | ||||||
| CodeSAU | -7.082 | |||||
| -9.781 | ||||||
| CodeSDN | -14.522 | |||||
| -20.054 | ||||||
| CodeSEN | -19.607 | |||||
| -28.833 | ||||||
| CodeSGP | 2.989 | |||||
| 4.127 | ||||||
| CodeSHN | 1.649 | |||||
| 2.278 | ||||||
| CodeSLB | -9.179 | |||||
| -12.676 | ||||||
| CodeSLE | -26.656 | |||||
| -38.842 | ||||||
| CodeSLV | -13.004 | |||||
| -19.405 | ||||||
| CodeSMR | 7.851 | |||||
| 10.842 | ||||||
| CodeSOM | -22.793 | |||||
| -31.478 | ||||||
| CodeSPM | 1.283 | |||||
| 1.772 | ||||||
| CodeSRB | 0.350 | |||||
| 0.484 | ||||||
| CodeSSD | -32.305 | |||||
| -44.613 | ||||||
| CodeSTP | -12.825 | |||||
| -17.712 | ||||||
| CodeSUR | -7.724 | |||||
| -10.668 | ||||||
| CodeSVK | 1.248 | |||||
| 1.858 | ||||||
| CodeSVN | 2.967 | |||||
| 4.097 | ||||||
| CodeSWE | 12.125 | |||||
| 20.715 | ||||||
| CodeSWZ | -17.833 | |||||
| -24.628 | ||||||
| CodeSXM | 1.456 | |||||
| 2.011 | ||||||
| CodeSYC | -2.540 | |||||
| -3.508 | ||||||
| CodeSYR | -6.222 | |||||
| -8.592 | ||||||
| CodeTCA | -2.498 | |||||
| -3.449 | ||||||
| CodeTCD | -25.139 | |||||
| -34.717 | ||||||
| CodeTGO | -18.274 | |||||
| -25.237 | ||||||
| CodeTHA | -6.683 | |||||
| -9.594 | ||||||
| CodeTJK | -9.586 | |||||
| -13.239 | ||||||
| CodeTKL | -2.896 | |||||
| -3.999 | ||||||
| CodeTKM | -9.273 | |||||
| -12.806 | ||||||
| CodeTLS | -25.462 | |||||
| -35.163 | ||||||
| CodeTON | -4.531 | |||||
| -6.258 | ||||||
| CodeTTO | -3.660 | |||||
| -5.451 | ||||||
| CodeTUN | -10.579 | |||||
| -15.691 | ||||||
| CodeTUR | -7.479 | |||||
| -10.736 | ||||||
| CodeTUV | -11.715 | |||||
| -16.178 | ||||||
| CodeTWN | 1.737 | |||||
| 2.399 | ||||||
| CodeTZA | -18.041 | |||||
| -24.915 | ||||||
| CodeUGA | -20.726 | |||||
| -30.480 | ||||||
| CodeUKR | -2.902 | |||||
| -4.482 | ||||||
| CodeURY | 3.538 | |||||
| 5.463 | ||||||
| CodeUSA | 6.788 | |||||
| 10.755 | ||||||
| CodeUZB | -6.485 | |||||
| -8.956 | ||||||
| CodeVAT | 5.421 | |||||
| 7.487 | ||||||
| CodeVCT | -4.028 | |||||
| -5.563 | ||||||
| CodeVEN | -7.592 | |||||
| -11.723 | ||||||
| CodeVGB | -0.199 | |||||
| -0.275 | ||||||
| CodeVIR | -1.276 | |||||
| -1.762 | ||||||
| CodeVNM | -4.801 | |||||
| -6.630 | ||||||
| CodeVUT | -8.572 | |||||
| -11.838 | ||||||
| CodeWLF | -6.533 | |||||
| -9.023 | ||||||
| CodeWSM | -5.549 | |||||
| -7.663 | ||||||
| CodeYEM | -18.081 | |||||
| -24.970 | ||||||
| CodeZAF | -12.397 | |||||
| -17.121 | ||||||
| CodeZMB | -17.243 | |||||
| -23.813 | ||||||
| CodeZWE | -15.370 | |||||
| -21.227 | ||||||
| Num.Obs. | 22601 | 22601 | 22601 | 22601 | 22601 | 22601 |
| R2 | 0.402 | 0.532 | 0.687 | 0.837 | 0.687 | 0.917 |
| R2 Adj. | 0.402 | 0.532 | 0.684 | 0.835 | 0.687 | 0.916 |
| AIC | 175543.8 | 170015.5 | 150684.4 | 135991.5 | 150973.3 | 131402.8 |
| BIC | 175567.9 | 170047.6 | 150700.4 | 136015.6 | 150997.4 | 133377.1 |
| Log.Lik. | -65455.395 | |||||
| F | 1011.897 | |||||
| RMSE | 11.76 | 10.40 | 6.78 | 4.90 | 6.83 | 4.38 |
loess_fit <- loess(life_expectancy ~ time_since_1950, data = df_panel, span = 0.3)
loess_df <- tibble(Year = df_panel$Year, fitted = predict(loess_fit)) %>% group_by(Year) %>% summarise(fitted = mean(fitted, na.rm = TRUE))
ggplot() +
geom_line(data = global_trend, aes(x = Year, y = mean_life), color = "#BD6262", linewidth = 1.1) +
geom_line(data = loess_df, aes(x = Year, y = fitted), color = "#74852A", linewidth = 1.2, linetype = "dashed") +
labs(title = "Loess Smoother vs Global Trend", y = "Life Expectancy") +
theme_minimal(base_family = "Cormorant Garamond") +
theme(panel.background = element_rect(fill = NA, color = NA),
plot.background = element_rect(fill = NA, color = NA),
panel.grid.major = element_line(color = "#422E12"),
panel.grid.minor = element_blank(),
text = element_text(color = "#1f2933"),
plot.title = element_text(color = "#1f2933",size = 20, face = "bold"),
axis.text = element_text(color = "#1f2933",size = 14, face = "bold"),
axis.title = element_text(color = "#1f2933",size = 14, face = "bold"))
diagnostic_objects <- list(pool = pool_diag, fe = fe_diag, mixed = mixed_diag, random_slopes = re_df)
save_objects <- list(models = model_list, comparison = comparison_df, diagnostics = diagnostic_objects, volatility = country_volatility, loess_fit = loess_fit)
save_objects
## $models
## $models$m1_pool
##
## Model Formula: life_expectancy ~ time_since_1950
##
## Coefficients:
## (Intercept) time_since_1950
## 55.83876 0.18164
##
##
## $models$m2_pool_quad
##
## Model Formula: life_expectancy ~ poly(time_since_1950, 2, raw = TRUE)
##
## Coefficients:
## (Intercept) poly(time_since_1950, 2, raw = TRUE)1
## 5.2245e+01 2.7011e-01
## poly(time_since_1950, 2, raw = TRUE)2
## 7.0559e-04
##
##
## $models$m3_fe
##
## Model Formula: life_expectancy ~ time_since_1950
##
## Coefficients:
## time_since_1950
## 0.2575
##
##
## $models$m4_fe_poly
##
## Model Formula: life_expectancy ~ poly(time_since_1950, 2, raw = TRUE)
##
## Coefficients:
## poly(time_since_1950, 2, raw = TRUE)1 poly(time_since_1950, 2, raw = TRUE)2
## 0.32505738 0.00069655
##
##
## $models$m5_re
##
## Model Formula: life_expectancy ~ time_since_1950
##
## Coefficients:
## (Intercept) time_since_1950
## 53.20279 0.25723
##
##
## $models$m6_mixed_ri
## Linear mixed model fit by REML ['lmerMod']
## Formula: life_expectancy ~ time_since_1950 + (1 | Code)
## Data: df_panel
## REML criterion at convergence: 152158.9
## Random effects:
## Groups Name Std.Dev.
## Code (Intercept) 9.436
## Residual 6.820
## Number of obs: 22601, groups: Code, 239
## Fixed Effects:
## (Intercept) time_since_1950
## 53.2008 0.2573
##
## $models$m7_mixed_rs
## Linear mixed model fit by REML ['lmerMod']
## Formula: life_expectancy ~ time_since_1950 + (time_since_1950 | Code)
## Data: df_panel
## REML criterion at convergence: 127297.8
## Random effects:
## Groups Name Std.Dev. Corr
## Code (Intercept) 10.4931
## time_since_1950 0.1232 -0.60
## Residual 3.8317
## Number of obs: 22601, groups: Code, 239
## Fixed Effects:
## (Intercept) time_since_1950
## 50.5730 0.3507
##
## $models$m8_spline_fe
##
## Call:
## lm(formula = life_expectancy ~ ns(time_since_1950, df = 6) +
## Code, data = df_panel)
##
## Coefficients:
## (Intercept) ns(time_since_1950, df = 6)1
## 37.72776 18.92888
## ns(time_since_1950, df = 6)2 ns(time_since_1950, df = 6)3
## 33.58481 33.12436
## ns(time_since_1950, df = 6)4 ns(time_since_1950, df = 6)5
## 49.60280 13.95899
## ns(time_since_1950, df = 6)6 CodeAFG
## 61.17829 -24.48011
## CodeAGO CodeAIA
## -24.41842 -0.42704
## CodeALB CodeAND
## -1.37321 7.84458
## CodeARE CodeARG
## -2.35038 0.36699
## CodeARM CodeASM
## -3.83829 -1.20234
## CodeATG CodeAUS
## 1.18903 10.02757
## CodeAUT CodeAZE
## 2.88577 -8.82985
## CodeBDI CodeBEL
## -22.17651 7.36849
## CodeBEN CodeBES
## -20.38343 3.19140
## CodeBFA CodeBGD
## -23.13315 -17.11926
## CodeBGR CodeBHR
## 0.57472 -2.42388
## CodeBHS CodeBIH
## -1.23732 -2.06521
## CodeBLM CodeBLR
## 2.37989 -4.17820
## CodeBLZ CodeBMU
## -4.56138 3.11110
## CodeBOL CodeBRA
## -15.26458 -9.04273
## CodeBRB CodeBRN
## 0.62817 -1.94099
## CodeBTN CodeBWA
## -16.95186 -12.81744
## CodeCAF CodeCAN
## -25.06203 8.59260
## CodeCHE CodeCHL
## 8.81965 -4.53349
## CodeCHN CodeCIV
## -8.02783 -21.39307
## CodeCMR CodeCOD
## -20.07316 -21.77882
## CodeCOG CodeCOK
## -14.94782 -5.59224
## CodeCOL CodeCOM
## -7.45009 -17.36966
## CodeCPV CodeCRI
## -8.04495 -2.46549
## CodeCUB CodeCUW
## -1.85532 0.68551
## CodeCYM CodeCYP
## -0.71535 1.75583
## CodeCZE CodeDEU
## 1.92601 4.61912
## CodeDJI CodeDMA
## -16.98751 -2.90768
## CodeDNK CodeDOM
## 9.91341 -9.37630
## CodeDZA CodeECU
## -12.90080 -4.86146
## CodeEGY CodeERI
## -13.40704 -21.14960
## CodeESH CodeESP
## -17.83680 1.89480
## CodeEST CodeETH
## 2.07611 -22.63116
## CodeFIN CodeFJI
## 7.11812 -7.15953
## CodeFLK CodeFRA
## -8.06978 7.51002
## CodeFRO CodeFSM
## 5.29331 -8.93023
## CodeGAB CodeGBR
## -12.89077 9.59561
## CodeGEO CodeGGY
## -2.60212 6.84936
## CodeGHA CodeGIB
## -17.62993 5.68927
## CodeGIN CodeGLP
## -23.54693 0.82582
## CodeGMB CodeGNB
## -21.43650 -23.56876
## CodeGNQ CodeGRC
## -20.60156 3.02304
## CodeGRD CodeGRL
## -1.35072 -7.51946
## CodeGTM CodeGUF
## -15.57997 -2.90893
## CodeGUM CodeGUY
## 0.47593 -9.20335
## CodeHKG CodeHND
## 5.19720 -11.47787
## CodeHRV CodeHTI
## 1.27925 -17.78701
## CodeHUN CodeIDN
## -0.20660 -12.99863
## CodeIMN CodeIND
## 1.29933 -15.29333
## CodeIRL CodeIRN
## 6.26864 -9.77017
## CodeIRQ CodeISL
## -9.55400 7.95586
## CodeISR CodeITA
## 5.45153 4.06638
## CodeJAM CodeJEY
## -2.66974 5.01317
## CodeJOR CodeJPN
## -4.71860 2.39673
## CodeKAZ CodeKEN
## -9.08610 -16.87525
## CodeKGZ CodeKHM
## -7.94336 -18.69644
## CodeKIR CodeKNA
## -11.75631 -5.18335
## CodeKOR CodeKWT
## -6.72690 -7.52284
## CodeLAO CodeLBN
## -17.53948 -2.12543
## CodeLBR CodeLBY
## -23.03379 -8.80296
## CodeLCA CodeLIE
## -3.24985 6.14189
## CodeLKA CodeLSO
## -5.47220 -18.19886
## CodeLTU CodeLUX
## 0.61294 5.24964
## CodeLVA CodeMAC
## 2.10674 3.70513
## CodeMAF CodeMAR
## 3.48190 -11.83247
## CodeMCO CodeMDA
## 7.90177 -5.10124
## CodeMDG CodeMDV
## -17.48945 -11.21371
## CodeMEX CodeMHL
## -8.10344 -9.43167
## CodeMKD CodeMLI
## -1.81172 -26.38567
## CodeMLT CodeMMR
## 4.44249 -16.27322
## CodeMNE CodeMNG
## 1.18241 -13.09423
## CodeMNP CodeMOZ
## -0.44931 -23.19993
## CodeMRT CodeMSR
## -14.67411 -2.20160
## CodeMTQ CodeMUS
## 2.18808 -6.63156
## CodeMWI CodeMYS
## -24.59148 -2.63857
## CodeMYT CodeNAM
## -2.05769 -15.00948
## CodeNCL CodeNER
## -4.15629 -24.42145
## CodeNGA CodeNIC
## -25.29056 -13.41492
## CodeNIU CodeNLD
## -5.04556 8.72278
## CodeNOR CodeNPL
## 12.11121 -16.92001
## CodeNRU CodeNZL
## -11.05427 5.53153
## CodeOMN CodeOWID_KOS
## -10.00797 -7.63398
## CodeOWID_USS CodeOWID_WRL
## -11.47812 -4.71046
## CodePAK CodePAN
## -16.59501 -3.24588
## CodePER CodePHL
## -8.50525 -7.21496
## CodePLW CodePNG
## -8.29705 -14.26747
## CodePOL CodePRI
## 0.61280 3.08462
## CodePRK CodePRT
## -9.38309 2.06265
## CodePRY CodePSE
## -7.59465 -7.14561
## CodePYF CodeQAT
## 0.50184 -0.01586
## CodeREU CodeROU
## 0.03908 -2.22018
## CodeRUS CodeRWA
## -6.04529 -19.19758
## CodeSAU CodeSDN
## -7.08223 -14.52153
## CodeSEN CodeSGP
## -19.60682 2.98851
## CodeSHN CodeSLB
## 1.64934 -9.17887
## CodeSLE CodeSLV
## -26.65562 -13.00377
## CodeSMR CodeSOM
## 7.85080 -22.79336
## CodeSPM CodeSRB
## 1.28276 0.35039
## CodeSSD CodeSTP
## -32.30483 -12.82503
## CodeSUR CodeSVK
## -7.72450 1.24771
## CodeSVN CodeSWE
## 2.96695 12.12508
## CodeSWZ CodeSXM
## -17.83292 1.45632
## CodeSYC CodeSYR
## -2.53994 -6.22175
## CodeTCA CodeTCD
## -2.49763 -25.13868
## CodeTGO CodeTHA
## -18.27406 -6.68330
## CodeTJK CodeTKL
## -9.58620 -2.89563
## CodeTKM CodeTLS
## -9.27277 -25.46177
## CodeTON CodeTTO
## -4.53134 -3.66041
## CodeTUN CodeTUR
## -10.57896 -7.47904
## CodeTUV CodeTWN
## -11.71469 1.73728
## CodeTZA CodeUGA
## -18.04077 -20.72622
## CodeUKR CodeURY
## -2.90234 3.53757
## CodeUSA CodeUZB
## 6.78786 -6.48506
## CodeVAT CodeVCT
## 5.42143 -4.02796
## CodeVEN CodeVGB
## -7.59180 -0.19944
## CodeVIR CodeVNM
## -1.27605 -4.80086
## CodeVUT CodeWLF
## -8.57190 -6.53335
## CodeWSM CodeYEM
## -5.54855 -18.08124
## CodeZAF CodeZMB
## -12.39725 -17.24296
## CodeZWE
## -15.37026
##
##
## $models$m9_gam
##
## Family: gaussian
## Link function: identity
##
## Formula:
## life_expectancy ~ s(time_since_1950, k = 12) + s(Code, bs = "re")
##
## Estimated degrees of freedom:
## 10.7 237.4 total = 249.12
##
## REML score: 65460.43
##
##
## $comparison
## # A tibble: 9 × 3
## model AIC BIC
## <chr> <dbl> <dbl>
## 1 m1_pool 175544. 175568.
## 2 m2_pool_quad 170015. 170048.
## 3 m3_fe 150684. 150700.
## 4 m4_fe_poly 135992. 136016.
## 5 m5_re 150973. 150997.
## 6 m6_mixed_ri 152167. 152199.
## 7 m7_mixed_rs 127310. 127358.
## 8 m8_spline_fe 131403. 133377.
## 9 m9_gam 129653. 131661.
##
## $diagnostics
## $diagnostics$pool
## # A tibble: 22,601 × 2
## fitted resid
## <pseries> <pseries>
## 1 55.83876 2.148437
## 2 56.02041 2.708293
## 3 56.20205 3.243449
## 4 56.38369 3.735706
## 5 56.56534 4.251262
## 6 56.74698 4.705218
## 7 56.92863 5.115874
## 8 57.11027 5.489930
## 9 57.29191 5.861686
## 10 57.47356 6.137042
## # ℹ 22,591 more rows
##
## $diagnostics$fe
## # A tibble: 22,601 × 2
## fitted resid
## <pseries> <pseries>
## 1 -9.398609 -2.76065320
## 2 -9.141113 -2.27664934
## 3 -8.883617 -1.81734548
## 4 -8.626121 -1.40094161
## 5 -8.368624 -0.96123775
## 6 -8.111128 -0.58313388
## 7 -7.853632 -0.24833002
## 8 -7.596136 0.04987384
## 9 -7.338640 0.34577771
## 10 -7.081144 0.54528157
## # ℹ 22,591 more rows
##
## $diagnostics$mixed
## # A tibble: 22,601 × 1
## resid
## <dbl>
## 1 -2.72
## 2 -2.23
## 3 -1.77
## 4 -1.36
## 5 -0.916
## 6 -0.537
## 7 -0.202
## 8 0.0963
## 9 0.392
## 10 0.592
## # ℹ 22,591 more rows
##
## $diagnostics$random_slopes
## # A tibble: 239 × 3
## Code intercept slope
## <chr> <dbl> <dbl>
## 1 ABW 11.7 -0.136
## 2 AFG -24.4 0.184
## 3 AGO -19.7 0.0547
## 4 AIA 7.17 -0.0229
## 5 ALB 2.90 0.0679
## 6 AND 19.0 -0.121
## 7 ARE -3.02 0.203
## 8 ARG 6.59 -0.00627
## 9 ARM 5.88 -0.0808
## 10 ASM 12.0 -0.177
## # ℹ 229 more rows
##
##
## $volatility
## # A tibble: 239 × 3
## Code volatility trend
## <fct> <dbl> <dbl>
## 1 ABW 4.83 0.212
## 2 AFG 12.0 0.539
## 3 AGO 10.4 0.405
## 4 AIA 7.17 0.328
## 5 ALB 9.71 0.421
## 6 AND 5.04 0.228
## 7 ARE 12.4 0.561
## 8 ARG 15.1 0.345
## 9 ARM 6.09 0.268
## 10 ASM 3.99 0.169
## # ℹ 229 more rows
##
## $loess_fit
## Call:
## loess(formula = life_expectancy ~ time_since_1950, data = df_panel,
## span = 0.3)
##
## Number of Observations: 22601
## Equivalent Number of Parameters: 11.85
## Residual Standard Error: 10.09