Simple Regression

Regression is the workhorse of applied statistics. "How much do sales rise per $1,000 of advertising?" "How much does each extra year of experience add to salary?" "What will this house sell for?" A fitted line answers all of them with a number — the slope — and, done properly, with how uncertain that number is.

Simple regression has one predictor. Everything in it carries straight over to Module 8's multiple regression, where the same ideas are written with matrices.

Least squares. For each point, the residual is the vertical gap between the actual y and the line's prediction ŷ. Ordinary least squares picks the intercept a and slope b that make the sum of squared residuals (SSE) as small as possible. Squaring treats over- and under-predictions alike and leads to a closed-form answer: b = Sxy/Sxx and a = ȳ − b·x̄. The line always passes through the point (x̄, ȳ), and its residuals always sum to zero.

Reading the coefficients. The slope is "the change in predicted y for a one-unit increase in x". The intercept is the prediction at x = 0, which only means something if x = 0 is sensible and inside your data.

R². SST is the total squared variation in y around its mean; SSE is what's left after the line. R² = 1 − SSE/SST is the fraction the line explains. With one predictor, R² is exactly r² from Module 6.

The residual standard error s = √(SSE/(n − 2)) is the typical size of a miss, in y's units. It divides by n − 2 because the line used up two parameters.

Inference on the slope. The slope from your sample is an estimate; its standard error is s/√Sxx. Spreading the x values out makes Sxx larger and the slope more certain. The t-test of "slope = 0" asks whether x is related to y at all, and the confidence interval gives the plausible range.

Two kinds of prediction interval. A confidence interval for the mean response at x₀ says where the line itself plausibly lies. A prediction interval for one new observation adds that observation's own scatter, so it's always wider and never shrinks to zero however much data you have. Both widen as x₀ moves away from x̄ — extrapolation is uncertain.

• A struct that carries its own sums. FitLine computes Sxx, x̄, SSE, and SST once and stores them in Fit, so every method afterwards — R², standard errors, intervals — is a one-line formula with no second pass over the data.
• Value receivers (func (f Fit) R2() float64) suit a small, read-only result. Nothing a method does can change the fit.
• Methods compose. SlopeCI calls SlopeSE, which calls ResidualSE. Each formula appears once, which is how the code matches the table above line for line.
• Two passes, not one. You need x̄ before you can form x − x̄. Folding everything into one pass with Σx², Σxy and so on is the shortcut formula from Module 1, with the same cancellation problem.

package main

import "fmt"

func main() {
	// Advertising spend ($k) and sales ($k) for five months (made-up numbers).
	x := []float64{1, 2, 3, 4, 5}
	y := []float64{3, 7, 9, 10, 16}

	var mx, my float64
	for i := range x {
		mx += x[i]
		my += y[i]
	}
	mx, my = mx/5, my/5
	fmt.Printf("x̄ = %.1f   ȳ = %.1f\n\n", mx, my)
	fmt.Println("   x     y   x−x̄   y−ȳ   (x−x̄)(y−ȳ)   (x−x̄)²")

	var sxy, sxx, sst float64
	for i := range x {
		dx, dy := x[i]-mx, y[i]-my
		fmt.Printf("%4.0f  %4.0f  %4.0f  %4.0f  %11.0f  %7.0f\n", x[i], y[i], dx, dy, dx*dy, dx*dx)
		sxy += dx * dy
		sxx += dx * dx
		sst += dy * dy
	}

	b := sxy / sxx // slope
	a := my - b*mx // intercept: the line goes through (x̄, ȳ)
	var sse float64
	for i := range x {
		e := y[i] - (a + b*x[i])
		sse += e * e
	}
	fmt.Printf("\nSxy = %.0f   Sxx = %.0f   slope b = Sxy/Sxx = %.1f   intercept a = ȳ − b·x̄ = %.1f\n", sxy, sxx, b, a)
	fmt.Printf("ŷ = %.1f + %.1f·x     SSE = %.1f   SST = %.1f   R² = 1 − SSE/SST = %.3f\n", a, b, sse, sst, 1-sse/sst)
}

It prints:

Every number in a regression output comes from these few sums. The slope says each extra $1k of advertising goes with about $2.9k more in sales; the line explains 93% of the month-to-month variation. Whether that slope is trustworthy from five months of data is what the standard error and t-test are for — and with n = 5 it's worth checking.

• Extrapolating. The line is only supported where you have data. Predicting sales at $50k of advertising from data between $1k and $5k assumes the straight line keeps going; it often doesn't.
• Reading the slope as causal. Regression measures association, just as correlation does. Advertising and sales might both follow the season.
• Judging a fit by R² alone. Anscombe's four data sets all have R² = 0.67. Plot the residuals: they should look like a structureless band around zero.
• Using a confidence interval for the mean when you need a prediction interval. "Where will this house's price be?" needs the wider prediction interval.
• Regressing y on x and then solving for x. The regression of x on y is a different line. If you need to predict x, fit that line.

