The five steps

Step 1 — State the hypotheses. Every test has exactly two complementary hypotheses. The null hypothesis (H₀) is the default claim — usually “no effect”, “no difference”, or “the parameter equals a specific value”. The alternative hypothesis (H₁) is what you would conclude if the null is implausible given the data.

Step 2 — Choose a significance level (α). Before collecting data, decide how often you are willing to reject a true null hypothesis by mistake. Convention is α = 0.05 (5%), but this is a social norm, not a law of nature. For high-stakes environmental decisions (e.g., declaring a site contaminated), you might use α = 0.01.

Step 3 — Compute the test statistic. Reduce the data to a single signal-to-noise ratio (t, F, z, W…) that measures how far the observed data are from what the null predicts.

Step 4 — Find the p-value. Ask: if H₀ were true, how often would we observe a test statistic at least as extreme as the one we computed? This probability is the p-value.

Step 5 — Make a decision and report. If p ≤ α: reject H₀ (“the data are inconsistent with no effect”). If p > α: fail to reject H₀ (“the data do not provide sufficient evidence against no effect”). Report the test statistic, degrees of freedom, p-value, and an effect size. Never report only “p < 0.05”.

Type I and Type II errors

No matter what we decide, we can be wrong in two distinct ways:

Step 1 — Visualise against the null value

Always plot the data and mark the hypothesised mean before running any test.

do_df <- data.frame(site = seq_along(do_sites), DO = do_sites)

ggplot(do_df, aes(x = "", y = DO)) +
  geom_hline(yintercept = mu_0, linetype = "dashed",
             colour = "#bf360c", linewidth = 0.9) +
  geom_boxplot(width = 0.35, fill = "#d4e6da", colour = "#1e3a2f",
               outlier.shape = NA, linewidth = 0.7) +
  geom_jitter(width = 0.12, size = 2.2, colour = "#4a7c5e", alpha = 0.7) +
  annotate("text", x = 1.28, y = mu_0 + 0.08,
           label = "EPA standard (7 mg/L)", colour = "#bf360c",
           size = 3.4, hjust = 0, fontface = "italic") +
  stat_summary(fun = mean, geom = "point", shape = 18,
               size = 4, colour = "#c8832a") +
  labs(x = NULL, y = "Dissolved oxygen (mg/L)",
       title = "DO at 22 Stream Sites vs. EPA Standard") +
  theme_minimal(base_size = 12) +
  theme(
    plot.title         = element_text(family = "serif", face = "bold", colour = "#1e3a2f"),
    axis.text.x        = element_blank(),
    panel.grid.major.x = element_blank()
  )

DO measurements at 22 stream sites. Each point is one site; the box shows median and IQR; the dashed red line marks the EPA standard (7 mg/L). Most sites fall below the standard.

Step 2 — Compute by hand, then verify

\[t = \frac{\bar{x} - \mu_0}{s / \sqrt{n}}\]

t_hand <- (mean(do_sites) - mu_0) / (sd(do_sites) / sqrt(n))
p_hand <- 2 * pt(t_hand, df = n - 1)   # two-tailed

cat("t (by hand) =", round(t_hand, 4), "\n")

## t (by hand) = -5.2707

cat("p (by hand) =", round(p_hand, 5), "\n")

## p (by hand) = 3e-05

t.test(do_sites, mu = mu_0)

## 
##  One Sample t-test
## 
## data:  do_sites
## t = -5.2707, df = 21, p-value = 3.172e-05
## alternative hypothesis: true mean is not equal to 7
## 95 percent confidence interval:
##  6.125230 6.620225
## sample estimates:
## mean of x 
##  6.372727

Step 3 — Visualise the t-statistic on the null distribution

The observed t-statistic (red line) on the t(21) null distribution. The amber shaded tails are the 5% rejection region. Our t falls well into the left tail — the sample mean is significantly below the EPA standard.

Step 4 — Report

