Correlation

Correlation is the first question you ask about two variables: do taller people weigh more, do ad budgets track sales, do two stocks move together? It's the number behind diversification in a portfolio, feature screening before a model, and the "r" that simple regression reports. Module 7's regression slope is the correlation rescaled into the data's units.

It's also where the most famous warnings in statistics live: correlation isn't causation, and one number can't describe a scatter plot.

Covariance averages the product of the two variables' deviations from their means. When both are above (or both below) their means at once, the product is positive; when one is up while the other is down, it's negative. Its size depends on the units, so on its own it's hard to read.

Pearson's r is covariance divided by both standard deviations, which strips out the units and pins it between −1 and 1. ±1 means the points lie exactly on a line; 0 means no linear relationship. Its square, r², is the share of one variable's variance a straight line through the other explains.

Testing r. With n points, t = r·√((n − 2)/(1 − r²)) follows a t distribution with n − 2 degrees of freedom when the true correlation is zero. Small samples need big correlations to be convincing: r = 0.6 isn't significant at 0.05 with 8 points, but r = 0.2 is with 100.

Spearman's ρ is Pearson's r computed on the ranks of the data. It asks "does y tend to rise when x rises?" without requiring a straight line, and a single wild value can only move one rank. A curve that always rises (like y = x²) has ρ = 1 but r < 1.

Correlation is not causation. Ice-cream sales and drownings correlate because summer drives both — a confounder. Two variables that both trend over time correlate no matter what. A correlation is a reason to investigate, not a conclusion.

• Check lengths before pairing slices. x[i] and y[i] must describe the same observation; return NaN (or an error) when len(x) != len(y) rather than silently using the shorter one or panicking.
• Ranking without losing your place: sort a slice of indices by the values (slices.SortStableFunc), then walk runs of equal values and assign each run the average rank. The result goes back by index, so it stays in input order.
• Reuse beats re-derivation. Pearson is covariance over two standard deviations, and each standard deviation is √Covariance(x, x). Spearman is Pearson(Ranks(x), Ranks(y)) — one line.
• A correlation matrix is a [][]float64 of size k×k. It's symmetric with 1s on the diagonal, so you could compute only half — but at this size, clarity wins.

package main

import (
	"fmt"
	"math"
)

// fit returns Pearson's r and the least-squares line y = a + b·x.
func fit(x, y []float64) (r, a, b float64) {
	n := float64(len(x))
	var mx, my float64
	for i := range x {
		mx += x[i]
		my += y[i]
	}
	mx, my = mx/n, my/n
	var sxy, sxx, syy float64
	for i := range x {
		sxy += (x[i] - mx) * (y[i] - my)
		sxx += (x[i] - mx) * (x[i] - mx)
		syy += (y[i] - my) * (y[i] - my)
	}
	b = sxy / sxx
	return sxy / math.Sqrt(sxx*syy), my - b*mx, b
}

func main() {
	// Anscombe's quartet (1973): four data sets built to share their statistics.
	x := []float64{10, 8, 13, 9, 11, 14, 6, 4, 12, 7, 5}
	sets := []struct {
		x, y []float64
	}{
		{x, []float64{8.04, 6.95, 7.58, 8.81, 8.33, 9.96, 7.24, 4.26, 10.84, 4.82, 5.68}},
		{x, []float64{9.14, 8.14, 8.74, 8.77, 9.26, 8.10, 6.13, 3.10, 9.13, 7.26, 4.74}},
		{x, []float64{7.46, 6.77, 12.74, 7.11, 7.81, 8.84, 6.08, 5.39, 8.15, 6.42, 5.73}},
		{[]float64{8, 8, 8, 8, 8, 8, 8, 19, 8, 8, 8},
			[]float64{6.58, 5.76, 7.71, 8.84, 8.47, 7.04, 5.25, 12.50, 5.56, 7.91, 6.89}},
	}
	for i, s := range sets {
		r, a, b := fit(s.x, s.y)
		fmt.Printf("set %d:  r = %.3f   y = %.2f + %.3f·x\n", i+1, r, a, b)
	}
}

It prints:

Four data sets, one correlation, one regression line — and they look nothing alike. Set 1 is a noisy straight line, the only one where r and the line tell the truth. Set 2 is a perfect curve. Set 3 is a perfect straight line knocked off course by one outlier (its Spearman ρ is 0.991). Set 4 has ten points at x = 8 and one at x = 19, which single-handedly decides the slope. The statistician Francis Anscombe built them in 1973 to make one point: always plot your data.

• Reading r = 0 as "unrelated". Pearson only sees straight lines. y = x² over −3…3 has r = 0 and a perfect relationship.
• Inferring causation. Without an experiment or a careful argument about confounders, a correlation only says the variables move together.
• Correlating trends. Any two series that both grow over time — population and internet use, say — correlate strongly. Correlate the year-to-year changes instead.
• Trusting r from a handful of points. With five points, r = 0.8 has a p-value of about 0.10. Report the p-value or a confidence interval alongside r.
• Pooling groups. Two groups that each show no correlation can show a strong one when combined, if their averages differ (and the reverse — Simpson's paradox).

Part A builds covariance, Pearson's r, and its significance test. Part B is ranking with ties and Spearman's ρ. Part C builds a correlation matrix for several variables at once.

go test ./06-correlation/

TestSpearman uses Anscombe's set 3: Pearson gives 0.816, Spearman 0.991. If your Ranks doesn't average ties, the third case (with ties in both variables) will catch it.

The files for this module

Everything below is also in the course download, in 06-correlation/. Work in exercises.go; the tests decide when you’re done.

// Package correlation is Module 6: covariance and correlation.
//
// Correlation measures how tightly two variables move together on a line
// (Pearson) or in the same order (Spearman). It summarizes a scatter plot in
// one number — which is exactly why you should always look at the plot.
package correlation

import "github.com/bozcode/go-statistics-and-regression/internal/special"

// The helpers below are given — you wrote them in Modules 1 and 4.

// tCDF returns P(T <= t) for Student's t with df degrees of freedom.
func tCDF(t, df float64) float64 { return special.TCDF(t, df) }

// mean returns the arithmetic mean of xs.
func mean(xs []float64) float64 {
	sum := 0.0
	for _, x := range xs {
		sum += x
	}
	return sum / float64(len(xs))
}

//go:build !solutions

package correlation

// ============================================================
// Part A — covariance and Pearson's r
// ============================================================

// Covariance returns the sample covariance Σ(x − x̄)(y − ȳ) / (n − 1):
// positive when the variables rise together, negative when one falls as the
// other rises. NaN if the lengths differ or n < 2.
func Covariance(x, y []float64) float64 {
	// TODO
	return 0
}

// Pearson returns the correlation coefficient r = cov(x, y) / (sx·sy),
// always between −1 and 1. NaN if either variable is constant (its standard
// deviation is 0) or the inputs are unusable.
func Pearson(x, y []float64) float64 {
	// TODO
	return 0
}

// PearsonTest returns r and the two-sided p-value for H₀: ρ = 0, using
// t = r·√((n − 2)/(1 − r²)) with n − 2 degrees of freedom and the given
// tCDF. If |r| is 1 the p-value is 0. Matches SciPy's pearsonr.
func PearsonTest(x, y []float64) (r, p float64) {
	// TODO
	return 0, 0
}

// ============================================================
// Part B — rank correlation
// ============================================================

// Ranks returns the rank of each value, 1 for the smallest. Tied values
// share the average of the ranks they span: [10 20 20 30] -> [1 2.5 2.5 4].
// Results are in the input's order. Matches SciPy's rankdata.
func Ranks(xs []float64) []float64 {
	// TODO
	return nil
}

// Spearman returns Spearman's rank correlation: Pearson's r computed on the
// ranks of x and the ranks of y. It measures whether the relationship is
// monotonic (always rising or always falling), not whether it's a straight
// line, and an outlier can only move it so far.
func Spearman(x, y []float64) float64 {
	// TODO
	return 0
}

// ============================================================
// Part C — many variables at once
// ============================================================

// Matrix returns the correlation matrix of several variables (each column
// is one variable's values): m[i][j] = Pearson(cols[i], cols[j]), with 1 on
// the diagonal.
func Matrix(cols [][]float64) [][]float64 {
	// TODO
	return nil
}

package correlation

import (
	"math"
	"slices"
	"testing"
)

// near reports whether a and b agree to about 9 significant digits.
// Expected values come from NumPy and SciPy.
func near(a, b float64) bool {
	return math.Abs(a-b) <= 1e-9*math.Max(1, math.Abs(b))
}

// Anscombe's quartet, sets 1 and 3 (they share x).
var (
	ansX  = []float64{10, 8, 13, 9, 11, 14, 6, 4, 12, 7, 5}
	ansY1 = []float64{8.04, 6.95, 7.58, 8.81, 8.33, 9.96, 7.24, 4.26, 10.84, 4.82, 5.68}
	ansY3 = []float64{7.46, 6.77, 12.74, 7.11, 7.81, 8.84, 6.08, 5.39, 8.15, 6.42, 5.73}
)

// heights (inches) and weights (pounds), made up
var (
	height = []float64{62, 65, 68, 70, 71, 73, 75}
	weight = []float64{118, 124, 150, 142, 160, 171, 168}
)

func TestCovariance(t *testing.T) {
	if got := Covariance(ansX, ansY1); !near(got, 5.501) {
		t.Errorf("Covariance(x, y1) = %v, want 5.501", got)
	}
	if got := Covariance(ansX, ansY3); !near(got, 5.497) {
		t.Errorf("Covariance(x, y3) = %v, want 5.497", got)
	}
	if got := Covariance(ansX, ansX); !near(got, 11) {
		t.Errorf("Covariance(x, x) = %v, want 11 (the variance of x)", got)
	}
	if got := Covariance([]float64{1, 2}, []float64{1, 2, 3}); !math.IsNaN(got) {
		t.Errorf("Covariance of different lengths = %v, want NaN", got)
	}
}

func TestPearson(t *testing.T) {
	if got := Pearson(ansX, ansY1); !near(got, 0.8164205163448395) {
		t.Errorf("Pearson(x, y1) = %v, want 0.8164205163448395", got)
	}
	if got := Pearson(height, weight); !near(got, 0.9478318558142701) {
		t.Errorf("Pearson(height, weight) = %v, want 0.9478318558142701", got)
	}
	if got := Pearson([]float64{1, 2, 3, 4}, []float64{8, 6, 4, 2}); !near(got, -1) {
		t.Errorf("Pearson of a perfect falling line = %v, want -1", got)
	}
	if got := Pearson([]float64{1, 2, 3}, []float64{5, 5, 5}); !math.IsNaN(got) {
		t.Errorf("Pearson with a constant variable = %v, want NaN", got)
	}
}

func TestPearsonTest(t *testing.T) {
	for _, tc := range []struct {
		x, y []float64
		r, p float64
	}{
		{ansX, ansY1, 0.8164205163448395, 0.0021696288730788152},
		{height, weight, 0.9478318558142701, 0.0011607019316604108},
		{[]float64{1, 2, 3, 4, 5}, []float64{3.2, 2.9, 2.1, 1.8, 1.0}, -0.9878291611472619, 0.0016088671217392764},
	} {
		r, p := PearsonTest(tc.x, tc.y)
		if !near(r, tc.r) || !near(p, tc.p) {
			t.Errorf("PearsonTest(%v, %v) = (%v, %v), want (%v, %v)", tc.x, tc.y, r, p, tc.r, tc.p)
		}
	}
}

func TestRanks(t *testing.T) {
	for _, tc := range []struct{ in, want []float64 }{
		{[]float64{10, 20, 20, 30}, []float64{1, 2.5, 2.5, 4}},
		{[]float64{3, 1, 4, 1, 5, 9, 2, 6, 5, 3, 5}, []float64{4.5, 1.5, 6, 1.5, 8, 11, 3, 10, 8, 4.5, 8}},
		{[]float64{7}, []float64{1}},
	} {
		if got := Ranks(tc.in); !slices.Equal(got, tc.want) {
			t.Errorf("Ranks(%v) = %v, want %v", tc.in, got, tc.want)
		}
	}
}

func TestSpearman(t *testing.T) {
	for _, tc := range []struct {
		x, y []float64
		want float64
	}{
		{ansX, ansY3, 0.990909090909091},                                 // one outlier drags Pearson to 0.816; ranks barely notice
		{[]float64{1, 2, 3, 4, 5, 6}, []float64{1, 4, 9, 16, 25, 36}, 1}, // a curve, but perfectly monotonic
		{[]float64{1, 2, 2, 3, 4}, []float64{10, 20, 30, 30, 50}, 0.9210526315789475},
	} {
		if got := Spearman(tc.x, tc.y); !near(got, tc.want) {
			t.Errorf("Spearman(%v, %v) = %v, want %v", tc.x, tc.y, got, tc.want)
		}
	}
}

func TestMatrix(t *testing.T) {
	age := []float64{30, 31, 29, 33, 35, 34, 36}
	want := [][]float64{
		{1, 0.9478318558142701, 0.8301425586780717},
		{0.9478318558142701, 1, 0.7174395873511835},
		{0.8301425586780717, 0.7174395873511835, 1},
	}
	got := Matrix([][]float64{height, weight, age})
	if len(got) != 3 {
		t.Fatalf("Matrix returned %d rows, want 3", len(got))
	}
	for i := range want {
		if len(got[i]) != 3 {
			t.Fatalf("row %d has %d entries, want 3", i, len(got[i]))
		}
		for j := range want[i] {
			if !near(got[i][j], want[i][j]) {
				t.Errorf("Matrix[%d][%d] = %v, want %v", i, j, got[i][j], want[i][j])
			}
		}
	}
}

//go:build solutions

package correlation

import (
	"math"
	"slices"
)

// ============================================================
// Part A — covariance and Pearson's r
// ============================================================

func Covariance(x, y []float64) float64 {
	if len(x) != len(y) || len(x) < 2 {
		return math.NaN()
	}
	mx, my := mean(x), mean(y)
	s := 0.0
	for i := range x {
		s += (x[i] - mx) * (y[i] - my)
	}
	return s / float64(len(x)-1)
}

func Pearson(x, y []float64) float64 {
	sx, sy := math.Sqrt(Covariance(x, x)), math.Sqrt(Covariance(y, y))
	if sx == 0 || sy == 0 {
		return math.NaN()
	}
	return Covariance(x, y) / (sx * sy)
}

func PearsonTest(x, y []float64) (r, p float64) {
	r = Pearson(x, y)
	if math.Abs(r) >= 1 {
		return r, 0
	}
	df := float64(len(x) - 2)
	t := r * math.Sqrt(df/(1-r*r))
	return r, 2 * tCDF(-math.Abs(t), df)
}

// ============================================================
// Part B — rank correlation
// ============================================================

func Ranks(xs []float64) []float64 {
	idx := make([]int, len(xs))
	for i := range idx {
		idx[i] = i
	}
	slices.SortStableFunc(idx, func(a, b int) int {
		switch {
		case xs[a] < xs[b]:
			return -1
		case xs[a] > xs[b]:
			return 1
		}
		return 0
	})
	out := make([]float64, len(xs))
	for i := 0; i < len(idx); {
		j := i
		for j < len(idx) && xs[idx[j]] == xs[idx[i]] {
			j++
		}
		avg := float64(i+1+j) / 2 // ranks i+1 … j, averaged
		for k := i; k < j; k++ {
			out[idx[k]] = avg
		}
		i = j
	}
	return out
}

func Spearman(x, y []float64) float64 {
	return Pearson(Ranks(x), Ranks(y))
}

// ============================================================
// Part C — many variables at once
// ============================================================

func Matrix(cols [][]float64) [][]float64 {
	m := make([][]float64, len(cols))
	for i := range cols {
		m[i] = make([]float64, len(cols))
		for j := range cols {
			if i == j {
				m[i][j] = 1
			} else {
				m[i][j] = Pearson(cols[i], cols[j])
			}
		}
	}
	return m
}

• Add a bootstrap confidence interval for r using Module 3's Bootstrap (resample pairs, not each variable separately).
• Add Kendall's τ, another rank correlation: count concordant versus discordant pairs.
• Practice: compute the correlation matrix of a few columns from a CSV you care about, then plot the strongest pair to see whether r told the truth.

Four questions. Answer each one in your head — or out loud — before you open the answer.

Commit your work before you move on to module 7, so git log shows your own progress through the course.