Standard errors and confidence intervals in R

statistics

uncertainty

ecology tutorial

What a standard error and a 95% confidence interval mean in R, how they differ from the standard deviation, and how to check interval coverage by simulation.

Author

Tidy Ecology

Published

2026-05-06

Modified

2026-09-27

Step 10 of the Tidy Ecology course, Part two: a first result, and what it is worth.

Updated 27 September 2026: a new section, How skewed is too skewed for n = 30, measures the coverage of the t interval at 30 and 100 plots for skewed counts and biomass beside Cochran’s skewness screen for the sample size.

Almost every ecological result is a number plus an admission that the number could have come out differently. You sampled 30 plots, not every plot; another 30 would have given another mean. The standard error and the confidence interval are the two tools for putting that uncertainty on the page. They are also two of the most misread quantities in the field, so this post builds them from a known population, checks them by simulation, and shows where the textbook formula quietly breaks.

We will use a population we control, so we always know the right answer. Picture per-plot biomass in grams, right-skewed the way biomass usually is, drawn from a gamma distribution with a true mean of 20 and a true standard deviation of 10.

library(ggplot2)

paper <- "#f5f4ee"; ink <- "#16241d"; forest <- "#275139"
gold  <- "#cda23f"; miss_red <- "#b5534e"; faint <- "#5d6b61"; line <- "#dad9ca"

theme_te <- theme_minimal(base_size = 12) +
  theme(panel.grid.minor = element_blank(),
        panel.grid.major = element_line(colour = line, linewidth = 0.3),
        plot.background = element_rect(fill = paper, colour = NA),
        panel.background = element_rect(fill = paper, colour = NA),
        axis.title = element_text(colour = ink),
        axis.text = element_text(colour = faint),
        plot.title = element_text(colour = ink, face = "bold"))

shape <- 4; rate <- 0.2
mu_true <- shape / rate            # 20
sd_true <- sqrt(shape) / rate      # 10

The standard error is not the standard deviation

Take one sample of 30 plots and measure the spread two ways.

set.seed(101)
n <- 30
x <- rgamma(n, shape, rate)

m  <- mean(x)
s  <- sd(x)            # standard deviation: spread of the data
se <- s / sqrt(n)      # standard error: spread of the mean

round(c(mean = m, sd = s, se = se), 2)

 mean    sd    se 
19.80  8.82  1.61

The sample mean is 19.8 grams. The standard deviation is 8.82, and the standard error is 1.61, about 5.5 times smaller. That gap is the whole point. The standard deviation answers “how variable are individual plots?” and stays roughly 10 no matter how many plots you measure. The standard error answers “how variable is my estimate of the mean?” and shrinks as the sample grows, because averaging cancels noise. Writing one where you mean the other is the most common error-bar mistake in ecology: standard deviation bars describe the organisms, standard error bars describe your confidence in the average.

The shrinking is exactly $\text{SE} = \sigma / \sqrt{n}$. Across a grid of sample sizes, the average sample standard deviation sits near the true 10 while the average standard error falls off with the square root of n.

set.seed(202)
ns <- c(5, 10, 20, 40, 80, 160, 320)
grid <- do.call(rbind, lapply(ns, function(nn) {
  reps <- 1500
  sds <- numeric(reps); ses <- numeric(reps)
  for (i in seq_len(reps)) {
    xi <- rgamma(nn, shape, rate)
    sds[i] <- sd(xi); ses[i] <- sd(xi) / sqrt(nn)
  }
  data.frame(n = nn, SD = mean(sds), SE = mean(ses))
}))

long <- rbind(
  data.frame(n = grid$n, value = grid$SD, quantity = "standard deviation"),
  data.frame(n = grid$n, value = grid$SE, quantity = "standard error")
)

ggplot(long, aes(n, value, colour = quantity)) +
  geom_line(linewidth = 1) + geom_point(size = 2) +
  scale_x_log10(breaks = ns) +
  scale_colour_manual(values = c("standard deviation" = forest,
                                 "standard error" = gold)) +
  labs(x = "sample size (n, log scale)", y = "average value (g)",
       colour = NULL,
       title = "Spread of the data versus spread of the estimate") +
  theme_te + theme(legend.position = "top")

A line chart on a log x axis. A nearly flat upper line near ten is the mean standard deviation; a lower line falling from about four to half a unit is the mean standard error. — Figure 1: Standard deviation stays near the true value of 10 as sample size grows, while the standard error shrinks with the square root of n. Spread of the data is not spread of the estimate.

Building a 95% confidence interval

A confidence interval turns the standard error into a range. For a mean it is the estimate plus or minus a multiplier from the t distribution times the standard error, where the multiplier accounts for the extra uncertainty of estimating the spread from the same small sample.

tcrit <- qt(0.975, df = n - 1)        # multiplier for 95%, df = n - 1
ci <- m + c(-1, 1) * tcrit * se
round(c(t_multiplier = tcrit, lower = ci[1], upper = ci[2]), 2)

t_multiplier        lower        upper 
        2.05        16.50        23.09

# the same interval, the easy way
as.numeric(t.test(x)$conf.int)

[1] 16.50317 23.09237

With 29 degrees of freedom the multiplier is 2.05, slightly above the 1.96 you would use if the spread were known exactly. The interval runs from 16.5 to 23.09 grams, and t.test() returns the identical bounds, so the formula and the built-in agree.

