Lab-9 (Multiple Linear Regression)

(Intro. to Data Science: IIT Delhi - Abu Dhabi)

In this lab, we will work with real datasets to study important ideas around multiple linear regression:

• fitting and interpreting a multiple linear regression model,
• implementing gradient descent for least squares,
• understanding the bias–variance tradeoff,
• using polynomial regression to capture nonlinearity,
• selecting predictor variables,
• evaluating models using train/test error and cross-validation.

We will use the Diabetes dataset from scikit-learn for much of the discussion:

• Number of observations: 442
• Number of predictors: 10
• Response: a quantitative measure of disease progression one year after baseline
• Predictors: age, sex, body mass index, blood pressure, and six blood serum measurements

This is a real medical dataset, and it is well-suited for multiple linear regression because:

• the response is numerical,
• there are multiple predictors,
• some predictors may be more useful than others,
• the relationship may not be perfectly linear.

Since our target variable is continuous, we will evaluate models using measures such as:

• Mean Squared Error (MSE),
• Root Mean Squared Error (RMSE),
• (R^2),
• cross-validation error.

import numpy as np
import pandas as pd
from matplotlib.pyplot import subplots
from sklearn.datasets import load_diabetes

# We have already seen some of the modules and functions below
# Please find out about those you have not seen before
from sklearn.model_selection import train_test_split, cross_val_score, KFold
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.metrics import mean_squared_error, r2_score
from sklearn.feature_selection import SequentialFeatureSelector

# Let us load the diabetes data and do a high-level inspection
diabetes = load_diabetes(as_frame=True)
X = diabetes.data.copy()
y = diabetes.target.copy()

df = X.copy()
df["target"] = y

print("Shape of X:", X.shape)
print("\nFeature names:")
print(list(X.columns))
print("\nFirst five rows:")
display(df.head())
print("\nSummary statistics:")
display(df.describe())

Shape of X: (442, 10)

Feature names:
['age', 'sex', 'bmi', 'bp', 's1', 's2', 's3', 's4', 's5', 's6']

First five rows:

	age	sex	bmi	bp	s1	s2	s3	s4	s5	s6	target
0	0.038076	0.050680	0.061696	0.021872	-0.044223	-0.034821	-0.043401	-0.002592	0.019907	-0.017646	151.0
1	-0.001882	-0.044642	-0.051474	-0.026328	-0.008449	-0.019163	0.074412	-0.039493	-0.068332	-0.092204	75.0
2	0.085299	0.050680	0.044451	-0.005670	-0.045599	-0.034194	-0.032356	-0.002592	0.002861	-0.025930	141.0
3	-0.089063	-0.044642	-0.011595	-0.036656	0.012191	0.024991	-0.036038	0.034309	0.022688	-0.009362	206.0
4	0.005383	-0.044642	-0.036385	0.021872	0.003935	0.015596	0.008142	-0.002592	-0.031988	-0.046641	135.0

Summary statistics:

	age	sex	bmi	bp	s1	s2	s3	s4	s5	s6	target
count	4.420000e+02	4.420000e+02	4.420000e+02	4.420000e+02	4.420000e+02	4.420000e+02	4.420000e+02	4.420000e+02	4.420000e+02	4.420000e+02	442.000000
mean	-2.511817e-19	1.230790e-17	-2.245564e-16	-4.797570e-17	-1.381499e-17	3.918434e-17	-5.777179e-18	-9.042540e-18	9.293722e-17	1.130318e-17	152.133484
std	4.761905e-02	4.761905e-02	4.761905e-02	4.761905e-02	4.761905e-02	4.761905e-02	4.761905e-02	4.761905e-02	4.761905e-02	4.761905e-02	77.093005
min	-1.072256e-01	-4.464164e-02	-9.027530e-02	-1.123988e-01	-1.267807e-01	-1.156131e-01	-1.023071e-01	-7.639450e-02	-1.260971e-01	-1.377672e-01	25.000000
25%	-3.729927e-02	-4.464164e-02	-3.422907e-02	-3.665608e-02	-3.424784e-02	-3.035840e-02	-3.511716e-02	-3.949338e-02	-3.324559e-02	-3.317903e-02	87.000000
50%	5.383060e-03	-4.464164e-02	-7.283766e-03	-5.670422e-03	-4.320866e-03	-3.819065e-03	-6.584468e-03	-2.592262e-03	-1.947171e-03	-1.077698e-03	140.500000
75%	3.807591e-02	5.068012e-02	3.124802e-02	3.564379e-02	2.835801e-02	2.984439e-02	2.931150e-02	3.430886e-02	3.243232e-02	2.791705e-02	211.500000
max	1.107267e-01	5.068012e-02	1.705552e-01	1.320436e-01	1.539137e-01	1.987880e-01	1.811791e-01	1.852344e-01	1.335973e-01	1.356118e-01	346.000000

