The arithmetic mean

The mean is the balance point of a distribution — the value that makes all deviations sum to zero. You already know how to compute it, but it is worth being precise about what it measures and when it misleads.

\[\bar{x} = \frac{1}{n} \sum_{i=1}^{n} x_i\]

Let’s use a real dataset: hourly PM2.5 concentrations (µg/m³) recorded at a Southern California monitoring station over one week.

# Weekly PM2.5 readings (µg/m³) — one value per day
pm25 <- c(8.2, 11.4, 9.7, 42.3, 10.1, 9.3, 8.8)
days <- c("Mon", "Tue", "Wed", "Thu", "Fri", "Sat", "Sun")

n         <- length(pm25)
pm25_mean <- mean(pm25)

With n = 7 days, the mean PM2.5 concentration is 14.3 µg/m³.

Thursday’s value (42.3 µg/m³) looks suspicious — possibly a smoke event. Notice how strongly it pulls the mean upward.

Daily PM2.5 values. The dashed red line shows the mean; the solid blue line shows the median. One high-smoke Thursday pulls the mean far above the typical week.

The median

The median is the middle value when data are sorted. It is resistant to extreme values — one wildfire day does not distort it the way it distorts the mean.

pm25_median <- median(pm25)

The median is 9.7 µg/m³ — a gap of 4.6 µg/m³ below the mean, entirely due to Thursday’s spike.

Checking both together

# summary() gives the five-number summary plus the mean
summary(pm25)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    8.20    9.05    9.70   14.26   10.75   42.30

Residuals: deviations from the mean

The spread of data is fundamentally about how far each observation is from the mean. These deviations are called residuals.

residuals_pm25 <- pm25 - pm25_mean
round(residuals_pm25, 1)

## [1] -6.1 -2.9 -4.6 28.0 -4.2 -5.0 -5.5

Residuals from the mean. Each bar shows how far that day’s PM2.5 was from the weekly mean. Notice the enormous positive residual on Thursday.

Notice that the residuals sum to zero — that is always true by definition of the mean:

round(sum(residuals_pm25), 10)   # essentially zero (floating-point rounding)

## [1] 0

Why we square the residuals

We cannot use the raw sum of residuals as a measure of spread (it is always zero). Instead, we square each residual, which makes all values positive and gives larger weight to observations far from the mean.

sq_resid <- residuals_pm25^2
round(sq_resid, 1)

## [1]  36.7   8.2  20.8 786.4  17.3  24.6  29.8

Sum of Squares, Variance, and SD

\[\text{SS} = \sum_{i=1}^{n}(x_i - \bar{x})^2\]

\[s^2 = \frac{\text{SS}}{n-1} \qquad \text{(variance — uses } n-1 \text{ for sample)}\]

\[s = \sqrt{s^2} \qquad \text{(standard deviation)}\]

We divide by \(n-1\) rather than \(n\) because we are estimating the population variance from a sample — using \(n\) would systematically underestimate it. This correction is called Bessel’s correction.

SS       <- sum(sq_resid)
variance <- SS / (n - 1)          # sample variance
sd_hand  <- sqrt(variance)        # standard deviation by hand

knitr::kable(
  data.frame(
    Statistic = c("Sum of Squares", "Variance", "SD (by hand)", "sd() check"),
    Value     = round(c(SS, variance, sd_hand, sd(pm25)), 2),
    Units     = c("µg²/m⁶", "µg²/m⁶", "µg/m³", "µg/m³")
  ),
  caption = "Variance decomposition for the weekly PM2.5 data."
)

Two distributions: same mean, different spread

This is the key insight. Here are two hypothetical monitoring stations with the same annual average PM2.5 but very different variability:

set.seed(42)
station_A <- rnorm(200, mean = 15, sd = 3)
station_B <- rnorm(200, mean = 15, sd = 10)

par(las = 1, mar = c(4, 4, 3, 1), mfrow = c(1, 2))

hist(station_A, breaks = 20, col = "#4a7c5e", border = "white",
     main = "Station A\nSD = 3 µg/m³", xlab = "PM2.5 (µg/m³)",
     xlim = c(-15, 55), ylim = c(0, 50))
abline(v = 15, col = "#c8832a", lty = 2, lwd = 2)
abline(v = 35, col = "red",     lty = 2, lwd = 1.5)

