Lab 05: Model Complexity and Variable Selection

BSAN 750 - Foundations of Machine Learning and AI

This notebook begins with a guided housing-price example and then asks you to repeat the workflow with a real diabetes dataset. You will compare simple and complex models, use AIC/BIC, and tune Ridge and Lasso with cross-validation.

Learning goals

1. Recognize overfitting by comparing training and unseen-data error.
2. Compare candidate regression models using AIC and BIC.
3. Apply stepwise regression and evaluate the final model.
4. Explain how Ridge and Lasso control complexity.
5. Use cross-validation to select alpha without using the testing outcomes.
6. Evaluate the selected model once on held-out testing data.
7. (Optional) Boostrap standard error and confidence interval.

1. Setup and load the housing data

Upload housing_prices.csv if Colab asks you to choose a file. One row represents one home sale. The target is sale_price. We will not use home_id as a predictor because it is only an identifier.

# Import the libraries used throughout the lab.
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.formula.api as smf

from sklearn.datasets import load_diabetes
from sklearn.linear_model import LinearRegression, Ridge, Lasso
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import train_test_split, cross_val_score, GridSearchCV
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler

pd.set_option('display.max_rows', 100)

print("Setup complete.")

# Read the CSV if it is already available; otherwise ask for a Colab upload.
try:
    housing = pd.read_csv('housing_prices.csv')
except FileNotFoundError:
    from google.colab import files
    uploaded = files.upload()
    housing = pd.read_csv('housing_prices.csv')

# Inspect dimensions, variables, missing values, and summary statistics.
print('Rows and columns:', housing.shape)
print('Variables:', housing.columns)
display(housing.head())

# check missing values and obtain summary statistics (skipping this due to long output)
display(housing.isna().sum().to_frame('missing'))
display(housing.describe().round(2))

YOUR TURN - Answer the following questions:

1. What is the target?
2. What does one row represent?
3. Which variables are meaningful candidate predictors?
4. Should you include home_id in subsequent analyses?

Type your answer here.

2. Create training and testing samples

The training sample is used to fit models and make model choices. The testing sample is held aside for final evaluation. We will use the same split for every candidate model so the comparisons are fair.

# change scale of Y to the unit of thousand
housing['sale_price'] = housing['sale_price']/1000

# Split the data into training and testing sets.
housing_train, housing_test = train_test_split(
    housing, test_size=0.20, random_state=751)

print('Training rows:', len(housing_train))
print('Testing rows:', len(housing_test))

3. Compare simple, practical, and overly complex models

We will compare three models:

• model1: two predictors
• model2: five predictors
• model3: eight predictors
• model4: all 25 predictors. Some of them could be noise.

## model specification
model1_spec = ('sale_price ~ square_feet + bedrooms')
print("model1:", model1_spec)
model2_spec = ('sale_price ~ square_feet + bedrooms + home_age '
    '+ distance_to_center + school_rating')
print("model2:", model2_spec)
model3_spec = (
    'sale_price ~ square_feet + bedrooms + home_age '
    '+ distance_to_center + school_rating '
    '+ days_on_market + annual_hoa_fee + property_tax_rate_pct '
)
print("model3:", model3_spec)
model4_predictors = [
    column for column in housing_train.columns
    if column not in ['sale_price', 'home_id']
]
model4_spec = 'sale_price ~ ' + ' + '.join(model4_predictors)
print("model4:", model4_spec)

# Fit the three models
model1 = smf.ols(model1_spec, data=housing_train).fit()
model2 = smf.ols(model2_spec, data=housing_train).fit()
model3 = smf.ols(model3_spec, data=housing_train).fit()
model4 = smf.ols(model4_spec, data=housing_train).fit()

print(model2.summary())

Use cross-validation inside the training sample

A single testing split is useful for final evaluation, but it should not be used repeatedly to choose among models. Five-fold cross-validation makes several training/validation splits inside the training sample. The testing sample stays untouched.

# Compare the four models using their formula specifications directly
model_specs = {
    'Model 1': model1_spec,
    'Model 2': model2_spec,
    'Model 3': model3_spec,
    'Model 4': model4_spec
}

# We also keep a dictionary of the already fitted statsmodels objects to get training RMSE, AIC, and BIC
fitted_models = {
    'Model 1': model1,
    'Model 2': model2,
    'Model 3': model3,
    'Model 4': model4
}

