---
title: "Homework 04"
pagetitle: "Homework 04"
body-classes: "lecture-page practice-page"
format:
  html:
    embed-resources: true
    html-math-method: mathml
    page-layout: full
    toc: true
    toc-location: body
    toc-title: "On this page"
    toc-depth: 2
    code-fold: true
    code-summary: "Show the workflow"
---

## Question 1: Signal strength and correlated substitutes

Consider a regression problem with $n=120$ observations and $p=10$ covariates. Create one simulated dataset as follows. Generate the first covariate $X_1$ as $n$ independent standard normal draws. Then form the second covariate as

$$
X_2 = 0.999X_1 + \sqrt{1-0.999^2}\,Z,
$$

where $Z$ is another vector of $n$ independent standard normal draws, so that $X_1$ and $X_2$ have population correlation $0.999$. Generate the remaining covariates $X_3,\ldots,X_{10}$ as independent standard normal vectors of length $n$.

Thus $X_1$ and $X_2$ are extremely highly correlated, while $X_3$ and $X_4$ are independent of one another and of the remaining covariates. Finally, generate the response as

$$
Y
=
0.15X_1
+0.15X_2
+0.30X_3
+0.15X_4
+\epsilon,
$$

where $\epsilon$ is a vector of $n$ independent standard normal errors, independent of the covariates. Both members of the correlated pair appear in the generating equation with the same coefficient, even though each one carries almost the same observed information as the other. The two independent signals $X_3$ and $X_4$ have coefficients that differ by a factor of two, and $X_5,\ldots,X_{10}$ have zero coefficients.

Use seed `43246` and independently repeat the complete data generation and fitting procedure 200 times. Within each repetition, center and standardize the covariates using divisor $n$ and center the response. Fit the lasso model with the `glmnet` package in R or `sklearn.linear_model.Lasso` in Python, without scaling the variables inside the fitting function, at

$$
\lambda\in\{0.08,0.12,0.18,0.24\}.
$$

Make sure you understand which argument specifies the scaling option in a lasso fit and how to give a specific penalty value; the default settings of the fitting function do not produce the fits required here.

Fit all four penalties to the same generated dataset before beginning the next repetition. Leave the intercept unpenalized and convert the fitted slopes back to the original covariate scale. Store the slopes in one array whose dimensions correspond to the 200 repetitions, 10 covariates, and four penalties. Treat a slope as nonzero when its absolute value exceeds $10^{-8}$. Use a tight convergence tolerance as well as a sufficiently large iteration limit, because the nearly identical covariates make coefficient allocation sensitive to numerical accuracy. In R, use `control = list(thresh = 1e-18, maxit = 1000000)`; in Python, use `tol=1e-10, max_iter=100000`. Verify that all 800 fits converge. R and Python use different random-number generators, so their exact numerical results need not agree.

a. For every covariate and penalty, report the selection frequency and the mean absolute fitted slope, where the mean includes fitted zeros. Draw one figure for each summary. Explain how increasing the penalty changes the fitted magnitudes and selection frequencies.

b. For the pair $(X_1,X_2)$, report the frequency of four mutually exclusive outcomes at each penalty: only the first covariate is selected, only the second is selected, both are selected, and neither is selected. Compare the two marginal selection frequencies and explain why they are similar. Describe how the four outcomes shift as the penalty increases, and explain why lasso can produce a sparse fitted rule that omits one member of the pair even though both members appear in the generating equation.

c. Compare the selection frequencies of $X_3$ and $X_4$ with the average selection frequency among $X_5,\ldots,X_{10}$. Explain how signal magnitude and the penalty account for the pattern. Use the correlated pair to explain why a fitted zero does not establish that the corresponding population coefficient is zero, and why a nonzero fitted coefficient does not establish that the variable is uniquely important.

## Question 2: Elastic net with an equal penalty mix

This question continues the simulation of Question 1. Generate the data exactly as in Question 1, with the same seed `43246`, so that the 200 repetitions produce the same datasets and the only change is the fitting method. Now fit the elastic net model with mixing parameter $\alpha=0.5$, which places equal weight on the $\ell_1$ and squared $\ell_2$ penalties. Use the `glmnet` package in R or `sklearn.linear_model.ElasticNet` in Python, with the same four penalties $\lambda\in\{0.08,0.12,0.18,0.24\}$ and the same fitting details as in Question 1. In R, additionally set `family = gaussian()` to avoid the response-rescaling convention of the default Gaussian solver and match the stated penalty mix. Make sure you understand which argument controls this mix in the package you use.

a. For every covariate and penalty, report the selection frequency, and draw the selection-frequency figure. Compare the frequencies with the lasso results from Question 1: which covariates change the most, and in which direction?

b. For the pair $(X_1,X_2)$, report the frequency of four mutually exclusive outcomes at each penalty: only the first covariate is selected, only the second is selected, both are selected, and neither is selected. Compare these frequencies with Question 1. Are the two covariates now properly selected together? Explain why the squared $\ell_2$ part of the penalty encourages the fit to keep both members of the pair.

