Why we are here

In class we built plots from NYC HANES: histograms, box plots, a bar chart, and scatterplots. Every one of those was the right plot, but we never slowed down to ask why it was the right plot.

That is the whole point of this session. Slide 13 from Week 3 asked:

What kind of graph should I use? Depends on what kind of data you have (discrete/categorical vs. continuous) and how many factors you want to compare at once.

We are going to turn that slide into a working decision rule you can apply to any dataset, and then practice it on new data.

The chart chooser

Almost every plot you will make in this course falls somewhere in this table. Find your variables, read across.

What you have	The question it answers	Go-to geom	Also consider
1 continuous	What is the distribution?	`geom_histogram()`	`geom_density()`, `geom_boxplot()`
1 categorical	How many in each group?	`geom_bar()`	`geom_col()` if you already have counts
1 continuous + 1 categorical	Does the distribution differ by group?	`geom_boxplot()`	`geom_violin()`, faceted histogram, overlaid density
2 continuous	Do they move together?	`geom_point()`	add `geom_smooth()`
2 categorical	How does composition differ?	`geom_bar(position = ...)`	`"dodge"`, `"fill"`, `"stack"`
time + 1 continuous	How does it change over time?	`geom_line()`	`geom_point()` for sparse data
3+ variables	All of the above, split	any of the above + `color =`	`facet_wrap()`

Two rules that will save you a lot of grief:

ggplot decides what to draw based on the type of your variable, not what it means to you. A month coded 5 is a number to R even if it is a category to you. This is exactly why we spent class time on set_value_labels() and mutate_if(is.labelled, to_factor) – that wrangling step is what makes your plots come out right.
One variable, one channel. Each variable gets exactly one visual channel (x, y, color, shape, facet). When you start doubling up, the reader stops being able to decode the plot.

Set Up

Load libraries

# Same stack as class, minus haven -- our data is already in R, so nothing to import

library("tidyverse")
library("labelled")
library("forcats")

Load data

We are using airquality, which ships with base R. No file path, no read_sas(): just type the name and it is there.

It is daily air quality measurements from New York City, May–September 1973: ozone, solar radiation, wind, and temperature. It is a good stand-in for NYC HANES because it has the same shape as the data you will meet all term: an environmental exposure measured on a continuous scale, some covariates, a numeric-coded categorical variable, and real missing data.

# note number of variables and observations
glimpse(airquality)

## Rows: 153
## Columns: 6
## $ Ozone   <int> 41, 36, 12, 18, NA, 28, 23, 19, 8, NA, 7, 16, 11, 14, 18, 14, …
## $ Solar.R <int> 190, 118, 149, 313, NA, NA, 299, 99, 19, 194, NA, 256, 290, 27…
## $ Wind    <dbl> 7.4, 8.0, 12.6, 11.5, 14.3, 14.9, 8.6, 13.8, 20.1, 8.6, 6.9, 9…
## $ Temp    <int> 67, 72, 74, 62, 56, 66, 65, 59, 61, 69, 74, 69, 66, 68, 58, 64…
## $ Month   <int> 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,…
## $ Day     <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18,…

# 153 observations (one per day, May 1 - Sept 30)
# 6 variables

Variable definitions, straight from the R help file (?airquality):

Variable	Meaning	Units
`Ozone`	Mean ozone, 1300–1500 hours at Roosevelt Island	ppb
`Solar.R`	Solar radiation, 0800–1200 hours at Central Park	Langleys
`Wind`	Mean wind speed, 0700 and 1000 hours at LaGuardia	mph
`Temp`	Maximum daily temperature at LaGuardia	degrees F
`Month`	Month	5–9
`Day`	Day of month	1–31

Prepare data

Same clean-up sequence as class: rename, then label, then convert to factors, then derive anything new we need.

# 1. rename to something you can type without looking

aq <- airquality %>%
  rename(ozone = Ozone,
         solar = Solar.R,
         wind  = Wind,
         temp  = Temp,
         month = Month,
         day   = Day)

