The hypothesis test

Testing whether r differs from zero uses a t-transformation:

\[t = \frac{r\sqrt{n-2}}{\sqrt{1-r^2}}, \quad df = n-2\]

Assumptions

R walkthrough — water quality dataset

We use a simulated dataset of 55 lake monitoring sites, measuring five water chemistry variables collected during summer stratification: surface temperature (°C), dissolved oxygen (mg/L), pH, total phosphorus (μg/L), and chlorophyll-a (μg/L).

set.seed(314)
n <- 55

temp_C  <- round(rnorm(n, 22, 3.2), 1)
do_mgl  <- round(11.2 - 0.29 * temp_C + rnorm(n, 0, 0.7), 2)
pH      <- round(7.4  + 0.04 * temp_C  + rnorm(n, 0, 0.25), 2)
tp_ugl  <- round(exp(2.8 - 0.07 * do_mgl + rnorm(n, 0, 0.45)), 1)
chla    <- round(exp(0.9 + 0.018 * tp_ugl + rnorm(n, 0, 0.5)), 2)

lakes <- data.frame(temp_C, do_mgl, pH, tp_ugl, chla)
head(lakes, 5)

##   temp_C do_mgl   pH tp_ugl chla
## 1   17.9   5.00 8.19   12.6 5.93
## 2   24.3   3.91 7.97    7.4 3.13
## 3   19.3   6.49 8.17   18.3 3.22
## 4   19.7   5.88 8.20   17.1 2.17
## 5   22.4   4.39 8.41   12.7 3.87

par(mar = c(4.5, 4.5, 3, 1.5), las = 1)

plot(lakes$temp_C, lakes$do_mgl,
     pch = 19, cex = 0.9,
     col = adjustcolor("#1e3a2f", 0.65),
     xlab = "Surface temperature (°C)",
     ylab = "Dissolved oxygen (mg/L)",
     main = "Temperature vs. dissolved oxygen — lake summer survey",
     cex.main = 0.91)

fit_td  <- lm(do_mgl ~ temp_C, data = lakes)
xseq    <- seq(min(lakes$temp_C), max(lakes$temp_C), length.out = 200)
ci      <- predict(fit_td, data.frame(temp_C = xseq), interval = "confidence")
polygon(c(xseq, rev(xseq)), c(ci[,"lwr"], rev(ci[,"upr"])),
        col = adjustcolor("#d4e6da", 0.55), border = NA)
lines(xseq, ci[,"fit"], col = "#c8832a", lwd = 2)

r_td <- cor(lakes$temp_C, lakes$do_mgl)
legend("topright",
       legend = c(sprintf("Pearson r = %.3f", r_td), "OLS fit ± 95% CI"),
       lty = c(NA, 1), pch = c(19, NA),
       col = c(adjustcolor("#1e3a2f", .65), "#c8832a"),
       lwd = c(NA, 2), bty = "n", cex = 0.82)

Surface temperature vs. dissolved oxygen at 55 lake monitoring sites. The amber line is the OLS regression line (shown for orientation); the shaded band is a 95% confidence interval around the mean. The negative relationship reflects the physical constraint that warmer water holds less dissolved gas.

cor.test(lakes$temp_C, lakes$do_mgl, method = "pearson")

## 
##  Pearson's product-moment correlation
## 
## data:  lakes$temp_C and lakes$do_mgl
## t = -10.78, df = 53, p-value = 5.705e-15
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.8968369 -0.7221960
## sample estimates:
##        cor 
## -0.8287144

par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
qqnorm(lakes$temp_C, main = "QQ — Temperature",
       col = "#1e3a2f", pch = 19, cex = 0.75, cex.main = 0.88)
qqline(lakes$temp_C, col = "#c8832a", lwd = 1.8)
qqnorm(lakes$do_mgl, main = "QQ — Dissolved oxygen",
       col = "#1e3a2f", pch = 19, cex = 0.75, cex.main = 0.88)
qqline(lakes$do_mgl, col = "#c8832a", lwd = 1.8)

Shapiro-Wilk test results and QQ plots confirm that both variables are approximately normally distributed, satisfying the key assumption for Pearson’s test.

par(mfrow = c(1,1))

shapiro.test(lakes$temp_C)

## 
##  Shapiro-Wilk normality test
## 
## data:  lakes$temp_C
## W = 0.97368, p-value = 0.2678

shapiro.test(lakes$do_mgl)

## 
##  Shapiro-Wilk normality test
## 
## data:  lakes$do_mgl
## W = 0.98431, p-value = 0.688

Confidence interval for r

R reports a 95% CI for ρ using Fisher’s z-transformation, which stabilises the variance of r across its range. The transformation is:

\[z = \tanh^{-1}(r) = \frac{1}{2}\ln\!\left(\frac{1+r}{1-r}\right), \quad \text{SE}(z) = \frac{1}{\sqrt{n-3}}\]

The CI is computed on the z scale and back-transformed to r. This CI is reported automatically by cor.test() in the output above.

z_x  <- scale(lakes$temp_C)   # (x - x̄) / sd(x)
z_y  <- scale(lakes$do_mgl)
r    <- sum(z_x * z_y) / (n - 1)

Why rank?

Ranks are bounded (1 to n), which makes ρ resistant to outliers and insensitive to the scale of measurement. Replacing raw values with their ranks is equivalent to a monotone transformation of the data — it preserves the order but removes information about the magnitude of gaps between values.

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

# Raw
plot(lakes$tp_ugl, lakes$chla,
     pch = 19, cex = 0.85, col = adjustcolor("#1e3a2f", 0.6),
     xlab = "Total phosphorus (μg/L)", ylab = "Chlorophyll-a (μg/L)",
     main = "Raw values", cex.main = 0.9)
abline(lm(chla ~ tp_ugl, data = lakes), col = "#c8832a", lwd = 1.8)

# Ranks
plot(rank(lakes$tp_ugl), rank(lakes$chla),
     pch = 19, cex = 0.85, col = adjustcolor("#2d5a3d", 0.6),
     xlab = "Rank of total phosphorus", ylab = "Rank of chlorophyll-a",
     main = "After rank transformation", cex.main = 0.9)
abline(lm(rank(chla) ~ rank(tp_ugl), data = lakes),
       col = "#c8832a", lwd = 1.8)

Left: raw total phosphorus vs. chlorophyll-a. The relationship is positive but clearly non-linear (the right tail of TP drives a disproportionate increase in chl-a). Right: the same data after rank transformation — the monotone relationship becomes close to linear on the rank scale. Pearson’s r on ranks is Spearman’s ρ.

par(mfrow = c(1,1))

# Pearson on raw (underestimates due to non-linearity)
cor.test(lakes$tp_ugl, lakes$chla, method = "pearson")

## 
##  Pearson's product-moment correlation
## 
## data:  lakes$tp_ugl and lakes$chla
## t = 1.4989, df = 53, p-value = 0.1398
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.06722891  0.44324798
## sample estimates:
##       cor 
## 0.2016651

# Spearman on ranks (captures monotone relationship more faithfully)
cor.test(lakes$tp_ugl, lakes$chla, method = "spearman")

## 
##  Spearman's rank correlation rho
## 
## data:  lakes$tp_ugl and lakes$chla
## S = 20204, p-value = 0.04525
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.2711446

The contrast between the two outputs illustrates the core trade-off: Pearson’s r is attenuated by the non-linear scaling between TP and chl-a, while Spearman’s ρ captures the full monotone trend.

Base R pairs()

pairs(lakes,
      labels  = c("Temp\n(°C)", "DO\n(mg/L)", "pH", "TP\n(μg/L)", "Chl-a\n(μg/L)"),
      pch     = 19,
      cex     = 0.55,
      col     = adjustcolor("#1e3a2f", 0.55),
      gap     = 0.4,
      cex.labels = 1.0,
      font.labels = 2)

Base R pairs() matrix for the five lake chemistry variables. Each off-diagonal cell shows the scatterplot for that variable pair. Scan column-by-column or row-by-row to identify which pairs have strong linear or monotone relationships, which are non-linear, and which are independent.

Enhanced pairs panel with r values and smooths

Adding the correlation coefficient directly to each panel and a loess smooth makes the matrix far more informative for exploratory work.

# ── Panel functions ────────────────────────────────────────────────────────
panel_scatter <- function(x, y, ...) {
  points(x, y, pch = 19, cex = 0.5,
         col = adjustcolor("#1e3a2f", 0.5))
  ok <- is.finite(x) & is.finite(y)
  if (sum(ok) > 4)
    lines(lowess(x[ok], y[ok]), col = "#c8832a", lwd = 1.8)
}

panel_r <- function(x, y, ...) {
  usr <- par("usr"); on.exit(par(usr))
  par(usr = c(0, 1, 0, 1))
  r   <- cor(x, y, use = "pairwise.complete.obs")
  cex_r <- 0.8 + 1.4 * abs(r)
  col_r <- if (r > 0.5)  "#1e3a2f" else
            if (r < -0.5) "#7a1e1e" else "#7a7a72"
  text(0.5, 0.5, sprintf("%.2f", r),
       cex = cex_r, col = col_r, font = 2)
}

