# Exploratory Data Analysis and Descriptive Statistics for Political Science Research in R Packages we will use:

``````library(tidyverse)      # of course
library(ggridges)       # density plots
library(GGally)         # correlation matrics
library(stargazer)      # tables
library(knitr)          # more tables stuff
library(kableExtra)     # more and more tables
library(ggstream)       # streamplots
library(bbplot)         # pretty themes
library(ggthemes)       # more pretty themes
library(ggside)         # stack plots side by side
library(forcats)        # reorder factor levels
``````

Before jumping into any inferentional statistical analysis, it is helpful for us to get to know our data. For me, that always means plotting and visualising the data and looking at the spread, the mean, distribution and outliers in the dataset.

Before we plot anything, a simple package that creates tables in the `stargazer` package. We can examine descriptive statistics of the variables in one table.

Click here to read this practically exhaustive cheat sheet for the stargazer package by Jake Russ. I refer to it at least once a week.

I want to summarise a few of the stats, so I write into the `summary.stat()` argument the number of observations, the mean, median and standard deviation.

The `kbl()` and `kable_classic()` will change the look of the table in R (or if you want to copy and paste the code into latex with the `type = "latex"` argument).

In HTML, they do not appear.

Choose the variables you want, put them into a` data.frame` and feed them into the `stargazer()` function

``````stargazer(my_df_summary,
covariate.labels = c("Corruption index",
"Civil society strength",
'Rule of Law score',
"Physical Integerity Score",
"GDP growth"),
summary.stat = c("n", "mean", "median", "sd"),
type = "html") %>%
kbl() %>%
kable_classic(full_width = F, html_font = "Times", font_size = 25)``````
 Statistic N Mean Median St. Dev. Corruption index 179 0.477 0.519 0.304 Civil society strength 179 0.670 0.805 0.287 Rule of Law score 173 7.451 7.000 4.745 Physical Integerity Score 179 0.696 0.807 0.284 GDP growth 163 0.019 0.020 0.032

Next, we can create a barchart to look at the different levels of variables across categories. We can look at the different regime types (from complete autocracy to liberal democracy) across the six geographical regions in 2018 with the `geom_bar()`.

``````my_df %>%
filter(year == 2018) %>%
ggplot() +
geom_bar(aes(as.factor(region),
fill = as.factor(regime)),
color = "white", size = 2.5) -> my_barplot``````

And we can add more theme changes

``````my_barplot + bbplot::bbc_style() +
theme(legend.key.size = unit(2.5, 'cm'),
legend.text = element_text(size = 15),
text = element_text(size = 15)) +
scale_fill_manual(values = c("#9a031e","#00a896","#e36414","#0f4c5c")) +
scale_color_manual(values = c("#9a031e","#00a896","#e36414","#0f4c5c")) ``````

This type of graph also tells us that Sub-Saharan Africa has the highest number of countries and the Middle East and North African (MENA) has the fewest countries.

However, if we want to look at each group and their absolute percentages, we change one line: we add `geom_bar(position = "fill")`. For example we can see more clearly that over 50% of Post-Soviet countries are democracies ( orange = electoral and blue = liberal democracy) as of 2018.

We can also check out the density plot of democracy levels (as a numeric level) across the six regions in 2018.

With these types of graphs, we can examine characteristics of the variables, such as whether there is a large spread or normal distribution of democracy across each region.

``````my_df %>%
filter(year == 2018) %>%
ggplot(aes(x = democracy_score, y = region, fill = regime)) +
geom_density_ridges(color = "white", size = 2, alpha = 0.9, scale = 2) -> my_density_plot``````

And change the graph theme:

``````my_density_plot + bbplot::bbc_style() +
theme(legend.key.size = unit(2.5, 'cm')) +
scale_fill_manual(values = c("#9a031e","#00a896","#e36414","#0f4c5c")) +
scale_color_manual(values = c("#9a031e","#00a896","#e36414","#0f4c5c")) ``````

Next, we can also check out Pearson’s correlations of some of the variables in our dataset. We can make these plots with the `GGally `package.