glimpse(aq)

## Rows: 153
## Columns: 6
## $ ozone <int> 41, 36, 12, 18, NA, 28, 23, 19, 8, NA, 7, 16, 11, 14, 18, 14, 34…
## $ solar <int> 190, 118, 149, 313, NA, NA, 299, 99, 19, 194, NA, 256, 290, 274,…
## $ wind  <dbl> 7.4, 8.0, 12.6, 11.5, 14.3, 14.9, 8.6, 13.8, 20.1, 8.6, 6.9, 9.7…
## $ temp  <int> 67, 72, 74, 62, 56, 66, 65, 59, 61, 69, 74, 69, 66, 68, 58, 64, …
## $ month <int> 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5…
## $ day   <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 1…

# 2. build a real date BEFORE month becomes a factor
#    (we need month as a number to do this, so order matters here)

aq <- aq %>%
  mutate(date = as.Date(paste(1973, month, day, sep = "-")))

# 3. labels -- month is coded 5-9, which is exactly the NYC HANES problem

aq <- aq %>%
  set_value_labels(month = c("May" = 5, "June" = 6, "July" = 7,
                             "August" = 8, "September" = 9))

Mutate_if

# change labelled to factors

aq <- aq %>%
  mutate_if(is.labelled, to_factor)

glimpse(aq)

## Rows: 153
## Columns: 7
## $ ozone <int> 41, 36, 12, 18, NA, 28, 23, 19, 8, NA, 7, 16, 11, 14, 18, 14, 34…
## $ solar <int> 190, 118, 149, 313, NA, NA, 299, 99, 19, 194, NA, 256, 290, 274,…
## $ wind  <dbl> 7.4, 8.0, 12.6, 11.5, 14.3, 14.9, 8.6, 13.8, 20.1, 8.6, 6.9, 9.7…
## $ temp  <int> 67, 72, 74, 62, 56, 66, 65, 59, 61, 69, 74, 69, 66, 68, 58, 64, …
## $ month <fct> May, May, May, May, May, May, May, May, May, May, May, May, May,…
## $ day   <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 1…
## $ date  <date> 1973-05-01, 1973-05-02, 1973-05-03, 1973-05-04, 1973-05-05, 197…

# month is now <fct> with levels in calendar order, not <int>
levels(aq$month)

## [1] "May"       "June"      "July"      "August"    "September"

Derive two new variables

We want a categorical version of our exposure (like pb_cat and hg_cat in class) and a simple two-level grouping variable.

aq <- aq %>%
  mutate(
    # EPA AQI breakpoints for ozone, in ppb
    ozone_cat = cut(ozone,
                    breaks = c(-Inf, 54, 70, 85, 105, Inf),
                    labels = c("Good", "Moderate", "Unhealthy (sensitive)",
                               "Unhealthy", "Very unhealthy")),
    hot_day = factor(if_else(temp >= 80, "Hot (80F+)", "Cooler (<80F)"),
                     levels = c("Cooler (<80F)", "Hot (80F+)"))
  )

table(aq$ozone_cat, useNA = "ifany")

## 
##                  Good              Moderate Unhealthy (sensitive) 
##                    83                     8                    13 
##             Unhealthy        Very unhealthy                  <NA> 
##                     5                     7                    37

Say this out loud in section. Those cut points are the EPA’s 8-hour ozone AQI breakpoints, and ozone here is a 2-hour midday average. So this is an illustrative category, not an official AQI. Real analyses live or die on details like this – get in the habit of stating the mismatch instead of quietly hoping nobody asks.

Know your types before you plot

Before you write a single ggplot() call, you should be able to fill in this table for your own data. It tells you which row of the chart chooser you are on.

tibble(
  variable = names(aq),
  class    = map_chr(aq, ~ class(.x)[1])
)

