Ch 7 - Exploratory Data Analysis

7.1 Introduction

This chapter will show you how to use visualisation and transformation to explore your data in a systematic way, a task that statisticians call exploratory data analysis, or EDA for short. EDA is an iterative cycle. You:

1. Generate questions about your data.
2. Search for answers by visualising, transforming, and modelling your data.
3. Use what you learn to refine your questions and/or generate new questions.

EDA is an important part of any data analysis, even if the questions are handed to you on a platter, because you always need to investigate the quality of your data. Data cleaning is just one application of EDA: you ask questions about whether your data meets your expectations or not. To do data cleaning, you’ll need to deploy all the tools of EDA: visualisation, transformation, and modelling.

7.1.1 Prerequisites

In this chapter we’ll combine what you’ve learned about dplyr and ggplot2 to interactively ask questions, answer them with data, and then ask new questions.

library(tidyverse)

7.2 Questions

Goal during EDA is to develop an understanding of your data.

EDA is fundamentally a creative process. And like most creative processes, the key to asking quality questions is to generate a large quantity of questions.

There is no rule about which questions you should ask to guide your research. However, two types of questions will always be useful for making discoveries within your data. You can loosely word these questions as:

1. What type of variation occurs within my variables?
2. What type of covariation occurs between my variables?

The rest of this chapter will look at these two questions.

• A variable is a quantity, quality, or property that you can measure.
• A value is the state of a variable when you measure it. The value of a variable may change from measurement to measurement.
• An observation is a set of measurements made under similar conditions (you usually make all of the measurements in an observation at the same time and on the same object). An observation will contain several values, each associated with a different variable. I’ll sometimes refer to an observation as a data point.
• Tabular data is a set of values, each associated with a variable and an observation. Tabular data is tidy if eavh value is placed in its own “cell”, each variable in its own column, and each observation in its own row.

7.3 Variation

Variation is the tendency of the values of a variable to change from measurement to measurement. You can see variation easily in real life; if you measure any continuous variable twice, you will get two different results. Every variable has its own pattern of variation, which can reveal interesting information. The best way to understand that pattern is to visualise the distribution of the variable’s values.

7.3.1 Visualising distributions

How you visualise the distribution of a variable will depend on whether the variable is categorical or continuous. A variable is categorical if it can only take one of a small set of values. In R, categorical variables are usually saved as factors or character vectors. To examine the distribution of a categorical variable, use a bar chart:

ggplot(data = diamonds) +
  geom_bar(mapping = aes(x = cut))

The height of the bars displays how many observations occurred with each x value. You can compute these values manually with dply::count():

diamonds %>%
  count(cut)

A variable is continuous if it can take any of an infinite set of ordered values. Numbers and date-times are two examples of continuous variables. To examine the distribution of a continuous variable, use a histogram:

ggplot(data = diamonds) +
  geom_histogram(mapping = aes(x = carat), binwidth = .5)

You can compute this by hand by combining dply::count() and ggplot2::cut_width():

diamonds %>%
  count(cut_width(carat, .5))

A histogram divides the x-axis into equally spaced bins and then uses the height of a bar to display the number of observations that fall in each bin. In the graph above, the tallest abar shows that almost 30,000 observations have a varat value between 0.25 and 0.75, which are the left and right edges of the bar.

You can set the width of the intervals in a histogram with the binwidth argument, which is measured in the units of the x variable. You should always explore a variety of binwidths when working with histograms, as different binwidths can reveal different patters. For example, here is how the graph above looks when we zoom into just the diamonds with a size of less than three carats and choose a small binwidth.

diamonds %>%
  filter(carat < 3) %>%
  ggplot(mapping = aes(x = carat)) +
  geom_histogram(binwidth = .1)

If you wish to overlay multiple histograms in the same plot, I recommend using geom_freqpoly() instead of geom_histogram(). geom_freqpoly() performs the same calculation as geom_histogram(), but instead of displaying the counts with bars, uses lines instead. It’s much easier to understand overlapping lines than bars.

diamonds %>%
  filter(carat < 3) %>%
  ggplot(mapping = aes(x = carat, color = cut)) +
  geom_freqpoly()

There are a few challenges with this type of plot, which we will come back to in visualising a categorical and a continuous variable. Now that we can visualise variation, how do we ask good questions and further explore? The key to asking good follow-up questions will be to rely on your curiousity (What do you want to learn more about?) as well as your skepticism (How could this be misleading?).

7.3.2 Typical values

In both bar charts and histograms, tall bars show the common values of a variable, and shorter bars show less-common values. Places that do not have bars reveal values that were not seen in your data. To turn this information into useful questions, look for anything unexpected:

• Which values are the most common? Why?
• Which values are rare? Why? Does that match your expectations?
• Can you see any unusual patterns? What might explain them?

As an example, the histogram below suggests several interesting questions:

• Why are there more diamonds at whole carats and common fractions of carats?
• Why are there more diamonds slightly to the right of each peak than there are slightly to the left of each peak?
• Why are there no diamonds bigger than 3 carats?

ggplot(data = smaller, mapping = aes(x = carat)) +
  geom_histogram(binwidth = .01)

Clusters of similar values suggest that subgroups exist in your data. To understand the subgroups, ask:

• How are the observations within each cluster similar to each other?
• How are the observations in separate clusters different from each other?
• How can you explain or describe the clusters?
• Why might the appearance of clusters be misleading?

