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(tmap)
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, 50)
## 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
## 6 Afghanistan AFG 2008 6.174382 Asia
## 7 Afghanistan AFG 2009 8.673145 Asia
## 8 Afghanistan AFG 2010 9.559737 Asia
## 9 Afghanistan AFG 2011 9.359386 Asia
## 10 Afghanistan AFG 2012 7.950485 Asia
## 11 Afghanistan AFG 2013 7.548227 Asia
## 12 Afghanistan AFG 2014 7.141338 Asia
## 13 Afghanistan AFG 2015 7.874765 Asia
## 14 Afghanistan AFG 2016 9.837913 Asia
## 15 Afghanistan AFG 2017 10.187634 Asia
## 16 Afghanistan AFG 2018 9.173828 Asia
## 17 Afghanistan AFG 2019 8.362843 Asia
## 18 Afghanistan AFG 2020 6.760211 Asia
## 19 Albania ALB 1992 13.806224 Europe
## 20 Albania ALB 1993 15.244532 Europe
## 21 Albania ALB 1994 16.086613 Europe
## 22 Albania ALB 1995 14.837486 Europe
## 23 Albania ALB 1996 12.825735 Europe
## 24 Albania ALB 1997 13.784090 Europe
## 25 Albania ALB 1998 18.750260 Europe
## 26 Albania ALB 1999 18.830282 Europe
## 27 Albania ALB 2000 20.769794 Europe
## 28 Albania ALB 2001 20.617838 Europe
## 29 Albania ALB 2002 21.121912 Europe
## 30 Albania ALB 2003 21.441029 Europe
## 31 Albania ALB 2004 22.503714 Europe
## 32 Albania ALB 2005 22.486690 Europe
## 33 Albania ALB 2006 23.188566 Europe
## 34 Albania ALB 2007 23.347008 Europe
## 35 Albania ALB 2008 24.228289 Europe
## 36 Albania ALB 2009 24.003004 Europe
## 37 Albania ALB 2010 23.277954 Europe
## 38 Albania ALB 2011 22.882267 Europe
## 39 Albania ALB 2012 22.238976 Europe
## 40 Albania ALB 2013 21.395224 Europe
## 41 Albania ALB 2014 23.933598 Europe
## 42 Albania ALB 2015 23.985617 Europe
## 43 Albania ALB 2016 23.537329 Europe
## 44 Albania ALB 2018 25.149136 Europe
## 45 Albania ALB 2019 24.942486 Europe
## 46 Albania ALB 2020 23.912548 Europe
## 47 Albania ALB 2021 24.741135 Europe
## 48 Angola AGO 1993 38.236233 Africa
## 49 Angola AGO 1994 41.451614 Africa
## 50 Angola AGO 1995 23.964706 Africa
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.
colors <- c(
"1-1" = "#FFF0F5",
"1-2" = "#f78fb6",
"1-3" = "#f73593",
"2-1" = "#a5e8cd",
"2-2" = "#a58fb6",
"2-3" = "#a53593",
"3-1" = "#40dba7",
"3-2" = "#408fa7",
"3-3" = "#403593"
)
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)
tm_shape(map.data) +
tm_polygons(
fill = tm_vars(
c("Tax.Median", "Tax.IQR"),
multivariate = TRUE
),
fill.scale = tm_scale_bivariate(
scale1 = tm_scale_categorical(),
scale2 = tm_scale_categorical(),
values = colors
),
fill.legend = tm_legend_bivariate(
xlab = "IQR",
ylab = "Median",
show = TRUE,
position = tm_pos_in("left", "bottom")
)
) +
tm_title(
"Median vs IQR Tax Revenues as Shares of GDP\nFor Years 2000 to 2020",
size = 1,
frame = FALSE,
position = tm_pos_out("center", "top", align.h = "center")
)
## [tip] Consider a suitable map projection, e.g. by adding `+ tm_crs("auto")`.
## Labels abbreviated by the first letters, e.g.: "1" => "1"
##
## This message is displayed once per session.
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)