Chapter 2
Describing data

Lectures 3–4 · centre, spread, shape, and pictures that tell the truth

Before anyone can predict from data or infer anything about the mechanism behind them, someone has to look at the data. That sounds trivial. It isn't. A table of 300 numbers tells the human eye almost nothing, and a single summary number can hide exactly the thing that matters. An average depth of three feet doesn't mean you can wade across the river.

This chapter gives you the tools for looking: numbers that summarize where data sit (the centre), how much they vary (the spread), and what pattern they make (the shape), and the pictures that show all three at once. Along the way we'll meet a geyser that keeps two timetables, a physicist who measured the speed of light slightly wrong, and a university that seemed to discriminate against women until someone looked more closely. Every named dataset is real and built into R, so you can explore them yourself.

Every summary in this chapter comes with three things:

1. a way to compute it by hand, step by step, so that you know exactly what it is (and can do it in an exam);
2. the R code that computes it, so that you can use it on real datasets with hundreds of values;
3. its strengths and weaknesses, because no single summary is best for every job.

The by-hand answers and the R answers always agree. If yours don't, that's a sign to look again at one of the steps.

Describing data is not an end in itself: it is the first step of both of the course's big goals. For prediction, a summary of past data is the raw material of every forecast: the centre of a distribution is our best single guess for the next value (Section 2.2 makes this precise), and the spread tells us how far off that guess is likely to be. For inference, the pictures and summaries of a sample are our first clues about the mechanism that produced it: two peaks hint at two hidden groups, a difference between boxplots hints at a real effect, and a pattern that reverses when the data are split warns us that something else is at work. Keep both questions in mind as you read: what would I predict for the next case? and what does this suggest about the process behind the data?

Lecture 3 — Shape, centre and spread

2.1 Distributions: the shape of data

A geyser with two timetables

Old Faithful, in Yellowstone National Park, is the most famous geyser in the world. Every hour or two it shoots water dozens of metres into the air, and crowds gather to watch. Park rangers post a prediction of the next eruption so visitors know when to arrive. How do they predict it?

R's faithful dataset records 272 eruptions: how long each lasted, and how long visitors then waited for the next one. Here are the eruption lengths.

summary(faithful$eruptions)       # minutes

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.600   2.163   4.000   3.488   4.454   5.100 

ggplot(faithful, aes(x = eruptions)) +
  geom_histogram(binwidth = 0.2, boundary = 1.6) +
  labs(x = "eruption length (minutes)", y = "number of eruptions")

The pattern of values in a dataset (which values occur, and how often) is its distribution (Section 1.4). A histogram is the standard picture of the distribution of a numerical variable.

Drawing a histogram

A histogram is built from a frequency table: a list of intervals and the number of observations in each.

The bars of a histogram touch, unlike the bars of a bar chart for categories (Section 2.5), because a numerical variable runs continuously along the number line. A gap in a histogram means "no observations here", which is information.

expt1 <- morley$Speed[morley$Expt == 1]
expt1
h <- hist(expt1, breaks = seq(600, 1100, by = 100), right = FALSE,
          main = NULL, xlab = "speed (km/s minus 299,000)", ylab = "frequency")
h$counts          # the frequencies, bin by bin

 [1]  850  740  900 1070  930  850  950  980  980  880 1000  980  930
[14]  650  760  810 1000 1000  960  960
[1] 1 2 4 9 4

A histogram depends on two choices you make: the bin width and the starting point. Too wide and every detail is smoothed away; too narrow and the picture is all noise. Try it on Old Faithful and three other datasets from this chapter.

There is no single correct bin width. A good habit is to try several, and to trust only the features (peaks, gaps, tails) that survive every reasonable choice.

Describing a distribution

When you describe a distribution, describe three things:

• shape: how many peaks, whether it is symmetric or skewed, and whether there are gaps or unusual values;
• centre: where a typical value lies;
• spread: how much the values vary.

Statisticians use a small vocabulary for shape. It's worth learning exactly, because "skewed to the right" is easy to get backwards.