The histogram below shows the length (in minutes) of 272 eruptions of the Old Faithful geyser in Yellowstone National Park. Eruptions appear to be clustered into two groups: those that were for about 2 minutes (short) and longer eruptions of about 4-5 minutes, but little in between.

ggplot(data = faithful, mapping = aes(x = eruptions)) +
  geom_histogram()

Many of the questions above will prompt you to explore a relationship between variables, for example, to see if the values of one variable can explain the behavior of another variable.

7.3.3 Unusual values

Outliers are observations that are unusual; data points that don’t seem to fit the pattern. Sometimes outliers are data entry errors; other times outliers suggest important new science. When you have a lot of data, outliers are sometimes difficult to see in a histogram. For example, take the distribution of the y variable from the diamonds dataset. The only evidence of outliers is the unusually wide limits on the x-axis.

ggplot(data = diamonds, mapping = aes(x = y)) +
  geom_histogram(binwidth = .5)

There are so many observations in the common bins that the rare bins are so short that you can’t see them. To make it easy to see the unusual values, we need to zoom to small values of the y-axis with coord_cartesian():

ggplot(data = diamonds, mapping = aes(x = y)) +
  geom_histogram() +
  coord_cartesian(ylim = c(0, 50))

(_coord_cartesian()__ also has an xlim() argument for when you need to zoom into the x-axis. ggplot2 also has xlim() and ylim() functions that work slightly different: they throw away the data outside the limits.)

This allows us to see that there are three unusual values: 0, ~30, and ~60. We pluck them out with dplyr:

unusual <- diamonds %>%
  filter(y < 3 | y > 20) %>%
  select(price, x, y, z) %>%
  arrange(y)

unusual

The y variable measures one of the three dimensions of these diamonds, in mm. We know that diamonds can’t have a width of 0mm, so these values must be incorrect. We might also suspect that measurements of 32mm and 59mm are implausible: those diamonds are over an inch long, but don’t cost hundreds of thousands of dollars!

It’s good practice to repeat your analysis with and without the outliers. If they have minimal effect on the results, and you can’t figure out why they’re there, it’s reasonable to replace them with missing values, and move on. However, if they have a substantial effect on your results, you shouldn’t drop them without justification. You’ll need to figure out what caused them (e.g. a data entry error) and disclose that you removed them in your write-up.

7.3.4 Exercises

1. Explore the distribution of each of the x, y, and z variables in diamonds. What do you learn? Think about a diamond and how you might decide which dimension is the length, width, and depth.

summary(select(filtered, x, y, z))

       x                y                z        
 Min.   : 3.730   Min.   : 3.680   Min.   :1.070  
 1st Qu.: 4.710   1st Qu.: 4.720   1st Qu.:2.910  
 Median : 5.700   Median : 5.710   Median :3.530  
 Mean   : 5.732   Mean   : 5.733   Mean   :3.539  
 3rd Qu.: 6.540   3rd Qu.: 6.540   3rd Qu.:4.040  
 Max.   :10.740   Max.   :10.540   Max.   :6.980  

filtered %>%
  filter(x < 10) %>%
  ggplot(mapping = aes(x = x)) +
  geom_histogram(binwidth = .01) +
  scale_x_continuous(breaks = 1:10)

filtered %>%
  filter(y < 10) %>%
  ggplot(mapping = aes(x = y)) +
  geom_histogram(binwidth = .01) +
  scale_x_continuous(breaks = 1:10)

filtered %>%
  filter(z < 10) %>%
  ggplot(mapping = aes(x = z)) +
  geom_histogram(binwidth = .01) +
  scale_x_continuous(breaks = 1:10)

2. Explore the distribution of price. Do you discover anything unusual or surprising? (Hint: Carefully think about the binwidth and make sure you try a wide range of values.)

summary(select(diamonds, price))

     price      
 Min.   :  326  
 1st Qu.:  950  
 Median : 2401  
 Mean   : 3933  
 3rd Qu.: 5324  
 Max.   :18823  

sd(diamonds$price)

[1] 3989.44

There is a large amount of variation between the price of the “cheapest” and most expensive diamonds. In fact, the most expensive diamond is greater than 3 standard deviations from the 3rd quartile.

ggplot(data = diamonds, mapping = aes(x = price)) +
  geom_histogram(binwidth = 1000)

ggplot(data = diamonds, mapping = aes(x = price)) +
  geom_histogram(binwidth = 100)

The distribution of price is similar to the distributions of x, y, and z. The data is right-skewed with most diamonds costing less than $5,000 and the largest cluster of diamonds costing less than $2,500. Beginning ~$2,500 the slope of the reduction becomes very predictable. Essentially, as the price increases about $1,000 the count of diamonds decreases by a proportional amount.

The graph with the narrower binwidth demonstrates some peculariaties. There appears to be a reversal of the trend for diamonds leading up to $5,000 and an oddly few amount of diamonds cost ~$1,250.

3. How many diamonds are .99 carat? How many are 1 carat? What do you think is the cause of the difference?

nrow(filter(diamonds, carat == .99))

[1] 23

nrow(filter(diamonds, carat == 1))

[1] 1558

I am not sure what the cause is? Maybe vanity? People want at least one carat for the diamond to be desirable? There is also likely a premium on full carat diamonds.

