Pearson vs Spearman correlation: which to use

Use Pearson’s r when the relationship between two continuous variables is roughly a straight line and no single point dominates the scatterplot. Use Spearman’s rho when the relationship rises or falls steadily but not in a straight line, when a variable is ordinal, or when extreme values are a concern. This post computes both on the same real data with confidence intervals, then shows a seeded example in which one stray point reverses the sign of Pearson’s r while Spearman’s rho keeps its sign.

What each coefficient measures

Pearson’s r is the covariance of two variables divided by the product of their standard deviations. It reaches -1 or +1 only when every point lies exactly on a straight line, so it answers one question: how well does a line describe this scatter? Because it uses the raw values, a point far from the others pulls on it hard.

Spearman’s rho (Spearman, 1904) is Pearson’s r computed on ranks. Each variable is ranked from smallest to largest, tied values share the average of their ranks, and the ordinary correlation is taken of the two rank columns. It reaches -1 or +1 whenever one variable always rises (or always falls) as the other rises, even along a curve. That property is called monotonic, and it is the whole difference in what the two coefficients claim.

Two consequences follow, and the rest of this post rests on them. First, a curved but steadily rising relationship scores a high rho and a lower r. Second, ranking throws away distances, so a very distant point counts no more than a slightly distant one. That is useful when the distance is a fluke and a loss when the distance is real information.

The data: Old Faithful eruptions

R ships the faithful data: 272 eruptions of the Old Faithful geyser, with the eruption duration in minutes (eruptions) and the waiting time to the next eruption in minutes (waiting). The question is simple: do longer eruptions come before longer waits? Each row is one eruption, so the design is two continuous measurements per unit, which is what a correlation is for. One caution applies throughout: eruptions of one geyser were recorded over time, and neighbouring eruptions may resemble each other, so the independence the standard tests assume is only approximate.

Plot the data before computing anything, because a coefficient hides the shape of the scatter (Anscombe, 1973). The faithful scatterplot is not one smooth cloud. It forms two groups, one of short eruptions followed by short waits and one of long eruptions followed by long waits, with few points between. That shape is exactly what makes the two coefficients disagree.

plot(faithful$eruptions, faithful$waiting,
     xlab = "Eruption duration (minutes)",
     ylab = "Wait to next eruption (minutes)")

What the engine ran

The deterministic engine ran Spearman’s correlation on the 272 pairs, from a fixed seed:

• Spearman’s rho = 0.778, n = 272, p < 0.001
• Engine assumption check: sample size adequate for the method
• The engine reports no confidence interval for this method, so the intervals below are base R

Both columns contain tied values, and R warns that it cannot compute an exact p-value with ties. The engine’s p-value therefore comes from the large-sample t approximation that R falls back to (R documentation for cor.test), which is reasonable at 272 pairs. Everything from here on, apart from this one result, is base R code and its printed output, not engine output. The example script on this page reruns the engine part so you can compare it with the figures above.

Both coefficients with intervals

Pearson’s r comes with a confidence interval from cor.test(), which uses Fisher’s z transformation (Fisher, 1921; R documentation for cor.test). That interval leans on the pairs being roughly bivariate normal, and the two-group shape here is not. So the code below also gives a percentile bootstrap interval for each coefficient: resample the 272 pairs with replacement 5000 times, recompute the coefficient each time, and take the 2.5th and 97.5th percentiles (Efron, 1979). The bootstrap also covers Spearman, for which cor.test() gives no interval (R documentation for cor.test).

boot_ci <- function(x, y, method, B = 5000) {
  idx <- replicate(B, sample(length(x), replace = TRUE))
  r <- apply(idx, 2, function(i) cor(x[i], y[i], method = method))
  quantile(r, c(0.025, 0.975))
}

d <- datasets::faithful
set.seed(20261025)
cor(d$eruptions, d$waiting)
cor.test(d$eruptions, d$waiting)$conf.int
boot_ci(d$eruptions, d$waiting, "pearson")
boot_ci(d$eruptions, d$waiting, "spearman")

