Regression Residual Analysis: Formulas, Diagnostic Plots, Thresholds and Interpretation

import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.stats.diagnostic import (
    linear_reset, het_breuschpagan, het_white
)
from statsmodels.stats.stattools import durbin_watson, jarque_bera
from scipy.stats import shapiro

predictors = ["G1", "G2", "studytime", "failures", "absences"]
work = data[["G3", *predictors]].dropna().copy()

X = sm.add_constant(work[predictors], has_constant="add")
fit = sm.OLS(work["G3"], X).fit()
infl = fit.get_influence()
frame = infl.summary_frame()

audit = work.reset_index(names="source_index")
audit["fitted"] = fit.fittedvalues.to_numpy()
audit["residual"] = fit.resid.to_numpy()
audit["studentized_deleted"] = infl.resid_studentized_external
audit["leverage"] = infl.hat_matrix_diag
audit["cooks_d"] = infl.cooks_distance[0]
audit["dffits"] = infl.dffits[0]
audit["covratio"] = infl.cov_ratio

n = len(audit)
p = X.shape[1]
audit["flag_studentized"] = audit["studentized_deleted"].abs() > 3
audit["flag_leverage"] = audit["leverage"] > 2 * p / n
audit["flag_cook"] = audit["cooks_d"] > 4 / n
audit["flag_dffits"] = audit["dffits"].abs() > 2 * np.sqrt(p / n)
audit["flag_covratio"] = (
    (audit["covratio"] < 1 - 3 * p / n) |
    (audit["covratio"] > 1 + 3 * p / n)
)
flag_cols = [c for c in audit if c.startswith("flag_")]
audit["flag_any"] = audit[flag_cols].any(axis=1)

reset = linear_reset(fit, power=2, use_f=True)
bp = het_breuschpagan(fit.resid, fit.model.exog)
white = het_white(fit.resid, fit.model.exog)
dw = durbin_watson(fit.resid)
jb = jarque_bera(fit.resid)
sw = shapiro(fit.resid)

Keep preprocessing and case identifiers explicit. The exact field names returned by an influence summary may differ by statsmodels version, so verify each column before publication.

vars <- c("G3", "G1", "G2", "studytime", "failures", "absences")
work <- na.omit(data[vars])
fit <- lm(G3 ~ G1 + G2 + studytime + failures + absences, data = work)

n <- nrow(work)
p <- length(coef(fit))

audit <- transform(
  work,
  fitted = fitted(fit),
  residual = resid(fit),
  studentized_deleted = rstudent(fit),
  leverage = hatvalues(fit),
  cooks_d = cooks.distance(fit),
  dffits = dffits(fit),
  covratio = covratio(fit)
)

audit$flag_studentized <- abs(audit$studentized_deleted) > 3
audit$flag_leverage <- audit$leverage > 2 * p / n
audit$flag_cook <- audit$cooks_d > 4 / n
audit$flag_dffits <- abs(audit$dffits) > 2 * sqrt(p / n)
audit$flag_covratio <- audit$covratio < (1 - 3*p/n) |
                       audit$covratio > (1 + 3*p/n)

plot(fit, which = 1)
qqnorm(resid(fit)); qqline(resid(fit))
hist(resid(fit), main = "Residual histogram", xlab = "Residual")

Additional packages can provide RESET and heteroscedasticity tests. Record package versions and do not assume that identically named residuals are calculated with identical small-sample conventions.

REGRESSION
  /DEPENDENT G3
  /METHOD=ENTER G1 G2 studytime failures absences
  /STATISTICS COEFF OUTS R ANOVA CI(95)
  /RESIDUALS DURBIN HISTOGRAM(ZRESID) NORMPROB(ZRESID)
  /CASEWISE PLOT(ZRESID) OUTLIERS(3)
  /SAVE PRED RESID ZRESID SRESID LEVER COOK.

Saved-variable names can vary. Inspect the active dataset after execution, retain the original case identifier and export the output viewer as the audit record.

Predicted G3:
=$B$2+$B$3*[@G1]+$B$4*[@G2]+$B$5*[@studytime]
+$B$6*[@failures]+$B$7*[@absences]

Raw residual:
=[@G3]-[@[Predicted G3]]

Studentized flag:
=ABS([@[Deleted Studentized Residual]])>3

Leverage flag:
=[@Leverage]>(2*6/649)

Cook flag:
=[@[Cook Distance]]>(4/649)

DFFITS flag:
=ABS([@DFFITS])>(2*SQRT(6/649))

COVRATIO flag:
=OR([@COVRATIO]<(1-3*6/649),[@COVRATIO]>(1+3*6/649))

Any flag:
=OR([@[Studentized Flag]],[@[Leverage Flag]],[@[Cook Flag]],
    [@[DFFITS Flag]],[@[COVRATIO Flag]])

Use full-precision diagnostic values in hidden or calculation columns and format the visible cells separately. Rounding before applying thresholds can change borderline flags.