## # A tibble: 9 × 2
##   variable  class  
##   <chr>     <chr>  
## 1 ozone     integer
## 2 solar     integer
## 3 wind      numeric
## 4 temp      integer
## 5 month     factor 
## 6 day       integer
## 7 date      Date   
## 8 ozone_cat factor 
## 9 hot_day   factor

ozone, solar, wind, temp are continuous. month, ozone_cat, hot_day are categorical. date is time. day is a number that is really an index – careful with that one.

The grammar, one layer at a time

Week 3 slides 4–10 showed the layers. Here they are as running code. Same plot, built up one line at a time.

# Layer 1: data only. ggplot knows what to plot with, but not what to draw.
ggplot(data = aq)

# Layer 2: aesthetic mapping. Now it knows the axes -- note the scales appear.
ggplot(data = aq, mapping = aes(x = temp, y = ozone))

# Layer 3: a geom. NOW something is drawn.
ggplot(data = aq, mapping = aes(x = temp, y = ozone)) +
  geom_point()

## Warning: Removed 37 rows containing missing values or values outside the scale range
## (`geom_point()`).

# Layer 4: everything else is decoration and refinement
ggplot(data = aq, mapping = aes(x = temp, y = ozone)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "lm", se = FALSE) +
  ggtitle("Ozone rises with temperature, NYC 1973") +
  xlab("Temperature (degrees F)") +
  ylab("Ozone (ppb)")

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 37 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 37 rows containing missing values or values outside the scale range
## (`geom_point()`).

That warning about removed rows is not a bug. ozone has 37 missing values, and ggplot is telling you it silently dropped those days. Read your warnings. If that number ever surprises you, go back to your data, not your plot.

ONE variable

One continuous variable: histogram

Question 1: What does the ozone distribution look like?

# Default bins -- ggplot picks 30 and warns you to choose deliberately
aq %>%
  ggplot(aes(x = ozone)) +
  geom_histogram()

## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.

## Warning: Removed 37 rows containing non-finite outside the scale range
## (`stat_bin()`).

# Name it if you want to reuse it
hist_ozone <- aq %>%
  ggplot(aes(x = ozone)) +
  geom_histogram(color = "white", fill = "steelblue", bins = 20) +
  ggtitle("Daily ozone, NYC May-Sept 1973") +
  xlab("Ozone (ppb)") +
  ylab("Number of days")
hist_ozone

## Warning: Removed 37 rows containing non-finite outside the scale range
## (`stat_bin()`).

# Too few bins: real structure gets smoothed away
aq %>%
  ggplot(aes(x = ozone)) +
  geom_histogram(color = "white", fill = "steelblue", bins = 5) +
  ggtitle("5 bins -- oversmoothed")

## Warning: Removed 37 rows containing non-finite outside the scale range
## (`stat_bin()`).

# Too many bins: noise looks like structure
aq %>%
  ggplot(aes(x = ozone)) +
  geom_histogram(color = "white", fill = "steelblue", bins = 60) +
  ggtitle("60 bins -- undersmoothed")

## Warning: Removed 37 rows containing non-finite outside the scale range
## (`stat_bin()`).

The point: bin width is an analytic choice, not a default. Change it, look at 2–3 versions, and pick the one that shows the shape without inventing detail. The honest thing to do in a paper is report the bin width you used.

This distribution is strongly right-skewed – most days are low, a few are very high. That is typical of environmental exposure data, and it is the single most common distribution shape you will meet in this field.

Same variable, three different geoms

# Histogram: shows the actual shape, bin by bin
aq %>% ggplot(aes(x = ozone)) +
  geom_histogram(bins = 20, fill = "steelblue", color = "white") +
  ggtitle("Histogram: shows shape")

## Warning: Removed 37 rows containing non-finite outside the scale range
## (`stat_bin()`).

# Density: smooth version, good for comparing groups later
aq %>% ggplot(aes(x = ozone)) +
  geom_density(fill = "steelblue", alpha = 0.4) +
  ggtitle("Density: smoothed shape")

