Regression Diagnostics

Least squares will fit a line through anything: a curve, a cloud with one wild point, predictors that are copies of each other. The output looks the same either way — coefficients, standard errors, a healthy R². Diagnostics are how you find out whether one point is steering the whole fit, whether the predictors are too entangled to separate, and whether the model will predict new data or has just memorized the old. Skipping them is how confident, wrong models get shipped.

Residuals are the evidence. If the model is right, the residuals look like structureless noise. A curve in a residual plot means a missing nonlinear term; a funnel shape means the spread changes with x; a few huge residuals mean outliers.

Leverage measures how unusual a point's predictor values are. It's the i-th diagonal entry of the hat matrix H = X(XᵀX)⁻¹Xᵀ (the matrix that turns y into ŷ). Leverages lie between 1/n and 1 and always sum to p, the number of parameters, so the average is p/n; above 2p/n or so is high. A high-leverage point pulls the line toward itself.

Standardized residuals. A raw residual is hard to judge, and a high-leverage point's residual is artificially small because the line bends toward it. Dividing each residual by its own standard error, s·√(1 − hᵢ), puts them on one scale: beyond ±2 is unusual, beyond ±3 is suspicious.

Cook's distance combines the two: how much would all the fitted values move if this point were left out? A point needs both leverage and a large residual to be influential. Values near 1 (or far above the rest) deserve investigation.

Multicollinearity, measured. The variance inflation factor for predictor j is 1/(1 − R²ⱼ), where R²ⱼ comes from regressing that predictor on all the others. VIF = 1 means it's unrelated to the rest; a VIF of 10 means its coefficient's variance is ten times what it would be otherwise, so its standard error is about three times larger.

In-sample fit flatters. Adding terms always lowers the error on the data you fitted. The honest test is data the model hasn't seen. k-fold cross-validation splits the data into k parts, fits on k − 1 of them, and measures error on the part held out — k times, so every point is predicted once without having been used to fit.

Overfitting. A flexible enough model (a high-degree polynomial, say) can chase the noise in the training data. Its training error keeps falling while its cross-validated error turns back up. The best model is the one with the lowest held-out error.

• Sum first, divide once. Accumulating a mean as m += x[i] / n adds n rounding errors: ten copies of 8 average to 7.999999999999999. In the worked example that turns "no spread at all" into a tiny nonzero spread and a meaningless slope. Add everything up, then divide.
• Float division by zero doesn't panic. 0.0/0.0 is NaN and 1.0/0.0 is +Inf. That's handy — a NaN slope is an honest answer — but check with math.IsNaN before a NaN quietly spreads through later arithmetic.
• Don't build what you only need the diagonal of. The hat matrix is n×n; leverage needs only n numbers, each xᵢᵀ·(XᵀX)⁻¹·xᵢ. For 100,000 rows that's the difference between 80 GB and 800 KB.
• Building training folds: collect the rows where i%k != f with append, then reuse the same fitting code. Interleaved folds (i % k) spread each fold across the whole range of x, which matters when the data is sorted.
• Given helpers are a seam. This module's fitOLS and xtxInverse live in types.go, wrapping internal/linalg. Swap in your Module 8 code if you like — the tests don't care.

package main

import "fmt"

// line fits y = a + b·x and also returns each point's leverage,
// h = 1/n + (x − x̄)²/Sxx: how far its x sits from the others.
func line(x, y []float64) (a, b float64, h []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 // sum first, divide once: x̄ is exactly 8, not 7.999999999999999
	var sxx, sxy float64
	for i := range x {
		sxx += (x[i] - mx) * (x[i] - mx)
		sxy += (x[i] - mx) * (y[i] - my)
	}
	b = sxy / sxx // 0/0 is NaN in Go's float64 — no panic
	a = my - b*mx
	for i := range x {
		h = append(h, 1/n+(x[i]-mx)*(x[i]-mx)/sxx)
	}
	return a, b, h
}

