Descriptive statistics

Describing distributions: centre and spread

Dr Massimiliano Canzi

Distributions

Did you know: tardigrades are animals about 0.5 mm long, with eight stubby legs. The German pastor J.A.E. Goeze first described them in 1773 and called them kleiner Wasserbär: “little water bear”.

Source: National Geographic

Why describe?

A data set with 344 penguins is too much to take in at once. A descriptive statistic summarises it in one number, such as a typical value or how much the values vary.

Like a map, a summary leaves most of the detail out. Today is about choosing a good summary, and noticing when it hides something important.

Today

What a distribution is: empirical and theoretical
Measures of central tendency: mean, median, mode
Measures of spread: standard deviation, MAD, range, quantiles
When each map misleads: skew, outliers, more than one peak
Grouped summaries and plots, building on {dplyr}

Distributions

A distribution shows how often each value, or range of values, occurs.

An empirical distribution describes the values we actually observed
A theoretical distribution is a mathematical shape, fully described by a few numbers called parameters

We compare the two all the time: does our sample look like a shape we know?

Plotting with `{ggplot2}`

Last week’s promise: {ggplot2} builds a plot in layers, added with +.

ggplot(data, aes(x = column)) +   # the data, and which column goes where
  geom_histogram() +              # how to draw it
  labs(x = "Axis label")          # labels

ggplot() and aes() say which data to use and which column goes on x, y or colour
A geom_*() draws the shapes: geom_histogram(), geom_point(), geom_boxplot()
labs() and theme() change the labels and the look

Our slides use a theme set up in the hidden first chunk of the .Rmd: your plots will look different, but show the same data.

An empirical distribution

ggplot(penguins, aes(x = flipper_length_mm)) +
  geom_histogram(binwidth = 3, fill = "#5aa0e0", colour = "#0b110e") +
  labs(x = "Flipper length (mm)", y = "Count")

Each bar counts the penguins whose flipper length falls in a 3 mm bin.

Histogram or density?

ggplot(penguins, aes(x = flipper_length_mm)) +
  geom_histogram(aes(y = after_stat(density)), binwidth = 3,
                 fill = "#5aa0e0", colour = "#0b110e", alpha = 0.5) +
  geom_density(colour = "#f2d45c", linewidth = 1.2) +
  labs(x = "Flipper length (mm)", y = "Density")

The density curve is a smoothed histogram. The area under it is 1.

The normal distribution

In a normal (or Gaussian) distribution, values are spread symmetrically around the mean, and most of them lie close to it: the familiar bell shape.

Two parameters describe it completely: the mean and the standard deviation.

set.seed(42)    # same random numbers every time we knit
normal <- rnorm(500, mean = 0, sd = 1)
head(round(normal, 2))

[1]  1.37 -0.56  0.36  0.63  0.40 -0.11

Sample vs theory

tibble(value = normal) |>
  ggplot(aes(x = value)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30,
                 fill = "#e08ab8", colour = "#0b110e", alpha = 0.6) +
  stat_function(fun = dnorm, colour = "#f2d45c", linewidth = 1.2) +
  labs(x = "Value", y = "Density")

Bars: our 500 random values. Line: the theoretical normal curve.

The 68–95–99.7 rule

In a normal distribution, about 68% of values lie within 1 SD of the mean, 95% within 2 SD, 99.7% within 3 SD.

Other distributions

Not every variable is normal. Three shapes you will meet in linguistic data:

Accuracy out of a fixed number of trials is binomial (lecture 12). Counts of events, such as pauses per minute, are often Poisson. Reaction times and durations are positive with a long right tail, closer to a log-normal or gamma shape.

Central tendency

Did you know: a protein found only in tardigrades, Dsup, protects DNA from radiation. Human cells given the gene for it showed about 40% less DNA damage after X-rays.

Source: Hashimoto et al. (2016), Nature Communications

The tools of the trade

Descriptive measures come in two kinds:

Central tendency	Spread
`mean()`	`sd()`, `var()`
`median()`	`mad()`, `IQR()`
mode	`range()`, `quantile()`

Always report one of each: a measure of centre and a measure of spread.

First: missing values

mean(penguins$bill_length_mm)

[1] NA

mean(penguins$bill_length_mm, na.rm = TRUE)

[1] 43.92193

NA means “unknown”. The mean of numbers that include an unknown is unknown, unless you tell R to drop them.

Three ways to drop missing values, three different data sets:

nrow(penguins)                               # all rows

[1] 344

nrow(drop_na(penguins, bill_length_mm))      # NA in bill length only

[1] 342

nrow(na.omit(penguins))                      # NA in any column