## Warning: Removed 37 rows containing non-finite outside the scale range
## (`stat_density()`).

# Boxplot: shows median, quartiles, and outliers -- but hides shape
aq %>% ggplot(aes(y = ozone)) +
  geom_boxplot() +
  ggtitle("Boxplot: shows summary and outliers")

## Warning: Removed 37 rows containing non-finite outside the scale range
## (`stat_boxplot()`).

Geom	Shows you	Hides from you
Histogram	The actual shape, gaps, multiple peaks	Exact summary statistics
Density	Shape, cleanly comparable across groups	That it is an estimate; invents smoothness
Boxplot	Median, IQR, outliers; compact	Shape entirely – a boxplot cannot show two peaks

A boxplot will happily draw the same box for a bell-shaped distribution and a two-humped one. When you only have room for one plot of one variable, use a histogram. Boxplots earn their keep when you are comparing many groups.

One categorical variable: bar chart

Question 2: How many days fell in each ozone category?

aq %>%
  ggplot(aes(x = ozone_cat)) +
  geom_bar(fill = "darkorange") +
  ggtitle("Days by ozone category") +
  xlab("EPA-style ozone category") +
  ylab("Number of days")

Two things worth knowing about bars:

# geom_bar COUNTS rows for you (stat = "count", the default)
aq %>% ggplot(aes(x = month)) + geom_bar() + ggtitle("geom_bar: counts rows")

# geom_col plots a number you already computed (this is stat = "identity")
aq %>%
  group_by(month) %>%
  summarise(mean_ozone = mean(ozone, na.rm = TRUE)) %>%
  ggplot(aes(x = month, y = mean_ozone)) +
  geom_col(fill = "darkorange") +
  ggtitle("geom_col: plots a value you supply") +
  ylab("Mean ozone (ppb)")

That is the stat = layer from slide 10, in practice. geom_bar() needs only x; geom_col() needs x and y.

Ordering matters

For an ordinal category (like ozone severity), keep the natural order. For a nominal category (no inherent order), sort by frequency – it makes the chart readable at a glance. This is what forcats is for.

aq %>%
  ggplot(aes(x = fct_infreq(ozone_cat))) +
  geom_bar(fill = "darkorange") +
  ggtitle("fct_infreq(): sorted by count") +
  xlab("Ozone category")

TWO variables

Continuous + categorical

Question 3: Is ozone worse in some months than others?

This is the row of the chart chooser with the most options, so let’s see all four.

# Option 1: boxplot -- the workhorse. Compact, comparable, shows outliers.
aq %>%
  ggplot(aes(x = month, y = ozone)) +
  geom_boxplot() +
  ggtitle("Boxplot") +
  ylab("Ozone (ppb)")

## Warning: Removed 37 rows containing non-finite outside the scale range
## (`stat_boxplot()`).

# Option 2: violin -- boxplot's shape-preserving cousin
aq %>%
  ggplot(aes(x = month, y = ozone)) +
  geom_violin(fill = "grey85") +
  geom_boxplot(width = 0.15) +
  ggtitle("Violin + boxplot") +
  ylab("Ozone (ppb)")

## Warning: Removed 37 rows containing non-finite outside the scale range
## (`stat_ydensity()`).
## Removed 37 rows containing non-finite outside the scale range
## (`stat_boxplot()`).

# Option 3: faceted histogram -- best when you care about the actual shapes
aq %>%
  ggplot(aes(x = ozone)) +
  geom_histogram(bins = 15, fill = "steelblue", color = "white") +
  facet_wrap(~ month) +
  ggtitle("Faceted histogram") +
  xlab("Ozone (ppb)")

## Warning: Removed 37 rows containing non-finite outside the scale range
## (`stat_bin()`).

# Option 4: overlaid density -- best for 2-3 groups you want directly compared
aq %>%
  ggplot(aes(x = ozone, color = month, fill = month)) +
  geom_density(alpha = 0.2) +
  ggtitle("Overlaid density") +
  xlab("Ozone (ppb)")

