1.1 Where Week 3 sits

Week 1 gave you chart types and categorical variables. Week 2 gave you populations, samples and the descriptive statistics that summarise discrete data. This week we handle continuous numerical variables, where every observation is effectively unique and counting distinct values no longer works.

1.2 What you will be able to do

By the end of this lecture and the accompanying notebook you should be able to:

1. Produce a stem-and-leaf plot and say what it preserves that a histogram discards.
2. Compute and interpret the mean, median, range, quantiles, IQR, variance and standard deviation in R, and explain which measure is appropriate for a given shape.
3. Build frequency and density histograms, and explain why unequal bin widths force you to use density.
4. Plot and read a cumulative relative frequency curve, including reading quantiles directly off it.
5. Construct boxplots, apply the 1.5 × IQR rule, and compare groups side by side.
6. Recognise when a numerical summary is hiding something a plot would reveal.

1.3 The datasets we will use

Everything in this lecture uses data that is already installed with R. Type the name and it is there — no files, no working directory, no import errors.

# Nothing to install, nothing to load.
str(airquality)
head(faithful)
?airquality        # every built-in dataset has a help page

1.4 Discrete vs continuous: why the tools change

In Week 2 you visualised discrete data with bar charts, because the variable took a manageable number of distinct values and you could count how often each one occurred.

A continuous variable does not work that way. Measure fuel use, ozone concentration or waiting time finely enough and almost every observation is unique, so every "bar" has a count of one. A bar chart of continuous data is just a picture of your measurement precision.

2.1 The stem-and-leaf idea

A stem-and-leaf plot splits every number into two parts. The leading digits become the stem; the final digit becomes the leaf. Leaves are then listed against their stem.

Turn the plot on its side and the rows form bars — it is a histogram whose bin widths are powers of ten. The difference is that no information is thrown away: you can read every original observation back out of the plot.

2.2 stem() in R

Base R provides stem(). It prints to the console rather than drawing a graphic, so it is text you can copy into a report.

stem(airquality$Temp)

  The decimal point is at the |

  56 | 0000
  58 | 0000
  60 | 000
  62 | 000
  64 | 0000
  66 | 0000000
  68 | 0000000
  70 | 0000
  72 | 00000000
  74 | 00000000
  76 | 0000000000000000
  78 | 000000000000
  80 | 0000000000000000
  82 | 0000000000000
  84 | 0000000000
  86 | 000000000000
  88 | 00000
  90 | 00000
  92 | 00000000
  94 | 00
  96 | 00

Temperature is recorded in whole degrees Fahrenheit, so every leaf is 0. R has chosen stems two units wide: the row 76 holds both 76° and 77°.

2.3 Reading the plot

Read the shape from the lengths of the rows, exactly as you would read bar lengths.

stem(airquality$Temp, scale = 0.5)  # fewer, wider stems
stem(airquality$Temp, scale = 2)    # more, narrower stems

Knowledge check — Section 2

3.1 Population and sample notation

Two sets of symbols describe the same quantities, depending on whether you have measured everyone or only a sample. Getting them straight now saves confusion later.

3.2 Mean and median

Finding the median by hand

Sort the data first, then:

• If $n$ is odd, the median is the $\frac{n+1}{2}$-th smallest value.
• If $n$ is even, it is the average of the $\frac{n}{2}$-th and $\left(\frac{n}{2}+1\right)$-th smallest values.

x <- c(31, 34, 35, 36, 38, 40, 41)
mean(x)      # 36.42857
median(x)    # 36

3.3 One outlier, two very different answers

Take seven observations and change only the largest one. The mean moves by more than fifty units; the median does not move at all.

a <- c(31, 34, 35, 36, 38, 40, 41)
b <- c(31, 34, 35, 36, 38, 40, 410)

mean(a);   median(a)     # 36.42857   36
mean(b);   median(b)     # 89.14286   36

