CS5100 Model Map ⌂ Modules
Stage 6 of 18Week 2 · Regressionabout 40 min

Multiple regression

Many inputs at once: one formula, $\hat\beta = (\mathbf X^T\mathbf X)^{-1}\mathbf X^T\mathbf y$, covers them all. Then the two tricks that let a linear model use categories and interactions.

Book: what to read for this stage

ISLP 3.2.1 (multiple regression, including the newspaper/radio correlation example): read. 3.2.2 (F-statistic, variable selection): skip. 3.3.1 qualitative predictors and 3.3.2 interactions: read. 3.3.3 potential problems: skim non-linearity, correlated errors, non-constant variance and collinearity; skip outliers, leverage and VIF.

Slide typo

Eqs. (48)-(49), the interaction model, should read \(Y = \beta_0 + (\beta_{TV} + \beta_3X_{Radio})X_{TV} + \beta_{Radio}X_{Radio} + \epsilon\): the interaction term needs its own coefficient \(\beta_3\).

1. Many inputs, one model

\[Y = \beta_0 + \beta_1X_1 + \dots + \beta_pX_p + \epsilon.\]

\(\beta_j\) is the average change in \(Y\) for a one-unit increase in \(X_j\), holding all the other inputs fixed.

Why not run separate simple regressions? They give no single prediction, and each one ignores the other inputs, which misleads when inputs are correlated. The slides’ example: alone, Newspaper advertising gets slope 0.055 and looks useful. With Radio in the model too, Newspaper’s slope drops to about \(-0.001\). Newspaper spend is correlated with Radio spend, so on its own it was taking credit for Radio’s effect.

A fitted model is \(\widehat{\text{Sales}} = 3 + 0.05\,\text{TV} + 0.19\,\text{Radio}\) (budgets in £1000s, sales in 1000s of units). Predict sales for TV = 100 and Radio = 20 (in 1000s of units).

\(3 + 0.05 \times 100 + 0.19 \times 20 = 3 + 5 + 3.8 = 11.8\), i.e. 11,800 units. Interpretation of 0.19: holding TV fixed, each extra £1000 on radio goes with about 190 more units sold, on average.

2. The matrix formula

Put a 1 in front of every input vector, \(\tilde x_i = (1, x_{i1}, \dots, x_{ip})\), so the intercept is just another coefficient. Stack these as the rows of the \(n \times (p + 1)\) design matrix \(\mathbf X\) (first column all ones). Then

\[\hat{\mathbf y} = \mathbf X\beta, \qquad \mathrm{RSS}(\beta) = (\mathbf y - \mathbf X\beta)^T(\mathbf y - \mathbf X\beta).\]

Setting the gradient \(\nabla\mathrm{RSS} = -2\mathbf X^T(\mathbf y - \mathbf X\beta)\) to zero gives the normal equations \(\mathbf X^T\mathbf X\hat\beta = \mathbf X^T\mathbf y\), so

\[\boxed{\hat\beta = (\mathbf X^T\mathbf X)^{-1}\mathbf X^T\mathbf y}.\]

The inverse exists when the columns of \(\mathbf X\) are linearly independent: no input is an exact combination of others, and \(n \ge p + 1\). For a \(2 \times 2\) matrix:

\[\begin{pmatrix}a & b\\ c & d\end{pmatrix}^{-1} = \frac{1}{ad - bc}\begin{pmatrix}d & -b\\ -c & a\end{pmatrix}.\]

In code, don’t invert: use np.linalg.solve(X.T @ X, X.T @ y) or np.linalg.lstsq(X, y).

Worked example: same answer as stage 5. \(x = (1, 2, 3, 4, 5)\), \(y = (2, 4, 5, 4, 5)\); rows of \(\mathbf X\) are \((1, x_i)\).

\(\mathbf X^T\mathbf X = \begin{pmatrix} n & \sum x_i\\ \sum x_i & \sum x_i^2\end{pmatrix} = \begin{pmatrix}5 & 15\\ 15 & 55\end{pmatrix}\) and \(\mathbf X^T\mathbf y = \begin{pmatrix}\sum y_i\\ \sum x_iy_i\end{pmatrix} = \begin{pmatrix}20\\ 66\end{pmatrix}\).

Determinant \(5 \times 55 - 15 \times 15 = 50\), so \(\hat\beta = \frac{1}{50}\begin{pmatrix}55 & -15\\ -15 & 5\end{pmatrix}\begin{pmatrix}20\\ 66\end{pmatrix} = \frac{1}{50}\begin{pmatrix}110\\ 30\end{pmatrix} = \begin{pmatrix}2.2\\ 0.6\end{pmatrix}\). ✓

For \(x = (0, 1, 2, 3)\), \(y = (1, 3, 2, 6)\): \(\mathbf X^T\mathbf X = \begin{pmatrix}4 & 6\\ 6 & 14\end{pmatrix}\) and \(\mathbf X^T\mathbf y = \begin{pmatrix}12\\ 25\end{pmatrix}\). Solve for \(\hat\beta_1\).