## Warning: Removed 37 rows containing non-finite outside the scale range
## (`stat_density()`).

# Option 5: show the raw data. With n = 153 you can afford to.
aq %>%
  ggplot(aes(x = month, y = ozone)) +
  geom_jitter(width = 0.2, alpha = 0.5) +
  ggtitle("Jittered points: every day is visible") +
  ylab("Ozone (ppb)")

## Warning: Removed 37 rows containing missing values or values outside the scale range
## (`geom_point()`).

How to choose:

2–3 groups, shape matters -> overlaid density
4–8 groups -> boxplot
Many groups, or a report where compactness wins -> boxplot
Shape matters and groups are few -> facets
Small n (say under ~200) -> show the points, alone or on top of a boxplot

Note facet_wrap() gives every panel the same axes, which is what makes the panels comparable. That is a feature, not a limitation.

Continuous + continuous: scatterplot

Question 4: What is the relationship between temperature and ozone?

temp_v_ozone <- aq %>%
  ggplot(aes(x = temp, y = ozone)) +
  geom_point(alpha = 0.5) +
  ggtitle("Temperature vs. ozone") +
  xlab("Temperature (degrees F)") +
  ylab("Ozone (ppb)")
temp_v_ozone

## Warning: Removed 37 rows containing missing values or values outside the scale range
## (`geom_point()`).

# Add a trendline. "lm" = straight line; "loess" = let the data bend it.
temp_v_ozone + geom_smooth(method = "lm", se = FALSE)

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 37 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 37 rows containing missing values or values outside the scale range
## (`geom_point()`).

temp_v_ozone + geom_smooth(method = "loess", se = TRUE)

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 37 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 37 rows containing missing values or values outside the scale range
## (`geom_point()`).

The loess curve bends upward at the high end. The straight line does not. That is a real finding, not a cosmetic difference: ozone accelerates with heat rather than tracking it proportionally. Fit the straight line first, then check whether it deserved to be straight.

Two scatterplot fixes worth memorizing

# Problem: right-skewed y makes the low values pile up at the bottom
# Fix: log scale
temp_v_ozone +
  scale_y_log10() +
  ggtitle("Same data, log-10 y axis")

## Warning: Removed 37 rows containing missing values or values outside the scale range
## (`geom_point()`).

# Problem: outliers stretch the axes so everything else is squashed
# Fix: coord_cartesian() zooms WITHOUT throwing away data
temp_v_ozone +
  coord_cartesian(ylim = c(0, 100)) +
  ggtitle("Zoomed with coord_cartesian()")

## Warning: Removed 37 rows containing missing values or values outside the scale range
## (`geom_point()`).

Use coord_cartesian(), not ylim(). ylim() deletes the points outside the range before any statistics run, so your trendline changes. coord_cartesian() just moves the camera.

Categorical + categorical

Question 5: Does the mix of ozone categories change across the season?

Three positions, three different questions. This trips people up constantly.

# position = "stack" (default): total height = total count
aq %>%
  filter(!is.na(ozone_cat)) %>%
  ggplot(aes(x = month, fill = ozone_cat)) +
  geom_bar(position = "stack") +
  ggtitle('position = "stack": how many days total, and of what kind')

# position = "dodge": side by side, easy to compare one category across months
aq %>%
  filter(!is.na(ozone_cat)) %>%
  ggplot(aes(x = month, fill = ozone_cat)) +
  geom_bar(position = "dodge") +
  ggtitle('position = "dodge": compare a single category across groups')

# position = "fill": every bar is 100%, shows composition only
aq %>%
  filter(!is.na(ozone_cat)) %>%
  ggplot(aes(x = month, fill = ozone_cat)) +
  geom_bar(position = "fill") +
  ggtitle('position = "fill": composition, ignoring group size') +
  ylab("Proportion of days")

