Regression Diagnostics and Model Checking

In plain language

After fitting a statistical model, regression diagnostics are the checks you run to make sure the model's assumptions are reasonable — for example, that the relationship really is roughly straight-line, that errors are not ballooning at one end, and that no single patient is dominating the result. Think of them as the quality-control step between "model ran" and "model is trustworthy enough to report." In large healthcare datasets (hundreds of thousands of patients), the same checks still apply, but the rules change slightly: formal statistical tests almost always reject even tiny problems, so analysts rely more on visual plots and ask "how big is the problem?" rather than "is the p-value below 0.05?"

What regression diagnostics are and why they matter

Regression diagnostics are the audit trail of a fitted model. After specifying and estimating a linear or generalized linear model, the analyst must ask: do the data actually support the assumptions built into the estimation procedure? An OLS linear model assumes the conditional mean of the outcome is linear in the predictors, that residuals are homoscedastic (constant variance), and that errors are independent across observations. A logistic regression model assumes the log-odds is linear in the predictors and that the fitted probabilities are calibrated to observed event rates. None of these assumptions is verified before fitting; they are checked after, using the residuals and influence measures the fitted model produces.

The critical distinction — one that separates rigorous HEOR analysis from data dredging — is that diagnostics are run on a prespecified model to assess its adequacy, not to iteratively modify the model until diagnostics improve. Running diagnostics, tweaking predictors, and re-running until the residual plot looks acceptable is specification search (sometimes called torturing the model until it confesses). It inflates type-I error and produces optimistically biased standard errors, because the model has been tuned to one dataset's noise patterns. Preregistration or a prespecified SAP is the antidote; diagnostics then serve as transparency about residual assumption violations rather than as inputs to model revision.

The diagnostics toolkit by model failure mode

Linearity — assessed with a residual-versus-fitted-value plot. If the conditional mean of the outcome is truly linear in the predictors, residuals should scatter randomly around zero across the range of fitted values, with no systematic curve. A U-shape or inverted-U-shape in the residual plot signals a missed quadratic term; a monotone trend signals a missing predictor or a need for log-transformation of the outcome or a predictor. The fix is model re-specification (add a squared term, a spline, or a log transformation), not a change of estimator. Restricted cubic splines (RCS) are the modern flexible response-surface tool: they add nonlinearity locally without imposing a global polynomial shape.

Heteroscedasticity — the pattern in the residual-versus-fitted plot is a funnel or fan shape: residuals grow larger as fitted values increase. This is endemic in health data — costs, utilization, and length of stay all fan out with disease severity. Heteroscedasticity does NOT bias the OLS coefficient estimates but DOES bias the conventional standard errors downward, producing confidence intervals that are too narrow. Two remedies exist: (1) use heteroscedasticity-robust (HC3) standard errors, which adjust the variance estimator without changing the coefficients; (2) re-specify the outcome distribution using a GLM with a variance function that matches the data (gamma or Tweedie for costs; Poisson or negative binomial for counts). In most RWE settings, both are applied: the GLM as the primary model and the OLS with robust SEs as a sensitivity check.

Non-normality of residuals — assessed with a quantile-quantile (Q-Q) plot of standardized residuals against the theoretical normal quantiles. Heavy tails appear as an S-curve; skewness appears as a curved departure at one end. At the sample sizes typical of claims and EHR research (n ≥ 5,000), the Central Limit Theorem protects coefficient estimates and their standard errors from non-normality of residuals — the asymptotic Gaussian approximation for inference is valid regardless of the residual shape. Non-normality of residuals therefore matters least at exactly the sample sizes where RWE analysts most often encounter it. The Q-Q plot remains useful for detecting outliers (isolated extreme points far from the reference line) rather than for deciding whether to abandon OLS.