[1] 333

`na.omit()` changes your result

mean(penguins$bill_length_mm, na.rm = TRUE)   # 342 penguins

[1] 43.92193

mean(na.omit(penguins)$bill_length_mm)        # 333 penguins

[1] 43.99279

na.omit() also dropped 9 penguins with a known bill length but unknown sex.

Rule of thumb: drop NAs only in the columns you are analysing, with na.rm = TRUE or drop_na(column).

`mean()`

The mean is the sum of all values divided by their number

x = (x₁ + x₂ + … + xₙ) / n

x <- c(4, 8, 15, 16, 23, 42)
sum(x) / length(x)

[1] 18

mean(x)

[1] 18

`mean()`

penguins |>
  summarise(mean_bill = mean(bill_length_mm, na.rm = TRUE))

# A tibble: 1 × 1
  mean_bill
      <dbl>
1      43.9

`median()`

The median is the middle value once all values are sorted from lowest to highest: half the values lie below it, half above.

median(c(3, 1, 5))        # sorted: 1 3 5

[1] 3

median(c(3, 1, 5, 7))     # sorted: 1 3 5 7 -> mean of 3 and 5

[1] 4

With an even number of values, the median is the mean of the two middle values.

`median()`

penguins |>
  summarise(mean_bill   = mean(bill_length_mm, na.rm = TRUE),
            median_bill = median(bill_length_mm, na.rm = TRUE))

# A tibble: 1 × 2
  mean_bill median_bill
      <dbl>       <dbl>
1      43.9        44.4

One number can hide a lot

The overall mean (43.9 mm) sits in a dip: few penguins actually have a bill that long. The distribution has three species in it.

Skewed data: reaction times

Reaction times cannot be negative, and some responses are very slow:

set.seed(3)
rt <- tibble(rt = rlnorm(2000, meanlog = log(600), sdlog = 0.35))
rt |> summarise(mean = mean(rt), median = median(rt))

# A tibble: 1 × 2
   mean median
  <dbl>  <dbl>
1  636.   602.

Skewed data: reaction times

A long right tail pulls the mean (dashed) away from the median (dotted).

Outliers

One participant fell asleep on one trial:

set.seed(7)
rts <- round(rlnorm(40, meanlog = log(600), sdlog = 0.25))
rts_sleepy <- c(rts, 5000)

tibble(data   = c("40 trials", "40 trials + one 5000 ms"),
       mean   = c(mean(rts), mean(rts_sleepy)),
       median = c(median(rts), median(rts_sleepy)))

# A tibble: 2 × 3
  data                     mean median
  <chr>                   <dbl>  <dbl>
1 40 trials                665.   626.
2 40 trials + one 5000 ms  770.   628

The mean moves by more than 100 ms. The median barely moves: it is robust.

Discussion

Here we know that the participant fell asleep. In real data, a 5,000 ms response usually comes with no explanation.

Would you remove it? How would you decide which responses are too slow?
Why does it matter whether you decide before or after looking at the results?

Mode

The mode is the most frequent value. It is used in categorical distributions.

With count() from last week:

penguins |> count(island) |> slice_max(n)

# A tibble: 1 × 2
  island     n
  <fct>  <int>
1 Biscoe   168

words <- tibble(word = c("the", "cat", "saw", "the", "dog", "and", "the", "bird"))
words |> count(word, sort = TRUE) |> head(3)

# A tibble: 3 × 2
  word      n
  <chr> <int>
1 the       3
2 and       1
3 bird      1

Mode: a trap

R has a function called mode(), but it does not compute the mode:

mode(c("Hello", "Hello", "Goodbye"))

[1] "character"

It returns the storage type. If you need a function, write your own:

find_mode <- function(x) {
  u <- unique(x)
  tab <- tabulate(match(x, u))
  u[tab == max(tab)]
}
find_mode(c("Hello", "Hello", "Goodbye", "Hello", "Penguin", "Dog", "Mountain"))

[1] "Hello"

Spread

Did you know: when their home dries out, tardigrades lose their body water, curl up and enter cryptobiosis, a deep sleep. They can revive years later, when water returns.

Sources: National Geographic; ScienceDaily

Same centre, different spread

Both speakers have a mean of 150 ms. Speaker B is far less consistent.

`sd()`

The standard deviation measures how spread out the values are around the mean: roughly, how far a typical value is from it.

s = √( Σ (xᵢ − x)² / (n − 1) )

x <- c(4, 8, 15, 16, 23, 42)
sqrt(sum((x - mean(x))^2) / (length(x) - 1))

[1] 13.49074

sd(x)

[1] 13.49074

