library(tidyverse)   # dplyr, tidyr, ggplot2, readr, ...
library(scales)      # axis label formatting
library(corrplot)    # correlation heatmap
library(broom)       # tidy regression output
library(knitr)       # tables

theme_set(theme_minimal(base_size = 12))

1 Downloading data and loading into R

Two Kaggle datasets are used:

pop_raw <- read_csv("world_population.csv", show_col_types = FALSE)
co2_raw <- read_csv("CO2_emission.csv", show_col_types = FALSE,
                    name_repair = "unique_quiet")

2 Exploratory data analysis (World Population Dataset)

2.1 Variable names

names(pop_raw)
##  [1] "Rank"                        "CCA3"                       
##  [3] "Country/Territory"           "Capital"                    
##  [5] "Continent"                   "2022 Population"            
##  [7] "2020 Population"             "2015 Population"            
##  [9] "2010 Population"             "2000 Population"            
## [11] "1990 Population"             "1980 Population"            
## [13] "1970 Population"             "Area (km²)"                 
## [15] "Density (per km²)"           "Growth Rate"                
## [17] "World Population Percentage"

2.2 Top 5 rows

head(pop_raw, 5)
Rank CCA3 Country/Territory Capital Continent 2022 Population 2020 Population 2015 Population 2010 Population 2000 Population 1990 Population 1980 Population 1970 Population Area (km²) Density (per km²) Growth Rate World Population Percentage
36 AFG Afghanistan Kabul Asia 41128771 38972230 33753499 28189672 19542982 10694796 12486631 10752971 652230 63.0587 1.0257 0.52
138 ALB Albania Tirana Europe 2842321 2866849 2882481 2913399 3182021 3295066 2941651 2324731 28748 98.8702 0.9957 0.04
34 DZA Algeria Algiers Africa 44903225 43451666 39543154 35856344 30774621 25518074 18739378 13795915 2381741 18.8531 1.0164 0.56
213 ASM American Samoa Pago Pago Oceania 44273 46189 51368 54849 58230 47818 32886 27075 199 222.4774 0.9831 0.00
203 AND Andorra Andorra la Vella Europe 79824 77700 71746 71519 66097 53569 35611 19860 468 170.5641 1.0100 0.00

2.3 Bottom 10 rows

tail(pop_raw, 10)
Rank CCA3 Country/Territory Capital Continent 2022 Population 2020 Population 2015 Population 2010 Population 2000 Population 1990 Population 1980 Population 1970 Population Area (km²) Density (per km²) Growth Rate World Population Percentage
43 UZB Uzbekistan Tashkent Asia 34627652 33526656 30949417 28614227 24925554 20579100 15947129 12011361 447400 77.3975 1.0160 0.43
181 VUT Vanuatu Port-Vila Oceania 326740 311685 276438 245453 192074 150882 118156 87019 12189 26.8061 1.0238 0.00
234 VAT Vatican City Vatican City Europe 510 520 564 596 651 700 733 752 1 510.0000 0.9980 0.00
51 VEN Venezuela Caracas South America 28301696 28490453 30529716 28715022 24427729 19750579 15210443 11355475 916445 30.8820 1.0036 0.35
16 VNM Vietnam Hanoi Asia 98186856 96648685 92191398 87411012 79001142 66912613 52968270 41928849 331212 296.4472 1.0074 1.23
226 WLF Wallis and Futuna Mata-Utu Oceania 11572 11655 12182 13142 14723 13454 11315 9377 142 81.4930 0.9953 0.00
172 ESH Western Sahara El Aaiún Africa 575986 556048 491824 413296 270375 178529 116775 76371 266000 2.1654 1.0184 0.01
46 YEM Yemen Sanaa Asia 33696614 32284046 28516545 24743946 18628700 13375121 9204938 6843607 527968 63.8232 1.0217 0.42
63 ZMB Zambia Lusaka Africa 20017675 18927715 16248230 13792086 9891136 7686401 5720438 4281671 752612 26.5976 1.0280 0.25
74 ZWE Zimbabwe Harare Africa 16320537 15669666 14154937 12839771 11834676 10113893 7049926 5202918 390757 41.7665 1.0204 0.20

2.4 Data types

tibble(variable = names(pop_raw),
       type     = map_chr(pop_raw, ~ class(.x)[1]))
