STAT151A Homework 3.

This homework is due on Gradescope on Friday October 11th at 9pm.

1 Interpretation of transforms of the response

Suppose I have data of the where $n=1,\dots 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.

2 Fit and regressors

Given a regression on $\mathit{X}$ with $P$ regressors and $N$ data points, and the corresponding $\mathit{Y}$, $\hat{\mathit{Y}}$, and $\hat{\epsilon }$, define the following quantities: $$\begin{array}{rlr}RSS:=& {\hat{\epsilon }}^{⊺}\hat{\epsilon }& \text{(Residual sum of squares)}\\ TSS:=& {\mathit{Y}}^{⊺}\mathit{Y}& \text{(Total sum of squares)}\\ ESS:=& {\hat{\mathit{Y}}}^{⊺}\hat{\mathit{Y}}& \text{(Explained sum of squares)}\\ R^2:=& \frac{ESS}{TSS}.\end{array}$$

1. Prove that $RSS+ESS=TSS$.
2. Express $R^2$ in terms of $TSS$ and $RSS$.
3. What is $R^2$ when we include no regressors? ($P=0$)
4. What is $R^2$ when we include $N$ linearly independent regressors? ($P=N$)
5. Can $R^2$ ever decrease when we add a regressor? If so, how?
6. Can $R^2$ ever stay the same when we add a regressor? If so, how?
7. Can $R^2$ ever increase when we add a regressor? If so, how?
8. Does a high $R^2$ mean the regression is useful? (You may argue by example.)
9. Does a low $R^2$ mean the regression is not useful? (You may argue by example.)

Solutions:

1. This follows from ${\hat{\mathit{Y}}}^{⊺}\hat{\epsilon }=\mathbf{0}$.
2. $R^2=\frac{TSS-RSS}{TSS}=1-\frac{RSS}{TSS}$
3. $R^2=0$ when we include no regressors
4. $R^2=1$ when we include $N$ linearly independent regressors?
5. No, it cannot, since you project onto the same or larger subspace.
6. Yes, if you add a regressor column that is colinear with the existing columns.
7. Yes, if you add a linearly independent regressor column.
8. No, you might overfit.
9. No, you might have low signal to noise ratio.

3 Prediction in the bodyfat example

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.

• Explain why the fitted error did not change.
• Explain why the coefficients did change.