hist(station_B, breaks = 20, col = "#8aad96", border = "white",
     main = "Station B\nSD = 10 µg/m³", xlab = "PM2.5 (µg/m³)",
     xlim = c(-15, 55), ylim = c(0, 50))
abline(v = 15, col = "#c8832a", lty = 2, lwd = 2)
abline(v = 35, col = "red",     lty = 2, lwd = 1.5)

Two stations with the same mean (15 µg/m³) but different standard deviations. Station B has much higher day-to-day variability — more frequent exceedance events even though the average is identical.

par(mfrow = c(1, 1))

knitr::kable(
  data.frame(
    Station         = c("Station A (SD = 3)", "Station B (SD = 10)"),
    Days_above_35   = c(sum(station_A > 35), sum(station_B > 35)),
    Total_days      = c(length(station_A), length(station_B)),
    Pct_exceedance  = round(100 * c(mean(station_A > 35), mean(station_B > 35)), 1)
  ),
  col.names = c("Station", "Days > 35 µg/m³", "Total days", "% exceedance"),
  caption = "Exceedance of the 24-hr NAAQS threshold (35 µg/m³) by station."
)

Shape and parameters

A normal distribution is completely described by two numbers: its mean (\(\mu\)) and standard deviation (\(\sigma\)).

\[X \sim \mathcal{N}(\mu,\ \sigma^2)\]

The normal distribution with mean = 0 and SD = 1 (the standard normal). The shaded regions show the 68–95–99.7 rule.

The 68–95–99.7 rule

This is the single most useful fact about the normal distribution:

• 68% of observations fall within 1 SD of the mean
• 95% fall within 2 SD of the mean
  
• 99.7% fall within 3 SD of the mean

# Verify with a large simulated normal sample
set.seed(7)
temps <- rnorm(10000, mean = 18.5, sd = 4.2)   # daily max °C, inland CA

within_1sd <- mean(abs(temps - mean(temps)) < sd(temps))
within_2sd <- mean(abs(temps - mean(temps)) < 2 * sd(temps))
within_3sd <- mean(abs(temps - mean(temps)) < 3 * sd(temps))

knitr::kable(
  data.frame(
    Range      = c("Within ±1 SD", "Within ±2 SD", "Within ±3 SD"),
    Theoretical = c("68.3%", "95.4%", "99.7%"),
    Observed    = paste0(round(c(within_1sd, within_2sd, within_3sd) * 100, 1), "%")
  ),
  caption = "Empirical check of the 68–95–99.7 rule (n = 10,000 simulated temperatures)."
)

Checking normality: the Q-Q plot

A quantile-quantile (Q-Q) plot compares your data to what a perfect normal distribution would look like. Points falling on the diagonal line = normality. Strong departures = skewness or heavy tails.

set.seed(3)
daily_temp <- rnorm(365, mean = 22, sd = 5.5)   # roughly normal

par(mfrow = c(1, 2), las = 1, mar = c(4, 4, 3, 1))

# Q-Q for PM2.5 (skewed by smoke event)
pm25_long <- c(rnorm(90, 10, 3), 42.3, rnorm(9, 10, 3))
qqnorm(pm25_long, main = "PM2.5 (skewed)", pch = 20, cex = 0.6,
       col = "#4a7c5e")
qqline(pm25_long, col = "#c8832a", lwd = 2)

# Q-Q for daily temperature (approximately normal)
qqnorm(daily_temp, main = "Daily Max Temp (≈ normal)", pch = 20, cex = 0.6,
       col = "#4a7c5e")
qqline(daily_temp, col = "#c8832a", lwd = 2)

Q-Q plots for two variables. Left: PM2.5 with a smoke event — the upper tail departs strongly from normality. Right: daily maximum temperature — much closer to the normal line.

par(mfrow = c(1, 1))

The Central Limit Theorem

If you take many samples from any distribution and compute their means, those means will be normally distributed — regardless of the shape of the original data. This remarkable result is the Central Limit Theorem (CLT).

The Central Limit Theorem in action. The original PM2.5 data are right-skewed (top left), but the distribution of 1000 sample means (n = 30 each) is approximately normal (bottom right).

Standard Error of the Mean (SEM)

The standard error measures how precisely a sample mean estimates the true population mean. It shrinks as sample size grows — with more data, your estimate becomes more reliable.