The distribution of carats shows that the carats peak at or near intervals of .5.

ggplot(data = diamonds, mapping = aes(x = carat)) +
  geom_histogram(binwidth = .1)

diamonds %>%
  filter(carat >= .9 & carat <= 1.1) %>%
  count(carat) %>%
  print(n = Inf)

4. Compare and contrast coord_cartesian() vs xlim() or ylim() when zooming in on a histogram. What happens if you leave binwidth unset? What happens if you try and zoom so only half a bar shows?

ggplot(data = diamonds, mapping = aes(x = x)) +
  geom_histogram()

ggplot(data = diamonds, mapping = aes(x = x)) +
  geom_histogram() +
  coord_cartesian(xlim = c(2, 10), ylim = c(0, 5000))

ggplot(data = diamonds) +
  geom_histogram(mapping = aes(x = x)) +
  xlim(3, 9) +
  ylim(0, 5000)

Because the coord_cartesian function is specified after the creation of the histogram, it zooms in on the coordinates, but does not change the shape of the histogram.

However, the x and y lim operators remove values prior to the creation of the histogram and thus influence the shape of the histogram that is produced.

7.4 Missing values

If you’ve encountered unusual values in your dataset, and simply want to move on to the rest of your analysis, you have two options.

1. Drop the entire row with the strange values:

diamonds2 <- diamonds %>%
  filter(between(y, 3, 20))

I don’t recommend this option because just because one measurement is invalid, doesn’t mean all the measurements are. Additionally, if you have low quality data, by time that you’ve applied this approach to every variable you might find that you don’t have any data left!

2. Instead, I recommend replacing the unusual values with missing values. The easiest way to do this is to use mutate() to replace the variable with a modified copy. You can use the ifelse() function to replace unusual values with NA:

diamonds2 <- diamonds %>%
  mutate(y = ifelse(y < 3 | y > 20, NA, y))

ifelse() has three arguments. The first argument test should be a logical vector. The result will contain the value of the second argument, yes, when test is TRUE, and the value of the third argument, no when it is false. Alternatively to ifelse, use dplyr::case_when(). case_when is particularly useful inside mutate when you want to create a new variable that relies on a complex combination of existing variables.

Like R, ggplot2 subscribes to the philosophy that missing values should never silently go missing. It’s not obvious where you should plot missing values, so ggplot2 doesn’t include them in the plot, but it does warn that they’ve been removed:

ggplot(data = diamonds2, mapping = aes(x = x, y = y)) +
  geom_point()

To suppress that warning, set na.rm = TRUE.

ggplot(data = diamonds2, mapping = aes(x = x, y = y)) +
  geom_point(na.rm = TRUE)

Other times you want to understand what makes observations with missing values different to observations with recorded values. For example, in nycflights13::flights, missing values in the dep_time variable indicate that the flight was cancelled. So you might want to compare the scheduled departure times for cancelled and non-cancelled times. You can do this by making a new variable with is.na().

flights %>%
  mutate(
    cancelled = is.na(dep_time),
    sched_hour = sched_dep_time %/% 100,
    sched_min = sched_dep_time %% 100,
    sched_dep_time = sched_hour + sched_min / 60
    ) %>%
  ggplot(mapping = aes(x = sched_dep_time)) +
  geom_freqpoly(mapping = aes(color = cancelled, binwidth = 1/4)) 

Ignoring unknown aesthetics: binwidth

7.4.1 Exercises

1. What happens to missing values in a histogram? What happens to missing values in a bar chart? Why is there a difference?

diamonds %>%
  mutate(y = ifelse(y < 3 | y > 20, NA, y)) %>%
  ggplot(mapping = aes(x = y)) +
  geom_histogram()

When producing a histogram, the histogram removes the missing values. This is likely because a histogram relies on continuous (numeric) data to plot the frequency of counts according to bins. When it does not know the value, it does not know which bin to put it in, so it removes it.

Alternatively,

diamonds %>%
  mutate(cut = ifelse(runif(n()) < .1, NA_character_, as.character(cut))) %>%
  ggplot(mapping = aes(x = cut)) +
  geom_bar()

In a bar chart, the graphic relies on categorical data and assumes NA is a character string indicating another category.

2. What does na.rm = TRUE do in mean() and sum().

diamonds %>%
  mutate(xyz = ifelse(z == 0, NA, z)) %>%
  select(x,y,z,xyz) %>%
  arrange(z) %>%
  summarise(
    total_xyz = sum(xyz, na.rm = TRUE),
    mean_xyz = mean(xyz, na.rm = TRUE)
  )

For this exercise, we just converted all of the “0s” in z to “NA” and then applied sums and means. When NA’s occur in the data, the sum and mean values are unsure how to calculate them so it returns values of NA. na.rm = TRUE removes the NA values before the calculation is performed though and then returns the sum and mean of the new column.

7.5 Covariation

If variation describes the behavior within a variable, covariation describes the behavior between variables. Covariation is the tendency for the values of two or more variables to vary together in a related way. The best way to spot covariation is to visualise the relationship between two or more variables. How you do that should again depend on the type of variables involved.

7.5.1 A categorical and continuous variable

It’s common to want to explore the distribution of a continuous variable broken down by a categorical variable, as in the previous frequency polygon. The default appearnace of geom_freqpoly() is not that useful for that sort of comparison because the height is given by the count. That means if one of the groups is much smaller than the other, it’s hard to see the differences in shape. For example, let’s explore how the price of a diamond varies with its quality:

ggplot(data = diamonds, mapping = aes(x = price)) +
  geom_freqpoly(mapping = aes(color = cut), binwidth = 500)

It’s hard to see the difference in distribution because the overall counts differ so much:

ggplot(data = diamonds) +
  geom_bar(mapping = aes(x = cut))

To make the comparison easier we need to swap what is displayed on the y-axis. Instead of displaying count, we’ll display density, which is the count standardised so that the area under each frequency polygon is one.

ggplot(data = diamonds, mapping = aes(x = price, y = ..density..)) +
  geom_freqpoly(mapping = aes(color = cut), binwidth = 500)

This graph is hard to interpert.

Another alternative to display the distribution of a continuous variable broken down by a categorical variable is the boxplot. A boxplot is a type of visual shorthand for a distribution of values that is popular among statisticians. Each boxplot consists of:

• A box that stretches from the 25th percentile of the distribution to the 75th percentile, a distance known as the intequartile range (IQR). In the middle of the box is a line that displays the median, i.e. 50th percentile, of the distribution. These three lines give you a sense of the spread of the distribution and whether or not the distribution is symmetric about the median or skewed to one side.

• Visual points that display observations that fall more than 1.5 times the IQR from either edge of the box. Theese outlying points are unusual so are plotted individually.

• A line (or whisker) that extends from each end of the box and goes to the farthest non-outlier point in the distribution.

Let’s take a look at the distribution of price by cut using geom_boxplot():

ggplot(data = diamonds) +
  geom_boxplot(mapping = aes(x = cut, y = price))

We see much less information about the distribution, but the boxplots are much more compact so we can more easily compare them (and fit more on one plot). It supports the counterintuitive finding that better quality diamonds are cheaper on average! In the exercises, you’ll be challenged to figure out why. *My assumption is that these graphs only tell you something about the cuts of the diamonds, but nothing about their size or color. It could be that there is a significantly larger portion of ideal diamonds that are relatively smaller than the other cuts.

cut is an ordered factor: fair is worse than good, which is worse than very good and so on. Many categorical variables don’t have such an intrisic order, so you might want to reorder them to make a more informative display. One way to do that is with the reorder() function.

For example, take the class variable in the mpg dataset. You might be interested to know how highway mileage varies across classes:

ggplot(data = mpg, mapping = aes(x = class, y = hwy)) +
  geom_boxplot()

To make the trend easier to see, we can reorder class based on the median value of hwy:

ggplot(data = mpg) +
  geom_boxplot(mapping = aes(x = reorder(class, hwy, FUN = median), y = hwy))

If you have long variable names, geom_boxplot() will work better if you flip it 90 degrees. You can do that with coord_flip().

ggplot(data = mpg) +
  geom_boxplot(mapping = aes(x = reorder(class, hwy, FUN = median), y = hwy)) +
  coord_flip()

7.5.1.1 Exercises

1. Use whay you’ve learned to improve the visualisation of the departure times of cancelled vs. non-cancelled flights.

cancelled_flights %>%
  group_by(year, month, day) %>%
  summarise(
    prop_cancelled = mean(cancelled)
  ) %>%
  ggplot() +
  geom_boxplot(mapping = aes(x = reorder(as.character(month), prop_cancelled, FUN = median), y = prop_cancelled))

`summarise()` has grouped output by 'year', 'month'. You can override using the `.groups` argument.

cancelled_flights %>%
  group_by(year, month) %>%
  summarise(
    prop_cancelled = mean(cancelled)) %>%
  arrange(prop_cancelled)

`summarise()` has grouped output by 'year'. You can override using the `.groups` argument.

cancelled_flights %>%
  mutate(
    sched_dep_hour = sched_dep_time %/% 100,
    sched_dep_min = sched_dep_time %% 100,
    sched_dep_time = sched_dep_hour + sched_dep_min / 60
  ) %>%
  ggplot() +
  geom_boxplot(mapping = aes(x = cancelled, y = sched_dep_time))

cancelled_flights %>%
  mutate(
    sched_dep_hour = sched_dep_time %/% 100,
    sched_dep_min = sched_dep_time %% 100,
    sched_dep_time = sched_dep_hour + sched_dep_min / 60
  ) %>%
  ggplot(mapping = aes(x = sched_dep_time, y = ..density..)) +
  geom_freqpoly(mapping = aes(color = cancelled))

The frequency polygon and the boxplot essentially display the same thing, that the majority of cancelled flights occur later in the day, ~ 15, or 3 pm. Note however that we have to supply a density argument to the frequency polygon graph in order to control for there being significantly more non-cancelled flights than cancelled flights.

2. What variable in the diamonds dataset is most important for predicting the price of a diamond? How is that variable correlated with cut? Why does the combination of those two relationships lead to lower quality diamonds being more expensive?

colnames(diamonds)

 [1] "carat"   "cut"     "color"   "clarity" "depth"   "table"   "price"   "x"      
 [9] "y"       "z"      

