Regression Performance Measures: Alternatives to MSE

Question

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)

Answer 1

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.

Answer 2

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.

1. 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.

2. 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.

3. 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())

Answer 3

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.