Most people learn linear regression as a function call and never recover from it. Calling fit is not the hard part, and it never was. Hard parts are knowing whether the number that came back means anything, and recognising the specific ways a linear model goes wrong while reporting excellent scores. Regression is the one model where you can see the entire mechanism, so this is the last chance in the series to build a real intuition before the models stop being inspectable.

Everything below was run against Python 3.10.12, NumPy 2.2.6 and pandas 2.3.3, with scikit-learn 1.7.2 as the installed version. Release 1.9.0 is the current one as of June 2026 and nothing used here changed between them. Numbers quoted are from the run, not from memory.

Where the churn project stands after Part 8

Last part we had a churn rate with an honest interval around it and a rule for comparing two models without kidding ourselves. What we did not have was a model. This part adds the first one, and deliberately picks a regression target rather than churn itself, because churn is a yes or no outcome and pointing a linear model at it produces predicted probabilities below zero. Part 10 handles that properly.

Instead the target here is monthly revenue per customer. That matters to the churn project directly: revenue at risk is churn probability multiplied by monthly charge, and finance will ask where the second number comes from long before they ask about the first. Subscription pricing is close to additive by construction, a base plus a line item per service, which makes it the honest case for a linear model rather than a forced one.

One planning note that shaped this whole part. Before pointing a fitter at a table where I do not know the answer, I build a table where I do. That habit cost me an hour to set up and has caught three genuine bugs since, and it is what the code below actually runs on.

Least squares in three lines of numpy

Linear regression assumes your target is a weighted sum of your features plus an intercept, and picks the weights that make the total squared error as small as possible. Squared error, not absolute error, and that choice is why the answer has a closed form. Differentiate the squared error, set the derivative to zero, and you get the normal equations: X transpose X times beta equals X transpose y. Solve for beta and you are done.

Geometrically, the fitted values are the projection of y onto the space your columns can reach, and residuals are what is left over, perpendicular to every column. That perpendicularity is why residuals sum to zero when you include an intercept, and it is worth holding on to, because it explains half the diagnostics later.

import numpy as np

# Xd already has a leading column of ones for the intercept
beta_ne = np.linalg.solve(Xd.T @ Xd, Xd.T @ y)      # normal equations
beta_ls, _, rank, sv = np.linalg.lstsq(Xd, y, rcond=None)  # what libraries use

print('max abs diff:', f'{np.abs(beta_ne - beta_ls).max():.3e}')

Output:

max abs diff: 2.842e-13

Both routes agree to thirteen decimal places, so on well conditioned data they are the same estimator. They stop agreeing when the data is badly conditioned, and that difference is the subject of a later section. Use lstsq in anything you ship. It goes through a singular value decomposition rather than inverting a matrix, which is slower and much harder to fool.

Calibration data with a known answer

Here is the practice I would push hardest in this whole part. Generate a table with the same shape and column types as your real one, but build the target yourself from coefficients you choose. Then fit it. If your pipeline cannot recover numbers you personally put in, it will not recover numbers you did not, and you have just found the bug for free instead of finding it in a review three weeks later.

My calibration table mirrors the subscription schema: an internet tier with three levels, six optional add-ons that only exist when internet exists, and a multi line flag. Prices are a base of 19.95, plus 25.10 for DSL or 45.30 for fibre, plus 5.00 for each of four support add-ons, plus 10.00 for each of two streaming add-ons, plus gaussian noise with a standard deviation of 2.00. Seven thousand and forty three rows, matching the extract the series uses.

X = pd.get_dummies(df.drop(columns='MonthlyCharges'),
                   columns=['InternetService'], drop_first=True).astype(float)
Xd = np.column_stack([np.ones(len(X)), X.values])

beta, _, rank, _ = np.linalg.lstsq(Xd, y, rcond=None)
resid = y - Xd @ beta
n, p = Xd.shape

R2  = 1 - (resid**2).sum() / ((y - y.mean())**2).sum()
s2  = (resid**2).sum() / (n - p)                 # unbiased error variance
se  = np.sqrt(np.diag(np.linalg.inv(Xd.T @ Xd)) * s2)

