Interpreting Regression Results in Python

Introduction to Regression Analysis

Regression analysis is a powerful statistical method used to examine the relationship between variables. It helps us understand how changes in one or more independent variables affect a dependent variable. In this presentation, we'll explore how to interpret regression results using Python, focusing on practical examples and actionable insights.

import pandas as pd
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split

# Generate sample data
np.random.seed(42)
X = np.random.rand(100, 1) * 10
y = 2 * X + 1 + np.random.randn(100, 1)

# Plot the data
plt.scatter(X, y)
plt.xlabel('Independent Variable')
plt.ylabel('Dependent Variable')
plt.title('Sample Data for Regression Analysis')
plt.show()

Simple Linear Regression

Simple linear regression models the relationship between two variables using a straight line. It's used when we have one independent variable and one dependent variable. Let's fit a simple linear regression model to our sample data and interpret the results.

model = LinearRegression()
model.fit(X, y)

# Print model coefficients
print(f"Intercept: {model.intercept_[0]:.2f}")
print(f"Slope: {model.coef_[0][0]:.2f}")

# Plot the regression line
plt.scatter(X, y)
plt.plot(X, model.predict(X), color='red')
plt.xlabel('Independent Variable')
plt.ylabel('Dependent Variable')
plt.title('Simple Linear Regression')
plt.show()

Interpreting Coefficients

The coefficients of a linear regression model provide valuable insights into the relationship between variables. The intercept represents the expected value of the dependent variable when all independent variables are zero. The slope (coefficient) indicates how much the dependent variable changes for a one-unit increase in the independent variable.

y_pred = model.predict(X)
residuals = y - y_pred

# Plot residuals
plt.scatter(X, residuals)
plt.axhline(y=0, color='red', linestyle='--')
plt.xlabel('Independent Variable')
plt.ylabel('Residuals')
plt.title('Residual Plot')
plt.show()

# Print mean of residuals
print(f"Mean of residuals: {np.mean(residuals):.4f}")

R-squared and Model Fit

R-squared (R²) is a statistical measure that represents the proportion of variance in the dependent variable explained by the independent variables. It ranges from 0 to 1, with higher values indicating a better fit. Let's calculate and interpret the R-squared value for our model.

r_squared = model.score(X, y)
print(f"R-squared: {r_squared:.4f}")

# Visualize model fit
plt.scatter(X, y)
plt.plot(X, model.predict(X), color='red')
plt.xlabel('Independent Variable')
plt.ylabel('Dependent Variable')
plt.title(f'Model Fit (R-squared: {r_squared:.4f})')
plt.show()

Multiple Linear Regression

Multiple linear regression extends simple linear regression to include multiple independent variables. This allows us to model more complex relationships and account for multiple factors influencing the dependent variable.

X_multi = np.column_stack((X, np.random.rand(100, 1) * 5))
y_multi = 2 * X_multi[:, 0] + 3 * X_multi[:, 1] + 1 + np.random.randn(100, 1)

# Fit multiple linear regression model
model_multi = LinearRegression()
model_multi.fit(X_multi, y_multi)

# Print coefficients
print(f"Intercept: {model_multi.intercept_[0]:.2f}")
print(f"Coefficient 1: {model_multi.coef_[0][0]:.2f}")
print(f"Coefficient 2: {model_multi.coef_[0][1]:.2f}")

Interpreting Multiple Regression Coefficients

In multiple regression, each coefficient represents the change in the dependent variable for a one-unit increase in the corresponding independent variable, holding all other variables constant. This concept is known as "ceteris paribus" in economics.

def plot_partial_effect(X, y, model, feature_index, feature_name):
    plt.scatter(X[:, feature_index], y)
    
    X_plot = np.linspace(X[:, feature_index].min(), X[:, feature_index].max(), 100).reshape(-1, 1)
    X_mean = X.mean(axis=0)
    X_pred = np.tile(X_mean, (100, 1))
    X_pred[:, feature_index] = X_plot.ravel()
    
    y_pred = model.predict(X_pred)
    plt.plot(X_plot, y_pred, color='red')
    
    plt.xlabel(feature_name)
    plt.ylabel('Dependent Variable')
    plt.title(f'Partial Effect of {feature_name}')
    plt.show()

# Plot partial effects for both features
plot_partial_effect(X_multi, y_multi, model_multi, 0, 'Feature 1')
plot_partial_effect(X_multi, y_multi, model_multi, 1, 'Feature 2')

Multicollinearity