Part A fits the line, predicts, and computes residuals. Part B measures the fit with R² and the residual standard error. Part C is inference: standard errors, the slope test, its confidence interval, and both intervals for a prediction.

go test ./07-simple-regression/

The expected values come from statsmodels' OLS: its coefficients, bse, tvalues, pvalues, conf_int, and get_prediction(...).summary_frame. Store everything you'll need in Fit and the methods in Parts B and C are short.

The files for this module

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

// Package simplereg is Module 7: simple linear regression.
//
// One predictor x, one response y, and the straight line ŷ = a + b·x that
// makes the squared vertical distances from the points as small as possible
// (ordinary least squares). Then: how good is the line, how sure are we of
// its slope, and how far off will a prediction be?
package simplereg

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

// Fit is a fitted line plus the sums every later calculation needs.
// FitLine fills in every field.
type Fit struct {
	Intercept float64 // a: the predicted y when x = 0
	Slope     float64 // b: the change in predicted y per unit of x
	N         int     // number of points
	MeanX     float64 // x̄
	Sxx       float64 // Σ(x − x̄)²
	SSE       float64 // Σ(y − ŷ)²: the squared error the line leaves
	SST       float64 // Σ(y − ȳ)²: the squared error with no line at all
}

// Interval is a range [Lo, Hi].
type Interval struct {
	Lo, Hi float64
}

// The helpers below are given — you built them in Module 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) }

// tQuantile returns the t with P(T <= t) = p.
func tQuantile(p, df float64) float64 { return special.TQuantile(p, df) }

//go:build !solutions

package simplereg

// ============================================================
// Part A — fitting the line
// ============================================================

// FitLine fits ŷ = a + b·x by least squares and fills in every field of
// Fit:
//
//	Sxx = Σ(x − x̄)²    Sxy = Σ(x − x̄)(y − ȳ)
//	b = Sxy / Sxx      a = ȳ − b·x̄
//	SSE = Σ(y − ŷ)²    SST = Σ(y − ȳ)²
//
// Target: O(n), two passes.
func FitLine(x, y []float64) Fit {
	// TODO
	return Fit{}
}

// Predict returns the fitted value a + b·x.
func (f Fit) Predict(x float64) float64 {
	// TODO
	return 0
}

// Residuals returns y[i] − ŷ[i] for each point: what the line gets wrong.
// For a least-squares line they always sum to zero.
func (f Fit) Residuals(x, y []float64) []float64 {
	// TODO
	return nil
}

// ============================================================
// Part B — how good is the line?
// ============================================================

// R2 returns the coefficient of determination, 1 − SSE/SST: the fraction
// of y's variation the line explains. For one predictor it equals r².
func (f Fit) R2() float64 {
	// TODO
	return 0
}

// ResidualSE returns s = √(SSE / (n − 2)), the typical size of a residual —
// "the line is usually off by about s". n − 2 because two parameters were
// estimated.
func (f Fit) ResidualSE() float64 {
	// TODO
	return 0
}

// ============================================================
// Part C — inference: how sure are we?
// ============================================================

// SlopeSE returns the standard error of the slope, s / √Sxx. Spreading the
// x values out (a bigger Sxx) pins the slope down better.
func (f Fit) SlopeSE() float64 {
	// TODO
	return 0
}

// InterceptSE returns the standard error of the intercept,
// s·√(1/n + x̄²/Sxx).
func (f Fit) InterceptSE() float64 {
	// TODO
	return 0
}

// SlopeTest tests H₀: slope = 0 ("x tells us nothing about y"):
// t = b / SE(b), df = n − 2, two-sided p = 2·tCDF(−|t|, df).
func (f Fit) SlopeTest() (t, p float64) {
	// TODO
	return 0, 0
}

// SlopeCI returns b ± t*·SE(b) with t* = tQuantile((1 + level)/2, n − 2).
func (f Fit) SlopeCI(level float64) Interval {
	// TODO
	return Interval{}
}

// MeanCI returns a confidence interval for the AVERAGE y at x0:
// ŷ ± t*·s·√(1/n + (x0 − x̄)²/Sxx). It's narrowest at x̄.
func (f Fit) MeanCI(x0, level float64) Interval {
	// TODO
	return Interval{}
}

// PredictionInterval returns an interval for ONE new observation at x0:
// ŷ ± t*·s·√(1 + 1/n + (x0 − x̄)²/Sxx). The extra 1 is the new point's own
// scatter, so it's always wider than MeanCI.
func (f Fit) PredictionInterval(x0, level float64) Interval {
	// TODO
	return Interval{}
}

