3
$\begingroup$

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)
$\endgroup$
14
  • 3
    $\begingroup$ OLS is, by design, the best at minimizing MSE. The parameters are the ones that minimize MSE, so no other parameter vector can beat it. Likewise, least absolute regression is designed to minimize MAE. It is not just that the optimization is similar to the performance assessment; the optimization is identical to the performance metric. // Are you doing out-of-sample testing? Then it can be a different story. // It is up to you to decide what you value for performance. If you want to minimize square loss, least absolute regression does not seek such a solution. $\endgroup$
    – Dave
    Jul 3 at 22:01
  • 1
    $\begingroup$ @Dave, 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. $\endgroup$
    – jeza
    Jul 5 at 8:34
  • 1
    $\begingroup$ OLS is designed to give the minimum square loss. That’s what “least squares” means. $\endgroup$
    – Dave
    Jul 5 at 11:55
  • 2
    $\begingroup$ @Dave No, OLS is best among unbiased estimators at minimizing MSE, but there exist many many biased estimators that have better MSE. Stein's example gives us ridge regression as an optimal ANOVA. Bayes optimal estimators generally do better than OLS $\endgroup$
    – AdamO
    Jul 6 at 14:09
  • 2
    $\begingroup$ I'm a bit confused. Can you share the code? One possible explanation is that your "multicollinearity" and "outliers" aren't as badly behaved as you are simulating them. It doesn't make sense to change the performance metric on which you are judging a model to be "good" because it is not giving you the results you expected. $\endgroup$
    – AdamO
    Jul 6 at 14:13
7
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ I tried different scenarios of outliers and multicollinearty but still, I have the same problem. $\endgroup$
    – jeza
    Jul 7 at 8:14
  • 1
    $\begingroup$ @jeza If everything has been done correctly, then it's not a problem; it's a result. But, your description doesn't match what I see when quickly comparing OLS to LAD. In my simulation, LAD can show a big advantage in test set MSE in the presence of strong outliers. Perhaps we're doing things differently. But, it's also worth checking your code to ensure that it's not an implementation issue. For example, the mean vs. median problem that AdamO pointed out is quite important, and also makes me wonder whether there might be mistakes in other places. I don't use R, so haven't looked at your code. $\endgroup$
    – user20160
    Jul 9 at 14:07
5
$\begingroup$

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.

  1. 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())
$\endgroup$
2
  • $\begingroup$ I edit this function but still have the same issue. 'mse=function(x){ mmm=rep(0,ncol(x)) for (i in 1:ncol(x)){ mmm[i]=median((x[,i])^2) } return(mmm) }' $\endgroup$
    – jeza
    Jul 7 at 8:03
  • $\begingroup$ @jeza you are really under prepared to tackle this problem. You need to break it down into smaller steps because you're in way too deep. I'm shocked you still didn't understand that you need to use the MEAN as in the function mean to calculate mean-squared error, and not the median. Alas this is just one problem, as I alluded to earlier. $\endgroup$
    – AdamO
    Jul 7 at 15:13
1
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.