variable type
Rank numeric
CCA3 character
Country/Territory character
Capital character
Continent character
2022 Population numeric
2020 Population numeric
2015 Population numeric
2010 Population numeric
2000 Population numeric
1990 Population numeric
1980 Population numeric
1970 Population numeric
Area (km²) numeric
Density (per km²) numeric
Growth Rate numeric
World Population Percentage numeric

2.5 Shape of the dataset

dim(pop_raw)
## [1] 234  17
cat("Rows:", nrow(pop_raw), " | Columns:", ncol(pop_raw), "\n")
## Rows: 234  | Columns: 17

2.6 Duplicate rows

n_dup <- sum(duplicated(pop_raw))
cat("Number of duplicated rows:", n_dup, "\n")
## Number of duplicated rows: 0
pop <- distinct(pop_raw)   # drop duplicates (if any)
cat("Rows after dropping duplicates:", nrow(pop), "\n")
## Rows after dropping duplicates: 234

The dataset contains 0 duplicate rows, so all 234 rows are kept.

For easier coding, the columns are renamed to short snake_case names:

pop <- pop %>%
  rename(rank = Rank, cca3 = CCA3, country = `Country/Territory`,
         capital = Capital, continent = Continent,
         pop_2022 = `2022 Population`, pop_2020 = `2020 Population`,
         pop_2015 = `2015 Population`, pop_2010 = `2010 Population`,
         pop_2000 = `2000 Population`, pop_1990 = `1990 Population`,
         pop_1980 = `1980 Population`, pop_1970 = `1970 Population`,
         area = `Area (km²)`, density = `Density (per km²)`,
         growth_rate = `Growth Rate`,
         world_pop_pct = `World Population Percentage`)
glimpse(pop)
## Rows: 234
## Columns: 17
## $ rank          <dbl> 36, 138, 34, 213, 203, 42, 224, 201, 33, 140, 198, 55, 9…
## $ cca3          <chr> "AFG", "ALB", "DZA", "ASM", "AND", "AGO", "AIA", "ATG", …
## $ country       <chr> "Afghanistan", "Albania", "Algeria", "American Samoa", "…
## $ capital       <chr> "Kabul", "Tirana", "Algiers", "Pago Pago", "Andorra la V…
## $ continent     <chr> "Asia", "Europe", "Africa", "Oceania", "Europe", "Africa…
## $ pop_2022      <dbl> 41128771, 2842321, 44903225, 44273, 79824, 35588987, 158…
## $ pop_2020      <dbl> 38972230, 2866849, 43451666, 46189, 77700, 33428485, 155…
## $ pop_2015      <dbl> 33753499, 2882481, 39543154, 51368, 71746, 28127721, 145…
## $ pop_2010      <dbl> 28189672, 2913399, 35856344, 54849, 71519, 23364185, 131…
## $ pop_2000      <dbl> 19542982, 3182021, 30774621, 58230, 66097, 16394062, 110…
## $ pop_1990      <dbl> 10694796, 3295066, 25518074, 47818, 53569, 11828638, 831…
## $ pop_1980      <dbl> 12486631, 2941651, 18739378, 32886, 35611, 8330047, 6560…
## $ pop_1970      <dbl> 10752971, 2324731, 13795915, 27075, 19860, 6029700, 6283…
## $ area          <dbl> 652230, 28748, 2381741, 199, 468, 1246700, 91, 442, 2780…
## $ density       <dbl> 63.0587, 98.8702, 18.8531, 222.4774, 170.5641, 28.5466, …
## $ growth_rate   <dbl> 1.0257, 0.9957, 1.0164, 0.9831, 1.0100, 1.0315, 1.0066, …
## $ world_pop_pct <dbl> 0.52, 0.04, 0.56, 0.00, 0.00, 0.45, 0.00, 0.00, 0.57, 0.…

2.7 Missing values per column

missing_tbl <- tibble(variable = names(pop),
                      n_missing = colSums(is.na(pop)))
missing_tbl
variable n_missing
rank 0
cca3 0
country 0
capital 0
continent 0
pop_2022 0
pop_2020 0
pop_2015 0
pop_2010 0
pop_2000 0
pop_1990 0
pop_1980 0
pop_1970 0
area 0
density 0
growth_rate 0
world_pop_pct 0 
cat("Total missing values:", sum(missing_tbl$n_missing), "\n")
## Total missing values: 0