Position	Answers	Loses
`"stack"`	How many days total, and what were they?	Hard to compare middle segments
`"dodge"`	Is July worse than June for this category?	Total per month
`"fill"`	What fraction of days were unhealthy?	How many days there were at all

The "fill" version is the honest one for comparing months here, since months have different numbers of non-missing days – but only if you also report the denominators somewhere. A 100% bar built on 8 observations looks identical to one built on 800.

Time + continuous: line plot

This is the nigerm96 example from slides 15–18, on our data.

# Line: the x axis is ordered and continuous, so connecting points is meaningful
aq %>%
  ggplot(aes(x = date, y = ozone)) +
  geom_line(color = "steelblue") +
  ggtitle("Daily ozone across the 1973 season") +
  xlab("Date") +
  ylab("Ozone (ppb)")

# Points + line together: shows where data actually exists
aq %>%
  ggplot(aes(x = date, y = ozone)) +
  geom_line(color = "grey60") +
  geom_point(size = 1) +
  geom_smooth(method = "loess", se = FALSE, color = "darkred") +
  ggtitle("Daily ozone with a seasonal trend") +
  xlab("Date") +
  ylab("Ozone (ppb)")

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 37 rows containing non-finite outside the scale range
## (`stat_smooth()`).

## Warning: Removed 37 rows containing missing values or values outside the scale range
## (`geom_point()`).

The gaps in the line are the missing days. A line plot is the one geom that visibly reveals missingness – the histogram and boxplot just quietly dropped those 37 days.

Only use geom_line() when the x axis has a meaningful order. Connecting categories with a line tells your reader “there is a path between these,” which is false for things like race or borough. If you cannot say what the space between two x values means, do not draw a line through it.

THREE or more variables

Now we are on slide 13’s “3 factors” case. Two ways to add a third variable: a new aesthetic or facets.

# Third variable as COLOR -- continuous variable gives a gradient legend
aq %>%
  ggplot(aes(x = temp, y = ozone, color = wind)) +
  geom_point(size = 2) +
  ggtitle("Third variable as a color gradient") +
  xlab("Temperature (degrees F)") +
  ylab("Ozone (ppb)")

## Warning: Removed 37 rows containing missing values or values outside the scale range
## (`geom_point()`).

# Third variable as COLOR, but categorical -- discrete legend
aq %>%
  ggplot(aes(x = temp, y = ozone, color = month)) +
  geom_point(size = 2) +
  ggtitle("Categorical color: discrete legend") +
  xlab("Temperature (degrees F)") +
  ylab("Ozone (ppb)")

## Warning: Removed 37 rows containing missing values or values outside the scale range
## (`geom_point()`).

# Third variable as SHAPE (works for categorical only, and only up to ~5 levels)
aq %>%
  ggplot(aes(x = temp, y = ozone, shape = hot_day)) +
  geom_point(size = 2, alpha = 0.6) +
  ggtitle("Categorical shape")

## Warning: Removed 37 rows containing missing values or values outside the scale range
## (`geom_point()`).

# Third variable as FACETS
aq %>%
  ggplot(aes(x = temp, y = ozone)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "lm", se = FALSE) +
  facet_wrap(~ month) +
  ggtitle("Faceted by month, with a trendline in each panel")

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 37 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 37 rows containing missing values or values outside the scale range
## (`geom_point()`).

Same color = argument, two completely different legends. ggplot read the class of the variable and picked a continuous or discrete scale accordingly. This is the single best argument for doing your factor conversion properly before you plot.

Four variables, and knowing when to stop

aq %>%
  ggplot(aes(x = temp, y = ozone, color = wind)) +
  geom_point(size = 2) +
  facet_wrap(~ month) +
  ggtitle("Temp vs. ozone, colored by wind, faceted by month") +
  xlab("Temperature (degrees F)") +
  ylab("Ozone (ppb)")

## Warning: Removed 37 rows containing missing values or values outside the scale range
## (`geom_point()`).