Explore the shapes below. Pay attention to where the mean and median sit for each one: the relationship between them is the quickest numerical clue to the shape.

Most quantities that can't go below zero but occasionally get very large (incomes, waiting times, pollution levels, drug concentrations) are right-skewed. Left skew is rarer, and usually appears when there is a natural upper limit that most values crowd against (scores on an easy test), or when a few values come out unusually low (Michelson's measurements).

Histograms: strengths and weaknesses

2.2 Centre: what is a typical value?

Mean and median

Suppose we have observations $x_1,x_2,\dots ,x_n$. The subscript $i$ in $x_i$ just labels the observations in the order they were recorded. When we sort them from smallest to largest, we write $x_{(1)}$$\le x_{(2)}$$\le \cdots$$\le x_{(n)}$: the brackets around the subscript mean "the $i$th smallest". There are two classic ways to describe the centre.

The two measure "typical" in different ways. The mean is the balance point of the data: put an equal weight at each observation on a see-saw, and the mean is where it balances. A value far from the others sits at the end of a long lever, so it has a lot of influence. The median is the halfway point: half the observations are at or below it and half at or above. It only cares about the order of the values, not how far out the extreme ones are.

conc  <- c(12.1, 9.8, 15.3, 11.0, 10.4)
mean(conc); median(conc)
sprayA <- InsectSprays$count[InsectSprays$spray == "A"]
sprayA
sort(sprayA)
mean(sprayA); median(sprayA)

[1] 11.72
[1] 11
 [1] 10  7 20 14 14 12 10 23 17 20 14 13
 [1]  7 10 10 12 13 14 14 14 17 20 20 23
[1] 14.5
[1] 14

two <- rbind(data.frame(data = "toxin (ng/mL)", x = conc),
             data.frame(data = "insects, spray A", x = sprayA))
lines <- aggregate(x ~ data, data = two, FUN = function(v) c(mean = mean(v), median = median(v)))
lines <- data.frame(data = lines$data, mean = lines$x[, "mean"], median = lines$x[, "median"])
two$stack <- ave(two$x, two$data, two$x, FUN = seq_along)   # stack repeated values
ggplot(two, aes(x = x, y = stack)) +
  geom_point(size = 3) +
  geom_vline(data = lines, aes(xintercept = mean)) +
  geom_vline(data = lines, aes(xintercept = median), linetype = "dashed") +
  geom_text(data = lines, aes(x = mean, y = 3.4, label = "mean"), hjust = -0.15, size = 3.8) +
  geom_text(data = lines, aes(x = median, y = 3.4, label = "median"), hjust = 1.15, size = 3.8) +
  facet_wrap(~ data, scales = "free_x", ncol = 1) +
  scale_y_continuous(NULL, breaks = NULL, limits = c(0.5, 3.8)) + labs(x = NULL)

For symmetric data the mean and median agree; for skewed data they don't. Which one should you report? That depends on the question you are answering.

oz <- airquality$Ozone
sum(!is.na(oz))                # number of days with a reading
mean(oz, na.rm = TRUE)
median(oz, na.rm = TRUE)

[1] 116
[1] 42.12931
[1] 31.5

ggplot(airquality, aes(x = Ozone)) +
  geom_histogram(binwidth = 10, boundary = 0, colour = "white", na.rm = TRUE) +
  geom_vline(xintercept = mean(oz, na.rm = TRUE)) +
  geom_vline(xintercept = median(oz, na.rm = TRUE), linetype = "dashed") +
  labs(x = "ozone (ppb)", y = "number of days")

The mode

For the insect counts in Example 2.3(b), the value 14 occurs three times and nothing else occurs more than twice, so the mode is 14. For glass type in fgl, the counts are WinF 70, WinNF 76, Veh 17, Con 13, Tabl 9 and Head 29, so the modal type is WinNF (building window, non-float). The mode is the only measure of centre that makes sense for a nominal variable: there is no "average glass type" and no "middle glass type". Base R has no function for the mode, but table() does the work:

table(sprayA)                       # how often each count occurs
names(which.max(table(sprayA)))     # the most frequent value
table(MASS::fgl$type)

sprayA
 7 10 12 13 14 17 20 23 
 1  2  1  1  3  1  2  1 
[1] "14"

 WinF WinNF   Veh   Con  Tabl  Head 
   70    76    17    13     9    29 

Resistance and the trimmed mean

To see the difference, imagine changing the largest value in a dataset to a billion. The mean shoots up towards a billion. The median doesn't move at all, because the largest value is still the largest, and the middle of the sorted list is unchanged.

Is there a compromise between the two? Yes: throw away the most extreme values at both ends, and average what is left.

Trimmed means are used where outliers are a real worry and the mean is still wanted: Olympic diving and gymnastics drop the highest and lowest judges' scores before averaging, which is a trimmed mean that stops one biased judge deciding a medal.

bac <- c(0.081, 0.079, 0.083, 0.080, 0.82)
mean(bac)
median(bac)
mean(bac, trim = 0.2)    # drops the lowest and highest 20%: one value each

[1] 0.2286
[1] 0.081
[1] 0.08133333

ggplot(data.frame(bac, y = c(-0.3, 0.3, 0, -0.15, 0.15)), aes(x = bac, y = y)) +
  geom_point(size = 3, alpha = 0.8) +
  geom_vline(xintercept = mean(bac)) +
  geom_vline(xintercept = c(median(bac), mean(bac, trim = 0.2)), linetype = "dashed") +
  annotate("text", x = mean(bac) + 0.01, y = -0.8, label = "mean", hjust = 0, size = 3.8) +
  annotate("text", x = median(bac) + 0.01, y = 0.8, label = "median and\ntrimmed mean", hjust = 0, size = 3.8) +
  scale_y_continuous(NULL, breaks = NULL, limits = c(-1, 1)) +
  labs(x = "blood alcohol reading (g/dL)")

The trimmed mean is not one summary but a whole family, with the mean and the median at its two ends.

• Trim 0%: nothing is removed, and the trimmed mean is the ordinary mean.
• Trim a little: the most extreme values go, and the result moves away from the mean.
• Trim almost 50% from each end: only the middle one or two values remain, and their average is the median. (R treats trim = 0.5 as the median.)

So choosing a trimming fraction is choosing how much to trust the extremes. The widget below makes the journey visible on the ozone data.

Choosing a measure of centre

Why the mean and median are "best" guesses

There is a deeper way to think about centre, and it connects directly to Chapter 1's airline. Suppose you must summarize the data with a single number $c$, and every data value charges you a penalty depending on how far it is from $c$. Which $c$ makes the total penalty smallest? The answer depends on the penalty rule. Try it before reading the theorem.

So the mean is the best single guess if big errors are especially costly (squared penalty: an error of 10 costs 100 times as much as an error of 1), the median is the best guess if every unit of error costs the same (absolute penalty), and the mode is the best guess if only an exact hit counts. This is the same idea as Chapter 1's airline: the best prediction depends on the cost of being wrong.

2.3 Spread: how much do the data vary?

Two drug-testing instruments both give the right answer on average. One gives readings within 1% of the truth; the other is sometimes 20% off. Their means are identical; only their spread tells them apart. In science, medicine and forensics, spread is often the more important number.

The range

