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 performance measure should reflect what you actually care about when making predictions. E.g. choose MSE if you think that the badness/cost of a prediction increases quadratically with its distance from the true value, or MAE if you think it increases linearly. Other factors might come into play as well, like choosing a common metric in order to compare your results to past work. In all cases, you should choose the performance measure (or multiple measures) ahead of time, based on your actual interests. Then, see how the results turn out. You shouldn't fiddle with the metric in order to get an answer you've already decided is true.
Some other considerations:

*

*The relative performance of the different regression methods can depend on the the fraction of points that are outliers, as well as how the outliers are distributed (e.g. how far they tend to lie from the inliers). It's worth exploring a range of these parameters in your code, because the results you find for a single setting might not apply to another. Robust regression methods would be expected to show a greater advantage as the frequency and severity of outliers increase.


*When measuring test set prediction error, it's worth thinking about whether the test set should contain outliers. If yes, you're asking about predictive performance on the same distribution as the training data (a 'contaminated' mix of inliers and outliers). If no, you're asking how well a model trained on a contaminated mix fits the inliers alone (which is presumably the distribution of interest). Both are valid approaches, but they address different questions. Outlier-free test sets would probably tend to show a greater difference between robust and non-robust regression methods.


*Besides prediction, you might also be interested in inference. I.e. how closely can you recover the true parameters of the (inlier) data generating process? For example, suppose the inliers are generated as $y_i = \beta^T x_i + \epsilon_i$ and you fit coefficients $\hat{\beta}$ to a training set containing a mix of inliers/outliers. The squared Euclidean distance $\|\beta - \hat{\beta}\|_2^2$ is one way to quantify how well the estimated coefficients match the true values.
