Chapter 4 Effective data visualization

4.5.1 The Mauna Loa CO2 data set

This data set contains the atmospheric concentration of carbon dioxide (CO2, in parts per million) at the Mauna Loa research station in Hawaii from the years 1959-1997. Question: Does the concentration of atmospheric CO2 change over time, and are there any interesting patterns to note?

# mauna loa carbon dioxide data 
co2_df <- read_csv("data/maunaloa.csv")
head(co2_df)

## # A tibble: 6 x 2
##   concentration  date
##           <dbl> <dbl>
## 1          315. 1959 
## 2          316. 1959.
## 3          316. 1959.
## 4          318. 1959.
## 5          318. 1959.
## 6          318  1959.

Since we are investigating a relationship between two variables (CO2 concentration and date), a scatter plot is a good place to start. Scatter plots show the data as individual points with x (horizonal axis) and y (vertical axis) coordinates. Here, we will use the date as the x coordinate and CO2 concentration as the y coordinate. When using the ggplot2 library, we create the plot object with the ggplot function; there are a few basic aspects of a plot that we need to specify:

• the data: the name of the dataframe object that we would like to visualize
  • here, we specify the co2_df dataframe
• the aesthetic mapping: tells ggplot how the columns in the dataframe map to properties of the visualization
  • to create an aesthetic mapping, we use the aes function
  • here, we set the plot x axis to the date variable, and the plot y axis to the concentration variable
• the geometric object: specifies how the mapped data should be displayed
  • to create a geometric object, we use a geom_* function (see the ggplot reference for a list of geometric objects)
  • here, we use the geom_point function to visualize our data as a scatterplot

There are many other possible arguments we could pass to the aesthetic mapping and geometric object to change how the plot looks. For the purposes of quickly testing things out to see what they look like, though, we can just go with the default settings:

co2_scatter <- ggplot(co2_df, aes(x = date, y = concentration)) + 
        geom_point() 
co2_scatter

Certainly the visualization shows a clear upward trend in the atmospheric concentration of CO2 over time. This plot answers the first part of our question in the affirmative, but that appears to be the only conclusion one can make from the scatter visualization. However, since time is an ordered quantity, we can try using a line plot instead using the geom_line function. Line plots require that the data are ordered by their x coordinate, and connect the sequence of x and y coordinates with line segments. Let’s again try this with just the default arguments:

co2_line <- ggplot(co2_df, aes(x = date, y = concentration)) + 
        geom_line() 
co2_line

Aha! There is another interesting phenomenon in the data: in addition to increasing over time, the concentration seems to oscillate as well. Given the visualization as it is now, it is still hard to tell how fast the oscillation is, but nevertheless, the line seems to be a better choice for answering the question than the scatter plot was. The comparison between these two visualizations illustrates a common issue with scatter plots: often the points are shown too close together or even on top of one another, muddling information that would otherwise be clear (overplotting).

Now that we have settled on the rough details of the visualization, it is time to refine things. This plot is fairly straightforward, and there is not much visual noise to remove. But there are a few things we must do to improve clarity, such as adding informative axis labels and making the font a more readable size. In order to add axis labels we use the xlab and ylab functions. To change the font size we use the theme function with the text argument:

co2_line <- ggplot(co2_df, aes(x = date, y = concentration)) + 
                   geom_line() +
                   xlab('Year') +
                   ylab('Atmospheric CO2 (ppm)') + 
                   theme(text = element_text(size = 18))
co2_line

Finally, let’s see if we can better understand the oscillation by changing the visualization a little bit. Note that it is totally fine to use a small number of visualizations to answer different aspects of the question you are trying to answer. We will accomplish this by using scales, another important feature of ggplot2 that allow us to easily transform the different variables and set limits. We scale the horizontal axis by using the scale_x_continuous function, and the vertical axis with the scale_y_continuous function. We can transform the axis by passing the trans argument, and set limits by passing the limits argument. In particular, here we will use the scale_x_continuous function with the limits argument to zoom in on just five years of data (say, 1990-1995):

co2_line <- ggplot(co2_df, aes(x = date, y = concentration)) + 
                   geom_line() +
                   xlab('Year') +
                   ylab('Atmospheric CO2 (ppm)') + 
                   scale_x_continuous(limits = c(1990, 1995)) +
                   theme(text = element_text(size = 18))
co2_line

Interesting! It seems that each year, the atmospheric CO2 increases until it reaches its peak somewhere around April, decreases until around late September, and finally increases again until the end of the year. In Hawaii, there are two seasons: summer from May through October, and winter from November through April. Therefore, the oscillating pattern in CO2 matches up fairly closely with the two seasons.

4.5.2 The island landmass data set

This data set contains a list of Earth’s land masses as well as their area (in thousands of square miles). Question: Are the continents (North / South America, Africa, Europe, Asia, Australia, Antarctica) Earth’s 7 largest landmasses? If so, what are the next few largest landmasses after those?

