Backwards stepwise regression code in R (using cross-validation as criteria)

Question

I want to write code that does backward stepwise selection using cross-validation as a criterion.

I have only started learning R a month ago and I have almost zero programming experience prior to that. Hence, I would appreciate any comments on the code.

(I understand that there are issues with the backward stepwise selection process itself.)

rm(list=ls())
library(forecast)
library(fpp)
set.seed(1)
x1<-c(1:500)*runif(20,min=0,max=200)
x2<-c(1:500)*runif(20,min=0,max=200)
library(caret)
data(credit)
score <- credit$score 
log.savings <-log(credit$savings+1)
log.income <-log(credit$income+1)
log.address<- log(credit$time.address+1)
log.employed <- log(credit$time.employed+1)
################################## 
#Input
full.model<- score ~ log.savings + log.income + log.address + log.employed + x1 +x2 
# 1. Start with the full model 
counter=1
full <- lm(full.model)
scopevars<- attr(full$terms, "predvars") # list of model predicotrs 
n <- length(full$residuals); # sample size
current_best <- full #what is the currenly the best model ? that would be our full model for now
while(T){ #process until we break the loop 
 best_cv <- CV(current_best)[1] #the CV of the current_best model
 temp <- summary(current_best) #summary output of the current_best model 
 p <- dim(temp$coefficients)[1] # current model's size
 rnames <- rownames(temp$coefficients) # list of terms in the current model
 compare <- matrix(NA,p,1) # create a matrix to store CV information 
 rownames(compare) <- rnames #name the matrix 
 #create a loop that drops one predictor at a time and calculates the CV of the regression when that predictor is dropped 
 #put all the information into "compare"
 #drop the variable that has the lowerst CV when it is dropped.
 for (i in 2:p) { 
 drop_var <- rnames[i] #variable to be deleted 
 f <- formula(current_best) #current formula 
 f <- as.formula(paste(f[2], "~", paste(f[3], drop_var, sep=" - "))) # modify the formula to drop the chosen variable (by subtracting it)
 new_fit <- lm(f) # fit the modified model 
 compare[i]<-CV(new_fit)[1]
 }
 remove_var <-rownames(compare)[which.min(compare)] #drop is the varibale that is associate with the lowerst CV 
 update_cv <- compare[which.min(compare)] 
 if(update_cv>best_cv) break #we should not continue if dropping a variable will not improve Cv 
 write(paste("--- Dropping",counter, remove_var , update_cv, "\n"), file="") # output the variables we are dropping. 
 f <- formula(current_best) #current formula
 f <- as.formula(paste(f[2], "~", paste(f[3], remove_var , sep=" - "))); # modify the formula to drop the chosen variable (by subtracting it)
 current_best <- lm(f)# we now have a new best model 
 counter=counter+1 #update our counter 
 next
} #end the while loop 
print(current_best)
CV(current_best)

Answer 1

Your code is good. I could not find important corrections to be made. I could find only some small cosmetic improvements - nothing serious.

For example. You do not need the caret package loaded for this task. You do not need to use next at the end of the while loop. It is good practice to keep your code in width of 80 characters. It improves readability a lot. See my code attached.

rm(list=ls())
gc()
library(forecast)
library(fpp)
set.seed(1)
x1 <- c(1:500) * runif(20, min=0, max=200)
x2 <- c(1:500) * runif(20, min=0, max=200)
data(credit)
score <- credit$score 
log.savings <- log(credit$savings+1)
log.income <- log(credit$income+1)
log.address <- log(credit$time.address+1)
log.employed <- log(credit$time.employed+1)
################################## 
#Input
full.model <- score ~ log.savings + log.income + log.address +
 log.employed + x1 +x2
# 1. Start with the full model 
counter <- 1L
full <- lm(full.model)
scopevars <- attr(full$terms, "predvars") # list of model predicotrs
n <- length(full$residuals) # sample size
current_best <- full #what is the currenly the best model?
 # that would be our full model for now
while(T) { #process until we break the loop
 best_cv <- CV(current_best)[1] #the CV of the current_best model
 temp <- summary(current_best) #summary output of the current_best model 
 p <- dim(temp$coefficients)[1] # current model's size
 rnames <- rownames(temp$coefficients) # list of terms in the current model
 compare <- matrix(NA,p,1) # create a matrix to store CV information 
 rownames(compare) <- rnames #name the matrix 
 # create a loop that drops one predictor at a time and calculates
 # the CV of the regression when that predictor is dropped 
 # put all the information into "compare"
 # drop the variable that has the lowest CV when it is dropped.
 for (i in 2:p) { 
 drop_var <- rnames[i] #variable to be deleted 
 f <- formula(current_best) #current formula 
 f <- as.formula(paste(f[2], "~", paste(f[3], drop_var, sep=" - ")))
 # modify the formula to drop the chosen variable (by subtracting it)
 new_fit <- lm(f) # fit the modified model 
 compare[i] <- CV(new_fit)[1]
 }
 remove_var <- rownames(compare)[which.min(compare)]
 # drop is the varibale that is associate with the lowest CV
 update_cv <- compare[which.min(compare)]
 if (update_cv > best_cv) break # we should not continue
 # if dropping a variable will not improve CV
 write(paste("--- Dropping", counter, remove_var , update_cv, "\n"), file="")
 # output the variables we are dropping
 f <- formula(current_best) #current formula
 f <- as.formula(paste(f[2], "~", paste(f[3], remove_var , sep=" - ")))
 # modify the formula to drop the chosen variable (by subtracting it)
 current_best <- lm(f) # we now have a new best model
 counter <- counter + 1L # update our counter
} #end the while loop 
print(current_best)
CV(current_best)