Robust methods and penalized regression

Question

Are penalized regression methods such as ridge or lasso sensitive to outliers? If so, what options are there in regards to robust methods for penalized regressions and are there any packages in R?

Answer 1

The majority of the development of Ridge and LASSO relates to estimation of OLS parameters. Recent work has expanded this to GLMs for exponential families, unified under the notion that it's the likelihood that's penalized.

In robust statistics, one views maximum likelihood as a special case of the general optimization problem for a general loss function: $\sum_i^n \rho(\mathbf{X}_i\beta - Y_i)$. Different loss functions can be defined to give robust estimators. For instance, $\rho(r) = \frac{r^2}{1+r^2}$ is locally quadratic but bounded above by 1, having the effect of "downweighting" large residuals. This is a good alternative to the OLS estimator which uses unbounded quadratic loss. With minimax estimation, one can often show these estimators are the "best" according to some arbitrarily defined risk functions.

I don't see any reason why the general minimax estimator would not be amenable to $\mathcal{L}_1$ or $\mathcal{L}_2$ penalties. As far as I know, however, substantive research hasn't been done with the problem of estimation in this case.

$$\hat{\beta} = \mbox{argmin}_\beta \left\{ \sum_i^n \rho(\mathbf{X}_i\beta - Y_i) - \lambda \| \beta\| \right\}$$

So play around with it with your own optimization functions and see what you find. Make good use of the R functions optim and nlm. I don't think there are any packages for this.

Answer 2

Sure, you can combine $l_1$ (or $l_2$) penalty with robust regression.

Consider for example Alfons et al. 2013 [0] which combines $l_1$ sparsity penalty with the LTS loss function (and a FastLTS like algorithm). Their Lasso-LTS estimator is defined as:

$$(1)\quad\hat{\pmb\beta}_{\text{LLTS}} = \arg\min_{\pmb\beta}\sum_{i=1}^h S(\pmb y - \pmb X\pmb\beta)_{(i)}+h\lambda||\pmb\beta||_1$$

where $n$ indicates the number of observations, $p$ the number of design variables, $S$ is a symmetric, smooth, positive loss function, $h$ an integer larger than $[n/2]+1$ and for any $n$ vector $\pmb y$, $S(\pmb y)$ is the $n$-vector obtained by applying $S$ to the entries of $\pmb y$ element-wise, $S(\pmb y)_i$ and $S(\pmb y)_{(i)}$ are respectively the $i$-th and the $i$-th largest entry of this vector so that $\sum_{i=1}^hS(\pmb y)_{(i)}$ is the sum of the $h$ smallest entries of $S(\pmb y)$.

In this notation, $h$ is a parameter that governs the desired robustness of the estimator to the presence of outliers in the data and $S$ is your usual loss function (for example squared loss). The robustness to outliers comes from the partial sum in $(1)$. Since $h>[n/2]+1$, the use of a partial sum prevents observations inconsistent with the multivariate pattern of the bulk of the data from influencing the fit (for a online, R code based treatment of robust regression see this tutorial, a more textbook-y treatment of robust estimation can be found in [2]).

Here is a link to a high quality R implementation of the Lasso-LTS method.

The Lasso-LTS method is robust in the sense of having a bounded loss function and a positive and high breakdown point. The (finite sample) breakdown point of an estimator is a pragmatic measure of its robustness to the presence of outliers in the data [1]. Informally, it is the smallest proportion of the original data that needs to be replaced by outliers to drive the estimates arbitrary far away from the values they would have had on the original data. Intuitively, the higher the breakdown point the more robust the estimator (to give you an idea, the breakdown point of the classical lasso is essentially the same as that of the univariate mean: 0). For the Lasso-LTS, the finite sample breakdown point is:

$$\varepsilon^*_n(\hat{\pmb\beta}_{\text{LLTS}})=\frac{n-h+1}{n}\approx0.5,$$ (for comparison this is essentially similar to the breakdown point of the univariate median) Typically the robustification of the fit penalty (the first term) of the total loss function renders it highly non convex and very complex (though smooth) in the sense of having a large number of local minimae. Moreover, these minimae are typically disconnected from one another (they don't say all lie on a manifold). This makes the search for a solution to $(1)$ a delicate affair. For this reasons, special algorithms have been developed to obtain (stochastic) approximations to the corresponding optimae. A big issue here revolves around data based procedures to pick good starting points. I wouldn't advice trying to rewrite your own version of them (or trying to reinvent the wheel).

• [0] A. Alfons, C. Croux, S. Gelper (2013). Sparse least trimmed squares regression for analyzing high-dimensional large data sets. The Annals of Applied Statistics, 7(1), 226-248. Ungated Link
• [1] Donoho, D. L., Breakdown properties of multivariate location estimators, Dept. Statistics, Harvard Univ. 1982. Ungated Link
• [2] Maronna R. A., Martin R. D. and Yohai V. J. (2006). Robust Statistics: Theory and Methods. Wiley, New York.