Multicollinearity occurs when independent variables in a regression model are highly correlated with each other. This can lead to unstable and unreliable coefficient estimates. Let's explore how to detect and address multicollinearity.

X_collinear = np.column_stack((X, X + np.random.randn(100, 1) * 0.1))

# Calculate correlation matrix
corr_matrix = np.corrcoef(X_collinear.T)

# Visualize correlation matrix
plt.imshow(corr_matrix, cmap='coolwarm')
plt.colorbar()
plt.xticks([0, 1], ['X1', 'X2'])
plt.yticks([0, 1], ['X1', 'X2'])
plt.title('Correlation Matrix')
plt.show()

# Calculate Variance Inflation Factor (VIF)
from statsmodels.stats.outliers_influence import variance_inflation_factor

vif_data = pd.DataFrame()
vif_data["Variable"] = ["X1", "X2"]
vif_data["VIF"] = [variance_inflation_factor(X_collinear, i) for i in range(X_collinear.shape[1])]
print(vif_data)

Handling Non-linear Relationships

Not all relationships between variables are linear. When faced with non-linear patterns, we can use polynomial regression or other non-linear transformations to capture these relationships more accurately.

X_nonlinear = np.linspace(0, 10, 100).reshape(-1, 1)
y_nonlinear = 2 * X_nonlinear**2 + X_nonlinear + 5 + np.random.randn(100, 1) * 10

# Fit linear and polynomial models
from sklearn.preprocessing import PolynomialFeatures

linear_model = LinearRegression().fit(X_nonlinear, y_nonlinear)
poly_features = PolynomialFeatures(degree=2)
X_poly = poly_features.fit_transform(X_nonlinear)
poly_model = LinearRegression().fit(X_poly, y_nonlinear)

# Plot results
plt.scatter(X_nonlinear, y_nonlinear)
plt.plot(X_nonlinear, linear_model.predict(X_nonlinear), color='red', label='Linear')
plt.plot(X_nonlinear, poly_model.predict(X_poly), color='green', label='Polynomial')
plt.xlabel('X')
plt.ylabel('y')
plt.title('Linear vs Polynomial Regression')
plt.legend()
plt.show()

Residual Analysis

Residual analysis is crucial for validating the assumptions of linear regression. By examining residuals, we can check for homoscedasticity, normality, and the presence of outliers or influential points.

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
model = LinearRegression().fit(X_train, y_train)
y_pred = model.predict(X_test)
residuals = y_test - y_pred

# Plot residuals vs predicted values
plt.scatter(y_pred, residuals)
plt.axhline(y=0, color='red', linestyle='--')
plt.xlabel('Predicted Values')
plt.ylabel('Residuals')
plt.title('Residuals vs Predicted Values')
plt.show()

# Q-Q plot for normality check
from scipy import stats
fig, ax = plt.subplots()
stats.probplot(residuals.ravel(), plot=ax)
ax.set_title("Q-Q plot")
plt.show()

Feature Selection

Feature selection is the process of choosing the most relevant independent variables for your model. This can improve model performance, reduce overfitting, and enhance interpretability.

# Generate sample data with multiple features
X_multi = np.random.rand(100, 5)
y_multi = 2*X_multi[:, 0] + 3*X_multi[:, 2] + np.random.randn(100)

# Perform Recursive Feature Elimination
model = LinearRegression()
rfe = RFE(estimator=model, n_features_to_select=2)
rfe = rfe.fit(X_multi, y_multi)

# Print selected features
print("Selected features:")
for i, selected in enumerate(rfe.support_):
    if selected:
        print(f"Feature {i+1}")

# Plot feature importance
plt.bar(range(1, 6), rfe.ranking_)
plt.xlabel('Feature')
plt.ylabel('Ranking')
plt.title('Feature Importance')
plt.show()

Cross-validation

Cross-validation helps assess how well our regression model generalizes to unseen data. It involves splitting the data into multiple subsets, training the model on some subsets, and evaluating it on others.

# Perform k-fold cross-validation
cv_scores = cross_val_score(LinearRegression(), X, y, cv=5)

print(f"Cross-validation scores: {cv_scores}")
print(f"Mean CV score: {cv_scores.mean():.4f}")
print(f"Standard deviation of CV scores: {cv_scores.std():.4f}")

# Visualize cross-validation results
plt.boxplot(cv_scores)
plt.title('Cross-validation Scores')
plt.ylabel('R-squared')
plt.show()

Real-life Example: Predicting House Prices

Let's apply our regression knowledge to a real-world scenario: predicting house prices based on various features such as size, number of bedrooms, and location.

