6.2 Residuals, Heteroscedasticity, and Influential Points: Diagnosing Which Failure Mode to Address

A successful fit means the software returned a result. A residual arc may suggest missing nonlinear structure; a funnel may suggest changing variance; a high-leverage point may influence the fitted line. These patterns are clues to investigate with the data-generation process, not instructions to delete inconvenient rows.

Learning Objectives

• Connect linear, independent, homoscedasticity, and distributional assumptions to specific conclusions;
• Explain residual-vs-fitted plots, Q-Q plots, leverage, and Cook's distance;
• Apply appropriate boundaries when using robust or clustered standard errors;
• Understand how multicollinearity affects coefficient stability rather than necessarily harming predictive performance.

1. Distinguish Prediction Conditions from Inference Conditions

Common assumptions:

1. The conditional mean form is correct, $E[\epsilon \mid X]=0$;
2. The observed or error dependence structure is properly addressed;
3. Whether constant conditional variance is required by the chosen method;
4. The error distribution is compatible with finite-sample inference methods;
5. The design matrix lacks complete multicollinearity;
6. The distribution of the training data matches that of the application setting.

Heteroskedasticity does not necessarily bias OLS coefficient estimates, but it can lead to incorrect standard errors, and OLS is no longer guaranteed to be the most efficient linear unbiased estimator. On the other hand, omitted confounding variables that cause $E[\epsilon \mid X]\ne 0$ can directly introduce bias into the coefficient estimates.

2. Residuals vs Fitted: Checking Mean and Variance Structure

Run the blocks in order with statsmodels 0.15, pandas, NumPy, seaborn and Matplotlib. The synthetic data deliberately contain a shared team effect, changing noise scale and two missing fields. The fit cohort is selected explicitly before fitting so every later diagnostic refers to the same rows.

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
import matplotlib.pyplot as plt
import seaborn as sns

rng = np.random.default_rng(12)
training = pd.DataFrame({
    "team_id": np.repeat(np.arange(20), 4),
    "team_size": rng.integers(2, 9, 80),
    "resources_used": rng.uniform(10, 50, 80),
    "mission_type": pd.Categorical(rng.choice(["survey", "repair"], 80)),
})
team_effect = rng.normal(0, 3, 20)
training["duration_minutes"] = (
    30 - 2 * training["team_size"] + .4 * training["resources_used"]
    + team_effect[training["team_id"]]
    + rng.normal(0, 1 + .04 * training["resources_used"], 80)
)
training.loc[2, "resources_used"] = np.nan
training.loc[5, "mission_type"] = np.nan
required = ["duration_minutes", "team_size", "resources_used", "mission_type", "team_id"]
included = training[required].notna().all(axis=1)
fit_rows = training.loc[included].copy()
model = smf.ols(
    "duration_minutes ~ team_size + resources_used + C(mission_type)",
    data=fit_rows, missing="raise",
).fit()
fitted, residuals = model.fittedvalues, model.resid
fig, ax = plt.subplots()
sns.scatterplot(x=fitted, y=residuals, alpha=0.5, ax=ax)
ax.axhline(0, color="black", linestyle="--")
ax.set_xlabel("Fitted duration (minutes)")
ax.set_ylabel("Residual (minutes)")
plt.show()

Observations:

• Arc-shaped patterns: indicate missing nonlinearity or interaction effects;
• Funnel shape: suggests variance changes with the mean;
• Grouped bands: may point to categorical variables, measurement discretization, or omitted subgroups;
• Large structured regions: could reflect temporal or spatial dependence or model specification errors.

Smooth lines can help identify systemic patterns, but residuals and fitted values share estimated parameters. Therefore, graphical interpretation must be contextualized with knowledge of the data generation process.

3. Q-Q Plots Are Not the Final Switch for Model Validity

A Q-Q plot compares the quantiles of standardized residuals against those of a reference distribution. Deviations in the tails may arise from heavy-tailed distributions, outliers, mixture populations, or incorrect variance specification.

import statsmodels.api as sm
influence = model.get_influence()
sm.qqplot(influence.resid_studentized_internal, line="45")
plt.show()

Normality of errors is primarily relevant for certain small-sample exact t/F inference; ordinary least squares (OLS) point estimation does not require either the response variable $Y$ or the residuals to be perfectly normal. Large-sample approximations also depend on independence, moment conditions, and design structure, sample size alone is not sufficient to justify ignoring clustering.

If the goal is prediction, the performance of the validation set and interval coverage are more directly relevant than whether residuals appear normally distributed.

4. How to Handle Heteroskedasticity Depends on the Objective

• For inference only: use heteroskedasticity-consistent covariance (HC series);
• When the variance structure is known: apply weighted least squares;
• When variance changes with the mean: transform the response or use an appropriate GLM;
• For prediction intervals: explicitly model the conditional variance or use quantile regression.

robust_model = model.get_robustcov_results(cov_type="HC3")
print(robust_model.summary())

Robust standard errors do not correct for nonlinearity, omitted variables, incorrect dependence structures, or bad data. They only adjust the covariance estimates and leave the original OLS fitted values unchanged.

HC3 inflates each squared residual contribution using $(1-h_{ii}{)}^{-2}$, where $h_{ii}$ is leverage. This addresses a finite-sample weakness of basic heteroskedastic covariance, but does not make correlated rows independent. High leverage near 1 is also a warning that the correction can be unstable.

5. Clustering and Time Dependence

Tasks within a team may share unobserved factors. Positive within-team dependence can make independence-based intervals too narrow. When treatment assignment or sampling occurs by team, use an inference method that respects that unit, such as clustered covariance or an appropriate hierarchical model:

clustered = model.get_robustcov_results(
    cov_type="cluster", groups=fit_rows["team_id"].to_numpy(),
    use_correction=True, df_correction=True, use_t=True,
)
print(pd.DataFrame({
    "ordinary_se": model.bse,
    "hc3_se": robust_model.bse,
    "cluster_se": clustered.bse,
}))

When the number of clusters is small, conventional asymptotic approximations may be unreliable and require small-sample corrections, randomization inference, or other specialized methods.

The group vector uses fit_rows, not the unfiltered training table: its length and order must match the fitted design matrix. This fixture has 20 groups, which does not by itself establish reliable coverage. The selected finite-sample corrections and t reference distribution are explicit; statsmodels covariance options document their behavior.

Temporal error components may also exhibit autocorrelation; methods like Newey–West/HAC rely on bandwidth selection and stationarity assumptions and cannot be treated as a one-size-fits-all fix.

6. Leverage, Residuals, and Influence Are Different

• Leverage: Whether a point is far from others in the $X$ space;
• Residual: How much the observed value differs from the fitted value;
• Influence: How the fit changes when that point is removed.

Cook's distance combines residual and leverage into a single diagnostic metric, it is a screening tool, not an automatic rule for deletion. 4/n serves as an empirical threshold, not a rigid rule that declares any point "wrong" upon exceeding it.

diagnostics = pd.DataFrame({
    "leverage": influence.hat_matrix_diag,
    "studentized_residual": influence.resid_studentized_internal,
    "cooks_d": influence.cooks_distance[0],
}, index=fit_rows.index)
print(diagnostics.sort_values("cooks_d", ascending=False).head())

Leverages are the diagonal of $H=X(X^{T}X{)}^{-1}{X}^T$ for a full-rank design. Their sum is the number of fitted coefficients, including the intercept. A remote predictor row can have small influence if its outcome lies close to the fitted surface.

Handling high-influence points:

1. Verify data collection and parsing processes;
2. Determine whether the point represents a valid member of the target population;
3. Assess the sensitivity of results to inclusion or exclusion;
4. Consider using robust regression or a more appropriate model;
5. Document the rules and how conclusions change.

Removing genuine outliers can cause the model to only apply to typical, non-extreme cases. It is essential to clearly define the new scope of applicability.

7. Collinearity Increases Coefficient Uncertainty

When columns are nearly linearly dependent, the model struggles to isolate individual contributions. This manifests as coefficient sensitivity to small data changes, large standard errors, and potential sign flips.

VIF:

$$VI{F}_j=\frac{1}{1-R_j^2},$$

where $R_j^2$ is derived from using all other features to predict the $j$th column. VIF > 5 or >10 are merely empirical thresholds and cannot replace judgment based on the specific context.

Collinear features may still jointly provide good predictive performance; the primary issue lies in the interpretability of individual coefficients and the stability of extrapolation. Removing features, creating composite indicators, redesigning sampling, or applying regularization each carries distinct semantic consequences.

8. Missing Patterns and Selection Also Require Diagnosis

Regression libraries often silently remove rows with missing data using the complete cases approach. First, compare the cases entering the model with those that are excluded:

The inclusion rule must cover every variable used by the formula, including categorical predictors and any transformation domain restrictions. An additional cluster-ID requirement also changes the analysis cohort. missing="raise" prevents a second silent deletion step after the explicit selection.

print({"input_rows": len(training), "fit_rows": len(fit_rows), "nobs": model.nobs})
print(training.groupby(included)["mission_type"].value_counts(dropna=False, normalize=True))
assert fit_rows.index.equals(model.model.data.row_labels)

If the excluded cases are related to the outcome or to the population characteristics, the target population of the complete-case model has been altered. You must return to examining the missing data mechanism and the estimand, rather than the final nobs alone.

9. Conduct Sensitivity Analysis After Diagnosis

Compare:

• Linear versus spline or nonlinear specifications;
• Ordinary, robust, and clustered standard errors;
• Models with and without influential outliers;
• Different approaches to handling missing data;
• Pre-specified combinations of covariates;
• Time or subgroup slicing.

If the direction of conclusions changes dramatically with reasonable specification choices, report this sensitivity, do not simply select the specification that yields the most attractive results.

Common Misconceptions

• Non-normal residuals rule out OLS: It’s the target variable and the inference method that matter, not the residuals alone.
• Robust standard errors fix model specification errors: They do not correct for misspecified mean structures or confounding.
• Cook’s D exceeding a threshold requires deletion: It only signals influential observations that warrant further investigation, not a justification for removal.
• Multicollinearity implies poor prediction: It primarily undermines the interpretability of coefficients, not the predictive accuracy of the model.

Exercise

1. Construct nonlinear, heteroskedastic, and high-leverage data sets, and compare residual plots.
2. Compare ordinary, HC3, and cluster-standard error estimates for the same model.
3. Remove observations with high Cook's D values, compare coefficient estimates before and after, and document your findings.
4. Assess whether complete-case deletion alters the population composition.

Summary

Diagnosis is not about assigning a simple "pass/fail" label to a model; it's about identifying which specific conclusions are at risk. Issues like mean misspecification, heteroscedasticity, correlated errors, influential points, and multicollinearity require distinct corrective actions and sensitivity reporting.

The next lesson integrates prediction workflows with regularization into cross-validation, ensuring that scaling, missing data handling, and hyperparameter selection do not inadvertently peek at the final test data.