Case Study: Exploratory Data Analysis in R - United Nations voting dataset

Data cleaning and summarizing with dplyr

We’ll learn how to clean and filter the United Nations voting dataset using the dplyr package, and how to summarize it into smaller, interpretable units.

votes <- readRDS("votes.rds")

The United Nations Voting Dataset

Filtering rows

The vote column in the dataset has a number that represents that country’s vote:

• 1 = Yes
• 2 = Abstain
• 3 = No
• 8 = Not present
• 9 = Not a member

One step of data cleaning is removing observations (rows) that you’re not interested in. In this case, you want to remove “Not present” and “Not a member”.

# Load the dplyr package
library(dplyr)

# Print the votes dataset
head(votes)

# Filter for votes that are "yes", "abstain", or "no"
votes_filtered <- votes %>% 
  filter(vote <= 3)
head(votes_filtered)

Adding a year column

The next step of data cleaning is manipulating your variables (columns) to make them more informative.

In this case, you have a session column that is hard to interpret intuitively. But since the UN started voting in 1946, and holds one session per year, you can get the year of a UN resolution by adding 1945 to the session number.

# Add another %>% step to add a year column
votes_filtered <- votes %>%
  filter(vote <= 3) %>% 
  mutate(year = session + 1945)
head(votes_filtered)

Adding a country column

The country codes in the ccode column are what’s called Correlates of War codes. This isn’t ideal for an analysis, since you’d like to work with recognizable country names.

You can use the countrycode package to translate. For example:

library(countrycode)

# Translate the country code 2
countrycode(2, "cown", "country.name")

[1] "United States"

# Translate multiple country codes
countrycode(c(2, 20, 40), "cown", "country.name")

[1] "United States" "Canada"        "Cuba"         

# Convert country code 100
countrycode(100, "cown", "country.name")

[1] "Colombia"

# Add a country column within the mutate: votes_processed
votes_processed <- votes %>%
  filter(vote <= 3) %>%
  mutate(year = session + 1945,
  country = countrycode(ccode, "cown", "country.name")) 

Some values were not matched unambiguously: 260

Grouping and summarizing

Summarizing the full dataset

In this analysis, we’re going to focus on “% of votes that are yes” as a metric for the “agreeableness” of countries.

We’ll start by finding this summary for the entire dataset: the fraction of all votes in their history that were “yes”. Note that within our call to summarize(), we can use n() to find the total number of votes and mean(vote == 1) to find the fraction of “yes” votes.

# Print votes_processed
head(votes_processed)

# Find total and fraction of "yes" votes
votes_processed %>%
  summarize(total = n(),
  percent_yes = mean(vote == 1))

Summarizing by year

The summarize() function is especially useful because it can be used within groups.

For example, we might like to know how much the average “agreeableness” of countries changed from year to year. To examine this, you can use group_by() to perform your summary not for the entire dataset, but within each year.

# Change this code to summarize by year
votes_processed %>%
  group_by(year) %>%
  summarize(total = n(),
            percent_yes = mean(vote == 1))

The group_by() function must go before your call to summarize() when you’re trying to perform your summary within groups.

Summarizing by country

We could instead summarize() within each country, which would let you compare voting patterns between countries.

# Summarize by country: by_country
by_country = votes_processed %>%
  group_by(country) %>%
  summarize(total = n(),
            percent_yes = mean(vote == 1))
head(by_country)

Sorting and filtering summarized data

Sorting by percentage of “yes” votes

Now that we ’ve summarized the dataset by country, we can start examining it and answering interesting questions.

For example, we might be especially interested in the countries that voted “yes” least often, or the ones that voted “yes” most often.

# You have the votes summarized by country
by_country <- votes_processed %>%
  group_by(country) %>%
  summarize(total = n(),
            percent_yes = mean(vote == 1))

# Print the by_country dataset
head(by_country)

# Sort in ascending order of percent_yes
head(by_country %>%
  arrange(percent_yes))

# Now sort in descending order
head(by_country %>%
  arrange(desc(percent_yes)))

Filtering summarized output