summary(diamonds)

     carat               cut        color        clarity          depth      
 Min.   :0.2000   Fair     : 1610   D: 6775   SI1    :13065   Min.   :43.00  
 1st Qu.:0.4000   Good     : 4906   E: 9797   VS2    :12258   1st Qu.:61.00  
 Median :0.7000   Very Good:12082   F: 9542   SI2    : 9194   Median :61.80  
 Mean   :0.7979   Premium  :13791   G:11292   VS1    : 8171   Mean   :61.75  
 3rd Qu.:1.0400   Ideal    :21551   H: 8304   VVS2   : 5066   3rd Qu.:62.50  
 Max.   :5.0100                     I: 5422   VVS1   : 3655   Max.   :79.00  
                                    J: 2808   (Other): 2531                  
     table           price             x                y                z         
 Min.   :43.00   Min.   :  326   Min.   : 0.000   Min.   : 0.000   Min.   : 0.000  
 1st Qu.:56.00   1st Qu.:  950   1st Qu.: 4.710   1st Qu.: 4.720   1st Qu.: 2.910  
 Median :57.00   Median : 2401   Median : 5.700   Median : 5.710   Median : 3.530  
 Mean   :57.46   Mean   : 3933   Mean   : 5.731   Mean   : 5.735   Mean   : 3.539  
 3rd Qu.:59.00   3rd Qu.: 5324   3rd Qu.: 6.540   3rd Qu.: 6.540   3rd Qu.: 4.040  
 Max.   :95.00   Max.   :18823   Max.   :10.740   Max.   :58.900   Max.   :31.800  
                                                                                   

ggplot(data = diamonds) +
  geom_histogram(mapping = aes(x = carat, fill = cut), binwidth = .1) +
  coord_cartesian(xlim = c(0, 3))

ggplot(data = diamonds) +
  geom_bar(mapping = aes(x = cut))

ggplot(data = diamonds, mapping = aes(x = color, y = price)) +
  geom_boxplot()

ggplot(data = diamonds, mapping = aes(x = clarity, y = price)) +
  geom_boxplot()

ggplot(data = diamonds, mapping = aes(x = x, y = price)) +
  geom_point() +
  xlim(c(3, 10))

ggplot(data = diamonds, mapping = aes(x = y, y = price)) +
  geom_point() +
  xlim(c(1, 12))

ggplot(data = diamonds, mapping = aes(x = z, y = price)) +
  geom_point() +
  xlim(c(1, 7.5))

Our previous hypothesis that ideal diamonds were underrepresented in the data was actually incorrect. In fact, it is the opposite. Ideal diamonds are the most common cut of data and fair diamonds are actually the smallest proportion.

ggplot(data = diamonds, mapping = aes(x = carat, y = price)) +
  geom_point() +
  geom_smooth()

The scatterplots of price against carat, x, y, and z are all very similar which leads me to believe that x, y, and z might be determinants of carat? In this case, we would expect that carat would be correlated with this variables.

ggplot(data = diamonds, mapping = aes(x = x, y = carat)) +
  geom_point() +
  geom_smooth(se = FALSE) +
  xlim(c(3, 10))

ggplot(data = diamonds, mapping = aes(x = y, y = carat)) +
  geom_point() +
  geom_smooth(se = FALSE) +
  xlim(c(1, 15))

ggplot(data = diamonds, mapping = aes(x = z, y = carat)) +
  geom_point() +
  geom_smooth(se = FALSE) +
  xlim(c(1, 10))

And yes in fact, carat and x, y, and z are all highly correlated (especially after filtering out previously identified likely outliers).

ggplot(data = diamonds) +
  geom_boxplot(mapping = aes(x = cut, y = carat))

The boxplot above shows that ideal cuts have the lowest median value of carat. This means that although ideal diamonds, all else equal, are more expensive, they generally are also much smaller, whereas fair diamonds are on average “the largest” diamonds in our dataset.

3. Installe the ggstance package, and create a horizontal boxplot. How does this compare to using coord_flip()?

ggplot(data = diamonds) +
  geom_boxplot(mapping = aes(x = cut, y = carat)) +
  coord_flip()

ggplot(data = diamonds) +
  geom_boxploth(mapping = aes(x = carat, y = cut))

They basically do the same thing, but the ggstance package assumes a horizontal graph. GGplots can do the same thing you just have to specify the axes or use orientation to specify.

4. One problem with boxplots is that they were developed in an era of much smaller datasets and tend to display a prohibitevly large number of “outlying values”. One approach to remedy this problem is the letter value plot. Install the lvplot package, and try using geom_lv() to display the distribution of price vs cut. What do you learn? How do you interpret the plots?

library(lvplot)
library(ggplot2)
ggplot(data = diamonds) +
  geom_lv(mapping = aes(x = cut, y = price))

What we see from the letter value plot is that while it works in the same way as a boxplot, it uses many more quartiles and thus gives a more accurate representation of the quantiles of the dataset. This is especially useful for larger datasets which should, in theory, have more (in absolute terms) outliers than smaller datasets.

5. Compare and contrast geom_violin() with a facetted geom_histogram(), or a colored geom_freqpoly(). What are the pros and cons of each method?

ggplot(data = diamonds) +
  geom_violin(mapping = aes(x = price, y = cut))

ggplot(data = diamonds) +
  geom_freqpoly(mapping = aes(x = price, y = ..density.., color = cut))

ggplot(data = diamonds) +
  geom_histogram(mapping = aes(x = price)) +
  facet_wrap(~cut, nrow = 5, ncol = 1)