var(x)       # the variance is the SD squared

[1] 182

`sd()`

penguins |>
  summarise(mean_bill = mean(bill_length_mm, na.rm = TRUE),
            sd_bill   = sd(bill_length_mm, na.rm = TRUE))

# A tibble: 1 × 2
  mean_bill sd_bill
      <dbl>   <dbl>
1      43.9    5.46

`mad()`

The median absolute deviation (MAD) is the robust partner of the SD: the median distance of the values from their median.

Step	Values
1. The data, sorted	1, 1, 2, 2, 4, 6, 9	median = 2
2. Distance of each value from 2	1, 1, 0, 0, 2, 4, 7
3. The same distances, sorted	0, 0, 1, 1, 2, 4, 7	MAD = 1

`mad()` in R

x <- c(1, 1, 2, 2, 4, 6, 9)
mad(x, constant = 1)

[1] 1

mad(x)

[1] 1.4826

By default, mad() multiplies the result by 1.4826. With normal data, this makes the MAD estimate the SD, so the two are comparable.

c(sd = sd(rts), mad = mad(rts))                    # 40 trials

      sd      mad 
187.2306 150.4839

c(sd = sd(rts_sleepy), mad = mad(rts_sleepy))      # + one 5000 ms trial

      sd      mad 
701.8433 158.6382

`mad()`

penguins |>
  summarise(median_bill = median(bill_length_mm, na.rm = TRUE),
            mad_bill    = mad(bill_length_mm, na.rm = TRUE))

# A tibble: 1 × 2
  median_bill mad_bill
        <dbl>    <dbl>
1        44.4     7.04

`range()`

range() returns the smallest and the largest value. The distance between them is what statisticians call the range.

range(penguins$bill_length_mm, na.rm = TRUE)

[1] 32.1 59.6

diff(range(penguins$bill_length_mm, na.rm = TRUE))

[1] 27.5

The range depends on two values only, the most extreme ones. One outlier changes it completely.

Quantiles / Quartiles

Quantiles cut the sorted values into groups of equal size. Quartiles cut them into four quarters:

quantile(penguins$bill_length_mm, na.rm = TRUE)

    0%    25%    50%    75%   100% 
32.100 39.225 44.450 48.500 59.600

IQR(penguins$bill_length_mm, na.rm = TRUE)     # 75% minus 25%

[1] 9.275

The median is the 50% quantile. The IQR holds the middle half of the data.

`summary()`

Most of today in one line:

summary(penguins$bill_length_mm)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
  32.10   39.23   44.45   43.92   48.50   59.60       2

It works on a whole data frame too: try summary(penguins).

Boxplots draw the quartiles

Box: 25% to 75% (the IQR). Line: the median. Whiskers: to the furthest point within 1.5 × IQR of the box. Dots: every penguin.

`group_by` and `summarise`

Did you know: in 2007, dried tardigrades flew into orbit on the FOTON-M3 mission and were exposed to open space: vacuum, cosmic rays and solar UV over 1,000 times stronger than on Earth. Most survived the vacuum, a few even the full UV, and the survivors reproduced normally.

Sources: Jönsson et al. (2008), Current Biology; ScienceDaily

Summaries per group

Until now, we have summarised one distribution at a time.

Usually we want the same summary for each group: bill_length_mm for each species separately.

The bad way to do it:

Adelie    <- filter(penguins, species == "Adelie")
Gentoo    <- filter(penguins, species == "Gentoo")
Chinstrap <- filter(penguins, species == "Chinstrap")
mean(Adelie$bill_length_mm, na.rm = TRUE)

[1] 38.79139

mean(Gentoo$bill_length_mm, na.rm = TRUE)

[1] 47.50488

mean(Chinstrap$bill_length_mm, na.rm = TRUE)

[1] 48.83382

The good way

Last week: group_by() + summarise(), or summarise(.by = )

penguins |>
  drop_na(bill_length_mm) |>
  summarise(n      = n(),
            mean   = mean(bill_length_mm),
            sd     = sd(bill_length_mm),
            median = median(bill_length_mm),
            IQR    = IQR(bill_length_mm),
            .by    = species)

# A tibble: 3 × 6
  species       n  mean    sd median   IQR
  <fct>     <int> <dbl> <dbl>  <dbl> <dbl>
1 Adelie      151  38.8  2.66   38.8  4   
2 Gentoo      123  47.5  3.08   47.3  4.25
3 Chinstrap    68  48.8  3.34   49.6  4.73

One row per species, one column per summary, no copies of the data.

Grouping by two variables

