[Free Template] Linear Regression in Python with Fit Diagnostics Matrix

This template has the code, explanation, sample output where you can just change the dataset and variable names only and get a complete regression report

Regression

The linear regression task fits a linear model to predict a single continuous dependent variable from one or more continuous or categorical predictor variables. This task produces statistics and graphs for interpreting the results.

Simple linear regression

In simple linear regression, there is only one independent variable in the model. The line equation will be as follows:

Y = a + b X

Where Y is the dependent variable that we want to predict its value, a is the intercept, b is the coefficient, and X is the independent variable. The error is ignored.

The null hypothesis is that there is no linear relationship between the variables. In other words, the value of b is zero. The alternative hypothesis is that there is a linear relationship, and b is not equal to zero.

H0: b = 0
Ha: b ^= 0

Let us see whether there is a linear relationship between oil prices per barrel and gold prices. The dataset is in my book’s GitHub account with the name: sp_oil_gold.xlsx, https://github.com/Apress/Learn-Data-Science-Using-Python/blob/main/Datasets/Chapter%207/sp_oil_gold.xlsx

There are several libraries that you can use in Python to perform linear regression. Two of the popular libraries are statsmodels and scikit-learn. The choice between using Ordinary Least Squares (OLS) from the statsmodels library and linear regression from scikit-learn (sklearn) depends on your specific needs and the context of your analysis. Both libraries offer linear regression functionalities, but they serve somewhat different purposes and have different features.

If you are primarily interested in statistical analysis, hypothesis testing, and detailed diagnostics, and you don’t need the model for broader machine learning tasks, statsmodels might be a better fit.

If you are more focused on predictive modeling, integration with other machine learning tools, and simplicity, scikit-learn is a solid choice.

In some cases, analysts use both libraries: statsmodels for initial exploration and hypothesis testing, and then scikit-learn for building a predictive model that can be easily integrated into a broader machine learning workflow. In our code, we shall include both scikit-learn and statsmodels, followed by diagnostic plots using seaborn to use the advantages of each one.

In the following example, we shall use the scikit-learn library for the model fitting and the fit diagnostics diagrams.

The steps to do the simple linear regression with its fit diagnostics:

1. Import the required libraries, load the dataset, and explore it.

2. Create the model.

3. Test the model and compute the predictions and accuracy score.

4. Plot five fit diagnostics for the oil price per barrel in the presidential election years:

a. Predicted vs. actual values.

b. Residuals vs. Fitted Values plot.

c. R-Student vs Predicted.

d. R-Student vs Leverage.

e. Cook’s D vs Observation.

The whole code is in the Example Code Folder at GitHub: https://github.com/Apress/Learn-Data-Science-Using-Python/blob/main/Example%20Code/Chapter%207.ipynb. However, in the Jupyter notebook, the code will be divided into cells matching the above steps. Now, we shall explain each cell to code each of the steps mentioned above, with the explanation of the fit diagnostics plots.

1. Import the Libraries, Load the Dataset, and Explore

In this step, we import the required libraries, load the dataset, and explore it. To explore the data more, we plot the multiple pairwise bivariate distributions of the variables.

Listing 7.1.1: Import the required libraries, load the dataset, and explore it

import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression
import matplotlib.pyplot as plt
import seaborn as sns
# importing r2_score module
from sklearn.metrics import r2_score
from sklearn.metrics import mean_squared_error
from sklearn.preprocessing import StandardScaler
import statsmodels.api as sm

df = pd.read_excel(’../Datasets/sp_oil_gold.xlsx’)
print(df.head())
df=df[[’oil price’,’gold_prices’]]
print(df.head())

# plot numerical data as pairs
plt.figure(0)
sns.pairplot(df);

In Listing 7.1.1, we start by importing nine libraries. Pandas, Numpy, and sklearn.preprocessing are for crunching the numbers and matrices. Matplotlib and Seaborn are for visualization. LinearRegression is a machine learning model. The sklearn.metrics libraries are for model evaluation. The statsmodels.api is to perform OLS and access influence statistics.

In Python, you can achieve the same task in multiple ways. Here, we loaded the data as a Data Frame using the pd.read_excel() function because the dataset is an Excel file. However, you could have loaded it as two Numpy arrays, as follows, because the dataset is small.

oil_prices = np.array([1.63, 1.45, 1.32, 1.82, 11.6, 35.52, 28.2, 14.24, 18.44, 20.29, 27.6, 36.05, 94.1, 109.45, 40.76])