The facetted histogram is the best at comparing the raw counts of the data. By this we mean it is very easy to see that “ideal” diamonds are the most common within the dataset. The frequency polygon is a little difficult to interpret, especially across categories that are similar; in this example premium and very good diamonds. It is however very useful in determining which cut has the highest density at a given price and tracking price change. It is for example, very easy to see how all prices converge as number of diamonds decreases. And the violin plot is ideal at seeing where there are sort of anomalies in the data. For example where there are “bulges” in the violin.

6. If you have a small dataset, it’s sometimes useful to use geom_jitter() to see the relationship between a continuous and categorical variable. The ggbeeswarm package provides a number of methods similar to geom_jitter(). List them and briefly describe what each one does.

library(ggbeeswarm)
?ggbeeswarm

The first two geom_beeswarm and geom_quasirandom produce plots that are a mix of point and violin. Functionally, the violins are produced using points and the difference between the two methods is whether to randomizing points within or across categories. Within each geom you can then specify randomization methods.

geom_beeswarm: works similarly to geom_jitter

geom_quasirandom: works similarly to geom_jitter, but it randomizes points within categories to reduce overplotting.

position_beeswarm: violin point-style plots to show overlapping points. x must be discrete

position_quasirandom: violin point-style plots to show overlapping points. x must be discrete

7.5.2 Two categorical variables

To visualise the covariation between categorical variables, you’ll need to count the number of observations for each combination. One way to do that is to rely on the built-in geom_count():

ggplot(data = diamonds) +
  geom_count(mapping = aes(x = cut, y = color))

The size of each circle in the plot displays how many observations occurred at each combination of values. Covariation will appear as a strong correlation between specific x values and specific y values.

Another approach is to compute the count with dplyr:

library(dplyr)
diamonds %>%
  group_by(cut, color) %>%
  summarise(
    cut_by_color = n()
  )

`summarise()` has grouped output by 'cut'. You can override using the `.groups` argument.

OR count

then visualise with geom_tile() and the fill aesthetic:

diamonds %>%
  count(color, cut) %>%
  ggplot() +
  geom_tile(mapping = aes(x = color, y = cut, fill = n))

If the categorical variables are unordered, you might want to use the seriation package to simultaneously reorder the rows and columns in order to more clearly reveal interesting patterns. For larger plots, you might want to try the d3heatmap or heatmaply packages, which create interactive plots.

7.5.2.1 Exercises

1. How could you rescale the count dataset above to more clearly show the distribution of cut within color, or color within cut?

diamonds %>%
  group_by(color) %>%
  count(color, cut) %>%
  mutate(
    prop = n/sum(n)
    ) %>%
  ggplot(mapping = aes(x = color, y = cut)) +
  geom_tile(aes(fill = prop))

diamonds %>%
  group_by(cut) %>%
  count(cut, color) %>%
  mutate(
    prop = n/sum(n)
  ) %>%
  ggplot(mapping = aes(x = color, y = cut)) +
  geom_tile(aes(fill = prop))

2. Use geom_tile() together with dplyr to explore how average flight delays vary by destination and month of year. What makes the plot difficult to read? How could you improve it?

library(nycflights13)
flights %>%
  filter(!is.na(arr_delay), arr_delay > 0) %>%
  group_by(dest, month) %>%
  mutate(avg_delay = mean(arr_delay)) %>%
  ggplot(mapping = aes(x = factor(month), y = reorder(dest, distance))) +
  geom_tile(mapping = aes(fill = avg_delay)) +
  labs(x = "Month", y = "Destination", fill = "Average Delay")

A few things make the plot difficult to read.

1.) There are too many destinations for all of the fields to fit comfortably on either axis. 2.) The month variable is contained within the data as an integer when it really is more of a nominal value with discrete values between 1:12.

We could fix these issues by grouping destinations by state and by converting month to a true discrete nominal variable. The latter issue is easier to fix, we can simply use factor(month).

3. Why is it slightly better to use aes(x = color, y = cut) rather than aes(x = cut, y = color) in he example above?

diamonds %>%
  group_by(color) %>%
  count(color, cut) %>%
  mutate(
    prop = n/sum(n)
    ) %>%
  ggplot(mapping = aes(x = color, y = cut)) +
  geom_tile(aes(fill = prop))

diamonds %>%
  group_by(color) %>%
  count(color, cut) %>%
  mutate(
    prop = n/sum(n)
    ) %>%
  ggplot(mapping = aes(x = cut, y = color)) +
  geom_tile(aes(fill = prop))

It is generally good practice to put the variable with more categories on the x-axis (horizontal) because it is easier to read. In addition, the flow of data is easier to understand when color is on the x-axis. Essentially, the graphic “waterfalls”. Larger quantities are at the top and the graph has a natural decreasing trend as it flows towards the lower-right quadrant. Inversely, when “cut” is on the x-axis, the understanding of the graph changes to mean that as we approach the bottom-right hand quadrant of the graph the values are increasing - which is somewhat more complicated to conceptualize.

7.5.3 Two continuous variables

You’ve already seen one great way to visualise the covariation between two continuous variables: draw a scatterplot with geom_point(). You can see covariation as a pattern in the points. For example, you can see an exponential relationship between the carat size and price of a diamond.

ggplot(data = diamonds, mapping = aes(x = carat, y = price)) +
  geom_point()

Scatterplots become less useful as the size of your dataset grows, because points begin to overplot, and pile up into areas of uniform black (as above). You’ve already seen one way to fix the problem: using the alpha aesthetic to add transparency.

