Prerequisites

This chapter focuses on ggplot2, one of the core packages in the tidyverse. To access the datasets, help pages, and functions used in this chapter, load the tidyverse by running this code:

library(tidyverse)
#> ── Attaching core tidyverse packages ──────────────── tidyverse 1.3.2.9000 ──
#> ✔ dplyr     1.0.99.9000     ✔ readr     2.1.3      
#> ✔ forcats   0.5.2.9000      ✔ stringr   1.5.0.9000 
#> ✔ ggplot2   3.4.0.9000      ✔ tibble    3.1.8      
#> ✔ lubridate 1.9.0           ✔ tidyr     1.2.1.9001 
#> ✔ purrr     1.0.1           
#> ── Conflicts ─────────────────────────────────────── tidyverse_conflicts() ──
#> ✖ dplyr::filter() masks stats::filter()
#> ✖ dplyr::lag()    masks stats::lag()
#> ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors

That one line of code loads the core tidyverse; packages which you will use in almost every data analysis. It also tells you which functions from the tidyverse conflict with functions in base R (or from other packages you might have loaded).

If you run this code and get the error message there is no package called 'tidyverse', you’ll need to first install it, then run library() once again.

install.packages("tidyverse")
library(tidyverse)

You only need to install a package once, but you need to reload it every time you start a new session.

In addition to tidyverse, we will also use the palmerpenguins package, which includes the penguins dataset containing body measurements for penguins on three islands in the Palmer Archipelago.

library(palmerpenguins)

Thepenguins data frame

You can test your answer with the penguins data frame found in palmerpenguins (a.k.a. palmerpenguins::penguins). A data frame is a rectangular collection of variables (in the columns) and observations (in the rows). penguins contains 344 observations collected and made available by Dr. Kristen Gorman and the Palmer Station, Antarctica LTERHorst AM, Hill AP, Gorman KB (2020). palmerpenguins: Palmer Archipelago (Antarctica) penguin data. R package version 0.1.0. https://allisonhorst.github.io/palmerpenguins/. doi: 10.5281/zenodo.3960218..

penguins
#> # A tibble: 344 × 8
#>   species island   bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
#>   <fct>   <fct>             <dbl>         <dbl>             <int>       <int>
#> 1 Adelie  Torgers…           39.1          18.7               181        3750
#> 2 Adelie  Torgers…           39.5          17.4               186        3800
#> 3 Adelie  Torgers…           40.3          18                 195        3250
#> 4 Adelie  Torgers…           NA            NA                  NA          NA
#> 5 Adelie  Torgers…           36.7          19.3               193        3450
#> 6 Adelie  Torgers…           39.3          20.6               190        3650
#> # … with 338 more rows, and 2 more variables: sex <fct>, year <int>

This data frame contains 8 columns. For an alternative view, where you can see all variables and the first few observations of each variable, use glimpse(). Or, if you’re in RStudio, run View(penguins) to open an interactive data viewer.

glimpse(penguins)
#> Rows: 344
#> Columns: 8
#> $ species           <fct> Adelie, Adelie, Adelie, Adelie, Adelie, Adelie, A…
#> $ island            <fct> Torgersen, Torgersen, Torgersen, Torgersen, Torge…
#> $ bill_length_mm    <dbl> 39.1, 39.5, 40.3, NA, 36.7, 39.3, 38.9, 39.2, 34.…
#> $ bill_depth_mm     <dbl> 18.7, 17.4, 18.0, NA, 19.3, 20.6, 17.8, 19.6, 18.…
#> $ flipper_length_mm <int> 181, 186, 195, NA, 193, 190, 181, 195, 193, 190, …
#> $ body_mass_g       <int> 3750, 3800, 3250, NA, 3450, 3650, 3625, 4675, 347…
#> $ sex               <fct> male, female, female, NA, female, male, female, m…
#> $ year              <int> 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2…

Among the variables in penguins are:

1. species: a penguin’s species (Adelie, Chinstrap, or Gentoo).

2. flipper_length_mm: length of a penguin’s flipper, in millimeters.

3. body_mass_g: body mass of a penguin, in grams.

To learn more about penguins, open its help page by running ?penguins.

Ultimate goal

Our ultimate goal in this chapter is to recreate the following visualization displaying the relationship between flipper lengths and body masses of these penguins, taking into consideration the species of the penguin.

