Multiple Regression Analysis

This tutorial provides a comprehensive introduction to Multiple Linear Regression (MLR) using Python, a versatile programming language with powerful libraries for statistical analysis. MLR is a widely used statistical technique that models the relationship between a dependent variable and multiple independent variables. This tutorial guides you through the fundamental concepts of MLR, its applications, and practical implementation in Python using libraries like NumPy, pandas, and statsmodels.

Introduction

Multiple Linear Regression is a statistical method used to examine the relationship between a dependent variable and two or more independent variables. It extends simple linear regression, which studies the relationship between a dependent variable and a single independent variable. In MLR, the goal is to understand how multiple independent variables collectively influence the variation in the dependent variable.

Suppose you have a dataset with \(n\) observations and \(k\) independent variables. The matrix notation for the multiple linear regression model is as follows:

1. Matrix of Independent Variables (Design Matrix):
   • Let \(X\) be the design matrix, with dimensions \(n \times (k+1)\).
   • The first column of \(X\) is typically a column of ones, representing the intercept term.
   • The remaining columns correspond to the values of the independent variables for each observation.

\[ X = \begin{bmatrix} 1 & x_{11} & x_{12} & \ldots & x_{1k} \\ 1 & x_{21} & x_{22} & \ldots & x_{2k} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{n1} & x_{n2} & \ldots & x_{nk} \end{bmatrix} \]

2. Vector of Coefficients:
   • Let \(\beta\) be the vector of coefficients, with dimensions \((k+1) \times 1\).
   • The elements of \(\beta\) represent the intercept and the coefficients for the independent variables.

\[ \beta = \begin{bmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \\ \vdots \\ \beta_k \end{bmatrix} \]

3. Vector of Dependent Variables:
   • Let \(Y\) be the vector of dependent variables, with dimensions \(n \times 1\).
   • Each element of \(Y\) corresponds to the dependent variable for each observation.

\[ Y = \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix} \]

4. Error Vector:
   • Let \(\varepsilon\) be the vector of errors (residuals), with dimensions \(n \times 1\).
   • Each element of \(\varepsilon\) represents the difference between the observed and predicted values.

\[ \varepsilon = \begin{bmatrix} \varepsilon_1 \\ \varepsilon_2 \\ \vdots \\ \varepsilon_n \end{bmatrix} \]

The multiple linear regression model can be expressed in matrix notation as:

\[ Y = X\beta + \varepsilon \]

This equation captures the relationship between the dependent variable \(Y\), the design matrix \(X\), the coefficient vector \(\beta\), and the error vector \(\varepsilon\). The goal of regression analysis is to estimate the coefficients in \(\beta\) that minimize the sum of squared errors (least squares method). This estimation can be expressed in matrix form as:

\[ \hat{\beta} = (X^TX)^{-1}X^TY \]

where \(\hat{\beta}\) represents the estimated coefficients. Matrix notation provides a concise and efficient way to represent and solve regression problems, especially with multiple variables.

Multiple Linear Regression Analysis from Scratch

Below is an example of performing multiple regression from scratch in Python with four predictors, including summary statistics and R² values:

import numpy as np
import pandas as pd
from scipy import stats

# Set seed for reproducibility
np.random.seed(123)

# Generate sample data
n = 100  # Number of observations
X = np.random.normal(size=(n, 4))  # Four predictors
X = pd.DataFrame(X, columns=['X1', 'X2', 'X3', 'X4'])
true_beta = np.array([2, 3, -1, 0.5, 2])  # True coefficients (intercept first)
X_with_intercept = np.c_[np.ones(n), X]  # Design matrix with intercept
error = np.random.normal(0, 2, n)  # Error term
y = X_with_intercept @ true_beta + error  # Dependent variable

# Calculate coefficients using normal equation
XTX = X_with_intercept.T @ X_with_intercept
beta_hat = np.linalg.inv(XTX) @ X_with_intercept.T @ y

# Calculate residuals and R-squared
y_hat = X_with_intercept @ beta_hat
residuals = y - y_hat
SS_residual = np.sum(residuals**2)
SS_total = np.sum((y - np.mean(y))**2)
R_squared = 1 - (SS_residual / SS_total)