There are no missing values in the World Population Dataset, so no imputation is required. (Missing values do exist in the CO2 dataset; they are handled in Section 5.)

2.8 Outlier detection with box plots

quant_vars <- pop %>% select(where(is.numeric), -rank) %>% names()

pop %>%
  select(all_of(quant_vars)) %>%
  pivot_longer(everything(), names_to = "variable", values_to = "value") %>%
  ggplot(aes(x = "", y = value)) +
  geom_boxplot(fill = "#8ecae6", outlier.colour = "#d62828", outlier.alpha = 0.7) +
  facet_wrap(~ variable, scales = "free_y", ncol = 5) +
  scale_y_continuous(labels = label_number(scale_cut = cut_short_scale())) +
  labs(title = "Box plots of the quantitative variables",
       subtitle = "Red points are outliers (beyond 1.5 × IQR)", x = NULL, y = NULL)

Number of outliers per variable using the 1.5 × IQR rule:

iqr_outliers <- function(x) {
  q <- quantile(x, c(.25, .75), na.rm = TRUE)
  iqr <- diff(q)
  sum(x < q[1] - 1.5 * iqr | x > q[2] + 1.5 * iqr, na.rm = TRUE)
}
outlier_tbl <- pop %>%
  summarise(across(all_of(quant_vars), iqr_outliers)) %>%
  pivot_longer(everything(), names_to = "variable", values_to = "n_outliers")
outlier_tbl
variable n_outliers
pop_2022 25
pop_2020 26
pop_2015 27
pop_2010 26
pop_2000 29
pop_1990 31
pop_1980 32
pop_1970 29
area 28
density 22
growth_rate 3
world_pop_pct 25

Which countries are the outliers? For example, population and density:

pop %>% arrange(desc(pop_2022)) %>% select(country, pop_2022) %>% head(5)
country pop_2022
China 1425887337
India 1417173173
United States 338289857
Indonesia 275501339
Pakistan 235824862 
pop %>% arrange(desc(density))  %>% select(country, density)  %>% head(5)
country density
Macau 23172.267
Monaco 18234.500
Singapore 8416.463
Hong Kong 6783.392
Gibraltar 5441.500

2.9 Dealing with missing values and outliers

  • Missing values: none in this dataset, so nothing to impute.
  • Outliers: every variable is strongly right-skewed. The outliers are not data-entry errors: they are real countries such as China and India (population), Russia and Canada (area), or Macau and Monaco (density). Deleting them would remove the most important countries from a study of population and CO2. They are therefore kept, and their influence is controlled by:
    1. using log10 scales in plots and regression models, and
    2. reporting Spearman rank correlation (robust to outliers) next to Pearson.

The box plots after a log10 transform show far fewer extreme points:

pop %>%
  select(all_of(quant_vars), -growth_rate) %>%
  pivot_longer(everything(), names_to = "variable", values_to = "value") %>%
  filter(value > 0) %>%
  ggplot(aes(x = "", y = log10(value))) +
  geom_boxplot(fill = "#b7e4c7", outlier.colour = "#d62828") +
  facet_wrap(~ variable, scales = "free_y", ncol = 5) +
  labs(title = "Box plots after log10 transformation", x = NULL, y = "log10(value)")

3 Generating a new variable: projected 2030 population

The exponential growth model is

\[N_t = P\,e^{rt}\]

where \(P\) is the population at the end of 2022, \(r\) is the annual growth rate and \(t = 2030 - 2022 = 8\) years.

Important: in this dataset Growth Rate is stored as a growth factor (e.g. 1.0257), not as a rate. The annual rate is therefore \(r = \text{Growth Rate} - 1\) (e.g. 0.0257 = 2.57 %).

t_years <- 2030 - 2022

pop <- pop %>%
  mutate(r = growth_rate - 1,
         pop_2030 = round(pop_2022 * exp(r * t_years)),
         change_2022_2030 = pop_2030 - pop_2022)

pop %>%
  arrange(desc(pop_2030)) %>%
  select(country, pop_2022, r, pop_2030, change_2022_2030) %>%
  head(15) %>%
  kable(format.args = list(big.mark = ","), digits = 4,
        caption = "Projected 2030 population (top 15 countries)")