Determinant \(56 - 36 = 20\). \(\hat\beta = \frac{1}{20}\begin{pmatrix}14 & -6\\ -6 & 4\end{pmatrix}\begin{pmatrix}12\\ 25\end{pmatrix} = \frac{1}{20}\begin{pmatrix}18\\ 28\end{pmatrix} = \begin{pmatrix}0.9\\ 1.4\end{pmatrix}\).

You have \(n = 50\) observations and keep adding inputs until \(p = 60\). What happens to \(\hat\beta = (\mathbf X^T\mathbf X)^{-1}\mathbf X^T\mathbf y\)?

With more parameters than observations there are infinitely many perfect fits. Stage 8’s ridge penalty is one fix.

3. Categories as inputs: dummy variables

To use “owns a house?” as an input, create \(x_i = 1\) for owners, 0 otherwise:

\[y_i = \beta_0 + \beta_1x_i + \epsilon_i:\quad \text{non-owners } \beta_0,\quad \text{owners } \beta_0 + \beta_1.\]

So \(\beta_0\) is the baseline group’s mean and \(\beta_1\) is the difference between the groups. With \(m\) levels, use \(m - 1\) dummies; the level without a dummy is the baseline. Using \(m\) dummies would make them sum to the intercept column (perfect collinearity), and \(\mathbf X^T\mathbf X\) would be singular. In pandas: pd.get_dummies(df['region'], drop_first=True).

Worked example. Mean balances: East 500, West 550, South 620. With East as the baseline: \(\beta_0 = 500\), \(\beta_{West} = 50\), \(\beta_{South} = 120\). With South as the baseline: \(\beta_0 = 620\), \(\beta_{East} = -120\), \(\beta_{West} = -70\). Predictions don’t change, only the interpretation.

Adding Age as well, \(\widehat{\text{balance}} = 300 + 50\,\text{West} + 120\,\text{South} + 4\,\text{Age}\) gives three parallel lines: same slope for age, different intercepts. A 40-year-old in the South: \(300 + 120 + 160 = 580\).

A predictor has 5 levels. How many dummy variables do you include alongside the intercept?

\(m - 1 = 4\). The fifth level is the baseline.

4. Interactions

Linear regression assumes additivity: the effect of \(X_1\) doesn’t depend on \(X_2\). If radio advertising makes TV advertising work better, add a product column:

\[Y = \beta_0 + \beta_1X_1 + \beta_2X_2 + \beta_3X_1X_2 + \epsilon.\]

Now the effect of one more unit of \(X_1\) is \(\beta_1 + \beta_3X_2\), which depends on \(X_2\). It is still linear in the \(\beta\)’s: \(X_1X_2\) is just another column of \(\mathbf X\), and the same least-squares formula fits it.

The model’s other assumptions (stage 7 has more): uncorrelated errors with constant variance. Check them with a plot of residuals against fitted values: a shapeless band is healthy, a U-shape means a missed non-linearity, a funnel means non-constant variance.

\(\hat Y = 10 + 3X_1 - 2X_2 + 1.5X_1X_2\). By how much does \(\hat Y\) change when \(X_2\) increases by 1, if \(X_1 = 2\)?

The effect of \(X_2\) is \(-2 + 1.5X_1 = -2 + 3 = 1\).

Stage test

In a multiple regression, what does \(\hat\beta_2 = -0.8\) mean?

“Holding the others fixed” is the key phrase. A coefficient is not a correlation.

\(\mathbf X^T\mathbf X = \begin{pmatrix}5 & 15\\ 15 & 55\end{pmatrix}\), \(\mathbf X^T\mathbf y = \begin{pmatrix}20\\ 66\end{pmatrix}\). Give \(\hat\beta_0\).

\(\frac{1}{50}(55 \times 20 - 15 \times 66) = \frac{110}{50} = 2.2\).

Why can’t you include a dummy for every one of the \(m\) levels of a categorical input, plus an intercept?

Perfect collinearity: no unique solution.

\(\widehat{\text{income}} = 20 + 15\,\text{Graduate} + 4\,\text{Age} + 0.5\,\text{Graduate}\times\text{Age}\) (in £1000s; Graduate is 0/1). Predict for a 30-year-old graduate.

\(20 + 15 \times 1 + 4 \times 30 + 0.5 \times 1 \times 30 = 20 + 15 + 120 + 15 = 170\). For graduates the age slope is \(4 + 0.5 = 4.5\); for non-graduates it is 4.

Alone, Newspaper has a positive slope for Sales; in a model with Radio it is about zero. The best explanation:

The same story appears in week 3 (the student paradox in logistic regression).

Is \(Y = \beta_0 + \beta_1X_1 + \beta_2X_2 + \beta_3X_1X_2 + \epsilon\) a linear model?

“Linear” refers to the parameters. Least squares fits it with the same formula.