Creating a ggplot

Let’s recreate this plot layer-by-layer.

With ggplot2, you begin a plot with the function ggplot(), defining a plot object that you then add layers to. The first argument of ggplot() is the dataset to use in the graph and So ggplot(data = penguins) creates an empty graph. This is not a very exciting plot, but you can think of it like an empty canvas you’ll paint the remaining layers of your plot onto.

ggplot(data = penguins)

Next, we need to tell ggplot() the variables from this data frame that we want to map to visual properties (aesthetics) of the plot. The mapping argument of the ggplot() function defines how variables in your dataset are mapped to visual properties of your plot. The mapping argument is always paired with the aes() function, and the x and y arguments of aes() specify which variables to map to the x and y axes. For now, we will only map flipper length to the x aesthetic and body mass to the y aesthetic. ggplot2 looks for the mapped variables in the data argument, in this case, penguins.

The following plots show the result of adding these mappings, one at a time.

ggplot(
  data = penguins,
  mapping = aes(x = flipper_length_mm)
)
ggplot(
  data = penguins,
  mapping = aes(x = flipper_length_mm, y = body_mass_g)
)

Our empty canvas now has more structure – it’s clear where flipper lengths will be displayed (on the x-axis) and where body masses will be displayed (on the y-axis). But the penguins themselves are not yet on the plot. This is because we have not yet articulated, in our code, how to represent the observations from our data frame on our plot.

To do so, we need to define a geom: the geometrical object that a plot uses to represent data. These geometric objects are made available in ggplot2 with functions that start with geom_. People often describe plots by the type of geom that the plot uses. For example, bar charts use bar geoms (geom_bar()), line charts use line geoms (geom_line()), boxplots use boxplot geoms (geom_boxplot()), and so on. Scatterplots break the trend; they use the point geom: geom_point().

The function geom_point() adds a layer of points to your plot, which creates a scatterplot. ggplot2 comes with many geom functions that each adds a different type of layer to a plot. You’ll learn a whole bunch of geoms throughout the book, particularly in #chp-layers.

ggplot(
  data = penguins,
  mapping = aes(x = flipper_length_mm, y = body_mass_g)
) +
  geom_point()
#> Warning: Removed 2 rows containing missing values (`geom_point()`).

Now we have something that looks like what we might think of as a “scatter plot”. It doesn’t yet match our “ultimate goal” plot, but using this plot we can start answering the question that motivated our exploration: “What does the relationship between flipper length and body mass look like?” The relationship appears to be positive, fairly linear, and moderately strong. Penguins with longer flippers are generally larger in terms of their body mass.

Before we add more layers to this plot, let’s pause for a moment and review the warning message we got:

Removed 2 rows containing missing values (geom_point()).