cv_rows = []

for name, formula in model_specs.items():
    # Extract predictors from the formula string (e.g., 'sale_price ~ x1 + x2')
    # Split by '~' and then by '+' to get the individual column names
    rhs = formula.split('~')[1]
    predictors_used = [col.strip() for col in rhs.split('+')]

    # Calculate Cross-Validation RMSE
    scores = cross_val_score(
        LinearRegression(),
        housing_train[predictors_used],
        housing_train['sale_price'],
        cv=5,
        scoring='neg_root_mean_squared_error'
    )

    # Get the already fitted statsmodels object to extract training metrics
    fit = fitted_models[name]
    #training_rmse = np.sqrt(np.mean(fit.resid**2))
    training_rmse = np.sqrt(fit.mse_resid)  # this is adjusted for degrees of freedom, i.e., divided by n-p-1

    cv_rows.append({
        'model': name,
        'model_size': len(predictors_used),
        'training_RMSE': training_rmse,
        'mean_CV_RMSE': -scores.mean(),
        'SD_CV_RMSE': scores.std(),
        'AIC': fit.aic,
        'BIC': fit.bic
    })

model_comparison = pd.DataFrame(cv_rows)
display(model_comparison.round(2))

YOUR TURN - Written response:

1. Which model has the smallest training RMSE?
2. Which model has the smallest CV RMSE?
3. What evidence suggests that the complex model is learning noise?
4. Which model do AIC and BIC prefer?

Type your answer here.

4. Stepwise selection using AIC and BIC

Previously, we manually specified several candidate models. Stepwise selection automates part of this process by repeatedly adding or removing predictors. At each step, it retains the change that produces the lowest AIC or BIC.

All stepwise selection must be conducted using the training sample. The testing sample remains untouched until the final models are evaluated.

## install a package for stepwise regression
%pip install step-criterion==0.2.2 -q
from step_criterion import step_criterion

## forward AIC model
forward_aic = step_criterion(
    data=housing_train,
    initial='sale_price ~ 1',
    scope={
        'lower': 'sale_price ~ 1',
        'upper': model4_spec
    },
    direction='forward',
    criterion='aic',
    #family=sm.families.Binomial(), ## if response is binary
    trace=0
)
print(forward_aic.model.summary())

## forward BIC model
forward_bic = step_criterion(
    data=housing_train,
    initial='sale_price ~ 1',
    scope={
        'lower': 'sale_price ~ 1',
        'upper': model4_spec
    },
    direction='forward',
    criterion='bic',
    #family=sm.families.Binomial(), ## if response is binary
    trace=0
)
print(forward_bic.model.summary())

print('AIC model size:', len(forward_aic.model.params) - 1)
print('AIC model:', forward_aic.model.model.formula)
print('BIC model size:', len(forward_bic.model.params) - 1)
print('BIC model:', forward_bic.model.model.formula)

stepwise_specs = {
    'Forward AIC': [col.strip() for col in forward_aic.model.model.formula.split('~')[1].split('+') if col.strip() != '1'],
    'Forward BIC': [col.strip() for col in forward_bic.model.model.formula.split('~')[1].split('+') if col.strip() != '1']
}

stepwise_cv_rows = []

for name, predictors_used in stepwise_specs.items():
    # Calculate Cross-Validation RMSE
    scores = cross_val_score(
        LinearRegression(),
        housing_train[predictors_used],
        housing_train['sale_price'],
        cv=5,
        scoring='neg_root_mean_squared_error'
    )

    fit = forward_aic.model if name == 'Forward AIC' else forward_bic.model
    training_rmse = np.sqrt(fit.mse_resid)

    stepwise_cv_rows.append({
        'model': name,
        'model_size': len(predictors_used),
        'training_RMSE': training_rmse,
        'mean_CV_RMSE': -scores.mean(),
        'SD_CV_RMSE': scores.std(),
        'AIC': fit.aic,
        'BIC': fit.bic
    })

stepwise_cv_comparison = pd.DataFrame(stepwise_cv_rows)
display(stepwise_cv_comparison.round(2))

YOUR TURN: 1. Which model is better? 2. Does BIC model pick less variables? Why?

So far, we have 4+2 models. Let’s put them together and compare.

# Combine the two cross-validation comparison tables (Models 1-4 and Stepwise models)
combined_cv_comparison = pd.concat(
    [model_comparison, stepwise_cv_comparison],
    ignore_index=True
)