The results, to three decimals:

• Pearson’s r = 0.901, Fisher interval 0.876 to 0.921, bootstrap interval 0.883 to 0.917
• Spearman’s rho = 0.778, bootstrap interval 0.719 to 0.821

The two intervals do not overlap. On these data the choice of coefficient changes the headline number by 0.12, and both statements are true of the same scatter. They describe different things: r says the points are tightly arranged around a line, rho says that eruption rank predicts wait rank less tightly.

The Fisher and bootstrap intervals for Pearson’s r agree closely here, 0.876 to 0.921 against 0.883 to 0.917, so the departure from bivariate normality did little harm to the interval at this sample size. At 272 pairs that is reassuring, but it is a property of this sample and not a guarantee. With 20 pairs the two methods can disagree much more, and a bootstrap interval from so few pairs is itself unstable, so say how many pairs you had.

Why Spearman is lower here

Pearson gives the gap between the two groups its full distance, and that gap is large compared with the spread inside each group. Spearman turns the gap into a step in the ranks and lets the scatter inside each group weigh more. Inside each group, the correlation is much weaker. Among the 98 eruptions shorter than 3.1 minutes, r is 0.351 and rho is 0.273. Among the other 174, r is 0.354 and rho is 0.283.

short <- d$eruptions < 3.1
within <- function(i) round(c(n = sum(i),
  r = cor(d$eruptions[i], d$waiting[i]),
  rho = cor(d$eruptions[i], d$waiting[i], method = "spearman")), 3)
within(short)    # n 98, r 0.351, rho 0.273
within(!short)   # n 174, r 0.354, rho 0.283

So most of the 0.901 comes from there being two groups, and a reader who takes it as “a very tight line” has overread it. Neither coefficient is wrong. Reporting only the larger one would mislead, and this is the case for plotting first and for naming the pattern in the text.

This also corrects a common rule of thumb. Spearman is often called the safe, conservative choice, and in this example it is lower, but not because it is safer: ranking changed what is being summarised. There is no general rule that rho is smaller or larger than r.

One outlier can flip Pearson

The second example is synthetic and seeded, so you can run it to the digit. Thirty values of x are drawn uniformly between 1 and 10 (observed range 1.2 to 9.6), and y is 2 plus 1.2 times x plus normal noise with a standard deviation of 2.5; x and y are rounded to one decimal. The relationship is positive and close to a line. Then one extra point is added at x = 30, y = -30, the kind of value that a unit mix-up or a typing slip produces.

fit_both <- function(x, y) {
  p <- cor.test(x, y)
  s <- suppressWarnings(cor.test(x, y, method = "spearman", exact = FALSE))
  set.seed(1)
  ci <- boot_ci(x, y, "spearman")
  c(pearson = p$estimate, p_lo = p$conf.int[1], p_hi = p$conf.int[2],
    spearman = s$estimate, s_lo = ci[[1]], s_hi = ci[[2]])
}

set.seed(20261025)
n <- 30
x <- round(runif(n, 1, 10), 1)
y <- round(2 + 1.2 * x + rnorm(n, 0, 2.5), 1)
fit_both(x, y)
fit_both(c(x, 30), c(y, -30))

Without the extra point (n = 30), Pearson’s r is 0.904 (95% interval 0.807 to 0.954) and Spearman’s rho is 0.883 (bootstrap interval 0.711 to 0.951). With it (n = 31), Pearson’s r is -0.575 (95% interval -0.772 to -0.277, p < 0.001) and Spearman’s rho is 0.706 (bootstrap interval 0.308 to 0.936).

One point out of 31 turned a strong positive correlation into a “statistically significant” negative one. Spearman kept its sign because the new point only has the largest x rank and the smallest y rank, and ranks are bounded: the point is as discordant as a point can be, but no more. Note what did not hold. Rho fell from 0.883 to 0.706 and its interval widened, so the outlier is not harmless to Spearman, only far less damaging.