For the toxin concentrations of Example 2.3, the range is $15.3$$-9.8$$=5.5$ ng/mL. The range is easy to compute and to explain, but it depends on only two observations, the most extreme ones, so a single outlier changes it completely, and it tends to grow as you collect more data. (R's range(x) returns the minimum and maximum; diff(range(x)) gives their difference.)

Variance and standard deviation

A better measure of spread uses every observation: how far, typically, is an observation from the mean? The distances from the mean, $x_i-\bar{x}$, are called deviations. We can't just average them, because they always add to zero (part 4 of the theorem in Section 2.2): the positive and negative ones cancel. So we square them first, making them all positive, average the squares, and then take a square root at the end to get back to the original units.

conc <- c(12.1, 9.8, 15.3, 11.0, 10.4)
conc - mean(conc)                   # the deviations
sum((conc - mean(conc))^2)          # sum of squared deviations
var(conc)
sd(conc)

[1]  0.38 -1.92  3.58 -0.72 -1.32
[1] 18.908
[1] 4.727
[1] 2.174167

d <- data.frame(sample = 1:5, conc = conc)
ggplot(d, aes(x = sample, y = conc)) +
  geom_hline(yintercept = mean(conc), linetype = "dashed") +
  geom_segment(aes(xend = sample, yend = mean(conc)), colour = "orange") +
  geom_point(size = 3) +
  labs(x = "sample", y = "toxin (ng/mL)")

Why divide by $n-1$ rather than $n$? The deviations are measured from $\bar{x}$, which was computed from the same data and therefore sits right in the middle of them. Theorem part 1 says $\bar{x}$ is the point that makes the sum of squared deviations as small as possible, so measuring from $\bar{x}$ gives a slightly smaller sum than measuring from the true population mean would. Dividing by $n-1$ instead of $n$ inflates the answer just enough to correct for this on average. Chapter 8 makes this precise: $s^2$ is an unbiased estimator of the population variance. For large $n$ the difference hardly matters.

Reading a standard deviation

How should you read a standard deviation? For many datasets, especially roughly symmetric, single-peaked ones, the following rule of thumb holds.

Here is the rule at work on the birth weights of Chapter 1. The darker band is one standard deviation either side of the mean; the lighter band, two.

bw <- MASS::birthwt$bwt; m <- mean(bw); s <- sd(bw)
c(within_1_sd = mean(abs(bw - m) <= s), within_2_sd = mean(abs(bw - m) <= 2 * s))
ggplot(data.frame(bw), aes(bw)) +
  annotate("rect", xmin = m - 2 * s, xmax = m + 2 * s, ymin = -Inf, ymax = Inf, fill = "steelblue", alpha = 0.12) +
  annotate("rect", xmin = m - s, xmax = m + s, ymin = -Inf, ymax = Inf, fill = "steelblue", alpha = 0.22) +
  geom_histogram(binwidth = 250, boundary = 0, colour = "white") +
  geom_vline(xintercept = m) +
  labs(x = "birth weight (g)", y = "number of babies")

within_1_sd within_2_sd 
  0.6772487   0.9629630 

For the 189 birth weights in birthwt (mean 2,945 g, standard deviation 729 g), one standard deviation either side of the mean is 2,215 g to 3,674 g, and 67.7% of babies fall in it; two standard deviations either side is 1,486 g to 4,403 g, and 96.3% fall in it. Very close to the rule. But for skewed or two-peaked data the rule can fail badly. The next result holds for every dataset, whatever its shape.

Chebyshev's guarantee is weak (75% instead of the rule of thumb's 95%), but it never fails. That's the trade-off between assuming something about the shape and assuming nothing.

Quantiles and quartiles

The standard deviation squares deviations, so it is even less resistant than the mean: one extreme value can inflate it enormously. Resistant measures of spread are built from quantiles instead.

The rule spreads the $n$ sorted values evenly from "0% of the way along" (the minimum, position 1) to "100% of the way along" (the maximum, position $n$). The 25% point is then a quarter of the way along, at position $1+(n-1)/4$, and if that falls between two observations we interpolate.

R's plain mad(x) multiplies the MAD by 1.4826. That constant makes the MAD comparable to the standard deviation for bell-shaped data, so that the two can be compared directly. Always say which version you are reporting.

quantile(sprayA)               # min, Q1, median, Q3, max
IQR(sprayA)
mad(sprayA, constant = 1)      # the plain MAD
mad(sprayA)                    # R's default: 1.4826 x MAD
sd(sprayA)                     # for comparison

   0%   25%   50%   75%  100% 
 7.00 11.50 14.00 17.75 23.00 
[1] 6.25
[1] 3.5
[1] 5.1891
[1] 4.719399