print("Unified Model Cross-Validation Comparison:")
display(combined_cv_comparison.round(2))

YOUR TURN: From above table, do you think there are duplicated models?

print('Model2:', model2.model.formula)
print('BIC model:', forward_bic.model.model.formula)

5. Ridge and Lasso regularization

Ridge and Lasso penalize large coefficients. Ridge usually shrinks every coefficient toward zero. Lasso can set some coefficients exactly to zero. Because the penalty acts on coefficient size, predictors should be standardized first.

A Pipeline ensures that scaling is learned separately inside each cross-validation training fold. GridSearchCV compares candidate values of alpha; larger values mean stronger regularization.

# Prepare the data
X_train = housing_train[model4_predictors]
y_train = housing_train['sale_price']

path_scaler = StandardScaler()
X_path = path_scaler.fit_transform(X_train)

# Setting a sequence of lambda (called alpha in python)
alpha_grid = np.exp(np.linspace(-5, 11, 30))
print(np.log(alpha_grid))

Fitting the ridge regression

We first demonstrate the solution path. Note: this part is NOT necessary in a real workflow.

# fitting ridge regression along the path of lambda
ridge_path = []

for alpha in alpha_grid:
    ridge_model = Ridge(alpha=alpha)
    ridge_model.fit(X_path, y_train)
    ridge_path.append(ridge_model.coef_)

ridge_path = np.asarray(ridge_path)

# showing the path figure (not needed in practice)
for index, predictor in enumerate(model4_predictors):
    plt.plot(
        np.log(alpha_grid),
        ridge_path[:, index],
        label=predictor
    )

plt.xlabel('log(alpha)')
plt.ylabel('Standardized coefficient')
plt.title('Ridge Coefficient Path')
plt.legend(bbox_to_anchor=(1.02, 1), fontsize=7)
plt.show()

Now, use CV to find the best value of \(\lambda\).

This part is what you need in a real workflow. The first cell below is configuration, while the second cell uses the configuration to actually fit the models.

# This part of code is a configuration.
# Create the Ridge pipeline.
ridge_pipeline = Pipeline([
    ('scale', StandardScaler()),
    ('model', Ridge())
])

# Set up five-fold cross-validation for Ridge.
# model__alpha refers to the alpha parameter of the pipeline's model step.
# GridSearchCV selects the alpha with the smallest mean CV RMSE.
ridge_search = GridSearchCV(
    ridge_pipeline,
    {'model__alpha': alpha_grid},
    cv=5,
    scoring='neg_root_mean_squared_error'
)

# This part of code is actual model fitting.
# After selecting alpha, GridSearchCV refits Ridge on the full training sample.
ridge_search.fit(X_train, y_train)

print('Ridge selected alpha:', round(ridge_search.best_params_['model__alpha'], 3))
print('Ridge CV RMSE:', round(-ridge_search.best_score_, 1))

The following cell is to show the CV plot. Not part of workflow in practice.

# pull the solution path from the grid search
ridge_cv = pd.DataFrame({
    'alpha': np.asarray(
        ridge_search.cv_results_['param_model__alpha'], dtype=float
    ),
    'CV_RMSE': -ridge_search.cv_results_['mean_test_score'],
    'SE': ridge_search.cv_results_['std_test_score'] / np.sqrt(5)
}).sort_values('alpha')

# Ridge cross-validation plot
ridge_alpha = ridge_search.best_params_['model__alpha']

plt.errorbar(
    np.log(ridge_cv['alpha']),
    ridge_cv['CV_RMSE'],
    yerr=ridge_cv['SE'],
    marker='o',
    capsize=3
)

plt.axvline(np.log(ridge_alpha), color='red', linestyle='--')
plt.xlabel('log(alpha)')
plt.ylabel('CV RMSE')
plt.title('Ridge Cross-Validation')
plt.show()

Fitting the LASSO regression

Similarly, the first part below is to show the LASSO solution path. It is not necessary in a real workflow.

# set a new lambda sequence for lasso
alpha_grid = np.exp(np.linspace(-5, 6, 30))

# fitting lasso path
lasso_path = []

for alpha in alpha_grid:
    lasso_model = Lasso(alpha=alpha, max_iter=50000)
    lasso_model.fit(X_path, y_train)
    lasso_path.append(lasso_model.coef_)