The new regressors are invertible linear transforms of the original regressors. This means $\hat{\mathit{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.

• Explain why the fitted error did not change.
• Explain why the difference coefficient was not estimated by 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 ${\mathit{X}}^{⊺}\mathit{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.

• Explain why we could we estimate the coefficient of the ratio of height to weight, but not the difference.
• Explain why the fitted error changed.
• It happened that, by including the regressor 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=({\text{Abdomen}}_n,{\text{Height}}_n,{\text{Weight}}_n{)}^{⊺}$ 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=\mathit{A}{x}_n$, where the matrix $\mathit{A}$ is chosen using a machine learning algorithm.

• Will this result produce a better fit to the data than simply regressing $y_n\sim {\mathit{x}}_n$? Why or why not?

Since $\mathit{A}{x}_n$ is a linear combination of ${\mathit{x}}_n$, it can never give a better fit than regressing on $x_n$ alone. (It can give a worse fit if $\mathit{A}$ is lower rank than ${\mathit{x}}_n$.) Searching over all possible $\mathit{A}$ will never give you a better fit than $y_n\sim {\mathit{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)$.

• Will this result produce a better fit to the data than simply regressing $y_n\sim {\mathit{x}}_n$? Why or why not?
• Do you think this result produce a useful prediction for new data? Why or why not?

Since $f({\mathit{x}}_n)$ can be a nonlinear function, you can produce a better fit than regressing on ${\mathit{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.

4 Leaving a single datapoint out of regression

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 $\mathit{u}$ and $\mathit{v}$ vectors the same length as $A$. Then

$$(A+\mathit{u}{\mathit{v}}^{⊺}{)}^{-1}=A^{-1}-\frac{A^{-1}\mathit{u}{\mathit{v}}^{⊺}{A}^{-1}}{1+{\mathit{v}}^{⊺}{A}^{-1}\mathit{u}}.$$

We will also use the definition of a “leverage score” from lecture $h_n:={\mathit{x}}_n^{⊺}({\mathit{X}}^{⊺}\mathit{X}{)}^{-1}{\mathit{x}}_n$. Note that $h_n=(\mathit{X}({\mathit{X}}^{⊺}\mathit{X}{)}^{-1}{\mathit{X}}^{⊺}{)}_{nn}$ is the $n$–th diagonal entry of the projection matrix $\underset{\mathit{X}}{\mathit{P}}$.

Let ${\hat{\mathit{\beta }}}_{-n}$ denote the estimate of $\hat{\beta }$ with the datapoint $n$ left out. Similarly, let ${\mathit{X}}_{-n}$ denote the $\mathit{X}$ matrix with row $n$ left out, and ${\mathit{Y}}_{-n}$ denote the $\mathit{Y}$ matrix with row $n$ left out.

a Prove that

$${\hat{\mathit{\beta }}}_{-n}=({\mathit{X}}_{-n}^{⊺}{\mathit{X}}_{-n}{)}^{-1}{\mathit{X}}_{-n}^{⊺}{\mathit{Y}}_{-n}=({\mathit{X}}^{⊺}\mathit{X}-{\mathit{x}}_n{\mathit{x}}_n^{⊺}{)}^{-1}({\mathit{X}}^{⊺}\mathit{Y}-{\mathit{x}}_n{y}_n)$$

Solution

This follows from ${\mathit{X}}^{⊺}\mathit{X}=\sum _{n=1}^N{\mathit{x}}_n{\mathit{x}}_n^{⊺}$ and ${\mathit{X}}^{⊺}\mathit{Y}=\sum _{n=1}^N{\mathit{x}}_n{y}_n$.

b

Using the Woodbury formula, derive the following expression: $$({\mathit{X}}^{⊺}\mathit{X}-{\mathit{x}}_n{\mathit{x}}_n^{⊺}{)}^{-1}=({\mathit{X}}^{⊺}\mathit{X}{)}^{-1}+\frac{({\mathit{X}}^{⊺}\mathit{X}{)}^{-1}{\mathit{x}}_n{\mathit{x}}_n^{⊺}({\mathit{X}}^{⊺}\mathit{X}{)}^{-1}}{1-h_n}$$

Solution

Direct application of the formula with $\mathit{u}={\mathit{x}}_n$ and $\mathit{v}=-{\mathit{x}}_n$ gives $$({\mathit{X}}^{⊺}\mathit{X}-{\mathit{x}}_n{\mathit{x}}_n^{⊺}{)}^{-1}=({\mathit{X}}^{⊺}\mathit{X}{)}^{-1}+\frac{({\mathit{X}}^{⊺}\mathit{X}{)}^{-1}{\mathit{x}}_n{\mathit{x}}_n^{⊺}({\mathit{X}}^{⊺}\mathit{X}{)}^{-1}}{1-{\mathit{x}}_n^{⊺}({\mathit{X}}^{⊺}\mathit{X}{)}^{-1}{\mathit{x}}_n}.$$

Then recognize the leverage score.

c

Combine (a) and (b) to derive the following explicit expression for ${\hat{\mathit{\beta }}}_{-n}$:

$${\hat{\mathit{\beta }}}_{-n}=\hat{\mathit{\beta }}-({\mathit{X}}^{⊺}\mathit{X}{)}^{-1}{\mathit{x}}_n\frac{1}{1-h_n}{\hat{\epsilon }}_n.$$

Solution

We have

$$({\mathit{X}}^{⊺}\mathit{X}-{\mathit{x}}_n{\mathit{x}}_n^{⊺}{)}^{-1}{\mathit{X}}^{⊺}\mathit{Y}=\hat{\beta }+\frac{({\mathit{X}}^{⊺}\mathit{X}{)}^{-1}{\mathit{x}}_n{\mathit{x}}_n^{⊺}\hat{\beta }}{1-h_n}.$$

and

$$({\mathit{X}}^{⊺}\mathit{X}-{\mathit{x}}_n{\mathit{x}}_n^{⊺}{)}^{-1}{\mathit{x}}_n{y}_n=({\mathit{X}}^{⊺}\mathit{X}{)}^{-1}{\mathit{x}}_n{y}_n+\frac{({\mathit{X}}^{⊺}\mathit{X}{)}^{-1}{\mathit{x}}_n{h}_n}{1-h_n}{y}_n.$$

Combining,

$$\begin{array}{rl}{\hat{\mathit{\beta }}}_{-n}=& \hat{\beta }+({\mathit{X}}^{⊺}\mathit{X}{)}^{-1}{\mathit{x}}_n(\frac{{\mathit{x}}_n^{⊺}\hat{\beta }}{1-h_n}-y_n-\frac{h_n}{1-h_n}{y}_n)\\ =& \hat{\beta }+({\mathit{X}}^{⊺}\mathit{X}{)}^{-1}{\mathit{x}}_n(\frac{1}{1-h_n}{\hat{y}}_n-(1+\frac{h_n}{1-h_n})y_n)\\ =& \hat{\beta }+({\mathit{X}}^{⊺}\mathit{X}{)}^{-1}{\mathit{x}}_n(\frac{1}{1-h_n}{\hat{y}}_n-\frac{1}{1-h_n}{y}_n)\\ =& \hat{\beta }-({\mathit{X}}^{⊺}\mathit{X}{)}^{-1}{\mathit{x}}_n\frac{1}{1-h_n}{\hat{\epsilon }}_n\end{array}$$

d

Letting ${\hat{y}}_{n,-n}={\mathit{x}}_n^{⊺}{\hat{\mathit{\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={\mathit{x}}_n^{⊺}{\hat{\mathit{\beta }}}_{-n}-{\mathit{x}}_n^{⊺}\hat{\mathit{\beta }}=-\frac{h_n}{1-h_n}{\hat{\epsilon }}_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 ${\mathit{x}}_n^{⊺}$ times each side of the result from c gives

$$\begin{array}{rl}{\hat{y}}_{n,-n}=& {\mathit{x}}_n^{⊺}{\hat{\mathit{\beta }}}_{-n}\\ =& {\mathit{x}}_n^{⊺}\hat{\beta }-{\mathit{x}}_n^{⊺}({\mathit{X}}^{⊺}\mathit{X}{)}^{-1}{\mathit{x}}_n\frac{1}{1-h_n}{\hat{\epsilon }}_n\\ =& {\hat{y}}_n-\frac{h_n}{1-h_n}{\hat{\epsilon }}_n.\end{array}$$

The result follows by subtracting ${\hat{y}}_n$ from both sides.