print(f'n={n} p={p} rank={rank} R2={R2:.4f} sigma={np.sqrt(s2):.3f}')

Output:

n=7043 p=10 rank=10 R2=0.9938 sigma=2.009

Sigma came back as 2.009 against a true noise level of 2.000, which is the first thing I check. If the estimated noise does not match the noise I injected, nothing downstream is trustworthy. Here is the full coefficient table against ground truth.

Reading coefficients without lying to yourself

Notice what the intercept is in that table. It is 45.05, not the 19.95 base price I injected. Nothing is wrong. Pandas dropped the DSL level alphabetically when I passed drop_first=True, so the baseline customer is a DSL subscriber, and the intercept is base plus DSL, which is 45.05 exactly. Fibre then reads as 20.14 because it is the gap between fibre and DSL rather than the fibre price. Every categorical coefficient is measured against whichever level got dropped, and if you cannot say which level that is, you cannot interpret a single row of your own output.

Second habit worth building: read every coefficient with the phrase holding the other columns fixed attached to it. Fibre being worth 20.14 more than DSL is a claim about two customers identical in every other column. When your columns cannot vary independently in the real world, that comparison describes customers who do not exist, and the number is arithmetic rather than insight. Standard errors carry the same caveat, and the interval logic from Part 8 applies directly here: a coefficient plus or minus two standard errors is the range the data actually supports.

Those t statistics in the 90 to 600 range are not something to celebrate either. They are that large because I generated data with almost no unexplained variation. On the real extract nothing looks like this, and a t statistic above 100 in production work is usually a sign that a column is leaking the target rather than predicting it.

Collinearity and a silent failure worth fearing

Drop the drop_first argument and keep every level of the categorical, and the three internet columns now sum to exactly one on every row, which is exactly what the intercept column already does. Your design matrix has eleven columns and a rank of ten. One column is redundant, and there is no longer a unique answer.

Xall = pd.get_dummies(df.drop(columns='MonthlyCharges'),
                      columns=['InternetService'], drop_first=False).astype(float)
Xa = np.column_stack([np.ones(n), Xall.values])

print('shape', Xa.shape, 'rank', np.linalg.matrix_rank(Xa))
print('cond %.3e' % np.linalg.cond(Xa))

np.linalg.inv(Xa.T @ Xa)          # does this raise?
b2, _, _, _ = np.linalg.lstsq(Xa, y, rcond=None)

Output:

shape (7043, 11) rank 10
cond 5.270e+14

const                          32.554
InternetService_DSL            12.498
InternetService_Fiber          32.634
InternetService_No            -12.577

Read that carefully, because the important thing is what is missing. No exception. np.linalg.inv did not raise LinAlgError: Singular matrix, because floating point rounding leaves the matrix technically invertible even though it is mathematically singular. Predictions stay perfect, R squared is unchanged, and every coefficient is now garbage. Fibre reads 32.63 against a true 45.30. The intercept is 32.55 against a true 45.05. A condition number of 5.27e14 is the only thing on that page that tells you.

Near collinearity is the more common version and it is worse, because the rank check passes. Two columns at a correlation of 0.98 will give you coefficients that swing wildly between refits while the model scores identically. If you only ever look at the score, you will never see it.

Residual plots that earn their place

R squared tells you how much variance you captured and nothing about whether the linear form was right. Residuals tell you the second thing, which is the one that decides whether you should keep going. Two plots are worth the effort and the rest usually are not: residuals against fitted values, which exposes curvature and changing spread, and a histogram of residuals, which exposes skew and outliers.

Splitting my fitted values into octiles, mean residual per bucket ran from -0.075 to +0.126 against a residual standard deviation of 2.009, so under 7 percent of a standard deviation of drift across the whole range. Flat, as it should be when the model matches how the data was made. On a real table you look for a smile shape, which means you needed a squared term, or a widening fan, which means the error grows with the prediction and you should probably model the logarithm instead.

One clarification that saves arguments. Normality of residuals is not required for the coefficients to be unbiased, and at 7043 rows it is not required for the standard errors either. It matters for small samples and for prediction intervals on individual customers. Plenty of teams throw away a perfectly usable model over a failed normality test on 50000 rows, where the test will reject any real dataset regardless.

