4.3 Checking a regression · Visual Statistics

4.3.1 The four standard plots

A regression assumes four things about the noise ε: it has no pattern left (the model’s shape is right), it is equally large everywhere, it is roughly normal, and no single case runs the show. (Independence, the fifth, comes in Section 4.3.7.) Four plots check them at a glance, and every statistics program draws them. Here they are for the diabetes model of Chapter 4.2:

Code

influence = diabetes_fit.get_influence()
standardised = influence.resid_studentized_internal
leverage = influence.hat_matrix_diag
fitted_values = diabetes_fit.fittedvalues.to_numpy()


def add_trend(ax: Axes, x: np.ndarray, y: np.ndarray) -> None:
    """A lowess curve through a scatter plot."""
    smooth = sm.nonparametric.lowess(y, x, frac=0.6)
    ax.plot(smooth[:, 0], smooth[:, 1], color=style.CORAL, ls="--", lw=1.5)


def plot_fitted(ax: Axes) -> None:
    """Residuals against fitted values."""
    ax.scatter(fitted_values, diabetes_fit.resid, s=7, color=style.BLUE, alpha=0.5)
    ax.axhline(0, color=style.INK, lw=1)
    add_trend(ax, fitted_values, diabetes_fit.resid.to_numpy())
    ax.set(title="Residuals vs fitted", xlabel="fitted value", ylabel="residual")


def plot_qq_residuals(ax: Axes) -> None:
    """QQ plot of the standardised residuals against the standard normal."""
    ordered = np.sort(standardised)
    expected = stats.norm.ppf((np.arange(1, len(ordered) + 1) - 0.5) / len(ordered))
    ax.plot([-3, 3], [-3, 3], color=style.INK, ls="--", lw=1)
    ax.scatter(expected, ordered, s=7, color=style.BLUE, alpha=0.5)
    ax.set(title="Normal QQ", xlabel="expected if normal", ylabel="standardised residual")


def plot_scale_location(ax: Axes) -> None:
    """√|standardised residual| against fitted values."""
    root = np.sqrt(abs(standardised))
    ax.scatter(fitted_values, root, s=7, color=style.BLUE, alpha=0.5)
    add_trend(ax, fitted_values, root)
    ax.set(title="Scale-location", xlabel="fitted value", ylabel="√|standardised residual|")


def plot_leverage(ax: Axes) -> None:
    """Standardised residuals against leverage, with Cook's distance contours."""
    ax.scatter(leverage, standardised, s=7, color=style.BLUE, alpha=0.5)
    p = len(diabetes_fit.params)
    h = np.linspace(0.005, 0.45, 200)  # wide enough to show where the contours are
    for d in (0.5, 1):
        bound = np.sqrt(d * p * (1 - h) / h)
        ax.plot(h, bound, color=style.CORAL, ls="--", lw=1)
        ax.plot(h, -bound, color=style.CORAL, ls="--", lw=1)
    ax.axhline(0, color=style.INK, lw=0.8)
    ax.set(
        title="Residuals vs leverage",
        xlabel="leverage",
        ylabel="standardised residual",
        ylim=(-4, 4),
        xlim=(0, h.max()),
    )


fig, axes = style.subplots(4, cols=2, height=3, width=9.5)
for ax, plot in zip(axes, [plot_fitted, plot_qq_residuals, plot_scale_location, plot_leverage]):
    plot(ax)
fig.tight_layout()
plt.show()

Figure 4.3.1: The standard diagnostic plots for the regression of disease progression on nine measurements. Top left: residuals against fitted values, with a smooth trend (dashed); a flat band means the shape is fine. Top right: QQ plot of the standardised residuals; a straight line means roughly normal noise. Bottom left: scale-location, the square root of the absolute standardised residuals against fitted values; a flat trend means equal noise. Bottom right: standardised residuals against leverage, with Cook’s distance contours at 0.5 and 1 (dashed); points beyond them would be influential.

This model passes: no curve, roughly normal residuals (slightly short tails), even noise, and no point near the Cook’s contours. Three quantities appear in these plots:

The leverage $h_i$ of case i measures how far its predictors are from the average predictors. With one predictor, $h_i = 1/n + (x_i - \bar x)^2 / \sum_j (x_j - \bar x)^2$ ; the average leverage is always (p + 1)/n.
A standardised residual divides each residual by its own standard error, $r_i = e_i / (\hat\sigma \sqrt{1 - h_i})$ . High-leverage points pull the line towards themselves, so their raw residuals look too small; this corrects for it, and puts all residuals on one scale where values beyond ±3 are rare.
Cook’s distance combines both: how far all the fitted values move when case i is left out, $D_i = \frac{r_i^2}{p + 1} \cdot \frac{h_i}{1 - h_i}$ . Large residual and high leverage makes a large D.

4.3.2 Residuals against fitted values

A regression assumes three things about the noise $\varepsilon$ : it has no pattern left (the shape of the model is right), it is equally large everywhere (constant variance), and the residuals do not depend on each other. One plot checks the first two: residuals against the fitted values.

Code

N = 200


def plot_residuals(
    ax: Axes,
    fitted: np.ndarray,
    residuals: np.ndarray,
    title: str,
    color: str,
    trend: bool = True,
) -> None:
    """Residuals against fitted values, with zero and (optionally) a smooth trend line."""
    ax.scatter(fitted, residuals, s=9, color=color, alpha=0.5)
    ax.axhline(0, color=style.INK, lw=1)
    if trend:
        smooth = sm.nonparametric.lowess(residuals, fitted, frac=0.5)
        ax.plot(smooth[:, 0], smooth[:, 1], color=style.INK, lw=1.5, ls="--")
    ax.set(title=title, xlabel="fitted value", ylabel="residual")