Neither is wrong. If you are budgeting total fuel, insurance payouts or total revenue, you need the mean, because the mean multiplied by $n$ gives the total. If you are describing a typical case, the median is usually the honest choice.

3.4 The same effect in real data

Daily ozone in airquality is strongly right-skewed: most days are low, a handful are very high. The mean sits well above the median as a result.

oz <- airquality$Ozone[!is.na(airquality$Ozone)]

length(oz)    # 116
mean(oz)      # 42.12931
median(oz)    # 31.5

The mean of 42.1 exceeds the median of 31.5 by a third. That gap is the skew, expressed numerically. For a symmetric distribution the two would nearly coincide.

3.5 summary() and the missing-value trap

summary() gives five numbers plus the mean in one call — and, crucially, tells you how many values are missing.

summary(airquality$Ozone)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's
   1.00   18.00   31.50   42.13   63.25  168.00      37

Thirty-seven of the 153 days have no ozone reading. Most R summary functions refuse to guess what to do about that:

mean(airquality$Ozone)                 # NA
mean(airquality$Ozone, na.rm = TRUE)   # 42.12931

Knowledge check — Section 3

4.1 Range, and why it is fragile

The range is the distance from the smallest observation to the largest. It is the simplest measure of spread and the least reliable.

oz <- airquality$Ozone[!is.na(airquality$Ozone)]

range(oz)               # 1 168   <- returns min and max, NOT the difference
max(oz) - min(oz)       # 167     <- the range as a single number

The range depends on exactly two observations — the two most extreme, which are also the two most likely to be errors, and the two whose values are least reproducible if you collected the data again. Every other observation is ignored.

We need a measure of spread that uses the middle of the data instead. That is what quantiles give us.

4.2 Quantiles

A quantile is a cut point that divides sorted data into groups of specified proportion. The p-th quantile is the value below which a proportion p of the data lies.

Quartiles split the data into four equal parts; deciles into ten; percentiles into a hundred. They are all the same idea at different resolutions.

4.3 The interquartile range

The interquartile range is the width of the middle half of the data:

Unlike the range, it is computed from the bulk of the data and is completely unaffected by however extreme the most extreme observations happen to be. Move the largest value from 168 to 168,000 and the IQR does not change.

oz <- airquality$Ozone[!is.na(airquality$Ozone)]

quantile(oz, 0.25)      # 18.00
quantile(oz, 0.75)      # 63.25
IQR(oz)                 # 45.25

# equivalently
quantile(oz, 0.75) - quantile(oz, 0.25)

So the central half of the ozone readings spans 45.25 ppb, from 18 to 63.25. That is a far more useful description of typical variability than the range of 167, which is driven entirely by one exceptional day.

4.4 How R actually computes a quantile

When $p$ does not land exactly on an observation, R interpolates between the two neighbouring values. Its default method computes a position

then takes the $h$-th value in the sorted data, interpolating linearly when $h$ is not a whole number.

x <- c(3, 7, 8, 12, 15, 21, 24, 30)
quantile(x, 0.25)       # 7.75  -- matches the hand calculation

Nine different quantile definitions exist and textbooks disagree about which to use; some define the position as $p(n+1)$ instead, giving 7.25 here. R implements all nine via the type argument. Small differences between software packages usually trace back to this choice, not to an error.

4.5 Quantiles in R

Called with no probabilities, quantile() returns the five-number summary. Pass a vector of probabilities to get exactly the cut points you want.

oz <- airquality$Ozone[!is.na(airquality$Ozone)]

quantile(oz)

    0%    25%    50%    75%   100%
  1.00  18.00  31.50  63.25 168.00

# any probabilities you like
quantile(oz, c(0.10, 0.90))

 10%  90%
 7.5 87.0

So 10% of days had ozone below 7.5 ppb and 10% were above 87 ppb. Note quantile() has no na.rm default — pass na.rm = TRUE if your vector still contains missing values.

4.6 Variance, step by step