Ridge regression as a trade of bias for variance

Ridge adds a penalty on the size of the coefficients, which stabilises the solve when columns are nearly duplicated. To show what you are buying and what you are paying, I built two features correlated at 0.9975 where only the first one genuinely drives the target with a weight of 3.0, and the second contributes nothing.

P = np.eye(3); P[0, 0] = 0        # never penalise the intercept
for lam in [1.0, 10.0, 100.0]:
    br = np.linalg.solve(Xc.T @ Xc + lam * P, Xc.T @ yy)
    print('ridge %6.1f b1=%.3f b2=%.3f' % (lam, br[1], br[2]))

Watch what heavy penalties do. At alpha 100 ridge has split the credit almost evenly between the real driver and the passenger, 1.715 against 1.268, and it is now confidently wrong about which feature matters. Predictions barely suffer, because the two columns are interchangeable for prediction purposes. Attribution is destroyed. If somebody is going to read your coefficients as a story about the business, a large alpha will make that story false while the scoreboard stays green.

My verdict: reach for ridge when you need predictions and have more columns than you can defend, and always standardise your features first, because the penalty is scale dependent and an unscaled feature gets penalised on the accident of its units. Avoid ridge when the coefficients themselves are the deliverable. Drop or combine the duplicated columns instead, which is a decision you can explain, unlike an alpha you tuned.

Prediction intervals against confidence intervals

Two intervals come out of a fitted regression and they get confused constantly, including by people who should know better. A confidence interval on the fitted value asks how well you know the average monthly charge for customers with this exact configuration. A prediction interval asks what range a single new customer will land in. Same model, same row, wildly different widths.

Arithmetic makes the gap obvious. Confidence interval width shrinks toward zero as your sample grows, because averages become knowable. Prediction interval width converges on the irreducible noise instead, which in my calibration run was a standard deviation of 2.009, so roughly plus or minus 3.94 at 95 percent no matter how many rows I add. Ten million customers would not narrow it, because individual customers vary and no amount of data removes that.

Practical consequence for the churn project: a revenue at risk figure aggregated over a segment can lean on confidence intervals and will be tight. A promise about what one named account will bill next month needs the prediction interval and will be embarrassingly wide. I have watched a forecast lose credibility because someone quoted the narrow interval for a single account question, and the account came in 6.00 outside it in the first month. Match the interval to the question being asked, and say out loud which one you used.

Setup I would use for a first regression

Fit with scikit-learn if the model is going into a pipeline, and refit the same specification in statsmodels when you have to explain it, because summary() gives you standard errors, t statistics and a multicollinearity note in one call and scikit-learn deliberately does not. Running both is not duplicated work, it is the cheapest review you will ever get, and the coefficients must agree to many decimal places or you have specified two different models.

Concretely, four things go in every regression script I write. Print the shape, the rank and the condition number of the design matrix before fitting. Name the dropped level for every categorical in a comment. Compare estimated residual noise against something you expect. Plot residuals before reporting a score to anybody. None of that is sophisticated and all of it has caught something real.

Part 10 keeps this design matrix and changes one thing: the target goes back to churn itself, and a link function keeps predictions between zero and one. Almost everything here carries over, including the collinearity trap, which behaves identically and is harder to spot. If your descriptive statistics are rusty before that, the analyst series covers the ground in statistics analysts actually use. Open a notebook, generate a table with coefficients you chose, and confirm you can get them back. If you cannot, you have learned something more valuable than a good score.

References

• statsmodels 0.14.6: Linear Regression, model classes and results classes
• scikit-learn: LinearRegression API reference
• scikit-learn: release notes for version 1.9
• NumPy: numpy.linalg.lstsq

About The Author

Dr. Pranay Jha

Dr. Pranay Jha is a Cloud and AI Consultant with 18+ years of experience in hybrid cloud, virtualization, and enterprise infrastructure transformation. He specializes in VMware technologies, multi-cloud strategy, and Generative AI solutions. He holds a PhD in Computer Applications with research focused on Cloud and AI, has published multiple research papers, and has been a VMware vExpert since 2016 and a VMUG Community Leader.

See author's posts