2
$\begingroup$

There is this Figure 7.14 in Elements of Statistical Learning:

enter image description here

Which I've been trying to make, and so far I've been able to only make the top-left plot "Prediction Error" using a code that would generate the 100 simulations as in figure 7.3 (see image).

enter image description here enter image description here

Figure 7.14 only requires the linear regression plot (top-right) in fig 7.3. The code is shown below (Note:case =2 is for the kNN plot which is not necessary):

# generate dataset
genX <- function(n = 80, p = 20){
  X = matrix(runif(n*p, 0, 1), ncol = p, nrow = n)
  return(X)
}
# generate response 
genY <- function(X, case = 1){
  n = nrow(X)
  Y = numeric(n)
  if (case == 1){ # for the left panel of fig. 7.3
    Y = sapply(X[, 1], function(x) ifelse(x <= 0.5, 0, 1))
  } 
  else {
    Y = apply(X[, 1:10], 1, function(x) ifelse(sum(x) > 5, 1, 0))  
  }
  return(Y)
}

## global parameters setting
ntest = 1000
percent = 0.75
B = 100 # the number of repetition

## case 2
seed = 1234
set.seed(seed)
X = genX()
Y = genY(X, case = 2)
X.test = genX(n = ntest)
Y.test = genY(X.test, case = 2)

## use leaps package to do best subset selection
library(leaps)
## predict test dataset by using the best subset model with size p
predict.regsub <- function(model, p, X.test){
  which = summary(model)$which
  coef.raw = coef(model, p)
  # construct coef vector
  if (length(coef.raw) == p+1)
  {
    coef.vec = numeric(1+ncol(X.test)) # include intercept
    coef.vec[1] = coef.raw[1]
    flag = 1
  }
  else
  {
    coef.vec = numeric(ncol(X.test))
    flag = 0
  }
  j = flag + 1 # point to raw coef
  for (i in c(1:ncol(X.test)) + flag){
    if (which[p, i]){
      coef.vec[i] = coef.raw[j]
      j = j + 1
    }
  }
  # for simplicity, consider intercept; and in fact, every regsubset models have intercept
  pred = apply(cbind(1, X.test), 1, function(x) sum(x*coef.vec))
  return(pred)
}

n = nrow(X)
## store all prediction
reg.pred.full = vector("list", 20)
for (i in 1:20){
  reg.pred.full[[i]] = matrix(nrow = nrow(X.test),
                              ncol = B)
}
for (i in 1:B){
  train.id = sample(n, n*percent)
  X.train = X[train.id, ]
  Y.train = Y[train.id]
  reg.sub = regsubsets(X.train, Y.train, nvmax = 20)
   for (j in 1:20){
     reg.pred.full[[j]][, i] = predict.regsub(reg.sub, j, X.test)
   }
}
## calculate bias2, variance, epe
reg.bias2.full = numeric(20)
reg.variance.full = numeric(20)
reg.epe.full = numeric(20)
cl.epe.full = numeric(20)

#apply 
#sapply
for (i in 1:20){
  bias2 = sapply(1:ntest, function(j) (mean(reg.pred.full[[i]][j, ]) - Y.test[j])^2)
  variance = apply(reg.pred.full[[i]], 1, function(x) var(x))
  # for regression
  epe = sapply(1:ntest, function(j) mean((reg.pred.full[[i]][j, ] - Y.test[j])^2))
  reg.variance.full[i] = mean(variance)
  reg.bias2.full[i] = mean(bias2)
  reg.epe.full[i] = mean(epe)
  }

Here, is where I begin to code the top-left plane of Figure 7.14.

#Prediction Error
error = matrix(0,nrow=20,ncol=100);
for (i in 1:20){
  error[i,] = sapply(1:100, function(j) mean((reg.pred.full[[i]][,j] - Y.test)^2))
}

Finally, this is where I got stuck. I am attempting the top-right panel for the 10-fold-CV but I do not know what to do next after I use part of the available data to fit the model. Same goes for the bottom plots.

#10-fold-cross-validation
n_folds<-10;
folds_i<-sample(rep(1:n_folds,length.out = B));
cv_err = matrix(0,nrow=n_folds,ncol=20);
for (k in 1:n_folds){
  test_i<-which(folds_i == k)
  #what's next??

}

#Leave-one-out CV

I know I have to fill the the cv_err to make a plot like the one for the Prediction error shown below,

#plots
plot(1:20,error[,1],"l",ylim=c(0,.4),xlab = "Subset size p", ylab = "Error", main = "Prediction Error");
colvec<-rainbow(99);
for(i in 2:100)
{
lines(1:20,error[,i],"l",col=colvec[i]);  
}

enter image description here

But I'm missing the point, as to how to use the 10 folds with my predicting models. Also, from the 10-fold, how to do the leave-one out and the bottom right panel. Is there anyone who could help me with this? Note: I have also posted the same question in StackOverflow in case this is more of a code issue rather than a stats issue.

$\endgroup$

0