The `ggpairs()` argument shows a scatterplot, a density plot and correlation matrix.

``````my_df %>%
filter(year == 2018) %>%
select(regime,
corruption,
civ_soc,
rule_law,
physical,
gdp_growth) %>%
ggpairs(columns = 2:5,
ggplot2::aes(colour = as.factor(regime),
alpha = 0.9)) +
bbplot::bbc_style() +
scale_fill_manual(values = c("#9a031e","#00a896","#e36414","#0f4c5c")) +
scale_color_manual(values = c("#9a031e","#00a896","#e36414","#0f4c5c"))``````

We can use the ggside package to stack graphs together into one plot.

There are a few arguments to add when we choose where we want to place each graph.

For example, `geom_xsideboxplot(aes(y = freedom_house), orientation = "y")` places a boxplot for the three Freedom House democracy levels on the top of the graph, running across the x axis. If we wanted the boxplot along the y axis we would write `geom_ysideboxplot`(). We add `orientation = "y"` to indicate the direction of the boxplots.

Next we indiciate how big we want each graph to be in the panel with `theme(ggside.panel.scale = .5)` argument. This makes the scatterplot take up half and the boxplot the other half. If we write .3, the scatterplot takes up 70% and the boxplot takes up the remainning 30%. Last we indicade `scale_xsidey_discrete()` so the graph doesn’t think it is a continuous variable.

We add Darjeeling Limited color palette from the Wes Anderson movie.

``````my_df %>%
filter(year == 2018) %>%
filter(!is.na(fh_number)) %>%
mutate(freedom_house = ifelse(fh_number == 1, "Free",
ifelse(fh_number == 2, "Partly Free", "Not Free"))) %>%
mutate(freedom_house = forcats::fct_relevel(freedom_house, "Not Free", "Partly Free", "Free")) %>%
ggplot(aes(x = freedom_from_torture, y = corruption_level, colour = as.factor(freedom_house))) +
geom_point(size = 4.5, alpha = 0.9) +
geom_smooth(method = "lm", color ="#1d3557", alpha = 0.4) +
geom_xsideboxplot(aes(y = freedom_house), orientation = "y", size = 2) +
theme(ggside.panel.scale = .3) +
scale_xsidey_discrete() +
bbplot::bbc_style() +
facet_wrap(~region) +
scale_color_manual(values= wes_palette("Darjeeling1", n = 3))
``````

The next plot will look how variables change over time.

We can check out if there are changes in the volume and proportion of a variable across time with the `geom_stream(type = "ridge")` from the `ggstream `package.

In this instance, we will compare urban populations across regions from 1800s to today.

``````my_df %>%
group_by(region, year) %>%
summarise(mean_urbanization = mean(urban_population_percentage, na.rm = TRUE)) %>%
ggplot(aes(x = year, y = mean_urbanization, fill = region)) +
geom_stream(type = "ridge") -> my_streamplot``````

``````  my_streamplot + ggthemes::theme_pander() +
theme(
legend.title = element_blank(),
legend.position = "bottom",
legend.text = element_text(size = 25),
axis.text.x = element_text(size = 25),
axis.title.y = element_blank(),
axis.title.x = element_blank()) +
scale_fill_manual(values = c("#001219",
"#0a9396",
"#e9d8a6",
"#ee9b00",
"#ca6702",
"#ae2012")) ``````

We can also look at interquartile ranges and spread across variables.

We will look at the urbanization rate across the different regions. The variable is calculated as the ratio of urban population to total country population.

Before, we will create a hex color vector so we are not copying and pasting the colours too many times.

``````my_palette <- c("#1d3557",
"#0a9396",
"#e9d8a6",
"#ee9b00",
"#ca6702",
"#ae2012")``````

We use the `facet_wrap(~year) `so we can separate the three years and compare them.