# Calculate standard errors and p-values
df = n - len(beta_hat)  # Degrees of freedom
sigma_sq_hat = SS_residual / df  # Residual variance estimate
var_cov = sigma_sq_hat * np.linalg.inv(XTX)  # Variance-covariance matrix
se = np.sqrt(np.diag(var_cov))  # Standard errors
t_values = beta_hat / se  # t-statistics
p_values = 2 * stats.t.sf(np.abs(t_values), df=df)  # p-values

# Create summary table
summary_df = pd.DataFrame({
    'Estimate': np.round(beta_hat, 4),
    'Std.Error': np.round(se, 4),
    't.value': np.round(t_values, 4),
    'p.value': np.where(p_values < 0.001, '<0.001', np.round(p_values, 4))
}, index=['(Intercept)', 'X1', 'X2', 'X3', 'X4'])

# Display results
print("Multiple Regression Results from Scratch")
print("=" * 40)
print(summary_df)
print(f"\nResidual standard error: {np.round(np.sqrt(sigma_sq_hat), 3)} on {df} degrees of freedom")
print(f"Multiple R-squared: {np.round(R_squared, 4)}")

Multiple Linear Regression Analysis in Python

In this section, we explore how to perform MLR in Python using the statsmodels library, interpret the results, and visualize the findings. We will also cover model diagnostics, performance evaluation, and relative importance of predictors.

Check and Install Required Packages

import importlib.util
import sys

# List of required packages
packages = ['numpy', 'pandas', 'statsmodels', 'scikit-learn', 'matplotlib', 'seaborn']

# Check and install missing packages
# Check and install missing packages
for package in packages:
    if not importlib.util.find_spec(package):
        try:
            import pip
            pip.main(['install', package])
        except ImportError:
            print(f"Failed to install {package}. Pip is not available.")

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
from statsmodels.stats.outliers_influence import variance_inflation_factor
from sklearn.metrics import r2_score

Data

# Load data
mf = pd.read_csv("https://github.com/zia207/Python_for_Beginners/raw/refs/heads/main/Data/gp_soil_data.csv")
print(mf.head())

Fit an MLR Model with Continuous Predictors

We aim to create an MLR model incorporating predictors such as DEM (Digital Elevation Model), Slope, MAT (Mean Annual Temperature), MAP (Mean Annual Precipitation), and NDVI (Normalized Difference Vegetation Index) to predict soil organic carbon (SOC).

# Create a DataFrame with selected variables
df = mf[['SOC', 'DEM', 'Slope', 'MAT', 'MAP', 'NDVI']]

# Fit an MLR model
X = df[['DEM', 'Slope', 'MAT', 'MAP', 'NDVI']]
X = sm.add_constant(X)  # Add intercept
y = df['SOC']
mlr_01 = sm.OLS(y, X).fit()

# Display summary
print(mlr_01.summary())

Create a Regression Summary Table

The statsmodels library provides a detailed summary of the regression model. For a custom summary table, you can extract relevant statistics:

# Create a summary table
summary_df = pd.DataFrame({
    'Estimate': mlr_01.params,
    'Std.Error': mlr_01.bse,
    't.value': mlr_01.tvalues,
    'p.value': mlr_01.pvalues
}).round(4)
print(summary_df)

Model Performance

To evaluate model performance, we can calculate metrics like R-squared and adjusted R-squared:

# Model performance
print(f"R-squared: {mlr_01.rsquared:.4f}")
print(f"Adjusted R-squared: {mlr_01.rsquared_adj:.4f}")

Comprehensive Visualization of Model Checks

To visualize model diagnostics, we can create plots for residuals, Q-Q plots, and fitted vs. residual plots:

# Residuals vs Fitted
plt.figure(figsize=(6.5, 4))
sns.scatterplot(x=mlr_01.fittedvalues, y=mlr_01.resid)
plt.axhline(0, color='red', linestyle='--')
plt.xlabel('Fitted Values')
plt.ylabel('Residuals')
plt.title('Residuals vs Fitted')
plt.show()

# Q-Q Plot
plt.figure(figsize=(6.5, 4))
stats.probplot(mlr_01.resid, dist="norm", plot=plt)
plt.title('Q-Q Plot')
plt.show()

Fit an MLR Model with Continuous and Categorical Predictors

