Lab 6: Generalized Linear Models

Course: INS-605: Data Analysis II
Lecturer: Sothea HAS, PhD

The jupyter notebook can be downloaded here.

# These are the libraries that you will need for this lab
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats, optimize
import statsmodels.api as sm
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score, mean_absolute_error, mean_poisson_deviance

1. Multiple Linear Regression: Boston Housing

Goal: predict MEDV (median home value, in $1000) and improve the model with feature engineering.

1.1. Data & Exploration

Run the cell below to load the data and keep the 7 features used in the lecture.

url = 'https://raw.githubusercontent.com/selva86/datasets/master/BostonHousing.csv'
try:
    data = pd.read_csv(url)
except Exception:                       # fallback if no internet access to GitHub
    from sklearn.datasets import fetch_openml
    data = fetch_openml(name='boston', version=1, as_frame=True).frame
data.columns = data.columns.str.upper()
data.head()
CRIM ZN INDUS CHAS NOX RM AGE DIS RAD TAX PTRATIO B LSTAT MEDV
0 0.00632 18.0 2.31 0 0.538 6.575 65.2 4.0900 1 296.0 15.3 396.90 4.98 24.0
1 0.02731 0.0 7.07 0 0.469 6.421 78.9 4.9671 2 242.0 17.8 396.90 9.14 21.6
2 0.02729 0.0 7.07 0 0.469 7.185 61.1 4.9671 2 242.0 17.8 392.83 4.03 34.7
3 0.03237 0.0 2.18 0 0.458 6.998 45.8 6.0622 3 222.0 18.7 394.63 2.94 33.4
4 0.06905 0.0 2.18 0 0.458 7.147 54.2 6.0622 3 222.0 18.7 396.90 5.33 36.2

A. Identify quantative and qualitative columns. Make sure they are in correct data type.

  • Print summary statistics of all inputs according to their types (detect problems such as invalid data).
  • Plot distribution of each feature. Notice features with outliers without handling them.
  • Print summary statistics and plot the distribution of MEDV.
# TODO: understand the inputs and target

B. Detect missing values, duplicaicates, and handle them.

# TODO: handle missing values and drop the duplicated data

C. Compute the Pearson and Spearman correlation of each feature with MEDV, and the gap \(|\rho_{\text{Spearman}}|-|\rho_{\text{Pearson}}|\).

corr = pd.DataFrame({
    'Pearson': None,     # TODO
    'Spearman': None     # TODO: df.corr(method='spearman')
}).drop('MEDV')
corr['gap'] = None       # TODO
corr.round(2)

D. In a markdown cell, briefly answer: which feature is the most correlated with MEDV? Which features have a clearly larger Spearman than Pearson, and what does it suggest?

D. Your answer here:

…

1.2. Feature Engineering

A. For each feature \(x\), compute the Pearson correlation of \(f(x)\) with MEDV for \(f\in\{x,\ \log x,\ \sqrt x,\ x^2\}\). Which transformations look promising?

transforms = {'x': lambda x: x, 'log(x)': np.log, 'sqrt(x)': np.sqrt, 'x^2': np.square}

table = {}
for c in cols:
    table[c] = {name: None for name, f in transforms.items()}   # TODO: Pearson corr. of f(df[c]) with MEDV
pd.DataFrame(table).T.round(3)

B. Create at least 3 new features in df_fe (transformation, polynomial, interaction, binning, …) inspired by the table above and your own intuition. Then check their correlation with MEDV.

Make sure there is no NaN/inf in the new columns.

df_fe = df.copy()
# TODO: add at least 3 new columns, e.g.
# df_fe['fe1'] = ...
# df_fe['fe2'] = ...

print(np.isfinite(df_fe).all().all())
df_fe.corr()['MEDV'].round(2)

1.3. Preprocessing: Split & Standardize

A. Complete prepare: split the data (80/20), then fit the scaler on the training set only and apply it to both sets.