9.1 Solving optimization problem

Before starting to use the Python linear regression modules, let's build our own function for fitting the data using the concepts related to solving a convex optimization problem that we discussed in class. For the first exercise, we shall use all 10 numerical predictor variables and fit the entire dataset. We shall predict using the following equation: $\hat{y}_i = \theta_0 + \theta_1 \cdot x_{1, i} + ... + \theta_{10} \cdot x_{10, i}$.

# We already have the predictors in the matrix X. 
# We need to add the intercept column with ones.
X = np.insert(X, 0, np.ones(X.shape[0]), axis=1)
# We shall use the closed-form solution that we derived in class
Z = X.T@X
theta = np.linalg.inv(Z)@X.T@y
theta = np.array(theta)

# Let us calculate the error
residuals = y - X@theta
MSE = np.mean(residuals**2)
print('MSE = ', MSE)
print('theta = ', theta)

MSE =  2859.69634758675
theta =  [ 152.13348416  -10.0098663  -239.81564367  519.84592005  324.3846455
 -792.17563855  476.73902101  101.04326794  177.06323767  751.27369956
   67.62669218]

Next, let's try to solve the optimization problem using gradient descent. We have seen how important scaling of the data is for gradient descent. So, we shall start with scaling the data. We shall use the StandardScaler function in the sklearn.preprocessing module to convert the attributes to standard units.

from sklearn.preprocessing import StandardScaler

# Standardize the training and test predictors for gradient descent
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)
# Let's not forget to set the intercept column back to all ones.
X_scaled[:, 0] = np.ones(X.shape[0])
X_scaled

array([[ 1.        ,  0.80050009,  1.06548848, ..., -0.05449919,
         0.41853093, -0.37098854],
       [ 1.        , -0.03956713, -0.93853666, ..., -0.83030083,
        -1.43658851, -1.93847913],
       [ 1.        ,  1.79330681,  1.06548848, ..., -0.05449919,
         0.06015558, -0.54515416],
       ...,
       [ 1.        ,  0.87686984,  1.06548848, ..., -0.23293356,
        -0.98564884,  0.32567395],
       [ 1.        , -0.9560041 , -0.93853666, ...,  0.55838411,
         0.93616291, -0.54515416],
       [ 1.        , -0.9560041 , -0.93853666, ..., -0.83030083,
        -0.08875225,  0.06442552]])

We shall create a function that implements the gradient descent algorithm.

def gradient_descent(Xmat, yvec, lr=0.05, n_iter=2000):
    n, p = Xmat.shape
    theta = np.zeros(p)
    history = []

    for _ in range(n_iter):
        gradient = (2 / n) * Xmat.T @ (Xmat @ theta - yvec)
        theta = theta - lr * gradient
        mse_loss = np.mean((yvec - Xmat @ theta)**2)
        history.append(mse_loss)

    return theta, history

We saw in class how sensitive GD is to learning rate. So, let us try a few learning rates and observe the convergence.

# Try several learning rates
learning_rates = [0.001, 0.01, 0.05, 0.1]
results = {}