Multicollinearity — arises when two or more predictors are strongly linearly related, which inflates the standard errors of their coefficients without biasing the estimates. The standard diagnostic is the variance inflation factor (VIF): for predictor j, run an auxiliary regression of x_j on all other predictors and compute R²_j; then VIF_j = 1 / (1 - R²_j). A VIF of 1 means no collinearity; VIF ≥ 10 is the traditional alarm threshold (some analysts use 5). A critical distinction: collinearity between control covariates (e.g., age and comorbidity score) inflates their SEs but does not bias the exposure coefficient if the exposure is not itself collinear with them. Collinearity becomes damaging when the exposure of interest is collinear with a control variable — this is the frailty-adjustment problem in claims data, where disease severity proxies are correlated with treatment assignment by design. The appropriate response to high VIF is not automatic deletion of one covariate (which induces omitted-variable bias) but rather domain-guided consolidation (combine correlated proxies into a single composite) or ridge/penalized regression when the research goal is prediction rather than causal inference.

Influence and outliers — a distinction matters between leverage (unusual predictor values) and influence (observations that materially shift the fitted coefficients). An observation can have high leverage but small residual (it fits the model well despite unusual predictors) or high residual but low leverage (a poor fit in a typical part of predictor space). Cook's distance combines both: Cook's D_i ≈ 1 is the rule-of-thumb threshold above which a single observation materially moves the fitted coefficients. In claims data, extremely high-cost observations (transplants, catastrophic illness events) often appear as influential; these are almost always real events, not data errors. The appropriate route is to the cost-outlier handling workflow (winsorization, two-part models) rather than silent deletion.

GLM-specific diagnostics — for logistic regression: Pearson and deviance residuals replace OLS residuals. The deviance residual for observation i is signed(y_i - p̂_i) × √(deviance contribution_i) and should be inspected for systematic patterns versus fitted probabilities. Overdispersion (variance exceeding the distributional assumption) is checked by the ratio of the Pearson chi-squared to its degrees of freedom; for binary logistic regression this ratio should be near 1 (values > 2 suggest clustering or unmodeled heterogeneity; consider GEE or mixed-effects logistic). Calibration — whether predicted probabilities match observed event rates — is assessed with the Hosmer-Lemeshow (HL) test (divide observations into deciles of predicted probability and compare observed vs expected event counts via chi-squared) and calibration plots. The HL test has well-known limitations at large n: it detects trivially small miscalibration with overwhelming power, making formal significance misleading. At claims scale (n ≥ 100,000), the calibration plot with LOESS smoother and the integrated calibration index (ICI) replace the HL p-value as the primary calibration evidence.

The RWE-scale theme: plots over tests

At n = 500,000 patient-years, every formal diagnostic test rejects. Shapiro-Wilk detects infinitesimal non-normality, Breusch-Pagan detects trivial heteroscedasticity, and the Hosmer-Lemeshow test detects calibration error of 0.001 probability units. These test rejections are statistically real but practically irrelevant. The modern approach at RWE scale is to replace test-based decisions with visual inspection of randomly sampled subsets (sample 10,000 rows for diagnostic plots on a dataset of 500,000) combined with effect-size-of-violation reasoning: how large is the detected violation, and does it meaningfully distort the coefficient or its confidence interval? Robust standard errors and GLM re-specification are preemptive defenses against the most common violations; diagnostics then confirm the defense was adequate.

Diagnostics are not specification search

The workflow corruption known as HARKing (Hypothesizing After Results are Known) has a regression analogue: running diagnostics, finding a violation, modifying the model to improve the diagnostic, and reporting the final model as if it had been prespecified. This inflates type-I error because the analyst has effectively performed multiple hypothesis tests on the same data while counting only one. The legitimate role of diagnostics is: (1) verify the prespecified model is adequate; (2) document residual violations transparently; (3) if a major violation is found, report a prespecified alternative specification (e.g., switch from OLS to gamma GLM for costs if the residual plot shows severe heteroscedasticity, as specified in the SAP). Using diagnostics to discover and then incorporate a new predictor, or to justify dropping a confounding variable because its presence creates collinearity, is specification search and belongs in an exploratory section clearly labeled as such.

Interpreting the output

Concrete diagnostic outputs from a claims-based OLS model of 30-day readmission cost (log-transformed) with predictors including age, comorbidity burden, prior utilization, and a binary treatment indicator:

VIF = 4 for the comorbidity burden index covariate. Formal interpretation: the variance of the comorbidity coefficient is 4 times what it would be if comorbidity were orthogonal to all other predictors — the standard error of that coefficient is doubled relative to a perfectly orthogonal design. Practical interpretation: a VIF of 4 is moderate and does not, by itself, demand remediation. Report the wider confidence interval honestly, note that age and comorbidity are correlated (as expected in any claims cohort), and resist the urge to drop comorbidity because it shares variance with age — omitting it would bias the treatment coefficient if comorbidity is also associated with treatment assignment.

Cook's distance D = 0.9 for one observation. Formal interpretation: D ≈ 1 is at the conventional threshold, meaning this single patient materially shifts the fitted coefficient vector if removed. Practical interpretation: investigate the record. In a cost model, a Cook's D near 1 almost always corresponds to a catastrophic-cost outlier — a $2 million transplant episode in a dataset where the 99th percentile is $150,000. This is a real healthcare event, not a data entry error. Route to the cost-outlier-handling workflow: re-run with and without winsorization at the 99th percentile; if results change materially, report both and discuss why the outlier exists clinically.

A funnel-shaped residual-versus-fitted plot. Formal interpretation: heteroscedasticity is present — the error variance increases with the fitted mean, violating OLS assumption (3). Conventional standard errors are negatively biased. Practical interpretation: switch to HC3-robust standard errors for the OLS model and compare to a gamma GLM with log link as the primary specification. In a typical claims cost analysis, the gamma GLM will show tighter confidence intervals because it correctly models the mean-variance relationship.

Hosmer-Lemeshow p = 0.02 at n = 200,000 in a logistic readmission model. Formal interpretation: the test rejects perfect calibration across the deciles of predicted probability. Practical interpretation: at n = 200,000 the HL test has >99% power to detect miscalibration of 0.001 probability units — a difference clinically indistinguishable from perfect calibration. Do not report "the model is poorly calibrated." Instead inspect the calibration plot: if the smoothed observed-probability curve closely tracks the 45-degree line across the range of predicted values, the model is well calibrated in practice. Report the ICI (integrated calibration index) as a quantitative summary and the calibration plot as the primary evidence; relegate the HL p-value to a footnote with appropriate context.

Pros, cons, and trade-offs

Pros of regression diagnostics: - Provide transparent documentation of assumption adequacy, supporting reproducible and auditable HEOR deliverables. - Catch model mis-specification before conclusions are drawn — a residual plot that reveals non-linearity is far cheaper to address before the analysis report is drafted. - Guide appropriate model corrections (robust SEs, GLM re-specification, splines) that improve inference validity without abandoning the prespecified approach. - Required by regulators (FDA, EMA) and HTA bodies for submitted regression models in HEOR dossiers; their absence is an audit flag.

Cons and limitations: - At large RWE n, formal diagnostic tests (Shapiro-Wilk, Breusch-Pagan, HL) have excessive power and reject for trivially small violations — they require interpretation in terms of practical effect size, not just p-values. - Diagnostic plots require visual judgment; two analysts may disagree on whether a residual plot shows "slight" or "material" non-linearity. - A clean diagnostic check does not guarantee valid causal inference — it verifies distributional assumptions, not identification assumptions (no unmeasured confounding, no positivity violations). A perfectly homoscedastic model on a confounded cohort still estimates a biased treatment effect. - Diagnostics invite specification search if the analytical culture does not enforce prespecification discipline.

When to use

Run regression diagnostics after every confirmatory regression model in an HEOR analysis. Specific situations: (1) any OLS model of a continuous outcome (LOS, costs, quality-of-life scores) — residual-vs-fitted and Q-Q plots plus VIF for all predictors; (2) any logistic regression model where calibration is a study deliverable (readmission models, risk stratification, propensity score models) — calibration plot and HL test with appropriate caveats at large n; (3) any model with highly correlated covariates (age + comorbidity, index year + calendar-time variables) — compute VIF before reporting coefficients; (4) cost models — Cook's distance and residual plots to identify influential episodes.

When NOT to use — and misuse modes

Dropping outliers because they are inconvenient. A Cook's D > 1 on a high-cost claims episode does not justify deletion. These are real healthcare costs — the influencing observation is the signal, not noise. Route to cost-outlier-handling methods (winsorization, two-part models, sensitivity analyses). Silently deleting high-influence observations without reporting produces results that underestimate true healthcare costs.

