22  ggplot2

The first graphs in Chapter 11 helped us explore relationships, distributions, and category shares. A figure for a research paper or presentation also needs deliberate decisions about its analysis sample, group distinctions, panel scales, labels, and saved dimensions. We now use the table and variable operations learned since that chapter to prepare those comparisons explicitly.

By the end of this chapter, you will be able to:

This chapter builds on the plotting grammar from Chapter 11, selecting observations, creating columns, and factors from Chapter 15. We use those tools to construct figures in stages. We return to the country data from the early manipulation chapters, then move to annual country records when the question requires a time series.

Load tidyverse for ggplot2 and the table operations, and tidystata for lifeexp. Copy the captured OWID file into your project’s examples folder. read_csv() restores that public table as owid. No GSS file is needed for this chapter.

library(tidyverse)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.2.1     ✔ readr     2.2.0
✔ forcats   1.0.1     ✔ stringr   1.6.0
✔ ggplot2   4.0.3     ✔ tibble    3.3.1
✔ lubridate 1.9.5     ✔ tidyr     1.3.2
✔ purrr     1.2.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(tidystata)

Attaching package: 'tidystata'

The following object is masked from 'package:dplyr':

    recode
owid <- read_csv("examples/life-expectancy.csv")
Rows: 21565 Columns: 4
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): Entity, Code
dbl (2): Year, Life expectancy

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

22.1 A relationship between countries

Do countries with higher income tend to have higher life expectancy? In lifeexp, gnppc records GNI per capita using the Atlas method, in current US dollars in 1998. lexp records life expectancy at birth in years in 1998. Each point will represent a country, so both measurements must be recorded to locate it. Check the familiar summary’s n column before plotting:

nrow(lifeexp)
[1] 68
summ(lifeexp, gnppc, lexp)
# A tibble: 2 × 8
  name  type      n unique   min   mean       sd   max
  <chr> <chr> <int>  <int> <dbl>  <dbl>    <dbl> <dbl>
1 gnppc dbl      67     64   170 9137.  11497.   46650
2 lexp  dbl      68     20    57   71.7     5.16    79

All 68 countries have life expectancy recorded, but only 67 have income recorded. Therefore, 67 countries can contribute points. The missing income is not zero.

ggplot() receives a data frame and returns a plot object. Inside it, aes() defines mappings: x = gnppc positions countries horizontally by income, and y = lexp positions them vertically by life expectancy. geom_point() adds a point for each country with both coordinates. The + operator adds a plotting layer; put it at the end of the line when continuing on the next line.

Set na.rm = TRUE in the layer to omit countries that cannot be positioned because a coordinate is missing. This choice avoids a warning about that explained omission; it neither replaces the missing value nor removes a row from lifeexp.

country_plot <- ggplot(lifeexp, aes(x = gnppc, y = lexp)) +
  geom_point(color = "#2B6487", na.rm = TRUE)

Assignment stores the plot without displaying it. The hexadecimal value specifies the book’s blue. Because color is outside aes(), it sets one appearance for every point. Only the variables inside aes() determine positions here.

labs() labelled axes in Chapter 11. Its additional subtitle argument puts a short line of context above the panel. Here it identifies the complete-pair sample. A subtitle belongs inside the plotted image; the automatic figure caption below it also records units and context when the figure is read in the book. Add the labels to the saved object, assign the updated result, then type its name to display it.

country_plot <- country_plot +
  labs(
    x = "GNI per capita (current US dollars, 1998)",
    y = "Life expectancy at birth (years, 1998)",
    subtitle = "Countries with both measurements recorded"
  )
country_plot
Scatterplot of country GNI per capita horizontally and life expectancy in years vertically. Life expectancy is generally higher among higher-income countries, with variation among countries at similar incomes.
Figure 22.1: GNI per capita and life expectancy in 1998. Each point is a country with both measurements; GNI uses the Atlas method in current US dollars.

Figure 22.1 suggests higher life expectancy at higher incomes, with substantial variation at the lower end of the income distribution. A horizontal distance measures a dollar difference; a vertical difference measures years. This is an association among countries in the dataset. It neither estimates an individual income effect nor shows that changing national income by itself causes a particular change in longevity.

22.2 Add regional context

Do countries with higher income tend to have higher life expectancy within the recorded regions? Return to the country relationship in Figure 22.1. gnppc records GNI per capita using the Atlas method, in current US dollars in 1998; lexp records life expectancy at birth in years in 1998. This time we prepare the complete country pairs in an object before plotting and distinguish regions by point shape.

The familiar filter() retains countries with both measurements: is.na() identifies a missing value, and ! reverses that result. The comma requires both conditions to hold. mutate() then converts the labelled region column into categories with readable labels. haven::as_factor() performs that conversion; the package prefix identifies the function without loading another package. The pipe passes each result to the following operation.