Projected 2030 population (top 15 countries)
country pop_2022 r pop_2030 change_2022_2030
India 1,417,173,173 0.0068 1,496,402,904 79,229,731
China 1,425,887,337 0.0000 1,425,887,337 0
United States 338,289,857 0.0038 348,731,782 10,441,925
Indonesia 275,501,339 0.0064 289,974,355 14,473,016
Pakistan 235,824,862 0.0191 274,757,645 38,932,783
Nigeria 218,541,212 0.0241 265,011,868 46,470,656
Brazil 215,313,498 0.0046 223,384,633 8,071,135
Bangladesh 171,186,372 0.0108 186,634,630 15,448,258
Ethiopia 123,379,924 0.0257 151,542,848 28,162,924
Russia 144,713,314 -0.0027 141,621,023 -3,092,291
Mexico 127,504,125 0.0063 134,095,029 6,590,904
Philippines 115,559,009 0.0147 129,980,092 14,421,083
DR Congo 99,010,212 0.0325 128,409,323 29,399,111
Egypt 110,990,103 0.0158 125,944,462 14,954,359
Japan 123,951,692 -0.0053 118,806,000 -5,145,692 
world_2022 <- sum(pop$pop_2022)
world_2030 <- sum(pop$pop_2030)
cat("World population 2022:", comma(world_2022), "\n")
## World population 2022: 7,973,413,042
cat("Projected world population 2030:", comma(world_2030), "\n")
## Projected world population 2030: 8,569,157,744

The model projects a world population of about 8.57 billion in 2030 (from 7.97 billion in 2022). India (1,496,402,904) overtakes China (1,425,887,337), because China’s growth rate is 0 (no growth), while India still grows at 0.68 % per year. Russia (r = -0.0027) is projected to shrink.

4 Value extraction and plots

4.1 Top 10 most populous countries in 2022

top10 <- pop %>% slice_max(pop_2022, n = 10) %>%
  select(cca3, country, continent, pop_2022)
top10 %>% kable(format.args = list(big.mark = ","))
cca3 country continent pop_2022
CHN China Asia 1,425,887,337
IND India Asia 1,417,173,173
USA United States North America 338,289,857
IDN Indonesia Asia 275,501,339
PAK Pakistan Asia 235,824,862
NGA Nigeria Africa 218,541,212
BRA Brazil South America 215,313,498
BGD Bangladesh Asia 171,186,372
RUS Russia Europe 144,713,314
MEX Mexico North America 127,504,125 
ggplot(top10, aes(x = reorder(country, pop_2022), y = pop_2022 / 1e6)) +
  geom_col(fill = "#219ebc") +
  geom_text(aes(label = comma(round(pop_2022 / 1e6, 1))), hjust = -0.1, size = 3.5) +
  coord_flip() +
  scale_y_continuous(expand = expansion(mult = c(0, 0.12))) +
  labs(title = "Top 10 most populous countries in 2022",
       x = NULL, y = "Population (millions)")

4.2 Population trend of the top 10 countries, 1990–2022

The dataset provides population for 1990, 2000, 2010, 2015, 2020 and 2022.

pop_trend <- pop %>%
  filter(cca3 %in% top10$cca3) %>%
  select(country, pop_1990, pop_2000, pop_2010, pop_2015, pop_2020, pop_2022) %>%
  pivot_longer(-country, names_to = "year", values_to = "population",
               names_prefix = "pop_", names_transform = as.integer)

ggplot(pop_trend, aes(year, population / 1e6, colour = reorder(country, -population))) +
  geom_line(linewidth = 1) +
  geom_point(size = 1.8) +
  scale_y_log10(labels = comma) +
  scale_x_continuous(breaks = c(1990, 2000, 2010, 2015, 2020, 2022)) +
  labs(title = "Population trend of the 10 most populous countries (1990-2022)",
       x = "Year", y = "Population (millions, log scale)", colour = "Country")

A log scale is used so that the smaller countries’ trends remain visible next to China and India. Growth since 1990:

pop %>%
  filter(cca3 %in% top10$cca3) %>%
  transmute(country, pop_1990, pop_2022,
            growth_pct = round(100 * (pop_2022 / pop_1990 - 1), 1)) %>%
  arrange(desc(growth_pct)) %>%
  kable(format.args = list(big.mark = ","))
