This R notebook explores the “Tax Revenues as Share of GDP” dataset from Our World in Data.

UNU-WIDER Government Revenue Dataset (2023) – with major processing by Our World in Data. “Taxes including social contributions – UNU-WIDER” [dataset]. UNU-WIDER, “Government Revenue Dataset (GRD) 2023” [original data]. Source: UNU-WIDER Government Revenue Dataset (2023) – with major processing by Our World In Data

library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.6
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ ggplot2   4.0.1     ✔ tibble    3.3.0
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.2
## ✔ purrr     1.2.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(sf)
## Linking to GEOS 3.13.1, GDAL 3.11.4, PROJ 9.7.0; sf_use_s2() is TRUE
library(biscale)
## Warning: package 'biscale' was built under R version 4.5.3
library(cowplot)
## Warning: package 'cowplot' was built under R version 4.5.3
## 
## Attaching package: 'cowplot'
## 
## The following object is masked from 'package:lubridate':
## 
##     stamp
tax.data <- read.csv(file = "tax-revenues-as-a-share-of-gdp-unu-wider.csv")
colnames(tax.data) <- c("Country", "Code", "Year", 
                        "Tax.Revenues.Shares.GDP", "Region")
head(tax.data, 5)
##       Country Code Year Tax.Revenues.Shares.GDP Region
## 1 Afghanistan  AFG 2003                2.512631   Asia
## 2 Afghanistan  AFG 2004                4.076170   Asia
## 3 Afghanistan  AFG 2005                4.668273   Asia
## 4 Afghanistan  AFG 2006                7.115892   Asia
## 5 Afghanistan  AFG 2007                6.061553   Asia
summary(tax.data$Tax.Revenues.Shares.GDP)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##  0.08581 11.24129 17.38397 19.84494 27.54285 60.94643
ggplot(
    data = tax.data
) +
    geom_histogram(
        mapping = aes(x= Tax.Revenues.Shares.GDP),
        bins = 30, 
        col = "navy", 
        fill = "royalblue"
    ) +
    labs(
        title = "Histogram of Tax Revenues",
        x = "% of GDP",
        y = "Frequency"
    ) +
    theme_bw() +
    theme(
        plot.title = element_text(hjust = 0.5)
    )

length(unique(tax.data$Code))
## [1] 191
# Count of records per year

ggplot(
    data = tax.data
) +
    geom_bar(
        stat = "count",
        mapping = aes(x = Year),
        col = "navy",
        fill = "royalblue"
    ) +
    labs(
        title = "Count of Countries by Year",
        x = "Year",
        y = "Countries"
    ) +
    theme_bw() +
    theme(
        plot.title = element_text(hjust = 0.5)
    )

ggplot(
    data = tax.data
) +
    geom_point(
        mapping = aes(x = Year, y = Tax.Revenues.Shares.GDP), 
        col = "royalblue"
    ) +
    geom_smooth(
        mapping = aes(x = Year, y = Tax.Revenues.Shares.GDP),
        col = "coral"
    ) +
    labs(
        title = "Tax Revenues by Years",
        x = "Year",
        y = "% of GDP"
    ) +
    theme_bw() +
    theme(
        plot.title = element_text(hjust = 0.5)
    )
## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'

# Find distribution of Tax Revenues % of GDP across years

fltr <- (tax.data$Year >= 2000 & tax.data$Year <= 2020)

ggplot(
    data = tax.data[fltr, ]
) +
    geom_boxplot(
        mapping = aes(x = as.factor(Year), y = Tax.Revenues.Shares.GDP),
        col = "navy",
        fill = "royalblue",
        staplewidth = 0.5,
        outlier.color = "red",
        median.color = "coral"
    ) +
    labs(
        title = "Distribution of Tax Revenues by Year",
        x = "Year",
        y = "% of GDP"
    ) +
    theme_bw() +
    theme(
        plot.title = element_text(hjust = 0.5),
        axis.text.x = element_text(angle = 90, vjust = 0.5)
    )

# Calculate distribution of Tax Revenues Shares as GDP by region

ggplot(
    data = tax.data
) +
    geom_boxplot(
        mapping = aes(x = as.factor(Region), y = Tax.Revenues.Shares.GDP),
        fill = "royalblue",
        col = "navy",
        outlier.color = "red",
        median.color = "coral",
        staplewidth = 0.5
    ) +
    labs(
        title = "Tax Revenues Shares as GDP Distribution by Region",
        x = "Region",
        y = "% of GDP"
    ) +
    theme_bw() +
    theme(
        plot.title = element_text(hjust = 0.5)
    )

fltr <- (tax.data$Year >= 2000 & tax.data$Year <= 2020)

tx.medians <- tax.data[fltr, ] %>%
    group_by(
        Code
    ) %>%
    summarise(
        Tax.Revenues.Median = median(Tax.Revenues.Shares.GDP)
    ) %>%
    arrange(
        desc(Tax.Revenues.Median)
    ) %>%
    head(50)


fltr <- (tax.data$Year >= 2000 & tax.data$Year <= 2020 &
             tax.data$Code %in% tx.medians$Code)

ggplot(
    data = tax.data[fltr, ]
) +
    geom_boxplot(
        mapping = aes(x = as.factor(Code), y = Tax.Revenues.Shares.GDP),
        fill = "royalblue",
        col = "navy",
        median.color = "coral",
        outlier.color = "red",
        staplewidth = 0.5
    ) +
    labs(
        title = "Distribution of Tax Revenues by Country",
        subtitle = "Top 50 Median Values (2000 - 2020)",
        x = "Country",
        y = "% of GDP"
    ) +
    theme_bw() +
    theme(
        plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5),
        axis.text.x = element_text(angle = 90, vjust = 0.5)
    )

