Regression Diagnostics

Regression diagnostics are NOT optional add-ons—they are fundamental to valid inference. A regression model is a hypothesis about the data-generating process, and diagnostics test whether that hypothesis holds. Without diagnostics, you're flying blind: your p-values may be wrong, your confidence intervals misleading, and your conclusions invalid. Every regression analysis MUST include diagnostics.

The consequences of ignoring diagnostics are severe: heteroskedasticity can inflate Type I error rates from 5% to 40%, influential outliers can reverse coefficient signs, and multicollinearity can make estimates unstable. This skill enforces rigorous diagnostic practices that protect against these failures.

"All models are wrong, but some are useful" - George Box. Diagnostics tell us HOW wrong our model is and whether it's still useful despite its wrongness.

FUNDAMENTAL TRUTH: Statistical software will happily compute invalid estimates. It's YOUR responsibility to verify assumptions hold. No diagnostic → no inference.

Complete Diagnostic Suite MUST generate all four core diagnostic plots: Residuals vs Fitted, Q-Q plot, Scale-Location, Residuals vs Leverage Each plot tests different assumptions. Missing any plot means missing potential violations (Belsley et al. 1980) Undetected violations lead to invalid inference, wrong p-values, misleading conclusions

Heteroskedasticity Testing MUST run both Breusch-Pagan and White tests, use robust SE if p MacKinnon & White (1985) show HC3 robust SE correct size distortions better than classical SE under heteroskedasticity Type I error rates can reach 40% instead of nominal 5% with uncorrected heteroskedasticity

Multicollinearity Assessment MUST compute VIF for all predictors, flag VIF > 10, recommend action if VIF > 5 Kutner et al. (2004) demonstrate VIF > 10 indicates serious multicollinearity requiring intervention Unstable estimates, inflated standard errors, sign reversals with minor data changes

Influence Diagnostics MUST compute Cook's distance, identify points > 4/n threshold, investigate high-leverage observations Cook & Weisberg (1982) show single influential points can dominate entire regression Results driven by outliers rather than general patterns, non-robust conclusions

Actionable Recommendations For EVERY violation detected, provide specific remediation (e.g., "use robust SE", "log transform", "remove variable X") Diagnostics without remediation are useless—practitioners need actionable guidance Known problems persist, analysis remains flawed despite awareness

When running regression diagnostics:

1. Generate diagnostic plots FIRST (visual inspection reveals patterns)
2. Run formal tests for each assumption
3. Compute influence measures for all observations
4. Check multicollinearity among predictors
5. Document ALL violations found
6. Provide specific remediation for each violation
7. Re-run diagnostics after any remediation
8. Report sensitivity to different specifications

# CRITICAL: Complete regression diagnostics template

