Regression plots
Seaborn has three functions dedicated to showing a relationship between two numeric variables together with a fitted trend: sns.regplot(), sns.lmplot() and sns.residplot(). They all do the same core work — fit a model behind the scenes and draw the result on top of a scatter — but they differ in scope and in what they let you control.
Before anything else, one sentence that should stay with you for the whole chapter:
That is not a flaw. It is a division of labour. Seaborn answers "does this relationship look linear, and how noisy is it?" — statsmodels, scipy or scikit-learn answer "what exactly is the model, and how confident am I in each coefficient?"
Throughout this chapter we use four built-in datasets:
import numpy as np import pandas as pd import matplotlib.pyplot as plt import seaborn as sns sns.set_theme(style="whitegrid") tips = sns.load_dataset("tips") # bill vs tip, continuous x and y penguins = sns.load_dataset("penguins") # body measurements, clean group structure titanic = sns.load_dataset("titanic") # binary outcome for logistic fits anscombe = sns.load_dataset("anscombe") # the classic four-dataset warning
19.1 regplot() and what the shaded band means
sns.regplot() is the axes-level function. It draws two things into a single Matplotlib Axes:
- a scatter of the raw observations;
- a fitted line with a shaded band around it.
sns.regplot(data=tips, x="total_bill", y="tip") plt.show()
The shaded band is a 95% confidence interval for the line itself — how uncertain the average relationship is. It is not a prediction interval for individual points, which is why most points fall outside it.
By default the model is ordinary least squares: a straight line y = a + b*x chosen to minimise the sum of squared vertical distances between the points and the line. The band is drawn by bootstrapping — seaborn resamples the rows many times, refits the line on each resample, and shows the middle 95% of the fitted lines.
The band is a confidence interval for the line, not for the points
This is the single most misread element in the whole seaborn library. Students look at a narrow band and conclude "the model predicts tips within ±0.2 dollars". That is wrong.
| Interpretation | Correct? | Why |
|---|---|---|
| 95% confidence interval for the mean value of y at each x | Yes | This is exactly what seaborn draws |
| Prediction interval for a single new observation | No | A prediction interval is much wider; it must include the residual noise, which the CI ignores |
| Standard deviation of y at that x | No | The spread of points is visible in the scatter itself, not in the band |
| A measure of how good the model is | No | The band shrinks as n grows even if the model is completely wrong |
The width of the band scales roughly with 1/√n. With 50,000 rows the band becomes a hairline even if the true relationship is a parabola and your straight line misses it entirely. Narrow band = "I know where the average line is"; it says nothing about "the average line is the right shape".
Imagine measuring the average height of students in each year group. With 1000 students per year you know each group's average height very precisely — that precision is the band. It does not mean every student is that height. Individual students still vary by 15 cm — that variation is the scatter cloud. The band describes your knowledge of the average, not the spread of individuals.
Controlling the band
# turn the band off entirely — faster and less misleading in dense plots sns.regplot(data=tips, x="total_bill", y="tip", ci=None) # use a different confidence level sns.regplot(data=tips, x="total_bill", y="tip", ci=68) # control the bootstrap count (default 1000); lower is faster, noisier sns.regplot(data=tips, x="total_bill", y="tip", n_boot=200)
Styling the two layers separately
regplot() draws a scatter and a line, so it accepts two dictionaries of keyword arguments that are forwarded to the underlying Matplotlib calls:
sns.regplot( data=tips, x="total_bill", y="tip", scatter_kws={"alpha": 0.4, "s": 25, "color": "#7a7a7a"}, line_kws={"color": "crimson", "linewidth": 2}, )
With more than a few hundred points the scatter turns into a solid blob and the line loses context. Set scatter_kws={"alpha": 0.3} so that density becomes visible: dark regions mean many overlapping observations.
Because it is axes-level, it composes
fig, axes = plt.subplots(1, 2, figsize=(10, 4), sharey=True) sns.regplot(data=tips, x="total_bill", y="tip", ax=axes[0]) sns.regplot(data=tips, x="size", y="tip", ax=axes[1]) axes[0].set_title("Tip vs bill") axes[1].set_title("Tip vs party size") fig.tight_layout()
You can also stack regplot() on top of another seaborn or Matplotlib call in the same Axes — for example a KDE background, or a second regression line from a filtered subset.
19.2 lmplot() and the regplot vs lmplot distinction
sns.lmplot() is the figure-level counterpart. Internally it calls regplot() once per facet, but it owns the whole figure and returns a FacetGrid.
g = sns.lmplot(data=tips, x="total_bill", y="tip") print(type(g)) # seaborn.axisgrid.FacetGrid g.set_axis_labels("Total bill ($)", "Tip ($)")
<class 'seaborn.axisgrid.FacetGrid'>
lmplot is figure-level, so it can fit a separate line per hue group and split into panels with col. regplot can do neither.
The two functions share almost all their statistical arguments — ci, order, logistic, lowess, robust, x_estimator, x_bins, truncate. What differs is who controls the figure.
| Aspect | regplot() |
lmplot() |
|---|---|---|
| Level | Axes-level | Figure-level |
| Returns | A Matplotlib Axes |
A FacetGrid |
Accepts ax= |
Yes | No — raises TypeError |
Faceting with col= / row= |
No | Yes |
hue= |
Supported in v0.12+ (was absent before) | Yes, long-standing |
| Size control | Via plt.subplots(figsize=...) |
Via height= and aspect= |
| Legend | Standard Matplotlib legend | Placed outside the axes by the grid |
| Data input | Accepts vectors or a DataFrame | Requires long-form data= with column names |
| Use when | Putting a fit inside a bigger composed figure | Comparing the same fit across subgroups |
TypeError: lmplot() got an unexpected keyword argument 'ax'. Figure-level functions create their own figure by design; they cannot be inserted into one you already made. If you need a regression inside a subplot grid you built yourself, use regplot(). This mirrors the pattern you saw with displot/histplot and catplot/boxplot.
# figure-level sizing: height is per-facet, in inches; width = height * aspect sns.lmplot(data=tips, x="total_bill", y="tip", col="time", height=4, aspect=1.1)
Getting to the Axes underneath a FacetGrid
g = sns.lmplot(data=tips, x="total_bill", y="tip", col="time") g.axes # 2-D numpy array of Axes g.axes.flat[0] # first Axes object g.figure # the Matplotlib Figure (v0.11.2+; g.fig is the older name) for ax in g.axes.flat: ax.axhline(3, color="grey", linestyle="--", linewidth=1)
19.3 Model variants: order=, logistic=, lowess=, robust=, x_estimator=
The default is a straight OLS line, but the same functions can fit several other models. Only one of order, logistic, lowess and robust may be active at a time.
| Argument | Model | Needs statsmodels | Use it when |
|---|---|---|---|
order=2, order=3 |
Polynomial least squares | No (uses numpy) | The cloud is visibly curved and you want a parametric curve |
lowess=True |
Locally weighted scatterplot smoothing | Yes | You want the shape with no assumption at all about its form |
robust=True |
Robust linear model (Huber-type down-weighting) | Yes | A few extreme points are dragging the OLS line |
logistic=True |
Logistic regression | Yes | y is binary 0/1 |
x_estimator=np.mean |
Still linear, but the scatter is collapsed | No | x takes few discrete values |
x_bins=8 |
Still linear, x grouped into bins for display | No | Continuous x, very noisy y |
Polynomial fits with order=
# the second Anscombe dataset is a perfect parabola with noise-free points d2 = anscombe.query("dataset == 'II'") fig, axes = plt.subplots(1, 2, figsize=(10, 4)) sns.regplot(data=d2, x="x", y="y", ci=None, ax=axes[0]) sns.regplot(data=d2, x="x", y="y", order=2, ci=None, ax=axes[1]) axes[0].set_title("order=1 (default) — wrong shape") axes[1].set_title("order=2 — fits perfectly")
With a binary outcome, logistic=True fits an S-shaped probability curve instead of a straight line. The raw points can only sit at 0 or 1; the curve gives the probability between them.
A degree-k polynomial can pass exactly through k+1 points. With 11 rows, order=5 will produce a wiggly curve that looks impressively tight and predicts nothing. Stay at order=2 unless there is a physical reason for more, and prefer lowess=True when you only want to see the shape.
Non-parametric shape with lowess=True
sns.regplot(data=tips, x="total_bill", y="tip", lowess=True, line_kws={"color": "crimson"})
LOWESS fits many tiny local regressions and stitches them together. It assumes nothing about global form, which makes it the honest first look at a relationship. Two caveats: seaborn does not draw a confidence band for lowess fits (the ci argument is ignored), and the curve can wander at the edges where few points support it.
Down-weighting outliers with robust=True
# Anscombe III: one outlier drags the OLS line off the true trend d3 = anscombe.query("dataset == 'III'") fig, axes = plt.subplots(1, 2, figsize=(10, 4)) sns.regplot(data=d3, x="x", y="y", ci=None, ax=axes[0]) sns.regplot(data=d3, x="x", y="y", robust=True, ci=None, ax=axes[1]) axes[0].set_title("OLS — pulled up by one point") axes[1].set_title("robust=True — follows the other ten")
Robust fitting is noticeably slow, because the bootstrap for the band refits the robust model 1000 times. Pass ci=None or lower n_boot while exploring.
Binary outcomes with logistic=True
When y is 0/1, a straight line is meaningless — it predicts values below 0 and above 1. Logistic regression fits an S-shaped curve interpreted as P(y = 1 | x).
sns.lmplot( data=titanic, x="age", y="survived", logistic=True, ci=95, height=4.5, aspect=1.3, y_jitter=0.03, scatter_kws={"alpha": 0.3}, )
Two practical details:
ymust be numeric 0/1, not the strings "yes"/"no".titanic["survived"]is already 0/1; a column likealivewould need(titanic["alive"] == "yes").astype(int).y_jitter=(andx_jitter=) shift the plotted points slightly so that the two horizontal stripes of 0s and 1s show their density. Jitter affects display only — the fit uses the original values.
Each bootstrap replicate refits a logistic regression. On a few thousand rows this is fine; on 100,000 rows it will hang. Use ci=None, or reduce n_boot, or fit the model once in statsmodels and plot the predicted curve yourself.
Discrete or very noisy x: x_estimator and x_bins
When x has only a handful of distinct values — tips["size"] takes 1 to 6 — the scatter is a set of vertical stripes and the eye cannot judge the trend. x_estimator replaces each stripe by a summary point plus an error bar:
sns.regplot(data=tips, x="size", y="tip", x_estimator=np.mean)
The error bars here are bootstrap confidence intervals of the estimator (the mean), controlled by the same ci argument. Crucially, the regression is still fitted on the raw observations, not on the collapsed means — only the display changes.
For continuous but noisy x, x_bins chops x into groups first:
# 8 equal-count bins; the line is still fitted on all raw points sns.regplot(data=tips, x="total_bill", y="tip", x_bins=8) # or give explicit bin centres sns.regplot(data=tips, x="total_bill", y="tip", x_bins=[10, 20, 30, 40, 50])
Other useful arguments
ArgumentEffect
truncate=True (default)The line stops at the data range instead of spanning the whole axis
fit_reg=FalseDraw only the scatter — turns lmplot into a faceted scatterplot
scatter=FalseDraw only the line and band — useful for overlaying several fits
marker= / markers=Point shape; markers=["o", "x"] distinguishes hue levels in print
seed=Fixes the bootstrap RNG so the band is reproducible
logx=TrueFits y = a + b*log(x); x must be positive
Getting the actual numbers
None of the above returns a coefficient. Three standard ways to obtain them:
from scipy import stats import statsmodels.api as sm d = tips.dropna(subset=["total_bill", "tip"]) # 1. numpy — coefficients only, highest power first slope, intercept = np.polyfit(d["total_bill"], d["tip"], deg=1) # 2. scipy — slope, intercept, r, p-value, standard error res = stats.linregress(d["total_bill"], d["tip"]) print(res.slope, res.intercept, res.rvalue ** 2, res.pvalue) # 3. statsmodels — the full table you would put in a report X = sm.add_constant(d["total_bill"]) model = sm.OLS(d["tip"], X).fit() print(model.summary())
0.1050245173843534 0.9202696135546726 0.4566165863516758 6.692470646864404e-34
OLS Regression Results
==============================================================================
Dep. Variable: tip R-squared: 0.457
Model: OLS Adj. R-squared: 0.454
Method: Least Squares F-statistic: 203.4
Date: Wed, 30 Sep 2026 Prob (F-statistic): 6.69e-34
Time: 12:32:10 Log-Likelihood: -350.54
No. Observations: 244 AIC: 705.1
Df Residuals: 242 BIC: 712.1
Df Model: 1
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 0.9203 0.160 5.761 0.000 0.606 1.235
total_bill 0.1050 0.007 14.260 0.000 0.091 0.120
==============================================================================
Omnibus: 20.185 Durbin-Watson: 2.151
Prob(Omnibus): 0.000 Jarque-Bera (JB): 37.750
Skew: 0.443 Prob(JB): 6.35e-09
Kurtosis: 4.711 Cond. No. 53.0
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
Plot first with seaborn to check the shape and spot outliers, then fit with statsmodels to report the numbers, then plot the residuals of that fitted model to validate it. Seaborn opens and closes the loop; statsmodels does the middle.
Overlaying a model you fitted yourself
xs = np.linspace(d["total_bill"].min(), d["total_bill"].max(), 200) ys = intercept + slope * xs ax = sns.scatterplot(data=d, x="total_bill", y="tip", alpha=0.4) ax.plot(xs, ys, color="crimson", linewidth=2, label=f"y = {intercept:.2f} + {slope:.3f}x") ax.legend()
This is the pattern to use whenever the fit must be reproducible, reportable, or shared with the modelling code — the line on the chart is then provably the same line you wrote in the text.
19.4 residplot() and reading residuals
A residual is what the model missed: residual = observed y − predicted y. sns.residplot() fits the regression, subtracts it, and plots the leftovers against x.
sns.residplot(data=tips, x="total_bill", y="tip") plt.axhline(0, color="grey", linestyle="--")
Residuals should be a formless cloud around zero. A curve means the model shape is wrong; a widening funnel means the variance grows with x and the confidence intervals cannot be trusted.
It is axes-level (accepts ax=), and it accepts order=, robust= and lowess= — the last one adds a smoother through the residuals, which makes leftover structure much easier to see:
sns.residplot(data=tips, x="total_bill", y="tip", lowess=True, line_kws={"color": "crimson"})
On the original scatter, a strong trend dominates your eye and hides small systematic errors. Subtracting the trend removes that dominance and rescales the y-axis to the size of the errors, so a curvature of 0.3 units becomes obvious instead of invisible next to a 10-unit slope.
What a healthy residual plot looks like
A formless cloud, symmetric about zero, with roughly constant vertical spread all the way across:
+2 | . . . . | . . . . . . . . 0 |------.---.-.----.---.------ | . . . . . . . -2 | . . . . +--------------------------- x GOOD: no trend, even spread
Pattern 1 — curvature (a missing non-linear term)
+2 | . . | . . 0 |------.------.-------------- | . . -2 | . +--------------------------- x BAD: a smile/frown shape. The straight line is too low at the ends and too high in the middle. Try order=2, log(x), or sqrt(y).
Pattern 2 — funnel / cone (heteroscedasticity)
+4 | . | . . . 0 |---.--.--.----.-----.------- | . . . -4 | . . +--------------------------- x BAD: spread grows with x. The line may be unbiased, but the standard errors and the CI band are wrong. Try log(y), sqrt(y), or robust standard errors.
Pattern 3 — outliers
| +8 | * |
|---|---|
| 0 | --.-.--.-.--.--.-.--.-.----- |
| -8 | * |
| +--------------------------- x | |
| BAD: one or two points far from the | |
| band of the rest. Investigate them | |
| individually — data error, a | |
| different population, or a real | |
| but rare case? |
Pattern 4 — a trend in the residuals
If residuals slope upward or downward across x, the fit is not actually a least-squares fit of the data you think it is — typically a sign that you plotted residuals from one model against a different x, or that a group variable was ignored. For a correctly fitted simple OLS, the residuals are uncorrelated with x by construction, so a visible slope always means something is off in the setup.
Same rule as before: you get a picture, not an array. To have the numbers, fit with statsmodels and read model.resid, or compute y - (intercept + slope * x) yourself. The picture is for diagnosis; the array is for tests such as Breusch-Pagan or Durbin-Watson.
model = sm.OLS(d["tip"], sm.add_constant(d["total_bill"])).fit() resid = model.resid # a real pandas Series fitted = model.fittedvalues # residuals vs fitted values — the standard diagnostic, not vs x ax = sns.scatterplot(x=fitted, y=resid, alpha=0.4) ax.axhline(0, color="grey", linestyle="--") ax.set(xlabel="Fitted values", ylabel="Residuals")
19.5 Fitting separate lines per group
A single global line is often the wrong summary because the population is a mixture of groups. Seaborn gives you three ways to split.
hue= — several lines in one Axes
sns.lmplot(data=penguins, x="bill_length_mm", y="bill_depth_mm", hue="species", height=5, aspect=1.2)
This is the penguins version of Simpson's paradox and it is worth running once with and once without hue. Across all penguins pooled together, longer bills go with shallower bills — a clear negative slope. Within each species, longer bills go with deeper bills — three clear positive slopes. The pooled line describes a difference between species (Gentoo have long shallow bills), not a relationship inside any of them.
# run these two side by side and compare the slope signs fig, axes = plt.subplots(1, 2, figsize=(11, 4.5), sharex=True, sharey=True) sns.regplot(data=penguins, x="bill_length_mm", y="bill_depth_mm", ax=axes[0], scatter_kws={"alpha": 0.4}) axes[0].set_title("Pooled: negative slope") for sp, sub in penguins.groupby("species"): sns.regplot(data=sub, x="bill_length_mm", y="bill_depth_mm", ax=axes[1], label=sp, scatter_kws={"alpha": 0.4}) axes[1].legend(title="species") axes[1].set_title("By species: three positive slopes")
Averaging over groups is like measuring the average temperature of a hotel by mixing water from the hot and cold taps of every room. The mixture tells you nothing about whether any individual tap runs hot or cold — and it can even suggest the opposite of what every room does.
col= and row= — one panel per group
# one column per lunch/dinner sns.lmplot(data=tips, x="total_bill", y="tip", col="time", height=4) # a full grid: rows = smoker, columns = time, colours = sex sns.lmplot(data=tips, x="total_bill", y="tip", row="smoker", col="time", hue="sex", height=3.5, aspect=1.1, ci=None, scatter_kws={"alpha": 0.5})
Faceting is better than hue when the groups overlap heavily — separate panels avoid a plate of spaghetti. hue is better when you need to compare slopes directly, because comparing angles across panels is harder than comparing them side by side.
lmplot shares x and y across facets by default (sharex=True, sharey=True). Keep it that way. If each panel rescales itself, a gentle slope and a steep one can end up looking identical, which defeats the entire comparison.
col_wrap= for many groups, and the anscombe reproduction
This is the single most important figure in the chapter:
sns.lmplot( data=anscombe, x="x", y="y", col="dataset", hue="dataset", col_wrap=2, height=3.5, ci=None, palette="muted", scatter_kws={"s": 60, "alpha": 0.8}, )
Four datasets, eleven points each. The fitted line is essentially y = 3.00 + 0.500x in all four, the correlation is 0.816 in all four, and the means and variances match to two decimals. The shapes could not be more different:
| Dataset | What it actually is | What the summary statistics say |
|---|---|---|
| I | A genuine noisy linear relationship | slope 0.500, r 0.816 — correct here |
| II | A clean parabola, no noise at all | slope 0.500, r 0.816 — badly wrong |
| III | A perfect straight line plus one outlier | slope 0.500, r 0.816 — the outlier tilts the line |
| IV | x is constant except for one point | slope 0.500, r 0.816 — the slope is defined by a single observation |
# confirm the identical statistics numerically for name, sub in anscombe.groupby("dataset"): r = stats.linregress(sub["x"], sub["y"]) print(name, round(r.slope, 3), round(r.intercept, 3), round(r.rvalue, 3))
I 0.5 3.0 0.816 II 0.5 3.001 0.816 III 0.5 3.002 0.816 IV 0.5 3.002 0.817
19.6 What these plots cannot tell you
Being explicit about the boundaries of a tool is part of teaching it.
1. No fit statistics come back
There is no .slope, no .r_squared, no p-value anywhere on the returned object. Seaborn deliberately does not expose the fit. If a report needs a coefficient, it must come from scipy.stats.linregress, numpy.polyfit or statsmodels — and then you should verify that the plotted line and the reported line are the same model, because lowess, robust and order=2 all draw something OLS did not produce.
2. Extrapolation is meaningless
The default truncate=True stops the line at the observed range for exactly this reason. Setting truncate=False extends the line across the whole axis, which invites the reader to read off a prediction in a region where no data exists. The model has no evidence there, and linearity that holds between 5 and 50 dollars need not hold at 500.
3. Correlation is not causation
A fitted line describes association. It does not tell you whether x causes y, y causes x, or a third variable causes both. Larger bills come with larger tips — but party size drives both, and neither plot nor slope can separate those stories.
4. A straight line will be drawn through anything
Ask for regplot on a parabola and you get a straight line with a confidence band, drawn with complete composure. Nothing warns you. Anscombe II is the demonstration. The check is the residual plot, not the fit plot.
5. Confounders are invisible — Simpson's paradox
A pooled trend can be the exact reverse of every within-group trend, as the penguins bill example shows. The defence is procedural: whenever you have a candidate grouping variable, refit with hue= or col= and compare the slope signs before you write down a conclusion.
6. Only one predictor is shown
These are bivariate pictures. Real models have many predictors, and the marginal relationship between x and y can differ in sign from the partial relationship after controlling for other variables. For that you need added-variable (partial regression) plots from statsmodels, not regplot.
7. The band assumes the model
The bootstrap CI is computed under the assumption that the chosen form is correct. If the form is wrong, the band is a precise interval around a wrong line. Precision and correctness are independent properties.
8. Nothing checks whether a regression makes sense at all
Seaborn will happily regress on an ordinal code, a year stored as an integer, or a categorical variable that pandas has encoded numerically. The function is not a statistician.
19.7 Gotchas
Restated because it is the number-one misreading. About 95% of the time the true mean line lies inside that band; individual points are expected to fall well outside it, and roughly half of them do. Seaborn cannot draw a prediction interval at all — for that use statsmodels: model.get_prediction(X).summary_frame() gives both mean_ci_* and obs_ci_* columns.
Figure-level function. Use regplot() for subplots, lmplot() for faceting. Corollary: plt.title("...") after lmplot titles only the last facet — use g.figure.suptitle("...") or g.set_titles("{col_name}").
With 11 or 30 rows, a quadratic or cubic will chase noise and the band will be misleadingly tight in the middle and explosive at the edges. Prefer lowess=True for shape exploration, and only adopt a polynomial when the residual plot of the linear fit shows genuine curvature.
Without it you get ImportError: The lowess method requires statsmodels (similar messages for the others). Install with pip install statsmodels. And expect slowness: robust and logistic fits are re-run n_boot times for the band, so a plot that took 0.2 s with the default can take 30 s or more. Set ci=None while exploring.
sns.regplot(data=tips, x="day", y="tip") either fails or silently regresses on the internal codes of the category, in which case the "slope" depends entirely on alphabetical order. If x is categorical, you want boxplot, barplot or pointplot — see Chapter 12 — not a regression line. An ordered ordinal variable is a grey area: a line can be defensible, but state the coding you used.
penguins has missing measurements and titanic["age"] has 177 missing values. Seaborn drops any row where x or y is NaN without warning, so the fit may use far fewer rows than len(df). Always check: df[["x", "y"]].dropna().shape[0]. This matters most when comparing facets — one panel can be fitted on 30 rows and another on 300, and the plot gives no hint.
Combinations like order=2, lowess=True or logistic=True, robust=True are not supported. Seaborn raises an error or silently ignores one of them depending on the version. Choose one model per call; draw two calls into the same Axes if you want to compare.
It downloads from the seaborn-data GitHub repository and fails with a URL error offline. Workarounds: cache once with sns.load_dataset("tips").to_csv("tips.csv", index=False) and later pd.read_csv("tips.csv"); or clone github.com/mwaskom/seaborn-data and pass sns.load_dataset("tips", data_home="path/to/seaborn-data"); or set the SEABORN_DATA environment variable to that folder.
Two runs of the same code produce very slightly different band edges. For figures that must be byte-identical across runs, pass seed=0, or use ci=None.
x_jitter and y_jitter exist purely to reveal overlapping points. The regression is always computed on the original values. Do not describe a jittered plot as if the noise were in the data.
19.8 Practice set
# Setup for all exercises import numpy as np, pandas as pd, matplotlib.pyplot as plt, seaborn as sns from scipy import stats sns.set_theme(style="whitegrid") tips = sns.load_dataset("tips") penguins = sns.load_dataset("penguins") titanic = sns.load_dataset("titanic") anscombe = sns.load_dataset("anscombe") # 1. Basic fit and the band. # Draw regplot of tip vs total_bill. Then draw it again with ci=None and # with ci=99. Write one sentence describing what the band represents and # one sentence saying what it does NOT represent. # 2. Axes-level composition. # Build a 1x3 subplot figure with plt.subplots and place three regplots in # it: tip vs total_bill, tip vs size, total_bill vs size. Give each a title. # Then try to do the same with lmplot and record the exact error message. # 3. Reproduce the Anscombe quartet. # Use lmplot with col="dataset", col_wrap=2, ci=None. Then use # scipy.stats.linregress inside a groupby loop to print slope, intercept # and r for each of the four datasets. Explain in two sentences why # identical numbers with different pictures is the whole argument for # plotting your data. # 4. Choose the right model. # For anscombe dataset "II", plot order=1, order=2 and lowess=True side by # side. For dataset "III", plot default OLS next to robust=True. In each # case say which fit you would report and why. # 5. Residual diagnosis. # Draw residplot(tips, x="total_bill", y="tip", lowess=True). Decide which # of the four patterns from 19.4 you see. Then repeat with # y=np.log(tips["tip"]) and state whether the funnel improved. # 6. Logistic regression. # Plot survived vs age on titanic with logistic=True, y_jitter=0.03 and # hue="sex". Read off the approximate age at which predicted survival # probability crosses 0.5 for each sex. Then check how many rows were # actually used: titanic[["age","survived"]].dropna().shape. # 7. Simpson's paradox. # Fit bill_depth_mm on bill_length_mm for all penguins pooled, and record # the sign of the slope with linregress. Then fit it separately per species # and record the three signs. Draw both versions and write the one-sentence # conclusion a careless analyst would have written from the pooled fit. # 8. From picture to report. # Fit tip ~ total_bill with statsmodels OLS. Print the slope, its 95% # confidence interval, the p-value and R-squared. Then reproduce the # seaborn plot manually: scatterplot + your own line from the fitted # coefficients, with the equation in the legend. Confirm the two lines # coincide.
19.9 Interview questions
Q: What exactly does the shaded region around a seaborn regression line represent?
A 95% bootstrap confidence interval for the regression line itself — that is, for the estimated mean of y at each value of x. It is not a prediction interval for individual observations and it is not a standard deviation of the data. It narrows as the sample size grows, so a thin band means the position of the average line is well estimated, not that the model fits well or predicts individual points accurately.
Q: When would you use regplot and when lmplot?
regplot is axes-level: it takes ax=, draws into a figure you already control, and composes with plt.subplots or other seaborn axes-level calls. lmplot is figure-level: it creates and owns its figure, returns a FacetGrid, supports col= and row= faceting, and rejects ax=. Use regplot inside a composed layout; use lmplot when the point is to compare the same fit across subgroups.
Q: Can you get the slope, intercept or R² from these functions?
No. Seaborn fits the model internally — through numpy for polynomial fits and through statsmodels for lowess, robust and logistic — but returns only drawing objects. To report numbers, use scipy.stats.linregress (slope, intercept, r, p, stderr), numpy.polyfit (coefficients only), or statsmodels.api.OLS(...).fit() for a full summary with confidence intervals and diagnostics.
Q: Why is the Anscombe quartet the standard example in this chapter?
Four eleven-point datasets share nearly identical means, variances, correlation (0.816) and fitted line (roughly y = 3 + 0.5x), yet one is genuinely linear, one is a clean parabola, one is a straight line ruined by a single outlier, and one has a constant x apart from one leverage point. Summary statistics cannot distinguish them; a plot distinguishes them instantly. It is the strongest available argument for looking at data before trusting a fit.
Q: What does a good residual plot look like, and name two bad patterns.
Good: a structureless cloud centred on zero with roughly constant vertical spread across the whole range of x or of the fitted values. Bad pattern one is curvature — a smile or frown shape, meaning the linear form is missing a non-linear term; the fix is a polynomial term or a transformation of x or y. Bad pattern two is a funnel, where spread widens with x — heteroscedasticity; the coefficients may still be unbiased but the standard errors and the confidence band are untrustworthy, so use a variance-stabilising transform such as log(y) or robust standard errors.
Q: Your scatter has a few extreme points. What are your options in seaborn?
First look at them: they may be data-entry errors, a different population, or genuine rare cases, and the response differs in each case. Within seaborn, robust=True fits a robust linear model that down-weights them, and comparing it with the default OLS line shows how much influence they carry. It requires statsmodels and is slow because of the bootstrap, so pair it with ci=None during exploration. Never delete points merely because they move the line.
Q: How do you plot a regression when the outcome is binary?
Use logistic=True with a numeric 0/1 y — for example sns.lmplot(data=titanic, x="age", y="survived", logistic=True). This fits a logistic regression and draws the S-shaped predicted probability curve instead of a straight line, so predictions stay inside [0, 1]. Add y_jitter so the two stripes of 0s and 1s reveal their density, and remember that jitter is display-only. It needs statsmodels and is slow on large data, so consider ci=None.
Q: What is Simpson's paradox and how do you check for it here?
It is the situation where an association present in pooled data reverses or disappears once the data are split by a relevant group variable. In penguins, bill length and bill depth are negatively associated overall but positively associated within every species — the pooled line is really describing differences between species. The check is to refit with hue= or col= set to each candidate grouping variable and compare the slope signs before writing any conclusion.
Q: What can go wrong silently in these plots?
Three things in particular. Rows with NaN in x or y are dropped without any warning, so the fit may rest on far fewer observations than the DataFrame contains — and facets can end up with wildly different sample sizes. A categorical x may be regressed on its internal integer codes, making the slope an artefact of alphabetical order. And a straight line will be drawn confidently through curved data with no complaint at all, which is why the residual plot, not the fit plot, is the actual diagnostic.
19.10 Summary
regplot()is axes-level (acceptsax=, composes into your own figures);lmplot()is figure-level (returns aFacetGrid, supportscol=/row=, rejectsax=). They share every statistical argument.- The shaded band is a 95% bootstrap confidence interval for the regression line — not a prediction interval, not a standard deviation, and not evidence that the model is correct. Remove it with
ci=None. - Both functions fit a model and neither returns numbers. Use
scipy.stats.linregress,numpy.polyfitorstatsmodelswhenever a slope, p-value or R² has to be reported. - Model variants:
order=for polynomials,lowess=Truefor assumption-free shape,robust=Trueagainst outliers,logistic=Truefor binary y,x_estimator=np.meanandx_bins=for discrete or noisy x. Only one model per call; the last three need statsmodels and are slow. residplot()removes the trend so that leftover structure becomes visible. A good residual plot is a formless cloud around zero; curvature, funnels and isolated outliers each point to a specific fix.hue=,col=androw=fit separate lines per group. Always try this before accepting a pooled trend — Simpson's paradox is real andpenguinsdemonstrates it in two lines of code.- The Anscombe quartet: four datasets with the same slope, intercept and correlation but four different shapes. Plot before you trust a fit.
- Known limits: no fit statistics, meaningless extrapolation, no causal claim, only one predictor, and a band that assumes the model form is right.
- Watch for silent failures: dropped NaNs, categorical x regressed on codes, high-order polynomials on small samples, and
sns.load_dataset()requiring internet.
Next: Chapter 20 covers multi-plot grids — FacetGrid, pairplot and jointplot.
Questions about this chapter?
The course, homework review and the other five books are in Uzbek. Questions in English are welcome on Telegram.