panel_hist <- function(x, ...) {
  usr <- par("usr"); on.exit(par(usr))
  h   <- hist(x, plot = FALSE, breaks = 10)
  breaks <- h$breaks
  par(usr = c(breaks[1], tail(breaks,1), 0, max(h$counts) * 1.25))
  rect(breaks[-length(breaks)], 0, breaks[-1], h$counts,
       col = "#d4e6da", border = "white")
}

pairs(lakes,
      labels      = c("Temp\n(°C)", "DO\n(mg/L)", "pH", "TP\n(μg/L)", "Chl-a\n(μg/L)"),
      lower.panel = panel_scatter,
      upper.panel = panel_r,
      diag.panel  = panel_hist,
      gap         = 0.35,
      cex.labels  = 1.0,
      font.labels = 2)

Enhanced scatterplot matrix with Pearson r in the upper triangle, loess smooths (amber) in the lower triangle, and variable histograms on the diagonal. The upper-triangle values are colour-coded: dark green = strong positive, red = strong negative, grey ≈ zero. Note the skewed distributions of TP and chl-a on the diagonal, reinforcing the case for Spearman correlation with those variables.

The correlation matrix in numbers

# Pearson correlation matrix — rounded
round(cor(lakes, method = "pearson"), 3)

##        temp_C do_mgl     pH tp_ugl   chla
## temp_C  1.000 -0.829  0.577  0.181 -0.048
## do_mgl -0.829  1.000 -0.587 -0.058  0.041
## pH      0.577 -0.587  1.000  0.161  0.151
## tp_ugl  0.181 -0.058  0.161  1.000  0.202
## chla   -0.048  0.041  0.151  0.202  1.000

# Spearman correlation matrix — better for skewed variables
round(cor(lakes, method = "spearman"), 3)

##        temp_C do_mgl     pH tp_ugl   chla
## temp_C  1.000 -0.835  0.538  0.145 -0.027
## do_mgl -0.835  1.000 -0.518 -0.033  0.006
## pH      0.538 -0.518  1.000  0.152  0.051
## tp_ugl  0.145 -0.033  0.152  1.000  0.271
## chla   -0.027  0.006  0.051  0.271  1.000

R_mat <- cor(lakes, method = "pearson")
vars_lab <- c("Temp", "DO", "pH", "TP", "Chl-a")
p  <- nrow(R_mat)

# Colour mapping
r_to_col <- function(r) {
  if (r > 0) adjustcolor("#1e3a2f", alpha.f = abs(r) * 0.85 + 0.08)
  else       adjustcolor("#7a1e1e", alpha.f = abs(r) * 0.85 + 0.08)
}

par(mar = c(1, 4, 4, 1))
plot(NA, xlim = c(0.5, p + 0.5), ylim = c(0.5, p + 0.5),
     xaxt = "n", yaxt = "n", xlab = "", ylab = "",
     main = "Pearson correlation heatmap — lake chemistry",
     cex.main = 0.91)

for (i in seq_len(p)) {
  for (j in seq_len(p)) {
    r_ij <- R_mat[i, j]
    rect(j - 0.48, p - i + 0.52, j + 0.48, p - i + 1.48,
         col    = if (i == j) "#ede8df" else r_to_col(r_ij),
         border = "white", lwd = 1.5)
    if (i != j) {
      text(j, p - i + 1, sprintf("%.2f", r_ij),
           cex = 0.82, col = if (abs(r_ij) > 0.5) "#fff" else "#1a1a18",
           font = if (abs(r_ij) > 0.6) 2 else 1)
    }
  }
}
axis(3, at = seq_len(p), labels = vars_lab, tick = FALSE, cex.axis = 0.88)
axis(2, at = seq_len(p), labels = rev(vars_lab), tick = FALSE,
     cex.axis = 0.88, las = 1)

Pearson correlation heatmap for the five lake variables. Colour saturation encodes strength; hue encodes direction (forest green = positive, red = negative). The upper triangle shows the r value. Read off-diagonal cells to identify clusters of co-varying variables — temperature, DO, and pH form one cluster; TP and chl-a form a second.

Redundancy and multicollinearity

Variables with |r| > 0.8 are near-redundant: they carry almost the same information. In a regression context this is multicollinearity, which inflates coefficient standard errors. In an exploratory context it tells you that the two variables may be proxies for the same underlying process.

# Which pairs are strongly correlated (|r| > 0.5)?
R_spear <- cor(lakes, method = "spearman")
strong  <- which(abs(R_spear) > 0.5 & upper.tri(R_spear), arr.ind = TRUE)
data.frame(
  var1 = vars_lab[strong[,1]],
  var2 = vars_lab[strong[,2]],
  rho  = round(R_spear[strong], 3)
)