func main() {
	// Anscombe's set 4: ten points at x = 8, and one at x = 19.
	x := []float64{8, 8, 8, 8, 8, 8, 8, 19, 8, 8, 8}
	y := []float64{6.58, 5.76, 7.71, 8.84, 8.47, 7.04, 5.25, 12.50, 5.56, 7.91, 6.89}

	a, b, h := line(x, y)
	fmt.Printf("all 11 points:      y = %.2f + %.3f·x\n", a, b)
	fmt.Print("leverage:          ")
	for _, v := range h {
		fmt.Printf(" %.2f", v)
	}
	fmt.Println()

	_, b, _ = line(append(x[:7:7], x[8:]...), append(y[:7:7], y[8:]...))
	fmt.Printf("without x = 19:     slope = %v (every x is 8: Sxy/Sxx = 0/0)\n", b)
}

It prints:

Ten of the eleven points have leverage 0.10; the one at x = 19 has leverage 1.00, the maximum possible, so the fitted line must pass exactly through it. Remove it and there's no slope to estimate at all. That's Anscombe's set 4: the same r and line as the other three sets, decided entirely by a single observation. A residual plot can't flag it either — its residual is exactly zero. Leverage can.

• Deleting influential points to improve the fit. Investigate them. A data-entry error should be fixed; a genuine extreme case may be the most informative observation you have, and removing it overstates your certainty.
• Checking only R² and p-values. They can't see curvature, outliers, or leverage. Plot the residuals against the fitted values every time.
• Judging predictive power on the training data. The training error of a model that has seen the data is optimistic, and more so the more flexible the model. Cross-validate, or hold out a test set and touch it once.
• Dropping a predictor just because its VIF is high. If you need its coefficient, collinearity is a problem; if you only need predictions, it often isn't. Know which question you're answering.
• Leaking information across folds. Anything learned from data — scaling factors, chosen features — has to be learned inside each training fold, or the held-out error is optimistic too.

Part A finds influential points: leverage, standardized residuals, and Cook's distance. Part B measures multicollinearity with variance inflation factors. Part C measures prediction error — RMSE, MAE — and uses k-fold cross-validation to choose the degree of a polynomial.

go test ./09-regression-diagnostics/

The expected values come from statsmodels' OLSInfluence and variance_inflation_factor, run on Module 8's house data. TestKFoldRMSE fits polynomials of degree 1 to 6 to a noisy quadratic: training error falls every time, and cross-validation should still pick degree 2.

The files for this module

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

// Package diagnostics is Module 9: checking a regression before you trust it.
//
// A fitted model always produces numbers. This module asks whether they
// deserve belief: which points are steering the fit, whether predictors are
// too entangled to separate, and how well the model predicts data it hasn't
// seen.
//
// Design matrices here are plain [][]float64 with a leading column of 1s,
// as in Module 8. The fitting helpers are given — they're the reference
// version of what you built there — so this module doesn't depend on yours.
package diagnostics

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

// designMatrix returns the rows [1, cols[0][i], cols[1][i], …].
func designMatrix(cols ...[]float64) [][]float64 {
	x := make([][]float64, len(cols[0]))
	for i := range x {
		x[i] = make([]float64, len(cols)+1)
		x[i][0] = 1
		for j, c := range cols {
			x[i][j+1] = c[i]
		}
	}
	return x
}

// xtxInverse returns (XᵀX)⁻¹. It panics if a predictor is an exact
// combination of the others.
func xtxInverse(x [][]float64) [][]float64 {
	inv, err := linalg.Inverse(linalg.Mul(linalg.Transpose(x), x))
	if err != nil {
		panic(err)
	}
	return inv
}

// fitOLS returns the least-squares coefficients b = (XᵀX)⁻¹Xᵀy and the
// residuals y − Xb.
func fitOLS(x [][]float64, y []float64) (coef, resid []float64) {
	coef = linalg.MulVec(xtxInverse(x), linalg.MulVec(linalg.Transpose(x), y))
	fitted := linalg.MulVec(x, coef)
	resid = make([]float64, len(y))
	for i := range y {
		resid[i] = y[i] - fitted[i]
	}
	return coef, resid
}

// rSquared returns R² for the least-squares fit of y on x.
func rSquared(x [][]float64, y []float64) float64 {
	_, resid := fitOLS(x, y)
	mean := 0.0
	for _, v := range y {
		mean += v / float64(len(y))
	}
	var sse, sst float64
	for i := range y {
		sse += resid[i] * resid[i]
		sst += (y[i] - mean) * (y[i] - mean)
	}
	return 1 - sse/sst
}

//go:build !solutions

package diagnostics

// ============================================================
// Part A — influential points
// ============================================================