Spread means distance from the centre. The obvious measure — average deviation from the mean — is useless, because the deviations always sum to exactly zero. Squaring them first solves that.

Take $x=2,4,4,4,5,5,7,9$, so $n=8$ and $\bar{x}=5$.

4.7 var() and sd()

Confirm the hand calculation, then use the built-in functions from here on.

x <- c(2, 4, 4, 4, 5, 5, 7, 9)
n <- length(x)

# the long way
d  <- x - mean(x)          # deviations
sum(d)                     # 0        <- always
sum(d^2)                   # 32
sum(d^2) / (n - 1)         # 4.571429

# the short way
var(x)                     # 4.571429
sd(x)                      # 2.13809
sqrt(var(x))               # 2.13809  <- sd is just the square root of the variance

4.8 Why divide by $n-1$?

Four reasons, in increasing order of depth:

1. It is the standard convention. Every downstream statistical procedure — confidence intervals, t-tests, regression — assumes the divisor was $n-1$.
2. Software agrees. R's var() and sd() use $n-1$. Excel's VAR.S does too, though VAR.P does not — a classic source of mismatched answers between the two.
3. It behaves sensibly at $n=1$. One observation tells you nothing about variability, and $0/0$ correctly yields an undefined answer rather than a falsely confident zero.
4. It is unbiased. Imagine drawing sample after sample from a population. Dividing by $n$ gives variances that are, on average, slightly too small, because the sample mean sits closer to your own sample than the true $\mu$ does. Dividing by $n-1$ corrects exactly this, so the long-run average of $s^2$ equals ${\sigma }^2$.

4.9 Interpreting the standard deviation

The standard deviation is in the same units as the data, which makes it directly interpretable as a typical distance from the mean. For data with a roughly bell-shaped distribution, three benchmarks are worth memorising.

4.10 Does the rule survive contact with real data?

Here is the proportion of observations actually falling within one standard deviation of the mean, for five built-in datasets. The bell-curve prediction is 68%.

v <- rivers
mean(abs(v - mean(v)) <= sd(v))     # 0.9007092

The rule works well for the two datasets closest to symmetric and fails badly for rivers, where a handful of enormous values inflate the standard deviation so much that 90% of rivers fall inside one sd of the mean.

4.11 Standard deviation or IQR?

Both measure spread. They answer slightly different questions and fail in different circumstances.

Because $\pm 1$ sd spans about two-thirds of the data while the IQR spans exactly one half, you should normally expect $2\times \text{sd}$ to exceed the IQR. It does for all five datasets on the previous slide.

4.12 Shifting and scaling

Changing units is a linear transformation, and its effect on each statistic is predictable.

Shifting slides the whole distribution along without stretching it, so spread is untouched. Scaling stretches it, and because variance is in squared units it scales by $k^2$.

tempF <- airquality$Temp
tempC <- (tempF - 32) * 5/9        # shift by -32, then scale by 5/9

mean(tempF)   # 77.88235      mean(tempC)   # 25.49020
sd(tempF)     #  9.46527      sd(tempC)     #  5.25848
var(tempF)    # 89.59133      var(tempC)    # 27.65165

sd(tempF) * 5/9         # 5.25848   <- matches sd(tempC)
var(tempF) * (5/9)^2    # 27.65165  <- matches var(tempC)

The −32 shift changed the mean but had no effect at all on the sd or the variance. Only the ×5/9 scaling did.

Knowledge check — Section 4

5.1 What a histogram does

A histogram divides the range of the data into intervals called bins, counts the observations in each, and draws a bar for each bin. Bars touch, because the underlying variable is continuous and the bins are adjacent intervals rather than separate categories.

A bar chart and a histogram look similar and mean different things. A bar chart's horizontal axis lists categories in whatever order you choose; a histogram's is a number line, so bar order and bar width both carry meaning.

5.2 hist() in R

oz <- airquality$Ozone[!is.na(airquality$Ozone)]

hist(oz,
     xlab = "Ozone (ppb)",
     main = "Ozone concentration, 116 days")

