quantile regression (cross validation in R) How can I do a cross validation for quantile regression in R or any other methods to compare between multiple regression and quantile regression?
 A: Ordinary linear regression uses MSE as loss function, while median regression uses mean absolute error as loss. Thus, you can expect the linear regression to do better than median regression regarding MSE (and vice versa for the mean absolute error).
Still, I have put together a couple of lines how to do something like this in R in a relatively generic way without meta learning packages like caret or mlr3. Creating random folds is simple, still I will use my splitTools package that can also do stratified or grouped splitting if you are working on real-world problems.
library(splitTools)
library(ggplot2)
library(quantreg)

# Function to calculate RMSE
rmse <- function(y, pred) {
  sqrt(mean((y - pred)^2))
}

# We will create a model for log-diamonds prices
dim(diamonds)
summary(diamonds)

dia <- data.frame(diamonds)
dia$log_price <- log(dia$price)

# Create fold indices
nfolds <- 5
folds <- create_folds(diamonds$price, 
                      k = nfolds,
                      type = "basic", # type = "stratified"
                      seed = 181)

# Fit the models on each fold and store rmse on out-of-fold
model_equation <- log_price ~ log(carat) + factor(color) + factor(cut) + factor(clarity)

rmse_lm <- rmse_qr <- c()

for (fold in folds) {
  fit_lm <- lm(model_equation, data = dia[fold, ])
  fit_qr <- rq(model_equation, data = dia[fold, ])
  
  pred_lm <- predict(fit_lm, newdata = dia[-fold, ])
  pred_qr <- predict(fit_qr, newdata = dia[-fold, ])
  
  rmse_lm[length(rmse_lm) + 1] <- rmse(dia[-fold, "log_price"], pred_lm)
  rmse_qr[length(rmse_qr) + 1] <- rmse(dia[-fold, "log_price"], pred_qr)
}

# Summarize cross-validation results
X <- data.frame(Model = c("LM", "QR"),
                RMSE = c(mean(rmse_lm), mean(rmse_qr)),
                sd = c(sd(rmse_lm), sd(rmse_qr)))

#  Model      RMSE          sd
#1    LM 0.1338296 0.001471568
#2    QR 0.1349292 0.001232364

ggplot(X, aes(x = Model, y = RMSE)) + 
  geom_bar(stat = "identity", fill = "darkred", width = 0.3) +
  geom_errorbar(aes(ymin = RMSE - sd, ymax = RMSE + sd, width = 0.1))


The two models show almost the same average performance on the 5 out-of-fold data slices. Still, as expected, linear regression wins. Since we get a performance score for each fold, we can get an impression on the variation of performance across folds. Here, it is very small.