// Leverage returns each observation's leverage h_i = x_iᵀ(XᵀX)⁻¹x_i, where
// x_i is row i of X: the diagonal of the "hat matrix". It measures how
// unusual a point's predictor values are, between 1/n and 1. Use the given
// xtxInverse; don't build the whole n×n hat matrix.
func Leverage(x [][]float64) []float64 {
	// TODO
	return nil
}

// StandardizedResiduals returns each residual divided by its own standard
// error, r_i = e_i / (s·√(1 − h_i)), with s = √(SSE/(n − p)). High-leverage
// points pull the line toward themselves, so their raw residuals look
// smaller than they should; this corrects for that. |r| > 2 deserves a look.
func StandardizedResiduals(x [][]float64, y []float64) []float64 {
	// TODO
	return nil
}

// CooksDistance returns Cook's distance for each observation,
// D_i = r_i²·h_i / (p·(1 − h_i)), with r_i the standardized residual and p
// the number of columns of X. It combines "unusual x" (leverage) with
// "badly fitted" (residual) into "how much the fit would move without this
// point". Values near 1 are a red flag.
func CooksDistance(x [][]float64, y []float64) []float64 {
	// TODO
	return nil
}

// ============================================================
// Part B — entangled predictors
// ============================================================

// VIF returns the variance inflation factor of each predictor column:
// regress column j on all the OTHER columns (with an intercept), take that
// R²_j, and return 1 / (1 − R²_j). A VIF of 5 means the coefficient's
// variance is 5 times what it would be with uncorrelated predictors; above
// about 5–10 is trouble. Use the given designMatrix and rSquared.
func VIF(cols [][]float64) []float64 {
	// TODO
	return nil
}

// ============================================================
// Part C — prediction error and overfitting
// ============================================================

// RMSE returns the root mean squared error √(Σ(a − p)²/n): the typical
// prediction error, in y's units, with big misses weighted heavily.
func RMSE(actual, predicted []float64) float64 {
	// TODO
	return 0
}

// MAE returns the mean absolute error Σ|a − p|/n: the typical miss, less
// sensitive to a few large ones than RMSE.
func MAE(actual, predicted []float64) float64 {
	// TODO
	return 0
}

// PolyFeatures returns the columns x, x², …, x^degree, so a polynomial can
// be fitted as a multiple regression: designMatrix(PolyFeatures(x, 3)...).
func PolyFeatures(x []float64, degree int) [][]float64 {
	// TODO
	return nil
}

// KFoldRMSE estimates how well a linear model on cols predicts NEW data.
// Put observation i in fold i % k. For each fold: fit on the other folds
// (use the given designMatrix and fitOLS on just those rows), predict the
// held-out rows, and collect their errors. Return the RMSE over all n
// held-out predictions.
func KFoldRMSE(cols [][]float64, y []float64, k int) float64 {
	// TODO
	return 0
}

package diagnostics

import (
	"math"
	"testing"
)

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

func slicesNear(t *testing.T, name string, got, want []float64, tol float64) {
	t.Helper()
	if len(got) != len(want) {
		t.Fatalf("%s has %d values, want %d", name, len(got), len(want))
	}
	for i := range want {
		if math.Abs(got[i]-want[i]) > tol*math.Max(1, math.Abs(want[i])) {
			t.Errorf("%s[%d] = %v, want %v", name, i, got[i], want[i])
		}
	}
}

// the 15 made-up house sales from Module 8
var (
	size  = []float64{14.2, 16.8, 12.5, 19.6, 22.1, 15.4, 18.3, 11.9, 24.6, 17.2, 13.8, 20.5, 16.1, 21.7, 14.9}
	beds  = []float64{3, 3, 2, 4, 4, 3, 3, 2, 5, 3, 2, 4, 3, 4, 3}
	age   = []float64{25, 12, 40, 8, 5, 30, 15, 45, 3, 20, 35, 10, 22, 6, 28}
	price = []float64{220.4, 252.2, 172.0, 304.6, 352.2, 215.8, 277.6, 149.8, 389.6, 246.4, 192.6, 323.0, 236.8, 329.2, 224.2}
)