Now the interpretation, where care is needed. A 95% confidence interval does not mean there is a 95% probability that the true mean lies inside this particular interval. The true mean is a fixed number; it is either in or out. The 95% is a property of the procedure: if you repeated the whole exercise many times, about 95% of the intervals you build would contain the truth. The confidence is in the method over the long run, not in any single interval. That distinction sounds pedantic until you check it.

Does the procedure deliver? Check coverage by simulation

Because we know the true mean is 20, we can build the interval thousands of times and simply count how often it catches the truth. A well-behaved 95% procedure should land near 95%.

coverage <- function(nn, nsim = 4000, seed) {
  set.seed(seed)
  hit <- 0L
  for (i in seq_len(nsim)) {
    xi <- rgamma(nn, shape, rate)
    lo_hi <- mean(xi) + c(-1, 1) * qt(0.975, nn - 1) * sd(xi) / sqrt(nn)
    if (lo_hi[1] <= mu_true && mu_true <= lo_hi[2]) hit <- hit + 1L
  }
  100 * hit / nsim
}

c(n30 = coverage(30, seed = 303), n5 = coverage(5, seed = 404))

  n30    n5 
94.55 93.50

At n = 30 the empirical coverage is 94.6%, close to the promised 95%. At n = 5 it slips to 93.5%. The t interval assumes the sampling distribution of the mean is roughly normal, which the central limit theorem delivers for moderate samples even from skewed data, but with only five skewed observations the approximation is not quite there and the interval is a touch too narrow. The lesson is not that the method is broken; it is that “95%” is earned by sample size and shape, and is worth verifying rather than assuming.

The figure below makes the long-run idea concrete: 100 independent samples, 100 intervals, each a horizontal segment. Five of them miss the true mean, close to the five in a hundred that a 95% procedure misses on average; another 100 samples would give a different count.

set.seed(20)
K <- 100
ci_df <- data.frame(id = seq_len(K), m = NA_real_, lo = NA_real_, hi = NA_real_)
for (i in seq_len(K)) {
  xi <- rgamma(n, shape, rate)
  se_i <- sd(xi) / sqrt(n)
  ci_df$m[i]  <- mean(xi)
  ci_df$lo[i] <- mean(xi) - qt(0.975, n - 1) * se_i
  ci_df$hi[i] <- mean(xi) + qt(0.975, n - 1) * se_i
}
ci_df$hit <- ci_df$lo <= mu_true & mu_true <= ci_df$hi
ci_df <- ci_df[order(ci_df$m), ]
ci_df$rank <- seq_len(K)

ggplot(ci_df, aes(y = rank)) +
  geom_vline(xintercept = mu_true, colour = ink, linewidth = 0.5) +
  geom_segment(aes(x = lo, xend = hi, yend = rank, colour = hit), linewidth = 0.5) +
  geom_point(aes(x = m, colour = hit), size = 0.7) +
  scale_colour_manual(values = c("TRUE" = forest, "FALSE" = miss_red),
                      labels = c("TRUE" = "contains 20", "FALSE" = "misses 20")) +
  labs(x = "biomass (g)", y = "sample (sorted by mean)", colour = NULL,
       title = "Most intervals catch the truth; a predictable few miss") +
  theme_te +
  theme(legend.position = "top",
        axis.text.y = element_blank(), panel.grid.major.y = element_blank())

One hundred horizontal interval segments ordered by their sample mean, forming a diagonal band, a vertical reference line at twenty, five segments highlighted at the top and bottom of the band because they fail to cross the line. — Figure 2: One hundred 95% confidence intervals from one hundred independent samples of n = 30. Each segment is one interval; the vertical line is the true mean. Five intervals miss it in this particular run, close to the 5% the procedure targets on average.

How skewed is too skewed for n = 30

The gamma above has a skewness of exactly 1 (a gamma’s skewness is 2 over the square root of its shape, and the shape is 4), and that is the population on which the check found n = 30 enough. Counts of seedlings, larvae or nests per plot are often far more skewed, because most plots hold a few and one or two hold many. Cochran (1977) gives a screen for when a normal-theory interval around a mean, such as this t interval, is adequate: n > 25 g1^2, where g1 is the skewness of the population. For this gamma the screen asks for more than 25 plots, and the familiar “at least 30” is the same screen for a skewness of 1.10. Boos and Hughes-Oliver (2000) review the role of skewness in how large n must be for t intervals to reach their stated coverage. The chunk runs the coverage check at 30 and 100 plots on six populations: this gamma, a lognormal with the same mean, Poisson counts with mean 2, and negative binomial counts (mean 5 for k = 1 and 0.3, mean 2 for k = 0.1) whose aggregation parameter k (size in rnbinom(), variance mu + mu^2 / k) runs from 1 down to 0.1. Each skewness is a closed form, and the chunk stops if one disagrees with the moments of its distribution.

sk_pop <- data.frame(
  id  = c("Poisson", "gamma, shape 4", "NB k = 1", "NB k = 0.3", "lognormal", "NB k = 0.1"),
  fam = c("pois", "gamma", "nb", "nb", "lnorm", "nb"),
  mu  = c(2, mu_true, 5, 5, mu_true, 2), k = c(NA, NA, 1, 0.3, NA, 0.1))