``````my_df %>%
filter(year == 1980 | year == 1990 | year == 2000)  %>%
ggplot(mapping = aes(x = region,
y = urban_population_percentage,
fill = region)) +
geom_jitter(aes(color = region),
size = 3, alpha = 0.5, width = 0.15) +
geom_boxplot(alpha = 0.5) + facet_wrap(~year) +
scale_fill_manual(values = my_palette) +
scale_color_manual(values = my_palette) +
coord_flip() +
bbplot::bbc_style()``````

If we want to look more closely at one year and print out the country names for the countries that are outliers in the graph, we can run the following function and find the outliers int he dataset for the year 1990:

``````is_outlier <- function(x) {
return(x < quantile(x, 0.25) - 1.5 * IQR(x) | x > quantile(x, 0.75) + 1.5 * IQR(x))
}``````

We can then choose one year and create a binary variable with the function

``````my_df_90 <- my_df %>%
filter(year == 1990) %>%
filter(!is.na(urban_population_percentage))

my_df_90\$my_outliers <- is_outlier(my_df_90\$urban_population_percentage)``````

And we plot the graph:

``````my_df_90 %>%
ggplot(mapping = aes(x = region, y = urban_population_percentage, fill = region)) +
geom_jitter(aes(color = region), size = 3, alpha = 0.5, width = 0.15) +
geom_boxplot(alpha = 0.5) +
geom_text_repel(data = my_df_90[which(my_df_90\$my_outliers == TRUE),],
aes(label = country_name), size = 5) +
scale_fill_manual(values = my_palette) +
scale_color_manual(values = my_palette) +
coord_flip() +
bbplot::bbc_style()
``````

In the next blog post, we will look at t-tests, ANOVAs (and their non-parametric alternatives) to see if the difference in means / medians is statistically significant and meaningful for the underlying population.

# Create density plots with ggridges package in R Packages we will need:

``````library(tidyverse)
library(ggridges)
library(ggimage)  # to add png images
library(bbplot)   # for pretty graph themes``````

We will plot out the favourability opinion polls for the three main political parties in Ireland from 2016 to 2020. Data comes from Louwerse and Müller (2020)

Before we dive into the ggridges plotting, we have a little data cleaning to do. First, we extract the last four “characters” from the date string to create a year variable.

I took this quick function from a StackOverflow response:

``````substrRight <- function(x, n){
substr(x, nchar(x)-n+1, nchar(x))}

polls_csv\$year <- substrRight(polls_csv\$Date, 4)``````

Next, pivot the data from wide to long format.

More information of pivoting data with dplyr can be found here. I tend to check it at least once a month as the arguments refuse to stay in my head.

I only want to take the main parties in Ireland to compare in the plot.

``````polls <- polls_csv %>%
select(year, FG:SF) %>%
pivot_longer(!year, names_to = "party", values_to = "opinion_poll")``````

