0

fold cross-validation method to compare the residuals 19 models fitted using forward selection. I was stuck at the last step. The code is like this:

library(ISLR)
summary(Hitters)  # Use the dataset of Hitters
Hitters = na.omit(Hitters)

library(leaps)

set.seed(11)
folds = sample(rep(1:10,length=nrow(Hitters)))  # used for cross validation later
table(folds)
cv.errors = matrix(NA,10,19)   
# store the errors from 10 validations, each contains an error for a model

# write a prediction function
predict.regsubsets = function(object,newdata,id,...){
  form = as.formula(object$call[[2]])   # extract the formula
  mat = model.matrix(form,newdata)      # extract the exploratory data
  coefi = coef(object,id=id)            # coefficients for the ith model
  return(mat[,names(coefi)]%*%coefi)    # manually get the predicted value
}
# write a function to extract the Mean of squared root of residuals
error = function(object,newdata,origin,num,...){
  pred = lapply(seq_along(1:num),function(x){predict.regsubsets(object,newdata,id=x)})
  sapply(pred,function(x){mean((x-origin)^2)})
}

# this gives error: $ operator is invalid for atomic vectors 
lapply(seq_along(1:10),function(X){
  best.fit = regsubsets(Salary~.,data=Hitters[folds!=X,],nvmax=19,method="forward")
  cv.errors[X,]=error(best.fit,newdata=Hitters[folds==X,],origin=Hitters$Salary[folds==X],num=19)
})

# this works well, except for being slow...
for(X in 1:10){
  best.fit = regsubsets(Salary~.,data=Hitters[folds!=X,],nvmax=19,method="forward")
  cv.errors[X,]=error(best.fit,newdata=Hitters[folds==X,],origin=Hitters$Salary[folds==X],num=19)
}

Thanks!

1 Answers1

0

This might be a scoping problem. Not tried it, but try the following:

lapply(1:10,function(X){
  best.fit = regsubsets(Salary~.,data=Hitters[folds!=X,],nvmax=19,method="forward")
  cv.errors[X,] <- error(best.fit,newdata=Hitters[folds==X,],origin=Hitters$Salary[folds==X],num=19)
})

The only change here is the assignment operator.

Also, not sure if this will solve your problem, but: Do you really need to use seq_along here? Just 1:10 should be enough.

Community
  • 1
  • 1
Shambho
  • 3,250
  • 1
  • 24
  • 37