Analysis

Set up

library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   4.0.0     ✔ tibble    3.3.0
✔ lubridate 1.9.4     ✔ tidyr     1.3.1
✔ purrr     1.1.0     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(ggplot2)
library(plotly)

Attaching package: 'plotly'

The following object is masked from 'package:ggplot2':

    last_plot

The following object is masked from 'package:stats':

    filter

The following object is masked from 'package:graphics':

    layout

# Load pre-cleaned data
load("data/clean.RData") 

Q1: Who dominates

What are the top 10 countries that have recorded the highest average Air Quality Index (AQI) during the entire 2020-2024 period, and how is this ranking distributed across their respective Region and Income Level?

First graph: a bar chart of top 10 highest AQI countries

# Step1: Find the top 10 countries that have recorded the highest avg air quality
df1 <- df_clean2 |>
  group_by(country) |>
  summarise(avg_air_quality = mean(air_quality, na.rm = TRUE)) |>
  arrange(-avg_air_quality)

df1

# A tibble: 25 × 2
   country     avg_air_quality
   <chr>                 <dbl>
 1 Kenya                 160. 
 2 Egypt                 157. 
 3 India                 156. 
 4 Nigeria               155. 
 5 Pakistan              155. 
 6 Vietnam               154. 
 7 Philippines           154. 
 8 Bangladesh            154. 
 9 Indonesia             153. 
10 Germany                71.8
# ℹ 15 more rows

# Step2: Draw the bar chart about top 10 highest AQI countries
p1 <- df1|>
  head(10) |>
  arrange(-avg_air_quality) |>
  mutate(country = factor(country, levels = country)) |>
  ggplot(aes(x = country,
            y = avg_air_quality,
             fill = country)) +
  geom_col(width = 0.88) +
  geom_text(aes(label = round(avg_air_quality, 1)), 
            hjust = -0.12, 
            size = 2.5,
            color = "black",
            fontface = "bold") +
  geom_text(aes(label = country), 
            hjust = 1.1, 
            size = 3.5,
            color = "white",
            fontface = "bold") +
  coord_flip() +
  scale_x_discrete(limits = rev) +
  labs(x= "Country",
       y = "Average Air Quality Index",
       title = "Top 10 Highest Air Quality Index Countries in the World during 2020~2024",
       caption = "Source: Kaggle | Author: Gu Raomeng") +
  scale_fill_manual(values = rev(c("#bdbdbd", "#b0b0b0", "#969696", "#8c8c8c", "#737373", "#6a6a6a", "#525252", "#252525", "#1a1a1a", "#000000"))) +
  theme_light() +
  theme(panel.grid = element_blank(), 
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        legend.position = "none",
        plot.title = element_text(face = "bold",
                                  size = 13,
                                  hjust = 0.5,
                                  margin = margin(b = 12)),
        axis.title.x = element_text(face = "bold",
                                    size = 12,
                                    color = "#000000",
                                    margin = margin(t = 12)),
        axis.title.y = element_text(face = "bold",
                                    size = 12,
                                    color = "#000000",
                                    margin = margin(r = 12)),
        plot.caption = element_text(margin = margin(t = 12)))
print(p1)

# Step3: save the graph
ggsave("images/p1.png", p1, width = 7, height = 4.2, dpi = 300)

Second graph: a world map but only focus on the top 10 countries mentioned in the first graph, this map will use different color to show the income level of each country

# how is this ranking distributed across their respective Region and Income Level
# Step1: make a table of 10 countries and their corresponding region and income level
top10_countries <- c("Kenya", "Egypt", "India", "Nigeria", "Pakistan", "Vietnam", "Philippines", "Bangladesh", "Indonesia", "Germany")
df2 <- df_clean2 |>
  select( region, income_level, country) |>
  filter(country %in% top10_countries) |>
  distinct()

glimpse(df2)

Rows: 10
Columns: 3
$ region       <chr> "South Asia", "Africa", "Europe", "Africa", "Southeast As…
$ income_level <chr> "Lower-Middle", "Lower-Middle", "High", "Lower-Middle", "…
$ country      <chr> "India", "Nigeria", "Germany", "Kenya", "Indonesia", "Pak…

# draw the world map of highlight these top 10 highest air quality index countries with their income level
# Step2: Load world map
library(rnaturalearth)
library(rnaturalearthdata)

Attaching package: 'rnaturalearthdata'

The following object is masked from 'package:rnaturalearth':

    countries110

library(sf)

Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUE

world <- ne_countries(scale = "medium", returnclass = "sf")
world_top10 <- world |>
  filter(name %in% df2$country) |>
  left_join(df2, by = c("name" = "country"))

# Step3: Plot map
map1 <- ggplot() +
  geom_sf(data = world, 
          fill = "gray90", 
          color = "white", 
          size = 0.1) +
  geom_sf(data = world_top10,
          aes(fill = income_level),
          color = "gray45",
          size = 0.3) +
  scale_fill_manual(values = c("High" = "#2E8B57",
                               "Upper-Middle" = "#FFA500",
                               "Lower-Middle" = "#CD5C5C", 
                               "Low" = "#8B0000"),
                    name = "Income Level",
                    labels = c("High", "Upper-Middle", "Lower-Middle", "Low")) +
  labs(title = "Top 10 Highest AQI Countries by Income Level (2020-2024)", 
       caption = "Source: Kaggle | Author: Gu Raomeng", 
       fill = "Income Level") +
  coord_sf(xlim = c(-180, 180), ylim = c(-60, 80), expand = FALSE) +
  theme_minimal() +
  theme(panel.grid = element_blank(), 
        axis.text.y = element_blank(),
        axis.text.x = element_blank(),
        axis.ticks.y = element_blank(),
        axis.ticks.x = element_blank(),
        legend.position = "right", 
        plot.title = element_text(face = "bold",
                                  size = 13,
                                  hjust = 0.5,
                                  margin = margin(b = 12)),
        axis.title.x = element_text(face = "bold",
                                    size = 12,
                                    color = "#000000",
                                    margin = margin(t = 12)),
        axis.title.y = element_text(face = "bold",
                                    size = 12,
                                    color = "#000000",
                                    margin = margin(r = 12)),
        plot.caption = element_text(margin = margin(t = 12)))
print(map1) 

# Step4: save the map into PNG
ggsave("images/map1.png", map1, width = 7, height = 4.2, dpi = 300)

Q2: How Air Quality Affects Respiratory Health

Is there a measurable linear correlation between the average Air Quality Index and the average Respiratory Disease Rate during 2020-2024?

Third graph: a scatter plot with total 25 countries’ data.

# Step1: calculate the avg AQI and avg respiratory disease rate
df_q2 <- df_clean2 |>
  select(country, air_quality, respiratory_disease) |>
  group_by(country) |>
  summarise(avg_air_quality = mean(air_quality, na.rm = TRUE),
            avg_respiratory_disease = mean(respiratory_disease, na.rm = TRUE)) |>
  ungroup()

glimpse(df_q2)

Rows: 25
Columns: 3
$ country                 <chr> "Argentina", "Australia", "Bangladesh", "Brazi…
$ avg_air_quality         <dbl> 68.83142, 66.05364, 153.80843, 65.65517, 67.40…
$ avg_respiratory_disease <dbl> 65.00881, 62.93180, 81.26667, 62.33831, 64.100…

# Step2: draw scatter plot
p_scatter <- ggplot(df_q2, aes(x = avg_air_quality, y = avg_respiratory_disease)) + 
  geom_point(aes(color = country), 
             size = 3.5,
             alpha = 0.7) +
  geom_smooth(method = "lm", 
              se = FALSE, 
              color = "red",
              size = 0.8) +
  labs(x = "Average Air Quality Index (AQI)",
       y = "Average Respiratory Disease Rate (%)",
       title = "The Relationship Between AQI and Respiratory Disease Rate (2020-2024)",
       caption = "Source: Kaggle | Author: Gu Raomeng") +
  theme_minimal() +
  theme(legend.position = "right",
        plot.title = element_text(face = "bold", 
                              size = 12, 
                              hjust = 0.5),
        axis.title.x = element_text(face = "bold",
                                    size = 10,
                                    margin = margin(t = 10)),
        axis.title.y = element_text(face = "bold",
                                    size = 10,
                                    margin = margin(r = 10)),
        plot.caption = element_text(margin = margin(t = 10)),)

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.

print(p_scatter)

`geom_smooth()` using formula = 'y ~ x'

p_interactive <- ggplotly(p_scatter)

`geom_smooth()` using formula = 'y ~ x'

library(htmlwidgets)
saveWidget(p_interactive, file = "out/p_interactive.html")

Q3: What Protects Respiratory Health

Explore the relationship between the average healthcare access index and average respiratory disease rate in three different income levels countries during 2020-2024.

Fourth graph: 2 line plot.