fltr <- (tax.data$Year >= 2000 & tax.data$Year <= 2020)

map.data <- tax.data[fltr, ] %>%
    group_by(Code) %>%
    summarise(
        Tax.Median = median(Tax.Revenues.Shares.GDP)
    )

country.data <- st_read(dsn = "../world-ash-ms.geojson") %>%
    .[.$iso_a3 != "ATA", c("iso_a3", "admin", "geometry")] %>%
    mutate(
        iso_a3 = ifelse(admin == "France", "FRA", iso_a3)
    )
## Reading layer `world-ash-ms' from data source 
##   `C:\Workspace\Data Science Programming\world-ash-ms.geojson' 
##   using driver `GeoJSON'
## Simple feature collection with 242 features and 169 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -180 ymin: -89.99893 xmax: 180 ymax: 83.59961
## Geodetic CRS:  WGS 84
map.data <- merge(country.data, map.data,
                  by.x = "iso_a3", by.y = "Code",
                  all.x = TRUE)
ggplot(
    data = map.data
) +
    geom_sf(
        mapping = aes(fill = Tax.Median)
    ) +
    scale_fill_viridis_c(option = "plasma") +
    labs(
        title = "Median Tax Revenues as Shares of GDP by Country",
        subtitle = "For Years 2000 to 2020",
        fill = "% of GDP"
    ) +
    theme_bw() +
    theme(
        plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5),
        axis.text = element_blank()
    )

# Countries with the highest range of distribution

calc_iqr <- function(cntry) {
    val = 0.0
    
    fltr <- (tax.data$Year >= 2000 & tax.data$Year <= 2020 
                   & tax.data$Code == cntry)
    
    x <- tax.data[fltr, ]
    s <- summary(x$Tax.Revenues.Shares.GDP)
    val <- (s[3] - s[1])
        
    returnValue(val)
}

fltr <- (tax.data$Year >= 2000 & tax.data$Year <= 2020)

map.data <- tax.data[fltr, ] %>%
    group_by(Code) %>%
    summarise(Country = first(Country))

map.data$Tax.IQR <- sapply(map.data$Code, FUN = function(x) {calc_iqr(x)})

map.data <- merge(country.data, map.data,
                  by.x = "iso_a3", by.y = "Code",
                  all.x = TRUE)

ggplot(
    data = map.data
) +
    geom_sf(
        mapping = aes(fill = Tax.IQR)
    ) +
    scale_fill_viridis_c(option = "plasma") +
    labs(
        title = "Inter-Quartile Range of Tax Revenues as Shares of GDP by Country",
        subtitle = "For Years 2000 to 2020",
        fill = "IQR"
    ) +
    theme_bw() +
    theme(
        plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5),
        axis.text = element_blank()
    )


Generate a multi-variate map of Median Tax Revenues by IQR Tax Revenues.

Generate the same multi-variate map using ggplot.

fltr <- (tax.data$Year >= 2000 & tax.data$Year <= 2020)

map.data <- tax.data[fltr, ] %>%
    group_by(Code) %>%
    summarise(
        Country = first(Country),
        Tax.Median = median(Tax.Revenues.Shares.GDP)
    )

map.data$Tax.IQR <- sapply(map.data$Code, FUN = function(x) {calc_iqr(x)})

map.data$Tax.IQR <- ntile(map.data$Tax.IQR, n = 3)
map.data$Tax.Median <- ntile(map.data$Tax.Median, n = 3)

map.data <- merge(country.data, map.data,
                  by.x = "iso_a3", by.y = "Code",
                  all.x = TRUE)

map.data <- bi_class(
    map.data, 
    x = Tax.Median, y = Tax.IQR, 
    style = "quantile", dim = 3
)
## Warning in classInt::classIntervals(.data[[var]], n = dim, style = style): var
## has missing values, omitted in finding classes
## Warning in classInt::classIntervals(.data[[var]], n = dim, style = style): n
## same as number of different finite values\neach different finite value is a
## separate class
## Warning in classInt::classIntervals(.data[[var]], n = dim, style = style): var
## has missing values, omitted in finding classes
## Warning in classInt::classIntervals(.data[[var]], n = dim, style = style): n
## same as number of different finite values\neach different finite value is a
## separate class
# Create the map
map <- ggplot() +
    geom_sf(
        data = map.data, 
        mapping = aes(fill = bi_class), 
        color = "black", 
        size = 0.1, 
        show.legend = FALSE
    ) +
    bi_scale_fill(
        pal = "DkBlue",
        dim = 3
    ) +
    labs(
        title = "Median vs IQR Tax Revenues as Shares of GDP", 
        subtitle = "For Years 2000 to 2020"
    ) +
    theme_bw() +
    theme(
        plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5),
        axis.text = element_blank()
    )

# Create the legend
legend <- bi_legend(
    pal = "DkBlue", 
    dim = 3, 
    xlab = "Median", 
    ylab = "IQR", 
    size = 8
)

# Combine map with legend
cowplot::ggdraw() +
    cowplot::draw_plot(map, 0, 0, 1, 1) +
    cowplot::draw_plot(legend, 0.0, 0.15, 0.2, 0.2)