country pop_1990 pop_2022 growth_pct
Nigeria 95,214,257 218,541,212 129.5
Pakistan 115,414,069 235,824,862 104.3
India 870,452,165 1,417,173,173 62.8
Bangladesh 107,147,651 171,186,372 59.8
Mexico 81,720,428 127,504,125 56.0
Indonesia 182,159,874 275,501,339 51.2
Brazil 150,706,446 215,313,498 42.9
United States 248,083,732 338,289,857 36.4
China 1,153,704,252 1,425,887,337 23.6
Russia 148,005,704 144,713,314 -2.2

4.3 CO2 emission trend of the top 10 countries, 1990–2019

The CO2 dataset contains a duplicated 2019 column (renamed 2019...35 on import), which is dropped. Values are CO2 emissions in metric tons per capita.

co2 <- co2_raw %>%
  select(-`2019...35`) %>%
  rename(country_name = `Country Name`, region = Region, `2019` = `2019...34`) %>%
  pivot_longer(`1990`:`2019`, names_to = "year", values_to = "co2_pc",
               names_transform = as.integer)

cat("Missing CO2 values:", sum(is.na(co2$co2_pc)), "of", nrow(co2), "\n")
## Missing CO2 values: 746 of 6450

The ten countries are matched with the population data using their ISO3 code. Names differ between the two datasets (for example Russia vs Russian Federation), so codes are more reliable.

co2_top10 <- co2 %>%
  filter(country_code %in% top10$cca3) %>%
  left_join(select(top10, cca3, country), by = c("country_code" = "cca3"))

cat("Missing values for the top 10:", sum(is.na(co2_top10$co2_pc)), "\n")
## Missing values for the top 10: 0
ggplot(co2_top10, aes(year, co2_pc, colour = reorder(country, -co2_pc))) +
  geom_line(linewidth = 1) +
  scale_x_continuous(breaks = seq(1990, 2019, 5)) +
  labs(title = "CO2 emissions of the 10 most populous countries (1990-2019)",
       x = "Year", y = "CO2 emissions (metric tons per capita)", colour = "Country")

co2_top10 %>%
  filter(year %in% c(1990, 2019)) %>%
  select(country, year, co2_pc) %>%
  pivot_wider(names_from = year, values_from = co2_pc, names_prefix = "co2_") %>%
  mutate(change_pct = round(100 * (co2_2019 / co2_1990 - 1), 1)) %>%
  arrange(desc(co2_2019)) %>%
  kable(digits = 2)
country co2_1990 co2_2019 change_pct
United States 19.41 14.67 -24.4
Russia 14.62 11.80 -19.3
China 1.91 7.61 297.3
Mexico 3.21 3.52 9.7
Indonesia 0.82 2.29 179.7
Brazil 1.33 2.06 54.9
India 0.65 1.80 178.5
Pakistan 0.55 0.88 60.5
Nigeria 0.76 0.57 -24.9
Bangladesh 0.11 0.56 398.4

5 Correlation analysis

5.1 Correlations between Area, Density, Growth Rate and World Population Percentage

cor_vars <- pop %>%
  select(Area = area, Density = density,
         `Growth Rate` = growth_rate, `World Pop %` = world_pop_pct)

cor_pearson  <- cor(cor_vars, method = "pearson")
cor_spearman <- cor(cor_vars, method = "spearman")

kable(round(cor_pearson, 3),  caption = "Pearson correlation matrix")
Pearson correlation matrix
Area Density Growth Rate World Pop %
Area 1.000 -0.063 -0.014 0.453
Density -0.063 1.000 -0.070 -0.028
Growth Rate -0.014 -0.070 1.000 -0.021
World Pop % 0.453 -0.028 -0.021 1.000 
kable(round(cor_spearman, 3), caption = "Spearman (rank) correlation matrix")
Spearman (rank) correlation matrix
Area Density Growth Rate World Pop %
Area 1.000 -0.584 0.356 0.841
Density -0.584 1.000 -0.222 -0.154
Growth Rate 0.356 -0.222 1.000 0.278
World Pop % 0.841 -0.154 0.278 1.000 
par(mfrow = c(1, 2))
corrplot(cor_pearson, method = "color", type = "upper", addCoef.col = "black",
         tl.col = "black", tl.srt = 45, col = COL2("RdBu", 200),
         title = "Pearson", mar = c(0, 0, 2, 0))