for lr in learning_rates:
    theta_hat, history = gradient_descent(X_scaled, y, lr=lr, n_iter=3000)
    results[lr] = (theta_hat, history)

fig, ax = subplots(figsize=(8, 5))
for lr, (_, history) in results.items():
    ax.plot(history, label=f"lr={lr}")
ax.set_xlabel("Iteration")
ax.set_ylabel("Training MSE")
ax.set_title("Gradient descent convergence for different learning rates")
ax.legend()

<matplotlib.legend.Legend at 0x1212ebf10>

# Let us print the MSEs for different learning rates
for lr in learning_rates:
    print('MSE for lr=', lr, 'is', results[lr][1][len(results[lr][1])-1], 'and theta = ',
         results[lr][0])

MSE for lr= 0.001 is 2880.4429471103517 and theta =  [151.75864179  -0.30200546 -11.25543845  25.11925483  15.3004405
  -3.32916817  -5.26002714  -9.74444376   5.52395046  22.22502029
   3.64080951]
MSE for lr= 0.01 is 2867.708933668896 and theta =  [152.13348416  -0.37608323 -11.29493935  24.98004303  15.33158009
 -15.96660941   5.43606419  -4.89225062   5.67190048  27.63822932
   3.29611669]
MSE for lr= 0.05 is 2859.8276494650745 and theta =  [152.13348416  -0.46334361 -11.39253135  24.75884834  15.41678469
 -34.90023404  20.47050174   3.56291498   8.06723443  34.69870704
   3.22689774]
MSE for lr= 0.1 is 2859.6971143759492 and theta =  [152.13348416  -0.47514436 -11.40577141  24.72901716  15.42843976
 -37.46752872  22.5076079    4.71113201   8.39492544  35.65529544
   3.21745503]

Now that we know that even a high learning rate of 0.1 works, we can use it for further analysis.

9.2 Training and test set

Until now, we have tried to fit the entire labeled dataset into our linear regression model. However, this does not serve our key agenda, which is to come up with a prediction model that works well on new observations that are not part of the training set. To get a sense of how our prediction model will behave on a new observation, we split the training data into a training set and a test set. This is not a new idea since we have done this for the k-NN model in the previous lab as well. We'll create an 80:20 train:test split, train the linear regression model using the training set, and then test it on the test set.

# Let us restart with the diabetes dataset
X = diabetes.data.copy()
y = diabetes.target.copy()

# Split the data and fit a baseline multiple linear regression model
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

# Let us scale the dataset using the standard scaler. Even though scaling was crucial for GD,
# scaling is not required when we use the linear regression function that is part of the
# sklearn.linear_model module. 
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)
X_train_scaled = np.insert(X_train_scaled, 0, np.ones(X_train_scaled.shape[0]), axis=1)
X_test_scaled = np.insert(X_test_scaled, 0, np.ones(X_test_scaled.shape[0]), axis=1)
y_train_np = y_train.to_numpy() if hasattr(y_train, "to_numpy") else np.asarray(y_train)
y_test_np = y_test.to_numpy() if hasattr(y_test, "to_numpy") else np.asarray(y_test)

# We shall train our model using GD as well as using the sklearn library
theta_gd, history_gd = gradient_descent(X_train_scaled, y_train_np, lr=0.1, n_iter=3000)
# Let us print the theta vector and the MSE over the test and training sets
print('theta obtained from GD = ', theta_gd)

y_train_pred = X_train_scaled@theta_gd
y_test_pred = X_test_scaled@theta_gd
train_mse = mean_squared_error(y_train, y_train_pred)
test_mse = mean_squared_error(y_test, y_test_pred)
train_r2 = r2_score(y_train, y_train_pred)
test_r2 = r2_score(y_test, y_test_pred)


print("Training MSE:", train_mse)
print("Test MSE:", test_mse)
print("Training R^2:", train_r2)
print("Test R^2:", test_r2)