c. Compare the average selection frequency among the noise covariates $X_5,\ldots,X_{10}$ with the lasso results at each penalty. Are the noise covariates still screened out as the penalty increases? At the same numerical value of $\lambda$, the elastic net applies a weaker $\ell_1$ threshold than the lasso; use this to explain any difference. State the penalties at which the elastic net selects both members of the pair while screening out most noise covariates.

## Question 3: Coordinate descent and the lasso path

This question asks you to implement coordinate descent yourself and use it to compute a complete lasso path. Create one simulated dataset with $n=100$ observations and $p=3$ covariates. Generate three covariates as independent standard normal vectors of length $n$, and generate the response as

$$
Y
=
1.5+X_1-2X_2+\epsilon,
$$

where $\epsilon$ is a vector of $n$ independent standard normal errors, independent of the covariates. The third covariate does not appear in the generating equation. Use seed `43247`. Center the response and center and standardize each covariate using divisor $n$. The intercept is fitted separately as the mean of the response and is not penalized.

The lasso objective is

$$
L_\lambda(\boldsymbol\beta)
=
\frac{1}{2n}
\left\lVert
\widetilde{\mathbf y}-\mathbf X\boldsymbol\beta
\right\rVert_2^2
+
\lambda\lVert\boldsymbol\beta\rVert_1,
$$

where $\widetilde{\mathbf y}$ is the centered response and $\mathbf X$ is the standardized covariate matrix. Coordinate descent updates one slope at a time while holding the other slopes fixed. For coordinate $j$, form the residual that leaves covariate $j$ out of the current fit, and compute the score $a_j$ as the mean of the elementwise product of covariate $j$ and that residual. Because each standardized column satisfies $\frac{1}{n}\mathbf x_j^{\mathsf T}\mathbf x_j=1$, the update is the soft-thresholding rule from the Week 4 lecture:

$$
\widehat\beta_j\leftarrow S(a_j,\lambda),
\qquad
S(a,\lambda)=\operatorname{sign}(a)(|a|-\lambda)_+.
$$

Follow this procedure:

1. Fix a value of $\lambda$ and loop over the indices $j=1,\ldots,p$ until the coefficients stop changing. For each $j$, form the residual without covariate $j$, compute its score, and save the soft-thresholded value as the new $\widehat\beta_j$ before moving to the next index.
2. Check convergence after each complete cycle: stop when the largest absolute coefficient change is below a small tolerance, such as $10^{-8}$. Record the fitted vector for that $\lambda$.
3. Then reduce $\lambda$ to the next grid value and start the next fit from the current coefficients (a warm start). Repeat for the whole grid

$$
\lambda\in\{2.5,2.0,1.5,1.0,0.6,0.4,0.25,0.15,0.08,0.03\}.
$$

a. Implement the procedure. Report the fitted coefficient vector at each grid value in a small table, and verify that every fit converged. Compute $\lambda_{\max}=\lVert\frac{1}{n}\mathbf X^{\mathsf T}\widetilde{\mathbf y}\rVert_\infty$ and check that all slopes are zero for grid values at or above it.

b. Plot the coefficient path: one curve per covariate, with the fitted slope on the vertical axis and $\log(\lambda)$ on the horizontal axis. Describe the order in which the covariates enter the model and explain that order using the sizes of the true coefficients.

c. The third covariate has a zero population coefficient. Does it enter the path at small penalties in your run? Compare the fit at the smallest grid value with a least-squares fit on the same standardized data, and explain why the two are close but not identical.

## Question 4: Tuning and comparing penalized regressions for diabetes prediction

The supplied `data/diabetes.csv` contains 442 observations of ten baseline covariates and the quantitative response `y`. The file `data/diabetes-split-folds.csv` assigns 353 observations to the training set and 89 observations to the final test set; it also supplies ten fold labels for the training observations.

Tune three penalized regressions (lasso, ridge, and elastic net with mixing parameter $\alpha=0.5$) for predicting `y` from the ten covariates. Using only the training observations, perform ten-fold cross-validation with the supplied folds for each method and select each penalty by the smallest mean cross-validation error. In R, use `cv.glmnet` with its default Gaussian solver and package-generated penalty grids. In Python, use `GridSearchCV` with 61 logarithmically spaced penalty values from $10^{-4}$ to $10^2$ for each method, placing `StandardScaler` inside a `Pipeline`. Make sure you understand which argument selects the penalty family. In every fit, estimate the covariate means and scales from the observations used for that fit, center the response, and leave the intercept unpenalized. The package-specific grids and penalty conventions can give different tuning results across languages.

Select the method with the smallest minimum cross-validation error. Produce the cross-validation plot for each of the three fits (in R, the default plot of a `cv.glmnet` object) and a table comparing the three selected penalties and their cross-validation errors. Then refit the chosen model on all training observations, evaluate the final test set once, and report the selected method, its penalty, and the test error.