corrplot(cor_spearman, method = "color", type = "upper", addCoef.col = "black",
         tl.col = "black", tl.srt = 45, col = COL2("RdBu", 200),
         title = "Spearman", mar = c(0, 0, 2, 0))

par(mfrow = c(1, 1))

Significance test of each pair (Pearson):

combn(names(cor_vars), 2, simplify = FALSE) %>%
  map_dfr(function(p) {
    ct <- cor.test(cor_vars[[p[1]]], cor_vars[[p[2]]])
    tibble(var1 = p[1], var2 = p[2], r = ct$estimate, p_value = ct$p.value)
  }) %>%
  kable(digits = 4)
var1 var2 r p_value
Area Density -0.0631 0.3363
Area Growth Rate -0.0140 0.8317
Area World Pop % 0.4533 0.0000
Density Growth Rate -0.0698 0.2880
Density World Pop % -0.0276 0.6740
Growth Rate World Pop % -0.0209 0.7501

5.2 Population (2022) vs CO2 emission (2019)

5.2.1 Merging the datasets

The task asks for a merge on country name. First, a merge on exact names:

co2_2019 <- co2 %>% filter(year == 2019) %>%
  select(country_name, country_code, region, co2_2019 = co2_pc)

by_name <- inner_join(pop, co2_2019, by = c("country" = "country_name"))
by_code <- inner_join(pop, co2_2019, by = c("cca3" = "country_code"))
cat("Matched on exact country name:", nrow(by_name), "\n")
## Matched on exact country name: 186
cat("Matched on ISO3 country code :", nrow(by_code), "\n")
## Matched on ISO3 country code : 214

Some countries are lost because the two datasets name them differently:

by_code %>% filter(country != country_name) %>%
  select(population_dataset = country, co2_dataset = country_name)
population_dataset co2_dataset
Bahamas Bahamas, The
Brunei Brunei Darussalam
Cape Verde Cabo Verde
DR Congo Congo, Dem. Rep.
Egypt Egypt, Arab Rep.
Gambia Gambia, The
Hong Kong Hong Kong SAR, China
Iran Iran, Islamic Rep.
Ivory Coast Cote d’Ivoire
Kyrgyzstan Kyrgyz Republic
Laos Lao PDR
Macau Macao SAR, China
Micronesia Micronesia, Fed. Sts.
North Korea Korea, Dem. People’s Rep.
Palestine West Bank and Gaza
Republic of the Congo Congo, Rep.
Russia Russian Federation
Saint Kitts and Nevis St. Kitts and Nevis
Saint Lucia St. Lucia
Saint Martin St. Martin (French part)
Saint Vincent and the Grenadines St. Vincent and the Grenadines
Sint Maarten Sint Maarten (Dutch part)
Slovakia Slovak Republic
South Korea Korea, Rep.
Syria Syrian Arab Republic
Turkey Turkiye
United States Virgin Islands Virgin Islands (U.S.)
Venezuela Venezuela, RB
Yemen Yemen, Rep.

To keep these countries, the final merge uses the ISO3 code, which identifies the same country in both datasets. Rows with no 2019 CO2 value are then removed. A second new variable is also created: total CO2 emissions = per-capita emission × population.

merged <- by_code %>%
  select(cca3, country, continent, pop_2022, co2_2019) %>%
  filter(!is.na(co2_2019)) %>%
  mutate(co2_total_mt = co2_2019 * pop_2022 / 1e6)   # million metric tons (approx.)

cat("Countries in merged dataset:", nrow(merged), "\n")
## Countries in merged dataset: 191
head(merged, 10)
cca3 country continent pop_2022 co2_2019 co2_total_mt
AFG Afghanistan Asia 41128771 0.1598244 6.5733800
ALB Albania Europe 2842321 1.6922483 4.8099129
DZA Algeria Africa 44903225 3.9776505 178.6093335
AND Andorra Europe 79824 6.4812174 0.5173567
AGO Angola Africa 35588987 0.7921371 28.1913559
ATG Antigua and Barbuda North America 93763 5.3544765 0.5020518
ARG Argentina South America 45510318 3.7406503 170.2381843
ARM Armenia Asia 2780469 2.0860607 5.8002270
AUS Australia Oceania 26177413 15.2382672 398.8984126
AUT Austria Europe 8939617 7.2939843 65.2054256