In the last exercise, we may have noticed that the country that voted least frequently, Zanzibar, had only 2 votes in the entire dataset. We certainly can’t make any substantial conclusions based on that data!

Typically in a progressive analysis, when we find that a few of your observations have very little data while others have plenty, you set some threshold to filter them out.

# Filter out countries with fewer than 100 votes
head(by_country %>%
  arrange(percent_yes) %>%
  filter(total >= 100))

Data visualization with ggplot2

Plotting a data over time

Line plotter

We’ll now use the ggplot2 package to turn our results into a visualization of the percentage of “yes” votes over time.

# Define by_year
by_year <- votes_processed %>%
  group_by(year) %>%
  summarize(total = n(),
            percent_yes = mean(vote == 1))

# Load the ggplot2 package
library(ggplot2)

# Create line plot
ggplot(by_year, aes(x = year, y = percent_yes)) +
  geom_line()

Scatter plotter

A line plot is one way to display this data. Weu could also choose to display it as a scatter plot, with each year represented as a single point. This requires changing the layer (i.e. geom_line() to geom_point()).

We can also add additional layers to your graph, such as a smoothing curve with geom_smooth().

# Change to scatter plot and add smoothing curve
ggplot(by_year, aes(year, percent_yes)) +
  geom_point() +
  geom_smooth()

Visualizing by country

Summarizing by year and country

We’re more interested in trends of voting within specific countries than we are in the overall trend. So instead of summarizing just by year, summarize by both year and country, constructing a dataset that shows what fraction of the time each country votes “yes” in each year.

# Group by year and country: by_year_country
by_year_country <- votes_processed %>%
  group_by(year, country) %>%
  summarize(total = n(),
            percent_yes = mean(vote == 1))

Plotting just the UK over time

Now that we have the percentage of time that each country voted “yes” within each year, you can plot the trend for a particular country. In this case, you’ll look at the trend for just the United Kingdom.

This will involve using filter() on your data before giving it to ggplot2.

# Start with by_year_country dataset
by_year_country <- votes_processed %>%
  group_by(year, country) %>%
  summarize(total = n(),
            percent_yes = mean(vote == 1))

# Print by_year_country
head(by_year_country)

# Create a filtered version: UK_by_year
UK_by_year <- by_year_country %>%
  filter(country == "United Kingdom")

# Line plot of percent_yes over time for UK only
ggplot(UK_by_year, aes(x = year, y = percent_yes)) +
  geom_line()

Plotting multiple countries

Plotting just one country at a time is interesting, but you really want to compare trends between countries. For example, suppose you want to compare voting trends for the United States, the UK, France, and India.

We’ll have to filter to include all four of these countries and use another aesthetic (not just x- and y-axes) to distinguish the countries on the resulting visualization. Instead, we’ll use the color aesthetic to represent different countries.

# Vector of four countries to examine
countries <- c("United States", "United Kingdom",
               "France", "India")

# Filter by_year_country: filtered_4_countries
filtered_4_countries <- by_year_country %>%
  filter(country %in% countries)

# Line plot of % yes in four countries
ggplot(filtered_4_countries, aes(x = year, y = percent_yes, color = country)) +
  geom_line()

Faceting by country

Faceting the time series

Now we’ll take a look at six countries. While in the previous exercise you used color to represent distinct countries, this gets a little too crowded with six.

Instead, we will facet, giving each country its own sub-plot. To do so, you add a facet_wrap() step after all of your layers.

# Vector of six countries to examine
countries <- c("United States", "United Kingdom",
               "France", "Japan", "Brazil", "India")

# Filtered by_year_country: filtered_6_countries
filtered_6_countries <- by_year_country %>%
  filter(country %in% countries)

# Line plot of % yes over time faceted by country
ggplot(filtered_6_countries, aes(x = year, y = percent_yes)) +
  geom_line() +
  facet_wrap(~country)

Faceting with free y-axis

All six graphs had the same axis limits. This made the changes over time hard to examine for plots with relatively little change.

Instead, we may want to let the plot choose a different y-axis for each facet.

# Vector of six countries to examine
countries <- c("United States", "United Kingdom",
               "France", "Japan", "Brazil", "India")

