# Regression Diagnostics > Comprehensive regression model diagnostics and assumption checking. Use when validating regression models, checking assumptions (linearity, homoskedasticity, normality, independence), detecting outliers/influence, testing multicollinearity, or when user mentions residuals, heteroskedasticity, VIF, Cook's distance, or diagnostic plots. - **Type:** Skill - **Install:** `agentstack add skill-sshtomar-claude-code-skills-social-science-regression-diagnostics` - **Verified:** Yes — security-reviewed for prompt injection and unsafe behavior - **Seller:** [sshtomar](https://agentstack.voostack.com/s/sshtomar) - **Installs:** 0 - **Category:** [Agent Skills](https://agentstack.voostack.com/c/agent-skills) - **Latest version:** 0.1.0 - **License:** MIT - **Upstream author:** [sshtomar](https://github.com/sshtomar) - **Source:** https://github.com/sshtomar/claude-code-skills-social-science/tree/main/skills/regression-diagnostics ## Install ```sh agentstack add skill-sshtomar-claude-code-skills-social-science-regression-diagnostics ``` Requires the [AgentStack CLI](https://agentstack.voostack.com/docs/cli). Works with Claude Code, Cursor, and any MCP-compatible agent. ## About 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 ```python # 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 ```python @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. ## Pricing - **Free** — Free ## Security capabilities Automated source analysis of v0.1.0 — what this tool can access: - **Network access:** no - **Filesystem access:** no - **Shell / process execution:** no - **Environment & secrets:** no - **Dynamic code execution:** no *"Yes" means the capability is present in the source — more access means more to trust, not that it is unsafe.* ## Versions - **0.1.0** — security scan: passed — Imported from the upstream source. ## Links - Listing page: https://agentstack.voostack.com/l/skill-sshtomar-claude-code-skills-social-science-regression-diagnostics - Seller: https://agentstack.voostack.com/s/sshtomar - Browse the marketplace: https://agentstack.voostack.com/browse --- Listed on AgentStack — the marketplace for AI agent skills and MCP servers. Every listing is security-reviewed. Creators keep 70%.