penguins |>
  drop_na(sex) |>
  group_by(species, sex) |>
  summarise(mean_bill = mean(bill_length_mm),
            median_bill = median(bill_length_mm))

`summarise()` has regrouped the output.
ℹ Summaries were computed grouped by species and sex.
ℹ Output is grouped by species.
ℹ Use `summarise(.groups = "drop_last")` to silence this message.
ℹ Use `summarise(.by = c(species, sex))` for per-operation grouping
  (`?dplyr::dplyr_by`) instead.

# A tibble: 6 × 4
# Groups:   species [3]
  species   sex    mean_bill median_bill
  <fct>     <fct>      <dbl>       <dbl>
1 Adelie    female      37.3        37  
2 Adelie    male        40.4        40.6
3 Chinstrap female      46.6        46.3
4 Chinstrap male        51.1        51.0
# ℹ 2 more rows

The message: the result is still grouped by species. Either add .groups = "drop", ungroup() it, or use .by = c(species, sex), which never leaves groups behind.

`group_by()` or `.by`?

Use group_by() when several steps in a row should work per group, and ungroup() when you are done, before you save the result.

Use .by when only one step needs the groups. Nothing stays grouped afterwards, so there is nothing to forget.

Many columns at once: `across()`

From last week: across(columns, function)

penguins |>
  summarise(across(ends_with("_mm"),
                   list(mean = \(x) mean(x, na.rm = TRUE),
                        sd   = \(x) sd(x, na.rm = TRUE))),
            .by = species) |>
  select(species, starts_with("bill_length"))

# A tibble: 3 × 3
  species   bill_length_mm_mean bill_length_mm_sd
  <fct>                   <dbl>             <dbl>
1 Adelie                   38.8              2.66
2 Gentoo                   47.5              3.08
3 Chinstrap                48.8              3.34

list(mean = ..., sd = ...) applies two functions to each column, and adds _mean and _sd to the names.

Plotting group summaries

bill_summary <- penguins |>
  drop_na(bill_length_mm) |>
  summarise(mean = mean(bill_length_mm), sd = sd(bill_length_mm), .by = species)

ggplot(drop_na(penguins, bill_length_mm), aes(x = species, y = bill_length_mm)) +
  geom_jitter(aes(colour = species), width = 0.2, alpha = 0.3) +
  geom_pointrange(data = bill_summary, colour = "#e3e8e1", linewidth = 1, size = 0.8,
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd)) +
  labs(x = NULL, y = "Bill length (mm)") + theme(legend.position = "none")

Points: every penguin. White: mean ± 1 SD per species.

Reporting descriptives

In R Markdown, numbers in your text can come straight from R:

gentoo <- filter(bill_summary, species == "Gentoo")

Gentoo bills were `r round(gentoo$mean, 1)` mm long on average
(SD = `r round(gentoo$sd, 1)`).

becomes

Gentoo bills were 47.5 mm long on average (SD = 3.1).

If the data change, the text updates itself, so you never copy numbers by hand.

Summary

A summary is a simplified map of a distribution
Centre: mean() (sensitive), median() (robust), mode with count()
Spread: sd() (sensitive), mad() and IQR() (robust), range(), quantile()
Always plot the data: skew, outliers and several peaks hide behind one number
Per group: summarise() with .by or group_by(); many columns with across()

Slides and code

These slides

mcanzi.de/ws26_statistics/03_descriptive_statistics.html

Code, data and all lectures

github.com/mcanzi/ws26_statistics

Descriptive statistics

Describing distributions: centre and spread

Dr Massimiliano Canzi

Distributions

Why describe?

Today

Distributions

Plotting with {ggplot2}

An empirical distribution

Histogram or density?

The normal distribution

Sample vs theory

The 68–95–99.7 rule

Other distributions

Central tendency

The tools of the trade

First: missing values

na.omit() changes your result

mean()

mean()

median()

median()

One number can hide a lot

Skewed data: reaction times

Skewed data: reaction times

Outliers

Discussion

Mode

Mode: a trap

Spread

Same centre, different spread

sd()

sd()

mad()

mad() in R

mad()

range()

Quantiles / Quartiles

summary()

Boxplots draw the quartiles

group_by and summarise

Summaries per group

The good way

Grouping by two variables

group_by() or .by?

Many columns at once: across()

Plotting group summaries

Reporting descriptives

Summary

Slides and code

Plotting with `{ggplot2}`

`na.omit()` changes your result

`mean()`

`mean()`

`median()`

`median()`

`sd()`

`sd()`

`mad()`

`mad()` in R

`mad()`

`range()`

`summary()`

`group_by` and `summarise`

`group_by()` or `.by`?

Many columns at once: `across()`