I am trying to recreate FIGURE 3.6 from Elements of Statistical Learning. The only information about the figure is included in the caption.
To recreate the forward stepwise line my process is as follows:
For 50 repetitions:
- Generate data as described
- Apply forward stepwise regression (via AIC) 31 times to add variables
- Calculate the absolute difference between each $\hat{\beta}$ and its corresponding ${\beta}$ and store results
The leaves me with a $50 \times 31$ matrix of these differences on which I can calculate the mean of column wise to produce the plot.
The above approach is incorrect but it is not clear to me what exactly it is supposed to be. I believe my issue is with the interpretation of the mean squared error on the Y axis. What exactly does the formula on the y axis mean? Is it just the kth beta being compared?
Code for reference
Generate data:
library('MASS')
library('stats')
library('MLmetrics')
# generate the data
generate_data <- function(r, p, samples){
corr_matrix <- suppressWarnings(matrix(c(1,rep(r,p)), nrow = p, ncol = p)) # ignore warning
mean_vector <- rep(0,p)
data = mvrnorm(n=samples, mu=mean_vector, Sigma=corr_matrix, empirical=TRUE)
coefficients_ <- rnorm(10, mean = 0, sd = 0.4) # 10 non zero coefficients
names(coefficients_) <- paste0('X', 1:10)
data_1 <- t(t(data[,1:10]) * coefficients_) # coefs by first 10 columns
Y <- rowSums(data_1) + rnorm(samples, mean = 0, sd = 6.25) # adding gaussian noise
return(list(data, Y, coefficients_))
}
Apply forward stepwise regression 50 times:
r <- 0.85
p <- 31
samples <- 300
# forward stepwise
error <- data.frame()
for(i in 1:50){ # i = 50 repititions
output <- generate_data(r, p, samples)
data <- output[[1]]
Y <- output[[2]]
coefficients_ <- output[[3]]
biggest <- formula(lm(Y~., data.frame(data)))
current_model <- 'Y ~ 1'
fit <- lm(as.formula(current_model), data.frame(data))
for(j in 1:31){ # j = 31 variables
# find best variable to add via AIC
new_term <- addterm(fit, scope = biggest)[-1,]
new_var <- row.names(new_term)[min(new_term$AIC) == new_term$AIC]
# add it to the model and fit
current_model <- paste(current_model, '+', new_var)
fit <- lm(as.formula(current_model), data.frame(data))
# jth beta hat
beta_hat <- unname(tail(fit$coefficients, n = 1))
new_var_name <- names(tail(fit$coefficients, n = 1))
# find corresponding beta
if (new_var_name %in% names(coefficients_)){
beta <- coefficients_[new_var_name]
}
else{beta <- 0}
# store difference between the two
diff <- beta_hat - beta
error[i,j] <- diff
}
}
# plot output
vals <-apply(error, 2, function(x) mean(x**2))
plot(vals) # not correct
Output:
sqrt(x**2)
is the same as the absolute value ofx
, at the end you are computing the mean absolute error. The mean squared error is computed by omitting thesqrt
call. $\endgroup$