Empirical Project 5 Working in R

R-specific learning objectives

In addition to the learning objectives for this project, in this section you will learn how to use loops to repeat specified tasks for a list of values (Note: this is an extension task so may not apply to all users).

Getting started in R

For this project you will need the following packages:

If you need to install any of these packages, run the following code:

install.packages(c("readxl","tidyverse","ineq","reshape2"))

You can import these libraries now, or when they are used in the R walk-throughs below.

library(readxl)
library(tidyverse)
library(ineq)
library(reshape2)

Part 5.1 Measuring income inequality

One way to visualize the income distribution in a population is to draw a Lorenz curve. This curve shows the entire population lined up along the horizontal axis from the poorest to the richest. The height of the curve at any point on the vertical axis indicates the fraction of total income received by the fraction of the population given by that point on the horizontal axis.

We will start by using income decile data from the Global Consumption and Income Project to draw Lorenz curves and compare changes in the income distribution of a country over time. Note that income here refers to market income, which does not take into account taxes or government transfers (see Section 5.9 of Economy, Society, and Public Policy for further details).

To answer the question below:

R walk-through 5.1 Importing an Excel file (either .xlsx or .xls format) into R

As we are dealing with an Excel file, we use the read_excel function from the readxl package. The file is called GCIPrawdata.xlsx, but before you import the file into R, open the datafile in Excel to understand its structure. You will see that the data is on one worksheet (which is convenient), and that the headings for the variables are in the third row. Hence we will use the skip = 2 option in the read_excel function to skip the first two rows.

library(tidyverse)  
library(readxl)  
decile_data <- read_excel("GCIPrawdata.xlsx", skip = 2)  

The data is now in a tibble (like a spreadsheet for R). Let’s look at the first few rows:

head(decile_data)  
## # A tibble: 6 x 14
##   Country      Year `Decile 1 Income` `Decile 2 Income` `Decile 3 Income`
##   <chr>       <dbl>             <dbl>             <dbl>             <dbl>
## 1 Afghanistan  1980               206               350               455
## 2 Afghanistan  1981               212               361               469
## 3 Afghanistan  1982               221               377               490
## 4 Afghanistan  1983               238               405               527
## 5 Afghanistan  1984               249               424               551
## 6 Afghanistan  1985               256               435               566
## # ... with 9 more variables: `Decile 4 Income` <dbl>, `Decile 5
## #   Income` <dbl>, `Decile 6 Income` <dbl>, `Decile 7 Income` <dbl>,
## #   `Decile 8 Income` <dbl>, `Decile 9 Income` <dbl>, `Decile 10
## #   Income` <dbl>, `Mean Income` <dbl>, Population <dbl>

As you can see, we have an entry (row) for every country and every year. The first entry (for Afghanistan in the Year 1980) is 206, and it is the value for the variable Decile 1 Income. This value indicates that the mean annual income of the poorest 10% in Afghanistan was the equivalent of 206 US Dollars (in 1980, adjusted using purchasing power parity). The mean income of the next richest 10% (those in the 11th to 20th percentiles for income) was 350.

To see the list of variables, we examine the structure of decile_data.

str(decile_data)  
## Classes 'tbl_df', 'tbl' and 'data.frame':  4799 obs. of 14 variables:
##  $ Country         : chr  "Afghanistan" "Afghanistan" "Afghanistan" "Afghanistan" ...
##  $ Year            : num  1980 1981 1982 1983 1984 ...
##  $ Decile 1 Income : num  206 212 221 238 249 256 268 243 223 202 ...
##  $ Decile 2 Income : num  350 361 377 405 424 435 457 414 380 344 ...
##  $ Decile 3 Income : num  455 469 490 527 551 566 594 539 493 447 ...
##  $ Decile 4 Income : num  556 574 599 644 674 692 726 659 603 547 ...
##  $ Decile 5 Income : num  665 686 716 771 806 828 869 788 722 654 ...
##  $ Decile 6 Income : num  793 818 854 919 961 ...
##  $ Decile 7 Income : num  955 986 1029 1107 1157 ...
##  $ Decile 8 Income : num  1187 1225 1278 1376 1438 ...
##  $ Decile 9 Income : num  1594 1645 1717 1848 1932 ...
##  $ Decile 10 Income: num  3542 3655 3814 4105 4291 ...
##  $ Mean Income     : num  1030 1063 1109 1194 1248 ...
##  $ Population      : num  13211412 12996923 12667001 12279095 11912510 ...

In addition to the country, year, and the ten income deciles we have mean income and the population.

  1. Choose two countries. You will be using their data, for 1980 and 2014, as the basis for your Lorenz curves. Use the country data you have selected to calculate cumulative income shares. (Remember that each decile represents 10% of the population.)