speed <- morley$Speed + 299000        # undo the coding: km/s
mean(speed); median(speed)
sd(speed); IQR(speed)
mean(speed) - 299734.5                # how far off is the average (in air)?

[1] 299852.4
[1] 299850
[1] 79.01055
[1] 85
[1] 117.9

ggplot(morley, aes(x = factor(Expt), y = Speed + 299000)) +
  geom_boxplot() +
  geom_hline(yintercept = 299734.5, linetype = "dashed") +
  labs(x = "experiment", y = "speed of light (km/s)")

Choosing a measure of spread

A good rule of thumb: report the mean with the standard deviation for roughly symmetric data, and the median with the IQR for skewed data or data with outliers.

Lecture 4 — Boxplots, honest pictures and Simpson's paradox

2.4 Position, outliers and boxplots

Is a birth weight of 2,500 g unusual? Is a glass fragment's refractive index consistent with window glass? Questions like these ask about one observation's position within a distribution. Quantiles (Section 2.3) answer one version: a baby "on the 10th percentile" for weight is lighter than 90% of babies of the same age. z-scores answer another.

z-scores

z-scores let us compare observations measured on different scales: is a student who scored 85 on an exam with mean 70 and SD 10 ($z=1.5$) further above average than one who scored 620 on an exam with mean 500 and SD 100 ($z=1.2$)? Yes. And z-scores don't depend on the units: converting grams to pounds changes $x$, $\bar{x}$ and $s$ by the same factor, which cancels. Here is the general fact.

For example, converting Celsius to Fahrenheit ($F=32+1.8C$) multiplies the mean by 1.8 and then adds 32, but only multiplies the standard deviation by 1.8: shifting every temperature up by 32 degrees doesn't make them any more spread out. The theorem also shows that z-scores always have mean 0 and standard deviation 1: the z-score is $y=a+bx$ with $b=1/s$ and $a=-\bar{x}/s$, so its mean is $-\bar{x}/s$$+\bar{x}/s$$=0$ and its SD is $s/s=1$.

library(MASS)
winf <- fgl$RI[fgl$type == "WinF"]
c(mean = mean(winf), sd = sd(winf), n = length(winf))
(6.0 - mean(winf)) / sd(winf)       # z-score of the jacket fragment
mean(winf >= 6.0)                    # fraction of window fragments at least this high

      mean         sd          n 
 0.7182857  2.2680969 70.0000000 
[1] 2.328699
[1] 0.01428571

ggplot(data.frame(RI = winf), aes(RI)) +
  geom_histogram(binwidth = 1, boundary = 0, colour = "white") +
  geom_vline(xintercept = 6.0, linetype = "dashed") +
  labs(x = "refractive index (coded), float window glass", y = "number of fragments")

Outliers

Why 1.5? It's a convention, chosen by John Tukey because it works well in practice: for bell-shaped data, it flags fewer than 1 observation in 100, so a flagged point is genuinely unusual, while still catching real outliers.

An outlier isn't automatically a mistake. It may be a data-entry error (Example 2.5), a measurement under unusual conditions, or the most interesting observation in the dataset. The rule only says: look at this one.

Boxplots

A boxplot (or box-and-whisker plot) draws the five-number summary and the outliers in one compact picture.

Work through the algorithm in the widget, one step at a time.

sprayC <- InsectSprays$count[InsectSprays$spray == "C"]
quantile(sprayC)
boxplot.stats(sprayC)
ggplot(data.frame(count = sprayC), aes(x = count, y = "")) +
  geom_boxplot(width = 0.4) +
  labs(x = "insects remaining (spray C)", y = NULL)

  0%  25%  50%  75% 100% 
 0.0  1.0  1.5  3.0  7.0 
$stats
[1] 0.0 1.0 1.5 3.0 4.0

$n
[1] 12

$conf
[1] 0.5877866 2.4122134

$out
[1] 7

What a boxplot can hide

A boxplot compresses a whole dataset into five numbers plus outliers. That makes it compact, and it also means very different datasets can have identical boxplots.