5.2.2 Visualisation

ggplot(merged, aes(pop_2022, co2_2019)) +
  geom_point(aes(colour = continent), alpha = 0.75, size = 2) +
  geom_smooth(method = "lm", se = TRUE, colour = "black") +
  geom_text(data = filter(merged, cca3 %in% top10$cca3), aes(label = country),
            size = 3, vjust = -0.8, check_overlap = TRUE) +
  scale_x_log10(labels = label_number(scale_cut = cut_short_scale())) +
  labs(title = "2022 population vs 2019 CO2 emissions per capita",
       x = "Population 2022 (log scale)", y = "CO2 (metric tons per capita)",
       colour = "Continent")

ggplot(merged, aes(pop_2022, co2_total_mt)) +
  geom_point(aes(colour = continent), alpha = 0.75, size = 2) +
  geom_smooth(method = "lm", se = TRUE, colour = "black") +
  geom_text(data = filter(merged, cca3 %in% top10$cca3), aes(label = country),
            size = 3, vjust = -0.8, check_overlap = TRUE) +
  scale_x_log10(labels = label_number(scale_cut = cut_short_scale())) +
  scale_y_log10(labels = comma) +
  labs(title = "2022 population vs total CO2 emissions (log-log)",
       x = "Population 2022 (log scale)", y = "Total CO2 (million metric tons, log scale)",
       colour = "Continent")

5.2.3 Correlation tests

ct_pearson      <- cor.test(merged$pop_2022, merged$co2_2019, method = "pearson")
ct_spearman     <- cor.test(merged$pop_2022, merged$co2_2019, method = "spearman", exact = FALSE)
ct_total_log    <- cor.test(log10(merged$pop_2022), log10(merged$co2_total_mt))

tibble(
  relationship = c("Population vs CO2 per capita (Pearson)",
                   "Population vs CO2 per capita (Spearman)",
                   "log Population vs log Total CO2 (Pearson)"),
  correlation = c(ct_pearson$estimate, ct_spearman$estimate, ct_total_log$estimate),
  p_value     = c(ct_pearson$p.value, ct_spearman$p.value, ct_total_log$p.value)
) %>% kable(digits = 4)
relationship correlation p_value
Population vs CO2 per capita (Pearson) 0.0145 0.8426
Population vs CO2 per capita (Spearman) -0.1103 0.1288
log Population vs log Total CO2 (Pearson) 0.8182 0.0000

5.2.4 Regression

m_pc    <- lm(co2_2019 ~ log10(pop_2022), data = merged)
m_total <- lm(log10(co2_total_mt) ~ log10(pop_2022), data = merged)

tidy(m_pc)    %>% kable(digits = 4, caption = "Model 1: CO2 per capita ~ log10(population)")
Model 1: CO2 per capita ~ log10(population)
term estimate std.error statistic p.value
(Intercept) 5.7780 2.4914 2.3192 0.0215
log10(pop_2022) -0.2437 0.3616 -0.6738 0.5013 
glance(m_pc)  %>% select(r.squared, adj.r.squared, p.value, nobs) %>% kable(digits = 4)
r.squared adj.r.squared p.value nobs
0.0024 -0.0029 0.5013 191 
tidy(m_total)   %>% kable(digits = 4, caption = "Model 2: log10(total CO2) ~ log10(population)")
Model 2: log10(total CO2) ~ log10(population)
term estimate std.error statistic p.value
(Intercept) -5.1042 0.3214 -15.8795 0
log10(pop_2022) 0.9127 0.0467 19.5630 0 
glance(m_total) %>% select(r.squared, adj.r.squared, p.value, nobs) %>% kable(digits = 4)
r.squared adj.r.squared p.value nobs
0.6694 0.6677 0 191

Interpretation.

  • Per capita: the correlation between population and CO2 per capita is very weak (Pearson r = 0.014, p = 0.843; Spearman ρ = -0.11, p = 0.129). Model 1 explains only 0.2 % of the variation. A country’s size does not tell us how much each person emits. That depends on income, industry and the energy mix (e.g. Qatar is small but has very high per-capita emissions, while Nigeria and Bangladesh are large but low-emitting).
  • Total emissions: total CO2 rises strongly with population (r = 0.818 on the log–log scale, R² = 0.669). The slope of 0.913 means that a 1 % larger population is associated with about 0.91 % more total CO2.