lasso_path = np.asarray(lasso_path)

# showing the lasso path
for index, predictor in enumerate(model4_predictors):
    plt.plot(
        np.log(alpha_grid),
        lasso_path[:, index],
        label=predictor
    )

plt.xlabel('log(alpha)')
plt.ylabel('Standardized coefficient')
plt.title('LASSO Coefficient Path')
plt.legend(bbox_to_anchor=(1.02, 1), fontsize=7)
plt.show()

Choose optimal \(\lambda\) via cross-validation. This part is what you need in a real workflow.

# Both of the following code are configuration, not actual model fitting
# Create the Lasso pipeline.
# max_iter is increased to help Lasso converge for all candidate alpha values.
lasso_pipeline = Pipeline([
    ('scale', StandardScaler()),
    ('model', Lasso(max_iter=50000))
])

# Set up the same five-fold cross-validation search for Lasso.
lasso_search = GridSearchCV(
    lasso_pipeline,
    {'model__alpha': alpha_grid},
    cv=5,
    scoring='neg_root_mean_squared_error'
)

# After selecting alpha, GridSearchCV refits Ridge on the full training sample.
lasso_search.fit(X_train, y_train)

print('Lasso selected alpha:', round(lasso_search.best_params_['model__alpha'], 3))
print('Lasso CV RMSE:', round(-lasso_search.best_score_, 1))

# pull the solution path from the grid search for Lasso
lasso_cv = pd.DataFrame({
    'alpha': np.asarray(
        lasso_search.cv_results_['param_model__alpha'], dtype=float
    ),
    'CV_RMSE': -lasso_search.cv_results_['mean_test_score'],
    'SE': lasso_search.cv_results_['std_test_score'] / np.sqrt(5)
}).sort_values('alpha')

# Lasso cross-validation plot
lasso_alpha = lasso_search.best_params_['model__alpha']

plt.errorbar(
    np.log(lasso_cv['alpha']),
    lasso_cv['CV_RMSE'],
    yerr=lasso_cv['SE'],
    marker='o',
    capsize=3
)

plt.axvline(np.log(lasso_alpha), color='red', linestyle='--')
plt.xlabel('log(alpha)')
plt.ylabel('CV RMSE')
plt.title('LASSO Cross-Validation')
plt.show()

# Examine the selected Lasso coefficients on the standardized scale.
lasso_coefficients = pd.DataFrame({
    'predictor': model4_predictors,
    'coefficient': lasso_search.best_estimator_.named_steps['model'].coef_
})
lasso_coefficients['selected'] = lasso_coefficients['coefficient'] != 0

display(
    lasso_coefficients
    .assign(abs_coefficient=lambda d: d['coefficient'].abs())
    .sort_values('abs_coefficient', ascending=False)
    .drop(columns='abs_coefficient')
    .round(3)
)
print('Number of nonzero coefficients:', lasso_coefficients['selected'].sum())

YOUR TURN - Written response:

1. How many variables are selected by lasso?
2. The above table shows the coefficients on the standardized scale. What does that mean? (Hint: Compare the coefficient scale to the early OLS models.)
3. Please vibe code to obtain the coefficients on the original scales.
4. Do you think Lasso model is the best one comparing to model 1-4 you had before? What would you do to answer this question?

Type your answer here.

Finally, we make prediction on the testing sample and compare different models.

# Define helper function for RMSE if not already loaded
def rmse(y_true, y_pred):
    return np.sqrt(mean_squared_error(y_true, y_pred))

# Helper to get statsmodels model size (excluding intercept)
def get_sm_size(model_obj):
    return len(model_obj.params) - 1

# Helper to get statsmodels training RMSE
def get_sm_train_rmse(model_obj):
    return np.sqrt(model_obj.mse_resid)

# Extract coefficients from the selected Ridge and Lasso models
ridge_coef = ridge_search.best_estimator_.named_steps['model'].coef_
lasso_coef = lasso_search.best_estimator_.named_steps['model'].coef_

