Regression analysis is a statistical method that examines the relationship between dependent and independent variables, allowing us to model and predict outcomes based on input data. SciPy’s stats module provides several useful options for performing different types of regression analysis. In this article, we will explore few of these methods.

Getting Started with SciPy Regression

First, let’s import the necessary modules:

from scipy import stats
import numpy as np

Basic Linear Regression Analysis

The linregress() function provides a comprehensive way to perform linear regression in SciPy. Here’s how to implement it:

# Create sample data
np.random.seed(42)
x = np.linspace(0, 10, 30)
y = 2 * x + 1 + np.random.normal(0, 1.5, 30)

# Perform linear regression
result = stats.linregress(x, y)

# Print results
print(f"Slope: {result.slope:.4f}")
print(f"Intercept: {result.intercept:.4f}")
print(f"R-squared: {result.rvalue**2:.4f}")
print(f"P-value: {result.pvalue:.4f}")
print(f"Standard Error: {result.stderr:.4f}")

Slope: 1.8471
Intercept: 1.4821
R-squared: 0.9514
P-value: 0.0000
Standard Error: 0.0789

Let’s break down these results:

1. Slope (1.8471): This indicates that for each one-unit increase in x, y increases by approximately 1.85 units. Since we generated our data with a true slope of 2, our model has captured this relationship quite well.
2. Intercept (1.4821): This is the predicted value of y when x equals 0. Our data was generated with a true intercept of 1, so our model’s estimate is reasonably close.
3. R-squared (0.9514): This value indicates that approximately 95.14% of the variance in y is explained by our model. This is a very strong fit, suggesting that our linear model captures the relationship between x and y very well.
4. P-value (0.0000): The extremely low p-value (< 0.05) indicates that our regression is statistically significant. We can reject the null hypothesis that there is no relationship between x and y.
5. Standard Error (0.0789): This represents the average distance that the observed values fall from the regression line. The small standard error suggests our slope estimate is precise.

To calculate confidence intervals for our regression line, we can use the standard error:

# Calculate 95% confidence intervals for the slope
confidence_level = 0.95
degrees_of_freedom = len(x) - 2  # n-2 degrees of freedom
t_value = stats.t.ppf((1 + confidence_level) / 2, degrees_of_freedom)

slope_ci_lower = result.slope - t_value * result.stderr
slope_ci_upper = result.slope + t_value * result.stderr

print(f"\n95% Confidence Interval for Slope: [{slope_ci_lower:.4f}, {slope_ci_upper:.4f}]")

95% Confidence Interval for Slope: [1.6855, 2.0088]

The confidence interval indicates we can be 95% confident that the true slope falls between 1.69 and 2.01, which nicely contains our true slope value of 2 that we used to generate the data.

Robust Regression with Siegel Slopes

Let’s explore how different regression methods handle outliers. We’ll use the exact same dataset from our previous example, but introduce a single outlier:

# Use the same data but add an outlier
np.random.seed(42)
x = np.linspace(0, 10, 30)
y = 2 * x + 1 + np.random.normal(0, 1.5, 30)
y[25] = 40  # Replace one value with an outlier

# Calculate Siegel slopes
slope, intercept = stats.siegelslopes(y, x)

# Compare with regular linear regression
regular_result = stats.linregress(x, y)

print("Siegel Slopes Results:")
print(f"Robust Slope: {slope:.4f}")
print(f"Robust Intercept: {intercept:.4f}")

print("\nRegular Linear Regression Results:")
print(f"Regular Slope: {regular_result.slope:.4f}")
print(f"Regular Intercept: {regular_result.intercept:.4f}")

Siegel Slopes Results:
Robust Slope: 1.8632
Robust Intercept: 1.3680

Regular Linear Regression Results:
Regular Slope: 2.1397
Regular Intercept: 0.7391

The Siegel slopes method provides estimates close to our true parameters (slope = 2, intercept = 1) despite the outlier. Notice how the regular linear regression slope is pulled higher (2.1397) by the outlier, while the Siegel slopes method maintains an estimate closer to our initial results without the outlier (1.8632). This is because Siegel slopes uses median regression principles, making it less sensitive to extreme values.

Theil-Sen Regression: An Alternative Robust Approach

The Theil-Sen estimator offers another robust approach to regression that differs from both regular linear regression and Siegel slopes:

# Perform Theil-Sen regression using the same data
slope, intercept, lo_slope, up_slope = stats.theilslopes(y, x, 0.95)

print("Theil-Sen Results:")
print(f"Slope: {slope:.4f}")
print(f"Intercept: {intercept:.4f}")
print(f"95% Confidence Interval for Slope: [{lo_slope:.4f}, {up_slope:.4f}]")

Theil-Sen Results:
Slope: 1.8573
Intercept: 0.7838
95% Confidence Interval for Slope: [1.6964, 2.0267]

While both Siegel slopes and Theil-Sen are robust methods, they use different approaches to achieve this robustness. Theil-Sen calculates slopes between all possible pairs of points and uses their median, while Siegel slopes uses repeated medians. In this case, both methods give us similar slopes (1.8632 and 1.8573), suggesting strong agreement about the true relationship in our data despite the outlier.

Choosing the Right Method: Comparison and Conclusions

Our analysis has demonstrated three different regression approaches in SciPy, each with its own strengths and use cases:

1. Linear Regression (linregress):
   • Best for clean data with no outliers
   • Provides comprehensive statistics (R-squared, p-values, standard errors)
   • Most efficient when assumptions about normal distribution are met
   • As we saw in our initial analysis, produced excellent results with clean data (slope ≈ 1.85)
2. Siegel Slopes (siegelslopes):
   • Resistant to outliers through median regression
   • Particularly useful when you suspect extreme values in your data
   • Computationally simpler than Theil-Sen
   • In our outlier example, maintained accuracy (slope ≈ 1.86) despite the extreme value
3. Theil-Sen (theilslopes):
   • Robust against outliers using pairwise slopes
   • Provides confidence intervals for slope estimates
   • In our analysis, produced similar robust results (slope ≈ 1.86) with added confidence bounds

When working with regression analysis, SciPy provides you with few options that lets you quickly adapt to different data scenarios. Before diving into analysis, it’s important to examine your data for outliers and consider which statistical measures you need—whether that’s comprehensive statistics from linear regression, robust estimates from Siegel slopes, or confidence intervals from Theil-Sen regression.

The strength of SciPy’s regression tools lies in their accessibility and ease of use. With just a few lines of code, you can switch between different regression methods to validate your findings or handle challenging data. This flexibility makes SciPy a good option for data scientists and analysts who need to quickly explore different regression approaches before settling on the most appropriate method for their specific dataset.