Question 1: Another Simulation Study

Here is a GPT-5 prompt I used to create some content:

How does ridge regression deal with correlated variables? Is it always better than Marrow’s C_p for prediction? Can you suggest a code with 20 variables to confirm your answer?

The following is the answer produced by GPT-5

Here’s a concise answer:

Ridge with correlated variables: Ridge regression shrinks coefficients toward zero in a way that distributes the effect across correlated predictors. Instead of picking one variable and ignoring the others, it tends to keep correlated variables together and stabilize their estimates. This reduces variance and avoids the instability of ordinary least squares.

Ridge vs. Mallows’ \(C_p\): Ridge is not always better. Ridge generally performs better when the signal is spread across many correlated predictors (dense signals). Mallows’ \(C_p\), used for subset selection, can perform better when only a few predictors truly matter (sparse signals).

In practice, which is better depends on the structure of the true underlying model. Would you like me to provide a small R code example (with 20 variables) that demonstrates both cases — dense vs. sparse signals — so you can see the difference in prediction performance?

Part a) [15 pts]

I would like to first consider a simulation study to understand the performance. And while doing this, let’s also compare the bias and variance of the estimations. Consider the following data generator, with \(p = 20\), but \(y\) depends on only the first two variables:

\[ y = x_1 + x_2 + \epsilon \]

where \(\epsilon \sim N(0,1)\) and \(x_1\) and \(x_2\) are generated from independent standard normal distributions, while the rest of the \(p-2\) variables are generated in the following way:

\[ x_j \sim \eta_j x_1 + \phi_j x_2 + \xi_j, \quad j = 3, \ldots, p \]

where \(\eta_j, \phi_j \sim \text{Uniform}(-1, 1)\) and \(\xi_j \sim N(0, 0.3^2)\) are all independently generated. Write your code to generate this data with \(n = 100\), and fit both the ridge regression (use leave-one-out CV and lambda.min to pick the best tuning parameter) and the best subset selection (using Mallows’ \(C_p\)). Generate another \(1000\) testing data, and evaluate the prediction performance of both methods.

  • Which variable(s) were selected by your Mallows’ \(C_p\) method?
  • Which variable was given the largest coefficient (based on absolute values) by your ridge regression?
  • Which method had the smaller test error?

Part b) [25 pts]

Now, let’s consider a simulation. For this question, we will repeat the above procedure 200 times, record their prediction errors and also evaluate the bias and variance of each of the estimations of the coefficients. To do this, we need to record the estimates of the coefficients for all variables for both methods, and store them. For the Marlow’s \(C_p\) method, if a variable is not selected, then we record its coefficient as 0. After all simulations, you should have a matrix of \(200 \times 20\) for each method. The intercept is not needed. Based on this matrix, you can calculate the bias and variance of the estimations of these coefficients for each variable by

  • Comparing the mean of these estimates to the true \(\beta\) values
  • Calculating the variance of these estimates over the 200 simulation runs

Answer the following question:

  • Provide a barplot() to visualize and compare the parameter estimates for both methods.
  • Looking at the bias result overall, what patterns do you see?
  • Looking at the variance result overall, what patterns do you see?
  • Which method has a smaller testing error?

Part c) [10 pts]

Based on your results (b)

  • Can you now compare your results to the answer provided by GPT-5? Do you agree or disagree with the statements made there? Please provide your reasoning.
  • Based on these simulation results, can you comment on the bias-variance trade-off for these two methods?

Question 2: Modeling High-Dimensional Data

We will use the golub dataset from the multtest package. This dataset contains 3051 genes from 38 tumor mRNA samples from the leukemia microarray study Golub et al. (1999). This package is not included in R, but on bioconductor. Install the latest version of this package from bioconductor, and read the documentation of this dataset to understand the data structure of golub and golub.cl.

Part a) [25 points]

We will not use this data for classification (the original problem). Instead, we will do a toy regression example to show how genes are highly correlated and could be used to predict each. Carry out the following tasks:

  • Perform marginal association test for each gene with the response golub.cl using mt.teststat(). Use t.equalvar (two sample \(t\) test with equal variance) as the test statistic.
  • Sort the genes by their p-values and select the top 100 genes
  • Construct a dataset with the top 10 genes and another one (call it \(X\)) with the remaining genes
  • Perform principal component analysis (PCA) on the top 100 genes and extract the first principal component, use this as the outcome \(y\). Be careful about the orientation of the data matrix.
  • Perform ridge regression with 19-fold cross-validation on \(X\) and the outcome \(y\). Does your model fit well? Can you provide detailed model fitting results to support your claim?
  • Fit ridge regression but use GCV as the criterion. Does your model fit well?

Part b) [5 points]

Based on your results, do you observe a “U” shaped curve? Can you comment on how is this related to the bias-variance trade-off?

Question 3: Linear Regression with Coordinate Descent

Recall the previous homework, we have a quadratic function for minimization. We know that analytic solution exist. However, in this example, let’s use coordinate descent to solve the problem. To demonstrate this, let’s consider the following simulated dataset, with design matrix \(x\) (without intercept) and response vector \(y\):

  set.seed(432)
  n <- 150
  x <- matrix(rnorm(n*2), n, 2)
  y <- 0.3 * x[, 1] + 0.5 * x[, 2] + rnorm(n)

We will consider a model without the intercept term. In this case, our objective function (of \(\beta_1\) and \(\beta_2\) for linear regression is to minimize the sum of squared residuals:

\[ f(\beta_1, \beta_2) = \frac{1}{n} \sum_{i=1}^n (y_i - \beta_1 x_{i1} - \beta_2 x_{i2})^2 \]

where \(x_{ij}\) represents the \(j\)th variable of the \(i\)th observation.

Part a) [10 points]

Write down the objective function in the form of

\[ f(x,y) = a \beta_1^2 + b \beta_2^2 + c \beta_1 \beta_2 + d \beta_1 + e \beta_2 + f \]

by specifying what are coefficients a, b, c, d, e, and f, using the simulated data. Calculate them in R, using vector operations rather than for-loops.

Part b) [10 points]

A coordinate descent algorithm essentially does two steps:

  1. Update \(\beta_1\) to its optimal value while keeping \(\beta_2\) fixed
  2. Update \(\beta_2\) to its optimal value while keeping \(\beta_1\) fixed

Write down the updating rules for \(\beta_1\) and \(\beta_2\) using the coordinate descent algorithm. Use those previously defined coefficients in your formula and write them in Latex. Implement them in a for-loop algorithm in R that iterates at most 100 times. Use the initial values \(\beta_1 = 0\) and \(\beta_2 = 0\). Decide your stopping criterion based on the change in \(\beta_1\) and \(\beta_2\). Validate your solution using the lm() function.