@app.cell
def complete_regression_diagnostics(model, df, cluster_var=None):
    """
    Comprehensive diagnostics are mandatory for valid inference.
    Every assumption violation has specific consequences and remediation.
    This function enforces systematic checking and clear reporting.
    """
    import matplotlib.pyplot as plt
    import numpy as np
    import pandas as pd
    from scipy import stats
    from statsmodels.stats.diagnostic import het_breuschpagan, het_white
    from statsmodels.stats.outliers_influence import variance_inflation_factor, OLSInfluence
    import os

    print("=" * 70)
    print("COMPREHENSIVE REGRESSION DIAGNOSTICS")
    print("=" * 70)

    # Extract residuals and fitted values
    fitted = model.fittedvalues
    residuals = model.resid
    standardized_resid = residuals / residuals.std()

    # 1. DIAGNOSTIC PLOTS (mandatory visualization)
    fig, axes = plt.subplots(2, 2, figsize=(14, 12))

    # Residuals vs Fitted (linearity + homoskedasticity)
    axes[0, 0].scatter(fitted, residuals, alpha=0.5, s=30)
    axes[0, 0].axhline(0, color='red', linestyle='--', linewidth=1.5)

    # Add loess smoother to detect patterns
    from scipy.interpolate import interp1d
    sorted_idx = np.argsort(fitted)
    window = max(10, len(fitted) // 20)
    smoothed = pd.Series(residuals[sorted_idx]).rolling(window, center=True).mean()
    axes[0, 0].plot(fitted[sorted_idx], smoothed, 'g-', linewidth=2, label='Loess smoother')

    axes[0, 0].set_xlabel("Fitted Values", fontsize=11)
    axes[0, 0].set_ylabel("Residuals", fontsize=11)
    axes[0, 0].set_title("Residuals vs Fitted\n(Check: Random scatter around zero)", fontsize=12)
    axes[0, 0].grid(True, alpha=0.3)
    axes[0, 0].legend()

    # Q-Q Plot (normality)
    stats.probplot(residuals, dist="norm", plot=axes[0, 1])
    axes[0, 1].set_title("Normal Q-Q Plot\n(Check: Points follow diagonal)", fontsize=12)
    axes[0, 1].grid(True, alpha=0.3)

    # Scale-Location (homoskedasticity)
    sqrt_abs_resid = np.sqrt(np.abs(standardized_resid))
    axes[1, 0].scatter(fitted, sqrt_abs_resid, alpha=0.5, s=30)

    # Add trend line
    z = np.polyfit(fitted, sqrt_abs_resid, 1)
    p = np.poly1d(z)
    axes[1, 0].plot(fitted, p(fitted), "r-", linewidth=2, label=f'Trend: slope={z[0]:.3f}')

    axes[1, 0].set_xlabel("Fitted Values", fontsize=11)
    axes[1, 0].set_ylabel("√|Standardized Residuals|", fontsize=11)
    axes[1, 0].set_title("Scale-Location\n(Check: Horizontal band, no trend)", fontsize=12)
    axes[1, 0].grid(True, alpha=0.3)
    axes[1, 0].legend()

    # Residuals vs Leverage (influence)
    influence = OLSInfluence(model)
    leverage = influence.hat_matrix_diag
    cooks_d = influence.cooks_distance[0]

    axes[1, 1].scatter(leverage, standardized_resid, alpha=0.5, s=30)
    axes[1, 1].axhline(0, color='red', linestyle='--', linewidth=1.5)

    # Mark influential points
    n = len(residuals)
    threshold = 4 / n
    influential_mask = cooks_d > threshold
    if influential_mask.any():
        axes[1, 1].scatter(leverage[influential_mask],
                          standardized_resid[influential_mask],
                          color='red', s=100, alpha=0.7,
                          label=f'Influential (Cook\'s D > {threshold:.3f})')

    # Add contour lines for Cook's distance
    x_leverage = np.linspace(0, max(leverage) * 1.1, 50)
    for cooks_level in [0.5, 1.0]:
        y_resid = np.sqrt(cooks_level * n / (x_leverage * (1 - x_leverage)))
        axes[1, 1].plot(x_leverage, y_resid, '--', color='gray', alpha=0.5)
        axes[1, 1].plot(x_leverage, -y_resid, '--', color='gray', alpha=0.5)

    axes[1, 1].set_xlabel("Leverage", fontsize=11)
    axes[1, 1].set_ylabel("Standardized Residuals", fontsize=11)
    axes[1, 1].set_title("Residuals vs Leverage\n(Check: No influential outliers)", fontsize=12)
    axes[1, 1].grid(True, alpha=0.3)
    axes[1, 1].legend()

    plt.suptitle("Regression Diagnostic Plots - ALL FOUR REQUIRED", fontsize=14, y=1.02)
    plt.tight_layout()

    os.makedirs("./images", exist_ok=True)
    plt.savefig("./images/regression_diagnostics_complete.png", dpi=144, bbox_inches="tight")

    # Save figure reference for inline display
    diagnostic_fig = plt.gcf()

    # 2. HETEROSKEDASTICITY TESTS (mandatory formal testing)
    print("\n" + "=" * 70)
    print("HETEROSKEDASTICITY TESTS")
    print("-" * 70)

    bp_stat, bp_pval, _, _ = het_breuschpagan(residuals, model.model.exog)
    white_stat, white_pval, _, _ = het_white(residuals, model.model.exog)

    print(f"Breusch-Pagan Test:")
    print(f"  Statistic: {bp_stat:.4f}, p-value: {bp_pval:.4f}")
    print(f"  Interpretation: {'VIOLATION - Heteroskedasticity detected' if bp_pval  10 else 'High' if vif > 5 else 'OK'
        })

    vif_df = pd.DataFrame(vif_data).sort_values('VIF', ascending=False)
    print(vif_df.to_string(index=False))

    severe_collinear = vif_df[vif_df['VIF'] > 10]
    high_collinear = vif_df[vif_df['VIF'] > 5]

    if len(severe_collinear) > 0:
        print("\nWARNING: SEVERE MULTICOLLINEARITY DETECTED:")
        for _, row in severe_collinear.iterrows():
            print(f"  {row['Variable']}: VIF = {row['VIF']:.2f}")
        print("\nREQUIRED ACTIONS:")
        print("  1. Remove one of the correlated variables")
        print("  2. Combine into index/principal component")
        print("  3. Use ridge regression for regularization")
    elif len(high_collinear) > 0:
        print("\nWARNING: Moderate multicollinearity detected. Monitor but may be acceptable.")

    # 4. INFLUENCE DIAGNOSTICS (identify problematic observations)
    print("\n" + "=" * 70)
    print("INFLUENTIAL OBSERVATIONS")
    print("-" * 70)

    n_influential = influential_mask.sum()
    high_leverage = leverage > (2 * (X.shape[1]) / n)
    n_high_leverage = high_leverage.sum()

    print(f"Cook's Distance threshold (4/n): {threshold:.4f}")
    print(f"Influential observations: {n_influential} ({100*n_influential/n:.1f}% of data)")
    print(f"High leverage observations: {n_high_leverage} ({100*n_high_leverage/n:.1f}% of data)")

    if n_influential > 0:
        influential_idx = np.where(influential_mask)[0]
        print(f"\nInfluential observation indices: {influential_idx[:10]}")  # Show first 10
        print("\nWARNING: REQUIRED ACTIONS:")
        print("  1. Investigate these observations for data errors")
        print("  2. Check if they represent valid but unusual cases")
        print("  3. Report results with AND without influential points")
        print("  4. Consider robust regression if many influential points")

    # 5. NORMALITY TESTS (for valid inference)
    print("\n" + "=" * 70)
    print("NORMALITY OF RESIDUALS")
    print("-" * 70)

    shapiro_stat, shapiro_pval = stats.shapiro(residuals[:5000])  # Shapiro-Wilk (limit 5000)
    jb_stat, jb_pval = stats.jarque_bera(residuals)

    print(f"Shapiro-Wilk Test: p-value = {shapiro_pval:.4f}")
    print(f"Jarque-Bera Test: p-value = {jb_pval:.4f}")

    if shapiro_pval  30, CLT often ensures valid inference despite non-normality")
        print("  - For small samples, consider bootstrap or transformation")
        print("  - Check for outliers causing non-normality")

    # 6. SUMMARY RECOMMENDATIONS
    print("\n" + "=" * 70)
    print("DIAGNOSTIC SUMMARY & REQUIRED ACTIONS")
    print("=" * 70)

    violations = []
    if bp_pval  0:
        violations.append("Severe multicollinearity → Remove/combine variables")
    if n_influential > 0:
        violations.append("Influential observations → Report sensitivity")
    if (shapiro_pval 

Well-behaved regression with no major violations

```python
@app.cell
def diagnose_clean_model():
    #Even with clean data, diagnostics are mandatory.
    # We verify assumptions hold before trusting inference.

    import pandas as pd
    import numpy as np
    import statsmodels.formula.api as smf

    # Simulated clean data for demonstration
    np.random.seed(42)
    n = 500
    df = pd.DataFrame({
        'x1': np.random.normal(0, 1, n),
        'x2': np.random.normal(0, 1, n),
        'x3': np.random.normal(0, 1, n),
    })

    # Well-specified model with homoskedastic errors
    df['y'] = 2 + 3*df['x1'] - 1.5*df['x2'] + 0.8*df['x3'] + np.random.normal(0, 2, n)

    # Estimate model
    model = smf.ols('y ~ x1 + x2 + x3', data=df).fit()

    print("MODEL SUMMARY")
    print("=" * 70)
    print(model.summary())

    # Run complete diagnostics
    vif_df, influential, cooks = complete_regression_diagnostics(model, df)

    # INTERPRETATION
    print("\n" + "=" * 70)
    print("INTERPRETATION FOR CLEAN MODEL")
    print("=" * 70)
    print("""
    This is an example of a well-behaved regression:

    1. Residuals vs Fitted: Random scatter around zero [OK]
       → Linear specification is appropriate

    2. Q-Q Plot: Points follow diagonal line [OK]
       → Residuals are approximately normal

    3. Scale-Location: Horizontal band, no trend [OK]
       → Variance is constant (homoskedastic)

    4. Residuals vs Leverage: No points outside Cook's distance contours [OK]
       → No influential outliers distorting results

    5. VIF values all  0.05 [OK]
       → Can use standard errors as computed

    CONCLUSION: Standard inference is valid. Report results as-is.
    """)

    return model, vif_df,

Even with "clean" data, running diagnostics is non-negotiable. This example shows what good diagnostics look like—use this as a reference for comparison.

Regression with heteroskedasticity requiring robust standard errors

@app.cell
def diagnose_heteroskedastic_model():
    #Heteroskedasticity is common in cross-sectional data.
    # We detect it, apply robust SE, and show the difference in inference.

    import pandas as pd
    import numpy as np
    import statsmodels.formula.api as smf

    np.random.seed(123)
    n = 500

    # Generate data with heteroskedastic errors
    df = pd.DataFrame({
        'income': np.random.lognormal(10, 1, n),  # Log-normal income
        'education': np.random.normal(12, 3, n),
        'experience': np.random.uniform(0, 40, n),
    })

    # Error variance increases with income (heteroskedasticity)
    error_sd = 0.1 * df['income']  # Variance proportional to income
    df['consumption'] = (
        1000 +
        0.7 * df['income'] +
        50 * df['education'] +
        20 * df['experience'] +
        np.random.normal(0, error_sd)
    )

    # Naive model (ignoring heteroskedasticity)
    model_naive = smf.ols('consumption ~ income + education + experience', data=df).fit()

    print("NAIVE MODEL (Standard Errors)")
    print("=" * 70)
    print(model_naive.summary())

    # Run diagnostics - will detect heteroskedasticity
    vif_df, influential, cooks = complete_regression_diagnostics(model_naive, df)

    # CORRECTED MODEL with robust standard errors
    print("\n" + "=" * 70)
    print("CORRECTED MODEL (HC3 Robust Standard Errors)")
    print("=" * 70)

    # HC3 is preferred for finite samples (MacKinnon & White 1985)
    model_robust = model_naive.get_robustcov_results(cov_type='HC3')
    print(model_robust.summary())

    # Compare standard errors
    comparison = pd.DataFrame({
        'Variable': model_naive.params.index,
        'Coef': model_naive.params.values,
        'SE (Naive)': model_naive.bse.values,
        'SE (Robust)': model_robust.bse.values,
        'SE Ratio': model_robust.bse.values / model_naive.bse.values,
        'p-val (Naive)': model_naive.pvalues.values,
        'p-val (Robust)': model_robust.pvalues.values,
    })

    print("\n" + "=" * 70)
    print("STANDARD ERROR COMPARISON")
    print("=" * 70)
    print(comparison.to_string(index=False))

    print("\n" + "=" * 70)
    print("KEY INSIGHTS")
    print("=" * 70)
    print("""
    1. Heteroskedasticity detected (p 

ALWAYS compare naive vs robust standard errors when heteroskedasticity is detected. The differences can be substantial and change conclusions about statistical significance.

Regression dominated by influential outliers requiring sensitivity analysis

```python
@app.cell
def diagnose_influential_outliers():
    #Influential observations can completely dominate results.
    # We identify them, assess impact, and report sensitivity.

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

    np.random.seed(789)
    n = 200

    # Generate mostly normal data
    df = pd.DataFrame({
        'x': np.random.normal(50, 10, n),
        'z': np.random.normal(0, 1, n),
    })
    df['y'] = 10 + 2*df['x'] + 5*df['z'] + np.random.normal(0, 10, n)

    # Add influential outliers
    outliers = pd.DataFrame({
        'x': [100, 105, 110],  # High leverage points
        'z': [0, 0, 0],
        'y': [50, 45, 40],  # Don't follow pattern
    })
    df = pd.concat([df, outliers],

…

## Source & license

This open-source skill is cataloged on AgentStack and links to its original source — we do not rehost the code.

- **Author:** [sshtomar](https://github.com/sshtomar)
- **Source:** [sshtomar/claude-code-skills-social-science](https://github.com/sshtomar/claude-code-skills-social-science)
- **License:** MIT

Install and usage instructions live in the source repository linked above.