The lesson is not “use Spearman and ignore outliers”. The lesson is to look at the point. If it is an error, correct it or remove it for a stated reason. If it is a real observation, the sign flip from Pearson is information: a single case is steering the linear fit, and the honest report says so. The post on checking outliers before a statistical test covers how to decide.

How to choose

Start from the design and the scatterplot, not from a normality test. Some software, Celeus included, screens for normality and suggests Spearman; treat that as a prompt to look at the plot, not as the decision.

• Both variables continuous, relationship close to a straight line, no point that dominates: Pearson, with its interval.
• Relationship consistently rising or falling but curved, for instance a saturating dose response: Spearman. A log transform followed by Pearson is an alternative if the transformed scale is meaningful.
• An ordinal variable, such as a single item on a five-point agreement scale: Spearman, because ordinal values carry order but not equal spacing. The Likert scale data analysis post covers the wider question.
• Extreme values you cannot explain: report both, show the scatterplot and say what you did with the point.
• A relationship that goes up and then down: neither coefficient describes it well, because both measure a one-directional trend. Model the curve instead.

Do not pick the coefficient by which p-value is smaller. Choose from the design before looking, then report the other as a sensitivity check if it matters.

Common mistakes

Switching to Spearman because Pearson came out low. A low r can mean no relationship, a relationship that is not one-directional, or a scatter that is simply wide. A U-shaped pattern lowers both coefficients, so switching changes nothing. Plot first, and switch only when the plot shows a steady curve or an ordinal variable.

Reading a coefficient as a slope. Neither r nor rho says how much y changes for one unit of x. That is the job of a regression coefficient, and the same r can go with a steep or a shallow line. The interpreting regression coefficients post picks that question up.

Comparing coefficients across samples with different spreads. When a sample covers only part of the range of a variable, the correlation shrinks even if the underlying relationship has not changed. The faithful data show a related effect: the correlation is 0.901 over all 272 eruptions and about 0.35 within each group, because most of the 0.901 comes from the gap between the groups rather than from the spread within either. A coefficient from a narrow sample and one from a wide sample are not directly comparable.

Squaring rho and calling it variance explained. The square of Pearson’s r is the share of variance in y that a straight line accounts for. The square of rho is a statement about ranks, not about the original values, so quote it as such or leave it out.

Dropping the inconvenient point. Removing an outlier because it spoils the result is a way of choosing the answer. Decide the rule for handling extreme values before you see the coefficient, write it down, and apply it to both.

The assumptions, stated plainly

Pearson’s coefficient is a descriptive number that needs no distribution at all. Its test and Fisher interval assume the pairs are roughly bivariate normal, and simulation work compares how the tests behave when that fails (Bishara and Hittner, 2012). Spearman’s test assumes no particular distribution.

Both tests assume the pairs are independent of one another, which repeated measures on the same unit break, and neither is a test of linearity or of normality. A significant result says only that the data are unlikely if the true association were zero. It does not say the relationship is strong, linear or causal. For that reading, see correlation versus causation.

Ties need a word. Tied values, which both faithful columns have, are given average ranks, and R then uses its large-sample t approximation for the p-value (R documentation for cor.test). That is a small matter at 272 pairs and a bigger one at ten.

Reporting the result

State the coefficient, the sample size, an interval and the method for it, and the shape you saw. For the first example, something like: “Eruption duration and waiting time were positively associated (Spearman’s rho = 0.778, 95% bootstrap interval 0.719 to 0.821, n = 272); the scatterplot showed two groups of eruptions.” Give the interval for the coefficient you chose and not only its p-value. The r versus R squared post explains how the two relate, and the wider modelling context sits in the regression analysis guide.

If you report both, say which was planned. Reporting only the one that looks better is a form of selective reporting, whichever way the difference runs. Everything above is reproducible from the seed and the faithful data, in the spirit of the reproducible statistical analysis guide and the advice in choosing a statistical test.

Where this goes next