A one-sample t-test indicated that the regional mean dissolved oxygen (M = 6.37 mg/L, SD = 0.56) was significantly below the EPA coldwater standard of 7 mg/L, t(21) = -5.27, p = 0.

Welch’s t-test (the safe default)

R’s t.test() uses Welch’s correction by default — it does not assume equal variances between groups, making it robust to heteroscedasticity. Only use var.equal = TRUE if you have strong evidence that population variances are equal (Levene’s test, or subject-matter knowledge).

(result_2t <- t.test(C_MgCha ~ treatment, data = carbon))

## 
##  Welch Two Sample t-test
## 
## data:  C_MgCha by treatment
## t = 2.2, df = 28.826, p-value = 0.03598
## alternative hypothesis: true difference in means between group Active replanting and group Passive recovery is not equal to 0
## 95 percent confidence interval:
##  0.359638 9.898948
## sample estimates:
## mean in group Active replanting  mean in group Passive recovery 
##                        41.96818                        36.83889

# Levene's test for equal variances (car package, or manual approach)
var.test(C_MgCha ~ treatment, data = carbon)

## 
##  F test to compare two variances
## 
## data:  C_MgCha by treatment
## F = 0.45695, num df = 21, denom df = 17, p-value = 0.09039
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.1757527 1.1347449
## sample estimates:
## ratio of variances 
##          0.4569535

Wilcoxon signed-rank test (one-sample / paired)

The non-parametric equivalent of the one-sample or paired t-test. Tests whether the median (not mean) equals a hypothesised value, or whether the median difference equals zero.

Environmental context: Heavy metal concentration (mg/kg) in a contaminated floodplain — small sample, right-skewed distribution.

set.seed(19)
# Heavy metal (Pb) — right-skewed, small n
pb_mg   <- c(12, 18, 145, 23, 9, 31, 88, 14, 27, 19, 52, 11)
pb_limit <- 40   # remediation threshold (mg/kg)
n_pb    <- length(pb_mg)

cat("n =", n_pb, "\n")

## n = 12

cat("Mean =", round(mean(pb_mg), 1), "mg/kg\n")

## Mean = 37.4 mg/kg

cat("Median =", median(pb_mg), "mg/kg\n")

## Median = 21 mg/kg

cat("Shapiro-Wilk p =", round(shapiro.test(pb_mg)$p.value, 4),
    "— normality", ifelse(shapiro.test(pb_mg)$p.value < 0.05,
                          "REJECTED", "not rejected"), "\n")

## Shapiro-Wilk p = 0.001 — normality REJECTED

pb_df <- data.frame(pb = pb_mg)

ggplot(pb_df, aes(x = "", y = pb)) +
  geom_hline(yintercept = pb_limit, linetype = "dashed",
             colour = "#bf360c", linewidth = 0.9) +
  geom_boxplot(width = 0.35, fill = "#d4e6da", colour = "#1e3a2f",
               linewidth = 0.7) +
  geom_jitter(width = 0.1, size = 2.5, colour = "#4a7c5e", alpha = 0.75) +
  annotate("text", x = 1.28, y = pb_limit + 5,
           label = "Remediation threshold (40 mg/kg)",
           colour = "#bf360c", size = 3.2, hjust = 0, fontface = "italic") +
  labs(x = NULL, y = "Lead concentration (mg/kg)",
       title = "Pb in Floodplain Soil — Small Skewed Sample") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(family = "serif", face = "bold", colour = "#1e3a2f"),
        axis.text.x = element_blank(),
        panel.grid.major.x = element_blank())

Pb concentration in a small, skewed sample. The single extreme value at 145 mg/kg violates the normality assumption that the t-test requires. The Shapiro-Wilk test confirms non-normality (p < 0.05).

# t-test: sensitive to the outlier, pulls the mean up
cat("── t-test (sensitive to outlier) ──────────────────\n")

## ── t-test (sensitive to outlier) ──────────────────

print(t.test(pb_mg, mu = pb_limit))

