eliminating outliers in MARS regression

Question

I using the regression method called MARS, in R is it called earth and is located in the package earth, in order to find the best regression model for my datat.

I know that this method is suitable for large data-sets, can handle NA and also decides which variables will be used and which not into the regression.

What I'm doing

After the regression is estimated, I detect the outliers using boxplot and then I eliminate from the data the observations which are extreme values and compute the model again.

I do this until maximum of grsq and rsq are found.

CODE

model <- earth(log(price) ~ ., data = data, weights = weights)
max_grsq <- round(model$grsq, digits = 4)
    max_rsq <- round(model$rsq, digits = 4)
min_diff <- abs(max_grsq - max_rsq)

while(!done) {
  residuals_abs <- abs(model$residuals)
      boxplot <- boxplot(residuals_abs, plot=F)
      indexes_to_remove <- c(which((residuals_abs > boxplot$stats[4]) == T), which((residuals_abs < boxplot$stats[2]) == T))

  if (length(indexes_to_remove) > 0) {
    data <- data[-indexes_to_remove, ]
    distances <- distances[-indexes_to_remove]
    weights <- (1/distances)/(sum(1/distances))
  }

  tempModel <- earth(log(price) ~ ., data = data, weights = weights)
  temp_grsq <- round(tempModel$grsq, digits = 4)
      temp_rsq <- round(tempModel$rsq, digits = 4)
  temp_diff <- abs(temp_grsq - temp_rsq)

  if ((temp_grsq > max_grsq && temp_rsq >= max_rsq) || (temp_grsq >= max_grsq && temp_rsq > max_rsq)) {
    model <- tempModel
    max_grsq <- temp_grsq
    max_rsq <- temp_rsq
    min_diff <- temp_diff
  } else {
    done = T
  }
 }

QUESTION

I'm not a statistician so I don't know any better way for removing the outliers.

• is my approach correct?
• should I use another approach?
• I know that there are bad outliers and good outliers (leverage points), how can I remove only the bad outliers?
• I'm using the semi-log form of the regression. because of the use of dummy variables I can't use the log-log form. Is there any other approach for data transformation? or should I standardize the data? x <- (x - x_min)/(x_max - x_min)

Does anyone has some hints?

Answer 1

My suggestion is independent of the software used.

We need to clarify whether the outliers are outliers in Y, the dependent variable, or outliers in the predictor, X. Outliers in the predictor, X, are easily handled with the vast numbers of transformations available that would reshape the PDF (probability density function) of X.

While I agree with Eric Farng that deleting outliers in Y is not recommended, I disagree with him that they should be deleted only after "careful consideration." In my opinion, one should never delete outliers since they contain important and useful information, that is, unless one can determine that these values are somehow "bad," illegal or fraudulent, etc. The alternatives to outlier deletion in Y are to leverage modeling methods that are robust to outliers.

Why am I opposed to a priori outlier deletions in Y? Let me use an example: after "careful consideration" (Eric Farng's wording, however one chooses to define this) a first set of outliers is deleted. Does this mean that you are done with deleting outliers? Probably not, since a second analysis would almost certainly reveal a new set of outliers relative to the new mean and standard deviation. What does one do with this new information? And how many passes of the data are required to completely scrub it of outliers? Clearly, this is a potentially endless process of outlier deletion that makes little or no sense.

Most importantly and even before one gets into the almost philosophical question of whether or not to delete outliers, it should be noted that MARS is one of the robust, non-parametric alternatives that exploits nonlinear, quantile functions of the relationship between X and Y. From a purely applied and practical point of view this means that MARS is highly robust, almost immune, to the presence of outliers in Y: by definition, deleting outliers in Y is unnecessary, even moot, when leveraging MARS.

Answer 2

From the R package documentation and from the original MARS paper, it looks like rsq and grsq are used for model selection within the package and it looks like your code is removing outliers until the fit of the model is maximized. This is usually not recommended. There are statistical tools to help identify potential outliers and outliers are removed only after careful consideration. Points that are mistakes (value of -1 for count data) or do not represent the population (a population of adults containing a child's data point) can be removed. As for other points, they may be an unlucky outlier. They may also be unlucky to be only one point and not more. So these are removed only after careful inspection.

In the earth package, residuals and leverage are easily available. As you say, not all points with high leverage are necessarily a problem. Cook's Distance tries to solve that issue by actually removing the point from the data and checking how any predicted values change. Unfortunately, it may be difficult to find a library to compute this and other outlier metrics that support the earth package. However, this package does support the plot(model) command which will also give some potential outliers.

For the log-log form, I believe you can skip the categorical variables if, for a variable with $n$ categories, they are being represented as $n-1$ binary/dummy variables. For example in a simple linear regression, transforming the categorical dummy variables [0, 1] into [0, 0.5] will only cause its coefficient to a double. Here is an example using MARS regression. a is the dependent variable. b, d, e are the independent variables where d and e are categorical variables. In the second example, d and e are changed from (0 and 1) to (0 and 0.5) and (-1, 1) respectively.

a <- rnorm(100)
b <- rexp(100) + a
c <- rgamma(100, shape=1) + 2*a
d <- sapply(sapply(round(b + c), function(x) min(1,x)), function(y) max(0, y))
e <- sapply(sapply(round(a - c), function(x) min(1,x)), function(y) max(0, y))
fit <- earth( a ~ b + d + e)

e2 <- e * 2
e2 <- e2 - 1
d2 <- d * 0.5
fit2 <- earth( a ~ b + d2 + e2)

summary(fit)
summary(fit2)
sum(predict(fit) - predict(fit2))

You can see from the output, the coefficients for d, e and the intercept have changed, but none of the other coefficients changed. In addition, the predictions are the same between models.

For reference:

$$ rsq = R^2\\ grsq = 1 - \frac{gcv}{gcv.null}\\ gcv = \frac{1}{N}\sum_{i=1}^{N}\frac{\big[y_i-\hat{f}_M(x_i)\big]^2}{\bigg[1 - \frac{C(M)}{N}\bigg]^2}\\ C(M) = trace(B(B^TB)^{-1}B^T) + 1 $$

$B$ is the "data" matrix. $gcv.null$ is the Generalized Cross Validation of the intercept only model.

Answer 3

To elaborate slightly on Eric Farng's comment that removing outliers until the fit of the model is maximized is not recommended:

The fundamental problem with tweaking the data until you get a good GRSq is that although you will build a model that gives a good fit to your selected data, your model will not give good predictions for future data --- because you aren't modelling the underlying distribution of the data. It's like a model that predicts stock market prices very accurately for historical data but is useless for prediction of future stock prices.