sk_pop$g1 <- with(sk_pop, ifelse(fam == "pois", 1 / sqrt(mu),     # closed-form skewness
  ifelse(fam == "gamma", 2 / sqrt(shape), ifelse(fam == "nb",
  (1 + 2 * mu / k) / sqrt(mu * (1 + mu / k)), (exp(1) + 2) * sqrt(exp(1) - 1)))))
sk_pop$rule_n <- floor(25 * sk_pop$g1^2) + 1                      # smallest n > 25 g1^2
sk_dens <- function(i, v) with(sk_pop[i, ], switch(fam, pois = dpois(v, mu),
  gamma = dgamma(v, shape, rate), nb = dnbinom(v, mu = mu, size = k),
  lnorm = dlnorm(v, log(mu) - 0.5, 1)))
sk_draw <- function(i, n_draw) with(sk_pop[i, ], switch(fam, pois = rpois(n_draw, mu),
  gamma = rgamma(n_draw, shape, rate), nb = rnbinom(n_draw, mu = mu, size = k),
  lnorm = rlnorm(n_draw, log(mu) - 0.5, 1)))
sk_mom <- function(i, f) if (sk_pop$fam[i] %in% c("pois", "nb")) {
  sum(f(0:20000) * sk_dens(i, 0:20000)) } else {
  integrate(function(v) f(v) * sk_dens(i, v), 0, Inf, rel.tol = 1e-10)$value }
for (sk_i in seq_len(nrow(sk_pop))) local({                   # closed forms vs moments
  m1 <- sk_mom(sk_i, identity); m2 <- sk_mom(sk_i, function(v) (v - m1)^2)
  m3 <- sk_mom(sk_i, function(v) (v - m1)^3); pp <- sk_pop[sk_i, ]
  stopifnot(abs(m1 / pp$mu - 1) < 1e-6, abs(m3 / m2^1.5 / pp$g1 - 1) < 1e-4,
            is.na(pp$k) || abs(m2 / (pp$mu + pp$mu^2 / pp$k) - 1) < 1e-6)
})
sk_cover <- function(i, nn, reps = 50000) {
  xx <- matrix(sk_draw(i, nn * reps), reps, nn); mu_i <- sk_pop$mu[i]
  mm <- rowMeans(xx); dev <- xx - mm; ss <- sqrt(rowSums(dev^2) / (nn - 1))
  half <- qt(0.975, nn - 1) * ss / sqrt(nn)
  g1_hat <- rowMeans(dev^3) / rowMeans(dev^2)^1.5             # NaN for an all-zero sample
  c(cover = 100 * mean(abs(mm - mu_i) <= half), below = 100 * mean(mm + half < mu_i),
    above = 100 * mean(mm - half > mu_i), g1_med = median(g1_hat, na.rm = TRUE),
    passes = 100 * mean(nn > 25 * g1_hat^2, na.rm = TRUE), r_mean_sd = cor(mm, ss))
}
set.seed(27092026); sk_30 <- t(sapply(seq_len(nrow(sk_pop)), sk_cover, nn = n))
sk_100 <- t(sapply(seq_len(nrow(sk_pop)), sk_cover, nn = 100))
cbind(sk_pop[, c("id", "mu", "rule_n")], g1 = round(sk_pop$g1, 2), round(sk_30[, 1:3], 1),
      cover_100 = round(sk_100[, 1], 1), g1_med = round(sk_30[, 4], 2), passes = round(sk_30[, 5]))

              id mu rule_n   g1 cover below above cover_100 g1_med passes
1        Poisson  2     13 0.71  94.7   3.7   1.7      95.0   0.53     89
2 gamma, shape 4 20     26 1.00  94.4   4.3   1.3      94.9   0.73     77
3       NB k = 1  5    101 2.01  92.8   6.5   0.8      94.2   1.39     28
4     NB k = 0.3  5    334 3.65  88.4  11.3   0.3      92.5   2.25      2
5      lognormal 20    957 6.18  88.4  11.3   0.2      91.8   2.01      9
6     NB k = 0.1  2   1001 6.33  80.1  19.9   0.1      88.7   3.21      0

At 30 plots the Poisson counts and the gamma pass the screen and cover 94.7% and 94.4%, a few tenths of a point short of 95% (each coverage here has a Monte Carlo standard error of about 0.1 points near 95% and under 0.2 for the lowest). The other four fail the screen at 30 and undercover: 92.8% for negative binomial counts with k = 1, 88.4% with k = 0.3 and 80.1% with k = 0.1. The misses fall almost all on one side. With k = 0.3, 11.3% of the intervals lie wholly below the true mean and 0.3% wholly above it: a sample that happens to miss the few large counts has a low mean and a low standard deviation at once (their correlation across these samples is 0.85), so its interval is both too low and too narrow.

sk_lab <- sprintf("%s (screen: n >= %d)", sk_pop$id, sk_pop$rule_n)
sk_plot <- data.frame(pop = factor(sk_lab, rev(sk_lab)), cover = c(sk_30[, 1], sk_100[, 1]),
                      plots = factor(rep(c("n = 30", "n = 100"), each = 6), c("n = 30", "n = 100")))
ggplot(sk_plot, aes(cover, pop)) +
  geom_vline(xintercept = 95, colour = ink, linewidth = 0.5) +
  geom_line(aes(group = pop), colour = line, linewidth = 1.4) +
  geom_point(aes(colour = plots), size = 2.8) + scale_x_continuous(limits = c(NA, 96)) +
  scale_colour_manual(values = c("n = 30" = miss_red, "n = 100" = forest)) +
  labs(x = "coverage of the nominal 95% t interval (%)", y = NULL, colour = NULL,
       title = "More skew, more misses at n = 30") +
  theme_te + theme(legend.position = "top", plot.title.position = "plot",
                   panel.grid.major.y = element_blank())