Now, we update the MLR model by adding categorical variables: NLCD (National Land Cover Database) and FRG (Forest Regimes).

# Create a DataFrame with selected variables
df_02 = mf[['SOC', 'DEM', 'Slope', 'MAT', 'MAP', 'NDVI', 'NLCD', 'FRG']]

# Drop rows with missing values (or impute as needed)
df_02 = df_02.dropna()

# Ensure categorical columns are strings
df_02['NLCD'] = df_02['NLCD'].astype(str)
df_02['FRG'] = df_02['FRG'].astype(str)

# Convert categorical variables to dummy variables
df_02 = pd.get_dummies(df_02, columns=['NLCD', 'FRG'], drop_first=True)

# Ensure all columns in X are numeric
X = df_02.drop('SOC', axis=1)
X = X.astype(float)  # Convert to float
X = sm.add_constant(X)  # Add intercept
y = df_02['SOC'].astype(float)  # Ensure y is numeric

# Check data types
print("Data types in X:\n", X.dtypes)

# Fit the MLR model
mlr_02 = sm.OLS(y, X).fit()

# Display summary
print(mlr_02.summary())

Model Performance

# Model performance
print(f"R-squared: {mlr_02.rsquared:.4f}")
print(f"Adjusted R-squared: {mlr_02.rsquared_adj:.4f}")

Model Comparison

To compare the performance of mlr_01 and mlr_02, we can compare their R-squared and adjusted R-squared values:

# Model comparison
comparison_df = pd.DataFrame({
    'Model': ['mlr_01', 'mlr_02'],
    'R-squared': [mlr_01.rsquared, mlr_02.rsquared],
    'Adjusted R-squared': [mlr_01.rsquared_adj, mlr_02.rsquared_adj]
}).round(4)
print(comparison_df)

Visualizing Model Predictions

To visualize the effect of a predictor (e.g., NDVI) on the dependent variable:

# Effect of a continuous variable (NDVI)
plt.figure(figsize=(6.5, 4))
sns.scatterplot(x=df['NDVI'], y=df['SOC'], label='Data')
x_range = np.linspace(df['NDVI'].min(), df['NDVI'].max(), 100)
X_pred = pd.DataFrame({
    'const': 1,
    'DEM': df['DEM'].mean(),
    'Slope': df['Slope'].mean(),
    'MAT': df['MAT'].mean(),
    'MAP': df['MAP'].mean(),
    'NDVI': x_range
})
y_pred = mlr_01.predict(X_pred)
plt.plot(x_range, y_pred, color='red', label='Prediction')
plt.xlabel('NDVI')
plt.ylabel('SOC')
plt.title('Effect of NDVI on SOC')
plt.legend()
plt.show()

Variance Inflation Factors (VIF)

VIF measures multicollinearity among predictors. A VIF greater than 5 or 10 indicates potential multicollinearity issues.

# Calculate VIF
X = df[['DEM', 'Slope', 'MAT', 'MAP', 'NDVI']]
vif_data = pd.DataFrame()
vif_data['Variable'] = X.columns
vif_data['VIF'] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]
print(vif_data)

Relative Importance

To calculate the relative importance of predictors, we can use the shap library or compute standardized coefficients:

# Standardized coefficients
X_std = (X - X.mean()) / X.std()
X_std = sm.add_constant(X_std)
mlr_std = sm.OLS((y - y.mean()) / y.std(), X_std).fit()
print("Standardized Coefficients:")
print(mlr_std.params[1:].round(4))  # Exclude intercept

Summary and Conclusion

This tutorial provided a comprehensive overview of Multiple Linear Regression using Python, a powerful tool for modeling relationships between multiple variables. We covered data preparation, model building, and evaluation using statsmodels. By fitting MLR models, we interpreted coefficients to understand predictor impacts on the response variable. Diagnostic plots and statistical tests validated model assumptions.

MLR is crucial in fields like finance, economics, and social sciences, enabling informed predictions and insights into complex variable interactions. To further enhance your models, explore advanced topics like feature engineering, interaction terms, and regularization. This tutorial serves as a foundation for applying MLR to your datasets, uncovering patterns and relationships that drive meaningful insights.

Resources

1. Statsmodels Documentation
2. Scikit-learn Documentation
3. Pandas Documentation
4. Seaborn Documentation