14 LDS Enclave Macroeconomy

14_macro_econ.r

# Part 5 of American Communities Project - GDP Analysis

library(tidyverse)
library(openxlsx)
library(scales)

# Outline:
# I. Load data and clean
# II. Create visualizations
# IIa. Rank communities by GDP with latest data
# IIb. Visualize GDP over time by community
# IIc. Rank GDP across countries and American Communities

# links to data: 
# county level gdp https://www.bea.gov/data/gdp/gdp-county-metro-and-other-areas
# population https://www.census.gov/data/tables/time-series/demo/popest/2020s-counties-total.html
# country gdp https://data.worldbank.org/indicator/NY.GDP.MKTP.CD

# population data
population_data_1 <- read.csv('./data/macroecon/co-est2010-alldata.csv')
population_data_2 <- read.csv('./data/macroecon/co-est2020-alldata.csv')
population_data_3 <- read.csv('./data/macroecon/co-est2024-alldata.csv')

population_data_1 <- population_data_1 |> 
    filter(COUNTY != 0) |> 
    mutate(fips = sprintf("%02d%03d", STATE, COUNTY)) |> 
    select(fips, tidyselect::contains('popestimate')) |> 
    select(-(POPESTIMATE2010))


population_data_2 <- population_data_2 |> 
    filter(COUNTY != 0) |> 
    mutate(fips = sprintf("%02d%03d", STATE, COUNTY)) |> 
    select(fips, tidyselect::contains('popestimate')) |> 
    select(-(POPESTIMATE2020))

population_data_3 <- population_data_3 |> 
    filter(COUNTY != 0) |> 
    mutate(fips = sprintf("%02d%03d", STATE, COUNTY)) |> 
    select(fips, tidyselect::contains('popestimate'))

population_data <- population_data_1 |> 
    inner_join(population_data_2, by = 'fips') |> 
    inner_join(population_data_3, by = 'fips') |> 
    mutate(fips = stringr::str_trim(fips))

#county gdp data
county_gdp_data <- read.csv('./data/macroecon/county_gdp.csv') |> 
    filter(Description == "Real GDP (thousands of chained 2017 dollars) ") |> 
    rename(fips = GeoFIPS) |>
    select(fips, tidyselect::contains('X')) |> 
    mutate(fips = stringr::str_trim(fips))

county_type <- openxlsx::read.xlsx('./data/macroecon/2023-Typology.xlsx') |> 
    mutate(fips = stringr::str_sub(GEO_ID, -5, -1)) |> 
    dplyr::rename(Typology = `2023.Typology`) |> 
    select(fips, Typology)

gdp_pop_data <- population_data |> 
    left_join(county_gdp_data, by = 'fips') |> 
    inner_join(county_type, by = 'fips') |> 
    mutate(across(starts_with("X")|starts_with("POPESTIMATE"), ~ as.numeric(.))) |> 
    mutate(across(starts_with("X"), ~ . *1000))

# Overall GDP Compared To American Communities

gdp_pop_data |> 
    group_by(Typology) |>
    summarize(across(starts_with("X"), sum, na.rm = TRUE)) |> 
    pivot_longer(
        cols = starts_with("X"),
        names_to = "year",
        values_to = "gdp"
    ) |> 
    mutate(year = as.numeric(str_remove(year, "X"))) |> 
    print(n=500)

