Homework 11
Robert Papshev
Part 1: Bears
bears <- read.table("/Users/robertpapshev/Downloads/bears.txt", header = TRUE) # Read the bears data from Downloads.
bears_indep <- bears[bears$Obs.No == 1, ] # Keep the first observation for each bear.
bears_indep$d <- as.integer(bears_indep$Sex == 2) # Code female as 1 and male as 0.
bears_indep$y <- sqrt(bears_indep$Weight) # Define the response.
predictors <- c("Head.L", "Head.W", "Neck.G", "Length", "Chest.G", "d") # Specify the six predictors.
full_bears <- lm(reformulate(predictors, response = "y"), data = bears_indep) # Fit the full model without interactions.
nrow(bears_indep) # Report the number of bears.## [1] 99I used the 99 first observations. There were no missing values in the response or six predictors.
1. Forward selection
forward <- lm(y ~ 1, data = bears_indep) # Start with the intercept-only model.
forward_steps <- data.frame() # Store the accepted additions.
repeat { # Continue until no remaining predictor meets the entry threshold.
candidates <- add1(forward, scope = formula(full_bears), test = "F") # Test each remaining addition.
candidates <- candidates[rownames(candidates) != "<none>", , drop = FALSE] # Keep only additions.
if (nrow(candidates) == 0) break # Stop if all predictors have entered.
j <- which.min(candidates[["Pr(>F)"]]) # Find the smallest partial F-test p-value.
if (candidates[["Pr(>F)"]][j] >= 0.10) break # Stop if the entry threshold is not met.
added <- rownames(candidates)[j] # Record the selected predictor.
forward_steps <- rbind(forward_steps, data.frame(Step = nrow(forward_steps) + 1, Added = added, F = candidates[["F value"]][j], p_value = candidates[["Pr(>F)"]][j])) # Save the step's test results.
forward <- update(forward, paste(". ~ . +", added)) # Add the selected predictor.
} # End forward selection.
forward_display <- forward_steps # Make a display copy of the step results.
forward_display$p_value <- format.pval(forward_display$p_value, digits = 4, eps = 1e-99) # Preserve small p-values.
knitr::kable(forward_display, digits = 4) # Print the selection sequence.| Step | Added | F | p_value |
|---|---|---|---|
| 1 | Chest.G | 1970.9364 | 2.968e-66 |
| 2 | Neck.G | 37.8245 | 1.773e-08 |
| 3 | Length | 29.0018 | 5.222e-07 |
summary(forward) # Print the final model estimates and tests.##
## Call:
## lm(formula = y ~ Chest.G + Neck.G + Length, data = bears_indep)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.4439 -0.3976 -0.0565 0.3903 1.8841
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.3737 0.4447 -12.08 < 2e-16 ***
## Chest.G 0.2550 0.0248 10.30 < 2e-16 ***
## Neck.G 0.2122 0.0384 5.53 2.8e-07 ***
## Length 0.0800 0.0149 5.39 5.2e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.663 on 95 degrees of freedom
## Multiple R-squared: 0.974, Adjusted R-squared: 0.973
## F-statistic: 1.2e+03 on 3 and 95 DF, p-value: <2e-16add1(forward, scope = formula(full_bears), test = "F") # Show the remaining entry tests.## Single term additions
##
## Model:
## y ~ Chest.G + Neck.G + Length
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 41.82 -77.31
## Head.L 1 0.8615 40.96 -77.38 1.977 0.163
## Head.W 1 0.0683 41.75 -75.48 0.154 0.696
## d 1 0.0192 41.80 -75.36 0.043 0.836I added Chest.G first (\(F=1970.9364\), \(p=2.9684\times10^{-66}\)), then
Neck.G (\(F=37.8245\),
\(p=1.7732\times10^{-8}\)), then
Length (\(F=29.0018\),
\(p=5.2224\times10^{-7}\)). The
smallest remaining entry p-value was 0.1630 for Head.L, so
I stopped.
The final fitted model is
\[ \widehat{\sqrt{\mathrm{Weight}}} = -5.37368 + 0.25499\,\mathrm{Chest.G} + 0.21225\,\mathrm{Neck.G} + 0.08005\,\mathrm{Length}. \]
backward <- full_bears # Start backward selection with all six predictors.
backward_steps <- data.frame() # Store the removed predictors.
repeat { # Continue until all remaining predictors meet the retention threshold.
candidates <- drop1(backward, test = "F") # Test each possible deletion.
candidates <- candidates[rownames(candidates) != "<none>", , drop = FALSE] # Keep only deletions.
if (nrow(candidates) == 0) break # Stop if no predictors remain.
j <- which.max(candidates[["Pr(>F)"]]) # Find the largest partial F-test p-value.
if (candidates[["Pr(>F)"]][j] <= 0.10) break # Retain predictors significant at 10%.
removed <- rownames(candidates)[j] # Record the predictor to remove.
backward_steps <- rbind(backward_steps, data.frame(Step = nrow(backward_steps) + 1, Removed = removed, p_value = candidates[["Pr(>F)"]][j])) # Save the deletion.
backward <- update(backward, paste(". ~ . -", removed)) # Remove the selected predictor.
} # End backward selection.
knitr::kable(backward_steps, digits = 4) # Print the deletion sequence.| Step | Removed | p_value |
|---|---|---|
| 1 | d | 0.8342 |
| 2 | Head.W | 0.4955 |
| 3 | Head.L | 0.1630 |
formula(backward) # Print the backward-selected model.## y ~ Neck.G + Length + Chest.GUsing a 10% removal threshold for the backward comparison, I removed
d, Head.W, and Head.L, in that
order. Forward and backward selection gave the same final model.
2. Best subsets by RSS
best_subsets <- regsubsets(y ~ Head.L + Head.W + Neck.G + Length + Chest.G + d, data = bears_indep, nvmax = 6, nbest = 1, method = "exhaustive") # Find the lowest-RSS model at each size.
subset_summary <- summary(best_subsets) # Extract the subset results.
best_terms <- lapply(1:6, function(k) names(coef(best_subsets, id = k))[-1]) # Extract each model's predictors.
best_models <- lapply(best_terms, function(v) lm(reformulate(v, response = "y"), data = bears_indep)) # Refit each best subset.
subset_table <- data.frame(Predictors = 1:6, Model = vapply(best_terms, paste, character(1), collapse = " + "), RSS = subset_summary$rss) # Summarize the best models.
knitr::kable(subset_table, digits = 4) # Display models and RSS values.| Predictors | Model | RSS |
|---|---|---|
| 1 | Chest.G | 76.0932 |
| 2 | Neck.G + Chest.G | 54.5860 |
| 3 | Neck.G + Length + Chest.G | 41.8193 |
| 4 | Head.L + Neck.G + Length + Chest.G | 40.9578 |
| 5 | Head.L + Head.W + Neck.G + Length + Chest.G | 40.7526 |
| 6 | Head.L + Head.W + Neck.G + Length + Chest.G + d | 40.7331 |
The forward-selected model is also the lowest-RSS three-predictor model, with RSS = 41.8193.
3-6. Adjusted R-squared, Mallows’ Cp, AIC, and BIC
I counted the intercept in the number of betas, so \(p=k+1\) for a model with \(k\) predictors. I used the lab formulas \(\mathrm{AIC}=n\log(\mathrm{RSS}/n)+2p\) and \(\mathrm{BIC}=n\log(\mathrm{RSS}/n)+p\log(n)\).
n <- nrow(bears_indep) # Count the observations used in all bear models.
p <- 2:7 # Count regression coefficients, including the intercept.
full_mse <- deviance(full_bears) / df.residual(full_bears) # Estimate error variance from the full model.
cp <- subset_summary$rss / full_mse - n + 2 * p # Compute Mallows' Cp.
aic <- n * log(subset_summary$rss / n) + 2 * p # Compute AIC using the lab formula.
bic <- n * log(subset_summary$rss / n) + log(n) * p # Compute BIC using the lab's stated log(n) penalty.
criteria <- data.frame(Predictors = 1:6, Betas = p, Adjusted_R2 = vapply(best_models, function(m) summary(m)$adj.r.squared, numeric(1)), Cp = cp, AIC = aic, BIC = bic) # Combine the selection criteria.
knitr::kable(criteria, digits = 6) # Display the criteria for all six models.| Predictors | Betas | Adjusted_R2 | Cp | AIC | BIC |
|---|---|---|---|---|---|
| 1 | 2 | 0.952610 | 76.86454 | -22.0529 | -16.8627 |
| 2 | 3 | 0.965650 | 30.28827 | -52.9389 | -45.1535 |
| 3 | 4 | 0.973407 | 3.45338 | -77.3143 | -66.9338 |
| 4 | 5 | 0.973678 | 3.50749 | -77.3752 | -64.3996 |
| 5 | 6 | 0.973528 | 5.04408 | -75.8724 | -60.3017 |
| 6 | 7 | 0.973253 | 7.00000 | -73.9198 | -55.7539 |
plot(p, criteria$Adjusted_R2, type = "b", pch = 19, xaxt = "n", xlab = "Number of betas (including intercept)", ylab = "Adjusted R-squared") # Plot adjusted R-squared against model size.
axis(1, at = p) # Label each model size.
Adjusted \(R^2\) is largest for
Head.L + Neck.G + Length + Chest.G (\(p=5\), adjusted \(R^2=0.973678\)).
plot(p, criteria$Cp, type = "b", pch = 19, xaxt = "n", xlab = "Number of betas (including intercept)", ylab = "Mallows' Cp") # Plot Cp against model size.
axis(1, at = p) # Label each model size.
abline(a = 0, b = 1, col = "red", lty = 2) # Add the reference line y = p.
legend("topright", legend = "Cp = p", col = "red", lty = 2, bty = "n") # Label the reference line.
The minimum \(C_p\) is 3.453379 for
Neck.G + Length + Chest.G (\(p=4\)), which is also close to the \(C_p=p\) line. I selected this model by
\(C_p\).
plot(p, criteria$AIC, type = "b", pch = 19, xaxt = "n", xlab = "Number of betas (including intercept)", ylab = "AIC") # Plot AIC against model size.
axis(1, at = p) # Label each model size.
AIC is smallest for Head.L + Neck.G + Length + Chest.G
(\(p=5\), AIC = -77.375182). The
three-predictor model is nearly tied (AIC = -77.314327).
plot(p, criteria$BIC, type = "b", pch = 19, xaxt = "n", xlab = "Number of betas (including intercept)", ylab = "BIC") # Plot BIC against model size.
axis(1, at = p) # Label each model size.
BIC is smallest for Neck.G + Length + Chest.G (\(p=4\), BIC = -66.933848).
Part 2: Prostate
1. Full model
data("prostate", package = "faraway") # Load all 97 observations.
full_prostate <- lm(lpsa ~ ., data = prostate) # Fit lpsa on all eight other variables.
summary(full_prostate) # Print coefficients, t-tests, and model fit statistics.##
## Call:
## lm(formula = lpsa ~ ., data = prostate)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.733 -0.371 -0.017 0.414 1.638
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.66934 1.29639 0.52 0.6069
## lcavol 0.58702 0.08792 6.68 2.1e-09 ***
## lweight 0.45447 0.17001 2.67 0.0090 **
## age -0.01964 0.01117 -1.76 0.0823 .
## lbph 0.10705 0.05845 1.83 0.0704 .
## svi 0.76616 0.24431 3.14 0.0023 **
## lcp -0.10547 0.09101 -1.16 0.2496
## gleason 0.04514 0.15746 0.29 0.7750
## pgg45 0.00453 0.00442 1.02 0.3089
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.708 on 88 degrees of freedom
## Multiple R-squared: 0.655, Adjusted R-squared: 0.623
## F-statistic: 20.9 on 8 and 88 DF, p-value: <2e-16The full model has \(R^2=0.6548\),
adjusted \(R^2=0.6234\), and residual
standard error 0.7084. Its overall test is \(F(8,88)=20.86\), \(p<2.2\times10^{-16}\). At the 5% level,
lcavol, lweight, and svi are
significant (\(p=2.11\times10^{-9}\),
0.00896, and 0.00233, respectively).
2. Ridge regression
ridge_grid <- seq(0, 10, length.out = 100) # Define 100 equally spaced penalties.
ridge_fit <- lm.ridge(lpsa ~ ., data = prostate, lambda = ridge_grid) # Fit ridge across the requested grid.
ridge_index <- which.min(ridge_fit$GCV) # Find the smallest GCV score.
ridge_lambda <- ridge_fit$lambda[ridge_index] # Extract the selected penalty.
ridge_coef <- coef(ridge_fit)[ridge_index, ] # Extract coefficients on the original predictor scales.
names(ridge_coef)[1] <- "(Intercept)" # Match the OLS intercept name.
knitr::kable(data.frame(Lambda = ridge_lambda, GCV = ridge_fit$GCV[ridge_index]), digits = 9) # Report the selected penalty and score.| Lambda | GCV | |
|---|---|---|
| 6.464646 | 6.46465 | 0.00551973 |
knitr::kable(cbind(OLS = coef(full_prostate), Ridge = ridge_coef), digits = 6) # Compare ridge and OLS estimates.| OLS | Ridge | |
|---|---|---|
| (Intercept) | 0.669337 | 0.489814 |
| lcavol | 0.587022 | 0.517343 |
| lweight | 0.454467 | 0.443616 |
| age | -0.019637 | -0.015545 |
| lbph | 0.107054 | 0.096111 |
| svi | 0.766157 | 0.696395 |
| lcp | -0.105474 | -0.042807 |
| gleason | 0.045142 | 0.060521 |
| pgg45 | 0.004525 | 0.003520 |
The minimum GCV is 0.005519732 at \(\lambda=6.464646\). Ridge retains all eight
predictors. Most slope magnitudes decrease; for example,
lcavol decreases from 0.587022 to 0.517343 and
lcp changes from -0.105474 to -0.042807. The
gleason coefficient increases from 0.045142 to
0.060521.
3. LASSO coefficient paths
x <- model.matrix(lpsa ~ ., data = prostate)[, -1] # Build the predictor matrix without an intercept column.
y <- prostate$lpsa # Extract the response.
lasso_fit <- glmnet(x, y, family = "gaussian", alpha = 1, standardize = TRUE) # Fit the LASSO coefficient path.
plot(lasso_fit, xvar = "lambda", label = TRUE) # Plot the coefficient paths as a function of the penalty.
As the penalty increases, the coefficients generally shrink toward
zero and predictors drop out. At sufficiently large penalties, every
slope is zero; lcavol remains nonzero at the largest
penalties before that point.
4. Cross-validated LASSO estimates
set.seed(462) # Fix the random fold assignment.
cv_lasso <- cv.glmnet(x, y, family = "gaussian", alpha = 1, standardize = TRUE, nfolds = 10, type.measure = "mse") # Select penalties by ten-fold cross-validation.
lasso_at_min <- glmnet(x, y, family = "gaussian", alpha = 1, standardize = TRUE, lambda = cv_lasso$lambda.min) # Fit at lambda.min using the lecture's separate-fit procedure.
lasso_at_1se <- glmnet(x, y, family = "gaussian", alpha = 1, standardize = TRUE, lambda = cv_lasso$lambda.1se) # Fit at lambda.1se using the same procedure.
lasso_min <- as.matrix(coef(lasso_at_min))[, 1] # Extract estimates from the lambda.min fit.
lasso_1se <- as.matrix(coef(lasso_at_1se))[, 1] # Extract estimates from the lambda.1se fit.
chosen_lambdas <- c(cv_lasso$lambda.min, cv_lasso$lambda.1se) # Collect the two selected penalties.
cv_indices <- match(chosen_lambdas, cv_lasso$lambda) # Locate their cross-validation results.
cv_table <- data.frame(Rule = c("lambda.min", "lambda.1se"), Lambda = chosen_lambdas, CV_MSE = cv_lasso$cvm[cv_indices], CV_SE = cv_lasso$cvsd[cv_indices]) # Summarize the selected penalties.
knitr::kable(cv_table, digits = 6) # Display penalties and cross-validation errors.| Rule | Lambda | CV_MSE | CV_SE |
|---|---|---|---|
| lambda.min | 0.035671 | 0.581824 | 0.095724 |
| lambda.1se | 0.251649 | 0.675031 | 0.133129 |
coefficient_table <- cbind(OLS = coef(full_prostate), Ridge = ridge_coef, LASSO_min = lasso_min, LASSO_1se = lasso_1se) # Combine all coefficient estimates.
knitr::kable(coefficient_table, digits = 6) # Display the coefficient comparison.| OLS | Ridge | LASSO_min | LASSO_1se | |
|---|---|---|---|---|
| (Intercept) | 0.669337 | 0.489814 | 0.601327 | 1.512110 |
| lcavol | 0.587022 | 0.517343 | 0.527788 | 0.446411 |
| lweight | 0.454467 | 0.443616 | 0.386505 | 0.083424 |
| age | -0.019637 | -0.015545 | -0.006887 | 0.000000 |
| lbph | 0.107054 | 0.096111 | 0.072558 | 0.000000 |
| svi | 0.766157 | 0.696395 | 0.596599 | 0.272045 |
| lcp | -0.105474 | -0.042807 | 0.000000 | 0.000000 |
| gleason | 0.045142 | 0.060521 | 0.000000 | 0.000000 |
| pgg45 | 0.004525 | 0.003520 | 0.002304 | 0.000000 |
With these folds, \(\lambda_{\min}=0.035671\) and \(\lambda_{1\mathrm{se}}=0.251649\). Their cross-validation MSEs are 0.581824 and 0.675031, respectively.
At lambda.min, LASSO retains lcavol,
lweight, age, lbph,
svi, and pgg45; it sets lcp and
gleason to zero. At lambda.1se, it retains
only lcavol, lweight, and svi,
with coefficients 0.446411, 0.083424, and 0.272045. These three slopes
are smaller than their OLS, ridge, and lambda.min
estimates. OLS and ridge retain all eight predictors.
5. Refit the predictors selected at lambda.1se
selected_terms <- setdiff(names(lasso_1se)[lasso_1se != 0], "(Intercept)") # Extract the nonzero LASSO predictors.
refit <- lm(reformulate(selected_terms, response = "lpsa"), data = prostate) # Refit the selected predictors by ordinary least squares.
summary(refit) # Print the refitted estimates and model statistics.##
## Call:
## lm(formula = reformulate(selected_terms, response = "lpsa"),
## data = prostate)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.7296 -0.4576 0.0281 0.4640 1.5701
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.2681 0.5435 -0.49 0.623
## lcavol 0.5516 0.0747 7.39 6.3e-11 ***
## lweight 0.5085 0.1502 3.39 0.001 **
## svi 0.6662 0.2098 3.18 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.717 on 93 degrees of freedom
## Multiple R-squared: 0.626, Adjusted R-squared: 0.614
## F-statistic: 52 on 3 and 93 DF, p-value: <2e-16refit_comparison <- data.frame(Term = names(coef(refit)), Full_OLS = coef(full_prostate)[names(coef(refit))], LASSO_1se = lasso_1se[names(coef(refit))], Refit_OLS = coef(refit)) # Compare coefficients retained in the refit.
knitr::kable(refit_comparison, digits = 6, row.names = FALSE) # Display the retained coefficient estimates.| Term | Full_OLS | LASSO_1se | Refit_OLS |
|---|---|---|---|
| (Intercept) | 0.669337 | 1.512110 | -0.268093 |
| lcavol | 0.587022 | 0.446411 | 0.551638 |
| lweight | 0.454467 | 0.083424 | 0.508541 |
| svi | 0.766157 | 0.272045 | 0.666158 |
The refitted model is
\[ \widehat{\mathrm{lpsa}} = -0.26809 + 0.55164\,\mathrm{lcavol} + 0.50854\,\mathrm{lweight} + 0.66616\,\mathrm{svi}. \]
Compared with the full model, the lcavol and
svi slopes decrease, while the lweight slope
increases. The refitted slopes are larger than the
lambda.1se estimates. The refit uses three predictors
instead of eight; its \(R^2\) is 0.6264
versus 0.6548, adjusted \(R^2\) is
0.6144 versus 0.6234, and residual standard error is 0.7168 versus
0.7084.