\[SE = \frac{s}{\sqrt{n}}\]

# Course time-use data (hours spent outside class per week)
# R1 = Section 1, R2 = Section 2
R1 <- c(4.5, 2, 10, 7, 10, 10, 6.5, 4, 4, 9, 13, 5, 4.5, 6, 8, 7.5, 4)
R2 <- c(10, 6, 8, 8, 8, 8, 4.5, 8, 3, 10, 10, 9.5, 11, 6, 9)

se <- function(x) sd(x, na.rm = TRUE) / sqrt(sum(!is.na(x)))

knitr::kable(
  data.frame(
    Section = c("Section 1", "Section 2"),
    n       = c(sum(!is.na(R1)), sum(!is.na(R2))),
    Mean    = round(c(mean(R1, na.rm = TRUE), mean(R2, na.rm = TRUE)), 1),
    SD      = round(c(sd(R1, na.rm = TRUE),   sd(R2, na.rm = TRUE)), 1),
    SE      = round(c(se(R1), se(R2)), 2)
  ),
  col.names = c("Section", "n", "Mean (hrs/wk)", "SD", "SE"),
  caption   = "Descriptive statistics for course time-use data by section."
)

Error bars: SD vs. SE — which to use?

This is one of the most misunderstood choices in data visualization.

SD vs SE error bars for the same data (course time-use, two sections). SD bars describe data spread; SE bars describe uncertainty about the mean. The difference grows with sample size.

Formula

\[CI_{95} = \bar{x} \pm t^*_{n-1} \times SE\]

where \(t^*_{n-1}\) is the critical value from the t-distribution with \(n-1\) degrees of freedom. For large samples (\(n > 30\)), \(t^* \approx 1.96\). For smaller samples, \(t^*\) is larger, making the CI wider to reflect greater uncertainty.

ci <- function(x, level = 0.95) {
  x   <- x[!is.na(x)]
  n   <- length(x)
  m   <- mean(x)
  s   <- sd(x)
  se  <- s / sqrt(n)
  t_star <- qt((1 + level) / 2, df = n - 1)
  c(mean  = m,
    lower = m - t_star * se,
    upper = m + t_star * se,
    t_star = t_star)
}

ci_R1 <- ci(R1)
ci_R2 <- ci(R2)

knitr::kable(
  data.frame(
    Section = c("Section 1", "Section 2"),
    Mean    = round(c(ci_R1["mean"],  ci_R2["mean"]),  1),
    CI_lo   = round(c(ci_R1["lower"], ci_R2["lower"]), 1),
    CI_hi   = round(c(ci_R1["upper"], ci_R2["upper"]), 1),
    t_star  = round(c(ci_R1["t_star"],ci_R2["t_star"]),2)
  ),
  col.names = c("Section", "Mean", "95% CI lower", "95% CI upper", "t*"),
  caption   = "95% confidence intervals for mean hours per week outside class."
)

Visualizing CIs vs. SD Bars

95% confidence intervals (left) vs. ±1 SD bars (right) for the same course time-use data. CIs are the preferred choice for inference: non-overlapping CIs suggest a statistically meaningful difference between groups.

What’s next?

With these descriptive tools in hand, the next modules build directly on them:

• t-test — are the two section means significantly different?
• ANOVA — extend the comparison to three or more groups
• Regression — describe the relationship between two continuous variables

ozone <- c(42, 38, 51, 44, 89, 46, 40)
mean(ozone); median(ozone)

ozone_clean <- ozone[ozone != 89]
mean(ozone_clean); median(ozone_clean)

n  <- length(ozone)
sd_oz <- sd(ozone)
se_oz <- sd_oz / sqrt(n)
ci_lower <- mean(ozone) - qt(0.975, df = n-1) * se_oz
ci_upper <- mean(ozone) + qt(0.975, df = n-1) * se_oz
c(lower = ci_lower, mean = mean(ozone), upper = ci_upper)

data(airquality)
hist(airquality$Ozone, breaks = 20,
     main = "Distribution of Ozone (ppb)", xlab = "Ozone (ppb)")
abline(v = mean(airquality$Ozone, na.rm=TRUE),   col="red",  lty=2)
abline(v = median(airquality$Ozone, na.rm=TRUE), col="blue", lty=1)