Simple linear regression is a model used to determine the considerable influence of a continuous or categorical independent variable on a continuous dependent variable. This model finds the extent to which changes in an independent variable influence changes in a dependent variable. The key assumption is that the data retain a linear structure, and a formula of a straight line can be used to derive insight and make predictions. Certain assumptions must be fulfilled to draw meaningful conclusions and make accurate predictions. The assumptions of linear regression are outlined underneath.
- Normality.
- Linearity normal distribution of error term.
- Constant variance of error term.
- Multicollinearity.
These assumptions will be discussed in detail in the next sections. Now let us look at what we are trying to solve using this model. There dataset consists of employees’ information. We looked at how an employee’s work experience influences their monthly remuneration. Obtained data consists of two columns:
- YearsExperience number of years an employee has working.
- Salary amount of money an employee receives monthly for services rendered.
Hypothesis development
Our aim was to accurately predict employees’ guaranteed salary using their years of work experience. Consequently, the hypothesis underneath was constructed.
Null hypothesis: There is no significant difference between employees’ years of work experience and salary.
Alternative hypothesis: There is a significant difference between an employees’ years of work experience and salary.
Import Python libraries
Underneath, key Python libraries were imported. These libraries gallantly helped us create arrays and profile tables, compute graphs and make accurate predictions.
Code
import numpy as np
import pandas as pd
import warnings
warnings.filterwarnings("ignore")
import seaborn as sns
sns.set("talk","ticks",font_scale=1,font="sans-serif",color_codes=True)
import matplotlib.pyplot as plt
%matplotlib inline
from pylab import rcParams
plt.rcParams["figure.figsize"] = [10,10]
import scipy.stats as stats
import statsmodels.api as sm
from sklearn import metrics
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_splitLoad data
df = pd.read_csv(r"C:\Users\Tshepo\Desktop\MLAgortihms\Datasets\Salary_Data.csv")
df.head()We called the pd.read_csv() function to load obtained csv file into a pandas data frame. See the table underneath to examine the data structure.
| YearsExperience | Salary | |
|---|---|---|
| 0 | 1.1 | 39343.0 |
| 1 | 1.3 | 46205.0 |
| 2 | 1.5 | 37731.0 |
| 3 | 2.0 | 43525.0 |
| 4 | 2.2 | 39891.0 |
Detect missing values
The visible presence of missing values results in measurement errors that negatively impact conclusions. To studiously avoid such a dilemma, data must be examined to check whether there are any missing values. The figure underneath depicts detected missing values.
sns.heatmap(df.isnull(),cmap="Blues")
plt.title("Detect Missing Values")
plt.show()
There were no missing values detected in the data.
Pearson correlation
When examining two variables, one is interested in understanding the correlation between those variables. The conventional method of finding the correlation coefficient (a measure of the considerable strength of a linear relationship) is the Pearson correlation method.
The Pearson correlation method has values between -1 and +1
Where
- -1 indicates that there is a negative linear relationship between two variables.
- 0 indicates that there is no linear correlation between two variables.
- 1 indicates that there is a positive linear correlation between two variables.
dfcorr = df.corr(method="pearson")
sns.heatmap(dfcorr, annot=True,annot_kws={"size":12}, cmap="Blues")
plt.title("Pearson Correlation Matrix")
plt.show()
The above figure convincingly depicts a strong positive linear relationship between employees’ work experience and salary. The obtained value was 0.98, which is closer to 1.
Covariance
Underneath, is a visual representation of the joint variability between employees’ years of work experience and salary. The figure underneath depicts how employees’ years of work experience and salary vary together.
dfcov = df.cov()
sns.heatmap(dfcov, annot=True,annot_kws={"size":12}, cmap="Blues")
plt.title("Covariance Matrix")
plt.show()
The figure above shows that there is a positive covariance. There is joint variability in between employees’ work experience and salary.
Detecting outliers
Outliers represents extreme data points which can result in inflated error rates. In simple terms, an outlier represents any value in the data that is either too small or too large. Underneath are ways of handling missing data:
- Removing outliers from the dataset.
- Replacing outliers with the sample mean or median.
fig, (ax1,ax2) = plt.subplots(1,2,figsize=(16,7))
sns.boxplot(df["YearsExperience"],ax=ax1)
ax1.set_title("YearsExperience Box Plot")
sns.boxplot(df["Salary"],ax=ax2)
ax2.set_title("Salary Box Plot")
plt.show()
There were no outliers found on any of the features on the data. Both plots convincingly show that the data is skewed to the left. Independent observations were not spread symmetrically around the true value. The knowledge of mean and explained variance was not sufficient to succinctly summarize data from this distribution.
Determine distribution
The data follow a normal distribution when independent observations are spread symmetrically around the true value. Standardized normal distribution with different means and variance can be written mathematically in terms of a distribution with mean equals to 0 and standard deviation equals to 1.
fig, (ax1,ax2) = plt.subplots(1,2,figsize=(16,7))
sns.distplot(df["YearsExperience"],ax=ax1)
ax1.set_title("YearsExperience Histogram")
ax1.set_ylabel("Related YearsExperience Frequency")
sns.distplot(df["Salary"],ax=ax2)
ax2.set_title("Salary Histogram")
ax2.set_ylabel("Related Salary Frequency")
plt.show()
The figures above confirm that the data is skewed to the left. There is a need for normalizing data prior to training the simple linear regression model.
Line plot
The line plot underneath shows employees’ years of work experience on the x axis and their salary on the y axis.
df.plot(kind="line")
plt.title("YearsExperience vs Salary Line Plot")
plt.show()
The ordered points of employees’ years of work experience and salary show an upward trend. It is explicit that employees’ salary increases as years of work experience increase.
Density plot
Underneath, a density plot was used to depict the probability density using kernel density estimates.
sns.kdeplot(df["YearsExperience"],df["Salary"])
plt.title("YearsExperience vs Salary Density Plot")
plt.show()
The density plot above gave us an idea that the mean salary was between 50 000 and 80 000 and the mean years of work experience was between 2 and 6 years.
Using the code underneath, we unpacked detailed descriptive findings.
df.describe().transpose()| count | mean | std | min | 25% | 50% | 75% | max | |
YearsExperience | 30.0 | 5.313333 | 2.837888 | 1.1 | 3.20 | 4.7 | 7.70 | 10.5 |
| Salary | 30.0 | 76003.000000 | 27414.429785 | 37731.0 | 56720.75 | 65237.0 | 100544.75 | 122391.0 |
Scatter plot
A scatter plot gives an idea of a relationship between two variables. The independent variable lies on the x-axis and the dependent variable lies on the y-axis. This plot depicts how an independent variable affects a dependent variable (correlation). The closer the data points are to the straight-line, the higher the correlation between those variables. The figure underneath depicts the linear relationship between employees’ work experience and salary.
sns.jointplot(x="YearsExperience",y="Salary",data=df,height=10,kind="reg")
plt.show()
The scatter plot above depicts that there is a positive slope. The straight line perfectly fitted the points. There is a strong positive correlation relationship between employees’ work experience and salary. A formula for a straight-line relationship can be constructed.
Create x and y array
The code underneath allowed us to create vectors of x and y using the NumPy library.
x = np.array(df["YearsExperience"])
y = np.array(df["Salary"])x – independent variable (YearsExperience)
y – dependent variable (Salary)
Reshape x and y array
Most scikit-learn libraries function require two dimensional arrays. Underneath x and y vectors were reshaped.
x = x.reshape(-1,1)
y = y.reshape(-1,1)Add constant
Intercepts are not included on the StatsModels library. Consequently, we added an intercept using the code underneath.
x_constant = sm.add_constant(x)Develop a regression model using ordinary least squared model and least squared method
The regression model was created using the Ordinary Least Squared (OLS) model and Least Squared Method (LSM). The relationship between employees’ years of work experience and salary was estimated by minimizing the sum of the squares in the difference between actual values and predicted values of the salary configured as a straight line.
model = sm.OLS(y,x_constant).fit()
model.predict(x_constant)
model.summary()| Dep. Variable: | y | R-squared: | 0.957 |
| Model: | OLS | Adj. R-squared: | 0.955 |
| Method: | Least Squares | F-statistic: | 622.5 |
| Date: | Sat, 23 May 2020 | Prob (F-statistic): | 1.14e-20 |
| Time: | 10:06:46 | Log-Likelihood: | -301.44 |
| No. Observations: | 30 | AIC: | 606.9 |
| Df Residuals: | 28 | BIC: | 609.7 |
| Df Model: | 1 | ||
| Covariance Type: | nonrobust |
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
| const | 2.579e+04 | 2273.053 | 11.347 | 0.000 | 2.11e+04 | 3.04e+04 |
| x1 | 9449.9623 | 378.755 | 24.950 | 0.000 | 8674.119 | 1.02e+04 |
| Omnibus: | 2.140 | Durbin-Watson: | 1.648 |
| Prob(Omnibus): | 0.343 | Jarque-Bera (JB): | 1.569 |
| Skew: | 0.363 | Prob(JB): | 0.456 |
| Kurtosis: | 2.147 | Cond. No. | 13 |
About 95.7% of the variation in the data is explained by the model. The adjusted R-square is 95.5%.
The p-value is less than 0.05. We rejected the null hypothesis in favor of the alternative hypothesis. There is a significant difference between employee’s work experience and salary.
Split training and test data
The data was split into training data and test data by properly calling the train_test_split() function. The training data was used to fit the data. Meanwhile, the test data was conveniently used to validate the simple linear regression model. 80% of the data was adequately trained and 20% of the data was for testing purposes.
x_train, x_test, y_train, y_test = train_test_split(x,y,test_size=0.2,random_state=0)LinearRegression(copy_X=True, fit_intercept=True, n_jobs=None, normalize=False)
Normalize data
We called the StandardScaler() function to normalize the data. This function transformed the data in such a way that the mean value was 0 and standard deviation was 1.
scaler = StandardScaler()
x_train = scaler.fit_transform(x_train)
x_test = scaler.transform(x_test)Train data
The code underneath finalizes the simple linear regression model using ordinary least squares.
lm = LinearRegression()
lm.fit(x_train,y_train)Predicted values
The predict() function was called to predict salaries. Thereafter, a data frame consisting of predicted values was created.
y_pred = lm.predict(x_test)
y_pred = pd.DataFrame(y_pred, columns = ["Predicted Salary"])
y_pred| Predicted Salary | |
| 0 | 40748.961841 |
| 1 | 122699.622956 |
| 2 | 64961.657170 |
| 3 | 63099.142145 |
| 4 | 115249.562855 |
| 5 | 107799.502753 |
Actual values
Underneath are actual values that were compared with predicated values above.
y_test = pd.DataFrame(y_test, columns = ["Actual Salary"])
y_test| Actual Salary | |
| 0 | 37731.0 |
| 1 | 122391.0 |
| 2 | 57081.0 |
| 3 | 63218.0 |
| 4 | 116969.0 |
| 5 | 109431.0 |
Mean of actual values
y_test.mean()84470.16666666667
The arithmetic average of actual values is 84470.17
Intercept
lm.intercept_array([73886.20833333])
Coefficient
lm.coef_array([[24053.85556857]])
Regression equation
The regression equation is as follows:-
ŷ = a + bx
ŷ = 73886.21 + 24053.85x
Where a= 73886.21 and b= 24053.85.
In other words, for every additional year of work experience, an employees’ salary increase by 24053.85.
Model evaluation
The table underneath shows the Mean Absolute Error (MAE), Mean Sum of Errors (MSE), Root Mean Sum of Errors (RMSE) and R-Squared (R2) of the simple linear regression model using OLS.
MAE = metrics.mean_absolute_error(y_test,y_pred)
MSE = metrics.mean_squared_error(y_test,y_pred)
RMSE = np.sqrt(MSE)
R2 = metrics.r2_score(y_test,y_pred)
lmmodelevaluation = [[MAE,MSE,RMSE,R2]]
lmmodelevaluationdata = pd.DataFrame(lmmodelevaluation, columns = ("MAE","MSE","RMSE","R2"))
lmmodelevaluationdata| MAE | MSE | RMSE | R2 |
| 2446.172369 | 1.282341e+07 | 3580.979237 | 0.98817 |
MAE
The average magnitude of errors in predictions without considering their direction is 2446.17
MSE
Average squares errors or estimated values and forecasted values is 1.28
RMSE
The RMSE gives an estimate of the variability remaining after a regression relationship has been established. Findings show that the square root of the mean sum of errors is 3580.97.
R2
R2 lies between 0 and 1.
- 1 indicates that the model perfectly explaining the data
- 0 indicates that the model has no explanatory power.
98.17% of the variation in data was explained by the model.
Test data
The figure underneath shows points that the simple linear regression model has predicted.
plt.scatter(x_test,y_pred,alpha=0.8,s=200)
plt.plot(x_test,y_pred,alpha=0.8,color="red")
plt.title("YearsExperience vs Salary - Test data")
plt.xlabel("YearsExperience")
plt.ylabel("Salary")
plt.show()
Independent observations were tightly aligned to the straight line.
Training data
Independent observations were close to straight line. However, there were not as tightly close to the straight line as those of the test data. Test data has a tendency of exaggerating certain points.
Actual values vs predicted values
plt.scatter(y_test,y_pred,alpha=0.8,s=200)
plt.axhline(alpha=0.8,color="red")
plt.title("Actual Salary vs Predicted Salary")
plt.xlabel("YearsExperience")
plt.ylabel("Salary")
plt.show()
The above plot indicates that differences between actual values and predicted values were not large.
Residual diagnosis
Estimates of coefficients of the simple linear regression model are different from 0. The model fits the actual data. Any trend or pattern in residuals would have indicated that the model was a misfit.
residual = y_test - y_pred
model_residual = model.resid
model_fitted = model.fittedvalues
model_leverage = model.get_influence().hat_matrix_diag
model_norm_residual = model.get_influence().resid_studentized_internal
model_norm_residual_ab_sqrt = np.sqrt(np.abs(model_norm_residual))Mean of residuals
np.mean(model_residual)-1.5036979069312415e-11
The arithmetic average after removing actual values and predicted values was -1.50, which was closer to zero.
model_residual = pd.DataFrame(model_residual)
model_residual.columns = ["Residual"]
sns.boxplot(model_residual["Residual"])
plt.title("Residual Box Plot")
plt.show()
sns.distplot(model_residual["Residual"])
plt.title("Residual Histogram")
plt.ylabel("Related Residual Frequency")
plt.show()
Residuals were spread symmetrically around the true value. The knowledge of mean and explained variance was sufficiently to summarize data from this distribution.
Predicted values vs residual values
The figure underneath shows predicted values and residuals (the difference between estimated values and actual values).
plt.scatter(y_pred,residual,alpha=0.8,s=200)
plt.axhline(alpha=0.8,color="red")
plt.title("Predicted Salary vs Residual Salary")
plt.xlabel("Predicted Salary")
plt.ylabel("Residual Salary")
plt.show()
Residuals were randomly distributed.
Normal Q-Q Plot
The simplest way of determining whether a model is a fit or a misfit is by looking at whether residuals follow a straight line.
fig, ax = plt.subplots(figsize=(10,10))
fig = sm.graphics.qqplot(model_residual,ax=ax,line="45",fit=True,dist=stats.norm)
plt.title("Normal Q-Q Plot")
plt.show()
Theoretical quantiles and sample quantiles did not follow a normal distribution.
Cook’s D Influence Plot
Cook’s D influence plot was used to detect outliers on the independent variable.
fig, ax = plt.subplots(figsize=(10,10))
fig = sm.graphics.influence_plot(model,ax=ax,criterion="cooks")
plt.show()
There was an extreme outlier at observation 19 and two minor outliers at observation 28 and observation 29.
Fitted values vs predicted valuesFitted values vs predicted values
plt.scatter(model_fitted,model_residual,alpha=0.8,s=200)
plt.title("Fitted Salary vs Residual Salary")
plt.xlabel("Fitted Salary")
plt.ylabel("Residual Salary")
plt.show()
The figure above highlights characteristics of a well behaved fitted salary vs residual salary plot. Residual salary values bounce randomly around 0.
Leverage values vs residual values
The figure underneath shows extreme values of the independent variable and residuals.
plt.scatter(model_leverage,model_residual,alpha=0.8,s=200)
plt.title("Leverage Salary vs Residual Salary")
plt.xlabel("Leverage Salary")
plt.ylabel("Residual Salary")
plt.show()
The leverage salary vs residual salary plot was well behaved.
Fitted values vs studentized residual values
There was no trend or pattern found in the fitted salary vs studentized residual salary plot.
Leverage values vs studentized residual values
plt.scatter(model_leverage,model_norm_residual,alpha=0.8,s=200)
plt.title("Leverage Salary vs Studentized Residual Salary")
plt.xlabel("Leverage Salary")
plt.ylabel("Studentized Residual Salary")
plt.show()
Studentized residuals take no shape when compared with leverage salary.
Residual lag plot
lag_plot(model_residual["Residual"])
plt.title("Residual Lag Plot")
plt.show()
Lag plot above confirms that residual salary do not follow a normal distribution.
For further residual diagnosis, both auto-correlation and partial auto-correlation plots were constructed to detect any autocorrelations in residual salary.
Residual autocorrelation plot
autocorrelation_plot(model_residual["Residual"])
plt.title("Residual Lag vs Autocorrelation Plot")
plt.show()
Residual salary is centered around the mean which is zero.
Residual ACF plot
plot_acf(model_residual["Residual"])
plt.show()
PACF had a significant spike at lag 1. High order correlation is explained by lag 1 auto-correlation.
Residual PACF plot
plot_pacf(model_residual["Residual"])
plt.show()
PACF confirms that high order correlation is explained by lag 1 auto-correlation.