package simplereg

import (
	"math"
	"testing"
)

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

// nearish allows for TQuantile's bisection (about 12 digits): 1e-8.
func nearish(a, b float64) bool {
	return math.Abs(a-b) <= 1e-8*math.Max(1, math.Abs(b))
}

// hours studied and exam score for 12 students (made-up data)
var (
	hours = []float64{1, 2, 2, 3, 4, 5, 5, 6, 7, 8, 9, 10}
	score = []float64{52, 55, 61, 58, 66, 70, 68, 75, 77, 84, 85, 91}
)

func TestFitLine(t *testing.T) {
	f := FitLine(hours, score)
	if !near(f.Intercept, 48.23131672597867) || !near(f.Slope, 4.245551601423493) {
		t.Errorf("line = %v + %v·x, want 48.23131672597867 + 4.245551601423493·x", f.Intercept, f.Slope)
	}
	if f.N != 12 || !near(f.MeanX, 5.166666666666667) || !near(f.Sxx, 93.66666666666666) {
		t.Errorf("N, MeanX, Sxx = %d, %v, %v; want 12, 5.1667, 93.6667", f.N, f.MeanX, f.Sxx)
	}
	if !near(f.SSE, 41.3523131672598) || !near(f.SST, 1729.6666666666665) {
		t.Errorf("SSE, SST = %v, %v; want 41.3523131672598, 1729.6666666666665", f.SSE, f.SST)
	}
	// a perfect line is recovered exactly
	g := FitLine([]float64{0, 1, 2, 3}, []float64{3, 5, 7, 9})
	if !near(g.Intercept, 3) || !near(g.Slope, 2) || !near(g.SSE, 0) {
		t.Errorf("FitLine on y = 3 + 2x gave %+v", g)
	}
}

func TestPredict(t *testing.T) {
	f := FitLine(hours, score)
	if got := f.Predict(5.5); !near(got, 71.58185053380788) {
		t.Errorf("Predict(5.5) = %v, want 71.58185053380788", got)
	}
	if got := f.Predict(f.MeanX); !near(got, 842.0/12) { // the line passes through (x̄, ȳ)
		t.Errorf("Predict(x̄) = %v, want ȳ = 842/12 = 70.1667", got)
	}
}

func TestResiduals(t *testing.T) {
	want := []float64{-0.47686832740216545, -1.7224199288256585, 4.2775800711743415, -2.9679715302491516,
		0.7864768683273553, 0.5409252669038551, -1.4590747330961449, 1.2953736654803691,
		-0.9501779359431168, 1.804270462633383, -1.4412811387901172, 0.3131672597863968}
	got := FitLine(hours, score).Residuals(hours, score)
	if len(got) != len(want) {
		t.Fatalf("Residuals returned %d values, want %d", len(got), len(want))
	}
	sum := 0.0
	for i := range want {
		if !near(got[i], want[i]) {
			t.Errorf("residual %d = %v, want %v", i, got[i], want[i])
		}
		sum += got[i]
	}
	if math.Abs(sum) > 1e-9 {
		t.Errorf("residuals sum to %v, want 0", sum)
	}
}

func TestR2(t *testing.T) {
	if got := FitLine(hours, score).R2(); !near(got, 0.9760923223160957) {
		t.Errorf("R2 = %v, want 0.9760923223160957", got)
	}
}

func TestResidualSE(t *testing.T) {
	if got := FitLine(hours, score).ResidualSE(); !near(got, 2.0335268173117313) {
		t.Errorf("ResidualSE = %v, want 2.0335268173117313", got)
	}
}

func TestStandardErrors(t *testing.T) {
	f := FitLine(hours, score)
	if got := f.SlopeSE(); !near(got, 0.21011515187367277) {
		t.Errorf("SlopeSE = %v, want 0.21011515187367277", got)
	}
	if got := f.InterceptSE(); !near(got, 1.234147077182429) {
		t.Errorf("InterceptSE = %v, want 1.234147077182429", got)
	}
}

func TestSlopeTest(t *testing.T) {
	tStat, p := FitLine(hours, score).SlopeTest()
	if !near(tStat, 20.20583267586547) || math.Abs(p-1.9415982727901893e-09) > 1e-9*1.9415982727901893e-09 {
		t.Errorf("SlopeTest = (%v, %v), want (20.20583267586547, 1.9415982727901893e-09)", tStat, p)
	}
	// no relationship at all: a flat cloud
	_, p = FitLine([]float64{1, 2, 3, 4, 5, 6}, []float64{4, 6, 5, 6, 4, 5}).SlopeTest()
	if p < 0.5 {
		t.Errorf("SlopeTest on a flat cloud: p = %v, want a large p-value", p)
	}
}