plot <- gdp_pop_data |> 
    group_by(Typology) |>
    summarize(across(starts_with("X"), sum, na.rm = TRUE)) |> 
    pivot_longer(
        cols = starts_with("X"),
        names_to = "year",
        values_to = "gdp"
    ) |> 
    mutate(year = as.numeric(str_remove(year, "X"))) %>%
    ggplot(aes(x = year, y = gdp, color = Typology)) +
    geom_line() +
    geom_label(data = . %>% filter(year == max(year)), aes(label = Typology), size = 1.5, nudge_x = -0.4) +
    labs(
        title = "GDP Over Time by Community Type",
        x = "Year",
        y = "GDP (in 2017 dollars)",
        caption = "@mormon_metrics\ndata: US Bureau of Economic Analysis"
    ) + 
    theme_minimal() +
    theme(
        axis.ticks = element_blank(),
        axis.title = element_text(size = 14),
        axis.text.y = element_text(size = 10),
        axis.text.x = element_text(size = 10, margin = margin(t = -5)),
        plot.background = element_rect(fill = "grey95"),  # Change entire plot background color here
        panel.background = element_blank(),
        text = element_text(face = "bold", family = "Cairo"),
        plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
        plot.title.position = "plot",
        plot.subtitle.position = "plot",
        legend.title = element_blank(),
        legend.background = element_blank(),
        legend.box.background = element_blank(),
        legend.key = element_blank(),
        panel.grid = element_blank(),
        legend.position = "none",
        panel.grid.major.y = element_line(color = "lightgrey", linetype = 'solid'),
        plot.caption = element_text(size = 8,family = "Cairo"))+
    scale_x_continuous(
        breaks = seq(2000,2020, 5)
    ) +
    scale_y_continuous(labels = scales::label_number(scale = 1e-12, suffix = "T"))

ggsave("./images/14_macro_1.png", plot, width = 1500, height = 2500, units = "px", dpi = 300)

# GDP Per Capita Over Time
plot <- gdp_pop_data |> 
    group_by(Typology) |>
    summarize(across(starts_with("X")|starts_with("POPESTIMATE"), sum, na.rm = TRUE)) |> 
    mutate(
        across(
            .cols = all_of(paste0("X", 2001:2023)),
            .fns = ~ .x / get(str_replace(cur_column(), "X", "POPESTIMATE")),
            .names = "gdppc{.col}"
        )
    ) |> 
    select(Typology, starts_with("gdppc")) |>
    pivot_longer(
        cols = starts_with("gdppc"),
        names_to = "year",
        values_to = "gdp_per_capita"
    ) |>
    mutate(year = as.numeric(str_remove(year, "gdppcX"))) %>%
    ggplot(aes(x = year, y = gdp_per_capita, color = Typology)) +
    geom_line() +
    geom_label(data = . %>% filter(year == max(year)), aes(label = Typology), size = 2, nudge_x = -1.3) +
    labs(
        title = "GDP Per Capita Over Time by Community Type",
        x = "Year",
        y = "GDP (in 2017 dollars)",
        caption = "@mormon_metrics\ndata: US Bureau of Economic Analysis"
    ) + 
    theme_minimal() +
    theme(
        axis.ticks = element_blank(),
        axis.title = element_text(size = 14),
        axis.text.y = element_text(size = 10),
        axis.text.x = element_text(size = 10, margin = margin(t = -5)),
        plot.background = element_rect(fill = "grey95"),  # Change entire plot background color here
        panel.background = element_blank(),
        text = element_text(face = "bold", family = "Cairo"),
        plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
        plot.title.position = "plot",
        plot.subtitle.position = "plot",
        legend.title = element_blank(),
        legend.background = element_blank(),
        legend.box.background = element_blank(),
        legend.key = element_blank(),
        panel.grid = element_blank(),
        legend.position = "none",
        panel.grid.major.y = element_line(color = "lightgrey", linetype = 'solid'),
        plot.caption = element_text(size = 8,family = "Cairo"))+
    scale_x_continuous(
        breaks = seq(2000,2020, 5)
    ) +
    scale_y_continuous(labels = scales::label_number(), breaks = seq(0, 100000, 10000), limits = c(0, 100000))

ggsave("./images/14_macro_2.png", plot, width = 1500, height = 2500, units = "px", dpi = 300)