## 
##  One Sample t-test
## 
## data:  pb_mg
## t = -0.2205, df = 11, p-value = 0.8295
## alternative hypothesis: true mean is not equal to 40
## 95 percent confidence interval:
##  11.62994 63.20340
## sample estimates:
## mean of x 
##  37.41667

cat("\n── Wilcoxon signed-rank (robust to outlier) ────────\n")

## 
## ── Wilcoxon signed-rank (robust to outlier) ────────

print(wilcox.test(pb_mg, mu = pb_limit, conf.int = TRUE))

## 
##  Wilcoxon signed rank exact test
## 
## data:  pb_mg
## V = 25, p-value = 0.3013
## alternative hypothesis: true location is not equal to 40
## 95 percent confidence interval:
##  15 70
## sample estimates:
## (pseudo)median 
##          24.75

Wilcoxon rank-sum test (two independent groups)

The non-parametric equivalent of the two-sample t-test. Also known as the Mann-Whitney U test. Tests whether one group tends to have larger values than the other.

set.seed(5)
# Species richness on degraded vs. restored plots — count data, right-skewed
restored  <- c(12, 18, 9, 24, 15, 11, 22, 8, 17, 20, 13, 27)
degraded  <- c(4,  7,  3, 11, 6,  2,  9,  5,  8,  4,  6,  3)

cat("Restored  median =", median(restored), "  mean =", round(mean(restored), 1), "\n")

## Restored  median = 16   mean = 16.3

cat("Degraded  median =", median(degraded), "  mean =", round(mean(degraded), 1), "\n")

## Degraded  median = 5.5   mean = 5.7

wilcox.test(restored, degraded, conf.int = TRUE)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  restored and degraded
## W = 139.5, p-value = 0.0001074
## alternative hypothesis: true location shift is not equal to 0
## 95 percent confidence interval:
##   5.999944 15.000018
## sample estimates:
## difference in location 
##               10.00006

rich_df <- data.frame(
  status    = rep(c("Restored", "Degraded"), c(length(restored), length(degraded))),
  richness  = c(restored, degraded)
)
rich_df$status <- factor(rich_df$status, levels = c("Degraded", "Restored"))

ggplot(rich_df, aes(status, richness, fill = status)) +
  geom_boxplot(width = 0.45, outlier.shape = NA, colour = "#1e3a2f",
               alpha = 0.8, linewidth = 0.6) +
  geom_jitter(aes(colour = status), width = 0.16, size = 2, alpha = 0.6) +
  scale_fill_manual(values   = c("Degraded" = "#c8832a", "Restored" = "#4a7c5e")) +
  scale_colour_manual(values = c("Degraded" = "#c8832a", "Restored" = "#4a7c5e")) +
  labs(x = NULL, y = "Species richness",
       title = "Plant Species Richness: Restored vs. Degraded Plots") +
  theme_minimal(base_size = 12) +
  theme(plot.title         = element_text(family = "serif", face = "bold",
                                          colour = "#1e3a2f"),
        legend.position    = "none",
        panel.grid.major.x = element_blank())

Species richness on restored vs. degraded plots. Count data are right-skewed with unequal variances — the Wilcoxon rank-sum test is preferred over the t-test here.

One-proportion z-test

Question: Does the true proportion of exceedances equal a specified value?

Environmental context: In a regulatory sample of 85 wells, 23 exceed the nitrate drinking-water limit (10 mg/L NO₃). Does this fraction differ from the statewide exceedance rate of 15%?

n_wells  <- 85
x_exceed <- 23
p_0      <- 0.15        # statewide rate (null)
p_hat    <- x_exceed / n_wells
cat("Observed proportion:", round(p_hat, 3), "\n")

## Observed proportion: 0.271

# Exact binomial test (preferred when np < 10)
binom.test(x_exceed, n_wells, p = p_0)

