Linear Regression From Scratch in NumPy: Assumptions and Residual Analysis

Your linear regression ran without errors. It printed a slope, an intercept, and a nice high score. So is the model right? Not necessarily. A line can fit and still be misleading, and the evidence is hiding in the leftovers, the residuals.

• The code ran
• The score looks high
• Is it trustworthy?

This is lesson forty. In lesson fourteen we derived linear regression and saw the normal equation and gradient descent. Today we build it as real code, test it, and then interrogate it with residual analysis. We will not re-derive the equation, and regularization waits for lesson seventeen.

• Build the class
• Test it
• Diagnose it

Here is the plan. First, a class with two fit methods, every line tied to a symbol. Second, tests against scikit-learn and numerical gradients. Third, the five assumptions and how to diagnose each. Finally, R squared, standard errors, and fixes when an assumption breaks.

• Build
• Test
• Diagnose
• Fix

Start with the symbols. We have n rows and p features. The matrix X has shape n by p, and y has shape n. To handle the intercept, we add a column of ones to X, so it becomes n by p plus one. Now the intercept is just one more weight, and the prediction is X times w.

• X is n by p
• Prepend a ones column
• Prediction is X w

We will use a tiny dataset all lesson. Four points: x equals one, two, three, four, and y equals two, four, five, eight. The means are two point five and four point seven five. The sum of centered cross products is nine point five, and the sum of squared x deviations is five. So the slope is one point nine, and the intercept is zero.

• x: 1, 2, 3, 4
• y: 2, 4, 5, 8
• Slope 1.9, intercept 0

Now the normal equation fit. The math says w equals the inverse of X transpose X, times X transpose y. The code starts with an array of ones stacked beside X, called A. Then gram equals A transpose times A, and rhs equals A transpose times y. Each line is one symbol from the formula.

• gram is X transpose X
• rhs is X transpose y
• One symbol, one line

Notice what the code does next. It does not call an inverse. It calls solve on gram and rhs. Solving the system directly is faster and more accurate than forming the inverse and multiplying. The inverse is a notation in the formula, not an instruction for the computer.

• Inverse is notation
• solve is faster
• solve is more accurate

There is a deeper problem, though. Forming X transpose X squares the condition number. If two features are nearly identical, the condition number might be ten to the sixth, and squaring makes it ten to the twelfth. You lose about half of your digits of precision before solving anything.

• cond of XtX = cond of X squared
• Nearly equal features hurt
• Digits vanish

The stable fix is QR decomposition. Write A as Q times R, where Q has orthonormal columns and R is upper triangular. The normal equation collapses to R w equals Q transpose y, solved by back substitution. We never form A transpose A. NumPy's lstsq goes further and uses the SVD, which even handles rank deficient inputs.

• A equals Q R
• Solve R w = Q^T y
• lstsq uses SVD

Our class therefore offers a method argument. Option one is normal, with solve. Option two is qr, using numpy linalg qr and back substitution. Option three is lstsq. After fitting, we store the weights, with the intercept first. Predict is then one line: add the ones column, and multiply by the weights.

• fit with a method switch
• w_ stores all weights
• predict is X times w

Now gradient descent. With mean squared error, the gradient of the loss with respect to w is two over n, times X transpose times the residual vector, where the residual is prediction minus target. Each step moves w by minus the learning rate times that gradient. The code is literally one line for the gradient and one for the update.

• grad is 2/n X^T (Xw - y)
• Update: w minus alpha grad
• Two lines of code

Let us do one step by hand. Start with weights at zero, so every prediction is zero. The gradient for the slope is two over four, times the negative sum of x times y, which is minus twenty eight point five. For the intercept it is minus nine point five. With a learning rate of zero point zero one, the slope becomes zero point two eight five, and the intercept zero point zero nine five.

• Gradient slope: -28.5
• Gradient intercept: -9.5
• New w 0.285, b 0.095

Gradient descent has its own pitfalls. If features have very different scales, the loss surface is a long narrow valley and steps zigzag. Standardizing the features turns it into a round bowl. And if the learning rate is too large, the loss grows each step instead of falling. Always plot the loss curve.

