Regression Performance Measures: Alternatives to MSE Preface
I am dealing with five Regression models (Ordinary Least square, Least Absolute regression, Huber Regression, MM Estimator, and Ridge Regression).
I try to check which model is more robust to outliers and multicollinearity simultaneously. For this reason, I created a simulation in R software with different scenarios of outliers and multicollinearity.
I used cross-validation to split my data into training and test samples.
My question as follows:
I used MSE, and MAE to check which model provides good results in terms of robustness and dealing with multicollinearity. However, I found that always Ordinary Least square gives the best result with MSE measure and Least Absolute regression with MAE measure.
I think this is due to the similarity of optimization of Ordinary Least square with MSE, and Least Absolute regression with MAE.
It means that I really need an alternative unbiased universal measure to check the performance of these five regression models.
Which Measure do you suggest to solve this problem?
Clarification
What I did is simply generating data from a normal distribution, then I introduce some artificial outliers and multicollinearity to this data. I split the data into training and test samples using cross-validation. I built the five models on training data and then check performance measures on testing data. We repeated this 1000 times and takes the average. OLS has the best results with MSE which is very strange because it supposes to be the worse when outliers and multicollinearity presents. The same for Least Absolute regression with MAE measure.
R-code
rm(list=ls())

library(MASS)
library(glmnet)

### Calling the important functions ###

# Mean Square meausre: MSE#
mse=function(x){
  mmm=rep(0,ncol(x))
  for (i in 1:ncol(x)){
    mmm[i]=mean((x[,i])^2)
  }
  return(mmm)
}

# Mean Absloute Deviation measure: MAD#
mad=function(x){
  mmm=rep(0,ncol(x))
  for (i in 1:ncol(x)){
    mmm[i]=mean(abs(x[,i]))
  }
  return(mmm)
}

# mean of the results ##
mee=function(x){
  mmm=rep(0,ncol(x))
  for (i in 1:ncol(x)){
    mmm[i]=mean((x[,i]))
  }
  return(mmm)
}