R walk-through 5.2 Calculating cumulative shares using the cumsum function

Here we have chosen China (a country that recently underwent enormous economic changes) and the US (a developed country).

sel_Year <- c(1980,2014)
sel_Country <- c("United States","China")
temp <- subset(decile_data, (decile_data$Country %in% sel_Country) & (decile_data$Year %in% sel_Year)) #  Select the data for the chosen country and years
temp
## # A tibble: 4 x 14
##   Country        Year `Decile 1 Income` `Decile 2 Income` `Decile 3 Incom~
##   <chr>         <dbl>             <dbl>             <dbl>            <dbl>
## 1 China          1980                79               113              146
## 2 China          2014               448               927             1440
## 3 United States  1980              3392              5820             7855
## 4 United States  2014              3778              6534             9069
## # ... with 9 more variables: `Decile 4 Income` <dbl>, `Decile 5
## #   Income` <dbl>, `Decile 6 Income` <dbl>, `Decile 7 Income` <dbl>,
## #   `Decile 8 Income` <dbl>, `Decile 9 Income` <dbl>, `Decile 10
## #   Income` <dbl>, `Mean Income` <dbl>, Population <dbl>

Before we calculate cumulative income shares, we need to calculate the total income using the mean income and the population size.

print("Total incomes are:")  
## [1] "Total incomes are:"  
total_income <-temp[,"Mean Income"]*temp[,"Population"]  
total_income  
##   Mean Income  
## 1 2.472624e+11  
## 2 6.609944e+12  
## 3 3.366422e+12  
## 4 6.401280e+12  

These numbers are very large, so for our purpose it is easier to assume that there is only one person in each decile, in other words the total income is 10 times the mean income. This simplification works because, by definition, each decile has exactly the same number of people (10% of the population).

We will be using the very useful cumsum function to calculate the cumulative income. To see what this function does, look at this simple example.

test <- c(2,4,10,22)  
cumsum(test)  
## [1] 2 6 16 38  

With this functionality in mind, we now calculate the cumulative income shares for China (1980).

decs_c80 <- unlist(temp[1,3:12])   # Pick the 10 deciles (Columns 3 to 12) in Row 1 (China, 1980)
                               # The unlist function transforms temp[1,3:12] from a tibble to simple vector with data which simplifies the calculations.
total_inc <- 10*unlist(temp[1,"Mean Income"])  # Give the total income, assuming a population of 10
cum_inc_share_c80 = cumsum(decs_c80)/total_inc
cum_inc_share_c80
##  Decile 1 Income  Decile 2 Income  Decile 3 Income  Decile 4 Income 
##       0.03134921       0.07619048       0.13412698       0.20436508 
##  Decile 5 Income  Decile 6 Income  Decile 7 Income  Decile 8 Income 
##       0.28769841       0.38492063       0.49841270       0.63174603 
##  Decile 9 Income Decile 10 Income 
##       0.79206349       0.99841270

We repeat the same process for China in 2014 and for the US in 1980 and 2014.

# For China, 2014  
decs_c14 <- unlist(temp[2,3:12]) # Go to Row 2 (China, 2014)
total_inc <- 10*unlist(temp[2,"Mean Income"]) # Give the total income, assuming a population of 10
cum_inc_share_c14 = cumsum(decs_c14)/total_inc  
  
# For the US, 1980  
decs_us80 <- unlist(temp[3,3:12]) # Select Row 3 (USA, 1980)  
total_inc <- 10*unlist(temp[3,"Mean Income"]) # Give the total income, assuming a population of 10
cum_inc_share_us80 = cumsum(decs_us80)/total_inc  
  
# For the US, 2014  
decs_us14 <- unlist(temp[4,3:12]) # Select Row 4 (USA, 2014)  
total_inc <- 10*unlist(temp[4,"Mean Income"]) # Give the total income, assuming a population of 10 
cum_inc_share_us14 = cumsum(decs_us14)/total_inc  
  1. Use the cumulative income shares to draw Lorenz curves for each country in order to visually compare the income distributions over time.

R walk-through 5.3 Drawing Lorenz curves

Let us plot the cumulative income shares for China (1980).

plot(cum_inc_share_c80, type = "l",col="blue",lwd = 2,
     ylab="Cumulative income share")  
abline(a=0,b=0.1,col="black",lwd=2) # Add the perfect equality line 
title("Lorenz curve, China, 1980") 

Lorenz curve, China, 1980.

Figure 5.1 Lorenz curve, China, 1980.

The blue line is the Lorenz curve. The Gini coefficient is the ratio of the area between the two lines and the total area under the black line. We will calculate that in the R walk-through 5.4.

Now we add the other Lorenz curves to the chart.

plot(cum_inc_share_c80, type = "l",col="blue", lty = 2, lwd=2,
     xlab = "Deciles", ylab="Cumulative income share")  