A dot plot with six rows, one per population, each labelled with the sample size the screen requires, from Poisson at 13 at the top to negative binomial with k of 0.1 at 1001 at the bottom. Each row has a red point for 30 plots and a green point for 100 plots joined by a grey bar, and a vertical line marks 95 per cent. The red points of the top two rows sit just left of the line, lower rows sit further left with the bottom red point at 80, and in every row the green point is closer to the line than the red one. — Figure 3: Coverage of the nominal 95% t interval for a mean at 30 and 100 plots, six populations ordered by skewness, with the smallest sample size Cochran’s screen allows in brackets. At 30 plots the two populations the screen allows come close to 95% and the four it does not fall short; every population covers better at 100.

Coverage keeps improving with n: at 100 plots the k = 0.3 counts cover 92.5% and the k = 0.1 counts 88.7%. Nothing breaks at a particular sample size; the screen says how many plots it takes before the shortfall is small, and 100 is still below the 334 and 1001 these two ask for. It is a screen, not a guarantee. The lognormal has almost the skewness of the k = 0.1 counts (6.18 against 6.33) yet covers 88.4% at 30 plots where they cover 80.1%, close to the far less skewed k = 0.3 counts (88.4%), so one moment does not fix the coverage.

The screen needs the skewness of the population, and 30 plots understate it: the median sample skewness is 2.25 for the k = 0.3 counts (population 3.65) and 2.01 for the lognormal (6.18). Applied to each sample’s own skewness (the moment estimator m3 / m2^1.5 in the chunk), the screen passes 28% of the samples of k = 1 counts, whose population fails it, and only 77% of the gamma samples, whose population passes it. For aggregated counts the skewness is set largely by k (the closed form is close to 2 / sqrt(k) when the mean is well above k), and k from a few dozen plots is poorly determined in its own right (Negative binomial k from small host samples). A t interval on log(count + 1) is no way round this for the mean count: it is an interval for the mean of the logged counts, a different quantity (Poisson and negative binomial GLMs in R covers the logging reflex). When counts fail the screen, the options include more plots or a method that models the count distribution, and any of them deserves the same coverage check.

Proportions need a better interval

Counts and proportions are everywhere in ecology: occupied sites, infected hosts, germinated seeds. The obvious interval, $\hat{p} \pm 1.96\sqrt{\hat{p}(1-\hat{p})/n}$, is the Wald interval, and it misbehaves badly near 0 and 1. Suppose 3 of 25 sites are occupied.

x_occ <- 3; n_site <- 25
phat <- x_occ / n_site

wald   <- phat + c(-1, 1) * qnorm(0.975) * sqrt(phat * (1 - phat) / n_site)
wilson <- as.numeric(prop.test(x_occ, n_site, correct = FALSE)$conf.int)

round(c(p_hat = phat, wald_lo = wald[1], wald_hi = wald[2],
        wilson_lo = wilson[1], wilson_hi = wilson[2]), 3)

    p_hat   wald_lo   wald_hi wilson_lo wilson_hi 
    0.120    -0.007     0.247     0.042     0.300

The estimate is 0.12. The Wald interval runs down to -0.007, a negative proportion, which is impossible. The Wilson interval from prop.test() stays inside [0, 1] at 0.042 to 0.300 and keeps its coverage near nominal even for small counts. For proportions, reach for the Wilson interval (or an exact method) and leave Wald for large samples well away from the boundaries.

What to carry forward

Report the standard error or a confidence interval, not the standard deviation, when the claim is about a mean or a model coefficient: readers want the precision of your estimate, not the variability of individual organisms. Read a confidence interval as a statement about the procedure, not a probability attached to the one interval you happened to compute. Treat the nominal level as a target to verify: small samples and skew can erode it, and a few lines of simulation will tell you by how much. For proportions, use an interval that respects the [0, 1] boundary. None of this requires a package beyond base R, and all of it travels directly from these synthetic plots to real field data.

References

Cumming G 2014 Psychological Science 25(1):7-29 (10.1177/0956797613504966)

Brown LD, Cai TT, DasGupta A 2001 Statistical Science 16(2):101-133 (10.1214/ss/1009213286)

Crawley MJ 2013 The R Book, 2nd edn. Wiley. ISBN 9780470973929

Cochran WG 1977 Sampling Techniques, 3rd edn. Wiley. ISBN 978-0-471-16240-7

Boos DD, Hughes-Oliver JM 2000 The American Statistician 54(2):121-128 (10.1080/00031305.2000.10474524)