Using the Hosmer-Lemeshow test as the sole calibration evidence at large n. At n ≥ 50,000, HL detects trivial miscalibration with near-certain power. A significant p-value does not mean the model is clinically miscalibrated. Always accompany HL with a calibration plot and a quantitative calibration metric (ICI, E:O ratio) that is interpretable on a clinical scale.

VIF-based automatic deletion of confounders. High VIF in a covariate means its coefficient has a wide confidence interval — not that the variable should be dropped. If the variable is a confounder (associated with both exposure and outcome), dropping it induces omitted-variable bias in the treatment coefficient. The bias from omitting a confounder almost always exceeds the efficiency loss from multicollinearity. The only situation where dropping a collinear predictor is warranted is when two variables measure exactly the same construct (e.g., both Charlson and Elixhauser comorbidity indices in the same model) — and even then, domain justification is required.

Treating diagnostics as a license for post-hoc model modification. Finding a QQ-plot deviation, adding a Box-Cox transformation, refitting, finding a residual pattern, adding a squared term, refitting — all without prespecification — is p-hacking dressed as rigor. If the diagnostic reveals a clear structural issue not anticipated in the SAP, document it as a protocol deviation, propose a prespecified sensitivity analysis, and do not revise the primary model silently.

Worked example

Scenario

A health outcomes analyst builds an OLS model to predict inpatient length of stay (LOS, in days) for five patients using two predictors: age (years) and a comorbidity index (0-5 scale). After fitting, she runs three diagnostics: (1) a VIF calculation to check whether age and comorbidity_index are too correlated to disentangle; (2) a residual computation by hand for the highest-LOS patient; (3) an assessment of which patient has the most influence on the fitted line. The dataset is small enough to work through by hand.

Dataset

Five patients with age, comorbidity index (higher = sicker), and observed LOS. Age and comorbidity are moderately correlated (older patients tend to have more comorbidities), which is typical in claims cohorts and motivates the VIF check.

Steps

• VIF for age: run the auxiliary regression of age on comorbidity_index. The R-squared from that regression is R² = 0.75 (age is moderately explained by comorbidity in this cohort). VIF for age = 1/(1 - R²) = 1/(1 - 0.75) = 1/0.25 = 4.0. A VIF of 4.0 means the standard error for the age coefficient is twice as wide as it would be if age were orthogonal to comorbidity_index (because sqrt(4.0) = 2.0). This is moderate — below the common threshold of 10 — and does not require remediation.

• Residual for P5: after fitting OLS on all five patients, the model predicts los_hat = 7 days for patient P5 (age=80, comorbidity=4). The observed LOS is 10 days. Residual = observed - predicted = 10 - 7 = 3. This is the largest residual in the dataset.

• Influence: P5 has both the highest predictor values (leverage) and the largest residual, which combine to give the highest Cook's distance in this small dataset. With n=5 and p=3 (intercept + 2 predictors), the 50th-percentile threshold for Cook's D is approximately F(0.5, p, n-p) = F(0.5, 3, 2). In practice, any Cook's D > 0.5 in this dataset warrants a sensitivity check. Re-running without P5 would shift the age coefficient noticeably, confirming that P5 is influential.

• Interpretation summary: the VIF = 4.0 signals moderate collinearity between age and comorbidity_index; report wider CIs for both predictors and resist deleting either. The residual of 3 days for P5 and its high influence suggest this patient's episode is unusually costly relative to model predictions — in a real claims dataset this would trigger investigation for a catastrophic-cost record or coding anomaly.

Result

VIF for age = 1/(1-0.75) = 1/0.25 = 4.0 (standard error of age coefficient is 2.0 times wider than under orthogonality, since sqrt(4.0) = 2.0). Residual for P5 = 10 - 7 = 3 (model underpredicts by 3 days). P5 is the most influential observation in the dataset by Cook's distance.

Runnable example

import numpy as np
import statsmodels.api as sm
from statsmodels.stats.outliers_influence import (
    variance_inflation_factor, OLSInfluence
)
from scipy import stats
import matplotlib.pyplot as plt

