Let's dive into creating diagnostic plots to check the crucial regression assumptions of linearity and homoscedasticity. These are often among the first assumptions to check, as their violation can lead to biased coefficients, incorrect standard errors, and unreliable p-values.

We'll primarily use residual plots for this. Residuals are the differences between the observed values and the values predicted by your regression model ($ei = yi - \hat{y}_i$). If your model is a good fit and assumptions are met, the residuals should be randomly scattered around zero.

We'll use Python with pandas, numpy, statsmodels, and matplotlib/seaborn.

Key Plots for Linearity & Homoscedasticity

• Residuals vs. Fitted Values Plot:


• Purpose: The primary plot for checking both linearity and homoscedasticity.


• What to look for (Good): A random cloud of points scattered evenly around the horizontal line at zero, with no discernible pattern. The spread of the residuals should be roughly constant across the range of fitted values.


• What to look for (Bad):


• Non-linearity: A curved pattern (e.g., U-shape, inverted U-shape) indicates that the relationship between predictors and the response is not linear. Your linear model is missing a non-linear component.


• Heteroscedasticity: A "fanning out" (cone shape) or "fanning in" (inverted cone shape) pattern indicates that the variance of the residuals is not constant. This means the model's errors are larger for some ranges of fitted values than others.


• Scale-Location Plot (or Spread-Location Plot):


• Purpose: Specifically designed to check for homoscedasticity more clearly than the Residuals vs. Fitted plot, by plotting the square root of the absolute standardized residuals against the fitted values.


• What to look for (Good): A horizontal line with randomly scattered points.


• What to look for (Bad): A non-horizontal line (e.g., upward or downward trend) indicates heteroscedasticity.

Python Implementation

Let's create some example data, run a regression, and then generate these plots. We'll show examples of both good and bad scenarios.

