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)