# Filtered by_year_country: filtered_6_countries
filtered_6_countries <- by_year_country %>%
  filter(country %in% countries)

# Line plot of % yes over time faceted by country
ggplot(filtered_6_countries, aes(year, percent_yes)) +
  geom_line() +
  facet_wrap(~ country, scales = "free_y")

More countries

# Add three more countries to this list
countries <- c("United States", "United Kingdom",
               "France", "Japan", "Germany", "Spain", "Ireland", "Russia", "China")

# Filtered by_year_country: filtered_countries
filtered_countries <- by_year_country %>%
  filter(country %in% countries)

# Line plot of % yes over time faceted by country
ggplot(filtered_countries, aes(year, percent_yes)) +
  geom_line() +
  facet_wrap(~ country, scales = "free_y")

Tidy modeling with broom

Linear regression

Linear regression on the United States

A linear regression is a model that lets us examine how one variable changes with respect to another by fitting a best fit line. It is done with the lm() function in R.

Here, we’ll fit a linear regression to just the percentage of “yes” votes from the United States.

# Percentage of yes votes from the US by year: US_by_year
US_by_year <- by_year_country %>%
  filter(country == "United States")

# Print the US_by_year data
head(US_by_year)

# Perform a linear regression of percent_yes by year: US_fit
US_fit <- lm(percent_yes~year, data = US_by_year)

# Perform summary() on the US_fit object
summary(US_fit)

Call:
lm(formula = percent_yes ~ year, data = US_by_year)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.222491 -0.080635 -0.008661  0.081948  0.194307 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) 12.6641455  1.8379743   6.890 8.48e-08 ***
year        -0.0062393  0.0009282  -6.722 1.37e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1062 on 32 degrees of freedom
Multiple R-squared:  0.5854,    Adjusted R-squared:  0.5724 
F-statistic: 45.18 on 1 and 32 DF,  p-value: 1.367e-07

What is the estimated slope of this relationship? -0.006

Finding the p-value of a linear regression

Not all positive or negative slopes are necessarily real. A p-value is a way of assessing whether a trend could be due to chance. Generally, data scientists set a threshold by declaring that, for example, p-values below .05 are significant.

What is the p-value of the relationship between year and percent_yes? 1.37e-07 It seems year is very significant when modeling percent_yes.

Tidy models with broom

Tidying a linear regression model

We’ll use the tidy() function in the broom package to turn that model into a tidy data frame.

# Load the broom package
library(broom)

# Call the tidy() function on the US_fit object
tidy(US_fit)

Combining models for multiple countries

One important advantage of changing models to tidied data frames is that they can be combined.

Earlier, we fit a linear model to the percentage of “yes” votes for each year in the United States. Now we’ll fit the same model for the United Kingdom and combine the results from both countries.

# Linear regression of percent_yes by year for US
US_by_year <- by_year_country %>%
  filter(country == "United States")
US_fit <- lm(percent_yes ~ year, US_by_year)

# Fit model for the United Kingdom
UK_by_year <- by_year_country %>%
  filter(country == "United Kingdom")
UK_fit <- lm(percent_yes ~ year, UK_by_year)

# Create US_tidied and UK_tidied
US_tidied <- tidy(US_fit)
UK_tidied <- tidy(UK_fit)

# Combine the two tidied models
bind_rows(US_tidied, UK_tidied)

Nesting for multiple models

Nesting a data frame

Right now, the by_year_country data frame has one row per country-vote pair. So that you can model each country individually, we’re going to “nest” all columns besides country, which will result in a data frame with one row per country. The data for each individual country will then be stored in a list column called data.

# Load the tidyr package
library(tidyr)
# Nest all columns besides country
head(by_year_country %>%
  nest(-country))

All elements of `...` must be named.
Did you want `data = c(year, total, percent_yes)`?

List columns

This “nested” data has an interesting structure. The second column, data, is a list, a type of R object that allows complicated objects to be stored within each row. This is because each item of the data column is itself a data frame. We can use nested$data to access this list column and double brackets to access a particular element. For example, nested$data[[1]] would give we the data frame with Afghanistan’s voting history (the percent_yes per year), since Afghanistan is the first row of the table.