import pandas as pd
import numpy as np
import statsmodels.api as sm
import statsmodels.formula.api as smf
import matplotlib.pyplot as plt
import seaborn as sns
# Set a style for better aesthetics
sns.set_style("whitegrid")
plt.rcParams['figure.figsize'] = (12, 8)
plt.rcParams['font.size'] = 12
# --- 1. Generate Synthetic Data for Demonstration ---
np.random.seed(42)
n_samples = 100
# CASE A: Data that generally meets assumptions (linear, homoscedastic)
X_linear = np.random.rand(n_samples) * 10
y_linear = 2 * X_linear + 5 + np.random.normal(0, 3, n_samples)
df_linear = pd.DataFrame({'X': X_linear, 'y': y_linear})
# CASE B: Data with Non-linearity
X_nonlinear = np.random.rand(n_samples) * 10
y_nonlinear = 0.5 * X_nonlinear**2 - 3 * X_nonlinear + 10 + np.random.normal(0, 3, n_samples) # Quadratic relationship
df_nonlinear = pd.DataFrame({'X': X_nonlinear, 'y': y_nonlinear})
# CASE C: Data with Heteroscedasticity
X_hetero = np.random.rand(n_samples) * 10
y_hetero = 2 * X_hetero + 5 + np.random.normal(0, X_hetero * 0.5, n_samples) # Error variance increases with X
df_hetero = pd.DataFrame({'X': X_hetero, 'y': y_hetero})
print("Synthetic Data Generated.\n")
# --- 2. Run OLS Regression for each case ---
# Model A: Good Case
model_linear = smf.ols('y ~ X', data=df_linear).fit()
print("Model A (Linear, Homoscedastic) Summary:\n", model_linear.summary().tables[1])
# Model B: Non-linear Case (we'll fit a linear model to it to show the violation)
model_nonlinear = smf.ols('y ~ X', data=df_nonlinear).fit()
print("\nModel B (Non-linear Data, Linear Fit) Summary:\n", model_nonlinear.summary().tables[1])
# Model C: Heteroscedastic Case
model_hetero = smf.ols('y ~ X', data=df_hetero).fit()
print("\nModel C (Heteroscedastic Data) Summary:\n", model_hetero.summary().tables[1])
# --- 3. Create Diagnostic Plots ---
def plot_regression_diagnostics(model, title_suffix=""):
    """
    Generates Residuals vs Fitted and Scale-Location plots for a given OLS model.
    """
    # Get residuals and fitted values
    residuals = model.resid
    fitted_values = model.fittedvalues
    
    # Standardize residuals for Scale-Location plot
    std_residuals = model.get_influence().resid_studentized_internal
    sqrt_abs_std_residuals = np.sqrt(np.abs(std_residuals))
    fig, axes = plt.subplots(1, 2, figsize=(14, 6))
    fig.suptitle(f'Regression Diagnostic Plots: {title_suffix}', fontsize=16)
    # Plot 1: Residuals vs. Fitted Values
    sns.scatterplot(x=fitted_values, y=residuals, ax=axes[0], alpha=0.7)
    axes[0].axhline(0, color='red', linestyle='--', linewidth=2)
    # Add a LOESS smoother to help visualize trends
    sns.regplot(x=fitted_values, y=residuals, lowess=True, 
                scatter=False, color='green', line_kws={'linestyle':'-', 'lw':2, 'alpha':0.8},
                ax=axes[0])
    axes[0].set_title('1. Residuals vs. Fitted Values')
    axes[0].set_xlabel('Fitted Values')
    axes[0].set_ylabel('Residuals')
    
    # Plot 2: Scale-Location Plot
    sns.scatterplot(x=fitted_values, y=sqrt_abs_std_residuals, ax=axes[1], alpha=0.7)
    # Add a LOESS smoother
    sns.regplot(x=fitted_values, y=sqrt_abs_std_residuals, lowess=True, 
                scatter=False, color='green', line_kws={'linestyle':'-', 'lw':2, 'alpha':0.8},
                ax=axes[1])
    axes[1].set_title('2. Scale-Location Plot (sqrt(|Standardized Residuals|) vs. Fitted)')
    axes[1].set_xlabel('Fitted Values')
    axes[1].set_ylabel('$\sqrt{|Standardized Residuals|}§§§CODEBLOCK0§§§#39;)
    
    plt.tight_layout(rect=[0, 0.03, 1, 0.95])
    plt.show()
# --- Generate Plots for each case ---
print("\n--- Visualizing Assumptions for Model A (Linear, Homoscedastic) ---")
plot_regression_diagnostics(model_linear, "Model A: Good Assumptions")
print("\n--- Visualizing Assumptions for Model B (Non-linear Data, Linear Fit) ---")
plot_regression_diagnostics(model_nonlinear, "Model B: Non-linearity Detected")
print("\n--- Visualizing Assumptions for Model C (Heteroscedastic Data) ---")
plot_regression_diagnostics(model_hetero, "Model C: Heteroscedasticity Detected")

Interpretation of the Plots

Let's look at what each set of plots tells us:

1. Model A: Good Assumptions (Linear, Homoscedastic)

• Residuals vs. Fitted Values:


• You'll see a random cloud of points centered around the horizontal red line at zero.


• The green LOESS line (Local Regression) will be very close to the red zero line, showing no systematic trend.


• The vertical spread of the points will be roughly consistent across the entire range of fitted values.


• Conclusion: Linearity and homoscedasticity appear to hold.


• Scale-Location Plot:


• The points will be randomly scattered around a horizontal green LOESS line.


• There will be no discernible pattern or trend in the spread.


• Conclusion: Homoscedasticity appears to hold.

2. Model B: Non-linearity Detected

• Residuals vs. Fitted Values:


• You'll observe a clear curved pattern (e.g., a U-shape or inverted U-shape) in the residuals. The green LOESS line will follow this curve, deviating significantly from the red zero line.


• This indicates that your linear model is not adequately capturing the true relationship between X and y.


• Conclusion: Linearity assumption is violated.


• Scale-Location Plot:


• This plot might also show some pattern, but its primary indicator is for heteroscedasticity. In this specific non-linear case, the variance might look somewhat constant, or it might also show some pattern if the non-linearity happens to interact with error variance. The primary evidence for non-linearity comes from the Residuals vs. Fitted plot.


• Conclusion: Heteroscedasticity might be okay, but linearity is definitely an issue.

3. Model C: Heteroscedasticity Detected

• Residuals vs. Fitted Values:


• The points will still be centered around zero (meaning the model isn't systematically biased on average), but the spread of the points will change across the range of fitted values.


• You'll likely see a "fanning out" (cone shape) or "fanning in" pattern. For example, if variance increases with X, the residuals will be tightly packed on the left and spread out on the right.


• The green LOESS line should still be close to zero (as the model is still unbiased in its predictions), but the spread around it is the key.


• Conclusion: Linearity might be fine, but the homoscedasticity assumption is violated.


• Scale-Location Plot:


• This plot will clearly show a non-horizontal trend in the green LOESS line (e.g., an upward slope if variance increases with fitted values).


• The spread of points around this non-horizontal line indicates that the variability of the residuals is not constant.


• Conclusion: Clear violation of homoscedasticity.

What's Next? (Hint for your next request)

These plots are powerful visual tools. After identifying potential violations, you'll often want to back them up with statistical tests. For instance:

• For Linearity: The visual check is often the most direct. If you suspect a specific non-linear form (e.g., quadratic), you can add polynomial terms to your model and check their significance.


• For Homoscedasticity:


• Breusch-Pagan Test


• White Test


• Goldfeld-Quandt Test