##   var1 var2    rho
## 1 Temp   DO -0.835
## 2 Temp   pH  0.538
## 3   DO   pH -0.518

Spurious correlation via a shared driver

Temperature drives both DO (physics) and pH (photosynthesis rates). The observed positive correlation between pH and DO may therefore be partly or entirely spurious — driven by their shared dependence on temperature rather than any direct chemical link.

Partial correlation removes the linear effect of the confounding variable:

# Partial correlation: pH ~ DO after removing shared variance with temperature
# Residualise both variables on temperature, then correlate residuals

res_pH <- residuals(lm(pH     ~ temp_C, data = lakes))
res_do <- residuals(lm(do_mgl ~ temp_C, data = lakes))

r_partial <- cor.test(res_pH, res_do, method = "pearson")
cat(sprintf(
  "Raw Pearson r(pH, DO)       = %.3f\n",
  cor(lakes$pH, lakes$do_mgl)
))

## Raw Pearson r(pH, DO)       = -0.587

cat(sprintf(
  "Partial r(pH, DO | temp)    = %.3f  (p = %.3f)\n",
  r_partial$estimate,
  r_partial$p.value
))

## Partial r(pH, DO | temp)    = -0.239  (p = 0.079)

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

plot(lakes$do_mgl, lakes$pH,
     pch = 19, cex = 0.85, col = adjustcolor("#1e3a2f", 0.6),
     xlab = "Dissolved oxygen (mg/L)", ylab = "pH",
     main = sprintf("Raw: r = %.2f", cor(lakes$do_mgl, lakes$pH)),
     cex.main = 0.9)
abline(lm(pH ~ do_mgl, data = lakes), col = "#c8832a", lwd = 1.8)

plot(res_do, res_pH,
     pch = 19, cex = 0.85, col = adjustcolor("#2d5a3d", 0.6),
     xlab = "DO residuals (temp removed)", ylab = "pH residuals (temp removed)",
     main = sprintf("Partial: r = %.2f (temp controlled)",
                    r_partial$estimate),
     cex.main = 0.9)
abline(lm(res_pH ~ res_do), col = "#c8832a", lwd = 1.8)

Left: the raw pH × DO scatterplot, which shows a positive association. Right: pH and DO residuals after removing the shared linear effect of temperature. The association weakens substantially, revealing that most of the raw correlation is attributable to temperature as a common driver.

par(mfrow = c(1,1))

Multiple testing in correlation matrices

A 5-variable matrix has \(\binom{5}{2} = 10\) unique pairs; a 10-variable matrix has 45. Testing all pairs at α = 0.05 produces approximately \(0.05 \times 10 = 0.5\) false positives per matrix by chance alone — nearly one spurious significant result for every two matrices examined.

# Adjust p-values for multiple comparisons (Holm method)
p_mat <- matrix(NA, p, p, dimnames = list(vars_lab, vars_lab))
for (i in 1:(p-1)) for (j in (i+1):p) {
  p_mat[i,j] <- cor.test(lakes[[i]], lakes[[j]],
                           method = "pearson")$p.value
}

raw_p  <- p_mat[upper.tri(p_mat)]
adj_p  <- p.adjust(raw_p, method = "holm")

result_df <- data.frame(
  pair    = apply(combn(vars_lab, 2), 2, paste, collapse = " × "),
  r       = round(cor(lakes, method = "pearson")[upper.tri(R_mat)], 3),
  p_raw   = round(raw_p, 4),
  p_holm  = round(adj_p, 4),
  sig_raw = ifelse(raw_p  < 0.05, "✓", ""),
  sig_adj = ifelse(adj_p  < 0.05, "✓", "")
)
result_df

##            pair      r  p_raw p_holm sig_raw sig_adj
## 1     Temp × DO -0.829 0.0000 0.0000       ✓       ✓
## 2     Temp × pH  0.577 0.0000 0.0000       ✓       ✓
## 3     Temp × TP -0.587 0.0000 0.0000       ✓       ✓
## 4  Temp × Chl-a  0.181 0.1870 1.0000                
## 5       DO × pH -0.058 0.6722 1.0000                
## 6       DO × TP  0.161 0.2414 1.0000                
## 7    DO × Chl-a -0.048 0.7280 1.0000                
## 8       pH × TP  0.041 0.7665 1.0000                
## 9    pH × Chl-a  0.151 0.2702 1.0000                
## 10   TP × Chl-a  0.202 0.1398 0.9788

Checklist before submitting a correlation analysis

What’s next?

• Simple Linear Regression — when X causes Y and you need a slope, not just an r
• Tests for Association — chi-square and Fisher’s exact for categorical variables
• Hypothesis Testing — the five-step framework underlying cor.test()