6 Comparing CO2 emissions across continents

As required, the 2019 CO2 emissions (metric tons per capita) in the merged dataset are summed by continent. The mean per country and the population-weighted mean are shown too, because a simple sum also depends on how many countries each continent has.

continent_co2 <- merged %>%
  group_by(continent) %>%
  summarise(n_countries   = n(),
            total_co2_pc  = sum(co2_2019),
            mean_co2_pc   = mean(co2_2019),
            weighted_co2_pc = weighted.mean(co2_2019, pop_2022),
            total_co2_mt  = sum(co2_total_mt)) %>%
  arrange(desc(total_co2_pc)) %>%
  mutate(rank = row_number(), .before = 1)

continent_co2 %>% kable(digits = 2, format.args = list(big.mark = ","))
rank continent n_countries total_co2_pc mean_co2_pc weighted_co2_pc total_co2_mt
1 Asia 46 290.55 6.32 4.32 20,247.94
2 Europe 42 246.67 5.87 6.97 5,179.65
3 North America 23 96.90 4.21 10.35 6,164.84
4 Africa 54 63.73 1.18 1.04 1,486.85
5 Oceania 14 54.78 3.91 10.12 446.64
6 South America 12 33.36 2.78 2.40 1,046.78 
ggplot(continent_co2, aes(reorder(continent, total_co2_pc), total_co2_pc, fill = continent)) +
  geom_col(show.legend = FALSE) +
  geom_text(aes(label = round(total_co2_pc, 2)), hjust = -0.1, size = 4) +
  coord_flip() +
  scale_y_continuous(expand = expansion(mult = c(0, 0.12))) +
  labs(title = "Sum of 2019 CO2 emissions by continent",
       x = NULL, y = "Sum of CO2 emissions (metric tons per capita)")

Answers (from the graph):

  1. First continent with the highest 2019 CO2 emission: Asia (290.55 metric tons per capita).
  2. Third continent: North America, with a total of 96.9 metric tons per capita.
  3. Last continent: South America, with a total of 33.36 metric tons per capita.

Note: summing per-capita values favours continents with many countries. The population-weighted mean column gives a fairer comparison of how much an average person emits on each continent.

7 Session information

sessionInfo()
## R version 4.6.1 (2026-06-24)
## Platform: aarch64-apple-darwin23
## Running under: macOS Tahoe 26.3.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.6/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: Africa/Kigali
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] knitr_1.52      broom_1.0.13    corrplot_0.95   scales_1.4.0   
##  [5] lubridate_1.9.5 forcats_1.0.1   stringr_1.6.0   dplyr_1.2.1    
##  [9] purrr_1.2.2     readr_2.2.0     tidyr_1.3.2     tibble_3.3.1   
## [13] ggplot2_4.0.3   tidyverse_2.0.0
## 
## loaded via a namespace (and not attached):
##  [1] sass_0.4.10        generics_0.1.4     lattice_0.22-9     stringi_1.8.9     
##  [5] hms_1.1.4          digest_0.6.39      magrittr_2.0.5     evaluate_1.0.5    
##  [9] grid_4.6.1         timechange_0.4.0   RColorBrewer_1.1-3 fastmap_1.2.0     
## [13] Matrix_1.7-5       jsonlite_2.0.0     backports_1.5.1    mgcv_1.9-4        
## [17] jquerylib_0.1.4    cli_3.6.6          rlang_1.3.0        crayon_1.5.3      
## [21] splines_4.6.1      bit64_4.8.6        withr_3.0.3        cachem_1.1.0      
## [25] yaml_2.3.12        otel_0.2.0         tools_4.6.1        parallel_4.6.1    
## [29] tzdb_0.5.0         vctrs_0.7.3        R6_2.6.1           stats4_4.6.1      
## [33] lifecycle_1.0.5    bit_4.6.0          vroom_1.7.1        pkgconfig_2.0.3   
## [37] pillar_1.11.1      bslib_0.12.0       gtable_0.3.6       glue_1.8.1        
## [41] xfun_0.62          tidyselect_1.2.1   rstudioapi_0.19.0  farver_2.1.2      
## [45] nlme_3.1-169       htmltools_0.5.9    rmarkdown_2.32     labeling_0.4.3    
## [49] compiler_4.6.1     S7_0.2.2