Data Analysis with Linear Models

In this article, we illustrate the IEBAE (pronounced as “eBay”, I/O, Exploration, Benchmark, Analysis, Evaluation) framework with a linear model example.

We use the Graduate Admission 2 dataset, which can be obtained from https://www.kaggle.com/mohansacharya/graduate-admissions. Our ultimate goal is to predict “Chance of Admit” from other variables.

We import some useful packages:

import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
import statsmodels.api as sm
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split, KFold
from sklearn.pipeline import Pipeline
from sklearn.linear_model import Ridge
from sklearn.ensemble import GradientBoostingRegressor

I/O

We can glimpse at the dataset using a tabular form with describe

df = pd.read_csv(“Admission_Predict_Ver1.1.csv”, index_col = “Serial No.”)
df.describe().transpose()

We can save the dataset using HDFS

df.to_hdf("data.h5", "linear_model/admission")

Exploration

The straight forward data exploration technique is pairplot .

sns.pairplot(data = df)

We can see that the chance of admit is highly correlated with GRE/TOFLE, and CGPA score.

Benchmark

Linear models are the de facto baseline models for regression problems. Therefore, we first construct a linear model to predict the Chance of Admit. Before our machine learning analysis, we need to

(1) determine the appropriate metric for evaluating the performance of our models;

(2) design a train-validate-test split strategy.

For the first question, we use the mean squared error between the prediction and true values.

def metric(y1, y2):
return np.sqrt(np.mean((y1-y2)**2))

For the second question, we can use the standard KFold cross-validation for model selection

kf = KFold(n_splits=10)

For model testing, we use the standard 80%-20% train test split

train_idx, test_idx = train_test_split(range(0, df.shape[0]), test_size = 0.2)

Because linear models have strong statistical interpretation, we use the linear regression tools from statsmodels , which provides more interpretable metrics.

sx = StandardScaler()
sy = StandardScaler()
X = sx.fit_transform(df.iloc[train_idx,:7])
y = sy.fit_transform(df.iloc[train_idx,7:8])
X = pd.DataFrame(X, index = df.index[train_idx], columns = df.columns[:7])
X = sm.add_constant(X) # can be omitted
res = sm.OLS(y, X).fit()
res.summary()

We can see that University Rating and SOP do not show statistical significance for predicting chance of admit.

y1 = res.predict(sm.add_constant(sx.transform(df.iloc[train_idx,:7])))
y2 = res.predict(sm.add_constant(sx.transform(df.iloc[test_idx,:7])))
y1 = sy.inverse_transform(y1)
y2 = sy.inverse_transform(y2)
sns.jointplot(x = "True", y = "Prediction", data = pd.DataFrame({"True": df.iloc[train_idx, 7], "Prediction": y1}))
sns.jointplot(x = "True", y = "Prediction", data = pd.DataFrame({"True": df.iloc[test_idx, 7], "Prediction": y2}))

The training and testing scores are 0.0603, 0.0566 respectively.

train_score = metric(y1, df.iloc[train_idx, 7])
test_score = metric(y2, df.iloc[test_idx, 7])

Analysis

We now consider a slightly more advanced (compared with linear regression) model: ridge regression. Because in ridge regression, we have one hyperparameter α, we need to do cross-validation.

def get_score(alpha = 1.0):
pipe = Pipeline([
("scale", StandardScaler()),
("ridge regression", Ridge(alpha = alpha))
])
kf = KFold(n_splits=10)
t_mse = []
v_mse = []
for tid, vid in kf.split(range(X_train.shape[0])):
xt = X_train.iloc[tid, :]
yt = y_train.iloc[tid]
xv = X_train.iloc[vid, :]
yv = y_train.iloc[vid]
pipe.fit(xt, yt)
yhat = pipe.predict(xv)
t_mse.append(np.sqrt(np.mean((yt - pipe.predict(xt))**2)))
v_mse.append(np.sqrt(np.mean((yhat - yv)**2)))
return np.mean(t_mse), np.std(t_mse)

We train with different hyperparameters, and collect the results:

ms = []
ss = []
a = np.linspace(0.0001, 1000.0, 50)
for alpha in a:
m, s = get_score(alpha)
ms.append(m)
ss.append(s)
ms = np.array(ms)
ss = np.array(ss)
plt.plot(a, ms, "g")
plt.fill_between(a, ms - ss, ms + ss, alpha = 0.5)
plt.xscale("log")

In this case, the smaller α, the better. We consider α=10

pipe = Pipeline([
("scale", StandardScaler()),
("ridge regression", Ridge(alpha = 10.0))
])
pipe.fit(X_train, y_train)
yhat = pipe.predict(X_test)
sns.jointplot(x = 'True', y = 'Prediction', data = pd.DataFrame({'True': y_test, 'Prediction': yhat}))

The training and testing scores are 0.0599 and 0.0590, respectively.

The next option is the Gradient Boosting Regression

pipe = Pipeline([
('scale', StandardScaler()),
('gbr', GradientBoostingRegressor())
])
pipe.fit(X_train, y_train)
y1 = pipe.predict(X_train)
y2 = pipe.predict(X_test)
sns.jointplot(x = 'True', y = 'Prediction', data = pd.DataFrame({'True': y_train, 'Prediction': y1}))
sns.jointplot(x = 'True', y = 'Prediction', data = pd.DataFrame({'True': y_test, 'Prediction': y2}))

The training and testing scores are 0.0405 and 0.0594, respectively.

Evaluation

The evaluation step aggregates all the evaluation metrics and presents different models using hbar .

dat = np.array([
[0.06031010091513192, 0.05659089069370097],
[0.05891919602076211, 0.05989847557550637]
])
d = pd.DataFrame(dat, index = ['Linear Regression', 'Ridge Regression'], columns = ['Training', 'Testing'])
d.plot.barh()

Although gradient boosting has a small training error, the testing error is large, indicating a certain degree of overfitting. The ridge regression has the smallest discrepancy between training and testing error, making it the most generalizable model among the three.

Get the Medium app

A button that says 'Download on the App Store', and if clicked it will lead you to the iOS App store
A button that says 'Get it on, Google Play', and if clicked it will lead you to the Google Play store