def prepare_data(data_in, features, test_size=0.2, seed=605):
    X_tr, X_te, y_tr, y_te = train_test_split(
        data_in[features], data_in['MEDV'], test_size=test_size, random_state=seed)
    scaler = StandardScaler()
    # TODO: fit the scaler on X_tr only
    Xs_tr = pd.DataFrame(None, columns=features, index=X_tr.index)   # TODO: scaled X_tr
    Xs_te = pd.DataFrame(None, columns=features, index=X_te.index)   # TODO: scaled X_te
    return Xs_tr, Xs_te, y_tr, y_te, scaler

B. In one sentence: why do we fit the scaler on the training set only?

B. Your answer here:

…

1.4. Model & Performance

A. Complete evaluate (RMSE, RMSE/Y_mean and \(R^2\)).

B. Fit three different LinearRegression models:

  • full_lm: using all columns of the original dataset.
  • sel_lm: using only the selected features from correlations + your gut feeling.
  • fe_lm: using all columns of df_fe except MEDV.
def evaluate(y_true, y_pred):
    rmse = None   # TODO
    r2 = None     # TODO
    return {'RMSE': None, 'R2': None, 'RMSE/mean(y)': None} # TODO: replace None with actual values

feat_base = None # TODO: full model
feat_sel = None # TODO: selected feature names
feat_fe = [c for c in df_fe.columns if c != 'MEDV']

results = {}
for name, feats in {'Full': feat_base, 'Selected': feat_sel, 'Engineered': feat_fe}.items():
    Xs_tr, Xs_te, ytr, yte, sc = prepare(df_fe, feats)
    lm = LinearRegression()
    lm.fit(None, None)   # TODO
    results[(name, 'train')] = evaluate(ytr, lm.predict(Xs_tr))
    results[(name, 'test')] = None   # TODO
pd.DataFrame(results).T.round(3)

C. Does feature engineering help on the test set? Is there any sign of overfitting?

C. Your answer here:

…

1.5. Test & Interpretation

A. Fit the OLS on the full standardized training data using statsmodels and read the summary.

# You should define a function to add the constant 1 to your input matrix:
add_c = lambda X: sm.add_constant(X, has_constant='add')   # add the intercept column

Xs_tr, Xs_te, ytr, yte, sc = prepare(df_fe, feat_base)
ols = None   # TODO: sm.OLS(y, X).fit()  (add the intercept to input here)
print(ols.summary())
coef_tab = pd.concat([
    ols.params.rename('coef'),
    ols.conf_int().rename(columns={0: '2.5%', 1: '97.5%'}),
    ols.pvalues.rename('p-value')], axis=1)
coef_tab.round(3)

B. Coefficients above are per 1 std. Convert them to per 1 raw unit (\(\beta_j^{\text{raw}}=\beta_j/\sigma_j\), with \(\sigma_j\) = sc.scale_).

beta_raw = None   # TODO
pd.DataFrame({'per 1 std': ols.params.iloc[1:], 'per 1 raw unit': beta_raw}).round(3)

C. In a markdown cell, briefly answer:

  1. Which coefficients are not significant at the 5% level?

  2. Interpret \(\beta_{\text{RM}}\) and \(\beta_{\text{LSTAT}}\) (in dollars, holding the others fixed).

  3. Can you see explain why we don’t want to use highly correlated features in our model?

  4. What is the null hypothesis of the \(F\)-test reported in the summary?

C. Your answer here:

  1. …

  2. …

  3. …

# (optional) verify your interpretation here

2. Poisson Regression: Student Absences

Data: Student Alcohol Consumption (Kaggle, student-mat.csv): 395 secondary-school students. Target: absences, the number of school absences (a count).

2.1. Data & Exploration

stu = pd.read_csv('student-mat.csv')
print(stu.shape)
stu.head()

A. Compute the mean, the variance and the proportion of zeros of absences, and plot its distribution.

