Regression is where statistics stops describing data and starts making claims about it. This guide covers fitting linear and logistic models with statsmodels, reading the output honestly, checking whether the model was appropriate in the first place, and turning predictions into numbers a stakeholder can act on.
Two models cover an enormous share of applied work. Linear regression answers “how much?”, where the thing you are predicting is a quantity: a price, a duration, a revenue figure. Logistic regression answers “will it?”, where the outcome is one of two states: churned or stayed, converted or bounced, approved or declined.
They look like different tools and they are almost the same tool. Same formula syntax, same fit step, same prediction workflow. The difference is a single transformation applied to the response, and once you see where that transformation sits, logistic regression stops being a separate subject and becomes linear regression wearing a coat.
This guide uses statsmodels rather than scikit-learn. Both fit the same models, and the difference is what they hand back. Scikit-learn gives you a predictor optimised for accuracy on new data. statsmodels gives you the full statistical report, with standard errors, p-values and diagnostics. When you want to understand a relationship rather than just predict it, statsmodelsis the right tool.
The data
Four small fictional tables carry every example. flats holds rental listings with rent_psm for monthly rent per square metre, n_cafes counting cafes within a ten minute walk, metro_distance_m in metres, and age_band as a categorical build era. advertsholds n_views and n_signups, both heavily right skewed. stores holds growth_2023 and growth_2024 for year on year sales. And members holds cancelled, a 1 or 0 flag, against months_since_joining.
import numpy as npimport pandas as pdimport seaborn as snsimport matplotlib.pyplot as pltfrom statsmodels.formula.api import ols, logit
ols fits linear models and logit fits logistic ones, both through the formula interface that section 2 explains.
1. Look at the data before you model it
Fitting first and plotting later is how people end up publishing a straight line through a curve.
sns.regplot(x="n_cafes", y="rent_psm", data=flats, ci=None, scatter_kws={"alpha": 0.5})plt.show()sns.displot(data=flats, x="rent_psm", col="age_band", bins=10)plt.show()
sns.regplot is the workhorse because it does two jobs in one call: it draws the scatter and overlays the fitted regression line, so you can see immediately whether a straight line is a defensible summary of what is in front of you. ci=None removes the confidence band, which declutters a busy plot. For a categorical predictor, displot with col= gives you one histogram per category side by side.
The alpha=0.5 matters more than it looks. With a few thousand overlapping points an opaque scatter collapses into a solid blob and you lose all sense of where the data actually concentrates. Translucent points restore that density information for free.
For a binary response, plot both candidate shapes on the same axes and let the picture make the argument.
sns.regplot(x="months_since_joining", y="cancelled", data=members, ci=None, line_kws={"color": "red"}) # straight linesns.regplot(x="months_since_joining", y="cancelled", data=members, ci=None, line_kws={"color": "blue"}, logistic=True) # S-curveplt.show()
logistic=True fits an S-shaped curve instead of a straight line, and line_kws styles the trend without touching the points. Plotting both together shows the straight line running past 1 and below 0, which a probability cannot do. Keep ci=None on the logistic call in particular, because its confidence band is computed by bootstrapping and is slow on large data.
2. Fitting a linear model
mdl = ols("rent_psm ~ n_cafes", data=flats).fit()print(mdl.params)print(mdl.summary())
OLS stands for ordinary least squares, which describes exactly what it does: it finds the line minimising the sum of squared vertical distances between the line and the points.
The formula string is R-style syntax, read as “response modelled by explanatory”. Whatever sits left of the tilde is what you are predicting and whatever sits right of it does the predicting. It looks unusual coming from scikit-learn’s X and y, but it pays off quickly: a second predictor is "rent_psm ~ n_cafes + metro_distance_m" and an interaction is "rent_psm ~ n_cafes * metro_distance_m", with no matrix assembly required.
The two-step pattern of ols(...) then .fit() is universal in statsmodels. The first call defines the model and the second estimates it. Afterwards, mdl.params gives the coefficients as a Series with the intercept first, and mdl.summary() prints the full report: coefficients, standard errors, p-values, R-squared and the F-statistic.
The fitted model is the line:
Where beta-zero is the intercept, the predicted rent when n_cafes is zero, and beta-one is the slope, the change in predicted rent per additional cafe.
Read that slope carefully. It is the average difference in rent between listings that differ by one cafe, holding nothing else constant because nothing else is in the model. It is not the amount rent would rise if someone opened a cafe. Regression coefficients describe associations in the data you fitted, and turning one into a causal claim requires an argument the data alone cannot supply.
3. Categorical predictors
statsmodels handles a text column automatically, and the default output is awkward to read.
mdl_cat = ols("rent_psm ~ age_band", data=flats).fit() # baseline codingmdl_cat = ols("rent_psm ~ age_band + 0", data=flats).fit() # group means
By default one category is absorbed into the intercept and every other coefficient reports a difference from that hidden baseline. It is mathematically correct and mentally expensive, because reading off the average rent for a category means finding the intercept and adding the right coefficient to it.
Adding + 0 removes the intercept, and each category coefficient then becomes that group’s mean directly, readable straight off mdl.params. For a single categorical predictor this is almost always the version you want to look at. If your categories are coded as numbers rather than text, wrap the column in C(...) to stop statsmodels treating them as a continuous variable.
One trap here is real and easy to fall into. Drop the intercept and the reported R-squared jumps, sometimes to something like 0.98, because it is now measured against zero rather than against the mean of the response. The model has not improved. Use + 0 to read group means, and go back to the standard form whenever you want to judge fit.
4. Making predictions
There is one prediction workflow and it has three steps.
explanatory_data = pd.DataFrame({"n_cafes": np.arange(0, 11)})predictions = mdl.predict(explanatory_data)prediction_data = explanatory_data.assign(rent_psm=predictions)
Build a DataFrame of the x-values you want predictions for, call predict to get a Series of predicted responses, then assignthem back as a new column. assign returns a new frame rather than modifying the original, which keeps the explanatory grid reusable.
The detail that catches people is that predict() wants a DataFrame whose column names match the ones used in the formula, even when you are predicting a single value. A bare list or array sometimes works and sometimes fails in a way that is hard to trace. Build the DataFrame every time and the problem never arises.
Then check your work visually.
sns.regplot(x="n_cafes", y="rent_psm", data=flats, ci=None)sns.scatterplot(x="n_cafes", y="rent_psm", data=prediction_data, color="red")plt.show()
The red points must sit exactly on the regression line, because that line is what generated them. Any red point off the line means a coding error somewhere: wrong model, wrong column, or a stale prediction frame. It is a thirty second check that catches a whole category of bugs, and it is worth doing every time.
5. Putting an interval on a prediction
A point prediction is one number offered as if it were the answer. Two different intervals can go around it, and confusing them is common.
pred = mdl.get_prediction(explanatory_data)frame = pred.summary_frame(alpha=0.05)frame[["mean", "mean_ci_lower", "mean_ci_upper", "obs_ci_lower", "obs_ci_upper"]]
mean_ci_lower and mean_ci_upper bound the average response at that x-value. They answer “across all flats with four nearby cafes, where is the true mean rent?” That interval reflects only the uncertainty in the fitted line, so it narrows as your sample grows and it is quite tight in the middle of your data.
obs_ci_lower and obs_ci_upper bound a single new observation. They answer “if I look at one more flat with four nearby cafes, what rent might it have?” That interval is always substantially wider, because it carries both the uncertainty in the line and the scatter of individual points around it, and it does not shrink towards zero however much data you collect.
Which one you want depends on the question. Forecasting a portfolio average wants the first. Advising one landlord on one property wants the second, and quoting the confidence interval there would understate the range by a large margin. alpha=0.05 gives 95% intervals and alpha=0.1 gives 90%.
Reporting a prediction without either is the regression equivalent of quoting a survey result with no margin of error.
6. Doing it by hand once
There is no machinery hidden inside predict(). It is the line equation and nothing else.
intercept = mdl.params["Intercept"]slope = mdl.params["n_cafes"]manual_predictions = intercept + slope * explanatory_data["n_cafes"]
That arithmetic reproduces exactly what mdl.predict() computes internally.
Index by name, not by position. You will see mdl.params[0] in older material, but positional indexing on a Series is deprecated in current pandas and raises a warning. More importantly, position is fragile: add a second predictor and params[1] silently starts referring to something else, while params["n_cafes"] cannot break that way.
7. Transforming when the relationship curves
Linear regression fits a straight line. When the scatter is clearly bent you have two options: use a different model, or bend the data until a straight line fits it.
# Transform the explanatory variableflats["sqrt_metro_distance"] = np.sqrt(flats["metro_distance_m"])mdl_dist = ols("rent_psm ~ sqrt_metro_distance", data=flats).fit()# Transform both sidesadverts["signups_4rt"] = adverts["n_signups"] ** 0.25adverts["views_4rt"] = adverts["n_views"] ** 0.25mdl_ads = ols("signups_4rt ~ views_4rt", data=adverts).fit()# Back-transform the predictions before reporting themprediction_data["n_signups"] = prediction_data["signups_4rt"] ** 4
Which transformation to reach for depends on how severe the curve is. A square root tames moderate right skew and squares back. A fourth root, written as ** 0.25, compresses much harder and suits count data spanning several orders of magnitude like ad impressions, and raising to the fourth power inverts it. A log handles relationships that are multiplicative rather than additive, where a predictor changes the response by a percentage rather than a fixed amount, and np.exp inverts it.
The step people forget is the last one. Fit on a fourth root, report the raw coefficients, and you are telling a stakeholder that an extra unit of fourth-root-views buys 0.42 fourth-root-signups, which means nothing to anyone. Predict on the transformed scale, then invert before the number leaves your notebook.
Two cautions. Log and square root are undefined at or below zero, so a column containing zeros needs np.log1p or a small constant added first, and be explicit about which you did. And back-transforming a prediction from a log model gives you the median of the response rather than the mean, because the transformation is not linear. For most reporting the distinction does not matter, and it is worth knowing the number is not quite what you might assume.
8. How well does the model fit
print(mdl.rsquared)print(mdl.rsquared_adj)rse = np.sqrt(mdl.mse_resid)
R-squared runs from 0 to 1 and reports the proportion of variation in the response the model accounts for. Zero means the model does no better than always guessing the mean, and one means every point sits exactly on the line.
Residual standard error is the more useful number day to day, because it is expressed in the units of the response:
If you are modelling rent in pounds per square metre and the RSE is 3.2, your predictions are typically off by about 3.2 pounds per square metre. That is a sentence anyone can act on. mse_resid carries the same information in squared units, which nobody can interpret directly, so the square root is not optional.
Use rsquared_adj rather than rsquared the moment you have more than one predictor. Plain R-squared can never decrease when you add a term, even a column of pure noise, so it will always favour the bigger model. Adjusted R-squared charges you for each term and will fall when a predictor is not earning its place.
And be sceptical of R-squared as a verdict at all. A high value on a model with the wrong shape is a well-fitted wrong answer, and the next section is how you find that out.
9. Diagnostics, and what to do about them
R-squared tells you how close the points are to the line. It says nothing about whether a line was the right thing to draw. Three plots answer that, and together they are the standard regression checkup.
# 1. Residuals against fitted values: should be a flat band centred on zerosns.regplot(x=mdl.fittedvalues, y=mdl.resid, ci=None, lowess=True)plt.axhline(y=0, linestyle="dotted", color="grey")plt.xlabel("Fitted values")plt.ylabel("Residuals")plt.show()# 2. Q-Q plot: are the residuals normally distributedfrom statsmodels.api import qqplotqqplot(data=mdl.resid, fit=True, line="45")plt.show()# 3. Scale-location: is the spread of residuals constantresiduals_std = mdl.get_influence().resid_studentized_internalsqrt_abs_resid = np.sqrt(np.abs(residuals_std))sns.regplot(x=mdl.fittedvalues, y=sqrt_abs_resid, ci=None, lowess=True)plt.show()
mdl.fittedvalues gives the model’s prediction for every training row and mdl.resid gives the gap between that and what actually happened. lowess=True overlays a locally smoothed curve, which is what makes any trend in the residuals visible rather than something you have to squint for.
Plotted against the fitted values, those residuals should look like structureless noise around zero. If the smoothed line curves, the model is missing a pattern that is plainly there in the data, which is a signal to transform a variable or add a term.
The Q-Q plot compares residual quantiles against a normal distribution, with line="45" drawing the reference the points should follow. Departures at the ends are common and usually tolerable at large sample sizes. Departures through the middle are not, and they undermine the p-values and confidence intervals in the summary, which assume normal residuals.
The scale-location plot checks whether the residual spread is constant across the range of predictions. resid_studentized_internal rescales the residuals so they are comparable, and taking the square root of the absolute value folds them to one side so the plot shows spread rather than direction. A flat band is what you want.
A fan shape, where errors grow as predictions grow, is heteroscedasticity, and this is the diagnostic most often reported and then ignored. The coefficients survive it. The standard errors do not, so every p-value in the summary becomes unreliable and the significance stars stop meaning what they appear to mean.
The fix is one argument:
mdl_robust = ols("rent_psm ~ n_cafes", data=flats).fit(cov_type="HC3")
cov_type="HC3" computes heteroscedasticity-consistent standard errors. The coefficients come back identical, because least squares still estimates them the same way, but the standard errors, t-statistics, p-values and confidence intervals are all recomputed in a form that stays valid when the residual variance is not constant. HC3 is the recommended variant for small and moderate samples, and HC1 is common in econometrics.
Robust standard errors cost you almost nothing when the residuals were well behaved anyway, so on any observational data where you have not specifically confirmed constant variance, they are a reasonable default rather than a remedy.
You will also see the shortcut sns.residplot(x="n_cafes", y="rent_psm", data=flats, lowess=True), which fits the simple regression internally and plots its residuals in one call. It is convenient, and note that it puts the explanatory variable on the x-axis rather than the fitted values. For a single predictor the two differ only by a rescaling of the axis so the shape is identical, but label the axis for what it is. With more than one predictor, fitted values are the only correct choice.
10. Leverage and influence
Not every row pulls on the line equally.
summary_info = mdl.get_influence().summary_frame()flats["leverage"] = summary_info["hat_diag"]flats["cooks_dist"] = summary_info["cooks_d"]print(flats.sort_values("cooks_dist", ascending=False).head())
get_influence().summary_frame() returns a DataFrame of per-row diagnostics, of which two matter.
hat_diag is leverage, measuring how unusual a row’s explanatory values are. A listing 40 minutes from the nearest station, when almost every other listing is within 10, has high leverage. That is potential rather than guilt: the point is positioned to move the line, but if it sits where the line already goes it moves nothing.
cooks_d is Cook’s distance, which combines position with outcome by measuring how much the fitted line would shift if you deleted that single row. High leverage plus a large residual gives a high Cook’s distance, and sorting descending on it surfaces exactly the rows worth opening up individually.
What to do about them is a judgement call rather than a rule. A high-influence point can be a data entry error, in which case fix or drop it. It can be a genuine but rare case, in which case it belongs in the model and its influence is real information. It can also be a sign the model is too simple to cover the full range of the data. The one thing not to do is delete points because removing them improves R-squared. That is not cleaning, it is fitting the data to the conclusion.
11. Regression to the mean
sns.regplot(x="growth_2023", y="growth_2024", data=stores, ci=None)plt.axline(xy1=(0, 0), slope=1, linewidth=2, color="green")plt.axis("equal")plt.show()
Measure anything twice and the extremes drift towards the middle. The stores with the wildest growth in 2023 will on average look more ordinary in 2024, and the worst performers will look less bad, even when nothing about the stores changed. An extreme result usually combines genuine performance with a helping of luck, and the luck does not repeat.
plt.axline(xy1=(0, 0), slope=1) draws the y equals x reference line, and plt.axis("equal") forces the same scale on both axes, without which that line is visually misleading. Points above it improved between periods and points below fell back. If regression to the mean is operating, the fitted line comes out flatter than the reference line, crossing it near the average and pulling both extremes inwards.
This is why last year’s top fund is a poor bet for next year, and why any intervention targeted specifically at the worst performers will appear to work. Select the bottom 10% and measure them again and they improve, whether or not you did anything at all. Any evaluation that selects on an extreme and lacks a control group is measuring this effect, at least in part.
12. Logistic regression
When the response is binary, a straight line is the wrong shape. It has no upper or lower bound, so it will confidently predict a 130% probability of cancellation for a long-standing member and a negative probability for a new one.
Logistic regression fits an S-curve bounded between 0 and 1 instead.
mdl_logit = logit("cancelled ~ months_since_joining", data=members).fit()print(mdl_logit.params)print(mdl_logit.summary())
The API is deliberately identical to ols: same formula string, same .fit(), same .predict(). The response must be binary, coded as 0 and 1 or as two categories.
What changes underneath is that coefficients are estimated by maximum likelihood rather than least squares, and the linear part of the model predicts log odds rather than the response itself:
Rearranging for the probability gives the S-curve:
Which is why mdl.params comes back on the log odds scale while .predict() returns probabilities.
Two things in the summary differ from ols. There is no R-squared, because the concept does not transfer, and mdl.prsquaredgives a pseudo R-squared that is useful for comparing models against each other but does not mean “proportion of variance explained”. And if .fit() reports perfect separation, some predictor splits the two classes cleanly, so the likelihood has no finite maximum and the coefficients run off towards infinity. That is usually a leaked variable encoding the answer rather than a genuine finding.
13. Four ways to express a logistic prediction
The same prediction sits on four scales, and the right one depends entirely on who is reading it.
explanatory_data = pd.DataFrame({"months_since_joining": np.arange(0, 41)})prediction_data = explanatory_data.assign( cancelled=mdl_logit.predict(explanatory_data) # 1. probability)prediction_data["most_likely"] = np.round(prediction_data["cancelled"]) # 2.prediction_data["odds"] = ( # 3. odds prediction_data["cancelled"] / (1 - prediction_data["cancelled"]))prediction_data["log_odds"] = np.log(prediction_data["odds"]) # 4. log odds
Probability is what you put in front of stakeholders. “A 70% chance of cancelling” needs no explanation.
Most likely outcome rounds at 0.5 to give a class label, which is what you need when the output drives a yes or no decision like whether to trigger a retention email.
Odds is the ratio of the probability of the event to the probability of it not happening:
/math \text{odds} = \frac{p}{1 – p}
Odds of 1 is an even chance and odds of 4 means four to one in favour. This is the framing anyone who has bet on anything already understands.
Log odds is the scale on which the model is genuinely linear. It is unintuitive to report and it is where the coefficients live, so a slope of 0.08 means each additional month adds 0.08 to the log odds of cancelling. Running a coefficient through np.exp() converts it into an odds ratio, which is the form most often quoted in published research.
Rounding at 0.5 is a choice, not a law. That threshold treats a false positive and a false negative as equally costly, which is rarely true. If missing a cancellation costs you a customer and a false alarm costs you one email, you should be cutting well below 0.5. Compare prediction_data["cancelled"] > 0.3 against the default and pick the threshold matching the actual cost of each mistake.
14. Confusion matrix and performance
conf_matrix = mdl_logit.pred_table()# Predicted 0 Predicted 1# Actual 0 [ TN FP ]# Actual 1 [ FN TP ]TN, FP = conf_matrix[0]FN, TP = conf_matrix[1]accuracy = (TP + TN) / (TP + TN + FP + FN)sensitivity = TP / (TP + FN) # recall, the true positive ratespecificity = TN / (TN + FP) # the true negative rate
pred_table() returns the two by two table of what happened against what the model said would happen, with rows as actual and columns as predicted. It cuts at 0.5 by default, and pred_table(threshold=0.3) moves that.
Four cells, four names. True negatives are correct nos, true positives are correct yeses, false negatives are the yeses you missed, and false positives are the yeses you invented. Three metrics come out of them:
Accuracy asks how often the model is right overall. Sensitivity asks how many of the real positives it caught. Specificity asks how many of the real negatives it correctly left alone.
Accuracy is the metric everyone reaches for and the one most likely to mislead. If 95% of your members never cancel, a model predicting “nobody cancels” for every single row scores 95% accuracy while being completely useless. Its sensitivity is zero, which is the number that exposes it. Whenever the classes are imbalanced, report sensitivity and specificity alongside accuracy or the headline figure will flatter a model that has learned nothing.
Which of the two matters more depends on the cost of each mistake. For fraud detection or a medical screen, a missed positive is far worse than a false alarm, so you optimise sensitivity and accept the false positives. For anything triggering an expensive or intrusive action, false positives cost real money and specificity matters more.
The matrix also plots, with each cell drawn as a rectangle whose area is proportional to its count, which shows at a glance which error type dominates:
from statsmodels.graphics.mosaicplot import mosaicmosaic(conf_matrix)plt.show()
15. Common pitfalls
Reading a slope as a cause. A coefficient describes an association in the data you fitted. Turning it into “if we change x, y will follow” needs an argument about how the data was generated, which the regression cannot provide.
Trusting R-squared without looking at residuals. A high R-squared on a curved relationship is a well-fitted wrong answer.
Comparing R-squared across a + 0 boundary. Dropping the intercept changes how it is calculated and inflates it.
Using plain R-squared with several predictors. It never decreases when you add a term, so it always prefers the bigger model. Use rsquared_adj.
Indexing coefficients by position. mdl.params[0] is deprecated in current pandas and breaks silently when the formula changes. Use the name.
Diagnosing heteroscedasticity and then doing nothing. fit(cov_type="HC3") fixes the standard errors in one argument.
Reporting a point prediction with no interval. And when you do add one, mean_ci is for the average response and obs_ciis for a single new observation. They are not interchangeable.
Forgetting to back-transform. A prediction from a fourth-root model stays on the fourth-root scale until you raise it to the fourth power.
Transforming through zero. np.log and np.sqrt are undefined at or below zero. Use np.log1p or shift the column.
Treating 0.5 as the natural threshold. It assumes false positives and false negatives cost the same.
Reporting accuracy on imbalanced classes. Always pair it with sensitivity and specificity.
Hitting perfect separation in logit. Look for a leaked column before you look for a fix.
17. The mental model
LINEAR REGRESSION the response is a quantity ols("y ~ x", data).fit() y_hat = b0 + b1 * x coefficients read directly in y units fit: rsquared_adj, sqrt(mse_resid) check: residuals vs fitted, Q-Q, scale-location if it fans: refit with cov_type="HC3" outliers: hat_diag for leverage, cooks_d for influence report: get_prediction, mean_ci for averages, obs_ci for one caseLOGISTIC REGRESSION the response is binary logit("y ~ x", data).fit() log(p / (1-p)) = b0 + b1 * x coefficients live on the log odds scale predict: probability -> round -> class label scales: probability -> odds -> log odds -> exp() for an odds ratio check: pred_table, then sensitivity and specificity, not accuracySHARED WORKFLOW plot the data -> is a line, or an S-curve, even plausible? fit the model -> ols or logit, same formula syntax build a grid -> DataFrame of explanatory values predict -> mdl.predict(grid), or get_prediction for intervals diagnose -> residuals for linear, confusion matrix for logistic report -> back-transform first, with an interval, in real units
Three habits separate a regression you can defend from one you cannot, and all of them are cheap. Plot before you fit, so you know whether the shape is right. Plot the residuals after you fit, so you know whether it stayed right, and act on what they show rather than noting it and moving on. And say out loud what the coefficient actually claims, because the moment you phrase it as a causal statement you will hear whether the data can support it.
Fit the model, then interrogate it. The summary table is where the work starts, not where it ends.
Thanks for reading, Andrei.
[…] Linear and Logistic Regression in Python […]