Basic cross-validation in R

Question

We would like to perform n-fold non-stratified cross-validation on a model, say an lm or glm.

We have a list of quality metrics we would like to evaluate, which we wrap in a function:

fn_rmse <- function(v_pred,v_true) sqrt(mean((v_true-v_pred)^2))
fn_mae <- function(v_pred,v_true) mean(abs(v_true-v_pred))
l_q <- list(fn_rmse, fn_mae, fn_smdape)
fn_q <- function(v_pred, v_true) {
  return(sapply(l_q, function(f)f(v_pred,v_true)))
}

Now we set up a cross-validation routine:

# Set up n-fold non-stratified cross-validation routines

# fn_fold: Run cross validation on a fold
# Input: 
#   df_train = training data
#   df_test = testing data
#   mdl = model to test
#   fn_q = quality metrics to try
# Output: list of qualities for this fold
fn_fold <- function(df_train, df_test, mdl, fn_q) {
    mdl_fold <- update(mdl, data=df_train)
    v_predictions <- predict(mdl_fold, newdata=df_test, type="response")[[1]]
    fmla <- formula(mdl$call["formula"])
    s_regressand <- all.vars(fmla[[2]])
    v_quality <- sapply(v_predictions, fn_q, v_true=df_test[[s_regressand]])
    return(v_quality)
}

# fn_cv: Perform n-fold cross validation
# Input:
#   df = sample
#   mdl = model
#   n_folds = number of folds
#   fn_q = quality metrics to try
# Output: ([# of metrics] x n_folds)-matrix containing cross validation results
fn_cv <- function(df, mdl, n_folds, fn_q) {
    v_folds <- cut(seq(1,nrow(df)), breaks=n_folds, labels=FALSE)
    v_folds <- sample(v_folds)
    v_results <- sapply(1:n_folds,
        function(n_f) fn_fold(
            df[-which(v_folds==n_f,arr.ind=TRUE),],
            df[which(v_folds==n_f,arr.ind=TRUE),],
            mdl, 
            fn_q)
        )
    return(v_results)
}

To run 10-fold cross-validation on a model, I call:

mtx_cv_results <- fn_cv(df=df_my_sample, 
                        mdl=mdl_my_model,
                        n_folds=10,
                        fn_q=fn_q)

This feels a bit like reinventing the wheel. Can anything be done to make this code nicer?

Answer 1

My advice is to use the caret-package in which you can easily define your own metrics. The use of CV is described here in a clear manner:

https://machinelearningmastery.com/how-to-estimate-model-accuracy-in-r-using-the-caret-package/

The ingredients to your metrics requirements should be deriveable from the output via str(model). Here's a small example that should help you to get to your results:

library(caret)
library(ModelMetrics)

# Generate some data
Xs<- matrix(rnorm(300*20), nrow = 300, ncol = 20)
Yvec<- rnorm(300)

# Define your metrics
maeSummary <- function (data,
                        lev = NULL,
                        model = NULL) {
  mae <- mae(data$obs, data$pred)  
  names(mae) <- "MAE"

  mse <- rmse(data$obs,data$pred)
  names(mse) <- "RMSE"

  out = c(mae,mse)
  out
}


# Initialize 
train= trainControl(method="CV", number=10, summaryFunction = maeSummary, returnResamp="all", savePredictions="all")

# Train
model = train(x= Xs, y = Yvec, trControl=train, method="lm")

# Overall results
print(model)

# Metrics per fold
model$resample

Good luck!