Four variables is roughly the ceiling for a plot someone can read in a talk. The science here is genuinely visible: the highest ozone days are hot and calm (dark points, top right), and they cluster in July and August. That is the photochemistry – ozone forms when sunlight cooks precursor pollutants on hot, stagnant days.

If you find yourself wanting a fifth variable, that is your data telling you it wants a model, not another aesthetic.

Same data, wrong chart

Three failure modes, each of which you will see in a published paper this year.

1. A category R thinks is a number

This is the class workflow’s whole reason for existing.

# WRONG: month left as an integer
airquality %>%
  ggplot(aes(x = Temp, y = Ozone, color = Month)) +
  geom_point(size = 2) +
  ggtitle("WRONG: month as a number -- a meaningless blue gradient")

## Warning: Removed 37 rows containing missing values or values outside the scale range
## (`geom_point()`).

# RIGHT: month converted to a factor first
aq %>%
  ggplot(aes(x = temp, y = ozone, color = month)) +
  geom_point(size = 2) +
  ggtitle("RIGHT: month as a factor -- five readable groups")

## Warning: Removed 37 rows containing missing values or values outside the scale range
## (`geom_point()`).

The first plot implies that “6.5” is a month. Nobody meant to say that, but the plot says it.

2. Bars where a distribution belongs

# WRONG: bar of means. Five numbers standing in for 116 observations.
aq %>%
  group_by(month) %>%
  summarise(mean_ozone = mean(ozone, na.rm = TRUE)) %>%
  ggplot(aes(x = month, y = mean_ozone)) +
  geom_col(fill = "grey60") +
  ggtitle("WRONG: bar of means hides everything") +
  ylab("Mean ozone (ppb)")

# RIGHT: same comparison, distribution intact
aq %>%
  ggplot(aes(x = month, y = ozone)) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(width = 0.2, alpha = 0.4) +
  ggtitle("RIGHT: boxplot + points shows spread, skew, and n") +
  ylab("Ozone (ppb)")

## Warning: Removed 37 rows containing non-finite outside the scale range
## (`stat_boxplot()`).

## Warning: Removed 37 rows containing missing values or values outside the scale range
## (`geom_point()`).

The bar chart cannot tell you that June has almost no high-ozone days while July has several extreme ones, or that the July mean is being dragged by a handful of observations. Bars are for counts. Distributions get boxplots.

3. Angles instead of lengths

# WRONG: a pie chart (a stacked bar in polar coordinates)
aq %>%
  filter(!is.na(ozone_cat)) %>%
  ggplot(aes(x = "", fill = ozone_cat)) +
  geom_bar(width = 1) +
  coord_polar(theta = "y") +
  ggtitle("WRONG: which is bigger, Unhealthy or Very unhealthy?") +
  xlab("") + ylab("")

# RIGHT: a bar chart
aq %>%
  filter(!is.na(ozone_cat)) %>%
  ggplot(aes(x = ozone_cat)) +
  geom_bar(fill = "darkorange") +
  ggtitle("RIGHT: lengths on a common baseline are easy to compare") +
  xlab("Ozone category") + ylab("Number of days")

People compare lengths against a shared baseline accurately, and angles poorly. That is the entire case against pie charts, and it is why ggplot makes you work to build one.

Making it presentation-ready

Everything above used ggtitle()/xlab()/ylab() to match class. labs() does all of it in one call, and adds a few slots you cannot get otherwise.

final_plot <- aq %>%
  ggplot(aes(x = temp, y = ozone, color = month)) +
  geom_point(size = 2, alpha = 0.8) +
  geom_smooth(method = "lm", se = FALSE, color = "black", linewidth = 0.5) +
  labs(
    title    = "Ozone rises steeply with temperature in NYC",
    subtitle = "Daily measurements, May-September 1973",
    x        = "Maximum daily temperature (degrees F)",
    y        = "Mean midday ozone (ppb)",
    color    = "Month",
    caption  = "Source: New York State Dept. of Conservation and NOAA, via R datasets"
  ) +
  theme_minimal() +
  theme(legend.position = "bottom")