• Unscaled features zigzag
• Standardize first
• Watch the loss curve

How do we know our gradient formula is right? A numerical gradient check. For each weight, nudge it up and down by a tiny epsilon, like ten to the minus five, compute the loss each time, and divide the difference by two epsilon. This centered difference should match the analytic gradient to many decimals.

• Nudge weight by epsilon
• Centered difference
• Compare with analytic

Next, compare against scikit-learn. Fit both on the same data, then check that the weights agree using allclose. Our normal, QR, and lstsq fits should match to about ten decimals, and gradient descent should match after enough iterations. If gradient descent is slightly off, increase iterations before blaming the code.

• Use allclose
• Direct methods match closely
• Gradient descent needs iterations

The fit is done, so now the question is trust. Predictions alone need few assumptions. But standard errors, confidence intervals, and p values rely on five. The relationship is linear. Errors are independent. Errors have constant variance. Errors are normal. And features are not strongly collinear. We check them through the residuals.

• Linearity
• Independence
• Constant variance
• Normal errors
• No multicollinearity

Residuals are the leftovers: actual minus predicted. In our toy fit they are zero point one, zero point two, minus zero point seven, and zero point four. They sum to zero, which always happens with an intercept. A good model leaves only random noise behind. Any pattern in the leftovers is a signal that the model missed something.

• Residual is y minus prediction
• Toy residuals sum to zero
• Patterns mean missed signal

The workhorse plot is residuals versus fitted values. You want a flat, formless cloud around zero. If you see a curve, like a smile or a frown, linearity is broken. The truth bends, and the line cannot follow. This is the first plot to make, every single time.

• Good: random cloud
• Curve: nonlinearity
• Make this plot first

The same plot reveals homoscedasticity, which means constant variance. If the cloud forms a funnel, narrow on the left and wide on the right, the noise grows with the prediction. This happens often with prices and incomes. Coefficients stay unbiased, but standard errors become wrong, so confidence intervals lie.

• Funnel means heteroscedasticity
• Coefficients stay unbiased
• Standard errors become wrong

For independence, plot residuals in time or row order. If positive residuals follow positive residuals, in long runs, errors are correlated. That is common in time series, and we will treat it in lesson forty six. Shuffled, independent data should look like random flicker around zero.

• Plot residuals in order
• Runs signal correlation
• Common in time series

For normal errors, use a Q-Q plot. Sort the residuals, and plot them against the quantiles a normal distribution would give. If the errors are normal, the points hug the diagonal line. Heavy tails make the ends curl away from it. Skewed errors bend the whole plot into a curve.

• Sort residuals
• Compare to normal quantiles
• Ends curling means heavy tails

A reassurance: with large samples, normality matters less for coefficients, thanks to the central limit theorem from lesson nine. It matters most for prediction intervals and for small datasets. So treat a mild Q-Q wobble as a note, not an emergency. A severe bend deserves a transform.

• CLT helps big samples
• Matters for intervals
• Mild wobble is fine

The fifth assumption, no multicollinearity. If two features carry nearly the same information, the model cannot decide which deserves the credit. Weights swing wildly and standard errors explode, even though predictions look fine. The tool is the variance inflation factor. It equals one over one minus R squared, from regressing that feature on all the others.

• Redundant features share credit
• Standard errors explode
• VIF is 1 over 1 minus R squared

Worked example. Suppose a feature can be predicted from the others with an R squared of zero point nine. The VIF is one over zero point one, which is ten. That means the variance of its coefficient is ten times larger than if it were uncorrelated, so the standard error is about three times bigger. A common rule of thumb flags values above five or ten.

• R squared 0.9 gives VIF 10
• Variance times 10
• Standard error about 3.2 times

Not every point has equal power. Leverage measures how unusual a point's features are. It comes from the diagonal of the hat matrix, which maps y to y hat. In our toy data, the leverage of x equals one and x equals four is zero point seven, and for the middle points it is zero point three. They sum to two, the number of parameters.