# All countries are nested besides country
nested <- by_year_country %>%
  nest(-country)

All elements of `...` must be named.
Did you want `data = c(year, total, percent_yes)`?

# Print the nested data for Brazil
nested$data[[7]]

Unnesting

The opposite of the nest() operation is the unnest() operation. This takes each of the data frames in the list column and brings those rows back to the main data frame.

We’ll learn how to fit a model in between these nesting and unnesting steps that makes this process useful.

# All countries are nested besides country
nested <- by_year_country %>%
  nest(-country)

All elements of `...` must be named.
Did you want `data = c(year, total, percent_yes)`?

# Unnest the data column to return it to its original form
head(nested %>%
  unnest(data))

Fitting multiple models

Performing linear regression on each nested dataset

Now that we’ve divided the data for each country into a separate dataset in the data column, we need to fit a linear model to each of these datasets.

The map() function from purrr works by applying a formula to each item in a list, where . represents the individual item. For example, you could add one to each of a list of numbers:

map(numbers, ~ 1 + .)

This means that to fit a model to each dataset, we can do:

map(data, ~ lm(percent_yes ~ year, data = .))

where . represents each individual item from the data column in by_year_country. Recall that each item in the data column is a dataset that pertains to a specific country.

# Load tidyr and purrr
library(tidyr)
library(purrr)

# Perform a linear regression on each item in the data column
head(by_year_country %>%
  nest(-country) %>%
  mutate(model = map(data, ~ lm(percent_yes~year, data = .))))

All elements of `...` must be named.
Did you want `data = c(year, total, percent_yes)`?

Tidy each linear regression model

We’ve now performed a linear regression on each nested dataset and have a linear model stored in the list column model. But we can’t recombine the models until we’ve tidied each into a table of coefficients. To do that, we’ll need to use map() one more time and the tidy() function from the broom package.

Recall that we can simply give a function to map() (e.g. map(models, tidy)) in order to apply that function to each item of a list.

# Add another mutate that applies tidy() to each model
head(by_year_country %>%
  nest(-country) %>%
  mutate(model = map(data, ~ lm(percent_yes ~ year, data = .))) %>%
  mutate(tidied = map(model, tidy)))

All elements of `...` must be named.
Did you want `data = c(year, total, percent_yes)`?

Unnesting a data frame

We now have a tidied version of each model stored in the tidied column. We want to combine all of those into a large data frame, similar to how you combined the US and UK tidied models earlier. Recall that the unnest() function from tidyr achieves this.

# Add one more step that unnests the tidied column
country_coefficients <- by_year_country %>%
  nest(-country) %>%
  mutate(model = map(data, ~ lm(percent_yes ~ year, data = .)),
         tidied = map(model, tidy)) %>%
         unnest(tidied)

All elements of `...` must be named.
Did you want `data = c(year, total, percent_yes)`?

# Print the resulting country_coefficients variable
head(country_coefficients)

Working with many tidy models

Filtering model terms

We currently have both the intercept and slope terms for each by-country model. We’re probably more interested in how each is changing over time, so we want to focus on the slope terms.

# Print the country_coefficients dataset
head(country_coefficients)

# Filter for only the slope terms
head(country_coefficients %>%
  filter(term == "year"))

Filtering for significant countries

Not all slopes are significant, and we can use the p-value to guess which are and which are not.

However, when we have lots of p-values, like one for each country, we run into the problem of multiple hypothesis testing, where we have to set a stricter threshold. The p.adjust() function is a simple way to correct for this, where p.adjust(p.value) on a vector of p-values returns a set that we can trust.

Here we’ll add two steps to process the slope_terms dataset: use a mutate to create the new, adjusted p-value column, and filter to filter for those below a .05 threshold.

# Filter for only the slope terms
slope_terms <- country_coefficients %>%
  filter(term == "year")

# Add p.adjusted column, then filter
head(slope_terms %>%
  mutate(p.adjusted = p.adjust(p.value)) %>%
  filter(p.adjusted < 0.05))