# ── Teaching dataset from worked_example ──
age  = np.array([45, 55, 65, 70, 80], dtype=float)
comorbidity = np.array([1, 2, 3, 3, 4], dtype=float)
los_obs = np.array([3, 4, 7, 9, 10], dtype=float)

# ── 1. Fit OLS ──
X = sm.add_constant(np.column_stack([age, comorbidity]))
model = sm.OLS(los_obs, X).fit()
print(model.summary())

# ── 2. Variance Inflation Factors ──
# VIF for each predictor column (skip the constant at index 0)
for i, col_name in enumerate(["age", "comorbidity_index"], start=1):
    vif = variance_inflation_factor(X, i)
    print(f"VIF({col_name}) = {vif:.2f}")
# Manual check matching worked_example:
# auxiliary R^2 of age ~ comorbidity_index = 0.75 -> VIF = 1/(1-0.75) = 1/0.25 = 4.0
r2_aux = np.corrcoef(age, comorbidity)[0, 1] ** 2
vif_manual = 1 / (1 - r2_aux)
print(f"VIF from auxiliary R^2={r2_aux:.4f}: {vif_manual:.2f}")

# ── 3. Residuals and Cook's distance ──
influence = OLSInfluence(model)
cooks_d, _ = influence.cooks_distance
residuals = model.resid
fitted    = model.fittedvalues
print(f"\nResiduals:     {residuals.round(2)}")
print(f"Cook's D:      {cooks_d.round(3)}")
# Patient P5: residual = 10 - fitted[4]
print(f"P5 residual = {los_obs[4]:.0f} - {fitted[4]:.2f} = {residuals[4]:.2f}")

# ── 4. Diagnostic plots (four-panel layout) ──
fig, axes = plt.subplots(2, 2, figsize=(10, 8))

# Panel 1: Residuals vs Fitted
axes[0, 0].scatter(fitted, residuals)
axes[0, 0].axhline(0, color="red", linestyle="--")
axes[0, 0].set(xlabel="Fitted values", ylabel="Residuals",
               title="Residuals vs Fitted")

# Panel 2: Q-Q plot of standardized residuals
std_resid = residuals / residuals.std()
sm.qqplot(residuals, line="s", ax=axes[0, 1])
axes[0, 1].set_title("Normal Q-Q")

# Panel 3: Scale-Location (sqrt(|std resid|) vs fitted)
sqrt_abs = np.sqrt(np.abs(std_resid))
axes[1, 0].scatter(fitted, sqrt_abs)
axes[1, 0].set(xlabel="Fitted values", ylabel="sqrt(|Std. residuals|)",
               title="Scale-Location")

# Panel 4: Residuals vs Leverage with Cook's D contours
leverage = influence.hat_matrix_diag
axes[1, 1].scatter(leverage, std_resid, c=cooks_d, cmap="Reds")
axes[1, 1].set(xlabel="Leverage", ylabel="Std. residuals",
               title="Residuals vs Leverage")
for label, hii, ri in zip(["P1","P2","P3","P4","P5"], leverage, std_resid):
    axes[1, 1].annotate(label, (hii, ri))

plt.tight_layout()
plt.savefig("regression_diagnostics.png", dpi=150)

# ── 5. Logistic calibration sketch (Hosmer-Lemeshow decile approach) ──
# For a fitted logistic model with predicted probabilities p_hat and outcomes y:
# (Replace p_hat and y with real model output)
np.random.seed(42)
n_test = 200
p_hat = np.random.beta(2, 5, n_test)        # synthetic predicted probs
y_test = (np.random.rand(n_test) < p_hat).astype(int)

# Decile-based calibration
deciles = np.quantile(p_hat, np.linspace(0, 1, 11))
obs_rates, pred_rates, counts = [], [], []
for lo, hi in zip(deciles[:-1], deciles[1:]):
    mask = (p_hat >= lo) & (p_hat < hi)
    if mask.sum() > 0:
        obs_rates.append(y_test[mask].mean())
        pred_rates.append(p_hat[mask].mean())
        counts.append(mask.sum())