fit5 <- function(x) {      # stretch x so its five-number summary is 0, 25, 50, 75, 100
  approx(quantile(x), c(0, 25, 50, 75, 100), xout = x)$y
}
n <- 101
three <- data.frame(
  set = rep(c("flat", "two peaks", "three clusters"), each = n),
  value = c(fit5(ppoints(n)),
            fit5(c(qnorm(ppoints(50), -2, 0.5), qnorm(ppoints(51), 2, 0.5))),
            fit5(c(qnorm(ppoints(34), 0, 0.3), qnorm(ppoints(33), 5, 0.3),
                   qnorm(ppoints(34), 10, 0.3)))))
aggregate(value ~ set, data = three, FUN = fivenum)

             set value.1 value.2 value.3 value.4 value.5
1           flat       0      25      50      75     100
2 three clusters       0      25      50      75     100
3      two peaks       0      25      50      75     100

ggplot(three, aes(x = value, y = set)) +
  geom_boxplot(width = 0.5) +
  labs(x = "value", y = NULL)

ggplot(three, aes(x = value)) +
  geom_histogram(binwidth = 5, boundary = 0, colour = "white") +
  facet_wrap(~ set) + scale_x_continuous(breaks = c(0, 50, 100)) +
  labs(x = "value", y = "count")

Side-by-side boxplots

Where boxplots really shine is in comparing a numerical variable across groups. Placed side by side on a common axis, they let you compare centres (the median lines), spreads (the box lengths), skewness (whether the median sits in the middle of its box, and whether one whisker is much longer) and outliers, for many groups at once.

ggplot(InsectSprays, aes(x = spray, y = count)) +
  geom_boxplot(outlier.shape = 1, outlier.size = 2.5) +
  geom_jitter(width = 0.12, height = 0, alpha = 0.3) +
  labs(x = "insecticide", y = "insects remaining")

tapply(InsectSprays$count, InsectSprays$spray, median)
tapply(InsectSprays$count, InsectSprays$spray, IQR)

   A    B    C    D    E    F 
14.0 16.5  1.5  5.0  3.0 15.0 
    A     B     C     D     E     F 
 6.25  5.00  2.00  1.25  2.25 10.00 

Boxplots: strengths and weaknesses

A good habit is to add the raw data as dots on top of a boxplot, as in Example 2.12, whenever there are few enough points to see.

2.5 Pictures that tell the truth

The right picture for the variable

The type of variable (Section 1.4) decides the picture. Every picture is good at some things and bad at others; this table collects them. (The ggplot2 functions are geom_histogram(), geom_boxplot(), geom_jitter() or geom_dotplot(), geom_bar(), geom_point(), geom_line(), and geom_qq() with geom_qq_line(), in the order of the table.)

Bar charts

A bar chart shows a categorical variable: one bar per category, with height equal to the count (or percentage). Unlike a histogram, the bars are separated by gaps, because the categories are separate, and for a nominal variable their order is up to you: sorting the bars by height usually makes the chart easier to read.

table(fgl$type)
ggplot(fgl, aes(x = reorder(type, type, function(t) -length(t)))) +
  geom_bar(width = 0.7) +
  labs(x = "glass type", y = "number of fragments")

 WinF WinNF   Veh   Con  Tabl  Head 
   70    76    17    13     9    29 

Bar charts: easy to read and compare; but bars must start at zero (Section "Honest and dishonest graphs" below), and pie charts showing the same information are much harder to judge by eye.

Scatterplots

A scatterplot shows the relationship between two numerical variables measured on the same individuals (units): one point per unit, with one variable on each axis.

ggplot(airquality, aes(x = Temp, y = Ozone)) +
  geom_point(alpha = 0.7, na.rm = TRUE) +
  labs(x = "maximum temperature (°F)", y = "ozone (ppb)")

sapply(anscombe, mean)
sapply(anscombe, sd)
sapply(1:4, function(i) cor(anscombe[[i]], anscombe[[i + 4]]))

      x1       x2       x3       x4       y1       y2       y3 