theta obtained from GD =  [153.73654391   1.7546429  -11.51146822  25.61024833  16.82786538
 -44.28638811  24.5138678    7.60530327  13.11647069  35.09884904
   2.3525401 ]
Training MSE: 2868.5501836770923
Test MSE: 2900.0593062621933
Training R^2: 0.5279193072034936
Test R^2: 0.4526281156302925

Let's now repeat the exercise using library calls to the LinearRegression function. Note that we need not scale the data.

# Let us restart with the diabetes dataset
X = diabetes.data.copy()
y = diabetes.target.copy()

# Split the data and fit a baseline multiple linear regression model
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

# Let us now fit the training data 
linreg = LinearRegression()
linreg.fit(X_train, y_train)

y_train_pred = linreg.predict(X_train)
y_test_pred = linreg.predict(X_test)

train_mse = mean_squared_error(y_train, y_train_pred)
test_mse = mean_squared_error(y_test, y_test_pred)
train_r2 = r2_score(y_train, y_train_pred)
test_r2 = r2_score(y_test, y_test_pred)

print("Training MSE:", train_mse)
print("Test MSE:", test_mse)
print("Training R^2:", train_r2)
print("Test R^2:", test_r2)


# Let us display the coefficients table to see the magnitude of influence of 
# each predictor on the response
coef_table = pd.DataFrame({
    "predictor": X.columns,
    "coefficient": linreg.coef_
}).sort_values("coefficient", key=np.abs, ascending=False)

display(coef_table)

Training MSE: 2868.549702835578
Test MSE: 2900.1936284934827
Training R^2: 0.5279193863361497
Test R^2: 0.45260276297191926

	predictor	coefficient
4	s1	-931.488846
8	s5	736.198859
2	bmi	542.428759
5	s2	518.062277
3	bp	347.703844
7	s4	275.317902
1	sex	-241.964362
6	s3	163.419983
9	s6	48.670657
0	age	37.904021

One thing we observe is that the R^2 for the test data is not as good as that of the training set. This indicates overfitting.

import matplotlib.pyplot as plt
# Residual analysis for the baseline linear model
residuals = y_test - y_test_pred

plt.figure(figsize=(4, 3))
plt.scatter(y_test, y_test_pred, alpha=0.7)
plt.xlabel("Actual target")
plt.ylabel("Predicted target")
plt.title("Actual vs Predicted (test set)")
plt.show()

plt.figure(figsize=(4, 3))
plt.scatter(y_test_pred, residuals, alpha=0.7)
plt.axhline(0, linestyle="--")
plt.xlabel("Predicted target")
plt.ylabel("Residual")
plt.title("Residuals vs Predicted (test set)")
plt.show()

9.3 Bias–variance tradeoff

In this section, we will explore the bias-variance tradeoff, i.e., the idea that:

• simpler models may have high bias and low variance,
• more flexible models may have low bias and high variance.

There are various ways to make the model more complex. One is considering a polynomial dependence rather than a linear dependence. We will use polynomial regression on a single predictor to make the tradeoff visible. The main ideas we explore in this section can also be considered later.

We know the theory for doing polynomial regression from scratch. If we are using two predictor variables $x_1, x_2$, and we do polynomial regression with degree 2, then we add three extra columns each for $x_1^2, x_2^2, (x_1\cdot x_2)$. We then find the coefficients corresponding to $intercept, x_1, x_2, x_1^1, x_2^2$, and $x_1\cdot x_2$ using the same techniques (i.e., closed form solution or GD). Instead of implementing this from scratch, we shall use the sklearn library called PolynomialFeatures that does precisely this. Here is a short write-up:

In scikit-learn, PolynomialFeatures is a preprocessing tool used to generate non-linear features by combining existing variables into polynomial and interaction terms. This allows linear models (like LinearRegression) to fit complex, non-linear data patterns.

PolynomialFeatures

Core Functionality: It transforms an input matrix into a new one containing all polynomial combinations of the features up to a specified degree.