# Step1: calculate all useful data
df4 <- df_clean2 |>
  select(country, healthcare, respiratory_disease, year, income_level) |>
  group_by(country, year,income_level) |>
  summarise(avg_healthcare = mean(healthcare, na.rm = TRUE),
            avg_respiratory_disease = mean(respiratory_disease, na.rm = TRUE)) |>
  ungroup()

`summarise()` has grouped output by 'country', 'year'. You can override using
the `.groups` argument.

glimpse(df4)

Rows: 125
Columns: 5
$ country                 <chr> "Argentina", "Argentina", "Argentina", "Argent…
$ year                    <dbl> 2020, 2021, 2022, 2023, 2024, 2020, 2021, 2022…
$ income_level            <chr> "Upper-Middle", "Upper-Middle", "Upper-Middle"…
$ avg_healthcare          <dbl> 61.54423, 61.51538, 61.88269, 61.40943, 61.315…
$ avg_respiratory_disease <dbl> 65.59615, 66.91346, 66.43077, 64.51509, 61.598…

# Step2: smaller the dataset with 3 explored countries, because United States and Kenya are the two extreme countries in Q2's plot, Argentina is a randomly upper-middle income level country
target_countries <- c("United States", "Argentina", "Kenya")

df_3countries <- df4 |>
  filter(country %in% target_countries) |>
  select(avg_healthcare, avg_respiratory_disease, year, country, income_level)

glimpse(df_3countries)

Rows: 15
Columns: 5
$ avg_healthcare          <dbl> 61.54423, 61.51538, 61.88269, 61.40943, 61.315…
$ avg_respiratory_disease <dbl> 65.59615, 66.91346, 66.43077, 64.51509, 61.598…
$ year                    <dbl> 2020, 2021, 2022, 2023, 2024, 2020, 2021, 2022…
$ country                 <chr> "Argentina", "Argentina", "Argentina", "Argent…
$ income_level            <chr> "Upper-Middle", "Upper-Middle", "Upper-Middle"…

# Step3: set colors for each countries and income levels
country_colors <- c("United States" = "#91cf60",
                    "Argentina" = "#ffeda0",
                    "Kenya" = "#fc8d59")
income_colors <- c("High" = "#2E8B57",
                   "Upper-Middle" = "#FFA500",
                   "Lower-Middle" = "#CD5C5C")

# Step4: Plot line graph(healthcare)
p_health <- ggplot(df_3countries, aes(x = year, y = avg_healthcare, group = country)) +
  geom_line(aes(color = country), linewidth = 1.5) +
  geom_point(aes(color = country), size = 3) +
  scale_color_manual(values = country_colors) +
  labs(
    title = "Healthcare Access Index Trends (2020-2024)",
    y = "Healthcare Access Index",
    x = "Year") +
  theme_minimal() +
  theme(legend.position = "right",
        plot.title = element_text(face = "bold",
                                  size = 12,
                                  hjust = 0.5),
        panel.grid.minor = element_blank(),
        axis.title.x = element_text(face = "bold",
                                    size = 10,
                                    margin = margin(t = 10)),
        axis.title.y = element_text(face = "bold",
                                    size = 10,
                                    margin = margin(r = 10)))

print(p_health)

# Step5: Plot line graph(respiratory disease rate)
p_respiratory <- ggplot(df_3countries, aes(x = year, y = avg_respiratory_disease, group = country)) +
  geom_line(aes(color = country), linewidth = 1.5) +
  geom_point(aes(color = country), size = 3) +
  scale_color_manual(values = country_colors) +
  labs(
    title = "Respiratory Disease Rate Trends (2020-2024)",
    y = "Respiratory Disease Rate (%)",
    x = "Year") +
  theme_minimal() +
  theme(legend.position = "right",
        plot.title = element_text(face = "bold",
                                  size = 12,
                                  hjust = 0.5),
        panel.grid.minor= element_blank(),
        axis.title.x = element_text(face = "bold",
                                    size = 10,
                                    margin = margin(t = 10)),
        axis.title.y = element_text(face = "bold",
                                    size = 10,
                                    margin = margin(r = 10)))

print(p_respiratory)

# Step6: combine 2 plots together
library(patchwork)

p_combined <- p_health / p_respiratory +
  plot_annotation(
    title = "Relationship between Healthcare Access and Respiratory Disease",
    subtitle = "With Countries in Three Income Level",
    caption = "Source: Kaggle | Author: Gu Raomeng",
    theme = theme(
      plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
      plot.subtitle = element_text(size = 13, hjust = 0.5, color = "gray25")))
print(p_combined)

# Step7: save the combined graph into correct size
ggsave("images/combined_plot.png", p_combined, 
       width = 10, height = 6, dpi = 300)