y = stu['absences']
mean_y = None     # TODO
var_y = None      # TODO
prop_zero = None  # TODO
print(f'mean = {mean_y:.2f} | variance = {var_y:.2f} | zeros = {prop_zero:.1%}')

# TODO: histogram of absences

B. In a markdown cell, briefly answer: why is a Poisson model a natural choice here? What do you notice when comparing the mean and the variance?

B. Your answer here:

…

2.2. Preprocessing

A. One-hot encode the categorical features (drop_first=True), then split (80/20) and standardize only the numeric columns (fit on train only).

num_cols = ['age', 'studytime', 'failures', 'famrel', 'freetime', 'goout', 'Dalc', 'Walc', 'health', 'G3']
cat_cols = ['sex', 'address', 'famsize', 'Pstatus', 'schoolsup', 'famsup', 'paid',
            'activities', 'internet', 'romantic', 'higher']

X = None   # TODO: numeric columns + one-hot encoded categorical columns (pd.get_dummies(..., drop_first=True, dtype=float))
y = stu['absences']

Xp_tr, Xp_te, yp_tr, yp_te = None, None, None, None   # TODO: train_test_split(X, y, test_size=0.2, random_state=605)

scaler_p = None   # TODO: fit a StandardScaler on the numeric columns of Xp_tr
Xp_tr_s, Xp_te_s = Xp_tr.copy(), Xp_te.copy()
Xp_tr_s[num_cols] = None   # TODO
Xp_te_s[num_cols] = None   # TODO

print(X.shape)
Xp_tr_s.head()

B. Why drop_first=True? Why don’t we standardize the dummy variables?

B. Your answer here:

…

2.3. Model

A. From scratch. With \(\widehat y_i=e^{x_i^T\beta}\), the Poisson negative log-likelihood (up to constants) is \[\text{Poisson loss}(\beta)=\frac1n\sum_{i=1}^n\left[e^{x_i^T\beta}-y_i\,x_i^T\beta\right],\qquad \nabla_\beta=\frac1n X^T\left(e^{X\beta}-y\right).\] Complete the loss and gradient, then minimize with scipy.optimize.minimize.

Xtr_np = add_c(Xp_tr_s).to_numpy(float)
ytr_np = yp_tr.to_numpy(float)

def poisson_loss(beta, X, y):
    z = X @ beta
    return None   # TODO: write possion loss here

def poisson_grad(beta, X, y):
    z = X @ beta
    return None   # TODO: what's its gradient?

beta0 = np.zeros(Xtr_np.shape[1])
beta0[0] = np.log(ytr_np.mean())          # start from the intercept-only solution
opt = optimize.minimize(poisson_loss, beta0, args=(Xtr_np, ytr_np), jac=poisson_grad, method='BFGS')
beta_scratch = opt.x
print(opt.success, opt.fun)

B. With statsmodels. Fit the Poisson GLM (sm.GLM with sm.families.Poisson() family) and compare with your solution.

pois = None   # TODO: to fit a Poisson GLM with statsmodels
print(pois.summary())
print('max |beta_scratch - beta_statsmodels| =', np.abs(beta_scratch - pois.params.to_numpy()).max())

2.4. Validation

A. Predict the expected counts on the test set and compare the model with the null baseline (always predict the training mean) using the mean Poisson deviance, MAPE and RMSE (compared to target mean):

\[\begin{align*}\text{Poisson Deviance }D&=2\sum_{i=1}^n\left[y_i\log\left(\frac{y_i}{\widehat{y}_i}\right)-(y_i-\widehat{y}_i)\right]\\ \text{MAE}/\overline{Y}&=\frac{1}{n\overline{Y}}\sum_{i=1}^n|y_i-\widehat{y}_i|\\ \text{RMSE}/\overline{Y}&=\frac{1}{n\overline{Y}}\sum_{i=1}^n(y_i-\widehat{y}_i)^2\end{align*}\]