gold_prices = np.array([35.27, 35.1, 39.31, 58.42, 124.74, 615, 361, 437, 343.82, 387.81, 279.11, 409.72, 871.96, 1668.98, 1250.74])

To verify that the data has been loaded corrected, we print the header and the first five rows of the dataset using df.head() function. We keep only the oil price and the gold_prices variables because this is a simple linear regression that needs only one independent variable, and we will drop the rest of the columns. We do that by statement:

df=df[[’oil price’,’gold_prices’]]

To explore the data, we print the pairwise bivariate distributions of the variables using sns.pairplot(df) function. To avoid the figures overlaying over each other, we enumerate them using plt.figure() function.

Figure 7.1: Output of the first cell, Listing 7.1.1

Figure 7.1 shows the output of Listing 7.1.1. The figure shows at first the first five rows from all the columns in the dataset. Then it shows the histogram of the oil price and the scatter plot of the oil price on the y-axis with the gold_price on the x-axis. In the second row of the figure distribution, the scatter plot is repeated again but with the gold_price on the y-axis and the oil price on the x-axis and the second row and column is the gold_price histogram. In both scatter plots, there are three clear outliers, which we will discuss the reason causing this phenomenon.

2. Create the Model

In this step, we specify the independent and target variables. Because the dataset is small, we shall skip the training and testing splitting step. We shall use the whole data set for training and testing. Then we create the model and fit the data to it. From the parameter estimates, we shall write the line equation. Plot the fitted line visually overlaid over the variables’ scatter plot.

Listing 7.1.2: Specify the independent and target variables

x=df[’gold_prices’].values.reshape(-1, 1)
y=df[’oil price’].values

# creating an object of LinearRegression class
LR = LinearRegression()

# fitting the data
LR.fit(x,y)

#parameter estimates
print(’intercept:’, LR.intercept_)
print(’slope:’, LR.coef_)
print(’equation of oil prices in the presedential election years=’,LR.intercept_,”+”,LR.coef_[0],”*gold_prices”)

#Draw the fitted line
plt.figure(2)
plt.scatter(x, y)
plt.xlabel(’Gold Prices’)
plt.ylabel(’Oil Price Per Barrel’)
plt.plot(x,LR.intercept_+LR.coef_*x)

In Listing 7.1.2, we start with changing the independent variable, x, from DataFrame to a Numpy by using the .values property. As we mentioned before, there are multiple ways to do the same task in Python. We can use the method DataFrame.to_numpy() instead. The output of this property is a 1D array.

However, the model fitting function requires a 2D array. So, we use the reshape(-1,1) method to do so. The syntax of the reshape function is .reshape(number of rows, number of columns). To include all the rows without specifying a certain value, we use -1. Since that this is a simple linear regression which means that we have one independent variable, so we have 1 as the number of columns.

For the target variable, y, we again change its type from DataFrame to Numpy by using the .values property.

Now, our variables are ready to be fitted into the model. We first create an object of LinearRegression class using LR = LinearRegression(), then we fit the variables using LR.fit(x,y).

From our model properties, we can get the parameter estimates and write the equation. We use them in plotting the line over the scatter plot of the two variables.

Figure 7.2: Output of Listing 7.1.2

Figure 7.2 shows that there is a linear relationship between oil and gold prices. The slope of the line = 0.05965. From the Parameter Estimate, the equation of the line is:

oil_price = 1.98831 + 0.05965 gold_prices

The plot consists of a scatter plot of data with the regression line overlaid over the scatter plot.

3. Test Model and Compute Predictions and Accuracy Score

Finally, we evaluate the model using R-squared, mean squared error, and root mean squared error.

Listing 7.1.3: Test the model and compute the predictions and accuracy score.

#Listing 7.1.3
#Test Model:
y_prediction = LR.predict(x)
print(’y_prediction=’,y_prediction)

# predicting the accuracy score
score=r2_score(y,y_prediction)
print(’r2 socre = ‘,score)
print(’mean_sqrd_error =’,mean_squared_error(y,y_prediction))
print(’root_mean_squared error =’,np.sqrt(mean_squared_error(y,y_prediction)))

In Listing 7.1.3, we computed the y_prediction values to calculate the accuracy of our model, which was implemented by the r2_score. It is a function inside sklearn. metrics module, where the value of r2_score varies between 0 and 100 percent.

Also, we compute the Mean squared error (MSE), the estimator that measures the average of the squares of errors.