## 
##  Exact binomial test
## 
## data:  x_exceed and n_wells
## number of successes = 23, number of trials = 85, p-value = 0.003581
## alternative hypothesis: true probability of success is not equal to 0.15
## 95 percent confidence interval:
##  0.1798891 0.3778595
## sample estimates:
## probability of success 
##              0.2705882

# Large-sample z-test (valid when np ≥ 10 and n(1-p) ≥ 10)
prop.test(x_exceed, n_wells, p = p_0, correct = FALSE)

## 
##  1-sample proportions test without continuity correction
## 
## data:  x_exceed out of n_wells, null probability p_0
## X-squared = 9.6943, df = 1, p-value = 0.001848
## alternative hypothesis: true p is not equal to 0.15
## 95 percent confidence interval:
##  0.1875962 0.3734196
## sample estimates:
##         p 
## 0.2705882

Two-proportion z-test

Question: Do two independent groups have the same proportion?

Environmental context: Compare the fraction of streams exceeding copper criteria upstream vs. downstream of a mining district.

# Upstream: 8 of 40 exceed criterion.  Downstream: 22 of 45.
prop.test(x = c(8, 22), n = c(40, 45), correct = FALSE)

## 
##  2-sample test for equality of proportions without continuity correction
## 
## data:  c(8, 22) out of c(40, 45)
## X-squared = 7.7387, df = 1, p-value = 0.005405
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.48045281 -0.09732497
## sample estimates:
##    prop 1    prop 2 
## 0.2000000 0.4888889

prop_df <- data.frame(
  location = c("Upstream", "Downstream"),
  x = c(8, 22), n = c(40, 45)
) |>
  mutate(
    p     = x / n,
    lower = mapply(function(xi, ni) binom.test(xi, ni)$conf.int[1], x, n),
    upper = mapply(function(xi, ni) binom.test(xi, ni)$conf.int[2], x, n)
  )
prop_df$location <- factor(prop_df$location, levels = c("Upstream", "Downstream"))

ggplot(prop_df, aes(location, p, fill = location)) +
  geom_col(width = 0.45, colour = "#1e3a2f", alpha = 0.85) +
  geom_errorbar(aes(ymin = lower, ymax = upper),
                width = 0.12, linewidth = 0.9, colour = "#1e3a2f") +
  scale_fill_manual(values = c("Upstream" = "#4a7c5e", "Downstream" = "#7a3a2f")) +
  scale_y_continuous(labels = scales::percent_format(), limits = c(0, 0.75)) +
  labs(x = NULL, y = "Proportion exceeding Cu criterion",
       title = "Stream Copper Exceedances Near Mining District") +
  theme_minimal(base_size = 12) +
  theme(plot.title         = element_text(family = "serif", face = "bold",
                                          colour = "#1e3a2f"),
        legend.position    = "none",
        panel.grid.major.x = element_blank())

Proportion of streams exceeding the copper criterion, with 95% Wilson confidence intervals. The upstream rate (~20%) and downstream rate (~49%) are clearly separated — consistent with the significant p-value.

Cohen’s d (for t-tests)

Cohen’s d expresses the difference in means as a multiple of the pooled standard deviation. It is interpretable across studies and variables.

\[d = \frac{\bar{x}_1 - \bar{x}_2}{s_{\text{pooled}}}\]

Cohen’s benchmarks: |d| = 0.2 small · 0.5 medium · 0.8 large.

cohens_d <- function(x, y) {
  nx <- length(x); ny <- length(y)
  sp <- sqrt(((nx - 1)*var(x) + (ny - 1)*var(y)) / (nx + ny - 2))
  (mean(x) - mean(y)) / sp
}

d_carbon <- cohens_d(active, passive)
d_rich   <- cohens_d(restored, degraded)

cat("Carbon (active vs. passive)  d =", round(d_carbon, 3), "\n")

## Carbon (active vs. passive)  d = 0.727

cat("Richness (restored vs. degraded) d =", round(d_rich, 3), "\n")

## Richness (restored vs. degraded) d = 2.271