9.000000 9.000000 9.000000 9.000000 7.500909 7.500909 7.500000 
      y4 
7.500909 
      x1       x2       x3       x4       y1       y2       y3 
3.316625 3.316625 3.316625 3.316625 2.031568 2.031657 2.030424 
      y4 
2.030579 
[1] 0.8164205 0.8162365 0.8162867 0.8165214

long <- data.frame(set = rep(paste("dataset", 1:4), each = 11),
                   x = unlist(anscombe[1:4]), y = unlist(anscombe[5:8]))
ggplot(long, aes(x, y)) +
  geom_point(size = 2.5) +
  facet_wrap(~ set)

Is it bell-shaped? The Gaussian Q-Q plot

Many methods later in the course assume data come from a bell-shaped (Gaussian) distribution. A histogram is a poor judge of this, especially for small samples, because its appearance depends so much on the bins. A better tool is the Gaussian quantile–quantile (Q-Q) plot. (Most books and software call the bell curve the normal distribution, and this plot a normal Q-Q plot; R's function is qqnorm(). These notes use the name Gaussian, for reasons explained in Section 5.3.)

The idea is simple. Sort the data. Then ask: if the data really came from a bell-shaped distribution, where would the smallest value be expected to sit? The second smallest? The median? A Q-Q plot puts these expected positions (from the Gaussian distribution, Chapter 5) on the horizontal axis and the actual sorted data on the vertical axis. If the data are bell-shaped, actual and expected agree up to a change of units, and the points lie close to a straight line. Curves and bends show how they depart:

• points curving upwards at the right end (and often at the left): the largest values are too large, a sign of right skew;
• points curving downwards at the left end: the smallest values are too small, left skew;
• an S shape, with both ends pulling away from the line in opposite directions: heavy tails (more extreme values at both ends than a bell would produce).

both <- rbind(data.frame(var = "birth weight (g)", value = birthwt$bwt),
              data.frame(var = "ozone (ppb)", value = na.omit(airquality$Ozone)))
ggplot(both, aes(sample = value)) +
  geom_qq() + geom_qq_line() +
  facet_wrap(~ var, scales = "free_y") +
  labs(x = "Gaussian quantiles", y = "data")

Honest and dishonest graphs

A graph is an argument, and it can mislead as easily as words. The most common trick is with the axis.

2.6 Two categorical variables and Simpson's paradox

Contingency tables

When both variables are categorical, we summarize them with a contingency table (a two-way table of counts): one row for each category of one variable, one column for each category of the other, and in each cell the number of units with that combination. Chapter 3 uses one for the Titanic. To compare groups, convert the counts to row or column percentages, so that groups of different sizes can be compared fairly.

overall <- margin.table(UCBAdmissions, c(1, 2))       # pool the departments
overall
round(prop.table(overall, 2), 3)                      # admission rate by sex
by_dept <- prop.table(UCBAdmissions, c(2, 3))["Admitted", , ]   # rate by sex and department
round(by_dept, 3)
margin.table(UCBAdmissions, c(2, 3))                  # number of applicants

          Gender
Admit      Male Female
  Admitted 1198    557
  Rejected 1493   1278
          Gender
Admit       Male Female
  Admitted 0.445  0.304
  Rejected 0.555  0.696
        Dept
Gender       A     B     C     D     E     F
  Male   0.621 0.630 0.369 0.331 0.277 0.059
  Female 0.824 0.680 0.341 0.349 0.239 0.070
        Dept
Gender     A   B   C   D   E   F
  Male   825 560 325 417 191 373
  Female 108  25 593 375 393 341

rates <- as.data.frame(as.table(by_dept))
names(rates) <- c("sex", "dept", "rate")
ggplot(rates, aes(x = dept, y = rate, shape = sex)) +
  geom_point(size = 3.5) +
  geom_hline(yintercept = c(0.445, 0.304), linetype = "dashed") +
  labs(x = "department", y = "admission rate", shape = NULL)