func TestLeverage(t *testing.T) {
	want := []float64{0.40412436677955377, 0.4382095760859276, 0.2453550135039821, 0.23334656462979964,
		0.24037885374411389, 0.18286061160829314, 0.26912918063697644, 0.3703888411947982,
		0.5794689014974526, 0.10064193948796983, 0.32593802598558574, 0.13024956468548918,
		0.09384674588656905, 0.18812331917416053, 0.19793849509932893}
	h := Leverage(designMatrix(size, beds, age))
	slicesNear(t, "Leverage", h, want, 1e-9)
	sum := 0.0
	for _, v := range h {
		sum += v
	}
	if !near(sum, 4) { // leverages always sum to p, the number of parameters
		t.Errorf("leverages sum to %v, want 4", sum)
	}
	// Anscombe's set 4: ten points at x = 8, one at x = 19
	h4 := Leverage(designMatrix([]float64{8, 8, 8, 8, 8, 8, 8, 19, 8, 8, 8}))
	if len(h4) != 11 || !near(h4[7], 1) || !near(h4[0], 0.1) {
		t.Errorf("Anscombe 4 leverages = %v, want 0.1 everywhere except 1 at index 7", h4)
	}
}

func TestStandardizedResiduals(t *testing.T) {
	want := []float64{1.483071512955371, -0.857737399574609, 0.8378562426912499, -0.5393593604393221,
		1.4929103601170326, -1.3137407991866639, 0.5979285100385263, -0.945979090806463,
		0.1866600903244587, -1.3870779133050453, 0.8183863917519304, 0.7426586954519439,
		-0.4244877052437035, -1.236866731649812, 0.8003068877277463}
	slicesNear(t, "StandardizedResiduals", StandardizedResiduals(designMatrix(size, beds, age), price), want, 1e-8)
}

func TestCooksDistance(t *testing.T) {
	want := []float64{0.3729268092118823, 0.14346839312507037, 0.05705993516610014, 0.022135981057478413,
		0.17632207423403143, 0.09655692051624647, 0.03291233527775299, 0.13160990044360404,
		0.012002589908987702, 0.05382550184786958, 0.08096401007347996, 0.020649054718565932,
		0.004665388391051025, 0.08862135591142416, 0.039516248617718235}
	slicesNear(t, "CooksDistance", CooksDistance(designMatrix(size, beds, age), price), want, 1e-8)
}

func TestVIF(t *testing.T) {
	// size, bedrooms and age are all entangled: every VIF is high
	slicesNear(t, "VIF", VIF([][]float64{size, beds, age}),
		[]float64{15.453470348060517, 10.217249396975733, 8.705220683741848}, 1e-8)
	// unrelated predictors: VIF 1
	v := VIF([][]float64{{1, 2, 3, 4}, {1, -1, -1, 1}})
	if len(v) != 2 || !near(v[0], 1) || !near(v[1], 1) {
		t.Errorf("VIF of uncorrelated columns = %v, want [1 1]", v)
	}
}

func TestErrors(t *testing.T) {
	a := []float64{3, 5.5, 2, 8}
	p := []float64{2.5, 6, 2, 5}
	if got := RMSE(a, p); !near(got, 1.541103500742244) {
		t.Errorf("RMSE = %v, want 1.541103500742244", got)
	}
	if got := MAE(a, p); !near(got, 1) {
		t.Errorf("MAE = %v, want 1", got)
	}
}

func TestPolyFeatures(t *testing.T) {
	cols := PolyFeatures([]float64{2, -1, 0.5}, 3)
	want := [][]float64{{2, -1, 0.5}, {4, 1, 0.25}, {8, -1, 0.125}}
	if len(cols) != 3 {
		t.Fatalf("PolyFeatures(…, 3) returned %d columns, want 3", len(cols))
	}
	for d := range want {
		slicesNear(t, "PolyFeatures column", cols[d], want[d], 1e-12)
	}
}

// curve is y = 1 + 2x − 3x² plus noise, at 20 points from −1 to 1.
func curve() (x, y []float64) {
	y = []float64{-3.79, -3.541, -2.349, -1.353, -1.343, -0.89, 0.186, 0.316, 0.199, 1.046,
		1.387, 1.121, 0.939, 1.57, 1.344, 0.932, 1.324, 0.629, 0.578, -0.31}
	x = make([]float64, 20)
	for i := range x {
		x[i] = (float64(i) - 9.5) / 9.5
	}
	return x, y
}