# islands data 
islands_df <- read_csv("data/islands.csv")
head(islands_df)

## # A tibble: 6 x 2
##   landmass      size
##   <chr>        <dbl>
## 1 Africa       11506
## 2 Antarctica    5500
## 3 Asia         16988
## 4 Australia     2968
## 5 Axel Heiberg    16
## 6 Baffin         184

Here, we have a list of Earth’s landmasses, and are trying to compare their sizes. The right type of visualization to answer this question is a bar plot, specified by the geom_bar function in ggplot2. However, by default, geom_bar sets the heights of bars to the number of times a value appears in a dataframe (its count); here we want to plot exactly the values in the dataframe, i.e., the landmass sizes. So we have to pass the stat = "identity" argument to geom_bar:

islands_bar <- ggplot(islands_df, aes(x = landmass, y = size)) + 
            geom_bar(stat = "identity")
islands_bar

Alright, not bad! This is definitely the right kind of visualization, as we can clearly see and compare sizes of landmasses. The major issues are that the sizes of the smaller landmasses are hard to distinguish, and that the names of the landmasses are obscuring each other as they have been squished into too little space. But remember that the question we asked was only about the largest landmasses; let’s make the plot a little bit clearer by keeping only the largest 12 landmasses. We do this using the top_n function. Then to help us make sure the labels have enough space, we’ll use horizontal bars instead of vertical ones. We do this using the coord_flip function, which swaps the x and y coordinate axes:

islands_top12 <- top_n(islands_df, 12, size)
islands_bar <- ggplot(islands_top12, aes(x = landmass, y = size)) + 
        geom_bar(stat = "identity") + 
        coord_flip()
islands_bar

This is definitely clearer now, and allows us to answer our question (“are the top 7 largest landmasses continents?”) in the affirmative. But the question could be made clearer from the plot by organizing the bars not by alphabetical order but by size, and to colour them based on whether or not they are a continent. In order to do this, we use mutate to add a column to the data regarding whether or not the landmass is a continent:

islands_top12 <- top_n(islands_df, 12, size)
continents <- c('Africa', 'Antarctica', 'Asia', 'Australia', 'Europe', 'North America', 'South America')
islands_ct <- mutate(islands_top12, is_continent = ifelse(landmass %in% continents, 'Continent', 'Other'))
head(islands_ct)

## # A tibble: 6 x 3
##   landmass    size is_continent
##   <chr>      <dbl> <chr>       
## 1 Africa     11506 Continent   
## 2 Antarctica  5500 Continent   
## 3 Asia       16988 Continent   
## 4 Australia   2968 Continent   
## 5 Baffin       184 Other       
## 6 Borneo       280 Other

In order to colour the bars, we add the fill argument to the aesthetic mapping. Then we use the reorder function in the aesthetic mapping to organize the landmasses by their size variable. Finally, we use the labs and theme functions to add labels, change the font size, and position the legend:

islands_bar <- ggplot(islands_ct, aes(x = reorder(landmass, size), y = size, fill = is_continent)) + 
                   geom_bar(stat="identity") +
                   labs(x = 'Landmass', y = 'Size (1000 square mi)', fill = 'Type') +
                   coord_flip() +
                   theme(text = element_text(size = 18), legend.position = c(0.75, 0.45))
islands_bar

This is now a very effective visualization for answering our original questions. Landmasses are organized by their size, and continents are coloured differently than other landmasses, making it quite clear that continents are the largest 7 landmasses.

4.5.3 The Old Faithful eruption / waiting time data set

This data set contains measurements of the waiting time between eruptions and the subsequent eruption duration (in minutes). Question: Is there a relationship between the waiting time before an eruption to the duration of the eruption?

# old faithful eruption time / wait time data
head(faithful)

##   eruptions waiting
## 1     3.600      79
## 2     1.800      54
## 3     3.333      74
## 4     2.283      62
## 5     4.533      85
## 6     2.883      55

Here again we are investigating the relationship between two quantitative variables (waiting time and eruption time). But if you look at the output of the head function, you’ll notice that neither of the columns are ordered. So in this case, let’s start again with a scatter plot:

faithful_scatter <- ggplot(faithful, aes(x = waiting, y = eruptions)) + 
            geom_point()  
faithful_scatter

We can see that the data tend to fall into two groups: one with a short waiting and eruption times, and one with long waiting and eruption times. Note that in this case, there is no overplotting: the points are generally nicely visually separated, and the pattern they form is clear. In order to refine the visualization, we need only to add axis labels and make the font more readable:

faithful_scatter <- ggplot(faithful, aes(x = waiting, y = eruptions)) + 
                         geom_point() +
                         labs(x = 'Waiting Time (mins)', y = 'Eruption Duration (mins)') +
                         theme(text = element_text(size = 18))
faithful_scatter

4.5.4 The Michelson speed of light data set

