Simple linear regression is a model used to determine the influence of a dependent variable on an independent variable. One can use this model to predict a continuous variable. Regression is based on a formula for a straight-line relationship between two or more variables. The closer the observations are to the line, the better the fit.
In this example, we compared several parametric methods (or linear regression methods), namely:
• Ordinary Least Squared Regression (which estimates a relationship between independent variable and dependent variable by minimizing the sum of the squares in the difference between the actual and predicted values).
• Ridge Regression (which allows for regularization of coefficients).
• RidgeCV Regression (which is a ridge regression with built-in cross-validation).
• Lasso Regression (which allows for regularization of coefficients using few parameters).
• Extreme Gradient Boost Regression (which allows gradient boosting of decision trees)
• Random Forest Regression (which uses multiple decision trees and Bootstrap Aggregation)
We determined whether changes in employees’ work experience influences changes in salary. Thereafter, we used employees’ work experience to predict salary. The hypothesis was follows:-
Null hypothesis: There is no significant difference between an employee’s years of work experience and salary.
Alternative hypothesis: There is a significant difference between an employee’s 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.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
import seaborn as sns
sns.set("talk","ticks",font_scale=1,font="sans-serif",color_codes=True)
import scipy.stats as stats
import statsmodels.api as sm
from sklearn import metrics
from sklearn.linear_model import LinearRegression, Lasso, Ridge, RidgeCV
from xgboost import XGBRegressor
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import cross_val_score, GridSearchCV, train_test_split
from sklearn.ensemble import RandomForestRegressor
import warnings
warnings.filterwarnings("ignore")
from pylab import rcParams
plt.rcParams["figure.figsize"] = [10,10]
df = pd.read_csv(r"C:\Users\Tshepo\Desktop\MLAgortihms\Datasets\Salary_Data.csv")
df.head()Load data
We called the pd.read_csv() function to load obtained csv file into a pandas data frame. See underneath table to examine the data structure.
df = pd.read_csv(r"C:\Users\Tshepo\Desktop\MLAgortihms\Datasets\Salary_Data.csv")
df.head()Detect missing values
The visible presence of missing values result in measurement errors that negatively impact conclusions. To studiously avoid such as a dilemma, data was carefully examined to unanimously determine whether there were any missing values. See underneath figure for a visual representation of the possible presence of missing values.
There were no missing values detected in the data.
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"])Above we created x and y arrays.
x – independent variable (YearsExperience)
y – dependent variable (Salary)
Reshape x and y array
Given that most scikit-lean libraries require a 1D array of the dependent variable to be shaped as a 2D array. The code underneath was used to reshape the created vectors.
x = x.reshape(-1,1)
y = y.reshape(-1,1)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)
x_train, x_val, y_train, y_val = train_test_split(x_train,y_train,test_size=0.2,random_state=0)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)Training data
Underneath we finalized different models, namely Ordinary Least Squared, Ridge, RidgeCV, Extreme Gradient Boost and Random Forest Regression.
lm = LinearRegression().fit(x_train,y_train)
ridge = Ridge().fit(x_train,y_train)
ridgecv = RidgeCV().fit(x_train,y_train)
lasso = Lasso().fit(x_train,y_train)
xgb = XGBRegressor().fit(x_train,y_train)
rfr = RandomForestRegressor().fit(x_train,y_train)Cross validation
def get_val_score(model):
scores = cross_val_score(model, x_train, y_train, scoring="r2")
print("CV mean: ", np.mean(scores))
print("CV std: ", np.std(scores))
print("\n")get_val_score(lm)Output
CV mean: 0.8799423889246292
CV std: 0.03213821058156567
get_val_score(lasso)CV mean: 0.8799429197435654
CV std: 0.032139173602634706
get_val_score(ridge)Output
CV mean: 0.8735101932639928
CV std: 0.041566258797297057
get_val_score(ridgecv)Output
CV mean: 0.880009188026896
CV std: 0.03238720294360753
get_val_score(xgb)Output
CV mean: 0.7548196430599866
CV std: 0.21279942011143785
get_val_score(rfr)Output
CV mean: 0.8992202034295529
CV std: 0.0444344358414017
Hyper-parameter tuning (Ordinary Least Squared)
param_gridOLS = {'fit_intercept':[True,False], 'normalize':[True,False], 'copy_X':[True, False]}
grid_modelOLS = GridSearchCV(estimator=lm, param_grid=param_gridOLS, n_jobs=-1)
grid_modelOLS.fit(x_train,y_train)Output
GridSearchCV(cv=’warn’, error_score=’raise-deprecating’, estimator=LinearRegression(copy_X=True, fit_intercept=True, n_jobs=None, normalize=False), iid=’warn’, n_jobs=-1, param_grid={‘copy_X’: [True, False], ‘fit_intercept’: [True, False], ‘normalize’: [True, False]}, pre_dispatch=’2*n_jobs’, refit=True, return_train_score=False, scoring=None, verbose=0)
print("Best score", grid_modelOLS.best_score_, "Best parameters", grid_modelOLS.best_params_)Output
Best score 0.8792397989835684 Best parameters {'copy_X': True, 'fit_intercept': True, 'normalize': True}
Hyper-parameter tuning (Ridge)
alpharidge = [0.0001,0.001,0.01,0.1,1,10,100,1000]
param_gridridge = dict(alpha=alpha1)
grid_modelridge = GridSearchCV(estimator=ridge, param_grid=param_gridridge, n_jobs=-1)
grid_modelridge .fit(x_train,y_train)Output
GridSearchCV(cv=’warn’, error_score=’raise-deprecating’, estimator=Ridge(alpha=1.0, copy_X=True, fit_intercept=True, max_iter=None, normalize=False, random_state=None, solver=’auto’, tol=0.001), iid=’warn’, n_jobs=-1, param_grid={‘alpha’: [0.0001, 0.001, 0.01, 0.1, 1, 10, 100, 1000]}, pre_dispatch=’2*n_jobs’, refit=True, return_train_score=False, scoring=None, verbose=0)
print("Best score", grid_modelridge.best_score_, "Best parameters", grid_modelridge.best_params_)Output
Best score 0.8792397321262045 Best parameters {'alpha': 0.0001}
Hyper-parameter tuning (lasso)
alphalasso = [0.0001,0.001,0.01,0.1,1,10,100,1000]
param_gridlasso = dict(alpha=alpha)
grid_modellasso = GridSearchCV(estimator=lasso, param_grid=param_gridlasso, n_jobs=-1)
grid_modellasso.fit(x_train,y_train)Output
GridSearchCV(cv=’warn’, error_score=’raise-deprecating’, estimator=Lasso(alpha=0.0001, copy_X=True, fit_intercept=True, max_iter=1000, normalize=False, positive=False, precompute=False, random_state=None, selection=’cyclic’, tol=0.0001, warm_start=False), iid=’warn’, n_jobs=-1, param_grid={‘alpha’: [0.0001, 0.001, 0.01, 0.1, 1, 10, 100, 1000]}, pre_dispatch=’2*n_jobs’, refit=True, return_train_score=False, scoring=None, verbose=0)
print("Best score", grid_modellasso.best_score_, "Best parameters", grid_modellasso.best_params_)Best score 0.8792397989069052 Best parameters {'alpha': 0.0001}
Hyper-parameter tuning (Random Forest Regression)
param_gridrfr = {"n_estimators": [2,4,6,8,10,12,14], "max_depth": [2,4,6,8], "min_samples_leaf": [1,2,3,4,5]}
grid_modelrfr = GridSearchCV(estimator=rfr, param_grid=param_gridrfr, n_jobs=-1)
grid_modelrfr.fit(x_train,y_train)Output
GridSearchCV(cv=’warn’, error_score=’raise-deprecating’, estimator=RandomForestRegressor(bootstrap=True, criterion=’mse’, max_depth=None, max_features=’auto’, max_leaf_nodes=None, min_impurity_decrease=0.0, min_impurity_split=None, min_samples_leaf=1, min_samples_split=2, min_weight_fraction_leaf=0.0, n_estimators=10, n_jobs=None, oob_score=False, random_state=None, verbose=0, warm_start=False), iid=’warn’, n_jobs=-1, param_grid={‘max_depth’: [2, 4, 6, 8], ‘min_samples_leaf’: [1, 2, 3, 4, 5], ‘n_estimators’: [2, 4, 6, 8, 10, 12, 14]}, pre_dispatch=’2*n_jobs’, refit=True, return_train_score=False, scoring=None, verbose=0)
print("Best score", grid_modelrfr.best_score_, "Best parameters", grid_modelrfr.best_params_)Output
Best score 0.9138512380928171 Best parameters {'max_depth': 4, 'min_samples_leaf': 1, 'n_estimators': 14}
Hyper-parameter tuning (Extreme Gradient Boost Regression)
param_gridxgb = {'nthread':[2,4,6,8,10], 'objective':['reg:linear'],'learning_rate': [0.0001,0.001,0.01,0.1,0.2,0.4,0.6,0.8], 'max_depth': [2,4,6,8],'min_child_weight': [2,4,6,8],'silent': [1],'subsample': [0.7],'colsample_bytree': [0.0001,0.001,0.01,0.1,0.2,0.4,0.6,0.8],'n_estimators': [2,4,6,8,10,12,14]}
grid_modelxgb = GridSearchCV(estimator=xgb, param_grid=param_gridxgb, n_jobs=-1)
grid_modelxgb.fit(x_train,y_train)GridSearchCV(cv=’warn’, error_score=’raise-deprecating’, estimator=XGBRegressor(base_score=0.5, booster=’gbtree’, colsample_bylevel=1, colsample_bynode=1, colsample_bytree=1, gamma=0, importance_type=’gain’, learning_rate=0.1, max_delta_step=0, max_depth=3, min_child_weight=1, missing=None, n_estimators=100, n_jobs=1, nthread=None, objective=’reg:linear’, random_sta… param_grid={‘colsample_bytree’: [0.0001, 0.001, 0.01, 0.1, 0.2, 0.4, 0.6, 0.8], ‘learning_rate’: [0.0001, 0.001, 0.01, 0.1, 0.2, 0.4, 0.6, 0.8], ‘max_depth’: [2, 4, 6, 8], ‘min_child_weight’: [2, 4, 6, 8], ‘n_estimators’: [2, 4, 6, 8, 10, 12, 14], ‘nthread’: [2, 4, 6, 8, 10], ‘objective’: [‘reg:linear’], ‘silent’: [1], ‘subsample’: [0.7]}, pre_dispatch=’2*n_jobs’, refit=True, return_train_score=False, scoring=None, verbose=0)
print("Best score", grid_modelxgb.best_score_, "Best parameters", grid_modelxgb.best_params_)Output
Best score 0.8837578939122723 Best parameters {‘colsample_bytree’: 0.0001, ‘learning_rate’: 0.4, ‘max_depth’: 2, ‘min_child_weight’: 2, ‘n_estimators’: 14, ‘nthread’: 2, ‘objective’: ‘reg:linear’, ‘silent’: 1, ‘subsample’: 0.7}
Configuring hyper-parameters for all regression models
lm = LinearRegression(copy_X= True, fit_intercept=True, normalize=True).fit(x_train,y_train)
ridge = Ridge(alpha=0.0001).fit(x_train,y_train)
ridgecv = RidgeCV().fit(x_train,y_train)
lasso = Lasso(alpha=0.0001).fit(x_train,y_train)
xgb = XGBRegressor(colsample_bytree = 0.0001, learning_rate=0.4, max_depth=2, min_child_weight=2, n_estimators=14, nthread=2, objective="reg:linear", silent=1, subsample=0.7).fit(x_train,y_train)
rfr = RandomForestRegressor(max_depth=4, min_samples_leaf=1, n_estimators=14).fit(x_train,y_train)Predictions
We passed on the model using the predict() function to make predictions. Thereafter, we created a data frame consisting of predicted values.
y_predlm = lm.predict(x_test)
y_predlm = pd.DataFrame(y_predlm, columns = ["Ordinary Least Squared"])
y_predridge = ridge.predict(x_test)
y_predridge = pd.DataFrame(y_predridge, columns = ["Ridge"])
y_predridgecv = ridgecv.predict(x_test)
y_predridgecv = pd.DataFrame(y_predridgecv, columns = ["RidgeCV"])
y_predlasso = lasso.predict(x_test)
y_predlasso = pd.DataFrame(y_predlasso, columns = ["Lasso"])
y_predxgb = xgb.predict(x_test)
y_predxgb = pd.DataFrame(y_predxgb, columns = ["Extreme Gradient Boost"])
y_predrfr = rfr.predict(x_test)
y_predrfr= pd.DataFrame(y_predrfr, columns = ["Random Forest"])
prediction1 = pd.concat([y_predlm,y_predridge],axis=1)
prediction2 = pd.concat([y_predridgecv,y_predlasso],axis=1)
prediction3 = pd.concat([prediction1,prediction2],axis=1)
prediction4 = pd.concat([y_predxgb,y_predrfr],axis=1)
prediction5 = pd.concat([prediction3,prediction4],axis=1)
prediction5| Ordinary Least Squared | Ridge | RidgeCV | Lasso | Extreme Gradient Boost | Random Forest | |
|---|---|---|---|---|---|---|
| 0 | 38773.916883 | 38774.105204 | 38961.253223 | 38773.917042 | 43122.152344 | 46893.619048 |
| 1 | 125228.013663 | 125227.746965 | 124962.710753 | 125228.013438 | 109790.335938 | 111207.928571 |
| 2 | 64317.172750 | 64317.226634 | 64370.774766 | 64317.172795 | 56123.136719 | 57209.486565 |
| 3 | 62352.306914 | 62352.371139 | 62416.196186 | 62352.306968 | 56123.136719 | 57209.486565 |
| 4 | 117368.550319 | 117368.324987 | 117144.396432 | 117368.550129 | 109790.335938 | 111207.928571 |
| 5 | 109509.086976 | 109508.903009 | 109326.082111 | 109509.086821 | 109790.335938 | 109360.928571 |
Actual
Underneath are actual values, which can be compared with predicated values above.
y_test = pd.DataFrame(y_test, columns = ["Actual"])
y_test| Actual | |
|---|---|
| 0 | 37731.0 |
| 1 | 122391.0 |
| 2 | 57081.0 |
| 3 | 63218.0 |
| 4 | 116969.0 |
| 5 | 109431.0 |
y_test.mean()Actual 84470.166667
dtype: float64
MAElm = metrics.mean_absolute_error(y_test,y_predlm)
MAEridge = metrics.mean_absolute_error(y_test,y_predridge)
MAEridgecv = metrics.mean_absolute_error(y_test,y_predridgecv)
MAElasso = metrics.mean_absolute_error(y_test,y_predlasso)
MAExgb = metrics.mean_absolute_error(y_test,y_predxgb)
MAErfr = metrics.mean_absolute_error(y_test,y_predrfr)
MAEmodel = [[MAElm,MAEridge,MAEridgecv,MAElasso,MAExgb,MAErfr]]
MAEmodeldata = pd.DataFrame(MAEmodel, columns = ("Ordinary Least Squared", "Ridge", "RidgeCV", "Lasso", "Extreme Gradient Boost", "Random Forest"), index=["MAE"]).transpose()
MSElm = metrics.mean_squared_error(y_test,y_predlm)
MSEridge = metrics.mean_squared_error(y_test,y_predridge)
MSEridgecv = metrics.mean_squared_error(y_test,y_predridgecv)
MSElasso = metrics.mean_squared_error(y_test,y_predlasso)
MSExgb = metrics.mean_squared_error(y_test,y_predxgb)
MSErfr = metrics.mean_squared_error(y_test,y_predrfr)
MSEmodel = [[MSElm,MSEridge,MSEridgecv,MSElasso,MSExgb,MSErfr]]
MSEmodeldata = pd.DataFrame(MSEmodel, columns = ("Ordinary Least Squared", "Ridge", "RidgeCV", "Lasso", "Extreme Gradient Boost", "Random Forest"), index=["MSE"]).transpose()
RMSElm = np.sqrt(MSElm)
RMSEridge = np.sqrt(MSEridge)
RMSEridgecv = np.sqrt(MSEridgecv)
RMSElasso = np.sqrt(MSElasso)
RMSExgb = np.sqrt(MSExgb)
RMSErfr = np.sqrt(MSErfr)
RMSEmodel = [[RMSElm,RMSEridge,RMSEridgecv,RMSElasso,RMSExgb,RMSErfr]]
RMSEmodeldata = pd.DataFrame(RMSEmodel, columns = ("Ordinary Least Squared", "Ridge", "RidgeCV", "Lasso", "Extreme Gradient Boost", "Random Forest"), index=["RMSE"]).transpose()
R2lm = metrics.r2_score(y_test,y_predlm)
R2ridge = metrics.r2_score(y_test,y_predridge)
R2ridgecv = metrics.r2_score(y_test,y_predridgecv)
R2lasso = metrics.r2_score(y_test,y_predlasso)
R2xgb = metrics.r2_score(y_test,y_predxgb)
R2rfr = metrics.r2_score(y_test,y_predrfr)
R2model = [[R2lm,R2ridge,R2ridgecv,R2lasso,R2xgb,R2rfr]]
R2modeldata = pd.DataFrame(R2model, columns = ("Ordinary Least Squared", "Ridge", "RidgeCV", "Lasso", "Extreme Gradient Boost", "Random Forest"), index=["R2"]).transpose()
modelevaluation1 = pd.concat([MAEmodeldata, MSEmodeldata],axis=1)
modelevaluation2 = pd.concat([RMSEmodeldata, R2modeldata],axis=1)
modelevaluation3 = pd.concat([modelevaluation1, modelevaluation2],axis=1)
modelevaluation3| MAE | MSE | RMSE | R2 | |
|---|---|---|---|---|
| Ordinary Least Squared | 2076.572280 | 1.040228e+07 | 3225.256599 | 0.990403 |
| Ridge | 2076.489277 | 1.040217e+07 | 3225.239535 | 0.990403 |
| RidgeCV | 2028.976146 | 1.032545e+07 | 3213.323730 | 0.990474 |
| Lasso | 2076.572209 | 1.040228e+07 | 3225.256584 | 0.990403 |
| Extreme Gradient Boost | 5597.090495 | 4.845970e+07 | 6961.299992 | 0.955293 |
| Random Forest | 5385.638889 | 4.638805e+07 | 6810.876974 | 0.957204 |
Conclusion
As seen on the above, RidgeCV is the best regression model for the following reasons:-
- There is a low average magnitude of errors in predictions without considering their direction (at 2028.976146).
- There is a low squares errors or estimated values and predicted values (at 1.032545e+07).
- 99.04% of the variation in data is explained by the model.
This comes as no surprise. RidgeCV is built-in with generalized cross-validation, which is a form of efficient Leave-One-Out cross-validation.
Given the above, we will proceed with our analysis using the RidgeCV Regression.
Test data
Underneath we plotted our test data to gain insight on predicted salary.
plt.scatter(x_test,y_predridgecv,alpha=0.8)
plt.plot(x_test,y_predridgecv,color="red",alpha=0.8)
plt.title("YearsExperience vs Salary - Test data")
plt.xlabel("YearsExperience")
plt.ylabel("Salary")
plt.show()
Observation are tightly aligned to the straight line.
Training data
plt.scatter(x_train,y_train,alpha=0.8)
plt.plot(x_test,y_predridgecv,color="red",alpha=0.8)
plt.title("YearsExperience vs Salary - Training data")
plt.xlabel("YearsExperience")
plt.ylabel("Salary")
plt.show()
Above we plotted the training data and found that observations are close to straight line – not as tight as the test data. The test data might have exaggerated some points.
Actual values vs Predicted values
Underneath we compared predicted values against actual values.
plt.scatter(y_test,y_predridgecv,alpha=0.8)
plt.axhline(color="red",alpha=0.8)
plt.title("Actual Salary vs Predicted Salary")
plt.xlabel("Actual Salary")
plt.ylabel("Predicted Salary")
plt.show()
The above plot shows that the difference between actual values and predicted values was not large.