Visualising Cohen’s d. Each panel shows two normal distributions with d = 0.2, 0.5, 0.8, and 1.4. The overlap shrinks as d grows — d = 0.2 means the distributions are nearly identical; d = 1.4 means a typical value in one group exceeds ~92% of values in the other.

r (for Wilcoxon tests)

For the Wilcoxon rank-sum, the effect size is: \[r = \frac{Z}{\sqrt{N}}\]

where Z is the normal approximation of the W statistic and N is the total sample size. Benchmarks: r = 0.1 small · 0.3 medium · 0.5 large.

w_result <- wilcox.test(restored, degraded, exact = FALSE)
Z_approx <- qnorm(w_result$p.value / 2)
r_effect <- abs(Z_approx) / sqrt(length(restored) + length(degraded))
cat("W =", w_result$statistic, "  Z ≈", round(Z_approx, 3),
    "  r =", round(r_effect, 3), "\n")

## W = 139.5   Z ≈ -3.873   r = 0.791

h (for proportions tests)

Cohen’s h standardises the difference between two proportions via the arcsine transformation, which stabilises variance.

\[h = 2 \arcsin\!\left(\sqrt{p_1}\right) - 2 \arcsin\!\left(\sqrt{p_2}\right)\]

cohens_h <- function(p1, p2) 2*asin(sqrt(p1)) - 2*asin(sqrt(p2))
h_copper <- cohens_h(22/45, 8/40)
cat("Cohen's h =", round(h_copper, 3),
    " (|h| benchmarks: 0.2 small, 0.5 medium, 0.8 large)\n")

## Cohen's h = 0.621  (|h| benchmarks: 0.2 small, 0.5 medium, 0.8 large)

P-values are uniform under the null

When H₀ is true, p-values are uniformly distributed between 0 and 1. There is nothing special about 0.05 except our convention to use it as a threshold.

set.seed(321)
n_trials <- 10000
p_null   <- replicate(n_trials, {
  A <- rnorm(20, 0, 1)
  B <- rnorm(20, 0, 1)     # both from same distribution — null TRUE
  t.test(A, B)$p.value
})
cat("Proportion p < 0.05 when H0 true:", mean(p_null < 0.05), "\n")

## Proportion p < 0.05 when H0 true: 0.0514

cat("(Expected: exactly 0.05 — these are false positives)\n")

## (Expected: exactly 0.05 — these are false positives)

Distribution of p-values from 10,000 t-tests when H₀ is exactly true. The distribution is flat (uniform) — every p-value between 0 and 1 is equally likely. The amber bar at the left (p < 0.05) contains exactly 5% of all tests — these are Type I errors. There is nothing special about values just below 0.05 vs. just above it.

How sample size and effect size jointly determine p

P-values from two-sample t-tests across varying sample sizes (x-axis) and true effect sizes (line colour). A tiny effect (d = 0.1) stays non-significant even with n = 200 per group. A large effect (d = 1.0) reaches significance with n ≈ 10. The p-value is a function of both effect size and sample size — it is not a measure of effect importance.

What’s next?

• Four Frameworks — ANOVA, regression, chi-square, logistic
• Signal Over Noise — permutation tests and the t/F null distributions
• Describing Data — mean, SD, SE, confidence intervals

above <- c(142, 138, 155, 149, 141)
below <- c(198, 210, 225, 189, 215, 202)
shapiro.test(above); shapiro.test(below)
result <- t.test(above, below)
result
# Cohen's d
pooled_sd <- sqrt(((length(above)-1)*var(above) + (length(below)-1)*var(below)) /
                   (length(above) + length(below) - 2))
d <- (mean(below) - mean(above)) / pooled_sd
cat("Cohen's d =", round(d, 2))

set.seed(7)
ctrl    <- c(2.1, 2.4, 2.0, 1.9, 2.3)
compost <- c(3.1, 3.4, 2.9, 3.2, 15.0)  # outlier at 15
t.test(compost, ctrl)
wilcox.test(compost, ctrl)