A raw residual is observed G3 minus fitted G3. Positive values indicate underprediction and negative values indicate overprediction. The residual retains the outcome unit, making it easy to explain but difficult to compare across observations with different leverage.

In an OLS model with an intercept, residuals sum to approximately zero and are orthogonal to the included predictor columns. Those algebraic facts do not prove that the residual pattern is random or that the model is correctly specified.

Caution: Do not call a negative residual a negative outcome effect; it is simply an overprediction under the fitted equation.

Internally standardized residuals use the model’s overall residual spread, including the focal case. Externally studentized residuals estimate the spread after deleting that case. The deleted version can reveal a severe outlier that partially masks itself by inflating the full-sample error estimate.

The worked Regression Residual Analysis uses |externally studentized residual| > 3 as a review rule. Ten cases cross it, and the maximum absolute value is 7.581732.

Caution: A threshold crossing requests verification; it does not establish fabrication, miscoding or exclusion eligibility.

The hat matrix H maps observed outcomes to fitted values. Its diagonal hᵢᵢ measures how far a predictor profile lies from the center of the design space. Leverage depends on X, not on the observed G3 value.

With p = 6 fitted coefficients and n = 649, the worked review line is 2p/n = 0.01848998. Forty-seven observations exceed it, and the maximum is 0.10328893.

Caution: Leverage is not an error measure. A perfectly recorded case can have high leverage because it represents a rare but valid predictor profile.

Cook’s distance summarizes the joint effect of deleting one case on the fitted coefficient vector. It increases when residual magnitude, leverage or both are substantial. The conventional 4/n line equals 0.00616333 in this analysis.

Twenty-nine cases exceed the line and the maximum Cook’s distance is 0.20287792. The line is intentionally sensitive in a large sample and should be treated as a screen.

Caution: Do not delete all observations above 4/n. That practice can remove a legitimate part of the population and manufacture apparent assumption compliance.

DFFITS measures the standardized change in an observation’s fitted value when that observation is deleted. It therefore focuses on prediction at the focal design point rather than the complete coefficient vector.

The worked threshold is 2√(p/n) = 0.19230177. Twenty-nine cases exceed it and the maximum absolute DFFITS is 1.15073682.

Caution: DFFITS and Cook’s distance are related but not interchangeable; report the quantity that matches the influence question.

DFBETAS quantify how much each coefficient changes, in standard-error units, when a case is omitted. They answer which specific slope is sensitive, whereas Cook’s distance summarizes joint influence.

A case can have moderate Cook’s distance yet strongly affect one coefficient if its predictor profile aligns with that term. DFBETAS are especially valuable when the research conclusion focuses on G2, failures or another named predictor.

Caution: Because no verified DFBETAS table was supplied in the attached results, do not invent case counts or maxima for this measure.

COVRATIO compares the determinant of the coefficient covariance matrix with and without a case. Values below one can indicate reduced precision after deletion; values above one can indicate increased precision, depending on the design and residual contribution.

The worked interval is 0.97226502 to 1.02773498, and 49 cases lie outside it. This is the largest individual flag count in the supplied analysis.

Caution: A COVRATIO flag does not necessarily imply a large change in prediction or coefficient direction.

A PRESS residual predicts each observation from a model fitted without that observation. It is closely related to leverage-adjusted residuals and provides a leave-one-out view of prediction error. PRESS can reveal optimistic fitted-sample performance that ordinary residuals conceal.

The current article reports fitted RMSE and case-deletion diagnostics but does not claim a verified PRESS statistic. A future prediction-focused extension should calculate PRESS or nested cross-validation without changing the primary diagnostic findings.

Caution: Do not describe R² = 0.850563 as future predictive accuracy.

This graph is the broad first check for curvature, changing spread and unusual errors. A smooth curve around zero suggests missing nonlinear mean structure. A funnel suggests heteroscedasticity. Integer-valued outcomes can produce bands without automatically violating the model.

In the worked analysis, the plot is interpreted with RESET p = 0.001759 and two highly significant variance tests.

Caution: Do not fit a complicated curve solely because a smoother is not perfectly horizontal in a large sample.

The Q–Q plot emphasizes quantile departure, while the histogram emphasizes overall shape. They are complementary because a histogram can hide tail behavior through bin choice and a Q–Q plot can make the center appear satisfactory while extremes depart sharply.

Shapiro–Wilk W = 0.756572 and Jarque–Bera = 10019.316610 agree with the strong lower-tail pattern.

Caution: Normality tests become highly sensitive as n grows; interpret effect size, tail severity and inferential purpose, not p-values alone.

The scale–location graph usually plots √|standardized residual| against fitted values. Removing the residual sign prevents positive and negative errors from canceling visually and makes variance trends easier to detect.

The worked plot is supported by Breusch–Pagan p = 1.391 × 10⁻⁵ and White p = 6.832 × 10⁻⁷.

Caution: A scale-location trend identifies a variance problem; it does not estimate correct weights by itself.

RESET adds powers or functions of the fitted values to test whether the original linear predictor may omit nonlinear structure. The worked F statistic is 9.868464 with p = 0.001759.