ggplot(data = diamonds) +
  geom_point(mapping = aes(x = carat, y = price), alpha = 1/100)

But using transparency can be challenging for very large datasets. Another solution is to use bin. Previously you used geom_histogram() and geom_freqpoly() to bin in one dimension. Now you’ll learn how to use geom_bin2d() and geom_hex() to bin in two dimensions.

geom_bin2d() and geom_hex() divide the coordinate plane into 2d bins and then use a fill color to display how many points fall into each bin. geom_bin2d() creates rectangular bins. geom_hex() creates hexagonal bins. You will need to install the hexbin package to use geom_hex().

library(hexbin)
ggplot(data = smaller) +
  geom_bin2d(mapping = aes(x = carat, y = price))

ggplot(data = smaller) +
  geom_hex(mapping = aes(x = carat, y = price))

Another option is to bin one continuous variable so it acts like a categorical variable. Then you can use one of the techniques for visualising the combination of a categorical and a continuous variable that you learned about. For example, you could bin carat and then for each group, displaying a boxplot:

ggplot(data = smaller, mapping= aes(x = carat, y = price)) +
  geom_boxplot(mapping = aes(group = cut_width(carat, 0.1)))

cut_width(x, width), as used above, divides x into bins of width width. By default, boxplots look roughly the same (apart from number of outliers( regardless of how many observations there are, so it’s difficult to tell that each boxplot summarises a different number of points. One way to show that is to make the width of the boxplot proportional to the number of points with varwidth = TRUE.))

ggplot(data = smaller, mapping= aes(x = carat, y = price)) +
  geom_boxplot(mapping = aes(group = cut_width(carat, 0.1)))

ggplot(data = smaller, mapping= aes(x = carat, y = price)) +
  geom_boxplot(mapping = aes(group = cut_width(carat, 0.1)), varwidth = TRUE)

Another approach is to display approximately the same number of points in each bin. That’s the job of cut_number():

ggplot(data = smaller, mapping = aes(x = carat, y = price)) +
  geom_boxplot(mapping = aes(group = cut_number(carat, 20)))

7.5.3.1 Exercises

1. Instead of summarising the conditional distribution with a boxplot, you could use a frequency polygon. What do you need to consider when using cut_width() vs cut_number()? How does that impact a visualisation of the 2d distribution of carat and price?

Both cut_width() and cut_number() split a variable into groups, but bins the values differently. cut_width bins the data into x number of bins all of about the same width if varwidth = TRUE is not specified. cut_number() divides the data into x number of bins where each bin has the same number of values within it.

ggplot(data = smaller, mapping = aes(x = price)) +
  geom_freqpoly(mapping = aes(color = cut_width(carat, .5)))

ggplot(data = smaller, mapping = aes(x = price)) +
  geom_freqpoly(mapping = aes(color = cut_number(carat, n = 5)))

When looking at graphs produced using these it is important to pay attention to the bins. Both will capture all of the data, but cut_width is better for telling you apprx. how many under a curve whereas cut_number is better at giving you an idea of how many within a range are necessary to form the curve.

Also, when using color remember that it is best practice to have no more than 8 categories because the colors become increasingly less distinct afterward.

2. Visualise the distribution of carat, partitioned by price.

ggplot(data = smaller) +
  geom_bin2d(mapping = aes(x = carat, y = price))

ggplot(data = diamonds, mapping = aes(x = carat, y = price)) +
  geom_boxplot(aes(y = cut_number(price, 10))) +
  labs(x = "carat", y = "price")

ggplot(data = diamonds, mapping = aes(x = cut_width(price, 2000, boundary = 0), y = carat)) +
  geom_boxplot(varwidth = TRUE) +
  coord_flip() +
  labs(x = "carat", y = "price")

3. How does the price distribution of very large diamonds compare to small diamonds Is it as you expect, or does it surprise you?

ggplot(data = diamonds, mapping = aes(x = cut_width(carat, 1, boundary = 0), y = price)) +
  geom_boxplot(varwidth = TRUE) +
  labs(x = "carat", y = "price")

ggplot(data = diamonds, mapping = aes(x = cut_interval(carat, n = 3), y = price)) +
  geom_boxplot(varwidth = TRUE) +
  labs(x = "carat", y = "price")

ggplot(data = diamonds, mapping = aes(x = cut_number(carat, 5), y = price)) +
  geom_boxplot() +
  labs(x = "carat", y = "price")

The first thing to note is that there are are far more “smaller” carats (carat sizes < 3) than there are “larger” carats (carat sizes > 3). In fact, only 40 of the diamonds in our dataset containing 53,940 observations have a carat size larger than 3. In fact, about 68% of diamonds in our dataset have sizes smaller than or equal to 1.

This is most evident by the boxplot using cut_number(). As you can see the smallest boxplot is able to achieve the same number of data points using only a range of [.2, .35] whereas the last box requires a range of [1.13, 5.01] to achieve the same number of diamonds.

There is a clear positive relationship between carat size and price. As either increases, so too does the other. There are several other factors to consider other than carat size, but what is perhaps most interesting is that the most significant increases occur when increasing a whole number in carat size (0:1, 1:2, 2:3) but the effect appears to sort of taper off after 3. This effect might reverse with more data points, but it is odd that there is not the price increase we would expect. Perhaps these “larger” carat diamonds are of “inferior” cut or color? Let’s check.