Figure 7.3: Output of Listing 7.1.3

Figure 7.3 shows the output of Listing 7.1.3.

4. Fit Diagnostics

In Listing 7.1.4, we create diagnostic plots to assess the model’s performance. In this code:

R-Student vs Predicted: This plot helps identify influential observations that may be driving the residuals.

R-Student vs Leverage: This plot helps identify observations with high leverage, which can disproportionately affect the model.

Cook’s D vs Observation: This plot helps identify influential observations that may have a large impact on the model parameters.

This part of the code uses statsmodels to fit an OLS model, and then it extracts influence statistics such as R-Student, leverage, and Cook’s D. The final part of the code creates diagnostic plots for R-Student vs Predicted, R-Student vs Leverage, and Cook’s D vs Observation using seaborn.

Listing 7.1.4: Fit Diagnostics Plots

#Listing 7.1.4
plt.figure(3)

# Plot diagnostic plots
plt.figure(figsize=(12, 5))
residuals = y - y_prediction

# Scatter plot of predicted vs. actual values
plt.subplot(1, 2, 1)
sns.scatterplot(x=y, y=y_prediction)
plt.title(’Predicted vs. Actual Values’)

# Residuals vs. Fitted Values plot
plt.subplot(1, 2, 2)
sns.scatterplot(x=LR.predict(x), y=residuals)
plt.axhline(y=0, color=’r’, linestyle=’--’)
plt.title(’Residuals vs. Fitted Values’)
plt.tight_layout()
plt.show() # Finalize and render the figure

# Use statsmodels to perform OLS and access influence statistics
X_ols = sm.add_constant(x)
ols_model = sm.OLS(y, X_ols).fit()

# Access influence statistics
influence = ols_model.get_influence()
r_student = influence.resid_studentized_internal
leverage = influence.hat_matrix_diag
cook_d = influence.cooks_distance

# Remove missing values
mask = ~np.isnan(r_student) & ~np.isnan(leverage)
r_student = r_student[mask]
leverage = leverage[mask]

# Plot R-Student vs Predicted
plt.figure(figsize=(12, 4))
plt.subplot(1, 3, 1)
sns.scatterplot(x=LR.predict(x), y=r_student)
plt.axhline(y=0, color=’r’, linestyle=’--’)
plt.title(’R-Student vs Predicted’)

# Plot R-Student vs Leverage
plt.subplot(1, 3, 2)
sns.scatterplot(x=leverage, y=r_student)
plt.axhline(y=0, color=’r’, linestyle=’--’)
plt.title(’R-Student vs Leverage’)

# Plot Cook’s D vs Observation
plt.subplot(1, 3, 3)
#print(len(df.index),len(cook_d))
sns.scatterplot(x=df.index, y=cook_d[0])
plt.xlim(0, 16) 
# Set y-axis limits to have values above and below zero
# Set x-axis ticks to increment by 1
plt.xticks(np.arange(min(df.index), 15, 1))
plt.axhline(y=0, color=’r’, linestyle=’--’)
plt.title(”Cook’s D vs Observation”)
plt.tight_layout()
plt.show()

These plots are based on statistical diagnostics and can provide insights into the influence of individual observations on the model. Adjust the code as needed for your specific analysis and requirements. Let us look at the output and explain them.

Figure 7.4: Output of the Fit Diagnostics Plots

Figure 7.4 shows the Fit Diagnostics for oil_price. The Dependent Variable (TASK) by the Predicted Value plot visualizes variability in the prediction. In this plot, the dots are random, and there is not a pattern that indicates a constant variance of the error. The outliers are clear in this plot.

Residual vs Fitted shows a random pattern of dots above and below the 0 line, which indicates an adequate model. Again, the three outliers are clear here.

The R-Student by Predicted Value plot shows a couple of outside the ±2 limits. A third dot lies in between them but is far away from the rest of the dots.

The R-Student by Leverage plot shows the outliers which have leverage on the calculation of the regression coefficients.

The Cook’s D plot is designed to identify outliers. Here, the three points are so clear with their values. They are at rows 12, 13, and 14. We shall explore these outliers and the reason behind them in the next section.

These diagnostic plots help assess the assumptions and performance of the linear regression model. They include checking for linearity, homoscedasticity (constant variance of residuals), and influential observations. The plots can provide insights into the model’s behavior and identify potential issues that might require further investigation or model refinement.

This snippet is from my book, Learn Data Science Using Python.