Example breakdown of PSA data on regions of the Philippines by GDP.

philippine-regional-gdp.R

if (!require("pacman")) install.packages("pacman")
pacman::p_load(
    "dplyr", "ggplot2", "readr", "reticulate",
    "sf", "showtext", "stringr"
)

# Utility function to save ggplots
save_image <- function(filename, plot, width = 7, height = 7) {
    ggsave(
        sprintf("%s.png", filename),
        plot,
        width = width, height = height,
        units = "in", dpi = 96
    ) 
}
# ggplot theme and text settings
plot_settings <- list(
    labs(
        caption = "Data: Philippine Statistics Authority | github.com/emordonez"
    ),
    theme_minimal(),
    theme(
        legend.position = "none",
        plot.margin = margin(10, 10, 10, 10, "pt"),
        text = element_text(family = "Open Sans"),
        plot.title = element_text(face = "bold"),
        plot.subtitle = element_text(face = "plain"),
        plot.caption = element_text(face = "italic"),
        axis.title.x = element_text(margin = margin(7, 0, 0, 0, "pt")),
        axis.title.y = element_text(margin = margin(0, 7, 0, 0, "pt"))
    )
)
font_add_google("Open Sans", "Open Sans")
showtext_auto(enable = TRUE)

# Read data
source_python("regional_data.py")
df <- format_data("./regional-gdp.csv", "./regional-population.csv")

# Compute regional shares of national aggregates
ph <- df %>%
    group_by(year) %>%
    dplyr::summarize(
        national_gdp = sum(gdp_constant_php), national_pop = sum(population)
    ) %>%
    mutate(national_gdp_pc = national_gdp / national_pop)

df <- merge(df, ph, by = "year") %>%
    mutate(
        pop_share = population / national_pop,
        gdp_pc_ratio = gdp_per_capita / national_gdp_pc,
        region = case_when(
            # Numbered regions have region name in parentheses
            startsWith(region, "Region") ~ region %>%
                str_extract(., "\\([^()]+\\)") %>%
                substring(2, nchar(.) - 1) %>%
                str_to_title(),
            region == "National Capital Region (NCR)" ~ "Metro Manila",
            region == "Cordillera Administrative Region (CAR)" ~ "Cordillera Region",
            region == "MIMAROPA Region" ~ "Mimaropa",
            region == "Bangsamoro Autonomous Region in Muslim Mindanao (BARMM)" ~ "Bangsamoro"
        )
    )


#'
#' PLOT 1: Philippine regions ranked by real GDP per capita
#' 
region_rankings <- df %>%
    # Rankings in 2000, 2010, and 2020
    filter(year == 2000 | year == 2010 | year == 2020) %>%
    select(year, region, gdp_per_capita) %>%
    group_by(year) %>%
    mutate(ranking = dense_rank(desc(gdp_per_capita))) %>%
    ungroup() %>%
    # Net movement up or down the rankings from 2000 to 2020
    group_by(region) %>%
    mutate(trend = case_when(
        ranking[year == 2020] < ranking[year == 2000] ~ "increase",
        ranking[year == 2020] > ranking[year == 2000] ~ "decrease",
        TRUE ~ "same"
    )) %>%
    ungroup()

rankings_2000 <- region_rankings %>% filter(year == 2000)
viz_region_rankings <- region_rankings %>%
    ggplot(aes(x = year, y = ranking, group = region, color = trend)) +
    geom_line() +
    geom_point() +
    scale_x_continuous(breaks = c(2000, 2010, 2020)) +
    scale_y_reverse(
        labels = rankings_2000$region,
        breaks = rankings_2000$ranking,
        sec.axis = dup_axis(labels = rankings_2000$ranking)
    ) +
    scale_color_manual(values = c("red", "blue", "black")) +
    labs(
        title = "Philippine regions ranked by real GDP per capita",
        subtitle = paste(
            "The fastest growing regions over the past 20 years have been",
            "\nNorthern Mindanao, Mimaropa, Western Visayas, and Caraga."
        ),
        x = NULL,
        y = NULL
    ) +
    plot_settings
save_image("1-region-rankings", viz_region_rankings)