mu_te = None                                         # TODO: predicted expected counts on the test set
mu_null = np.full(len(yp_te), yp_tr.mean())          # baseline

def poisson_scores(y_true, mu):
    return {
        'Poisson deviance': None,   # TODO: mean_poisson_deviance(y_true, mu)
        'MAE/Mean(Y)': None,                # TODO
        'RMSE/Mean(Y)': None,               # TODO
    }

pd.DataFrame({'Null (mean only)': poisson_scores(yp_te, mu_null),
              'Poisson GLM': poisson_scores(yp_te, mu_te)}).T.round(3)

B. Calibration: split the test predictions into 5 quantile bins and compare the mean predicted vs. mean observed count in each bin.

calib = pd.DataFrame({'pred': np.asarray(mu_te), 'obs': yp_te.to_numpy()})
calib['bin'] = pd.qcut(calib['pred'], 5, duplicates='drop')
calib_tab = None   # TODO: mean of 'pred' and 'obs' within each bin (groupby)

m = calib_tab.to_numpy().max()
plt.plot(calib_tab['pred'], calib_tab['obs'], 'o-', color='#345a8b')
plt.plot([0, m], [0, m], 'k--')
plt.xlabel('Mean predicted'); plt.ylabel('Mean observed'); plt.title('Calibration (test)'); plt.show()
calib_tab.round(2)

C. Check the dispersion: \(\displaystyle\widehat\phi=\frac{\chi^2_{\text{Pearson}}}{\text{df}_{\text{resid}}}\) on the training set using pois.pearson_chi2 and pois.df_resid (should be \(\approx1\) for a Poisson model).

dispersion = None   # TODO
print(f'dispersion = {dispersion:.2f}')

D. In a markdown cell, briefly answer: does the model beat the baseline? What does \(\widehat\phi\) tell you about the reliability of the standard errors and p-values?

D. Your answer here:

…

2.5. Test & Interpretation

A. Build a table of incidence-rate ratios \(\text{IRR}_j=e^{\widehat\beta_j}\) with their 95% CI (\(e^{\text{CI}(\beta_j)}\)) and p-values, and draw a forest plot.

ci = pois.conf_int()
irr = pd.DataFrame({'IRR': None,      # TODO: exp of coefficients
                    'lower': None,    # TODO: exp of lower CI bound
                    'upper': None,    # TODO
                    'p-value': pois.pvalues}).drop('const').sort_values('p-value')
display(irr.round(3))

fig, ax = plt.subplots(figsize=(6, 7))
ax.errorbar(irr['IRR'], range(len(irr)), xerr=None,   # TODO: [IRR - lower, upper - IRR]
            fmt='o', capsize=3, color='#345a8b')
ax.axvline(1, color='gray', linestyle='--')
ax.set_yticks(range(len(irr))); ax.set_yticklabels(irr.index); ax.invert_yaxis()
ax.set_xlabel('IRR (95% CI)'); plt.tight_layout(); plt.show()

B. In a markdown cell, interpret the two most significant effects as “% change in the expected number of absences” (numeric features: per +1 std; dummies: vs. the reference level). Is any result surprising?

B. Your answer here:

…

# (optional) verify your interpretation here

C. Test the global significance with a likelihood-ratio test against the intercept-only model: \[\text{LR}=D_{\text{null}}-D_{\text{model}}\ \sim\ \chi^2_{p}\quad\text{under }H_0:\beta_1=\dots=\beta_p=0,\] where \(D\) is the deviance and \(p\) the number of slopes. Report the p-value and conclude.

lr_stat = None   # TODO: pois.null_deviance - pois.deviance
df_lr = None     # TODO: number of slopes or coefficients excluding the intercept
p_lr = None      # TODO: stats.chi2.sf(lr_stat, df_lr)
print(f'LR = {lr_stat:.1f}, df = {df_lr}, p-value = {p_lr:.3g}')

C. Your answer here:

…