abline(a=0,b=0.1,col="black",lwd=2) # Add the perfect equality line  
lines(cum_inc_share_c14,col="green",lty = 1,lwd=2) # lty = 2 = solid line  
lines(cum_inc_share_us80,col="red", lty = 2,lwd=2) # lty = 1 = dashed line  
lines(cum_inc_share_us14,col="orange", lty = 1,lwd=2)  
title("Lorenz curves, China and the US (1980 and 2014)")  
legend("topleft", legend=c("China, 1980", "China, 2014", "US, 1980", "US, 2014"),  
col=c("blue", "green","red", "orange"), lty=2:1, lwd=2,cex=1.2)  

Lorenz curves, China and the US (1980 and 2014).

Figure 5.2 Lorenz curves, China and the US (1980 and 2014).

As the chart shows, the income distribution has changed more clearly for China than for the US.

  1. Using your Lorenz curves:

A rough way to compare income distributions is to use a summary measure such as the Gini coefficient. The Gini coefficient ranges from 0 (complete equality) to 1 (complete inequality). It is calculated by dividing the area between the Lorenz curve and the perfect equality line, by the total area underneath the perfect equality line. Intuitively, the further away the Lorenz curve is from the perfect equality line, the more unequal the income distribution is, and the higher the Gini coefficient will be.

To calculate the Gini coefficient you can either use a Gini coefficient calculator, or calculate it directly in R as shown in R walk-through 5.4.

  1. Calculate the Gini coefficient for each of your Lorenz curves. You should have four coefficients in total. Label each Lorenz curve with its corresponding Gini coefficient, and check that the coefficients are consistent with what you see in your charts.

R walk-through 5.4 Calculating Gini coefficients

In Section 5.8 of Economy, Society, and Public Policy you can learn that the Gini coefficient is graphically represented by dividing the area between the perfect equality line and the Lorenz curve by the total area under the perfect equality line. You could calculate this area by hand, by decomposing the area under the Lorenz curve into rectangles and triangles, but as with so many problems, someone else has already figured out how to do that and has provided R users with a package (ineq) to do this task very easily.

library(ineq) # Load the ineq library

g_c80<- Gini( decs_c80 ) # The decile mean incomes from R walk-through 5.3 are used.
g_c14<- Gini( decs_c14 )
g_us80<- Gini( decs_us80 )
g_us14<- Gini( decs_us14 )
paste("Gini coefficients")
## [1] "Gini coefficients"
paste("China - 1980: ", round(g_c80,2), ", 2014: ", round(g_c14,2))
## [1] "China - 1980: 0.29 , 2014: 0.51"
paste("United States - 1980: ", round(g_us80,2), ", 2014: ",round( g_us14,2))
## [1] "United States - 1980: 0.34 , 2014: 0.4"