Notice that there are now only 61 countries with significant trends.

Sorting by slope

Now that we’ve filtered for countries where the trend is probably not due to chance, we may be interested in countries whose percentage of “yes” votes is changing most quickly over time. Thus, we want to find the countries with the highest and lowest slopes; that is, the estimate column.

# Filter by adjusted p-values
filtered_countries <- country_coefficients %>%
  filter(term == "year") %>%
  mutate(p.adjusted = p.adjust(p.value)) %>%
  filter(p.adjusted < .05)

# Sort for the countries increasing most quickly
head(filtered_countries %>%
  arrange(estimate))

# Sort for the countries decreasing most quickly
head(filtered_countries %>%
  arrange(desc(estimate)))

Joining and tidying

##Joining datasets

We’ll now combine votes_processed dataset with the new descriptions dataset, which includes topic information about each country, so that we can analyze votes within particular topics.

To do this, we’ll make use of the inner_join() function from dplyr.

# Load dplyr package
library(dplyr)

# Print the votes_processed dataset
head(votes_processed)

# Print the descriptions dataset
descriptions <- readRDS("descriptions.rds")
head(descriptions)

# Join them together based on the "rcid" and "session" columns
votes_joined <- votes_processed %>%
  inner_join(descriptions, by= c("rcid" , "session"))

Column `rcid` has different attributes on LHS and RHS of joinColumn `session` has different attributes on LHS and RHS of join

Filtering the joined dataset

There are six columns in the descriptions dataset (and therefore in the new joined dataset) that describe the topic of a resolution:

1. me: Palestinian conflict

2. nu: Nuclear weapons and nuclear material

3. di: Arms control and disarmament

4. hr: Human rights

5. co: Colonialism

6. ec: Economic development

Each contains a 1 if the resolution is related to this topic and a 0 otherwise

# Filter for votes related to colonialism
votes_col <- votes_joined %>%
  filter(co == 1)
head(votes_col)

Visualizing colonialism votes

Earlier exercise, we graphed the percentage of votes each year where the US voted “yes”. Now we’ll create that same graph, but only for votes related to colonialism.

# Load the ggplot2 package
library(ggplot2)

# Filter, then summarize by year: US_co_by_year
US_co_by_year <- votes_joined %>%
  filter(country == "United States", co == 1) %>%
  group_by(year) %>%
  summarize(percent_yes = mean(vote==1))

# Graph the % of "yes" votes over time
ggplot(US_co_by_year, aes(x = year, y = percent_yes)) +
  geom_line()

Tidy data

Using gather to tidy a dataset

In order to represent the joined vote-topic data in a tidy form so we can analyze and graph by topic, we need to transform the data so that each row has one combination of country-vote-topic. This will change the data from having six columns (me, nu, di, hr, co, ec) to having two columns (topic and has_topic).

# Load the tidyr package
library(tidyr)

# Gather the six me/nu/di/hr/co/ec columns
votes_gathered <- votes_joined %>%
  gather(topic, has_topic, me:ec)
head(votes_gathered)

# Perform gather again, then filter
votes_gathered <- votes_joined %>%
  gather(topic, has_topic, me:ec) %>%
  filter(has_topic == 1)

Recoding the topics

There’s one more step of data cleaning to make this more interpretable. Right now, topics are represented by two-letter codes

So that we can interpret the data more easily, recode the data to replace these codes with their full name. We can do that with dplyr’s recode() function, which replaces values with ones you specify:

 example <- c("apple", "banana", "apple", "orange")
 recode(example,
        apple = "plum",
        banana = "grape")

# Replace the two-letter codes in topic: votes_tidied
votes_tidied <- votes_gathered %>%
  mutate(topic = recode(topic,
                        "me" = "Palestinian conflict",
                        "nu" = "Nuclear weapons and nuclear material",
                        "di" = "Arms control and disarmament",
                        "hr" = "Human rights",
                        "co" = "Colonialism",
                        "ec" = "Economic development"))

Summarize by country, year, and topic