• Hat matrix diagonal
• Edges have more pull
• Leverages sum to parameters

High leverage alone is not harmful. A point that sits on the trend line barely changes the fit. Trouble comes when high leverage meets a large residual. Cook's distance combines them: residual squared over the number of parameters times the residual variance, times leverage over one minus leverage, squared.

• Leverage times residual
• Influence needs both
• Measures change in fit

Let us compute it. The residual variance is point seven divided by two, which is zero point three five. For the last point, residual zero point four, leverage zero point seven. Residual squared over two times zero point three five is zero point two two nine. Leverage over one minus leverage squared is seven point seven eight. Multiply: Cook's distance is about one point eight. The point with the biggest residual only scores point four three.

• Point 4: about 1.78
• Point 3: about 0.43
• Edge point has the most pull

What do you do with an influential point? Do not delete it just because it is inconvenient. First check for a data entry error. Then refit with and without it and compare the coefficients. If one point swings your conclusions, report that honestly. Common cutoffs are one, or four over n, but treat them as prompts to look, not rules.

• Check for data errors
• Refit with and without
• Report honestly

Now fit quality. R squared is one minus the sum of squared errors over the total sum of squares. In our toy data, the errors sum to zero point seven, and the total is eighteen point seven five. So R squared is one minus zero point zero three seven, about zero point nine six three. It is the fraction of variance in y that the model explains.

• SSE is 0.70
• SST is 18.75
• R squared is 0.963

But R squared never goes down when you add features, even useless ones. Adjusted R squared fixes that by penalizing parameters. It equals one minus one minus R squared, times n minus one over n minus p minus one. For us, n is four and p is one, so we get one minus zero point zero three seven times one point five, which is zero point nine four four.

• R squared always rises
• Adjusted penalizes features
• Toy value: 0.944

Next, uncertainty in the weights. The residual variance is the sum of squared errors over n minus p minus one, which is zero point three five. The variance of the weights is that times the inverse of X transpose X. For the slope this is zero point three five over five, so the standard error is zero point two six four six.

• s squared is SSE over n-p-1
• Var(w) uses inverse of XtX
• Slope SE is 0.265

The t statistic is the estimate over its standard error: one point nine over zero point two six five, about seven point two. With only two degrees of freedom, the critical t value for a ninety five percent interval is about four point three, giving roughly plus or minus one point one four around the slope. Small samples mean wide intervals, which is honest.

• t equals 1.9 over 0.265
• About 7.2
• Wide interval from n of 4

Now the fixes. For a funnel shape or right skewed errors, transform the target. Taking the log of y turns multiplicative noise into additive noise and often straightens exponential growth. Just remember that predictions come back on the log scale, so exponentiate them to return to the original units.

• Log fixes funnels
• Straightens growth curves
• Exponentiate predictions

For curvature, add polynomial terms. Include x squared as a new column. The model is still linear, because it is linear in the weights, not in x. The residual plot with a U-shape should flatten. But higher degrees can overfit and create near collinear columns, so center your features first.

• Add an x squared column
• Still linear in weights
• Center to cut collinearity

Other fixes match other problems. For multicollinearity, drop or combine redundant features. For correlated errors, use time series models. For heteroscedasticity with correct mean, robust standard errors or weighted least squares help. And for stubborn outliers, investigate their source. Always refit and re-plot the residuals afterward, because a fix is only proven when the diagnostics improve.

• Curve: polynomial terms
• Funnel: transform or weights
• Collinear: drop or combine

Let us recap. We built a LinearRegression class with solve, QR, lstsq, and gradient descent, and validated it with scikit-learn and numerical gradients. We checked five assumptions with residual plots, Q-Q plots, VIF, leverage, and Cook's distance. We measured fit with R squared, adjusted R squared, and standard errors.

• Build and test the class
• Diagnose five assumptions
• Measure and fix

The big lesson: a good score is not a good model until the residuals agree. In lesson forty one we turn to classification and derive the gradient of logistic regression from scratch.

• Residuals tell the truth
• Lesson 41 is next
• Subscribe to continue