final_plot

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 37 rows containing non-finite outside the scale range
## (`stat_smooth()`).

## Warning: Removed 37 rows containing missing values or values outside the scale range
## (`geom_point()`).

labs(color = "Month") renames the legend – that one is worth memorizing, because a legend titled ozone_cat is the fastest way to make a figure look unfinished.

To save a plot for a homework write-up:

ggsave("ozone_by_temp.png", plot = final_plot,
       width = 7, height = 4.5, dpi = 300)

Cheat sheet

Pick your row:

Your variables	Plot
1 continuous	histogram (density, boxplot)
1 categorical	bar chart
continuous by categorical	boxplot (violin, faceted histogram, overlaid density)
continuous by continuous	scatterplot + smooth
categorical by categorical	bar with `position = "dodge"` / `"fill"`
continuous over time	line
add one more	`color =`, `shape =`, or `facet_wrap()`

Then check yourself:

Are my categorical variables actually factors?
Did I read the “removed N rows” warning?
Is my bin width or smoother a choice I can defend?
Does every axis have units in its label?
Can a reader who has never seen this data read the title and know what they are looking at?

Practice

Try these before looking at the solutions. Use aq.

What is the distribution of wind speed? Is it skewed like ozone?
Do hot days and cooler days differ in their solar radiation? Pick a geom and justify it.
Is the wind–ozone relationship different in July than in May?
One plot, three variables: does the temperature–ozone relationship depend on whether the day was windy?

Solutions

# 1. One continuous variable -> histogram. Wind is roughly symmetric, unlike ozone.
aq %>%
  ggplot(aes(x = wind)) +
  geom_histogram(bins = 20, fill = "steelblue", color = "white") +
  labs(title = "Wind speed is close to symmetric", x = "Wind (mph)", y = "Days")

# 2. Continuous by categorical, 2 groups -> boxplot or overlaid density.
#    Two groups and shape matters, so density is defensible; boxplot is safer.
aq %>%
  ggplot(aes(x = hot_day, y = solar)) +
  geom_boxplot() +
  geom_jitter(width = 0.15, alpha = 0.4) +
  labs(title = "Solar radiation by temperature group",
       x = NULL, y = "Solar radiation (Langleys)")

## Warning: Removed 7 rows containing non-finite outside the scale range
## (`stat_boxplot()`).

## Warning: Removed 7 rows containing missing values or values outside the scale range
## (`geom_point()`).

# 3. Continuous by continuous, split by a category -> scatterplot + facets
aq %>%
  filter(month %in% c("May", "July")) %>%
  ggplot(aes(x = wind, y = ozone)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "lm", se = FALSE) +
  facet_wrap(~ month) +
  labs(title = "Wind suppresses ozone more steeply in July",
       x = "Wind (mph)", y = "Ozone (ppb)")

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 10 rows containing non-finite outside the scale range
## (`stat_smooth()`).

## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).

# 4. Three variables -> two on the axes, one on color.
#    Bin wind first so the color legend is readable.
aq %>%
  mutate(windy = if_else(wind >= median(wind), "Windier", "Calmer")) %>%
  ggplot(aes(x = temp, y = ozone, color = windy)) +
  geom_point(alpha = 0.7, size = 2) +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title = "Heat drives ozone, but only when the air is still",
       x = "Temperature (degrees F)", y = "Ozone (ppb)", color = NULL)

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 37 rows containing non-finite outside the scale range
## (`stat_smooth()`).

## Warning: Removed 37 rows containing missing values or values outside the scale range
## (`geom_point()`).

Question 4 is the payoff. The two lines are not parallel: on calm days ozone climbs sharply with temperature, and on windy days the same heat produces much less ozone. That is one plot carrying an interaction – and it is the visual version of the modeling you will do later in the term.

TA Session 3: Choosing the Right Plot in ggplot2

Companion to Week 3 – Visualization using ggplot

TA Session

07/25/2026