# Create lists of specs for OLS models
ols_models = {
    'Model 1': (model1, [col.strip() for col in model1_spec.split('~')[1].split('+')]),
    'Model 2': (model2, [col.strip() for col in model2_spec.split('~')[1].split('+')]),
    'Model 3': (model3, [col.strip() for col in model3_spec.split('~')[1].split('+')]),
    'Model 4': (model4, model4_predictors),
    'Forward AIC': (forward_aic.model, [col.strip() for col in forward_aic.model.model.formula.split('~')[1].split('+') if col.strip() != '1']),
    'Forward BIC': (forward_bic.model, [col.strip() for col in forward_bic.model.model.formula.split('~')[1].split('+') if col.strip() != '1'])
}

# Collect metrics for all models
comparison_rows = []

for name, (model_obj, predictors) in ols_models.items():
    comparison_rows.append({
        'model': name,
        'model_size': len(predictors),
        'training_RMSE': get_sm_train_rmse(model_obj),
        'testing_RMSE': rmse(housing_test['sale_price'], model_obj.predict(housing_test))
    })

# Append Ridge
comparison_rows.append({
    'model': 'Ridge',
    'model_size': np.sum(np.abs(ridge_coef) > 1e-8),
    'training_RMSE': rmse(y_train, ridge_search.predict(X_train)),
    'testing_RMSE': rmse(housing_test['sale_price'], ridge_search.predict(housing_test[model4_predictors]))
})

# Append Lasso
comparison_rows.append({
    'model': 'Lasso',
    'model_size': np.sum(np.abs(lasso_coef) > 1e-8),
    'training_RMSE': rmse(y_train, lasso_search.predict(X_train)),
    'testing_RMSE': rmse(housing_test['sale_price'], lasso_search.predict(housing_test[model4_predictors]))
})

all_models_comparison = pd.DataFrame(comparison_rows)
display(all_models_comparison.round(2))

YOUR TURN - Written response:

1. Compare Ridge, Lasso, and the unregularized complex model. Which would you recommend from this exercise? Explain why.
2. Read the code in the above cell and try to understand the coding logic.

Type your answer here.

6. Optional: bootstrap coefficient and selection stability

When we train a LASSO model on slightly different samples of our data, does it always select the exact same predictors?

To find out, we can use bootstrap resampling (sampling our rows with replacement) and refit the model many times.

To conduct full boostrap, we retune the best penalty value using cross-validation within each individual resample to capture the full uncertainty of our modeling process.

## This process will take a while (~ 3 mins)
rng = np.random.default_rng(750)
B = 100

# Arrays to store coefficients and the tuned alpha for each resample
estimates_full = np.empty((B, len(model4_predictors)))
selected_alphas = np.empty(B)

print("Starting full bootstrap with retuned alpha...\n")

for b in range(B):
    # Show progress
    if (b + 1) % 10 == 0 or b == 0:
        print(f"Running bootstrap sample {b + 1}/{B}...")

    rows_b = rng.integers(0, len(X_train), size=len(X_train))
    X_boot = X_train.iloc[rows_b]
    y_boot = y_train.iloc[rows_b]

    # We fit a fresh GridSearchCV on this bootstrap sample
    lasso_pipeline_b = Pipeline([
        ('scale', StandardScaler()),
        ('model', Lasso(max_iter=50000))
    ])

    grid_b = GridSearchCV(
        lasso_pipeline_b,
        {'model__alpha': alpha_grid},
        cv=5,
        scoring='neg_root_mean_squared_error'
    )

    grid_b.fit(X_boot, y_boot)

    # Record the best alpha and scaled coefficients converted back to original units
    selected_alphas[b] = grid_b.best_params_['model__alpha']
    estimates_full[b] = grid_b.best_estimator_.named_steps['model'].coef_ / grid_b.best_estimator_.named_steps['scale'].scale_

print("\nFull bootstrap complete!")

# Create stability DataFrame for full bootstrap
stability_full = pd.DataFrame({
    "Bootstrap SD (Full CV)": estimates_full.std(axis=0, ddof=1),
    "Nonzero freq (Full CV)": (estimates_full != 0).mean(axis=0)
}, index=model4_predictors)

# Display summary statistics of the chosen alphas
print(f"Average selected alpha: {selected_alphas.mean():.3f} (SD: {selected_alphas.std(ddof=1):.3f})")

# Plot selection frequencies
ax = stability_full["Nonzero freq (Full CV)"].sort_values().plot.barh(figsize=(8, 6))
ax.set(xlim=(0, 1.05), xlabel="Fraction of bootstrap fits with a nonzero coefficient",
       title="LASSO selection stability (with inner CV tuning)")
plt.tight_layout()
plt.show()