3. Going Further (Homework)

The sections below are for you to complete after class.

3.1. MLR: Diagnostics & Do the New Features Really Help?

A. For the engineered OLS model, plot residuals vs. fitted values and a QQ-plot of the residuals. Are the assumptions (\(e_i\sim\mathcal N(0,\sigma^2)\): centered, constant variance, normal) reasonable?

Xf_tr, Xf_te, ytr, yte, _ = prepare(df_fe, feat_fe)
ols_f = sm.OLS(ytr, add_c(Xf_tr)).fit()
fitted, resid = ols_f.fittedvalues, ols_f.resid

fig, ax = plt.subplots(1, 2, figsize=(10, 4))
ax[0].scatter(None, None, color='#345a8b', alpha=.6)   # TODO: fitted vs. residuals
ax[0].axhline(0, color='r', linestyle='--'); ax[0].set_xlabel('Fitted'); ax[0].set_ylabel('Residual')
# TODO: QQ-plot of the residuals with sm.qqplot(resid, line='s', ax=ax[1])
plt.tight_layout(); plt.show()

A. Your answer here:

…

B. The baseline model is nested in the engineered one. Use the \(F\)-test (compare_f_test) to test whether the new features jointly improve the fit: \(H_0\): all coefficients of the new features are \(0\).

Xb_tr, _, _, _, _ = prepare(df_fe, feat_base)
ols_b = sm.OLS(ytr, add_c(Xb_tr)).fit()

F, p_val, df_diff = None, None, None   # TODO: ols_f.compare_f_test(ols_b)
print(f'F = {F:.2f}, df diff = {df_diff:.0f}, p-value = {p_val:.3g}')

B. Your answer here:

…

3.2. Poisson: Fixing Overdispersion

A. Quasi-Poisson: refit with .fit(scale='X2') (the covariance is multiplied by \(\widehat\phi\)). Compare the standard errors and the number of significant features with the classical Poisson.

Xc_tr, Xc_te = add_c(Xp_tr_s), add_c(Xp_te_s)
quasi = None   # TODO: sm.GLM(...).fit(scale='X2')

print('SE ratio (quasi / Poisson):', (quasi.bse / pois.bse).round(3).unique(), '| sqrt(dispersion):', round(np.sqrt(dispersion), 3))
print('# significant (5%) - Poisson:', (pois.pvalues.drop('const') < .05).sum(),
      '| quasi-Poisson:', (quasi.pvalues.drop('const') < .05).sum())

B. Negative Binomial regression allows \(\mathbb V[Y]>\mathbb E[Y]\). Fit sm.NegativeBinomial(y, X).fit(maxiter=500, disp=0), compare its test performance with the Poisson GLM (reuse poisson_scores), and report the estimated dispersion parameter alpha.

nb = None   # TODO: sm.NegativeBinomial(y, X).fit(maxiter=500, disp=0)
mu_nb = None   # TODO: predicted means on the test set (as an array)

print(f"alpha = {nb.params['alpha']:.3f}")
pd.DataFrame({'Poisson': poisson_scores(yp_te, mu_te),
              'NegBin': poisson_scores(yp_te, mu_nb)}).T.round(3)

C. In a markdown cell: how do the conclusions change after accounting for overdispersion (which effects stay significant)? Which model would you report?

C. Your answer here:

…

Wrap-up

In this lab you:

  • Explored, engineered features, standardized (without leakage), fitted and tested an MLR, and interpreted its coefficients (Section 1).
  • Derived the Poisson loss and minimized it by hand, matched it with statsmodels, validated the model on a test set and interpreted incidence-rate ratios (Section 2).
  • (At home) Checked diagnostics, used the \(F\)-test for nested models, corrected for overdispersion, and built a bootstrap CI for a GLM coefficient (Section 3).

The core idea: choose the distribution and the link that match the target, estimate \(\beta\) by maximum likelihood, then always validate, test, and check the assumptions.