# 2023 ranking of GDP per capita
plot <- gdp_pop_data |> 
    group_by(Typology) |>
    summarize(across(starts_with("X")|starts_with("POPESTIMATE"), sum, na.rm = TRUE)) |> 
    mutate(
        across(
            .cols = all_of(paste0("X", 2001:2023)),
            .fns = ~ .x / get(str_replace(cur_column(), "X", "POPESTIMATE")),
            .names = "gdppc{.col}"
        )
    ) |> 
    select(Typology, starts_with("gdppc")) |>
    pivot_longer(
        cols = starts_with("gdppc"),
        names_to = "year",
        values_to = "gdp_per_capita"
    ) |>
    mutate(year = as.numeric(str_remove(year, "gdppcX"))) |> 
    select(year, gdp_per_capita, Typology) |> 
    filter(year == 2023) |>
    ggplot(aes(x = reorder(Typology, gdp_per_capita), y = gdp_per_capita)) +
    geom_col(fill = "steelblue") +
    coord_flip() +
    labs(
        title = "GDP Per Capita by Community Type (2023)",
        x = "Community Type",
        y = "GDP Per Capita (in 2017 dollars)",
        caption = "@mormon_metrics\ndata: US Bureau of Economic Analysis"
    ) +
    geom_text(
        aes(label = paste0("$",format(round(gdp_per_capita),big.mark=",", scientific = FALSE))),
        family = "Cairo",
        fontface = "bold",
        nudge_y = 14000
    ) +
    theme_minimal() +
    theme(
        axis.ticks = element_blank(),
        axis.title = element_text(size = 14),
        axis.text.y = element_text(size = 10),
        axis.text.x = element_text(size = 10, margin = margin(t = -5)),
        plot.background = element_rect(fill = "grey95"),
        panel.background = element_blank(),
        text = element_text(face = "bold", family = "Cairo"),
        plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
        plot.title.position = "plot",
        plot.subtitle.position = "plot",
        legend.title = element_blank(),
        legend.background = element_blank(),
        legend.box.background = element_blank(),
        legend.key = element_blank(),
        panel.grid = element_blank(),
        legend.position = "none",
        panel.grid.major.y = element_line(color = "lightgrey", linetype = 'solid'),
        plot.caption = element_text(size = 8,family = "Cairo")) +
    scale_y_continuous(labels = scales::label_number(scale = 1e-3, suffix = "K"), breaks = seq(0, 120000, 20000), limits = c(0, 120000))

ggsave("./images/14_macro_3.png", plot, width = 1500, height = 1500, units = "px", dpi = 300)

# Add in countries and do GDP per capita and compare countries to LDS Enclave

country_gdp_data <- read.csv("./data/macroecon/country_gdp.csv") |> 
    select(Country.Name, Country.Code, Indicator.Name, Indicator.Code, X2023)

country_gdppc <- read.csv("./data/macroecon/gdp-per-capita-worldbank.csv")

country_gdp_data |> 
    filter(!is.na(Country.Code)) |>
    select(Country.Name, X2023) |>
    bind_rows(data.frame(Country.Name = "LDS Enclave", X2023 = 242975363000*1.26) #converting to 2023 dollars
        ) |> 
    mutate(
        X2023 = log(X2023),
        Community = ifelse(Country.Name == "LDS Enclave", "LDS Enclave", "Other")
    ) |> 
    drop_na() |> 
    mutate(pctl = percent_rank(X2023)) |> 
    filter(Country.Name == "LDS Enclave")

plot <- country_gdp_data |> 
    filter(!is.na(Country.Code)) |>
    select(Country.Name, X2023) |>
    bind_rows(data.frame(Country.Name = "LDS Enclave", X2023 = 242975363000*1.26) #converting to 2023 dollars
        ) |> 
    mutate(
        X2023 = log(X2023),
        Community = ifelse(Country.Name == "LDS Enclave", "LDS Enclave", "Other")
    ) |> 
    drop_na() |> 
    ggplot(aes(x = reorder(Country.Name, X2023), y = X2023, fill = Community)) +
    geom_bar(stat = "identity") +
    coord_flip() +
    theme_minimal() +
    theme(
        axis.ticks = element_blank(),
        axis.text.y = element_text(size = 6),
        axis.title = element_text(size = 14),
        plot.background = element_rect(fill = "grey95"),
        panel.background = element_blank(),
        text = element_text(family = "Cairo"),
        plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
        plot.title.position = "plot",
        legend.title = element_blank(),
        legend.background = element_blank(),
        legend.box.background = element_blank(),
        legend.key = element_blank(),
        panel.grid = element_blank(),
        legend.position = "top",
        panel.grid.major.y = element_line(color = "lightgrey", linetype = 'solid'),
        plot.caption = element_text(size = 8,family = "Cairo")) +
        labs(
            title = "LDS Enclave GDP Higher Than 75% Of Countries...",
            x = "Country",
            y = "Log GDP (in 2023 dollars)",
            caption = "@mormon_metrics\ndata: World Bank, US Bureau of Economic Analysis"
        )