Always label the axes and give the plot a title. xlab, ylab and main take character strings; the default labels are the R expression you passed in, which is meaningless to anyone reading your report.

5.3 Frequency and density

hist() can plot two different things on the vertical axis.

hist(oz, freq = FALSE, xlab = "Ozone (ppb)", main = "")

# verify the total area is 1
h <- hist(oz, plot = FALSE)
sum(h$density * diff(h$breaks))     # 1

5.4 Unequal bins: frequency misleads

Suppose we bin ozone finely at the low end and coarsely at the high end. The 50–100 bin is five times wider than its neighbours, so it accumulates a large count purely by being wide.

brk <- c(0, 10, 20, 30, 40, 50, 100, 200)

par(mfrow = c(1, 2))
hist(oz, breaks = brk, freq = TRUE,  xlab = "Ozone", main = "(a) frequency")
hist(oz, breaks = brk, freq = FALSE, xlab = "Ozone", main = "(b) density")

In panel (a) the highlighted 50–100 bin looks like the second-largest group in the dataset. In panel (b) the same bin is one of the shortest, because its 27 observations are spread thinly across a 50-unit interval.

5.5 Choosing the number of bins

The breaks argument controls bin count. R treats a single number as a suggestion and rounds to convenient boundaries.

par(mfrow = c(1, 3))
hist(faithful$waiting, breaks = 5,  main = "too few")
hist(faithful$waiting, breaks = 12, main = "about right")
hist(faithful$waiting, breaks = 60, main = "too many")

With five bins the two peaks are smeared into one broad lump. With twelve they are unmistakable. With sixty, random noise in individual bins starts to look like structure.

5.6 The vocabulary of shape

Four descriptions cover most of what you will need to say about a distribution.

Skew is named for the direction of the tail, not the location of the peak. rivers is right-skewed because its long thin tail extends to the right, even though its tallest bar is on the left.

5.7 When the summary hides the shape

faithful$waiting records the gap in minutes between eruptions of the Old Faithful geyser. Its numerical summary looks entirely unremarkable.

summary(faithful$waiting)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
   43.0    58.0    76.0    70.9    82.0    96.0

5.8 Heavy tails and the log transform

When one tail is extreme, a histogram on the original scale becomes a single tall bar next to a strip of near-empty space. Plotting the logarithm of the values spreads them out.

par(mfrow = c(1, 2))
hist(rivers,        xlab = "length (miles)",       main = "original scale")
hist(log10(rivers), xlab = "log10(length)",         main = "log scale")

On the log scale the distribution is close to symmetric, and structure that was compressed into the first bar becomes visible. This is why financial returns, city populations, incomes and file sizes are so often plotted on log axes.

The transformed values are logarithms, so a bar at 2.6 means $10^{2.6}\approx 400$ miles. Label such axes carefully or your reader will misread them by orders of magnitude.

Knowledge check — Section 5

6.1 The idea

Instead of asking how many observations fall in an interval, ask how many fall at or below a value. Plot that proportion against the value and you get the empirical cumulative distribution function.

The curve starts at 0 on the left, ends at 1 on the right, and never decreases. It steps up by $1/n$ at every observation, which is why it looks like a staircase.

6.2 ecdf() in R

oz <- airquality$Ozone[!is.na(airquality$Ozone)]

plot(ecdf(oz),
     verticals = TRUE,
     main = "Cumulative relative frequency",
     xlab = "Ozone (x, ppb)",
     ylab = "Proportion less than or equal to x",
     ylim = c(0, 1))

ecdf() returns a function, so you can also evaluate it directly rather than only plotting it:

F <- ecdf(oz)
F(50)     # 0.7068966  -> about 71% of days were at or below 50 ppb
F(100)    # 0.9396552  -> about 94% were at or below 100 ppb

6.3 Reading quantiles off the curve

The ECDF makes quantiles geometric. To find the median, start at 0.5 on the vertical axis, travel right to the curve, then drop to the horizontal axis.