umar <- function(R,n,sig,p,po,py,fx,fy){
  #' where 'R is the level of multicollinearity between 0 and 1'#
  #' "n" is the sample size
  #' "sig" is the error vatiance
  #' "p" is the number of explanaitory variable
  #' 'po' is percentage outlier in x direction
  #'  'py' is percentage outlier in y direction
  #' 'fx' is magnitude of outlier in x direction
  #' 'fy' is magnitude of outlier in y direction'#
  #' RR' is the number of replication 
  
  RR=20      
  set.seed(123)
  
  OP2=NULL
  OP3=NULL
  
  #explanatory vriables
  
  x=matrix(0,nrow=n,ncol=p)
  W <-matrix(rnorm(n*(p+1),mean=0,sd=1), n, p+1)  
  for (i in 1:n){
    for (j in 1:p){
      x[i,j] <- sqrt(1-R^2)*W[i,j]+(R)*W[i,p+1];      # Introduce multicollinearity
    }    
  }
  
  b=eigen(t(x)%*%x)$vec[,1]
  
  #Invoking outlier
  rep1=sample(1:n, size=po*n, replace=FALSE)
  x[rep1,2]=fx*max(x[,2])+x[rep1,2]     # the point of outlier
  for (i in 1:RR){
    u=rnorm(n,0,sig)
    y=x%*%b+u
    rep2=sample(1:n, size=py*n, replace=FALSE)
    y[rep2]=fy*max(y)+y[rep2]
    
    dat=data.frame(y,x)
    n=nrow(dat)
    
    # K-fold Cross validation
    #Create k equally size folds
    
    k=3 # number of folds
    folds <- cut(seq(1,n),breaks=k,labels=FALSE)
    
    mols=matrix(0,nrow= k);
    mM=matrix(0,nrow= k);mMM=matrix(0,nrow= k);
    mrls=matrix(0,nrow= k);mrm=matrix(0,nrow= k);mrmm=matrix(0,nrow= k);

    mols2=matrix(0,nrow= k);
    mM2=matrix(0,nrow= k);mMM2=matrix(0,nrow= k)
    mrls2=matrix(0,nrow= k);mrm2=matrix(0,nrow= k);mrmm2=matrix(0,nrow= k);
    
    #Perform 3 fold cross validation
    
    for(i in 1:k){
      #Segement your data by fold using the which() function 
      testIndexes <- which(folds==i,arr.ind=TRUE)
      testData <- dat[testIndexes, ]
      trainData <- dat[-testIndexes, ]
      xtr=as.matrix(trainData[,-1])
      ytr=trainData[,1]
      xte=as.matrix(testData[,-1])
      yte=testData[,1]
      
      mest=rlm(ytr~xtr,psi=psi.huber,k2=1.345,maxit=1000)$coefficients  # Huber Regression 
      
      mmest=rlm(ytr~xtr,method="MM",maxit = 1000)$coefficients  # MM Estimators 
      
      ols=lm(ytr~xtr)$coefficients     # OLS Regression 
      
      nxtr=model.matrix(~xtr)



      ridge.fit.cv <- cv.glmnet(nxtr, ytr, alpha = 0, standardize = FALSE, intercept = TRUE)
      ridge.fit.lambda <- ridge.fit.cv$lambda.1se
      
      I=diag(1,ncol(nxtr))
      ridols=solve(t(nxtr)%*%nxtr+ridge.fit.lambda*I)%*%(t(nxtr)%*%nxtr)%*%ols  # Ridge Regression 
      mrls[i]=mean(yte-cbind(1,xte)%*%ridols)^2
      ridM=solve(t(nxtr)%*%nxtr+ridge.fit.lambda*I)%*%(t(nxtr)%*%nxtr)%*%mest # Ridge Huber 
      mrm[i]=mean(yte-cbind(1,xte)%*%ridM)^2
      ridMM=solve(t(nxtr)%*%nxtr+ridge.fit.lambda*I)%*%(t(nxtr)%*%nxtr)%*%mmest # Ridge MM
      mrmm[i]=mean(yte-cbind(1,xte)%*%ridMM)^2


      mols[i]=mean(yte-cbind(1,xte)%*%ols)^2
      mM[i]=mean(yte-cbind(1,xte)%*%mest)^2
      mMM[i]=mean(yte-cbind(1,xte)%*%mmest)^2
      
      mrls2[i]=mean(abs(yte-cbind(1,xte)%*%ridols))
      mrm2[i]=mean(abs(yte-cbind(1,xte)%*%ridM))
      mrmm2[i]=mean(abs(yte-cbind(1,xte)%*%ridMM))
      mols2[i]=mean(abs(yte-cbind(1,xte)%*%ols))
      mM2[i]=mean(abs(yte-cbind(1,xte)%*%mest))
      mMM2[i]=mean(abs(yte-cbind(1,xte)%*%mmest))
      
    }
    
    res1=cbind(mols,mM,mMM,
               mrls,mrm,mrmm)
    
    res3=cbind(mols2,mM2,mMM2,
               mrls2,mrm2,mrmm2)
    
    op2=mse(res1)
    OP2=cbind(OP2,op2)
    op3=mad(res3)
    OP3=cbind(OP3,op3)
    
  }
  
  MSE=mee(t(OP2))
  MAD=mee(t(OP3))
  
  
  
  nam=c("OLS","M","MM","Ridge-OLS","Ridge-M","Ridge-MM")
  
  data.frame(nam,R,n,sig,p,po,py,fx,fy,MAD,MSE)
}


results=NULL
R=c(0.999)
n=c(100)
sig=c(5)
p=c(3)
po=c(0.2)
py=c(0.2)
fx=c(5)
fy=c(5)

for(i in 1:length(R)){
  for(j in 1:length(n)){
    for(k in 1:length(sig)){
      for(l in 1:length(p)){
        for(m in 1:length(po)){
          for(nn in 1:length(py)){
            for(o in 1:length(fx)){
              for(pp in 1:length(fy)){
                results=rbind(results,umar(R=R[i],n=n[j],sig=sig[k],p=p[l],
                                           po=po[m],py=py[nn],fx=fx[o],fy=fy[pp]))
              }
            }
          }
        }
      }
    }
  }
}

View(results)


The need for an alternative for the MSE is often discussed in the context of scoring rules, when one looks for an alternative for the Brier score (which is MSE in its way: $B=\frac{1}{N}\sum_{i}{(y_i-\hat{\pi}_i)^2}$. One common alternative is the log scoring rule, $L=\frac{1}{N}\sum_{i}{(y_i log(\hat{\pi}_i) + (1-y_i)log(1-\hat{\pi}_i))}$.
I'm not very sure about it, but let me offer you a direction:
Denote $p_i=\frac{|\hat{y}_i-y_i|}{|y_i|}$ as the relative error in percentage, then take $\frac{1}{N}\sum_{i}{p_i}$ as a measure of error.
Again, not sure where will it lead you, but it might help as it takes relative errors (as percentages) and then their log, this is a rather different scale and will be less sensitive to outliers.