This data set contains measurements of the speed of light (in kilometres per second with 299,000 subtracted) from the year 1879 for 5 experiments, each with 20 consecutive runs. Question: Given what we know now about the speed of light (299,792.458 kilometres per second), how accurate were each of the experiments?

# michelson morley experimental data
head(morley)

##     Expt Run Speed
## 001    1   1   850
## 002    1   2   740
## 003    1   3   900
## 004    1   4  1070
## 005    1   5   930
## 006    1   6   850

In this experimental data, Michelson was trying to measure just a single quantitative number (the speed of light). The data set contains many measurements of this single quantity. To tell how accurate the experiments were, we need to visualize the distribution of the measurements (i.e., all their possible values and how often each occurs). We can do this using a histogram. A histogram helps us visualize how a particular variable is distributed in a data set by separating the data into bins, and then using vertical bars to show how many data points fell in each bin. To create a histogram in ggplot2 we will use the geom_histogram geometric object, setting the x axis to the Speed measurement variable; and as we did before, let’s use the default arguments just to see how things look:

morley_hist <- ggplot(morley, aes(x = Speed)) + 
        geom_histogram() 
morley_hist

This is a great start. However, we cannot tell how accurate the measurements are using this visualization unless we can see what the true value is. In order to visualize the true speed of light, we will add a vertical line with the geom_vline function, setting the xintercept argument to the true value. There is a similar function, geom_hline, that is used for plotting horizontal lines. Note that vertical lines are used to denote quantities on the horizontal axis, while horizontal lines are used to denote quantities on the vertical axis.

morley_hist <- ggplot(morley, aes(x = Speed)) + 
                 geom_histogram() +
                 geom_vline(xintercept = 792.458, linetype = "dashed", size = 1.0)
morley_hist

We also still cannot tell which experiments (denoted in the Expt column) led to which measurements; perhaps some experiments were more accurate than others. To fully answer our question, we need to separate the measurements from each other visually. We can try to do this using a coloured histogram, where counts from different experiments are stacked on top of each other in different colours. We create a histogram coloured by the Expt variable by adding it to the fill aesthetic mapping. We make sure the different colours can be seen (despite them all sitting on top of each other) by setting the alpha argument in geom_histogram to 0.5 to make the bars slightly translucent:

morley_hist <- ggplot(morley, aes(x = Speed, fill = factor(Expt))) + 
                 geom_histogram(position = "identity", alpha = 0.5) +
                 geom_vline(xintercept = 792.458, linetype = "dashed", size = 1.0)
morley_hist

Unfortunately, the attempt to separate out the experiment number visually has created a bit of a mess. All of the colours are blending together, and although it is possible to derive some insight from this (e.g., experiments 1 and 3 had some of the most incorrect measurements), it isn’t the clearest way to convey our message and answer the question. Let’s try a different strategy of creating multiple separate histograms on top of one another.

In order to create a plot in ggplot2 that has multiple subplots arranged in a grid, we use the facet_grid function. The argument to facet_grid specifies the variable(s) used to split the plot into subplots. It has the syntax vertical_variable ~ horizontal_variable, where veritcal_variable is used to split the plot vertically, horizontal_variable is used to split horizontally, and . is used if there should be no split along that axis. In our case we only want to split vertically along the Expt variable, so we use Expt ~ . as the argument to facet_grid.

morley_hist <- ggplot(morley, aes(x = Speed, fill = factor(Expt))) + 
                 geom_histogram(position = "identity") +
                 facet_grid(Expt ~ .) +
                 geom_vline(xintercept = 792.458, linetype = "dashed", size = 1.0)
morley_hist

The visualization now makes it quite clear how accurate the different experiments were with respect to one another. There are two finishing touches to make this visualization even clearer. First and foremost, we need to add informative axis labels using the labs function, and increase the font size to make it readable using the theme function. Second, and perhaps more subtly, even though it is easy to compare the experiments on this plot to one another, it is hard to get a sense for just how accurate all the experiments were overall. For example, how accurate is the value 800 on the plot, relative to the true speed of light? To answer this question we’ll use the mutate function to transform our data into a relative measure of accuracy rather than absolute measurements:

morley_rel <- mutate(morley, relative_accuracy = 100*( (299000 + Speed) - 299792.458 ) / (299792.458))
morley_hist <- ggplot(morley_rel, aes(x = relative_accuracy, fill = factor(Expt))) + 
                 geom_histogram(position = "identity") +
                 facet_grid(Expt ~ .) +
                 geom_vline(xintercept = 0, linetype = "dashed", size = 1.0) + 
                 labs(x = 'Relative Accuracy (%)', y = '# Measurements', fill = 'Experiment ID') + 
                 theme(text = element_text(size = 18))
morley_hist

Wow, impressive! These measurements of the speed of light from 1879 had errors around 0.05% of the true speed. This shows you that even though experiments 2 and 5 were perhaps the most accurate, all of the experiments did quite an admirable job given the technology available at the time period.