#'
#' PLOT 2: Regional GDP per capita relative to the national average
#' 

# Using a 2011 PhilGIS shapefile
# 2011 region boundaries are coterminous with post-2018 region boundaries,
#   but some regions are named differently
download.file(
    "https://drive.google.com/uc?export=download&id=0B-JUh-O7Rg6YZmRwVEF5aWFsWVU",
    destfile = "Regions.7z"
)
system("7z e -oRegions Regions.7z")
regions_shp <- read_sf("./Regions/", "Regions")

# Match region names
regions_shp <- regions_shp %>%
    mutate(
        REGION = case_when(
            # Numbered regions end in " (Region ___)"
            grepl("(Region", REGION, fixed = TRUE) ~
                gsub(" \\([^)]*\\)", "", REGION, perl = TRUE) %>%
                str_to_title(),
            REGION == "Metropolitan Manila" ~ "Metro Manila",
            REGION == "Cordillera Administrative Region (CAR)" ~ "Cordillera Region",
            REGION == "Autonomous Region of Muslim Mindanao (ARMM)" ~ "Bangsamoro"
        )
    )

gdp_pc_ratios_2020 <- df %>%
    filter(year == 2020) %>%
    select(region, gdp_pc_ratio, national_gdp_pc) %>%
    merge(y = regions_shp, by.x = "region", by.y = "REGION")
gdp_pc_ratios_2020$ratio_bins <- cut(
    gdp_pc_ratios_2020$gdp_pc_ratio,
    breaks = c(
        0, 0.5, 0.6, 0.7, 0.8, 0.9, 1, Inf
    ),
    labels = c(
        "< 50", "50-60", "60-70", "70-80", "80-90", "90-100", "> 100"
    )
)

ph_gdp_pc_2020 <- tail(ph$national_gdp_pc, n = 1)
manila_ratio_2020 <- tail(
    filter(gdp_pc_ratios_2020, region == "Metro Manila")$gdp_pc_ratio,n = 1
)
northern_mindanao_ratio_2020 <- tail(
    filter(gdp_pc_ratios_2020, region == "Northern Mindanao")$gdp_pc_ratio,
    n = 1
)

viz_gdp_pc_ratios_2020 <- ggplot(data = gdp_pc_ratios_2020, aes(geometry = geometry)) +
    geom_sf(aes(fill = ratio_bins), color = "transparent", show.legend = "polygon") +
    geom_sf(fill = "transparent", color = "black", size = 0.25) +
    annotate(
        "text",
        x = 118.5,
        y = 14.5,
        label = sprintf(
            "Metro Manila: %0.0f%%",
            manila_ratio_2020 * 100
        ),
        size = 4
    ) +
    annotate(
        "text",
        x = 121,
        y = 8.75,
        label = sprintf(
            "Northern Mindanao: %0.0f%%",
            northern_mindanao_ratio_2020 * 100
        ),
        size = 4
    ) +
    scale_fill_manual(
        # There are no values in what should be the fourth bin
        values = c(
            "#003f5c", "#374c80", "#7a5195", #"#bc5090",
            "#ef5675", "#ff764a", "#ffa600"
        ),
        name = "Percent of national average",
        guide = guide_legend(
            keyheight = unit(8, units = "mm"),
            keywidth = unit(8, units = "mm"),
            title.position = "top",
            label.position = "right",
            nrow = 8
        )
    ) +    
    labs(
        title = "Regional GDP per capita relative to the national average",
        subtitle = paste(
            "The most productive regions per capita are those surrounding",
            "Manila and \nDavao. National GDP per capita in 2020 was",
            sprintf("%0.2f", ph_gdp_pc_2020),
            "thousand pesos at \nconstant 2018 prices."
        ),
        x = NULL,
        y = NULL
    ) +
    plot_settings +
    theme(
        legend.position = c(0, 0.975),
        legend.justification = c(0, 1),
        axis.text.x = element_blank(),
        axis.text.y = element_blank()
    )
save_image("2-regional-gdp-ratios-2020", viz_gdp_pc_ratios_2020, height = 10)