plot(ecdf(oz), verticals = TRUE, main = "", xlab = "Ozone (ppb)",
     ylab = "Proportion <= x")

q <- quantile(oz, c(0.25, 0.5, 0.75))
abline(h = c(0.25, 0.5, 0.75), lty = 2, col = "red")   # horizontal guides
abline(v = q,                  lty = 2, col = "red")   # vertical guides

The steepness of the curve is itself informative: the sharp rise between 10 and 40 ppb tells you most days are concentrated there, while the near-flat stretch beyond 100 tells you high-ozone days are rare.

6.4 Two plots in one figure

par(mfrow = c(rows, cols)) splits the plotting area into a grid, and subsequent plots fill it row by row. Stacking a histogram above its ECDF lets you match features between the two.

par(mfrow = c(2, 1))                       # 2 rows, 1 column

hist(oz, freq = FALSE, xlim = c(0, 180),
     xlab = "", ylab = "Density", main = "Ozone: density and cumulative views")

plot(ecdf(oz), verticals = TRUE, xlim = c(0, 180),
     xlab = "Ozone (x, ppb)", ylab = "Proportion <= x", main = "")

par(mfrow = c(1, 1))                       # reset when finished

Where the histogram has a tall bar, the ECDF is steep. Where the histogram is flat, the ECDF is nearly horizontal. They are two views of the same information.

Knowledge check — Section 6

7.1 The five-number summary

A boxplot condenses a distribution into five numbers, then draws them.

The box therefore spans the IQR and contains the central half of the data. The position of the median line within the box indicates skew: pushed towards the bottom means a longer upper tail.

7.2 Anatomy of a boxplot

Reading this plot: the median sits well below the centre of the box, and the upper whisker is much longer than the lower one. Both features say the same thing — the distribution is right-skewed, exactly as the histogram in Section 3 showed.

7.3 The 1.5 × IQR rule

Whiskers do not simply run to the minimum and maximum. R first computes two fences:

Each whisker then extends only as far as the most extreme observation still inside its fence. Anything beyond is drawn as an individual point.

boxplot.stats(oz)$out     # 135 168

7.4 boxplot() in R

boxplot(oz, ylab = "Ozone (ppb)", main = "Daily ozone, New York 1973")

# horizontal often reads better under a histogram
boxplot(oz, horizontal = TRUE, xlab = "Ozone (ppb)")

# the underlying numbers
boxplot.stats(oz)

$stats
[1]   1.0  18.0  31.5  63.5 122.0

$n
[1] 116

$out
[1] 135 168

The $stats vector holds the five plotted values in order: lower whisker, Q1, median, Q3, upper whisker. $out lists the flagged points and $n the sample size.

The Q3 reported here is 63.5, while quantile(oz, 0.75) gave 63.25. boxplot() uses hinges rather than R's default quantile definition — the same nine-definitions issue from Section 4.4. The difference is never material, but it is worth recognising rather than hunting for a bug.

7.5 Comparing groups side by side

This is what boxplots are genuinely best at. The formula interface y ~ g reads as "y broken down by g" and produces one box per group.

boxplot(weight ~ feed, data = chickwts,
        ylab = "Weight (g)", xlab = "Feed supplement",
        main = "Chick weight by feed type")

Casein and sunflower produce the heaviest chicks; horsebean is clearly worst. Sunflower also has the tightest box, meaning the most consistent results — a comparison of spread that a bar chart of group means would have concealed entirely.

7.6 What boxplots hide

A boxplot shows five numbers. Anything not captured by those five numbers is invisible — and that includes the number of peaks.

Both panels show faithful$waiting. The boxplot is entirely well behaved: a symmetric-looking box, whiskers reaching the extremes, and not a single flagged outlier. The histogram beside it shows two separated groups.

Knowledge check — Section 7

8.1 Which plot, and when

8.2 Function reference

8.3 Where this goes next