Prediction and Model Selection

Problem 1

Given the following small dataset: \[\mathbf{y} = \begin{bmatrix}3\\5\\7\\11\\13\end{bmatrix}\] \[\mathbf{X} = \begin{bmatrix}1&1\\1&2\\1&3\\1&4\\1&5\end{bmatrix}\]

a) Fit the linear regression model \(y = \beta_0 + \beta_1 x + \varepsilon\) and calculate the fitted values \(\hat{y}_i\) and residuals \(\hat{\varepsilon}_i\) for each observation.

b) Calculate the hat matrix \(\mathbf{H} = \mathbf{X}(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\) and extract the diagonal elements \(h_{ii}\) (the leverages) for each observation. You may use R for the matrix calculations.

c) Calculate the LOOCV error using the shortcut formula: \[\text{CV} = \frac{1}{n}\sum_{i=1}^n \left(\frac{y_i - \hat{y}_i}{1 - h_{ii}}\right)^2\]

d) To verify your answer, fit the model 5 times, each time leaving out one observation. For each left-out observation \(i\):

• Fit the model on the remaining \(n-1\) observations
  
• Predict the left-out observation to get \(\hat{y}_{-i}\)
  
• Calculate \((y_i - \hat{y}_{-i})^2\)

Calculate the LOOCV error as \(\frac{1}{5}\sum_{i=1}^5 (y_i - \hat{y}_{-i})^2\). Does this match your answer from part (c)?

e) Why is the shortcut formula useful? When would computing LOOCV by actually refitting the model \(n\) times become computationally prohibitive?

Problem 2

Consider a dataset where the true relationship is quadratic: \(y = 1 + 2x - 0.5x^2 + \varepsilon\) where \(\varepsilon \sim N(0, 1)\).

Generate data with \(n = 60\) observations and \(x\) uniformly distributed between 0 and 5.

set.seed(1)
n <- 60
x <- runif(n, 0, 5)
y <- 1 + 2*x - 0.5*x^2 + rnorm(n)

a) Fit polynomial models of degrees 1 through 6:

• model1 <- lm(y ~ poly(x, 1))
  
• model2 <- lm(y ~ poly(x, 2))
  
• … and so on through degree 6

For each model, calculate:

• AIC using AIC(model)
  
• LOOCV error using the shortcut formula
  
• Adjusted \(R^2\) from summary(model)$adj.r.squared

Create a table showing all three criteria for each polynomial degree.

b) Which degree is selected by each criterion? Do they agree?

c) Create 3 plots showing how each criterion changes with polynomial degree.

d) The true model is quadratic (degree 2). Do the selection criteria correctly identify this? If any criterion selects a higher degree, explain why this might happen even when we know the true model.

e) Generate a test dataset of 100 new observations from the same true model. Calculate the true prediction error (MSE on test data) for each of the 6 models. Which degree predicts best on new data? How does this compare to what the selection criteria chose?

Problem 3

The ridge regression objective function is:

\[\min_{\beta} (\mathbf{y} - \mathbf{X}\beta)^T(\mathbf{y} - \mathbf{X}\beta) + \lambda\beta^T\beta\]

a) To find the minimum, take the derivative with respect to \(\beta\) and set it equal to zero. Solve for \(\hat{\beta}_{\text{ridge}}\).

b) Derive the variance for \(\hat{\beta}_{\text{ridge}}\).

Problem 4

Generate data with highly correlated predictors:

set.seed(1)
n <- 50
x1 <- rnorm(n)
x2 <- x1 + rnorm(n, sd = 0.01)  # Highly correlated with x1
x3 <- rnorm(n)  # Independent
X <- cbind(scale(x1), scale(x2), scale(x3))
beta_true <- c(3, -2, 1.5)
y <- X %*% beta_true + rnorm(n, sd = 1)

a) Calculate \(\hat\beta\) using the OLS and examine the coefficients. Calculate the standard errors of the coefficient estimates. Do the coefficients have high standard errors? Why might this be?

b) Calculate the correlation matrix of the predictors. Which predictors are highly correlated?

c) Calculate the Variance Inflation Factor (VIF) for each predictor using:

\[\text{VIF}_j = \frac{1}{1 - R^2_j}\]

where \(R^2_j\) is from regressing \(x_j\) on all other predictors. Which predictor(s) have high VIF values?

d) Using \(\lambda=1\) fit a ridge regression model using the same data. Calculate the standard errors of the coefficients – how do they compare to the OLS errors calculated in part (a)?

Problem 5

Consider a regression setting where we observe data \((x_i, y_i)\) for \(i = 1, \ldots, n\) and the true relationship is:

\[y_i = f(x_i) + \varepsilon_i\]

where \(\varepsilon_i \sim N(0, \sigma^2)\) are independent errors and \(f(x)\) is the true regression function.

Suppose we fit a model \(\hat{f}(x)\) using our training data. For a new test point \(x_0\), the expected prediction error is:

\[\text{EPE}(x_0) = E[(y_0 - \hat{f}(x_0))^2]\]

where the expectation is over both the training data (which affects \(\hat{f}\)) and the new observation \(y_0\).

a) Show that the expected prediction error can be decomposed as:

\[\text{EPE}(x_0) = \sigma^2 + \text{Bias}^2[\hat{f}(x_0)] + \text{Var}[\hat{f}(x_0)]\]

where:

• \(\text{Bias}[\hat{f}(x_0)] = E[\hat{f}(x_0)] - f(x_0)\)
  
• \(\text{Var}[\hat{f}(x_0)] = E[(\hat{f}(x_0) - E[\hat{f}(x_0)])^2]\)

Hint: Start by adding and subtracting \(E[\hat{f}(x_0)]\) and \(f(x_0)\) inside the squared term, then expand.

b) Generate data to illustrate this decomposition empirically:

set.seed(1)

# Generate single test point
x0 <- 0.5
sigma <- 1

# Simulate 1000 training datasets
n_sims <- 1000
n_train <- 20
true_f <- x0 + x0^2 + x0^3

predictions <- numeric(n_sims)

for(i in 1:n_sims) {
  x_train <- runif(n_train, 0, 1)
  y_train <- x_train + x_train^2 + x_train^3 + rnorm(n_train, sd = sigma)
  
  # Fit polynomial model of degree 9
  model <- lm(y_train ~ poly(x_train, 9))
  predictions[i] <- predict(model, newdata = data.frame(x_train = x0))
}

Calculate the three components at \(x_0 = 0.5\):

• Irreducible error: \(\sigma^2\)
  
• Squared bias: \((E[\hat{f}(x_0)] - f(x_0))^2\) (estimate \(E[\hat{f}(x_0)]\) as mean(predictions))
  
• Variance: \(\text{Var}[\hat{f}(x_0)]\) (estimate as var(predictions))

Verify that their sum approximately equals the mean squared error: \(\frac{1}{1000}\sum_{i=1}^{1000}(y_0^{(i)} - \hat{f}^{(i)}(x_0))^2\) where you generate new test observations.

c) Repeat part (b) for polynomial degrees 1, 2, 5, and 9 Create a plot showing how bias squared, variance, and total expected prediction error change with model complexity. Explain the bias-variance tradeoff you observe: why does bias decrease and variance increase as we add more polynomial terms? Which degree would you choose and why?

Problem 6

Generate data from a simple linear model:

set.seed(1)
n <- 40
x <- runif(n, 0, 10)
y <- 5 + 2*x + rnorm(n)
model <- lm(y ~ x)

a) For a new observation at \(x_{\text{new}} = 6\), calculate:

• The point prediction \(\hat{y}_{\text{new}}\)
  
• A 95% confidence interval for \(E[y \mid x = 6]\) (the mean response)
  
• A 95% prediction interval for a new observation \(y_{\text{new}}\) at \(x = 6\).

b) Explain the difference between these two intervals. Why is the prediction interval wider than the confidence interval?

c) The variance of the prediction error is: \[\text{Var}(y_{\text{new}}-\hat{y}_{\text{new}}) = \sigma^2\left(1 + x_{\text{new}}^T(\mathbf{X}^T\mathbf{X})^{-1}x_{\text{new}}\right)\]

Show your work to derive why this is true.