This homework is due on Gradescope on Friday October 11th at 9pm.
Suppose I have data of the where \(n = 1, \ldots N\) indexes households, \(y_n\) is the expenditure on food in a time period (so that \(y_n > 0\) for all \(n\)). Suppose that households are randomly selected to either receive food stamps (for which \(x_n = 1\)) or to not receive food stamps (for which \(x_n = 0\)). Also, suppose we measure household income \(z_n\).
(a)
Suppose we regress \(y_n \sim \beta_0 + \beta_1 x_n + \beta_2 z_n\). Let \(\hat{f}(x_, z) = \hat{\beta}_0 + \hat{\beta}_1 x+ \hat{\beta}_2 z\) denote the fit.
Using this regression, we might estimate the effect of food stamps on food expenditure by \(\hat{f}(1, z) - \hat{f}(0, z)\). How does this estimate depend on \(z\)?
(b)
Suppose we regress \(\log y_n \sim \gamma_0 + \gamma_1 x_n + \gamma_2 z_n\). Let \(\hat{g}(x_, z) = \hat{\gamma}_0 + \hat{\gamma}_1 x+ \hat{\gamma}_2 z\).
Using this regression, we might estimate the effect of food stamps on food expenditure by \(\exp(\hat{g}(1, z)) - \exp(\hat{g}(0, z))\). How does this estimate depend on \(z\)?
(c)
The regressions in (a) and (b) make different implicit assumptions about how food stamps affect consumption for a particular household. State these assumptions in ordinary language.
Solutions
(a) It doesn’t depend on \(z\):
\[ \hat{f}(1, z) - \hat{f}(0, z) = \hat{\beta}_0 + \hat{\beta}_1 + \hat{\beta}_2 z- (\hat{\beta}_0 + \hat{\beta}_2 z) = \hat{\beta}_1. \]
(b) It is proportional to \(\exp(\hat{\gamma}_2 z)\):
\[ \exp(\hat{g}(1, z)) - \exp(\hat{g}(0, z)) = \exp(\hat{\gamma}_0 + \hat{\gamma}_1 + \hat{\gamma}_2 z) - \exp(\hat{\gamma}_0 + \hat{\gamma}_2 z) = (\exp(\hat{\gamma}_1) - 1) \exp(\hat{\gamma}_0) \exp(\hat{\gamma}_2 z) \]
(c)
In the first model, food stamps have an additive effect on expenditure. In the second model (with the log response), food stamps have a multiplicative effect on expenditure.
Given a regression on \(\boldsymbol{X}\) with \(P\) regressors and \(N\) data points, and the corresponding \(\boldsymbol{Y}\), \(\hat{\boldsymbol{Y}}\), and \(\hat{\varepsilon}\), define the following quantities: \[ \begin{aligned} RSS :={}& \hat{\varepsilon}^\intercal\hat{\varepsilon}& \textrm{(Residual sum of squares)}\\ TSS :={}& \boldsymbol{Y}^\intercal\boldsymbol{Y}& \textrm{(Total sum of squares)}\\ ESS :={}& \hat{\boldsymbol{Y}}^\intercal\hat{\boldsymbol{Y}}& \textrm{(Explained sum of squares)}\\ R^2 :={}& \frac{ESS}{TSS}. \end{aligned} \]
Solutions:
This exercise will use the bodyfat example from the datasets. Suppose we’re interested in predicting bodyfat, which is difficult to measure precisely, with other variables which are easier to measure: Height, Weight, and Abdomen circumference.
If we do so, we get the following sum of squared errors:
reg <- lm(bodyfat ~ Abdomen + Height + Weight, bodyfat_df)
print(reg$coefficients)(Intercept) Abdomen Height Weight
-36.6147193 0.9515631 -0.1270307 -0.1307606 print(sprintf("Error: %f", mean(reg$residuals^2)))[1] "Error: 19.456161"(a)
Noting that Heigh, Weight, and Abdomen are on different scales, your colleague suggests that you might get a better fit by normalizing them. But when you do, here’s what happened:
bodyfat_df <- bodyfat_df %>%
mutate(height_norm=(Height - mean(Height)) / sd(Height),
weight_norm=(Weight - mean(Weight)) / sd(Weight),
abdomen_norm=(Abdomen - mean(Abdomen)) / sd(Abdomen))
reg_norm <- lm(bodyfat ~ abdomen_norm + height_norm + weight_norm, bodyfat_df)
print(reg_norm$coefficients) (Intercept) abdomen_norm height_norm weight_norm
19.150794 10.260778 -0.465295 -3.842945 print(sprintf("Error: %f", mean(reg_norm$residuals^2)))[1] "Error: 19.456161"Our coefficients changed, but our fitted error didn’t change at all.
The new regressors are invertible linear transforms of the original regressors. This means \(\hat{\boldsymbol{Y}}\) does not change, but the coefficients change.
(b)
Chastened, your colleague suggests that maybe it’s the difference between normalized height and weight that would helps us predict. After all, it makes sense that height should only matter relative to weight, and vice versa. So they run the regresson on the difference:
bodyfat_df <- bodyfat_df %>%
mutate(hw_diff = height_norm - weight_norm)
reg_diff <- lm(bodyfat ~ abdomen_norm + height_norm + weight_norm + hw_diff, bodyfat_df)
print(reg_diff$coefficients) (Intercept) abdomen_norm height_norm weight_norm hw_diff
19.150794 10.260778 -0.465295 -3.842945 NA print(sprintf("Error: %f", mean(reg_diff$residuals^2)))[1] "Error: 19.456161"Now, our fitted error didn’t change at all, but our difference coefficient wasn’t even estimated.
R.The new regressors are again invertible linear transforms of the original regressors. The difference is a linear combination of height_norm and weight_norm, so \(\boldsymbol{X}^\intercal\boldsymbol{X}\) is not invertible, and \betavhat is not uniquely defined. R reports NA for coefficients that it cannot estimate.
(c)
Finally, your colleague suggests regressing instead on the ratio of weight to height. Here are the results:
bodyfat_df <- bodyfat_df %>%
mutate(hw_ratio = height_norm / weight_norm)
reg_ratio <- lm(bodyfat ~ abdomen_norm + height_norm + weight_norm + hw_ratio, bodyfat_df)
print(reg_ratio$coefficients) (Intercept) abdomen_norm height_norm weight_norm hw_ratio
19.149475485 10.271751693 -0.478683698 -3.848561820 0.009057317 print(sprintf("Error: %f", mean(reg_ratio$residuals^2)))[1] "Error: 19.423560"Our fitted error is different this time, and we could estimate this coefficient.
hw_ratio, the fitted error decreased. Your colleague tells you that had it been a bad regressor, the error would have increased. Are they correct?You can estimate the ratio because it is a nonlinear transformation of the regressors. With this new regressor, the fit can change. Including another regressor, no matter poorly associated with the response, can never make the error increase. Your colleague is incorrect.
(d)
Let \(x_n = (\textrm{Abdomen}_n, \textrm{Height}_n, \textrm{Weight}_n)^\intercal\) denote our set of regressors.
Your colleague suggests a research project where you improve your fit by regressing \(y_n \sim z_n\) for new regressors \(z_n\) of the form \(z_n = \boldsymbol{A}x_n\), where the matrix \(\boldsymbol{A}\) is chosen using a machine learning algorithm.
Since \(\boldsymbol{A}x_n\) is a linear combination of \(\boldsymbol{x}_n\), it can never give a better fit than regressing on \(x_n\) alone. (It can give a worse fit if \(\boldsymbol{A}\) is lower rank than \(\boldsymbol{x}_n\).) Searching over all possible \(\boldsymbol{A}\) will never give you a better fit than \(y_n \sim \boldsymbol{x}_n\). This is not a good idea for a research project.
(e)
Finally, your colleague suggests a research project where you again regression \(y_n \sim z_n\), but now you let \(z_n = f(x_n)\) for any function \(f\), where you use a neural network to find the best fit to the data over all possible functions \(f(x_n)\).
Since \(f(\boldsymbol{x}_n)\) can be a nonlinear function, you can produce a better fit than regressing on \(\boldsymbol{x}_n\) alone. In fact, the best possible fit among all possible functions is a highly expressive function which fits the data perfectly with zero error. You don’t need a neural network to identify this. Having zero error does not mean the fit is useful for anything — it likely overfits the training data and has poor predictive power on new data. This is not a good idea for a research project.
This homework problem derives a closed-form expression for the effect of leaving a datapoint out of the regression.
We will use the following result, known as the Woodbury formula (but also many other names, including the Sherman-Morrison-Woodbury formula). Let \(A\) denote an invertible matrix, and \(\boldsymbol{u}\) and \(\boldsymbol{v}\) vectors the same length as \(A\). Then
\[ (A + \boldsymbol{u}\boldsymbol{v}^\intercal)^{-1} ={} A^{-1} - \frac{A^{-1} \boldsymbol{u}\boldsymbol{v}^\intercal A^{-1}}{1 + \boldsymbol{v}^\intercal A^{-1} \boldsymbol{u}}. \]
We will also use the definition of a “leverage score” from lecture \(h_n := \boldsymbol{x}_n^\intercal(\boldsymbol{X}^\intercal\boldsymbol{X})^{-1} \boldsymbol{x}_n\). Note that \(h_n = (\boldsymbol{X}(\boldsymbol{X}^\intercal\boldsymbol{X})^{-1} \boldsymbol{X}^\intercal)_{nn}\) is the \(n\)–th diagonal entry of the projection matrix \(\underset{\boldsymbol{X}}{\boldsymbol{P}}\).
Let \(\hat{\boldsymbol{\beta}}_{-n}\) denote the estimate of \(\hat{\beta}\) with the datapoint \(n\) left out. Similarly, let \(\boldsymbol{X}_{-n}\) denote the \(\boldsymbol{X}\) matrix with row \(n\) left out, and \(\boldsymbol{Y}_{-n}\) denote the \(\boldsymbol{Y}\) matrix with row \(n\) left out.
a Prove that
\[ \hat{\boldsymbol{\beta}}_{-n} = (\boldsymbol{X}_{-n}^\intercal\boldsymbol{X}_{-n})^{-1} \boldsymbol{X}_{-n}^\intercal\boldsymbol{Y}_{-n} = (\boldsymbol{X}^\intercal\boldsymbol{X}- \boldsymbol{x}_n \boldsymbol{x}_n^\intercal)^{-1} (\boldsymbol{X}^\intercal\boldsymbol{Y}- \boldsymbol{x}_n y_n) \]
Solution
This follows from \(\boldsymbol{X}^\intercal\boldsymbol{X}= \sum_{n=1}^N\boldsymbol{x}_n \boldsymbol{x}_n^\intercal\) and \(\boldsymbol{X}^\intercal\boldsymbol{Y}= \sum_{n=1}^N\boldsymbol{x}_n y_n\).
b
Using the Woodbury formula, derive the following expression: \[ (\boldsymbol{X}^\intercal\boldsymbol{X}- \boldsymbol{x}_n \boldsymbol{x}_n^\intercal)^{-1} = (\boldsymbol{X}^\intercal\boldsymbol{X})^{-1} + \frac{(\boldsymbol{X}^\intercal\boldsymbol{X})^{-1} \boldsymbol{x}_n \boldsymbol{x}_n^\intercal(\boldsymbol{X}^\intercal\boldsymbol{X})^{-1} }{1 - h_n} \]
Solution
Direct application of the formula with \(\boldsymbol{u}= \boldsymbol{x}_n\) and \(\boldsymbol{v}= -\boldsymbol{x}_n\) gives \[ (\boldsymbol{X}^\intercal\boldsymbol{X}- \boldsymbol{x}_n \boldsymbol{x}_n^\intercal)^{-1} = (\boldsymbol{X}^\intercal\boldsymbol{X})^{-1} + \frac{(\boldsymbol{X}^\intercal\boldsymbol{X})^{-1} \boldsymbol{x}_n \boldsymbol{x}_n^\intercal(\boldsymbol{X}^\intercal\boldsymbol{X})^{-1} }{1 - \boldsymbol{x}_n^\intercal(\boldsymbol{X}^\intercal\boldsymbol{X})^{-1}\boldsymbol{x}_n}. \]
Then recognize the leverage score.
c
Combine (a) and (b) to derive the following explicit expression for \(\hat{\boldsymbol{\beta}}_{-n}\):
\[ \hat{\boldsymbol{\beta}}_{-n} = \hat{\boldsymbol{\beta}}- (\boldsymbol{X}^\intercal\boldsymbol{X})^{-1} \boldsymbol{x}_n \frac{1}{1 - h_n} \hat{\varepsilon}_n. \]
Solution
We have
\[ (\boldsymbol{X}^\intercal\boldsymbol{X}- \boldsymbol{x}_n \boldsymbol{x}_n^\intercal)^{-1} \boldsymbol{X}^\intercal\boldsymbol{Y}= \hat{\beta}+ \frac{(\boldsymbol{X}^\intercal\boldsymbol{X})^{-1} \boldsymbol{x}_n \boldsymbol{x}_n^\intercal\hat{\beta}}{1 - h_n}. \]
and
\[ (\boldsymbol{X}^\intercal\boldsymbol{X}- \boldsymbol{x}_n \boldsymbol{x}_n^\intercal)^{-1} \boldsymbol{x}_n y_n = (\boldsymbol{X}^\intercal\boldsymbol{X})^{-1} \boldsymbol{x}_n y_n + \frac{(\boldsymbol{X}^\intercal\boldsymbol{X})^{-1} \boldsymbol{x}_n h_n }{1 - h_n} y_n. \]
Combining,
\[ \begin{aligned} \hat{\boldsymbol{\beta}}_{-n} ={}& \hat{\beta}+ (\boldsymbol{X}^\intercal\boldsymbol{X})^{-1} \boldsymbol{x}_n \left( \frac{\boldsymbol{x}_n^\intercal\hat{\beta}}{1 - h_n} - y_n - \frac{h_n }{1 - h_n} y_n \right) \\={}& \hat{\beta}+ (\boldsymbol{X}^\intercal\boldsymbol{X})^{-1} \boldsymbol{x}_n \left( \frac{1 }{1 - h_n} \hat{y}_n - \left(1 + \frac{h_n }{1 - h_n} \right)y_n \right) \\={}& \hat{\beta}+ (\boldsymbol{X}^\intercal\boldsymbol{X})^{-1} \boldsymbol{x}_n \left( \frac{1 }{1 - h_n} \hat{y}_n - \frac{1 }{1 - h_n} y_n \right) \\={}& \hat{\beta}- (\boldsymbol{X}^\intercal\boldsymbol{X})^{-1} \boldsymbol{x}_n \frac{1 }{1 - h_n} \hat{\varepsilon}_n \end{aligned} \]
d
Letting \(\hat{y}_{n,-n} = \boldsymbol{x}_n^\intercal\hat{\boldsymbol{\beta}}_{-n}\) denote the estimate of \(y_n\) after deleting the \(n\)–th observation. Using (c), derive the following explicit expression the change in \(\hat{y}_n\) upon deleting the \(n\)–th observation:
\[ \hat{y}_{n,-n} - \hat{y}_n = \boldsymbol{x}_n^\intercal\hat{\boldsymbol{\beta}}_{-n} - \boldsymbol{x}_n^\intercal\hat{\boldsymbol{\beta}}= -\frac{h_n}{1 - h_n} \hat{\varepsilon}_n. \]
This shows that the effect of deleting observation \(n\) on \(\hat{y}_n\) is large only if both the residual and leverage score is large.
Solution Taking \(\boldsymbol{x}_n^\intercal\) times each side of the result from c gives
\[ \begin{aligned} \hat{y}_{n,-n} ={}& \boldsymbol{x}_n^\intercal\hat{\boldsymbol{\beta}}_{-n} \\={}& \boldsymbol{x}_n^\intercal\hat{\beta}- \boldsymbol{x}_n^\intercal(\boldsymbol{X}^\intercal\boldsymbol{X})^{-1} \boldsymbol{x}_n \frac{1 }{1 - h_n} \hat{\varepsilon}_n \\={}& \hat{y}_n - \frac{h_n}{1 - h_n} \hat{\varepsilon}_n. \end{aligned} \]
The result follows by subtracting \(\hat{y}_n\) from both sides.