AskStats

Answer

What goes wrong if you sort x and y separately before a regression?

Inspired by a question on Cross Validated ·

regressioncorrelation

The short answer

Sorting x and y independently throws away the one thing a regression studies: which y belongs with which x. The sorted data line up smallest with smallest and largest with largest, so the correlation is always positive and usually close to 1, even when the variables are unrelated or negatively related. The high R-squared and tiny p-value are artefacts of the sorting, and the fitted line predicts new observations worse than simply using the mean.

What sorting actually does to the data

A regression of y on x asks how y changes as x changes within the same unit: the same person, plot, day or machine. That information lives entirely in the pairing. Sorting the x column and the y column separately keeps both lists of values but deletes the pairing, and replaces it with an artificial one: the smallest x is matched with the smallest y, the second smallest with the second smallest, and so on.

What you are left with is a plot of the quantiles of y against the quantiles of x. That is exactly a Q-Q plot of the two distributions. It tells you whether x and y have similar shapes (a straight Q-Q line means one is a shifted and rescaled version of the other), but it says nothing about whether they are related.

See it in R and Python

The code below simulates two datasets of 50 observations: one where x and y are unrelated, and one where y falls as x rises (true correlation -0.5). It fits a regression to each before and after sorting, then uses the second dataset's two fitted lines to predict 10,000 new observations from the same process.

R

set.seed(185507)
n <- 50

# Case 1: x and y unrelated (true correlation 0)
x <- rnorm(n); y <- rnorm(n)
# Case 2: y falls as x rises (true correlation about -0.5)
x2 <- rnorm(n); y2 <- -0.5 * x2 + rnorm(n, sd = sqrt(0.75))

summ <- function(x, y) {
  f <- lm(y ~ x)
  c(r = cor(x, y), slope = coef(f)[[2]], R2 = summary(f)$r.squared,
    p = coef(summary(f))[2, 4])
}
round(rbind(
  unrelated_original = summ(x, y),
  unrelated_sorted   = summ(sort(x), sort(y)),
  negative_original  = summ(x2, y2),
  negative_sorted    = summ(sort(x2), sort(y2))
), 4)

# Out-of-sample prediction for case 2: new data from the same process
xn <- rnorm(10000); yn <- -0.5 * xn + rnorm(10000, sd = sqrt(0.75))
f_orig <- lm(y2 ~ x2)
f_sort <- lm(sy ~ sx, data = data.frame(sx = sort(x2), sy = sort(y2)))
mse <- function(pred) mean((yn - pred)^2)
round(c(
  original_fit = mse(coef(f_orig)[1] + coef(f_orig)[2] * xn),
  sorted_fit   = mse(coef(f_sort)[1] + coef(f_sort)[2] * xn),
  just_mean    = mse(rep(mean(y2), length(xn)))
), 3)

Python

import numpy as np
from scipy import stats

rng = np.random.default_rng(185507)
n = 50

# Case 1: x and y unrelated (true correlation 0)
x, y = rng.normal(size=n), rng.normal(size=n)
# Case 2: y falls as x rises (true correlation about -0.5)
x2 = rng.normal(size=n)
y2 = -0.5 * x2 + rng.normal(0, np.sqrt(0.75), n)

def summ(name, x, y):
    f = stats.linregress(x, y)
    print(f"{name:19s} r={f.rvalue:7.4f} slope={f.slope:7.4f} "
          f"R2={f.rvalue**2:.4f} p={f.pvalue:.4f}")

summ("unrelated_original", x, y)
summ("unrelated_sorted", np.sort(x), np.sort(y))
summ("negative_original", x2, y2)
summ("negative_sorted", np.sort(x2), np.sort(y2))

# Out-of-sample prediction for case 2: new data from the same process
xn = rng.normal(size=10000)
yn = -0.5 * xn + rng.normal(0, np.sqrt(0.75), 10000)
f_orig = stats.linregress(x2, y2)
f_sort = stats.linregress(np.sort(x2), np.sort(y2))
mse = lambda pred: np.mean((yn - pred) ** 2)
print(f"original_fit={mse(f_orig.intercept + f_orig.slope * xn):.3f} "
      f"sorted_fit={mse(f_sort.intercept + f_sort.slope * xn):.3f} "
      f"just_mean={mse(np.full(xn.size, y2.mean())):.3f}")
                         r   slope     R2      p
unrelated_original  0.0529  0.0463 0.0028 0.7152
unrelated_sorted    0.9495  0.8315 0.9016 0.0000
negative_original  -0.4597 -0.5090 0.2113 0.0008
negative_sorted     0.9932  1.0999 0.9865 0.0000
original_fit   sorted_fit    just_mean 
       0.762        3.388        1.029

The output above is from one seeded run of the R code. Python uses a different random number generator, so its figures differ slightly, but the pattern is the same.

Reading the results

The p-values after sorting are not merely too small; they have no meaning at all. The test assumes each row is an independent observation of a real (x, y) pair. Sorted rows are neither: each one depends on the ranks of every other value in the sample.

Why sorting can seem to give better regressions

The fit statistics printed for the sorted data are computed on the sorted data, so they only measure how straight the Q-Q plot is. Since almost any two roughly similar distributions give a nearly straight Q-Q plot, sorting will almost always look like an improvement by R-squared, residual error or p-value. None of those numbers tell you how well the model predicts a new unit whose x you know and whose y you do not.

The honest check is out-of-sample: fit each approach on part of the data and score it on held-out pairs that were never sorted. As the simulation shows, the sorted fit loses badly whenever the true relationship is weak or negative, and at best matches the ordinary fit when the relationship is already very strong and positive.

When pairing sorted values is legitimate

Matching quantiles is a real technique when the goal is to compare or map distributions, not to model a relationship between paired measurements. A Q-Q plot compares two samples' shapes. Equipercentile equating maps scores on one test form onto the scale of another by matching percentiles, and quantile mapping corrects the distribution of one data source to match another. None of these claims that a unit with a high x will have a high y. For a regression, keep the rows together, and check the data before modelling (our data audit checklist covers this). Our page on reading Q-Q plots explains what a sorted-against-sorted plot does show.

How to report a Pearson correlation in APA style (7th edition)

Report the correlation or regression computed on the original pairs, never on separately sorted values. Italicize r and p, give the degrees of freedom (n - 2) in parentheses, and drop the leading zero because r cannot exceed 1. Using the negative-relationship example above (n = 50):

"There was a moderate negative correlation between x and y, r(48) = -.46, p < .001. In a simple linear regression, each one-unit increase in x was associated with a 0.51-unit decrease in y, b = -0.51."

For a null result, report the exact p-value, for example "r(48) = .05, p = .715". Where space allows, add a 95% confidence interval for r in square brackets, [LL, UL].

Related tools and guides

More answered questions

Working with your own data?

General answers only go so far. Send us your situation and we'll reply by email within two business days.

Ask your question

Written by AskStats with AI assistance. This is general information, not advice for your specific data or study. When the results matter, check your approach with a qualified statistician.