I went online and found the logos for the three main parties (sorry, Labour) and saved them in the working directory I have for my RStudio. That way I can call the file with the prefix “~/**.png” rather than find the exact location they are saved on the computer.

``````polls %>%
filter(party == "FF" | party == "FG" | party == "SF" ) %>%
mutate(image = ifelse(party=="FF","~/ff.png",
ifelse(party=="FG","~/fg.png", "~/sf.png"))) -> polls_three``````

Now we are ready to plot out the density plots for each party with the `geom_density_ridges()` function from the `ggridges` package.

We will add a few arguments into this function.

We add an `alpha = 0.8` to make each density plot a little transparent and we can see the plots behind.

The `scale = 2` argument pushes all three plots togheter so they are slightly overlapping. If `scale =1`, they would be totally separate and 3 would have them overlapping far more.

The `rel_min_height = 0.01` argument removes the trailing tails from the plots that are under 0.01 density. This is again for aesthetics and just makes the plot look slightly less busy for relatively normally distributed densities

The `geom_image` takes the images and we place them at the beginning of the x axis beside the labels for each party.

Last, we use the `bbplot` package BBC style ggplot theme, which I really like as it makes the overall graph look streamlined with large font defaults.

``````polls_three %>%
ggplot(aes(x = opinion_poll, y = as.factor(party))) +
geom_density_ridges(aes(fill = party),
alpha = 0.8,
scale = 2,
rel_min_height = 0.01) +
ggimage::geom_image(aes(y = party, x= 1, image = image), asp = 0.9, size = 0.12) +
facet_wrap(~year) +
bbplot::bbc_style() +
scale_fill_manual(values = c("#f2542d", "#edf6f9", "#0e9594")) +
theme(legend.position = "none") +
labs(title = "Favourability Polls for the Three Main Parties in Ireland", subtitle = "Data from Irish Polling Indicator (Louwerse & Müller, 2020)")``````

# Comparing mean values across OECD countries with ggplot Packages we will need:

``````library(tidyverse)
library(magrittr) # for pipes
library(ggrepel) # to stop overlapping labels
library(ggflags)
library(countrycode) # if you want create the ISO2C variable``````

I came across code for this graph by Tanya Shapiro on her github for #TidyTuesday.

Her graph compares Dr. Who actors and their average audience rating across their run as the Doctor on the show. So I have very liberally copied her code for my plot on OECD countries.

That is the beauty of TidyTuesday and the ability to be inspired and taught by other people’s code.

I originally was going to write a blog about how to download data from the `OECD` R package. However, my attempts to download the data leads to an unpleasant looking error and ends the donwload request.

I will try to work again on that blog in the future when the package is more established.

So, instead, I went to the OECD data website and just directly downloaded data on level of trust that citizens in each of the OECD countries feel about their governments.

Then I cleaned up the data in excel and used `countrycode()` to add ISO2 and country name data.

Click here to read more about the `countrycode()` package.

First I will only look at EU countries. I tried with all the countries from the OECD but it was quite crowded and hard to read.

I add region data from another dataset I have. This step is not necessary but I like to colour my graphs according to categories. This time I am choosing geographic regions.

``````my_df %<>%
mutate(region = case_when(
e_regiongeo == 1 ~ "Western",
e_regiongeo == 2  ~ "Northern",
e_regiongeo == 3  ~ "Southern",
e_regiongeo == 4  ~ "Eastern"))``````

To make the graph, we need two averages:

1. The overall average trust level for all countries (`avg_trust`) and
2. The average for each country across the years (`country_avg_trust`),
``````my_df %<>%
mutate(avg_trust = mean(trust, na.rm = TRUE)) %>%
group_by(country_name) %>%
mutate(country_avg_trust = mean(trust, na.rm = TRUE)) %>%
ungroup()``````

When we plot the graph, we need a few `geom` arguments.

Along the x axis we have all the countries, and `reorder` them from most trusting of their goverments to least trusting.

We will color the points with one of the four geographic regions.

We use` geom_jitter()` rather than `geom_point()` for the different yearly trust values to make the graph a little more interesting.

I also make the sizes scaled to the year in the `aes()` argument. Again, I did this more to look interesting, rather than to convey too much information about the different values for trust across each country. But smaller circles are earlier years and grow larger for each susequent year.

The `geom_hline()` plots a vertical line to indicate the average trust level for all countries.

We then use the `geom_segment()` to horizontally connect the country’s individual average (the `yend `argument) to the total average (the `y` arguement). We can then easily see which countries are above or below the total average. The `x` and `xend` argument, we supply the `country_name `variable twice.

Next we use the `geom_flag()`, which comes from the `ggflags` package. In order to use this package, we need the ISO 2 character code for each country in lower case!

Click here to read more about the `ggflags` package.

``````my_df %>%
ggplot(aes(x = reorder(country_name, trust_score), y = trust_score, color = as.factor(region))) +
geom_jitter(aes(color = as.factor(region), size = year), alpha = 0.7, width = 0.15) +
geom_hline(aes(yintercept = avg_trust), color = "white", size = 2)+
geom_segment(aes(x = country_name, xend = country_name, y = country_avg_trust, yend = avg_trust), size = 2, color = "white") +
ggflags::geom_flag(aes(x = country_name, y = country_avg_trust, country = iso2), size = 10) +
coord_flip() +
scale_color_manual(values = c("#9a031e","#fb8b24","#5f0f40","#0f4c5c")) -> my_plot``````

Last we change the aesthetics of the graph with all the theme arguments!

``````my_plot +
theme(panel.border = element_blank(),
legend.position = "right",
legend.title = element_blank(),
legend.text = element_text(size = 20),
legend.background = element_rect(fill = "#5e6472"),
axis.title = element_blank(),
axis.text = element_text(color = "white", size = 20),
text= element_text(size = 15, color = "white"),
panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank(),
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
legend.key = element_rect(fill = "#5e6472"),
plot.background = element_rect(fill = "#5e6472"),
panel.background = element_rect(fill = "#5e6472")) +
guides(colour = guide_legend(override.aes = list(size=10))) ``````

And that is the graph.

We can see that countries in southern Europe are less trusting of their governments than in other regions. Western countries seem to occupy the higher parts of the graph, with France being the least trusting of their government in the West.

There is a large variation in Northern countries. However, if we look at the countries, we can see that the Scandinavian countries are more trusting and the Baltic countries are among the least trusting. This shows they are more similar in their trust levels to other Post-Soviet countries.

Next we can look into see if there is a relationship between democracy scores and level of trust in the goverment with a `geom_point()` scatterplot

The `geom_smooth()` argument plots a linear regression OLS line, with a standard error bar around.

We want the labels for the country to not overlap so we use the `geom_label_repel()` from the `ggrepel` package. We don’t want an `a` in the legend, so we add `show.legend = FALSE` to the arguments

``````
my_df %>%
filter(!is.na(trust_score)) %>%
ggplot(aes(x = democracy_score, y = trust_score)) +
geom_smooth(method = "lm", color = "#a0001c", size = 3) +
geom_point(aes(color = as.factor(region)), size = 20, alpha = 0.6) +
geom_label_repel(aes(label = country_name, color = as.factor(region)), show.legend = FALSE, size = 5) +
scale_color_manual(values = c("#9a031e","#fb8b24","#5f0f40","#0f4c5c")) -> scatter_plot``````

And we change the theme and add labels to the plot.

``````scatter_plot + theme(panel.border = element_blank(),
legend.position = "bottom",
legend.title = element_blank(),
legend.text = element_text(size = 20),
legend.background = element_rect(fill = "#5e6472"),
text= element_text(size = 15, color = "white"),

legend.key = element_rect(fill = "#5e6472"),
plot.background = element_rect(fill = "#5e6472"),
panel.background = element_rect(fill = "#5e6472")) +
guides(colour = guide_legend(override.aes = list(size=10)))  +
labs(title = "Democracy and trust levels",
y = "Democracy score",
x = "Trust level of respondents",
caption="Data from OECD") ``````

We can filter out the two countries with low democracy and examining the consolidated democracies.

# Lollipop plots with ggplot2 in R Packages we will need:

``````library(tidyverse)
library(rvest)
library(ggflags)
library(countrycode)
library(ggpubr)``````

We will plot out a lollipop plot to compare EU countries on their level of income inequality, measured by the Gini coefficient.

A Gini coefficient of zero expresses perfect equality, where all values are the same (e.g. where everyone has the same income). A Gini coefficient of one (or 100%) expresses maximal inequality among values (e.g. for a large number of people where only one person has all the income or consumption and all others have none, the Gini coefficient will be nearly one).

To start, we will take data on the EU from Wikipedia. With rvest package, scrape the table about the EU countries from this Wikipedia page.

Click here to read more about the `rvest` pacakge

With the `gsub()` function, we can clean up the different variables with some regex. Namely delete the footnotes / square brackets and change the variable classes.

``````eu_site <- read_html("https://en.wikipedia.org/wiki/Member_state_of_the_European_Union")

eu_tables <- eu_site %>% html_table(header = TRUE, fill = TRUE)

eu_members <- eu_tables[]

eu_members %<>% janitor::clean_names()  %>%
filter(!is.na(accession))

eu_members\$iso3 <- countrycode::countrycode(eu_members\$geo, "country.name", "iso3c")

eu_members\$accession <- as.numeric(gsub("([0-9]+).*\$", "\\1",eu_members\$accession))

eu_members\$name_clean <- gsub("\\[.*?\\]", "", eu_members\$name)

eu_members\$gini_clean <- gsub("\\[.*?\\]", "", eu_members\$gini)``````

Next some data cleaning and grouping the year member groups into different decades. This indicates what year each country joined the EU. If we see clustering of colours on any particular end of the Gini scale, this may indicate that there is a relationship between the length of time that a country was part of the EU and their domestic income inequality level. Are the founding members of the EU more equal than the new countries? Or conversely are the newer countries that joined from former Soviet countries in the 2000s more equal. We can visualise this with the following mutations:

``````eu_members %>%
filter(name_clean != "Totals/Averages") %>%
mutate(gini_numeric = as.numeric(gini_clean)) %>%
mutate(accession_decades = ifelse(accession < 1960, "1957", ifelse(accession > 1960 & accession < 1990, "1960s-1980s", ifelse(accession == 1995, "1990s", ifelse(accession > 2003, "2000s", accession))))) -> eu_clean ``````

To create the lollipop plot, we will use the `geom_segment()` functions. This requires an `x` and `xend `argument as the country names (with the `fct_reorder()` function to make sure the countries print out in descending order) and a `y` and `yend` argument with the gini number.

All the countries in the EU have a gini score between mid 20s to mid 30s, so I will start the y axis at 20.

We can add the flag for each country when we turn the ISO2 character code to lower case and give it to the country argument.

Click here to read more about the `ggflags` package

``````eu_clean %>%
ggplot(aes(x= name_clean, y= gini_numeric, color = accession_decades)) +
geom_segment(aes(x = forcats::fct_reorder(name_clean, -gini_numeric),
xend = name_clean, y = 20, yend = gini_numeric, color = accession_decades), size = 4, alpha = 0.8) +
geom_point(aes(color = accession_decades), size= 10) +
geom_flag(aes(y = 20, x = name_clean, country = tolower(iso_3166_1_alpha_2)), size = 10) +
ggtitle("Gini Coefficients of the EU countries") -> eu_plot``````

Last we add various theme changes to alter the appearance of the graph

``````eu_plot +
coord_flip() +
ylim(20, 40) +
theme(panel.border = element_blank(),
legend.title = element_blank(),
axis.title = element_blank(),
axis.text = element_text(color = "white"),
text= element_text(size = 35, color = "white"),
legend.text = element_text(size = 20),
legend.key = element_rect(colour = "#001219", fill = "#001219"),
legend.key.width = unit(3, 'cm'),
legend.position = "bottom",
panel.grid.major.y = element_line(linetype="dashed"),
plot.background = element_rect(fill = "#001219"),
panel.background = element_rect(fill = "#001219"),
legend.background = element_rect(fill = "#001219") )``````

We can see there does not seem to be a clear pattern between the year a country joins the EU and their level of domestic income inequality, according to the Gini score.

Of course, the Gini coefficient is not a perfect measurement, so take it with a grain of salt.

Another option for the lolliplot plot comes from the `ggpubr `package. It does not take the familiar aesthetic arguments like you can do with `ggplot2 `but it is very quick and the defaults look good!

``````eu_clean %>%
ggdotchart( x = "name_clean", y = "gini_numeric",
sorting = "descending",
rotate = TRUE,
dot.size = 10,
y.text.col = TRUE,
ggtheme = theme_pubr()) +
ggtitle("Gini Coefficients of the EU countries") +
theme(panel.border = element_blank(),
legend.title = element_blank(),
axis.title = element_blank(),
axis.text = element_text(color = "white"),
text= element_text(size = 35, color = "white"),
legend.text = element_text(size = 20),
legend.key = element_rect(colour = "#001219", fill = "#001219"),
legend.key.width = unit(3, 'cm'),
legend.position = "bottom",
panel.grid.major.y = element_line(linetype="dashed"),
plot.background = element_rect(fill = "#001219"),
panel.background = element_rect(fill = "#001219"),
legend.background = element_rect(fill = "#001219") )``````