country_points <- lifeexp |>
  filter(!is.na(gnppc), !is.na(lexp)) |>
  mutate(region = haven::as_factor(region))
nrow(lifeexp)
[1] 68
nrow(country_points)
[1] 67

The graph uses 67 of 68 countries, the same complete pairs as the preceding scatterplot. Saving that selection makes the analysis sample available for checking and reuse. It neither replaces the missing measurement nor changes lifeexp.

22.2.1 Map a variable, set an appearance

Map shape = region inside aes() to distinguish the recorded regions. The mapping uses each country’s category and produces a legend. All points retain the same blue because shape and the labelled legend carry the regional distinction. As in the preceding plot, the color argument outside aes() sets one appearance for every point.

region_plot <- ggplot(
  country_points,
  aes(x = gnppc, y = lexp, shape = region)
) +
  geom_point(color = "#2B6487")

Assignment stores this draft plot. We can refine it by adding layers to its name rather than rebuilding the data and mappings. labs() labels the axes, the shape legend, and a subtitle. theme_bw() supplies a white background and visible panel structure; it changes presentation without changing observations or calculations. Assign the updated result and print it.

region_plot <- region_plot +
  labs(
    x = "GNI per capita (current US dollars, 1998)",
    y = "Life expectancy at birth (years, 1998)",
    shape = "Region",
    subtitle = "Countries with both measurements recorded"
  ) +
  theme_bw()
region_plot
Country scatterplot with GNI per capita horizontally and life expectancy in years vertically. Life expectancy is generally higher among higher-income countries, with variation among countries at similar incomes. Point shapes and a legend identify regions.
Figure 22.2: GNI per capita and life expectancy in 1998. Each point is a country with both measurements; GNI uses the Atlas method in current US dollars. Shapes distinguish recorded regions.

Figure 22.2 describes an association among countries in this dataset. The points suggest higher life expectancy at higher incomes, with substantial variation at the lower end of the income distribution. A point’s horizontal distance measures a dollar difference on this scale; a vertical difference measures years. The graph neither estimates an individual income effect nor shows that changing national income by itself causes a particular change in longevity.

A quoted constant inside aes(), such as color = "blue", has a different meaning from the same argument outside it. Inside, ggplot2 treats "blue" as a category label and chooses a display color through a scale, usually with a legend. To request blue points, set color outside aes() as above. To distinguish regions by an appearance, map the unquoted variable region.

NoteOptional: compare income ratios

The linear income axis gives equal space to equal dollar differences. Much of the variation among lower-income countries occupies a small part of that axis. A logarithmic income scale offers another view of these same observations, with spacing determined by ratios instead of dollar differences.

A logarithm requires positive inputs. The recorded incomes among these complete country pairs are finite; check their minimum before changing the scale. min() returns the smallest value:

min(country_points$gnppc)
[1] 170

The minimum is 170 dollars, above zero. No further country exclusions are needed. scale_x_log10() replaces the horizontal scale with a base-ten logarithmic scale; it returns a scale to add to the existing plot object.

log_income_plot <- region_plot +
  scale_x_log10()
log_income_plot
Country scatterplot with income on a logarithmic horizontal axis and life expectancy in years vertically. Regional point shapes are unchanged. Lower-income countries are more spread out than in the linear view.
Figure 22.3: GNI per capita and life expectancy in 1998, with a logarithmic income axis. The same complete country pairs and regional shapes are retained; income labels remain current US dollars.

The axis labels still denote dollars. Equal horizontal distances now represent equal income ratios: 100 to 1,000 dollars has the same spacing as 1,000 to 10,000 dollars. A tenfold increase is one unit on the transformed axis, although the labels display the original currency values. Life expectancy remains on its original scale in years. The new view changes spacing, not recorded incomes, the 67 contributing countries, or the limits of a country-level association. The original region_plot remains available.

To inspect countries below 10,000 dollars more closely, change the display window. coord_cartesian() sets the visible coordinates; xlim receives the lower and upper horizontal limits in a vector.

zoomed_income_plot <- region_plot +
  coord_cartesian(xlim = c(0, 10000))
zoomed_income_plot
The regional country scatterplot zoomed to incomes from zero to ten thousand dollars. Higher-income points lie outside the displayed window; the original complete-pair plot data are retained.
Figure 22.4: The country income comparison zoomed to 0–10,000 current US dollars in 1998. Countries beyond this horizontal window are cropped; the plot retains all complete country pairs.

The plot still contains 67 country pairs, with 19 beyond the window. A cropped point remains part of the plot’s input even though it is hidden. Compare selecting rows:

income_window <- country_points |>
  filter(gnppc <= 10000)
nrow(country_points)
[1] 67
nrow(income_window)
[1] 48