ggsave("./images/14_macro_4.png", plot, width = 1500, height = 4000, units = "px", dpi = 300)

names_to_filter_out <- c(
    "East Asia and Pacific (WB)",
    "Europe and Central Asia (WB)",
    "Latin America and Caribbean (WB)",
    "Middle East and North Africa (WB)",
    "North America (WB)",
    "South Asia (WB)",
    "Sub-Saharan Africa (WB)",
    "European Union (27)",
    "High-income countries",
    "Low-income countries",
    "Lower-middle-income countries",
    "Middle-income countries",
    "Upper-middle-income countries"
)

country_gdppc |> 
    filter(!is.na(Code)) |>
    rename(gdppc = GDP.per.capita..PPP..constant.2021.international...) |>
    filter(Year == 2023 & !Entity %in% names_to_filter_out) |>
    select(Entity, gdppc) |>
    drop_na() |> 
    bind_rows(data.frame(Entity = "LDS Enclave", gdppc = 62862*1.08) #converting to 2021 dollars
    ) |> 
    mutate(
        gdppc = log(gdppc),
        Community = ifelse(Entity == "LDS Enclave", "LDS Enclave", "Other")
    ) |> 
    drop_na() |> 
    mutate(pctl = percent_rank(gdppc)) |> 
    filter(Entity == "LDS Enclave")

plot <- country_gdppc |> 
    filter(!is.na(Code)) |>
    rename(gdppc = GDP.per.capita..PPP..constant.2021.international...) |>
    filter(Year == 2023 & !Entity %in% names_to_filter_out) |>
    select(Entity, gdppc) |>
    drop_na() |>
    bind_rows(data.frame(Entity = "LDS Enclave", gdppc = 62862*1.08) #converting to 2021 dollars
    ) |> 
    mutate(
        Community = ifelse(Entity == "LDS Enclave", "LDS Enclave", "Other"),
        gdppc = round(gdppc)
    ) |> 
    drop_na() |> 
    ggplot(aes(x = reorder(Entity, gdppc), y = gdppc, fill = Community)) +
    geom_bar(stat = "identity") +
    coord_flip() +
    theme_minimal() +
    theme(
        axis.ticks = element_blank(),
        axis.text.y = element_text(size = 6),
        axis.title = element_text(size = 14),
        plot.background = element_rect(fill = "grey95"),
        panel.background = element_blank(),
        text = element_text(family = "Cairo"),
        plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
        plot.title.position = "plot",
        legend.title = element_blank(),
        legend.background = element_blank(),
        legend.box.background = element_blank(),
        legend.key = element_blank(),
        panel.grid = element_blank(),
        legend.position = "top",
        panel.grid.major.y = element_line(color = "lightgrey", linetype = 'solid'),
        plot.caption = element_text(size = 8,family = "Cairo")) +
        labs(
            title = "... And GDP Per Capita Higher Than 90%",
            x = "Country",
            y = "GDP Per Capita (2021 international dollars)",
            caption = "@mormon_metrics\ndata: World Bank, US Bureau of Economic Analysis"
        ) +
    scale_y_continuous(labels = scales::dollar_format(scale = 1))

ggsave("./images/14_macro_5.png", plot, width = 1500, height = 4000, units = "px", dpi = 300)