• Example: For input features $[a, b]$ at degree=2, it generates $[1, a, b, a^2, ab, b^2]$.
• Expansion: The number of features grows polynomially with the number of inputs and exponentially with the degree.

Key Parameters:

• degree (default=2): The maximum power for the features. Can be an integer or a tuple (min_degree, max_degree).
• interaction_only (default=False): If set to True, only products of distinct features are produced (e.g., $ab$ is kept, but $a^2$ is excluded).
• include_bias (default=True): Includes a column of ones, which acts as the intercept term in a linear model.

We will use the predictor bmi and s1 as the predictor variables and the diabetes progression score as the response. In the code below, we observe train/test MSE as a function of the polynomial degree. We observe that the train MSE decreases with degree whereas, the test MSE increases.

# Polynomial regression on two predictors
Xb = X[["bmi", "s1"]]
Xb_train, Xb_test, yb_train, yb_test = train_test_split(
    Xb, y, test_size=0.3, random_state=42
)

degrees = range(1, 8)
poly_results = []

for d in degrees:
    model = Pipeline([
        ("poly", PolynomialFeatures(degree=d, include_bias=False)),
        ("scaler", StandardScaler()),
        ("linreg", LinearRegression())
    ])
    model.fit(Xb_train, yb_train)
    yb_train_pred = model.predict(Xb_train)
    yb_test_pred = model.predict(Xb_test)

    poly_results.append({
        "degree": d,
        "train_mse": mean_squared_error(yb_train, yb_train_pred),
        "test_mse": mean_squared_error(yb_test, yb_test_pred)
    })

poly_df = pd.DataFrame(poly_results)
display(poly_df)

plt.figure(figsize=(6, 4))
plt.plot(poly_df["degree"], poly_df["train_mse"], marker="o", label="Train MSE")
plt.plot(poly_df["degree"], poly_df["test_mse"], marker="o", label="Test MSE")
plt.xlabel("Polynomial degree")
plt.ylabel("MSE")
plt.title("Bias–variance tradeoff via polynomial regression on BMI")
plt.legend()
plt.show()

	degree	train_mse	test_mse
0	1	3890.245986	3832.733230
1	2	3868.095804	3898.145924
2	3	3827.598717	3800.167325
3	4	3775.686334	3924.791041
4	5	3742.001920	4367.414162
5	6	3675.593276	4518.595731
6	7	3600.815373	6018.642387

9.4 Cross-validation for model selection

A single train/test split can be misleading. We now use cross-validation to select model complexity more reliably. This should not be a new concept since we covered this in the previous lab.

We will use the same bmi, s1 predictor, and compute 5-fold cross-validation MSE for polynomial degrees: $1,2,3,4,5,6,7,8,9,10$.

# Polynomial regression on two predictors
Xb = X[["bmi", "s1"]]
# Let us keep aside a test set that we can use after model selection 
# using cross-validation
Xb_train, Xb_test, yb_train, yb_test = train_test_split(
    Xb, y, test_size=0.2, random_state=42
)

# Now, we will run the cross-validation on the training set
cv_results = []

for d in range(1, 11):
    model = Pipeline([
        ("poly", PolynomialFeatures(degree=d, include_bias=False)),
        ("scaler", StandardScaler()),
        ("linreg", LinearRegression())
    ])

    scores = cross_val_score(
        model,
        Xb_train,
        yb_train,
        cv = KFold(n_splits=5, shuffle=True, random_state=42),
        scoring="neg_mean_squared_error"
    )

    cv_results.append({
        "degree": d,
        "cv_mse_mean": -scores.mean(),
        "cv_mse_std": scores.std()
    })

cv_df = pd.DataFrame(cv_results)
display(cv_df)

plt.figure(figsize=(6, 4))
plt.plot(cv_df["degree"], cv_df["cv_mse_mean"], marker="o")
plt.xlabel("Polynomial degree")
plt.ylabel("Mean 5-fold CV MSE")
plt.title("Cross-validation for polynomial degree selection")
plt.show()

	degree	cv_mse_mean	cv_mse_std