func TestKFoldRMSE(t *testing.T) {
	x, y := curve()
	if len(PolyFeatures(x, 1)) != 1 {
		t.Fatal("PolyFeatures(x, 1) should return one column — write PolyFeatures first")
	}
	want := []float64{1.1094386854745353, 0.2999916032680977, 0.3230957759122377,
		0.4121703436271352, 0.42893131710023574, 0.518756851228933}
	best, bestDeg := math.Inf(1), 0
	for deg := 1; deg <= 6; deg++ {
		got := KFoldRMSE(PolyFeatures(x, deg), y, 5)
		if math.Abs(got-want[deg-1]) > 1e-7 {
			t.Errorf("KFoldRMSE(degree %d) = %v, want %v", deg, got, want[deg-1])
		}
		if got < best {
			best, bestDeg = got, deg
		}
	}
	if bestDeg != 2 {
		t.Errorf("cross-validation picked degree %d; the data is a quadratic", bestDeg)
	}
}

//go:build solutions

package diagnostics

import "math"

// ============================================================
// Part A — influential points
// ============================================================

func Leverage(x [][]float64) []float64 {
	inv := xtxInverse(x)
	h := make([]float64, len(x))
	for i, row := range x {
		for a := range row {
			for b := range row {
				h[i] += row[a] * inv[a][b] * row[b]
			}
		}
	}
	return h
}

func StandardizedResiduals(x [][]float64, y []float64) []float64 {
	_, e := fitOLS(x, y)
	h := Leverage(x)
	n, p := len(x), len(x[0])
	sse := 0.0
	for _, v := range e {
		sse += v * v
	}
	s := math.Sqrt(sse / float64(n-p))
	r := make([]float64, n)
	for i := range e {
		r[i] = e[i] / (s * math.Sqrt(1-h[i]))
	}
	return r
}

func CooksDistance(x [][]float64, y []float64) []float64 {
	r := StandardizedResiduals(x, y)
	h := Leverage(x)
	p := float64(len(x[0]))
	d := make([]float64, len(r))
	for i := range r {
		d[i] = r[i] * r[i] * h[i] / (p * (1 - h[i]))
	}
	return d
}

// ============================================================
// Part B — entangled predictors
// ============================================================

func VIF(cols [][]float64) []float64 {
	out := make([]float64, len(cols))
	for j := range cols {
		var others [][]float64
		for k := range cols {
			if k != j {
				others = append(others, cols[k])
			}
		}
		out[j] = 1 / (1 - rSquared(designMatrix(others...), cols[j]))
	}
	return out
}

// ============================================================
// Part C — prediction error and overfitting
// ============================================================

func RMSE(actual, predicted []float64) float64 {
	s := 0.0
	for i := range actual {
		d := actual[i] - predicted[i]
		s += d * d
	}
	return math.Sqrt(s / float64(len(actual)))
}

func MAE(actual, predicted []float64) float64 {
	s := 0.0
	for i := range actual {
		s += math.Abs(actual[i] - predicted[i])
	}
	return s / float64(len(actual))
}

func PolyFeatures(x []float64, degree int) [][]float64 {
	cols := make([][]float64, degree)
	for d := range cols {
		cols[d] = make([]float64, len(x))
		for i, v := range x {
			cols[d][i] = math.Pow(v, float64(d+1))
		}
	}
	return cols
}

func KFoldRMSE(cols [][]float64, y []float64, k int) float64 {
	n := len(y)
	x := designMatrix(cols...)
	sq := 0.0
	for f := 0; f < k; f++ {
		var trainX [][]float64
		var trainY []float64
		for i := 0; i < n; i++ {
			if i%k != f {
				trainX = append(trainX, x[i])
				trainY = append(trainY, y[i])
			}
		}
		coef, _ := fitOLS(trainX, trainY)
		for i := f; i < n; i += k {
			pred := 0.0
			for j, v := range x[i] {
				pred += coef[j] * v
			}
			sq += (y[i] - pred) * (y[i] - pred)
		}
	}
	return math.Sqrt(sq / float64(n))
}

• Add externally studentized residuals, which use an s computed without point i, so an outlier can't hide by inflating s.
• Plot residuals against fitted values as text (# marks in a grid) and look for curves and funnels.
• Practice: take any regression you've run, compute its Cook's distances, and refit without the most influential point. How much did the coefficients move?

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 10, so git log shows your own progress through the course.