A significant result indicates misspecification risk but does not identify the omitted variable, interaction or transformation. It can also react to heteroscedasticity or other departures.

Caution: Do not add arbitrary polynomial terms until RESET becomes nonsignificant.

The Breusch–Pagan test examines whether squared residual variation is associated with the supplied design. The statistic is 30.129021 with p = 1.391 × 10⁻⁵, providing evidence against constant variance.

The test is most interpretable when the auxiliary variance relationship is plausibly linear in the predictors.

Caution: A significant Breusch–Pagan test does not mean the OLS point estimates should automatically be discarded.

The White test allows a broader auxiliary variance pattern, often including squares and cross-products of regressors. The worked statistic is 66.454172 with p = 6.832 × 10⁻⁷.

Agreement with Breusch–Pagan strengthens the conclusion that residual variance is not constant, although the White test can consume many degrees of freedom in smaller datasets.

Caution: Do not infer the exact variance function from the White test alone.

Durbin–Watson = 1.857347 is reasonably close to two for the retained row order. It primarily targets first-order residual autocorrelation and is meaningful only when adjacent rows represent a relevant sequence.

Cross-sectional educational data may be clustered by school or classroom even when arbitrary file order produces an acceptable Durbin–Watson value.

Caution: Do not write that independence is proven because Durbin–Watson is near two.

G3 has a bounded practical range, while a linear model can predict any real number. Near a lower or upper boundary, residual variance can compress and the conditional mean can become nonlinear. A cluster of zero grades can create a pronounced negative tail when predictor histories imply higher expected performance.

Caution: Do not choose a new model solely to remove difficult observations; the alternative must match the outcome scale and estimand.

HC0 through HC3 covariance estimators leave OLS coefficients and fitted values unchanged but adjust coefficient uncertainty for nonconstant variance. HC3 gives stronger leverage correction and is often preferred for finite-sample sensitivity reporting.

The attached results establish heteroscedasticity but do not provide verified HC3 intervals, so this article recommends the comparison without inventing numerical robust results.

Caution: Robust standard errors do not repair nonlinearity, omitted variables, clustering or extreme prediction error.

Weighted least squares changes the estimation criterion so observations with larger modeled variance receive less weight. It can improve efficiency when a defensible variance function is known or estimated reliably.

Unlike robust covariance, WLS changes coefficients, fitted values and residuals. A complete Regression Residual Analysis must therefore be rerun for the weighted model.

Caution: Data-dependent weights can overfit the same residual pattern they are intended to correct.

The significant RESET result motivates a structured search for mean-form improvements. Transformations may linearize a relationship; polynomial terms represent smooth curvature; interactions allow one predictor’s association to vary with another; splines provide flexible local shape.

Caution: Repeatedly testing many forms until diagnostics look favorable inflates researcher degrees of freedom.

Robust regression changes the loss function so extreme residuals have less influence. It can reveal whether the central coefficient pattern persists when tail observations receive reduced weight.

Robust regression estimates a different fitted solution from OLS and should be presented as a sensitivity or alternative estimand, not as a hidden cleanup step.

Caution: Downweighting a valid rare subgroup can conceal meaningful heterogeneity.

Observations nested within classrooms, schools, families or repeated students can share unmodeled variation. Ordinary residual plots may look acceptable while standard errors remain too small because independence is violated at the cluster level.

Caution: Durbin–Watson does not diagnose arbitrary cluster dependence.

With 649 observations, small departures can become statistically detectable. At the same time, the very large Jarque–Bera value, ten extreme studentized residuals and agreement between two heteroscedasticity tests indicate that the current concerns are not based on one marginal p-value.

Caution: Do not treat every p < .05 as requiring a different model, and do not dismiss every significant result as mere sample-size sensitivity.

A sensitivity analysis asks whether conclusions change under a documented alternative, such as correcting a confirmed data error or temporarily omitting a small set of the most influential verified observations. It is not a license to select the cleanest result.

The primary model in this article retains all 649 complete cases. The 66 flagged observations form a review list, not a deletion list.

Caution: Never remove all flagged cases as a single automatic batch.

Residual diagnostics evaluate the fitted model on the analyzed sample. Cross-validation and external validation evaluate performance on observations not used to estimate the same coefficients. A model can have attractive residual plots yet generalize poorly, or have imperfect residual normality while predicting well.

Caution: Do not use the current fitted R² or RMSE as evidence of out-of-sample performance.

Software can differ in definitions of standardized residuals, deleted residuals, leverage-adjusted scales, case indexing and small-sample test formulas. Agreement requires matching the equation, rows, intercept, precision and definitions.

Caution: Do not average conflicting outputs or assume one program is correct because it is more familiar.

A publishable Regression Residual Analysis should preserve the raw-data identifier, analysis formula, missing-data rule, software version, coefficient table, case-level diagnostic ledger, threshold formulas, plots, formal tests and every sensitivity decision. This record protects against selective interpretation and makes future correction possible.

Caution: Do not overwrite the original output with a revised model or cleaned dataset.