0	1	3.909985e+03	2.120624e+02
1	2	3.962487e+03	1.911091e+02
2	3	4.033647e+03	3.412430e+02
3	4	5.204778e+03	2.347937e+03
4	5	5.758134e+03	3.246469e+03
5	6	2.837413e+04	3.280303e+04
6	7	1.607620e+06	3.187779e+06
7	8	6.565720e+05	1.150962e+06
8	9	3.363120e+05	5.796066e+05
9	10	3.942730e+09	7.876418e+09

We observe that the the MSE is almost the same for $d=1, 2, 3$ and then there is a significant increase. So, we can pick $d=3$ for evaluation on the test set.

model = Pipeline([
        ("poly", PolynomialFeatures(degree=3, include_bias=False)),
        ("scaler", StandardScaler()),
        ("linreg", LinearRegression())
    ])

model.fit(Xb_train, yb_train)
yb_train_pred = model.predict(Xb_train)
yb_test_pred = model.predict(Xb_test)
print("test_mse: ", mean_squared_error(yb_test, yb_test_pred))

test_mse:  3915.878902705903

9.5 Feature Selection

An important question when doing multiple linear regression is:

Which predictors should we include?

Too few predictors may cause underfitting. Too many may add noise and hurt interpretability. We will try a simple forward sequential feature selection procedure. This is a greedy procedure that starts with no predictors (in forward selection) and greedily adds predictors based on the score (MSE or R^2) while doing cross-validation. A short write-up on the sklearn function is given below:

SequentialFeatureSelector (SFS) in scikit-learn is a greedy feature selection transformer that builds a feature subset by either adding or removing one feature at a time based on a model's cross-validation score.

Core Mechanisms:

• Forward Selection: Starts with zero features and iteratively adds the "best" feature—the one that most improves the model's cross-validated performance—until the target number is reached.
• Backward Selection: Starts with all features and iteratively removes the "least important" feature—the one whose removal least impacts or most improves the score—until the target number remains.

The code below is intuitive enough to make the usage of SequentialFeatureSelector clear.

# Let us restart with the diabetes dataset
X = diabetes.data.copy()
y = diabetes.target.copy()

# Split the data and fit a baseline multiple linear regression model
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

# Forward feature selection
base_estimator = LinearRegression()

sfs3 = SequentialFeatureSelector(
    base_estimator,
    n_features_to_select=3,
    direction="forward",
    #scoring="neg_mean_squared_error",
    scoring = "r2",
    cv=5
)
sfs5 = SequentialFeatureSelector(
    base_estimator,
    n_features_to_select=5,
    direction="forward",
    #scoring="neg_mean_squared_error",
    scoring = "r2",
    cv=5
)

sfs3.fit(X_train, y_train)
sfs5.fit(X_train, y_train)

selected_3 = X.columns[sfs3.get_support()].tolist()
selected_5 = X.columns[sfs5.get_support()].tolist()

print("Selected 3 predictors:", selected_3)
print("Selected 5 predictors:", selected_5)

def fit_and_evaluate(X_train_sub, X_test_sub, y_train, y_test):
    model = LinearRegression()
    model.fit(X_train_sub, y_train)
    pred = model.predict(X_test_sub)
    return {
        "test_mse": mean_squared_error(y_test, pred),
        "test_r2": r2_score(y_test, pred)
    }

all_results = fit_and_evaluate(X_train, X_test, y_train, y_test)
res3 = fit_and_evaluate(X_train[selected_3], X_test[selected_3], y_train, y_test)
res5 = fit_and_evaluate(X_train[selected_5], X_test[selected_5], y_train, y_test)

comparison = pd.DataFrame([
    {"model": "All 10 predictors", **all_results},
    {"model": "Selected 3 predictors", **res3},
    {"model": "Selected 5 predictors", **res5},
])

display(comparison)