We’re seeing this message because there are two penguins in our dataset with missing body mass and flipper length values and ggplot2 has no way of representing them on the plot. You don’t need to worry about understanding the following code yet (you will learn about it in #chp-data-transform), but it’s one way of identifying the observations with NAs for either body mass or flipper length.

penguins |>
  select(species, flipper_length_mm, body_mass_g) |>
  filter(is.na(body_mass_g) | is.na(flipper_length_mm))
#> # A tibble: 2 × 3
#>   species flipper_length_mm body_mass_g
#>   <fct>               <int>       <int>
#> 1 Adelie                 NA          NA
#> 2 Gentoo                 NA          NA

Like R, ggplot2 subscribes to the philosophy that missing values should never silently go missing. This type of warning is probably one of the most common types of warnings you will see when working with real data – missing values are a very common issue and you’ll learn more about them throughout the book, particularly in #chp-missing-values. For the remaining plots in this chapter we will suppress this warning so it’s not printed alongside every single plot we make.

Adding aesthetics and layers

Scatterplots are useful for displaying the relationship between two variables, but it’s always a good idea to be skeptical of any apparent relationship between two variables and ask if there may be other variables that explain or change the nature of this apparent relationship. Let’s incorporate species into our plot and see if this reveals any additional insights into the apparent relationship between flipper length and body mass. We will do this by representing species with different colored points.

To achieve this, where should species go in the ggplot call from earlier? If you guessed “in the aesthetic mapping, inside of aes()”, you’re already getting the hang of creating data visualizations with ggplot2! And if not, don’t worry. Throughout the book you will make many more ggplots and have many more opportunities to check your intuition as you make them.

ggplot(
  data = penguins,
  mapping = aes(x = flipper_length_mm, y = body_mass_g, color = species)
) +
  geom_point()

When a variable is mapped to an aesthetic, ggplot2 will automatically assign a unique value of the aesthetic (here a unique color) to each unique level of the variable (each of the three species), a process known as scaling. ggplot2 will also add a legend that explains which values correspond to which levels.

Now let’s add one more layer: a smooth curve displaying the relationship between body mass and flipper length. Before you proceed, refer back to the code above, and think about how we can add this to our existing plot.

Since this is a new geometric object representing our data, we will add a new geom: geom_smooth().

ggplot(
  data = penguins,
  mapping = aes(x = flipper_length_mm, y = body_mass_g, color = species)
) +
  geom_point() +
  geom_smooth()

We have successfully added smooth curves, but this plot doesn’t look like the plot from #sec-ultimate-goal, which only has one curve for the entire dataset as opposed to separate curves for each of the penguin species.

When aesthetic mappings are defined in ggplot(), at the global level, they’re inherited by each of the subsequent geom layers of the plot. However, each geom function in ggplot2 can also take a mapping argument, which allows for aesthetic mappings at the local level. Since we want points to be colored based on species but don’t want the smooth curves to be separated out for them, we should specify color = species for geom_point() only.

ggplot(
  data = penguins,
  mapping = aes(x = flipper_length_mm, y = body_mass_g)
) +
  geom_point(mapping = aes(color = species)) +
  geom_smooth()

Voila! We have something that looks very much like our ultimate goal, though it’s not yet perfect. We still need to use different shapes for each species of penguins and improve labels.

It’s generally not a good idea to represent information using only colors on a plot, as people perceive colors differently due to color blindness or other color vision differences. Therefore, in addition to color, we can also map species to the shape aesthetic.

ggplot(
  data = penguins,
  mapping = aes(x = flipper_length_mm, y = body_mass_g)
) +
  geom_point(mapping = aes(color = species, shape = species)) +
  geom_smooth()

Note that the legend is automatically updated to reflect the different shapes of the points as well.

And finally, we can improve the labels of our plot using the labs() function in a new layer. Some of the arguments to labs() might be self explanatory: title adds a title and subtitle adds a subtitle to the plot. Other arguments match the aesthetic mappings, x is the x-axis label, y is the y-axis label, and color and shape define the label for the legend.

ggplot(
  data = penguins,
  mapping = aes(x = flipper_length_mm, y = body_mass_g)
) +
  geom_point(aes(color = species, shape = species)) +
  geom_smooth() +
  labs(
    title = "Body mass and flipper length",
    subtitle = "Dimensions for Adelie, Chinstrap, and Gentoo Penguins",
    x = "Flipper length (mm)",
    y = "Body mass (g)",
    color = "Species",
    shape = "Species"
  )

We finally have a plot that perfectly matches our “ultimate goal”!

Exercises

1. How many rows are in penguins? How many columns?

2. What does the bill_depth_mm variable in the penguins data frame describe? Read the help for ?penguins to find out.

3. Make a scatterplot of bill_depth_mm vs. bill_length_mm. Describe the relationship between these two variables.

4. What happens if you make a scatterplot of species vs bill_depth_mm? Why is the plot not useful?

5. Why does the following give an error and how would you fix it?

   ggplot(data = penguins) + 
     geom_point()

6. What does the na.rm argument do in geom_point()? What is the default value of the argument? Create a scatterplot where you successfully use this argument set to TRUE.

7. Add the following caption to the plot you made in the previous exercise: “Data come from the palmerpenguins package.” Hint: Take a look at the documentation for labs().

8. Recreate the following visualization. What aesthetic should bill_depth_mm be mapped to? And should it be mapped at the global level or at the geom level?

   

9. Run this code in your head and predict what the output will look like. Then, run the code in R and check your predictions.

   ggplot(
     data = penguins,
     mapping = aes(x = flipper_length_mm, y = body_mass_g, color = island)
   ) +
     geom_point() +
     geom_smooth(se = FALSE)

10. Will these two graphs look different? Why/why not?

    ggplot(
      data = penguins,
      mapping = aes(x = flipper_length_mm, y = body_mass_g)
    ) +
      geom_point() +
      geom_smooth()
    
    ggplot() +
      geom_point(
        data = penguins,
        mapping = aes(x = flipper_length_mm, y = body_mass_g)
      ) +
      geom_smooth(
        data = penguins,
        mapping = aes(x = flipper_length_mm, y = body_mass_g)
      )

A categorical variable

A variable is categorical if it can only take one of a small set of values. To examine the distribution of a categorical variable, you can use a bar chart. The height of the bars displays how many observations occurred with each x value.

ggplot(penguins, aes(x = species)) +
  geom_bar()

In bar plots of categorical variables with non-ordered levels, like the penguin species above, it’s often preferable to reorder the bars based on their frequencies. Doing so requires transforming the variable to a factor (how R handles categorical data) and then reordering the levels of that factor.

ggplot(penguins, aes(x = fct_infreq(species))) +
  geom_bar()

You will learn more about factors and functions for dealing with factors (like fct_infreq() shown above) in #chp-factors.

A numerical variable

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

ggplot(penguins, aes(x = body_mass_g)) +
  geom_histogram(binwidth = 200)
ggplot(penguins, aes(x = body_mass_g)) +
  geom_density()

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 bar shows that 39 observations have a body_mass_g value between 3,500 and 3,700 grams, which are the left and right edges of the bar.

penguins |>
  count(cut_width(body_mass_g, 200))
#> # A tibble: 19 × 2
#>   `cut_width(body_mass_g, 200)`     n
#>   <fct>                         <int>
#> 1 [2.7e+03,2.9e+03]                 7
#> 2 (2.9e+03,3.1e+03]                10
#> 3 (3.1e+03,3.3e+03]                23
#> 4 (3.3e+03,3.5e+03]                38
#> 5 (3.5e+03,3.7e+03]                39
#> 6 (3.7e+03,3.9e+03]                37
#> # … with 13 more rows

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 patterns. In the plots below a binwidth of 20 is too narrow, resulting in too many bars, making it difficult to determine the shape of the distribution. Similarly, a binwidth of 2,000 is too high, resulting in all data being binned into only three bars, and also making it difficult to determine the shape of the distribution.

ggplot(penguins, aes(x = body_mass_g)) +
  geom_histogram(binwidth = 20)
ggplot(penguins, aes(x = body_mass_g)) +
  geom_histogram(binwidth = 200)
ggplot(penguins, aes(x = body_mass_g)) +
  geom_histogram(binwidth = 2000)

Exercises

1. Make a bar plot of species of penguins, where you assign species to the y aesthetic. How is this plot different?

2. How are the following two plots different? Which aesthetic, color or fill, is more useful for changing the color of bars?

   ggplot(penguins, aes(x = species)) +
     geom_bar(color = "red")
   
   ggplot(penguins, aes(x = species)) +
     geom_bar(fill = "red")

3. What does the bins argument in geom_histogram() do?

4. Make a histogram of the carat variable in the diamonds dataset. Experiment with different binwidths. What binwidth reveals the most interesting patterns?

A numerical and a categorical variable

To visualize the relationship between a numerical and a categorical variable we can use side-by-side box plots. A boxplot is a type of visual shorthand for a distribution of values that is popular among statisticians. As shown in #fig-eda-boxplot, each boxplot consists of:

• A box that stretches from the 25th percentile of the distribution to the 75th percentile, a distance known as the interquartile 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. These 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 body mass by species using geom_boxplot():

ggplot(penguins, aes(x = species, y = body_mass_g)) +
  geom_boxplot()

Alternatively, we can make frequency polygons with geom_freqpoly(). geom_freqpoly() performs the same calculation as geom_histogram(), but instead of displaying the counts with bars, it uses lines instead. It’s much easier to understand overlapping lines than bars of geom_histogram(). There are a few challenges with this type of plot, which we will come back to in #sec-cat-num on exploring the covariation between a categorical and a numerical variable.

ggplot(penguins, aes(x = body_mass_g, color = species)) +
  geom_freqpoly(binwidth = 200, linewidth = 0.75)

We’ve also customized the thickness of the lines using the linewidth argument in order to make them stand out a bit more against the background.

We can also use overlaid density plots, with species mapped to both color and fill aesthetics and using the alpha aesthetic to add transparency to the filled density curves. This aesthetic takes values between 0 (completely transparent) and 1 (completely opaque). In the following plot it’s set to 0.5.

ggplot(penguins, aes(x = body_mass_g, color = species, fill = species)) +
  geom_density(alpha = 0.5)

Note the terminology we have used here:

• We map variables to aesthetics if we want the visual attribute represented by that aesthetic to vary based on the values of that variable.
• Otherwise, we set the value of an aesthetic.

Two categorical variables

We can use segmented bar plots to visualize the distribution between two categorical variables. In creating this bar chart, we map the variable we want to divide the data into first to the x aesthetic and the variable we then further want to divide each group into to the fill aesthetic.

Below are two segmented bar plots, both displaying the relationship between island and species, or specifically, visualizing the distribution of species within each island. The plot on the left shows the frequencies of each species of penguins on each island and the plot on the right shows the relative frequencies (proportions) of each species within each island (despite the incorrectly labeled y-axis that says “count”). The relative frequency plot, created by setting position = "fill" in the geom is more useful for comparing species distributions across islands since it’s not affected by the unequal numbers of penguins across the islands. Based on the plot on the left, we can see that Gentoo penguins all live on Biscoe island and make up roughly 75% of the penguins on that island, Chinstrap all live on Dream island and make up roughly 50% of the penguins on that island, and Adelie live on all three islands and make up all of the penguins on Torgersen.

ggplot(penguins, aes(x = island, fill = species)) +
  geom_bar()
ggplot(penguins, aes(x = island, fill = species)) +
  geom_bar(position = "fill")

Two numerical variables

So far you’ve learned about scatterplots (created with geom_point()) and smooth curves (created with geom_smooth()) for visualizing the relationship between two numerical variables. A scatterplot is probably the most commonly used plot for visualizing the relationship between two variables.

ggplot(penguins, aes(x = flipper_length_mm, y = body_mass_g)) +
  geom_point()

Three or more variables

One way to add additional variables to a plot is by mapping them to an aesthetic. For example, in the following scatterplot the colors of points represent species and the shapes of points represent islands.

ggplot(penguins, aes(x = flipper_length_mm, y = body_mass_g)) +
  geom_point(aes(color = species, shape = island))

However adding too many aesthetic mappings to a plot makes it cluttered and difficult to make sense of. Another way, which is particularly useful for categorical variables, is to split your plot into facets, subplots that each display one subset of the data.

To facet your plot by a single variable, use facet_wrap(). The first argument of facet_wrap() is a formulaHere “formula” is the name of the type of thing created by ~, not a synonym for “equation”., which you create with ~ followed by a variable name. The variable that you pass to facet_wrap() should be categorical.

ggplot(penguins, aes(x = flipper_length_mm, y = body_mass_g)) +
  geom_point(aes(color = species, shape = species)) +
  facet_wrap(~island)

You will learn about many other geoms for visualizing distributions of variables and relationships between them in #chp-layers.

Exercises

1. Which variables in mpg are categorical? Which variables are continuous? (Hint: type ?mpg to read the documentation for the dataset). How can you see this information when you run mpg?

2. Make a scatterplot of hwy vs. displ using the mpg data frame. Next, map a third, numerical variable to color, then size, then both color and size, then shape. How do these aesthetics behave differently for categorical vs. numerical variables?

3. In the scatterplot of hwy vs. displ, what happens if you map a third variable to linewidth?

4. What happens if you map the same variable to multiple aesthetics?

5. Make a scatterplot of bill_depth_mm vs. bill_length_mm and color the points by species. What does adding coloring by species reveal about the relationship between these two variables?

6. Why does the following yield two separate legends? How would you fix it to combine the two legends?

   ggplot(
     data = penguins,
     mapping = aes(
       x = bill_length_mm, y = bill_depth_mm, 
       color = species, shape = species
     )
   ) +
     geom_point() +
     labs(color = "Species")

   

Exercises

1. Run the following lines of code. Which of the two plots is saved as mpg-plot.png? Why?

   ggplot(mpg, aes(x = class)) +
     geom_bar()
   ggplot(mpg, aes(x = cty, y = hwy)) +
     geom_point()
   ggsave("mpg-plot.png")

2. What do you need to change in the code above to save the plot as a PDF instead of a PNG?