Kendall’s tau is a third rank-based coefficient. It counts concordant and discordant pairs rather than correlating ranks, its numbers sit on a different scale from rho, and it is sometimes preferred when there are many ties. Do not compare its magnitude with r or rho. If you also want a matrix of many pairs, with adjustment for testing them all, see the correlation matrix in R post.

The reference pages for Pearson correlation and Spearman correlation document how each method is run in the engine, including how its p-value is computed, the interval Pearson reports and the checks each reports.

# Results

- Test: Spearman correlation
- n: 272
- statistic: 744700
- p-value: < 0.001
- effect size: Spearman's rho = 0.778
- seed: 20261025

## Assumption checks
- sample size: passed=TRUE (adequate n)

## Power
- Power is not calculated for this test. The matching parametric test gives a rough guide.

# Celeus - reproduce.R  (run from this package folder; bundled data in data/)
# Test: Spearman correlation  |  engine 0.31.8  |  catalog c39-d3cd
# Research context: not provided for this run.
# Reproduces the recorded result deterministically.

set.seed(20261025)
.data_path <- "data/pearson-vs-spearman.csv"
if (!file.exists(.data_path)) stop("Run reproduce.R from the package folder (the bundled data is in data/).")
data <- readr::read_csv(.data_path, show_col_types = FALSE)  # needs the readr package
clean_cells <- function (df) 
{
    ctrl_collapse <- "[\r\n\t\v\f]+"
    zero_width_chars <- intToUtf8(c(8203L, 65279L), multiple = TRUE)
    uni_space_chars <- intToUtf8(c(160L, 8239L, 8199L), multiple = TRUE)
    clean_text <- function(x) {
        out <- x
        for (ch in zero_width_chars) out <- gsub(ch, "", out, 
            fixed = TRUE)
        out <- gsub(ctrl_collapse, " ", out)
        for (ch in uni_space_chars) out <- gsub(ch, " ", out, 
            fixed = TRUE)
        trimws(gsub(" +", " ", out))
    }
    changed_cols <- character(0)
    cells_cleaned <- 0L
    retyped <- character(0)
    for (j in seq_along(df)) {
        col <- df[[j]]
        if (!is.character(col)) 
            next
        cleaned <- clean_text(col)
        n_changed <- sum(!is.na(col) & !is.na(cleaned) & col != 
            cleaned, na.rm = TRUE)
        if (n_changed > 0L) {
            changed_cols <- c(changed_cols, names(df)[j])
            cells_cleaned <- cells_cleaned + n_changed
        }
        df[[j]] <- cleaned
        non_blank <- cleaned[!is.na(cleaned) & nzchar(cleaned)]
        if (length(non_blank) && all(grepl("^[+-]?(0|[1-9][0-9]*)([.][0-9]+)?$", 
            non_blank))) {
            parsed <- suppressWarnings(as.numeric(cleaned))
            was_blank <- is.na(cleaned) | !nzchar(cleaned)
            lost <- (is.na(parsed) | is.infinite(parsed)) & !was_blank
            if (!any(lost)) {
                finite <- parsed[!is.na(parsed)]
                all_whole <- all(finite == floor(finite))
                in_range <- all(abs(finite) <= .Machine$integer.max)
                df[[j]] <- if (all_whole && in_range) 
                  as.integer(parsed)
                else parsed
                retyped <- c(retyped, names(df)[j])
            }
        }
    }
    names(df) <- clean_text(names(df))
    if (length(changed_cols)) {
        attr(df, "ingest_clean_columns") <- changed_cols
        attr(df, "ingest_clean_cells_cleaned") <- cells_cleaned
    }
    if (length(retyped)) 
        attr(df, "ingest_clean_retyped") <- retyped
    df
}
data <- clean_cells(data)
# Complete-case analysis: drop rows missing any analysis value.
.acols <- c("eruptions", "waiting")
data <- data[stats::complete.cases(data[, .acols, drop = FALSE]), , drop = FALSE]
result <- cor.test(data[["eruptions"]], data[["waiting"]], method = "spearman")
print(result)