Selected 3 predictors: ['bmi', 'bp', 's5']
Selected 5 predictors: ['sex', 'bmi', 'bp', 's3', 's5']

	model	test_mse	test_r2
0	All 10 predictors	2900.193628	0.452603
1	Selected 3 predictors	2891.037211	0.454331
2	Selected 5 predictors	2811.385100	0.469365

# Try polynomial regression on a selected subset
subset = selected_3  # or choose your own list of predictors
poly2_model = Pipeline([
    ("poly", PolynomialFeatures(degree=2, include_bias=False)),
    ("linreg", LinearRegression())
])

poly2_model.fit(X_train[subset], y_train)
poly2_pred = poly2_model.predict(X_test[subset])

print("Selected subset:", subset)
print("Degree-2 polynomial test MSE:", mean_squared_error(y_test, poly2_pred))
print("Degree-2 polynomial test R^2:", r2_score(y_test, poly2_pred))

Selected subset: ['bmi', 'bp', 's5']
Degree-2 polynomial test MSE: 2824.8994726603664
Degree-2 polynomial test R^2: 0.4668141633634232

9.6 Ridge and Lasso

We saw in class that Ridge and Lasso regression allow us to obtain a good bias-variance tradeoff in multiple linear regression using a single parameter $\alpha$, which is called the regularization parameter. Let us check the effectiveness of these techniques for the diabetes dataset.

from sklearn.linear_model import RidgeCV, LassoCV

# Let us restart with the diabetes dataset
X = diabetes.data.copy()
y = diabetes.target.copy()

# Split the data and fit a baseline multiple linear regression model
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

alphas = np.logspace(-4, 4, 100)
cv = KFold(n_splits=5, shuffle=True, random_state=42)

ridge_model = Pipeline([
    ("scaler", StandardScaler()),
    ("ridge", RidgeCV(alphas=alphas, cv=cv))#, scoring="neg_mean_squared_error"))
])

lasso_model = Pipeline([
    ("scaler", StandardScaler()),
    ("lasso", LassoCV(alphas=alphas, cv=cv, random_state=42, max_iter=10000))
])

ridge_model.fit(X_train, y_train)
lasso_model.fit(X_train, y_train)

ridge_pred = ridge_model.predict(X_test)
lasso_pred = lasso_model.predict(X_test)

print("Ridge best alpha:", ridge_model.named_steps["ridge"].alpha_)
print("Ridge RMSE:", np.sqrt(mean_squared_error(y_test, ridge_pred)))
print("Ridge R^2:", r2_score(y_test, ridge_pred))

print("Lasso best alpha:", lasso_model.named_steps["lasso"].alpha_)
print("Lasso RMSE:", np.sqrt(mean_squared_error(y_test, lasso_pred)))
print("Lasso R^2:", r2_score(y_test, lasso_pred))

Ridge best alpha: 0.9111627561154896
Ridge RMSE: 53.782216913587
Ridge R^2: 0.45404982842366015
Lasso best alpha: 0.521400828799969
Lasso RMSE: 53.44045229205649
Lasso R^2: 0.4609663753064137

9.7 Categorical predictors and one-hot encoding

As discussed in class, we must be careful when using categorical predictors. Performing multiple linear regression after assigning arbitrary numerical values to the categories may not give a good predictor. The solution is to use one-hot encoding that effectively allows us to consider the influence of different category classes on the response separately. For this discussion, we shall consider a housing dataset (please download) that has the "city name" as a categorical predictor variable.

housing = pd.read_csv('sfhousing.csv')
data = housing[['city', 'bsqft', 'price']]

# We shall consider data only from four cities
four_cities = ["Berkeley", "Albany", "Piedmont", "Richmond"]
# Retain data for only the above four cities
data1 = data[data['city'].isin(four_cities)]
# Some price entries are non numeric. Changing the column to numeric
data1['price'] = pd.to_numeric(data1['price'], errors='coerce')
# Drop the rows with undefined entries
data2 = data1.dropna().reset_index(drop=True)
# Drop the rows with extremely large value of bsqft
data3 = data2[data2['bsqft'] < 8000].reset_index(drop=True)
data3.shape