# Hosmer-Lemeshow chi-squared statistic
hl_chi2 = sum(
    n * (o - e)**2 / (e * (1 - e))
    for n, o, e in zip(counts, obs_rates, pred_rates)
    if e > 0 and e < 1
)
hl_df = len(obs_rates) - 2
hl_p = 1 - stats.chi2.cdf(hl_chi2, df=hl_df)
print(f"\nHL chi2 = {hl_chi2:.2f}, df = {hl_df}, p = {hl_p:.4f}")
print("At large n: inspect calibration plot; HL p-value alone is misleading.")

# ── Teaching dataset ──
df <- data.frame(
  patient_id = c("P1","P2","P3","P4","P5"),
  age        = c(45, 55, 65, 70, 80),
  comorbidity= c(1, 2, 3, 3, 4),
  los_obs    = c(3, 4, 7, 9, 10)
)

# ── 1. Fit OLS ──
lm_fit <- lm(los_obs ~ age + comorbidity, data = df)
summary(lm_fit)

# ── 2. Four-panel diagnostic plot ──
# Panels: (1) Residuals vs Fitted, (2) Normal Q-Q,
#         (3) Scale-Location, (4) Residuals vs Leverage
par(mfrow = c(2, 2))
plot(lm_fit)    # base R, labels outliers by row number automatically
par(mfrow = c(1, 1))

# ── 3. Variance Inflation Factors ──
# install.packages("car") if not already installed
library(car)
vif(lm_fit)     # reports VIF per predictor; > 10 is the conventional alarm threshold
# Manual check: auxiliary regression of age on comorbidity
aux_r2 <- summary(lm(age ~ comorbidity, data = df))$r.squared
cat(sprintf("Auxiliary R^2 = %.4f  -> VIF = 1/(1-R^2) = 1/%.4f = %.2f\n",
            aux_r2, 1 - aux_r2, 1 / (1 - aux_r2)))

# ── 4. Influence measures (Cook's D, leverage, DFFITS) ──
infl <- influence.measures(lm_fit)
print(infl)     # flags observations exceeding Cook's D / leverage thresholds with *
# Extract Cook's distance directly
cooks_d <- cooks.distance(lm_fit)
cat("\nCook's distances:\n")
print(round(cooks_d, 4))
cat(sprintf("P5 residual = %.0f - %.2f = %.2f\n",
            df$los_obs[5], fitted(lm_fit)[5], residuals(lm_fit)[5]))

# ── 5. Logistic regression diagnostics ──
# Synthetic binary outcome for demonstration
set.seed(42)
n <- 500
logistic_df <- data.frame(
  age    = rnorm(n, 65, 10),
  comor  = sample(0:5, n, replace = TRUE),
  y      = rbinom(n, 1, 0.25)
)
logit_fit <- glm(y ~ age + comor, data = logistic_df, family = binomial)
summary(logit_fit)

# VIF for logistic predictors
vif(logit_fit)

# Deviance residuals vs fitted probabilities
p_hat <- fitted(logit_fit)
dev_resid <- residuals(logit_fit, type = "deviance")
plot(p_hat, dev_resid, xlab = "Fitted probability", ylab = "Deviance residual",
     main = "Deviance residuals vs fitted (logistic)")
abline(h = 0, lty = 2)

# Calibration plot
cal_data <- data.frame(p_hat = p_hat, y = logistic_df$y)
cal_data$decile <- cut(p_hat, breaks = quantile(p_hat, seq(0, 1, 0.1)),
                       include.lowest = TRUE, labels = FALSE)
cal_summary <- aggregate(cbind(p_hat, y) ~ decile, data = cal_data, FUN = mean)
plot(cal_summary$p_hat, cal_summary$y,
     xlab = "Mean predicted probability", ylab = "Observed event rate",
     main = "Calibration plot (decile means)", xlim = c(0,1), ylim = c(0,1))
abline(0, 1, lty = 2, col = "red")   # 45-degree line of perfect calibration

# Hosmer-Lemeshow test (with large-n caveat: inspect calibration plot first)
# install.packages("ResourceSelection") if needed
library(ResourceSelection)
hl_test <- hoslem.test(logistic_df$y, p_hat, g = 10)
print(hl_test)
cat("At large n, HL p-value is unreliable. Prefer the calibration plot and ICI.\n")