func TestSlopeCI(t *testing.T) {
	f := FitLine(hours, score)
	for _, tc := range []struct {
		level  float64
		lo, hi float64
	}{
		{0.95, 3.7773858681427663, 4.71371733470422},
		{0.90, 3.8647260573387903, 4.6263771455081955},
	} {
		if got := f.SlopeCI(tc.level); !nearish(got.Lo, tc.lo) || !nearish(got.Hi, tc.hi) {
			t.Errorf("SlopeCI(%v) = %+v, want {%v %v}", tc.level, got, tc.lo, tc.hi)
		}
	}
}

func TestIntervalsAtX0(t *testing.T) {
	f := FitLine(hours, score)
	for _, tc := range []struct {
		x0             float64
		meanLo, meanHi float64
		predLo, predHi float64
	}{
		{5.5, 70.2645926619159, 72.89910840569986, 66.86327563856246, 76.3004254290533},
		{12, 95.72174383568341, 102.63412805043777, 93.4792479691428, 104.87662391697837}, // extrapolating: wider
	} {
		m, p := f.MeanCI(tc.x0, 0.95), f.PredictionInterval(tc.x0, 0.95)
		if !nearish(m.Lo, tc.meanLo) || !nearish(m.Hi, tc.meanHi) {
			t.Errorf("MeanCI(%v) = %+v, want {%v %v}", tc.x0, m, tc.meanLo, tc.meanHi)
		}
		if !nearish(p.Lo, tc.predLo) || !nearish(p.Hi, tc.predHi) {
			t.Errorf("PredictionInterval(%v) = %+v, want {%v %v}", tc.x0, p, tc.predLo, tc.predHi)
		}
	}
}

//go:build solutions

package simplereg

import "math"

// ============================================================
// Part A — fitting the line
// ============================================================

func FitLine(x, y []float64) Fit {
	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 sxx, sxy, sst float64
	for i := range x {
		sxx += (x[i] - mx) * (x[i] - mx)
		sxy += (x[i] - mx) * (y[i] - my)
		sst += (y[i] - my) * (y[i] - my)
	}
	f := Fit{Slope: sxy / sxx, N: len(x), MeanX: mx, Sxx: sxx, SST: sst}
	f.Intercept = my - f.Slope*mx
	for i := range x {
		e := y[i] - f.Predict(x[i])
		f.SSE += e * e
	}
	return f
}

func (f Fit) Predict(x float64) float64 {
	return f.Intercept + f.Slope*x
}

func (f Fit) Residuals(x, y []float64) []float64 {
	out := make([]float64, len(x))
	for i := range x {
		out[i] = y[i] - f.Predict(x[i])
	}
	return out
}

// ============================================================
// Part B — how good is the line?
// ============================================================

func (f Fit) R2() float64 {
	return 1 - f.SSE/f.SST
}

func (f Fit) ResidualSE() float64 {
	return math.Sqrt(f.SSE / float64(f.N-2))
}

// ============================================================
// Part C — inference: how sure are we?
// ============================================================

func (f Fit) SlopeSE() float64 {
	return f.ResidualSE() / math.Sqrt(f.Sxx)
}

func (f Fit) InterceptSE() float64 {
	return f.ResidualSE() * math.Sqrt(1/float64(f.N)+f.MeanX*f.MeanX/f.Sxx)
}

func (f Fit) SlopeTest() (t, p float64) {
	df := float64(f.N - 2)
	t = f.Slope / f.SlopeSE()
	return t, 2 * tCDF(-math.Abs(t), df)
}

func (f Fit) SlopeCI(level float64) Interval {
	half := tQuantile((1+level)/2, float64(f.N-2)) * f.SlopeSE()
	return Interval{f.Slope - half, f.Slope + half}
}

func (f Fit) MeanCI(x0, level float64) Interval {
	se := f.ResidualSE() * math.Sqrt(1/float64(f.N)+(x0-f.MeanX)*(x0-f.MeanX)/f.Sxx)
	half := tQuantile((1+level)/2, float64(f.N-2)) * se
	y := f.Predict(x0)
	return Interval{y - half, y + half}
}

func (f Fit) PredictionInterval(x0, level float64) Interval {
	se := f.ResidualSE() * math.Sqrt(1+1/float64(f.N)+(x0-f.MeanX)*(x0-f.MeanX)/f.Sxx)
	half := tQuantile((1+level)/2, float64(f.N-2)) * se
	y := f.Predict(x0)
	return Interval{y - half, y + half}
}

• Add Leverage(x0) = 1/n + (x₀ − x̄)²/Sxx and notice it's the quantity under the square roots above. Module 9 generalizes it.
• Fit a line through the origin (no intercept): b = Σxy / Σx². When is that justified?
• Practice: regress your monthly spending on your monthly income, and work out the slope's 95% confidence interval.

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