--- title: "Standard errors and confidence intervals in R" description: "What a standard error and a 95% confidence interval mean in R, how they differ from the standard deviation, and how to check interval coverage by simulation." date: "2026-05-06 15:00" date-modified: "2026-09-27" categories: [statistics, uncertainty, ecology tutorial, R] image: thumbnail.png image-alt: "A column of horizontal confidence intervals around a vertical line marking the true mean, most crossing it and a few missing." --- *Updated 27 September 2026: a new section, How skewed is too skewed for n = 30, measures the coverage of the t interval at 30 and 100 plots for skewed counts and biomass beside Cochran's skewness screen for the sample size.* Almost every ecological result is a number plus an admission that the number could have come out differently. You sampled 30 plots, not every plot; another 30 would have given another mean. The standard error and the confidence interval are the two tools for putting that uncertainty on the page. They are also two of the most misread quantities in the field, so this post builds them from a known population, checks them by simulation, and shows where the textbook formula quietly breaks. We will use a population we control, so we always know the right answer. Picture per-plot biomass in grams, right-skewed the way biomass usually is, drawn from a gamma distribution with a true mean of 20 and a true standard deviation of 10. ```{r} #| label: setup #| message: false #| warning: false library(ggplot2) paper <- "#f5f4ee"; ink <- "#16241d"; forest <- "#275139" gold <- "#cda23f"; miss_red <- "#b5534e"; faint <- "#5d6b61"; line <- "#dad9ca" theme_te <- theme_minimal(base_size = 12) + theme(panel.grid.minor = element_blank(), panel.grid.major = element_line(colour = line, linewidth = 0.3), plot.background = element_rect(fill = paper, colour = NA), panel.background = element_rect(fill = paper, colour = NA), axis.title = element_text(colour = ink), axis.text = element_text(colour = faint), plot.title = element_text(colour = ink, face = "bold")) shape <- 4; rate <- 0.2 mu_true <- shape / rate # 20 sd_true <- sqrt(shape) / rate # 10 ``` ## The standard error is not the standard deviation Take one sample of 30 plots and measure the spread two ways. ```{r} #| label: one-sample set.seed(101) n <- 30 x <- rgamma(n, shape, rate) m <- mean(x) s <- sd(x) # standard deviation: spread of the data se <- s / sqrt(n) # standard error: spread of the mean round(c(mean = m, sd = s, se = se), 2) ``` The sample mean is `r round(m, 2)` grams. The standard deviation is `r round(s, 2)`, and the standard error is `r round(se, 2)`, about `r round(s/se, 1)` times smaller. That gap is the whole point. The standard deviation answers "how variable are individual plots?" and stays roughly 10 no matter how many plots you measure. The standard error answers "how variable is my estimate of the mean?" and shrinks as the sample grows, because averaging cancels noise. Writing one where you mean the other is the most common error-bar mistake in ecology: standard deviation bars describe the organisms, standard error bars describe your confidence in the average. The shrinking is exactly $\text{SE} = \sigma / \sqrt{n}$. Across a grid of sample sizes, the average sample standard deviation sits near the true 10 while the average standard error falls off with the square root of n. ```{r} #| label: fig-sd-se #| fig-cap: "Standard deviation stays near the true value of 10 as sample size grows, while the standard error shrinks with the square root of n. Spread of the data is not spread of the estimate." #| fig-alt: "A line chart on a log x axis. A nearly flat upper line near ten is the mean standard deviation; a lower line falling from about four to half a unit is the mean standard error." #| fig-width: 7 #| fig-height: 4.2 set.seed(202) ns <- c(5, 10, 20, 40, 80, 160, 320) grid <- do.call(rbind, lapply(ns, function(nn) { reps <- 1500 sds <- numeric(reps); ses <- numeric(reps) for (i in seq_len(reps)) { xi <- rgamma(nn, shape, rate) sds[i] <- sd(xi); ses[i] <- sd(xi) / sqrt(nn) } data.frame(n = nn, SD = mean(sds), SE = mean(ses)) })) long <- rbind( data.frame(n = grid$n, value = grid$SD, quantity = "standard deviation"), data.frame(n = grid$n, value = grid$SE, quantity = "standard error") ) ggplot(long, aes(n, value, colour = quantity)) + geom_line(linewidth = 1) + geom_point(size = 2) + scale_x_log10(breaks = ns) + scale_colour_manual(values = c("standard deviation" = forest, "standard error" = gold)) + labs(x = "sample size (n, log scale)", y = "average value (g)", colour = NULL, title = "Spread of the data versus spread of the estimate") + theme_te + theme(legend.position = "top") ``` ## Building a 95% confidence interval A confidence interval turns the standard error into a range. For a mean it is the estimate plus or minus a multiplier from the t distribution times the standard error, where the multiplier accounts for the extra uncertainty of estimating the spread from the same small sample. ```{r} #| label: ci-mean tcrit <- qt(0.975, df = n - 1) # multiplier for 95%, df = n - 1 ci <- m + c(-1, 1) * tcrit * se round(c(t_multiplier = tcrit, lower = ci[1], upper = ci[2]), 2) # the same interval, the easy way as.numeric(t.test(x)$conf.int) ``` With 29 degrees of freedom the multiplier is `r round(tcrit, 2)`, slightly above the 1.96 you would use if the spread were known exactly. The interval runs from `r round(ci[1], 2)` to `r round(ci[2], 2)` grams, and `t.test()` returns the identical bounds, so the formula and the built-in agree. Now the interpretation, where care is needed. A 95% confidence interval does not mean there is a 95% probability that the true mean lies inside this particular interval. The true mean is a fixed number; it is either in or out. The 95% is a property of the procedure: if you repeated the whole exercise many times, about 95% of the intervals you build would contain the truth. The confidence is in the method over the long run, not in any single interval. That distinction sounds pedantic until you check it. ## Does the procedure deliver? Check coverage by simulation Because we know the true mean is 20, we can build the interval thousands of times and simply count how often it catches the truth. A well-behaved 95% procedure should land near 95%. ```{r} #| label: coverage coverage <- function(nn, nsim = 4000, seed) { set.seed(seed) hit <- 0L for (i in seq_len(nsim)) { xi <- rgamma(nn, shape, rate) lo_hi <- mean(xi) + c(-1, 1) * qt(0.975, nn - 1) * sd(xi) / sqrt(nn) if (lo_hi[1] <= mu_true && mu_true <= lo_hi[2]) hit <- hit + 1L } 100 * hit / nsim } c(n30 = coverage(30, seed = 303), n5 = coverage(5, seed = 404)) ``` At n = 30 the empirical coverage is `r round(coverage(30, seed = 303), 1)`%, close to the promised 95%. At n = 5 it slips to `r round(coverage(5, seed = 404), 1)`%. The t interval assumes the sampling distribution of the mean is roughly normal, which the central limit theorem delivers for moderate samples even from skewed data, but with only five skewed observations the approximation is not quite there and the interval is a touch too narrow. The lesson is not that the method is broken; it is that "95%" is earned by sample size and shape, and is worth verifying rather than assuming. The figure below makes the long-run idea concrete: 100 independent samples, 100 intervals, each a horizontal segment. Five of them miss the true mean, close to the five in a hundred that a 95% procedure misses on average; another 100 samples would give a different count. ```{r} #| label: fig-coverage #| fig-cap: "One hundred 95% confidence intervals from one hundred independent samples of n = 30. Each segment is one interval; the vertical line is the true mean. Five intervals miss it in this particular run, close to the 5% the procedure targets on average." #| fig-alt: "One hundred horizontal interval segments ordered by their sample mean, forming a diagonal band, a vertical reference line at twenty, five segments highlighted at the top and bottom of the band because they fail to cross the line." #| fig-width: 7 #| fig-height: 5 set.seed(20) K <- 100 ci_df <- data.frame(id = seq_len(K), m = NA_real_, lo = NA_real_, hi = NA_real_) for (i in seq_len(K)) { xi <- rgamma(n, shape, rate) se_i <- sd(xi) / sqrt(n) ci_df$m[i] <- mean(xi) ci_df$lo[i] <- mean(xi) - qt(0.975, n - 1) * se_i ci_df$hi[i] <- mean(xi) + qt(0.975, n - 1) * se_i } ci_df$hit <- ci_df$lo <= mu_true & mu_true <= ci_df$hi ci_df <- ci_df[order(ci_df$m), ] ci_df$rank <- seq_len(K) ggplot(ci_df, aes(y = rank)) + geom_vline(xintercept = mu_true, colour = ink, linewidth = 0.5) + geom_segment(aes(x = lo, xend = hi, yend = rank, colour = hit), linewidth = 0.5) + geom_point(aes(x = m, colour = hit), size = 0.7) + scale_colour_manual(values = c("TRUE" = forest, "FALSE" = miss_red), labels = c("TRUE" = "contains 20", "FALSE" = "misses 20")) + labs(x = "biomass (g)", y = "sample (sorted by mean)", colour = NULL, title = "Most intervals catch the truth; a predictable few miss") + theme_te + theme(legend.position = "top", axis.text.y = element_blank(), panel.grid.major.y = element_blank()) ``` ## How skewed is too skewed for n = 30 The gamma above has a skewness of exactly 1 (a gamma's skewness is 2 over the square root of its shape, and the shape is 4), and that is the population on which the check found n = 30 enough. Counts of seedlings, larvae or nests per plot are often far more skewed, because most plots hold a few and one or two hold many. Cochran (1977) gives a screen for when a normal-theory interval around a mean, such as this t interval, is adequate: `n > 25 g1^2`, where `g1` is the skewness of the population. For this gamma the screen asks for more than 25 plots, and the familiar "at least 30" is the same screen for a skewness of `r sprintf("%.2f", sqrt(30 / 25))`. Boos and Hughes-Oliver (2000) review the role of skewness in how large n must be for t intervals to reach their stated coverage. The chunk runs the coverage check at 30 and 100 plots on six populations: this gamma, a lognormal with the same mean, Poisson counts with mean 2, and negative binomial counts (mean 5 for k = 1 and 0.3, mean 2 for k = 0.1) whose aggregation parameter k (`size` in `rnbinom()`, variance `mu + mu^2 / k`) runs from 1 down to 0.1. Each skewness is a closed form, and the chunk stops if one disagrees with the moments of its distribution. ```{r} #| label: skew-coverage sk_pop <- data.frame( id = c("Poisson", "gamma, shape 4", "NB k = 1", "NB k = 0.3", "lognormal", "NB k = 0.1"), fam = c("pois", "gamma", "nb", "nb", "lnorm", "nb"), mu = c(2, mu_true, 5, 5, mu_true, 2), k = c(NA, NA, 1, 0.3, NA, 0.1)) sk_pop$g1 <- with(sk_pop, ifelse(fam == "pois", 1 / sqrt(mu), # closed-form skewness ifelse(fam == "gamma", 2 / sqrt(shape), ifelse(fam == "nb", (1 + 2 * mu / k) / sqrt(mu * (1 + mu / k)), (exp(1) + 2) * sqrt(exp(1) - 1))))) sk_pop$rule_n <- floor(25 * sk_pop$g1^2) + 1 # smallest n > 25 g1^2 sk_dens <- function(i, v) with(sk_pop[i, ], switch(fam, pois = dpois(v, mu), gamma = dgamma(v, shape, rate), nb = dnbinom(v, mu = mu, size = k), lnorm = dlnorm(v, log(mu) - 0.5, 1))) sk_draw <- function(i, n_draw) with(sk_pop[i, ], switch(fam, pois = rpois(n_draw, mu), gamma = rgamma(n_draw, shape, rate), nb = rnbinom(n_draw, mu = mu, size = k), lnorm = rlnorm(n_draw, log(mu) - 0.5, 1))) sk_mom <- function(i, f) if (sk_pop$fam[i] %in% c("pois", "nb")) { sum(f(0:20000) * sk_dens(i, 0:20000)) } else { integrate(function(v) f(v) * sk_dens(i, v), 0, Inf, rel.tol = 1e-10)$value } for (sk_i in seq_len(nrow(sk_pop))) local({ # closed forms vs moments m1 <- sk_mom(sk_i, identity); m2 <- sk_mom(sk_i, function(v) (v - m1)^2) m3 <- sk_mom(sk_i, function(v) (v - m1)^3); pp <- sk_pop[sk_i, ] stopifnot(abs(m1 / pp$mu - 1) < 1e-6, abs(m3 / m2^1.5 / pp$g1 - 1) < 1e-4, is.na(pp$k) || abs(m2 / (pp$mu + pp$mu^2 / pp$k) - 1) < 1e-6) }) sk_cover <- function(i, nn, reps = 50000) { xx <- matrix(sk_draw(i, nn * reps), reps, nn); mu_i <- sk_pop$mu[i] mm <- rowMeans(xx); dev <- xx - mm; ss <- sqrt(rowSums(dev^2) / (nn - 1)) half <- qt(0.975, nn - 1) * ss / sqrt(nn) g1_hat <- rowMeans(dev^3) / rowMeans(dev^2)^1.5 # NaN for an all-zero sample c(cover = 100 * mean(abs(mm - mu_i) <= half), below = 100 * mean(mm + half < mu_i), above = 100 * mean(mm - half > mu_i), g1_med = median(g1_hat, na.rm = TRUE), passes = 100 * mean(nn > 25 * g1_hat^2, na.rm = TRUE), r_mean_sd = cor(mm, ss)) } set.seed(27092026); sk_30 <- t(sapply(seq_len(nrow(sk_pop)), sk_cover, nn = n)) sk_100 <- t(sapply(seq_len(nrow(sk_pop)), sk_cover, nn = 100)) cbind(sk_pop[, c("id", "mu", "rule_n")], g1 = round(sk_pop$g1, 2), round(sk_30[, 1:3], 1), cover_100 = round(sk_100[, 1], 1), g1_med = round(sk_30[, 4], 2), passes = round(sk_30[, 5])) ``` At 30 plots the Poisson counts and the gamma pass the screen and cover `r sprintf("%.1f", sk_30[1, "cover"])`% and `r sprintf("%.1f", sk_30[2, "cover"])`%, a few tenths of a point short of 95% (each coverage here has a Monte Carlo standard error of about `r sprintf("%.1f", 100 * sqrt(0.95 * 0.05 / 50000))` points near 95% and under `r sprintf("%.1f", 100 * sqrt(0.8 * 0.2 / 50000))` for the lowest). The other four fail the screen at 30 and undercover: `r sprintf("%.1f", sk_30[3, "cover"])`% for negative binomial counts with k = 1, `r sprintf("%.1f", sk_30[4, "cover"])`% with k = 0.3 and `r sprintf("%.1f", sk_30[6, "cover"])`% with k = 0.1. The misses fall almost all on one side. With k = 0.3, `r sprintf("%.1f", sk_30[4, "below"])`% of the intervals lie wholly below the true mean and `r sprintf("%.1f", sk_30[4, "above"])`% wholly above it: a sample that happens to miss the few large counts has a low mean and a low standard deviation at once (their correlation across these samples is `r sprintf("%.2f", sk_30[4, "r_mean_sd"])`), so its interval is both too low and too narrow. ```{r} #| label: fig-skew-coverage #| fig-cap: "Coverage of the nominal 95% t interval for a mean at 30 and 100 plots, six populations ordered by skewness, with the smallest sample size Cochran's screen allows in brackets. At 30 plots the two populations the screen allows come close to 95% and the four it does not fall short; every population covers better at 100." #| fig-alt: "A dot plot with six rows, one per population, each labelled with the sample size the screen requires, from Poisson at 13 at the top to negative binomial with k of 0.1 at 1001 at the bottom. Each row has a red point for 30 plots and a green point for 100 plots joined by a grey bar, and a vertical line marks 95 per cent. The red points of the top two rows sit just left of the line, lower rows sit further left with the bottom red point at 80, and in every row the green point is closer to the line than the red one." #| fig-width: 7 #| fig-height: 4.2 sk_lab <- sprintf("%s (screen: n >= %d)", sk_pop$id, sk_pop$rule_n) sk_plot <- data.frame(pop = factor(sk_lab, rev(sk_lab)), cover = c(sk_30[, 1], sk_100[, 1]), plots = factor(rep(c("n = 30", "n = 100"), each = 6), c("n = 30", "n = 100"))) ggplot(sk_plot, aes(cover, pop)) + geom_vline(xintercept = 95, colour = ink, linewidth = 0.5) + geom_line(aes(group = pop), colour = line, linewidth = 1.4) + geom_point(aes(colour = plots), size = 2.8) + scale_x_continuous(limits = c(NA, 96)) + scale_colour_manual(values = c("n = 30" = miss_red, "n = 100" = forest)) + labs(x = "coverage of the nominal 95% t interval (%)", y = NULL, colour = NULL, title = "More skew, more misses at n = 30") + theme_te + theme(legend.position = "top", plot.title.position = "plot", panel.grid.major.y = element_blank()) ``` Coverage keeps improving with n: at 100 plots the k = 0.3 counts cover `r sprintf("%.1f", sk_100[4, "cover"])`% and the k = 0.1 counts `r sprintf("%.1f", sk_100[6, "cover"])`%. Nothing breaks at a particular sample size; the screen says how many plots it takes before the shortfall is small, and 100 is still below the `r sk_pop$rule_n[4]` and `r sk_pop$rule_n[6]` these two ask for. It is a screen, not a guarantee. The lognormal has almost the skewness of the k = 0.1 counts (`r sprintf("%.2f", sk_pop$g1[5])` against `r sprintf("%.2f", sk_pop$g1[6])`) yet covers `r sprintf("%.1f", sk_30[5, "cover"])`% at 30 plots where they cover `r sprintf("%.1f", sk_30[6, "cover"])`%, close to the far less skewed k = 0.3 counts (`r sprintf("%.1f", sk_30[4, "cover"])`%), so one moment does not fix the coverage. The screen needs the skewness of the population, and 30 plots understate it: the median sample skewness is `r sprintf("%.2f", sk_30[4, "g1_med"])` for the k = 0.3 counts (population `r sprintf("%.2f", sk_pop$g1[4])`) and `r sprintf("%.2f", sk_30[5, "g1_med"])` for the lognormal (`r sprintf("%.2f", sk_pop$g1[5])`). Applied to each sample's own skewness (the moment estimator `m3 / m2^1.5` in the chunk), the screen passes `r sprintf("%.0f", sk_30[3, "passes"])`% of the samples of k = 1 counts, whose population fails it, and only `r sprintf("%.0f", sk_30[2, "passes"])`% of the gamma samples, whose population passes it. For aggregated counts the skewness is set largely by k (the closed form is close to `2 / sqrt(k)` when the mean is well above k), and k from a few dozen plots is poorly determined in its own right ([Negative binomial k from small host samples](../negative-binomial-k-from-small-samples/)). A t interval on `log(count + 1)` is no way round this for the mean count: it is an interval for the mean of the logged counts, a different quantity ([Poisson and negative binomial GLMs in R](../glm-count-data-abundance/) covers the logging reflex). When counts fail the screen, the options include more plots or a method that models the count distribution, and any of them deserves the same coverage check. ## Proportions need a better interval Counts and proportions are everywhere in ecology: occupied sites, infected hosts, germinated seeds. The obvious interval, $\hat{p} \pm 1.96\sqrt{\hat{p}(1-\hat{p})/n}$, is the Wald interval, and it misbehaves badly near 0 and 1. Suppose 3 of 25 sites are occupied. ```{r} #| label: prop-ci x_occ <- 3; n_site <- 25 phat <- x_occ / n_site wald <- phat + c(-1, 1) * qnorm(0.975) * sqrt(phat * (1 - phat) / n_site) wilson <- as.numeric(prop.test(x_occ, n_site, correct = FALSE)$conf.int) round(c(p_hat = phat, wald_lo = wald[1], wald_hi = wald[2], wilson_lo = wilson[1], wilson_hi = wilson[2]), 3) ``` The estimate is `r round(phat, 2)`. The Wald interval runs down to `r round(wald[1], 3)`, a negative proportion, which is impossible. The Wilson interval from `prop.test()` stays inside [0, 1] at `r sprintf("%.3f", wilson[1])` to `r sprintf("%.3f", wilson[2])` and keeps its coverage near nominal even for small counts. For proportions, reach for the Wilson interval (or an exact method) and leave Wald for large samples well away from the boundaries. ## What to carry forward Report the standard error or a confidence interval, not the standard deviation, when the claim is about a mean or a model coefficient: readers want the precision of your estimate, not the variability of individual organisms. Read a confidence interval as a statement about the procedure, not a probability attached to the one interval you happened to compute. Treat the nominal level as a target to verify: small samples and skew can erode it, and a few lines of simulation will tell you by how much. For proportions, use an interval that respects the [0, 1] boundary. None of this requires a package beyond base R, and all of it travels directly from these synthetic plots to real field data. ## References Cumming G 2014 Psychological Science 25(1):7-29 (10.1177/0956797613504966) Brown LD, Cai TT, DasGupta A 2001 Statistical Science 16(2):101-133 (10.1214/ss/1009213286) Crawley MJ 2013 The R Book, 2nd edn. Wiley. ISBN 9780470973929 Cochran WG 1977 Sampling Techniques, 3rd edn. Wiley. ISBN 978-0-471-16240-7 Boos DD, Hughes-Oliver JM 2000 The American Statistician 54(2):121-128 (10.1080/00031305.2000.10474524) ## Related tutorials - [Bootstrap confidence intervals in R](../bootstrap-confidence-intervals/) - [The delta method: standard errors for derived values](../the-delta-method/) - [What a p-value tells you, and what it does not](../what-a-p-value-tells-you/) - [Maximum likelihood by hand with optim() in R](../maximum-likelihood-with-optim/)