Now we make the same line chart (copy and paste the code from R walk-through 5.3, but use the text command to label curves with their respective Gini coefficients.

plot(cum_inc_share_c80, type = "l",col="blue", lty = 2, lwd=2,
     xlab = "Deciles", ylab="Cumulative income share")   
abline(a=0,b=0.1,col="black", lwd=2) # Add the perfect equality line
lines(cum_inc_share_c14,col="green",lty = 1,lwd=2)  # lty = 2 = solid line
lines(cum_inc_share_us80,col="red", lty = 2,lwd=2)  # lty = 1 = dashed line
lines(cum_inc_share_us14,col="orange", lty = 1,lwd=2)
title("Lorenz curves, China and the US (1980 and 2014)")
legend("topleft", legend=c("China, 1980", "China, 2014", "US, 1980", "US, 2014"),
       col=c("blue", "green","red", "orange"), lty=2:1,lwd=2, cex=1.2)
text(8.5, 0.78,round(g_c80,digits = 3))
text(9.4, 0.6,round(g_c14,digits = 3))
text(5.7, 0.38,round(g_us80,digits = 3))
text(6.4, 0.3,round(g_us14,digits = 3))

Lorenz curves, China and the US (1980 and 2014).

Figure 5.3 Lorenz curves, China and the US (1980 and 2014).

The Gini coefficients have increased, confirming what we already saw from the Lorenz curves that in both countries the income distribution has become more unequal.

Extension R walk-through 5.5 Calculating Gini coefficients for all countries and all years using a loop

In this extention walk-though, we show you how to calculate the Gini coefficient for all countries and years in your dataset.

This sounds like a tedious task, and indeed if we were to use the same method as before it would be mind-numbing. However, we have a powerful programming language at hand, and this is the time to use it.

Here we use a very useful programming tool you may not have come across yet, which is loops. Let’s start with a very simple case: printing the values for i^2 for values of i=1,...,10.

for (i in seq(1,10)){
  print(i^2)
}
## [1] 1
## [1] 4
## [1] 9
## [1] 16
## [1] 25
## [1] 36
## [1] 49
## [1] 64
## [1] 81
## [1] 100

In the above command, seq(1,10) creates a vector of data (1,2,3,…,10). The command for (i in seq(1,10)) defines the variable i initially as 1, then performs all the commands that are between the curly brackets for each value of i (typically these commands will involve the variable i). Here our command prints the value of i^2 for each value of i.

Now we use loops to complete our task. We begin by creating a new variable in our dataset, gini, which we initially set to 0 for all country-year combinations.

decile_data$gini <- 0

Now we use a loop to run through all the rows in our dataset and for each row we will repeat the Gini coefficient calculation from R walk-through 5.4 then save the resulting value in the gini variable we created.

noc <- nrow(decile_data) # Give us the number of rows in decile_data

for (i in seq(1,noc)){
  decs_i <- unlist(decile_data[i,3:12]) # Go to Row I to get the decile data
  decile_data$gini[i] <- Gini( decs_i )
}

With this code, we calculated 4,799 Gini coefficients. We now look at the summary stats for the gini variable.

summary(decile_data$gini)
## Min.   1st Qu. Median Mean   3rd Qu. Max.
## 0.1791 0.3470  0.4814 0.4617 0.5700  0.7386

The average Gini coefficient is 0.46, the maximum is 0.74, and the minimum 0.18. Let’s look at these extreme cases.

First we will look at the extremely equal income distributions:

temp <- subset(decile_data, decile_data$gini < 0.20, select = c("Country","Year","gini"))
temp
## # A tibble: 17 x 3
##    Country          Year  gini
##    <chr>           <dbl> <dbl>
##  1 Bulgaria         1987 0.191
##  2 Czech Republic   1985 0.195
##  3 Czech Republic   1986 0.194
##  4 Czech Republic   1987 0.192
##  5 Czech Republic   1988 0.191
##  6 Czech Republic   1989 0.194
##  7 Czech Republic   1990 0.196
##  8 Czech Republic   1991 0.199
##  9 Slovak Republic  1985 0.195
## 10 Slovak Republic  1986 0.194
## 11 Slovak Republic  1987 0.193
## 12 Slovak Republic  1988 0.192
## 13 Slovak Republic  1989 0.193
## 14 Slovak Republic  1990 0.194
## 15 Slovak Republic  1991 0.195
## 16 Slovak Republic  1992 0.196
## 17 Slovak Republic  1993 0.179

These correspond to eastern European countries before the fall of communism.

Now the most unequal countries:

temp <- subset(decile_data, decile_data$gini > 0.73, select = c("Country","Year","gini"))
temp
## # A tibble: 27 x 3
##    Country       Year  gini
##    <chr>        <dbl> <dbl>
##  1 Burkina Faso  1980 0.738
##  2 Burkina Faso  1981 0.738
##  3 Burkina Faso  1982 0.738
##  4 Burkina Faso  1983 0.738
##  5 Burkina Faso  1984 0.738
##  6 Burkina Faso  1985 0.738
##  7 Burkina Faso  1986 0.738
##  8 Burkina Faso  1987 0.738
##  9 Burkina Faso  1988 0.738
## 10 Burkina Faso  1989 0.739
## # ... with 17 more rows

Extension R walk-through 5.6 Plotting time-series of Gini coefficients, using ggplot

In this extension walk-through, we show you how to make time series plots (time on the horizontal axis, the variable of interest on the vertical axis) with Gini coefficients for a list of countries of your choice.

There are many ways to plot data in R, one being the standard plotting function we used in previous walk-throughs. Another (and perhaps more beautiful) way is to use the ggplot function, which is part of the tidyverse package we loaded earlier. Our dataset is already in a format which the ggplot function can easily use.

First we select a small list of countries. As an example, we have chosen four anglophone countries: the UK, the US, Ireland, and Australia.

temp_data <- subset(decile_data, Country %in% c("United Kingdom","United States","Ireland","Australia"))

Now we plot the data using ggplot.

ggplot(temp_data,aes(x =Year, y=gini, color=Country)) +
  geom_line(size=1) +
  theme_bw() +
  ggtitle("Gini coefficients for anglophone countries") # Add a title

Time series plots of Gini coefficients for anglophone countries.

Figure 5.4 Time series plots of Gini coefficients for anglophone countries.

We asked the ggplot function to use the decile_data dataframe/tibble, with Year on the horizontal axis and gini on the vertical axis. The color option indicates which variable we use to separate the data (Country). The first line of code sets up the chart, and the + geom_line(size=1) then instructs R to draw lines. (See what happens if you replace + geom_line(size=1) with + geom_point(size=1).)

ggplot assumes that the different lines you want to show are identified through the different values in one variable (here, the Country variable). If your data is formatted differently, for example, if you have one variable for the Gini of each country, then in order to use ggplot you will first have to transform the dataset. Doing so is beyond the scope of this task, however you can find a worked example online, such as ‘R TSplots’.1

The ggplot set of commands are extremely powerful, and if you want to produce a variety of different charts, you may want to read more about that package, for example, see a Harvard R tutorial or an R statistics tutorial for great examples including code.

Now we will look at other measures of income inequality and see how they can be used along with the Gini coefficient to summarize a country’s income distribution. Instead of summarizing the entire income distribution like the Gini coefficient does, we can take the ratio of incomes at two points in the distribution. For example, the 90/10 ratio takes the ratio of the top 10% of incomes (Decile 10) to the lowest 10% of incomes (Decile 1). A 90/10 ratio of 5 means that the richest 10% earns 5 times more than the poorest 10%. The higher the ratio, the higher the inequality between these two points in the distribution.

  1. Look at the following ratios:

    • 90/10 ratio = ratio of the Decile 10 income to the Decile 1 income
    • 90/50 ratio = ratio of the Decile 10 income to the Decile 5 income (the median)
    • 50/10 ratio = ratio of the Decile 5 income (the median) to the Decile 1 income.

We will now compare these summary measures (ratios and the Gini coefficient) for a larger group of countries, using OECD data. The OECD has annual data for different ratio measures of income inequality for 42 countries around the world, and has an interactive chart function that plots this data for you.

Go to the OECD website to access the data. You will see a chart similar to Figure 5.5, which shows data for 2015. The countries are ranked from smallest to largest Gini coefficient on the horizontal axis, and the vertical axis gives the Gini coefficient.

  1. Compare summary measures of inequality:

OECD countries ranked according to their Gini coefficient.

Figure 5.5 OECD countries ranked according to their Gini coefficient.

The Gini coefficient and the ratios we have used are common measures of inequality, but there are other ways to measure income inequality.

  1. Go to the ‘Income Inequality’ section of the Our world in data website, and choose two other measures of income inequality that you find interesting.

Part 5.2 Measuring other kinds of inequality

There are many ways to measure income inequality, but income inequality is only one dimension of inequality within a country. To get a more complete picture of inequality within a country, we need to look at other areas in which there may be inequality in outcomes. We will explore two particular areas:

First, we will look at how researchers have measured inequality in health-related outcomes. Besides income, health is an important aspect of wellbeing because it determines how long an individual will be alive to enjoy his or her income. If two people had the same annual income throughout their lives, but one person had a much shorter life than the other, we might say that the distribution of wellbeing is unequal, despite annual incomes being equal.

As with income, inequality in life expectancy can be measured using a Gini coefficient. In the study ‘Mortality inequality’, researcher Sam Peltzman (2009) estimated Gini coefficients for life expectancy based on the distribution of total years lived (life-years) across people born in a given year (birth cohort). If everybody born in a given year lived the same number of years, then the total years lived would be divided equally among these people (perfect equality). If a few people lived very long lives but everybody else lived very short lives, then there would be a high degree of inequality (Gini coefficient close to 1).

We will now look at mortality inequality Gini coefficients for 10 countries around the world. First, download the data:

Import the data into R and investigate the structure of the data as explained in R walk-through 5.7.

R walk-through 5.7 Importing .csv files into R

Before importing, save the csv file in your working directory.

health_in <- read.csv("inequality-of-life-as-measured-by-mortality-gini-coefficient-1742-2002.csv") # Open the csv file from the working directory
str(health_in)
## 'data.frame':  320 obs. of 4 variables:
##  $ Entity    : Factor w/ 10 levels "Brazil","England and Wales",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Code      : Factor w/ 10 levels "","BRA","DEU",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ Year      : int  1892 1897 1902 1907 1912 1917 1922 1927 1932 1937 ...
##  $ X.percent.: num  0.566 0.557 0.547 0.482 0.494 ...

The variable Entity is the country and the variable X.percent is the health Gini. Let’s change these variable names to make them more intuitive for our analysis.

names(health_in)[1] <- "Country" # Country is the first variable.
names(health_in)[4] <- "HGini" # Health Gini is the fourth variable.

There is another quirk in the data that you may not have noticed in this initial data inspection: All countries have a short code (Code), except for England and Wales for which that field is empty (or formally ""). Let’s change the blanks to ENW.

levels(health_in$Code)[levels(health_in$Code)==""] <- "ENW"

Tip

The way this code works may seem a little mysterious, and you may find it difficult to remember the code for this step. However, an Internet search for ‘R renaming one factor level’ (recall that Code is a factor variable) will show you many ways to achieve this (including that shown above). Often you will find answers on stackoverflow.com, where experienced coders provide useful help.

  1. Using the mortality inequality data:

R walk-through 5.8 Creating line graphs with ggplot

As shown in R walk-through 5.7, the data is already formatted so that we can use ggplot directly, in other words we have only one variable for the mortality Gini (HGini), and we can separate the data by country using one variable (Country).

temp_data <- subset(health_in, Year > 1951) # Select all data after 1951

ggplot(temp_data,aes(x =Year, y=HGini, color=Country)) +
  geom_line(size=1) +
  labs(y="Mortality inequality Gini coefficient") +
  scale_color_brewer(palette="Paired") + # Change the colour palette
  theme_bw() +
  ggtitle("Mortality inequalities") # Add a title

Mortality inequality Gini coefficients (1952–2002).

Figure 5.6 Mortality inequality Gini coefficients (1952–2002).

  1. Now compare the Gini coefficients in the first year of your line chart (1952) with the last year (2002).

R walk-through 5.9 Drawing a column chart with sorted values

Plot a column chart for 1952

First we use subset to extract the data for 1952 only.

temp_data <- subset(health_in, Year == 1952)    # Select all data for 1952
temp_data <- temp_data[order(temp_data$HGini),] # Reorder the rows in temp_data according to the values of HGini, from smallest to largest

temp_data
##               Country Code Year     HGini
## 279            Sweden  SWE 1952 0.1194045
## 46  England and Wales  ENW 1952 0.1319542
## 310     United States  USA 1952 0.1471329
## 138           Germany  DEU 1952 0.1572112
## 86             France  FRA 1952 0.1605238
## 228             Spain  ESP 1952 0.1985371
## 184             Japan  JPN 1952 0.2021728
## 206            Russia  RUS 1952 0.2237161
## 161             India  IND 1952 0.3978703
## 13             Brazil  BRA 1952 0.4103805

The rows are now ordered according to HGini, in ascending order. Let’s use ggplot again.

ggplot(temp_data, aes(x=Code, y=HGini)) +
  geom_bar(stat="identity", width=.5, fill="tomato3") +
  theme_bw() +
  labs(title="Mortality Gini coefficients (1952)",
    caption="source: https://ourworldindata.org/health-inequality",y="Mortality inequality Gini coefficient")

Mortality Gini coefficients (1952).

Figure 5.7 Mortality Gini coefficients (1952).

Unfortunately, the columns are not ordered correctly, because when the horizontal axis variable (here, Code) is a factor, then ggplot uses the ordering of the factor levels, which is:

levels(temp_data$Code)
## [1] "ENW" "BRA" "DEU" "ESP" "FRA" "IND" "JPN" "RUS" "SWE" "USA"

A blog post from Data Se provides the following code for ‘R geom_bar change order’, and uses the reorder function to reorder the horizontal axis variable (Code) according to the HGini value.

ggplot(temp_data, aes(x=reorder(Code,HGini), y=HGini)) +
  geom_bar(stat="identity", width=.5, fill="tomato3") +
  coord_cartesian(ylim=c(0,0.45)) +
  theme_bw() +
  labs(title="Mortality Gini coefficients (1952)", x="Country",
    caption="source: https://ourworldindata.org/health-inequality",y="Mortality inequality Gini coefficient")

Mortality Gini coefficients (1952).

Figure 5.8 Mortality Gini coefficients (1952).

Plot a column chart for 2002

We want to compare this ranking with the ranking of 2002. First we extract the relevant data again.

temp_data <- subset(health_in, Year == 2002) # Select all data for 2002 

ggplot(temp_data, aes(x=reorder(Code,HGini), y=HGini)) +  
  geom_bar(stat="identity", width=.5, ylim = c(0,0.45),fill="tomato3") +  
  coord_cartesian(ylim=c(0,0.45)) + # Adjust the ylim (vertical axis scale) to ensure comparability with 1952
  theme_bw() +  
  labs(title="Mortality Gini coefficients (2002)", x="Country",  
    caption="source: https://ourworldindata.org/health-inequality",y="Mortality inequality Gini coefficient")   

Mortality Gini coefficients (2002).

Figure 5.9 Mortality Gini coefficients (2002).

It is fairly easy to plot the data for both years in the same chart.

temp_data <- subset(health_in, Year %in% c("1952","2002")) # Select all data for 1952 and 2002
temp_data$Year <- factor(temp_data$Year)

ggplot(temp_data, aes(x=reorder(Code,HGini), y=HGini, fill=Year)) +
  geom_bar(position="dodge", stat="identity") +
  theme_bw() +
  labs(title="Mortality Gini coefficients (1952 and 2002)", x="Country",
    caption="source: https://ourworldindata.org/health-inequality",y="Mortality inequality Gini coefficient")

Mortality Gini coefficients (1952 and 2002).

Figure 5.10 Mortality Gini coefficients (1952 and 2002).

Now the country ordering is in terms of the average HGini, rather than HGini in 1952 (which would have made comparisons easier).

Other measures of health inequality, such as those used by the World Health Organization (WHO), are based on access to healthcare, affordability of healthcare, and quality of living conditions. Choose one of the following measures of health inequality to answer Question 3:

To download the data for your chosen measure:

  1. For your chosen measure:

R walk-through 5.10 Drawing a column chart with sorted values

For this walk-through, we downloaded the ‘access to essential medicines’ data, as explained above. Here we saved it as WHO access to essential medicines.csv. Looking at the spreadsheet in Excel, you can see that the actual data starts in row three, meaning there are two header rows. So let’s skip the first row when uploading it.

med_access <- read.csv("WHO access to essential medicines.csv",skip = 1)
 str(med_access)
## 'data.frame':  38 obs. of 3 variables:
##  $ Country     : Factor w/ 38 levels "Afghanistan",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ X2007.2013  : num  94 42.9 86.7 76.7 72.1 58.3 13.3 90.7 31.3 33.3 ...
##  $ X2007.2013.1: num  81.1 43.2 31.9 0 87.1 46.7 15.5 86.7 21.2 100 ...

The second and third variables have lost their labels. From the spreadsheet you know that they are:

  • median availability of selected generic medicines (%) – Private
  • median availability of selected generic medicines (%) – Public.

Let’s change the names of these variables to make working with them easier:

names(med_access)[2] <- "Private_Access"
names(med_access)[3] <- "Public_Access"

To find details about these variables, click the ‘Metadata’ button on the website to find the following explanation:

A standard methodology has been developed by WHO and Health Action International (HAI). Data on the availability of a specific list of medicines are collected in at least four geographic or administrative areas in a sample of medicine dispensing points. Availability is reported as the percentage of medicine outlets where a medicine was found on the day of the survey.

Before we produce charts of the data we shall look at summary statistics of the access variable.

summary(med_access)
##                              Country   Private_Access   Public_Access   
##  Afghanistan                     : 1   Min.   :  2.80   Min.   :  0.00  
##  Bahamas                         : 1   1st Qu.: 54.62   1st Qu.: 39.67  
##  Bolivia (Plurinational State of): 1   Median : 70.15   Median : 55.95  
##  Brazil                          : 1   Mean   : 65.97   Mean   : 58.25  
##  Burkina Faso                    : 1   3rd Qu.: 86.70   3rd Qu.: 82.50  
##  Burundi                         : 1   Max.   :100.00   Max.   :100.00  
##  (Other)                         :32                    NA's   :2

On average, private sector patients have better access to essential medication.

From the summary statistics for the Public_Access variable, you can see that there are two missing observations. Here, we will keep these observations because leaving them in doesn’t affect the following analysis.

med_access <- med_access[complete.cases(med_access),]

There are a number of interesting aspects to look at. We shall produce a bar chart comparing the private and public access in countries.

med_access$Country <- reorder(med_access$Country,med_access$Private_Access) # Reorders according to values of private access (largest to smallest)
library(reshape2) # This is required for the melt function.
med_access_melt <- melt(med_access) # Rearrange the data for ggplot
# This creates a dataframe with three columns
# Country = Country name
# value = access percentage (either Private_Access or Public_Access).
# variable = indicates whether a row refers to Public_Access or Private_Access. 

ggplot(med_access_melt, aes(x=Country, y=value, fill = variable)) +
  geom_bar(position="dodge", stat="identity") +
  scale_fill_discrete(name  ="Access",labels=c("Private sector", "Public sector")) +
  theme(axis.text.x=element_text(angle=90,hjust=1,vjust=0.5)) +
  theme_bw() +
  labs(title="Access to essential medication", x="Country",y="Percent of patients with access to essential medication") +
  coord_flip()  # Flip axis to make country labels readable

Access to essential medication.

Figure 5.11 Access to essential medication.

Let’s find the extreme values. There are two countries where public sector patients have access to all (100%) essential medications.

med_access[med_access$Public_Access == 100,]
##               Country Private_Access Public_Access
## 10       Cook Islands           33.3           100
## 30 Russian Federation          100.0           100

Let’s see which countries provide 0% access to essential medication for people in the public sector.

med_access[med_access$Public_Access == 0,]
##   Country Private_Access Public_Access
## 4 Brazil  76.7           0

Since an individual’s income and available options in later life partly depend on their level of education, inequality in educational access or attainment can lead to inequality in income and other outcomes. We will focus on the aspect of gender inequality in educational attainment, using data from the Our world in data website, to make our own comparisons between countries and over time. Choose one of the following measures to answer Question 4:

To download the data for your chosen measure:

  1. For your chosen measure:

R walk-through 5.11 Using line and bar charts to illustrate changes in time

Import data and plot a line chart

First we import the data into R.

data_prim <- read.csv("OWID-gender-gap-in-primary-education.csv")  # Open the csv file from the working directory
str(data_prim)
## 'data.frame':  8780 obs. of 4 variables:
##  $ Entity                                         : Factor w/ 250 levels "Afghanistan",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Code                                           : Factor w/ 207 levels "","ABW","AFG",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ Year                                           : int  1970 1971 1972 1973 1974 1975 1976 1977 1978 1980 ...
##  $ Primary.education..pupils....female.....female.: num  14.1 13.7 14 14.5 14.4 ...

The data is now in the dataframe data_prim. The variable of interest (percentage of female enrolment) has a very long name so we will shorten it to PFE.

names(data_prim)[4] <- "PFE"

As usual, ensure that you understand the definition of the variables you are using. In the Our world in data website, look at the ‘Sources’ tab underneath the graph for a definition:

Percentage of female enrollment is calculated by dividing the total number of female students at a given level of education by the total enrolment at the same level, and multiplying by 100.

This definition implies that if the primary-school-age population was 50% male and 50% female and all children were enrolled in school, the percentage of female enrolment would be 50.

Before choosing ten countries, we check which countries are in the dataset using the unique(data_prim$Entity) command. Here we only show the first few countries using the head() command.

head(unique(data_prim$Entity))
## [1] Afghanistan         Albania             Algeria            
## [4] Andorra             Angola              Antigua and Barbuda
## 250 Levels: Afghanistan Albania Algeria Andorra ... Zimbabwe

You can find nearly all the countries in the world in this list (plus some sub- and supra-country entities, like OECD countries, which explains why the variable wasn’t initially called ‘country’).

Plot a line chart for a selection of countries

We now make a selection of ten countries. (You can of course make a different selection, but ensure that you get the spelling right as R is unforgiving!).

temp_data<- subset(data_prim, Entity %in% c("Albania","China","France", "India","Japan", "Switzerland", "United Arab Emirates", "United Kingdom","Zambia","United States"))

Now we plot the data, similar to what we did earlier.

ggplot(temp_data,aes(x =Year, y=PFE, color=Entity)) +
  geom_line(size = 1) +     # size = 1 sets the line thickness.
  theme_bw() +   # Remove grey background
  scale_colour_brewer(palette="Paired") +   # Change the set of colours used 
  scale_colour_discrete(name  ="Country") +
  ylab("Percentage (%)") +   # Set the vertical axis label
  ggtitle("Female pupils as a percentage of total enrolment in primary education") # Add a title

Female pupils as a percentage of total enrolment in primary education.

Figure 5.12 Female pupils as a percentage of total enrolment in primary education.

Plot a column chart with sorted values

To calculate the change in the value of this measure between 1980 and 2010 for each country chosen, we have to manipulate the data so that we have one entry (row) for each entity (or country), but two different variables for the percentage of female enrolment PFE (one for each year).

temp_data_80 <- subset(temp_data, Year == "1980") # Select all data for 1980
names(temp_data_80)[4] <- "PFE_80" # Rename variable to include year
temp_data_10 <- subset(temp_data, Year == "2010") # Select all data for 2010
names(temp_data_10)[4] <- "PFE_10" # Rename variable to include year
temp_data2 <- merge(temp_data_80,temp_data_10,by=c("Entity"))

Have a look at temp_data2, which now contains two variables for every country, PFE_80 and PFE_10. It also has multiple variables for Year (Year.x and Year.y) and Code (Code.x and Code.y), but that is a minor issue and you could delete one of them.

Now we can calculate the difference.

temp_data2$dPFE <- temp_data2$PFE_10 - temp_data2$PFE_80

You could plot a separate chart for each year and check the order, but here we show how to create one chart with the data from both years.

ggplot(temp_data2, aes(x=reorder(Code.x,dPFE), y=dPFE)) +
  geom_bar(stat="identity",fill="tomato3") +
  labs(title="Change (%) in female pupils’ share of total enrolment in primary education",
       x="Country", y="Percentage change (%)",
       caption="source: https://ourworldindata.org/educational-mobility-inequality") +
  theme_bw()

Change in percentage of female enrolment in primary school from 1980 to 2010.

Figure 5.13 Change in percentage of female enrolment in primary school from 1980 to 2010.

It is apparent that some countries saw very little or no change (the countries that already had very high PFE). The countries with initially low female participation have significantly improved.

  1. University of Manchester’s Econometric Computing Learning Resource (ECLR). 2018. ‘R TSplots’. Updated 26 July 2016.