Now that we have topic as an additional variable, you can summarize the votes for each combination of country, year, and topic (e.g. for the United States in 2013 on the topic of nuclear weapons.)

# Print votes_tidied
head(votes_tidied)

# Summarize the percentage "yes" per country-year-topic
by_country_year_topic <- votes_tidied %>%
  group_by(country, year, topic) %>%
  summarize(total = n(),
  percent_yes = mean(vote == 1)) %>%
  ungroup()

# Print by_country_year_topic
head(by_country_year_topic)

Visualizing trends in topics for one country

We can now visualize the trends in percentage “yes” over time for all six topics side-by-side. Here, we’ll visualize them just for the United States.

# Load the ggplot2 package
library(ggplot2)

# Filter by_country_year_topic for just the US
US_by_country_year_topic <- by_country_year_topic %>%
 filter(country == "United States")

# Plot % yes over time for the US, faceting by topic
ggplot(US_by_country_year_topic, aes(x = year, y = percent_yes)) +
  geom_line() +
  facet_wrap(~topic)

Tidy modeling by topic and country

Nesting by topic and country

We constructed a linear model for each country by nesting the data in each country, fitting a model to each dataset, then tidying each model with broom and unnesting the coefficients. The code looked something like this:

country_coefficients <- by_year_country %>%
  nest(-country) %>%
  mutate(model = map(data, ~ lm(percent_yes ~ year, data = .)),
         tidied = map(model, tidy)) %>%
  unnest(tidied)

Now, we’ll again be modeling change in “percentage” yes over time, but instead of fitting one model for each country, we’ll fit one for each combination of country and topic.

# Print by_country_year_topic
head(by_country_year_topic)

# Fit model on the by_country_year_topic dataset
country_topic_coefficients <- by_country_year_topic %>%
  nest(-country, -topic) %>%
  mutate(model = map(data, ~ lm(percent_yes ~ year, data = .)),
         tidied = map(model, tidy)) %>%
  unnest(tidied)

All elements of `...` must be named.
Did you want `data = c(year, total, percent_yes)`?essentially perfect fit: summary may be unreliableessentially perfect fit: summary may be unreliableessentially perfect fit: summary may be unreliableessentially perfect fit: summary may be unreliableessentially perfect fit: summary may be unreliableessentially perfect fit: summary may be unreliableessentially perfect fit: summary may be unreliableessentially perfect fit: summary may be unreliableessentially perfect fit: summary may be unreliableessentially perfect fit: summary may be unreliableessentially perfect fit: summary may be unreliableessentially perfect fit: summary may be unreliableessentially perfect fit: summary may be unreliableessentially perfect fit: summary may be unreliable

# Print country_topic_coefficients
head(country_topic_coefficients)

We can ignore the warning messages in the console for now.

Interpreting tidy models

Now we have both the slope and intercept terms for each model. Just as you did in the last chapter with the tidied coefficients, we’ll need to filter for only the slope terms.

We’ll also have to extract only cases that are statistically significant, which means adjusting the p-value for the number of models, and then filtering to include only significant changes.

# Create country_topic_filtered
country_topic_filtered <- country_topic_coefficients %>%
  filter(term == "year") %>%
  mutate(p.adjusted = p.adjust(p.value)) %>%
  filter(p.adjusted < .05)

Steepest trends by topic

Which combination of country and topic has the steepest downward trend?

head(country_topic_filtered %>% 
  arrange(estimate))

Vanuatu (an island nation in the Pacific Ocean) sharply changed its pattern of voting on the topic of Palestinian conflict.

Checking models visually

Let’s examine this country’s voting patterns more closely. Recall that the by_country_year_topic dataset contained one row for each combination of country, year, and topic. You can use that to create a plot of Vanuatu’s voting, faceted by topic.

# Create vanuatu_by_country_year_topic
vanuatu_by_country_year_topic <- by_country_year_topic %>%
  filter(country == "Vanuatu")

# Plot of percentage "yes" over time, faceted by topic
ggplot(vanuatu_by_country_year_topic, aes(x = year, y = percent_yes)) +
  geom_line() +
  facet_wrap(~topic)