Filtering returns a smaller table of 48 countries and changes the inputs to subsequent calculations. Zooming only changes what is displayed. Neither command replaces country_points; the full region_plot remains available.

22.4 Save the figure

ggsave() writes a plot to a file. The filename extension selects the image format; here .png requests a raster image. Specify plot = trend_plot so the saved figure is independent of whichever plot happened to be displayed most recently. width, height, and units set its physical dimensions, and dpi specifies the raster resolution in dots per inch.

ggsave(
  "life-expectancy-trends.png",
  plot = trend_plot,
  width = 8, height = 5, units = "in", dpi = 300
)

The file is saved in the current project folder at eight by five inches and 300 dpi. Open the exported image and check that labels, panel headings, and units are readable at the size where it will appear. Changing output dimensions can change the space available for text even when the plotted values are unchanged. Keep the code and source data with the figure so you can reproduce or revise it.

22.5 Review

  • Prepare a figure’s analysis sample explicitly and retain the source data unchanged.
  • A mapping inside aes() varies by values and can produce a legend; a constant outside it sets one appearance.
  • Add labels and a theme to a saved plot to refine its presentation without changing its observations.
  • geom_line() connects ordered observations, while facet_wrap() places groups in labelled panels.
  • Check annual coverage before drawing lines, and use common scales when comparing levels.
  • coord_cartesian() changes the display window while retaining plot inputs; filtering changes the data supplied to later calculations.
  • ggsave() exports a named plot with explicit dimensions and resolution.
  • Optional: scale_x_log10() changes income spacing from equal differences to equal ratios; check positive support first.

The next chapter, Chapter 23, uses explicit statistical comparisons to ask how sample differences relate to uncertainty. The graph remains useful for checking what those comparisons summarize.

22.6 Exercises

  1. A researcher writes geom_point(aes(color = "blue")) expecting every point to use that blue. Explain the result and create a corrected plot object using the complete country pairs. Then create a region-shape version. Keep both objects without changing lifeexp.
  2. Replace China with Canada in the annual trend comparison. Retain 2000–2023, use common-scale panels, and save the explicit plot object as canada-us-trends.png, at eight by five inches and 300 dpi. Check that both entities have the requested annual coverage before interpreting the lines.

Worked answers

Start a fresh session. These assignments create plot objects; type an object’s name in the console to display it.

library(tidyverse)
library(tidystata)
answer_countries <- lifeexp |>
  filter(!is.na(gnppc), !is.na(lexp)) |>
  mutate(region = haven::as_factor(region))
answer_blue_plot <- ggplot(answer_countries, aes(x = gnppc, y = lexp)) +
  geom_point(color = "#2B6487") +
  labs(x = "GNI per capita (current US dollars, 1998)",
       y = "Life expectancy at birth (years, 1998)") +
  theme_bw()
answer_region_plot <- ggplot(
  answer_countries, aes(x = gnppc, y = lexp, shape = region)
) +
  geom_point(color = "#2B6487") +
  labs(x = "GNI per capita (current US dollars, 1998)",
       y = "Life expectancy at birth (years, 1998)", shape = "Region") +
  theme_bw()

Inside aes(), the quoted word would label one mapped category and could produce a legend; it would not directly set the requested color. Outside aes(), the color applies to every point. The shape mapping varies by the recorded region. Both objects use the same 67 complete country pairs.

Start from the same saved public CSV in a fresh session. count(Entity) counts retained annual records within each entity, providing a coverage check for the chosen window.

library(tidyverse)
answer_owid <- read_csv("examples/life-expectancy.csv")
Rows: 21565 Columns: 4
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): Entity, Code
dbl (2): Year, Life expectancy

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

answer_trends <- answer_owid |>
  filter(
    Entity %in% c("Canada", "United States"),
    Year >= 2000, Year <= 2023,
    !is.na(`Life expectancy`)
  ) |>
  rename(lexp = `Life expectancy`) |>
  arrange(Entity, Year)
answer_coverage <- answer_trends |>
  count(Entity)
answer_coverage
# A tibble: 2 × 2
  Entity            n
  <chr>         <int>
1 Canada           24
2 United States    24

answer_trend_plot <- ggplot(answer_trends, aes(x = Year, y = lexp)) +
  geom_line(color = "#2B6487") +
  facet_wrap(~ Entity) +
  labs(x = "Year", y = "Life expectancy at birth (years)") +
  theme_bw()
ggsave(
  "canada-us-trends.png", plot = answer_trend_plot,
  width = 8, height = 5, units = "in", dpi = 300
)

The coverage output shows 24 annual records for Canada and 24 for the United States, covering the requested window. The shared axes permit direct comparison of levels. The exported file contains the named comparison even if another graph was displayed afterward; its annual lines do not imply exact within-year dates.