Let's do a scatter plot of price versus bsqft and note that the price crucially depends on the city in addition to the bsqft.

B = data3[data3['city'] == 'Berkeley']
L = data3[data3['city'] == 'Albany']
P = data3[data3['city'] == 'Piedmont']
R = data3[data3['city'] == 'Richmond']

fig, ax = subplots(nrows=2, ncols=2, figsize = (10, 8))
ax[0, 0].scatter(B['bsqft'], B['price'])
ax[0, 0].set_title('Berkeley')
ax[0, 1].scatter(L['bsqft'], L['price'])
ax[0, 1].set_title('Albany')
ax[1, 0].scatter(P['bsqft'], P['price'])
ax[1, 0].set_title('Piedmont')
ax[1, 1].scatter(R['bsqft'], R['price'])
ax[1, 1].set_title('Richmond')

Text(0.5, 1.0, 'Richmond')

The above scatter plots suggest that the city does indeed play an important part in determining the price. However, the ordering of the city in which they determine the price is not known. So, it is better to treat the city as a categorical predictor and use one-hot encoding for our linear regression. Please read the manual and review the code below to understand the mechanism.

from sklearn.preprocessing import OneHotEncoder
import numpy as np

# Initialize and fit
encoder = OneHotEncoder(sparse_output=False)
encoded_data = encoder.fit_transform(data3[['city']])

# Get the labels (e.g., 'category_column_value1', 'category_column_value2')
column_labels = encoder.get_feature_names_out(['city'])

# Create a labeled DataFrame
encoded_df = pd.DataFrame(encoded_data, columns=column_labels)
data4 = pd.concat([data3, encoded_df], axis=1)

# We can now work with the following X, y for our linear regression
y = data4['price']
X = data4[['bsqft', 'city_Albany', 'city_Berkeley', 'city_Piedmont', 'city_Richmond']]
X

	bsqft	city_Albany	city_Berkeley	city_Piedmont	city_Richmond
0	1178.0	1.0	0.0	0.0	0.0
1	1736.0	1.0	0.0	0.0	0.0
2	1736.0	1.0	0.0	0.0	0.0
3	1521.0	1.0	0.0	0.0	0.0
4	754.0	1.0	0.0	0.0	0.0
...	...	...	...	...	...
14030	3086.0	0.0	0.0	0.0	1.0
14031	2079.0	0.0	0.0	0.0	1.0
14032	2079.0	0.0	0.0	0.0	1.0
14033	2732.0	0.0	0.0	0.0	1.0
14034	2079.0	0.0	0.0	0.0	1.0

14035 rows × 5 columns

9.8 Project on Multiple Linear Regression (100 points)

You are asked to use all the concepts covered in this lab to design a predictor for a smaller version (with 30000 points) of the above housing dataset (also see this link). This dataset has been created for you in the code cell below.

X = housing[['county', 'city', 'zip', 'street', 'br', 'lsqft', 'bsqft', 'year', 'date']]
y = housing['price']

# create a smaller subset of 30,000 examples
X_small, _, y_small, _ = train_test_split(
    X, y, train_size=30000, random_state=42
)
# For this project you are asked to consider (X_small, y_small) as the entire labeled dataset available to you.
# Please feel free to create as many cells as required for this project

One of the aspects on which your classifier will be graded is R^2 (on a randomly chosen subset of the housing dataset). You also need to show and explain all your work on how you decided on your final predictor (including any data cleaning operations). Your understanding of various issues will determine your score on this minor project. You may create as many code/markdown cells to explain your work clearly. At the end, you should create a predictor named myPredictor(X) that takes as input a dataframe $X$ with features ('county', 'city', 'zip', 'street', 'br', 'lsqft', 'bsqft', 'year', 'date') and outputs a price prediction vector.