#'
#' PLOT 3: Real GDP per capita vs. share of national population
#' 
ph_pop_shares_2020 <- df %>%
    filter(year == 2020) %>%
    select(region, gdp_per_capita, pop_share, national_gdp, national_pop) %>%
    arrange(gdp_per_capita) %>%
    # Cumulative share of population with hack to get the plot width right
    mutate(
        pop_share = pop_share * 100,
        cum_share = cumsum(pop_share),
        loc = cum_share - pop_share / 2
    )

# IMF's 2018 average exchange rate for PHP to USD
xtr <- 52.661
# 2020 national aggregates
ph_gdp_2020 <- tail(ph$national_gdp, n = 1)
ph_pop_2020 <- tail(ph$national_pop, n = 1)
ph_gdp_pc_2020 <- tail(ph$national_gdp_pc, n = 1)

viz_ph_pop_shares_2020 <- ph_pop_shares_2020 %>%
    ggplot(aes(x = loc, y = gdp_per_capita, width = pop_share)) +
    geom_bar(
        aes(fill = region), alpha = 0.75,
        stat = "identity", position = "identity"
    ) +
    geom_hline(yintercept = ph_gdp_pc_2020) +
    scale_x_continuous(
        name = "Cumulative share of population (%)",
        breaks = seq(0, 100, 10),
        labels = as.character(seq(0, 100, 10)),
        sec.axis = sec_axis(
            ~ . * 1,
            breaks = ph_pop_shares_2020$loc,
            labels = ph_pop_shares_2020$region,            
        ),
        expand = c(0.01, 0.25)
    ) +
    scale_y_continuous(
        expand = c(0.01, 0.25)
    ) +
    annotate(
        "text",
        x = 50,
        y = 275,
        label = sprintf("National GDP per capita, %0.0f", ph_gdp_pc_2020),
        size = 4
    ) +
    labs(
        title = "Real GDP per capita vs. share of national population",
        subtitle = paste(
            "The population of the Philippines in 2020 was",
            sprintf("%0.2f million,", ph_pop_2020 / 1e6),
            "and total GDP at \nconstant 2018 prices was",
            sprintf("%0.2f trillion PHP", ph_gdp_2020 / 1e9),
            sprintf("(%0.0f billion USD, 1 USD = %0.2f PHP).",
                ph_gdp_2020 / 1e6 / xtr, xtr
            )
        ),
        x = NULL,
        y = "Real GDP per capita, thousand PHP (2018)"
    ) +
    plot_settings +
    theme(
        plot.caption = element_text(hjust = 0)
    ) +
    coord_flip()
save_image("3-rgdp-vs-population", viz_ph_pop_shares_2020)

regional-gdp.csv

regional-population.csv

regional_data.py

import pandas as pd
pd.options.mode.chained_assignment = None


def format_data(gdp_csv, population_csv):
    def read_data(filename, indicator):
        value_name = ''
        if indicator == 'gdp':
            value_name = 'gdp_constant_php'
        elif indicator == 'pop':
            value_name = 'population'
        else:
            return

        df = pd.read_csv(filename)
        df = pd.melt(
            df,
            id_vars='Region',
            value_vars=df.columns[1:],
            var_name='year',
            value_name=value_name
        )
        df.rename(columns={'Region': 'region'}, inplace=True)

        if indicator == 'gdp':
            df['year'] = df.year.apply(lambda x: x[-4:])
        df['year'] = pd.to_numeric(df.year, downcast='integer')

        return df


    gdp = read_data(gdp_csv, 'gdp')
    population = read_data(population_csv, 'pop')
    
    df = gdp.merge(population, on=['region', 'year'])
    df['region'] = df.region.apply(lambda x: x[2:])
    df['gdp_per_capita'] = df.gdp_constant_php / df.population

    cols = ['gdp_constant_php', 'population', 'gdp_per_capita']
    growth_cols = ['gdp_growth', 'pop_growth', 'gdp_pc_growth']
    df_list = []
    for region in set(df.region):
        region_df = df.loc[df.region == region]
        for col, growth_col in zip(cols, growth_cols):
            region_df[growth_col] = region_df[col].pct_change().fillna(0)
        df_list.append(region_df)
    
    return pd.concat(df_list, ignore_index=True)