def residuals_of_line(x: np.ndarray, y: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    """Fit a straight line; return its fitted values and residuals."""
    fit = sm.OLS(y, sm.add_constant(x)).fit()
    return fit.fittedvalues, fit.resid


x = rng.uniform(0, 10, N)
fig, axes = style.subplots(3, height=3, width=10.5)
plot_residuals(
    axes[0], *residuals_of_line(x, 2 + 0.5 * x + rng.normal(0, 1, N)), "Fine", style.BLUE
)
plot_residuals(
    axes[1],
    *residuals_of_line(x, 0.12 * (x - 5) ** 2 + 0.5 * x + rng.normal(0, 0.5, N)),
    "Curved",
    style.GOLD,
)
plot_residuals(
    axes[2],
    *residuals_of_line(engel.income.to_numpy(), engel.foodexp.to_numpy()),
    "Fanning out",
    style.CORAL,
    trend=False,
)
fig.tight_layout()
plt.show()

Figure 4.3.2: Residuals against fitted values for three regressions. Left: simulated data where a line is right: a structureless band. Middle: simulated data with a curve: the residuals bend, the line has the wrong shape. Right: food spending against income: the residuals fan out, the noise grows with income.

A curve in the residuals means the model’s shape is wrong: add a curved term, transform a variable, or use a more flexible model (Chapter 6.5). A fan means the noise changes size. The line itself is then still fine on average, but its standard errors and intervals are off: too narrow where the noise is large.

4.3.3 Is each effect straight? Partial residual plots

Residuals against fitted values show that the shape is wrong, not where. With several predictors, a partial residual plot (component plus residual plot) checks one predictor at a time: add the predictor’s own fitted contribution, $\hat\beta_j x_j$ , back onto the residuals and plot them against $x_j$ . If the predictor’s effect really is straight, the points scatter around the line; a curve in them shows the shape the model is missing.

Code

def plot_component(ax: Axes, name: str, label: str) -> None:
    """Residuals plus the predictor's own contribution, against the predictor."""
    x = diabetes[name].to_numpy()
    partial = diabetes_fit.resid.to_numpy() + diabetes_fit.params[name] * x
    ax.scatter(x, partial, s=7, color=style.BLUE, alpha=0.4)
    grid = np.array([x.min(), x.max()])
    ax.plot(
        grid,
        diabetes_fit.params[name] * grid
        - np.mean(diabetes_fit.params[name] * x)
        + partial.mean(),
        color=style.CORAL,
    )
    add_trend(ax, x, partial)
    ax.set(title=label, xlabel=label, ylabel="partial residual")


fig, axes = style.subplots(3, height=2.9, width=10.5)
for ax, (name, label) in zip(
    axes, (("bmi", "body mass index"), ("bp", "blood pressure"), ("s5", "s5 (triglycerides)"))
):
    plot_component(ax, name, label)
fig.tight_layout()
plt.show()

Figure 4.3.3: Partial residual plots for three predictors of the diabetes model. The red line is each predictor’s straight-line effect; the dashed curve is a smooth fit through the partial residuals. Blood pressure looks straight; body mass index bends slightly upward at its high end, and s5 flattens at its low end. Both departures are small, but worth checking with a curved term (Chapter 6.5).

4.3.4 Uneven noise: logs, weights and robust errors

Food spending fans out because households differ in proportion to their income: a rich household’s spending varies by more money, but by a similar percentage. That is the multiplying world of the log-normal (Chapter 2.1), and logs turn it back into adding.

Code

log_fit = smf.ols("np.log(foodexp) ~ np.log(income)", engel).fit()
TICKS = [300, 500, 1000, 2000, 4000]


def plot_log_scatter(ax: Axes) -> None:
    """Food against income, both on log scales, with the fitted line."""
    ax.scatter(engel.income, engel.foodexp, s=9, color=style.BLUE, alpha=0.5)
    grid = np.geomspace(engel.income.min(), engel.income.max(), 50)
    ax.plot(grid, np.exp(log_fit.predict(pd.DataFrame({"income": grid}))), color=style.CORAL)
    ax.set_xscale("log")
    ax.set_yscale("log")
    ax.set_xticks(TICKS, [f"{t:,}" for t in TICKS])
    ax.set_yticks(TICKS[:-1], [f"{t:,}" for t in TICKS[:-1]])
    ax.minorticks_off()
    style.note(ax, f"slope {log_fit.params.iloc[1]:.2f}", x=0.03)
    ax.set(
        title="Both on log scales",
        xlabel="income (log scale)",
        ylabel="food spending (log scale)",
    )


fig, axes = style.subplots(2, height=3.2, width=9.5)
plot_log_scatter(axes[0])
plot_residuals(axes[1], log_fit.fittedvalues, log_fit.resid, "Residuals", style.BLUE)
axes[1].set(xlabel="fitted value (log food spending)", ylabel="residual (log)")
fig.tight_layout()
plt.show()

Figure 4.3.4: Left: log food spending against log income: a straight line with even scatter. Right: its residuals have no fan.

On log scales the slope means something new: an elasticity. A slope of 0.86 says that 1% more income goes with 0.86% more food spending. Less than 1%, so the food share falls as income rises: Engel’s law again, now as a single number.

A fan can also be tested. The Breusch-Pagan test regresses the squared residuals on the predictors: if the predictors can explain how large the residuals are, the noise is not constant. Its statistic is n times the R² of that second regression; with constant noise it follows a chi-square distribution with one degree of freedom per predictor.

Code

raw_fit = smf.ols("foodexp ~ income", engel).fit()


def breusch_pagan(ax_left: Axes, ax_right: Axes, fit, x: np.ndarray, xlabel: str) -> None:
    """Squared residuals against the predictor, and n·R² against its chi-square yardstick."""
    squared = fit.resid.to_numpy() ** 2
    helper = sm.OLS(squared, sm.add_constant(x)).fit()
    ax_left.scatter(x, squared, s=8, color=style.BLUE, alpha=0.5)
    xs = np.array([x.min(), x.max()])
    ax_left.plot(xs, helper.params[0] + helper.params[1] * xs, color=style.CORAL, lw=1.6)
    ax_left.set(xlabel=xlabel, ylabel="squared residual")
    statistic = len(x) * helper.rsquared
    style.null_panel(ax_right, statistic, dist=stats.chi2(1), tail="right", xlabel="n·R²")
    ax_right.set_xlim(0, max(16, min(statistic * 1.1, 30)))
    if statistic > 30:  # too far out to draw: point to it instead
        ax_right.annotate(
            "observed value\nfar to the right →",
            xy=(0.97, 0.45),
            xycoords="axes fraction",
            ha="right",
            color=style.CORAL,
            fontsize=9,
        )


fig, axes = plt.subplots(2, 2, figsize=(10, 5.6))
breusch_pagan(axes[0, 0], axes[0, 1], raw_fit, engel.income.to_numpy(), "income")
breusch_pagan(axes[1, 0], axes[1, 1], log_fit, np.log(engel.income.to_numpy()), "log income")
axes[0, 0].set_title("Raw scale")
axes[1, 0].set_title("Log scales")
fig.tight_layout()
plt.show()

Figure 4.3.5: The Breusch-Pagan test, step by step, for food spending on income (top) and log food spending on log income (bottom). Left: each household’s squared residual against income, with the line that predicts it. For the raw data the squared residuals climb steeply with income; on the log scale the climb is small. Right: the statistic, n·R² of that line, against the chi-square distribution with 1 degree of freedom that constant noise would give. The raw model’s statistic is off the scale of the curve (p ≈ 10⁻²⁵); the log model’s lies beyond nearly all of it too (p = 0.0002), because with 235 households even a small leftover fan is detected.

For food spending on income it gives p ≈ 10⁻²⁵. For the log-log model it gives p = 0.0002: the fan is mostly gone, as the plot shows, but with 235 households the test still detects what is left. As with normality tests (Chapter 3.8), look at the size of the problem in the plot, not only at the p-value.

When logs do not fit the question, there are two other fixes:

Weighted least squares gives each point a weight inversely proportional to its noise variance, so the noisy points (here, rich households) count less. If the standard deviation grows in proportion to income, the weights are 1/income².
Robust standard errors (also called sandwich or heteroscedasticity-consistent standard errors) keep the ordinary line and estimate each coefficient’s uncertainty in a way that allows the noise to change size.

Code

ols = smf.ols("foodexp ~ income", engel).fit()
wls = smf.wls("foodexp ~ income", engel, weights=1 / engel.income**2).fit()
robust = ols.get_robustcov_results("HC3")
boot_slopes = [
    smf.ols("foodexp ~ income", engel.sample(len(engel), replace=True, random_state=rng))
    .fit()
    .params.income
    for _ in range(2_000)
]


def plot_two_lines(ax: Axes) -> None:
    """The data with the ordinary and the weighted line."""
    ax.scatter(engel.income, engel.foodexp, s=8, color=style.GREY, alpha=0.5)
    grid = pd.DataFrame({"income": [300, 5000]})
    ax.plot(
        grid.income,
        ols.predict(grid),
        color=style.CORAL,
        label=f"ordinary: slope {ols.params.income:.3f}",
    )
    ax.plot(
        grid.income,
        wls.predict(grid),
        color=style.BLUE,
        label=f"weighted: slope {wls.params.income:.3f}",
    )
    ax.legend(fontsize=8, loc="upper left")
    ax.set(title="Two lines", xlabel="household income", ylabel="food spending")


def plot_intervals(ax: Axes) -> None:
    """The slope's interval by four methods."""
    rows = [
        ("ordinary", *ols.conf_int().loc["income"]),
        ("robust (HC3)", *robust.conf_int()[1]),
        ("weighted", *wls.conf_int().loc["income"]),
        ("bootstrap", *np.percentile(boot_slopes, [2.5, 97.5])),
    ]
    for i, (name, low, high) in enumerate(rows):
        ax.plot([low, high], [i, i], color=style.CORAL if i == 0 else style.BLUE, lw=3)
    ax.set_yticks(range(len(rows)), [r[0] for r in rows])
    ax.invert_yaxis()
    ax.set(title="95% intervals for the slope", xlabel="food spending per unit of income")


fig, axes = style.subplots(2, height=3.2, width=9.5)
plot_two_lines(axes[0])
plot_intervals(axes[1])
fig.tight_layout()
plt.show()

Figure 4.3.6: Left: the ordinary least squares line (red) and the weighted least squares line (blue, weights 1/income²) for food spending; the weighted line is pulled less by the noisy rich households. Right: the slope’s 95% interval from ordinary standard errors, robust (HC3) standard errors, the weighted fit, and the case bootstrap of Chapter 3.3. The ordinary interval is too narrow for these fanning data; the robust and bootstrap intervals agree that the slope is much less certain. The weighted interval is narrow too, but for a good reason: it leans on the precisely measured low-income households.

4.3.5 Outliers, leverage and influence

Not all points pull equally on a line. A point far out along the x axis has leverage: like a long lever, a small push there tilts the whole line. A point with high leverage and a large residual has influence: removing it changes the fit.

Code

with_dc = smf.ols("murder ~ poverty", crime).fit()
without_dc = smf.ols("murder ~ poverty", crime.drop("District of Columbia")).fit()
cooks = with_dc.get_influence().cooks_distance[0]


def plot_lines(ax: Axes) -> None:
    """The states, the line with DC and the line without it."""
    ax.scatter(crime.poverty, crime.murder, s=14, color=style.BLUE, alpha=0.6)
    dc = crime.loc["District of Columbia"]
    ax.annotate(
        "DC", (dc.poverty, dc.murder), xytext=(6, -3), textcoords="offset points", fontsize=9
    )
    grid = pd.DataFrame({"poverty": [8, 23]})
    ax.plot(
        grid.poverty,
        with_dc.predict(grid),
        color=style.CORAL,
        label=f"with DC: slope {with_dc.params.poverty:.2f}",
    )
    ax.plot(
        grid.poverty,
        without_dc.predict(grid),
        color=style.INK,
        ls="--",
        label=f"without: slope {without_dc.params.poverty:.2f}",
    )
    ax.legend(fontsize=8, loc="upper left")
    ax.set(title="One point", xlabel="poverty (%)", ylabel="murders per 100,000")


def plot_cooks(ax: Axes) -> None:
    """Cook's distance for every state, in order."""
    ax.vlines(range(len(cooks)), 0, cooks, color=style.BLUE, lw=1.5)
    worst = int(np.argmax(cooks))
    ax.annotate(
        "DC", (worst, cooks[worst]), xytext=(6, -3), textcoords="offset points", fontsize=9
    )
    ax.set(
        title="Cook's distance", xlabel="state (alphabetical)", ylabel="influence on the fit"
    )


fig, axes = style.subplots(2, height=3.2, width=9.5)
plot_lines(axes[0])
plot_cooks(axes[1])
fig.tight_layout()
plt.show()

Figure 4.3.7: Murder rate against poverty in the 50 states and Washington DC. Left: the line with DC (red) and without it (dashed). Right: Cook’s distance, a measure of how much each point moves the fit; DC dwarfs every state.

DC is not a state but a city, with the murder rate of a city. One point out of 51 raises the slope by about a third. Influential points are not wrong by definition, and they should not be deleted to get a nicer answer. But they must be known, and the result should be reported with and without them. A Cook’s distance well above the others, or above 1, deserves a look.

A different answer to outliers is a method that does not let any single point dominate. Robust regression replaces the squared residuals of least squares with a loss that grows only linearly for large residuals (Huber’s loss), so far-off points get less weight. The fit is found by iteratively down-weighting the points with the largest residuals.

Code

huber = smf.rlm("murder ~ poverty", crime, M=sm.robust.norms.HuberT()).fit()


def plot_robust(ax: Axes) -> None:
    """Ordinary, robust, and without-DC lines."""
    ax.scatter(crime.poverty, crime.murder, s=14, color=style.GREY, alpha=0.7)
    grid = pd.DataFrame({"poverty": [8, 23]})
    ax.plot(
        grid.poverty,
        with_dc.predict(grid),
        color=style.CORAL,
        label=f"least squares: {with_dc.params.poverty:.2f}",
    )
    ax.plot(
        grid.poverty,
        huber.predict(grid),
        color=style.BLUE,
        label=f"robust (Huber): {huber.params.poverty:.2f}",
    )
    ax.plot(
        grid.poverty,
        without_dc.predict(grid),
        color=style.INK,
        ls="--",
        label=f"without DC: {without_dc.params.poverty:.2f}",
    )
    ax.legend(fontsize=8, loc="upper left")
    ax.set(title="Three lines", xlabel="poverty (%)", ylabel="murders per 100,000")


def plot_weights(ax: Axes) -> None:
    """Final weights of the robust fit, one bar per state."""
    weights = huber.weights.to_numpy()
    ax.vlines(range(len(weights)), 0, weights, color=style.BLUE, lw=1.5)
    worst = int(np.argmin(weights))
    ax.annotate(
        "DC", (worst, weights[worst]), xytext=(6, 4), textcoords="offset points", fontsize=9
    )
    ax.set(
        title="Robust weights",
        xlabel="state (alphabetical)",
        ylabel="weight in the fit",
        ylim=(0, 1.1),
    )


fig, axes = style.subplots(2, height=3.2, width=9.5)
plot_robust(axes[0])
plot_weights(axes[1])
fig.tight_layout()
plt.show()

Figure 4.3.8: Murder rate against poverty, all 51 points including DC. Ordinary least squares (red) is pulled up by DC; robust regression with Huber’s loss (blue) gives DC little weight and lands close to the line without DC (dashed). Right: the weight each state ends up with in the robust fit; DC gets almost none.

4.3.6 Are the residuals normal?

Tests and intervals for regression assume normal noise. This matters less than people think: by the central limit theorem, slopes are close to normal anyway once there are a few dozen points. The QQ plot of the standardised residuals (Section 4.3.1) shows whether the tails are wild enough to worry about; a Shapiro-Wilk test on the residuals (Chapter 3.8) adds little to it. With small samples and heavy tails, use the bootstrap (Chapter 3.3) instead. If the outcome is a count, a share or yes/no, the noise cannot be normal, and the models of Chapters 4.4 and 4.6 are the right ones.

4.3.7 Residuals in time

The last assumption is that residuals are independent. Data collected over time often break it: what is high today tends to be high tomorrow.

Code

weekly = co2.co2
years = (weekly.index - weekly.index[0]).days.to_numpy() / 365.25
trend = sm.OLS(weekly.to_numpy(), sm.add_constant(years)).fit()


def plot_trend(ax: Axes) -> None:
    """Weekly CO₂ with a straight line."""
    ax.plot(weekly.index, weekly, color=style.BLUE, lw=0.8)
    ax.plot(weekly.index, trend.fittedvalues, color=style.CORAL)
    ax.set(title="A line through time", xlabel="year", ylabel="CO₂ (ppm)")


def plot_time_residuals(ax: Axes) -> None:
    """Residuals in time order, with the correlation between neighbours."""
    ax.plot(weekly.index, trend.resid, color=style.GOLD, lw=0.8)
    ax.axhline(0, color=style.INK, lw=1)
    neighbours = np.corrcoef(trend.resid[:-1], trend.resid[1:])[0, 1]
    style.note(
        ax, f"correlation with\nthe week before: {neighbours:.2f}", x=0.35, color=style.INK
    )
    ax.set(title="Residuals in time", xlabel="year", ylabel="residual (ppm)")


fig, axes = style.subplots(2, height=3.2, width=9.5)
plot_trend(axes[0])
plot_time_residuals(axes[1])
fig.tight_layout()
plt.show()

Figure 4.3.9: A straight line through 44 years of weekly CO₂. Left: the data and the line. Right: the residuals in time order: a curve and a yearly wave. Neighbouring residuals are nearly equal, so the data hold far less independent information than 2,284 points suggest.

The Durbin-Watson statistic puts a number on it:

$\text{DW} = \frac{\sum_{t=2}^{n} (e_t - e_{t-1})^2}{\sum_{t=1}^{n} e_t^2} \approx 2(1 - r_1)$

where $r_1$ is the correlation of each residual with the one before. Independent residuals give about 2; strongly positively correlated ones give values near 0. For the CO₂ line it is 0.03.

Code

residual = trend.resid
dw = float(np.sum(np.diff(residual) ** 2) / np.sum(residual**2))
null_dw = np.array(
    [np.sum(np.diff(e) ** 2) / np.sum(e**2) for e in rng.normal(size=(2000, len(residual)))]
)

fig, axes = style.subplots(2, height=3.2, width=10)
axes[0].scatter(residual[:-1], residual[1:], s=3, color=style.GOLD, alpha=0.5)
axes[0].set(xlabel="residual the week before", ylabel="residual this week", title="Neighbours")
style.null_panel(axes[1], dw, samples=null_dw, tail="left", xlabel="Durbin-Watson statistic")
axes[1].set(title="If the residuals were independent", xlim=(0, 2.3))
fig.tight_layout()
plt.show()

Figure 4.3.10: The Durbin-Watson statistic. Left: each residual of the CO₂ line against the residual of the week before; they lie almost on a line (r₁ near 1). Right: DW for 2,000 simulated sets of independent residuals of the same length (grey), clustered tightly around 2. The CO₂ line’s DW of 0.03 is far below anything independent residuals produce.

Ordinary standard errors count every point as new information. With residuals this alike, they are much too small, and p-values too impressive. Time series need their own models (Chapter 9.1). The same goes for clusters: pupils in the same class, or repeated measurements of one person (Chapter 4.7).

4.3.8 In one picture

Check	Plot	If it fails
right shape	residuals against fitted; partial residual plots	add curved terms, transform, or a flexible model
constant noise	scale-location plot; Breusch-Pagan test	logs, weighted least squares, robust standard errors
influential points	residuals against leverage; Cook’s distance	report with and without them; robust regression
normal noise	QQ plot of standardised residuals	usually fine; bootstrap if small and heavy-tailed
independence	residuals in time or by group; Durbin-Watson	time series or mixed models

The residuals should look like noise: no curve, no fan, no order.
Logs turn multiplying effects into adding ones; the slope becomes an elasticity.
Know your influential points, and report what they do.

Next: regression for counts and yes/no outcomes.