diamonds %>%
  filter(carat >= 3) %>%
  select(carat, cut, color, price) %>%
  ggplot(mapping = aes(x = carat, y = price)) +
  geom_point()

diamonds %>%
  filter(carat >= 3) %>%
  select(carat, cut, color, price) %>%
  group_by(cut, color) %>%
  mutate(
    count = n()
  ) %>%
  ungroup() %>%
  mutate(
    prop = count/sum(count)
  ) %>%
  ggplot(mapping = aes(x = cut, y = color)) +
  geom_count()

diamonds %>%
  filter(carat >= 3) %>%
  select(carat, cut, color, price) %>%
  group_by(cut, color) %>%
  mutate(
    count = n()
  ) %>%
  ungroup() %>%
  mutate(
    prop = count/sum(count)
  ) %>%
  ggplot(mapping = aes(x = carat, y = cut_number(price, 5))) +
  geom_boxplot() +
  labs(x = "carat", y = "price")

4. Combine two of the techniques you’ve learned to visualise the combined distribution of cut, carat, and price.

ggplot(data = diamonds, mapping = aes(x = carat, y = price)) +
  geom_hex() +
  facet_wrap(~cut)

ggplot(data = diamonds, mapping = aes(x = cut_number(carat, 5), y = price)) +
  geom_boxplot(aes(color = cut))

ggplot(data = diamonds, mapping = aes(x = cut, y = price)) +
  geom_boxplot(aes(color = cut_number(carat, 5)))

5. Two dimensional plots reveal outliers that are not visible in one dimensional plots. For example, some points in the plot below have an unusual combination of x and y values, which makes the points outliers even though their x and y values appear normal when examined separately.

ggplot(data = diamonds) +
  geom_point(mapping = aes(x = x, y = y)) +
  coord_cartesian(xlim = c(4,11), ylim = c(4,11))

Why is a scatterplot a better display than a binned plot for this case?

A scatterplot is a better display for this case because there is a strong relationship between x and y. Were we to bin the plots, values belonging to x may mistakenly be categorized as outliers even though they fit strongly with the bivariate relationship.

7.6 Patterns and models

Patterns in your data provide clues about relationships. If a systematic relationship exists between two variables it will appear as a pattern in the data. If you spot a pattern, ask yourself:

• Could this pattern be due to coincidence (i.e. random chance)?
• How can you describe the relationship implied by the pattern?
• How strong is the relationship implied by the pattern?
• What other variables might affect the relationship?
• Does the relationship change if you look at individual subgroups of the data?

A scatterplot of Old Faithful eruption lengths versus the wait time between eruptions shows a pattern: longer wait times are associated with longer eruptions. The scatterplot also displays the two clusters that we noticed above.

ggplot(data = faithful) +
  geom_point(mapping = aes(x = eruptions, y = waiting))

Patterns provide one of the most useful tools for data scientists because they reveal covariation. If you think of variation as a phenomenon that creates uncertainty, covariation is a phenomenon that reduces it. If two variables covary, you can use the values of one variable to make better predictions about the values of the second. If the covariation is due to a causal relationship (a special case), then you can use the value of one variable to control the value of the second.

Models are a tool for extracting patterns out of data. For example, consider the diamonds data. It’s hard to understand the relationship between cut and price, because cut and carat, and carat and price are tightly related. It’s possible to use a model to remove the very strong relationship between price and carat so we can explore the subtleties that remain. The following code fits a model that predicts price from carat and then computes the residuals (the difference between the predicted value and the actual value). The residuals give us a view of the price of the diamond, once the effect of carat has been recmoved.

library(modelr)

mod <- lm(log(price) ~ log(carat), data = diamonds)

diamonds2 <- diamonds %>%
  add_residuals(mod) %>%
  mutate(resid = exp(resid))

ggplot(data = diamonds2) +
  geom_point(mapping = aes(x = carat, y = resid))

Once you’ve removed the strong relationship between carat and price, you can see what you expect in the relationship between cut and price: relative to their size, better quality diamonds are more expensive.

ggplot(data = diamonds2) +
  geom_boxplot(mapping = aes(x = cut, y = resid))

7.7 ggplot2 calls

As we move on from these introductory chapters, we’ll transition to a more concise expression of ggplot2 code. So far we’ve been very explicit, which is helpul when you are learning:

ggplot(data = faithful, mapping = aes(x = eruptions)) +
  geom_freqpoly(binwidth = .25)

Typically, the first one or two arguments to a function are so important that you should know them by heart. The first two arguments to ggplot() are data and mapping, and the first two arguments to aes() are x and y. In the remainder of the book, we won’t supply those names. That saves typing, and, by reducing the amount of boilerplate, makes it easier to see what’s different between plots. That’s a really important programming concern that we’ll come back in functions.

Rewriting the previous plot more concisely yields:

ggplot(faithful, aes(eruptions)) +
  geom_freqpoly(binwidth = .25)

Sometimes we’ll trun the end of a pipeline of ata transformation into a plot. Watch for the transition from %>% to +. I wish this transition wasn’t necessary but unfortunately ggplot2 was created before the pipe was discovered.

diamonds %>%
  count(cut, clarity) %>%
  ggplot(aes(clarity, cut, fill = n)) +
  geom_tile()