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)


There's a lot to unpack here, and pardon for the answer being somewhat frank, but I think there's a lot to learn in general regarding the approach and methods. Overall I think there are some good attempts here, but what's lacking is an understanding of why the methods were developed,  how to test them in different scenarios, and how to report the output. In short, the issue has nothing to do with whether MSE is the right accuracy metric, it likely is.

*

*You're exploring models built for other purposes in new ways, so you should think like a scientist and be prepared to be proven wrong about your hypotheses. For instance, the ridge regression was developed with a "large P" in mind. Changing metrics to get answers you like is not going to give you reliable (or reproducible) answers.


*Your code needs major work. We probably can't do that here. But right off the bat, you have this simple function written incorrectly.

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


It is mean squared error... not median squared error, assuming the x matrix is actually errors. There are other issues too, mainly that the code is really inefficient because you're "growing" vectors and matricies as a C programmer might do, rather than vectorizing the whole thing.


*I have a snippet of code on my end demonstrating indeed that while ridge beats OLS when $p$ is large relative to $n$, it doesn't handle collinearity any better than OLS. In fact, it is somewhat worse even for very correlated predictors. The issue with $p$ large never had to do with the collinearity, although it's guaranteed the issue will likely present itself. The issue has to do with the Neyman Scott problem: that is that when the number of free parameters to be estimated is large relative to $n$, likelihood based inference becomes unreliable. This was observed by Stein when he developed the Ridge estimator as a biased alternative to the ANOVA. Reference here: What is the problem in the Neyman-Scott problem?
In scenario 1, Ridge does 50% better than OLS at minimizing MSE as expected. However, if we focus on just 2 collinear predictors, there's no practical difference between the two. (if anything, Ridge does a bit worse).
library(glmnet)

set.seed(123)

## case 1: large p
do.one <- function() {
n <- 2000
p <- 500
b <- rnorm(p) ## vector of coefficients
x <- matrix(rnorm(n*p), n, p)
y <- x%*%b + rnorm(n)

## train and validate indicators
t <- sort(sample(seq(n), n/2))
v <- setdiff(seq(n), t)

m1 <- lm(y ~ x, subset=t)
m2 <- cv.glmnet(x=x,y=y, subset=t, alpha=0)

yhat1 <- cbind(1, x[v, ]) %*% coef(m1)
yhat2 <- predict(m2, newx = x[v, ])

mse1 <- var(c(yhat1 - y[v]))
mse2 <- var(c(yhat2 - y[v]))
c(mse1, mse2)
}

out1 <- replicate(100, do.one())

## case 2: small p, collinear x
do.two <- function() {
n <- 2000
p <- 2
b <- rnorm(p) ## vector of coefficients
Sigma <- matrix(c(1, 0.9, 0.9, 1), 2, 2)
x <- matrix(rnorm(n*p), n, p) %*% chol(Sigma)
y <- x%*%b + rnorm(n)

## train and validate indicators
t <- sort(sample(seq(n), n/2))
v <- setdiff(seq(n), t)

m1 <- lm(y ~ x, subset=t)
m2 <- cv.glmnet(x=x,y=y, subset=t, alpha=0)

yhat1 <- cbind(1, x[v, ]) %*% coef(m1)
yhat2 <- predict(m2, newx = x[v, ])

mse1 <- var(c(yhat1 - y[v]))
mse2 <- var(c(yhat2 - y[v]))

c(mse1, mse2)
}

out2 <- replicate(100, do.two())