np.random.seed(42)
n_samples = 1000
size = np.random.uniform(1000, 5000, n_samples)
bedrooms = np.random.randint(1, 6, n_samples)
location_score = np.random.uniform(1, 10, n_samples)
age = np.random.uniform(0, 50, n_samples)

# Create price with some noise
price = 100000 + 100 * size + 25000 * bedrooms + 50000 * location_score - 1000 * age
price += np.random.normal(0, 50000, n_samples)

# Combine features
X = np.column_stack((size, bedrooms, location_score, age))
y = price

# Split data and fit model
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
model = LinearRegression().fit(X_train, y_train)

# Print coefficients and R-squared
print("Coefficients:")
print(f"Size: ${model.coef_[0]:.2f} per sq ft")
print(f"Bedrooms: ${model.coef_[1]:.2f} per bedroom")
print(f"Location Score: ${model.coef_[2]:.2f} per point")
print(f"Age: ${model.coef_[3]:.2f} per year")
print(f"\nR-squared: {model.score(X_test, y_test):.4f}")

# Predict price for a sample house
sample_house = np.array([[2500, 3, 7, 15]])
predicted_price = model.predict(sample_house)[0]
print(f"\nPredicted price for sample house: ${predicted_price:.2f}")

Real-life Example: Analyzing Factors Affecting Plant Growth

In this example, we'll use regression analysis to understand how various environmental factors influence plant growth. This application demonstrates the versatility of regression in fields such as agriculture and ecology.

np.random.seed(42)
n_samples = 500
sunlight = np.random.uniform(2, 12, n_samples)  # hours of sunlight per day
water = np.random.uniform(100, 500, n_samples)  # ml of water per day
temperature = np.random.uniform(15, 35, n_samples)  # average temperature in Celsius
soil_quality = np.random.uniform(1, 10, n_samples)  # soil quality score

# Create growth rate with some noise
growth_rate = 0.5 + 0.1 * sunlight + 0.001 * water + 0.05 * temperature + 0.2 * soil_quality
growth_rate += np.random.normal(0, 0.5, n_samples)

# Combine features
X = np.column_stack((sunlight, water, temperature, soil_quality))
y = growth_rate

# Split data and fit model
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
model = LinearRegression().fit(X_train, y_train)

# Print coefficients and R-squared
print("Coefficients:")
print(f"Sunlight: {model.coef_[0]:.4f} cm/day per hour of sunlight")
print(f"Water: {model.coef_[1]:.4f} cm/day per ml of water")
print(f"Temperature: {model.coef_[2]:.4f} cm/day per degree Celsius")
print(f"Soil Quality: {model.coef_[3]:.4f} cm/day per soil quality point")
print(f"\nR-squared: {model.score(X_test, y_test):.4f}")

# Predict growth rate for a sample plant environment
sample_environment = np.array([[8, 300, 25, 7]])
predicted_growth = model.predict(sample_environment)[0]
print(f"\nPredicted growth rate for sample environment: {predicted_growth:.2f} cm/day")

# Visualize the effect of sunlight on growth rate
plt.scatter(X_test[:, 0], y_test, alpha=0.5)
plt.plot(X_test[:, 0], model.predict(X_test), color='red')
plt.xlabel('Sunlight (hours/day)')
plt.ylabel('Growth Rate (cm/day)')
plt.title('Effect of Sunlight on Plant Growth Rate')
plt.show()

Additional Resources

For those interested in diving deeper into regression analysis and its applications in Python, here are some valuable resources:

1. "An Introduction to Statistical Learning" by Gareth James, Daniela Witten, Trevor Hastie, and Robert Tibshirani - A comprehensive guide to statistical learning methods, including regression analysis.
2. "Python for Data Analysis" by Wes McKinney - An excellent resource for learning how to use Python for data manipulation and analysis, including regression techniques.
3. Scikit-learn documentation (https://scikit-learn.org/stable/modules/linear_model.html) - The official documentation for scikit-learn's linear models, including various regression techniques.
4. "Applied Predictive Modeling" by Max Kuhn and Kjell Johnson - A practical guide to predictive modeling techniques, with a focus on real-world applications.
5. ArXiv paper: "A Survey of Deep Learning Techniques for Neural Machine Translation" by Shuohang Wang and Jing Jiang (https://arxiv.org/abs/1703.01619) - While not directly about regression, this paper provides insights into advanced modeling techniques that can be applied to various prediction tasks.

These resources will help you expand your knowledge of regression analysis and its implementation in Python, from basic concepts to advanced techniques and real-world applications.