Robust regression for highly corrupted response by shifting outliers

Abstract

Outlying observations are often disregarded at the sacrifice of degrees of freedom or downsized via robust loss functions (e.g., Huber's loss) to reduce the undesirable impact on data analysis. In this article, we treat the outlying status of each observation as a parameter and propose a penalization method to automatically adjust the outliers. The proposed method shifts the outliers towards the fitted values, while preserve the non-outlying observations. We also develop a generally applicable algorithm in the iterative fashion to estimate model parameters and demonstrate the connection with the maximum likelihood based estimation procedure in the case of least squares estimation. We establish asymptotic property of the resulting parameter estimators under the condition that the proportion of outliers does not vanish as sample size increases. We apply the proposed outlier adjustment method to ordinary least squares and lasso-type penalization procedure and demonstrate its empirical value via numeric studies. Furthermore, we study applicability of the proposed method to two robust estimators, Huber's robust estimator and Huberized lasso, and demonstrate its noticeable improvement of model fit in the presence of extremely large outliers.

Keywords

case-specific parameter extreme outliers Huber's estimator robust lasso robust linear model

1 Introduction

Extensive research has been conducted to deal with outliers that can impede proper summarization of the overall tendency in a data set. The existing methods range from two-step approaches of first identifying the outliers and then reducing its effect (Kosinski 1999; Pena and Yohai 1999) to robust regression procedures (Bloomfield and Steiger 1983; Huber 1973; Beaton and Tukey 1974; Huber and Ronchetti 2009; Chambers and Heathcote 1981). Along the direction of adopting an explicit error distribution to account for the outlier data, Lange et al. (1989) used a heavy-tailed t-distribution to accommodate the outliers. Hadi and Simonoff (1993) provide a good literature review on the topic of outlier detection and robust regressions. More recently, Lee et al.(2012) introduced a penalization method with the parameters indicating the outlying status of all the observations, called ‘case-specific parameters’, and a lasso penalty (Tibshirani 1996) for the case-specific parameters to identify and downsize the outliers in an automated fashion. Case-specific parameters are devised to deal with each observation and should be penalized to avoid over-saturation. The method of Lee et al. (2012) with case-specific parameter turned out to be identical to Huber's estimator (Huber 1973) and Huberized lasso (Rosset and Zhu 2007), depending on the incorporation of lasso penalty. (See Section 3.1 of Lee et al. (2012) for the details.)

In this work, we are particularly interested in robust modelling of a data set containing massive outliers. Specifically, we propose a new penalty function for the case-specific parameters which are tailored to identify the influential outliers and to effectively reduce their impact on the fitted values. For the asymptotic behaviour of the proposed method, we consider the non-standard situation where the proportion of outlying samples remains non-zero as sample size increases.

We start with a classical linear regression model,

y_{i} = x_{i}^{'} β + ε_{i},

(1.1)

where

ε_{i} s

are independently and identically distributed (iid) with mean 0 and a finite variance

σ^{2}

, and

x_{i}

is a length-p vector for the ith observation. The

ε_{i}

is typically assumed to be normally distributed, which can be easily violated with the presence of large outliers (or leverage points). To address this, we use an additional parameter

γ_{i}

to capture the outlying observations and make the normality assumption more reasonable. Accordingly, the model is modified as

y_{i} = x_{i}^{'} β + γ_{i} + ε_{i}, i = 1, . . ., n,

(1.2)

where

γ_{i}

s are the case-specific parameters (Lee et al. 2012). It takes non-zero values for the outlying observations, and zero for the non-outlying observations. Note that

γ_{i} + ε_{i}

represents the deviation from the true regression line. To make

γ_{i}

and

ε_{i}

identifiable, we impose a constraint of setting

ε_{i}

to zero when its corresponding

γ_{i}

is non-zero. It means that

γ_{i}

is defined to be

y_{i} - x_{i}^{'} β

, if

| y_{i} - x_{i}^{'} β |

is greater than a pre-specified threshold

λ_{γ}

. When

| y_{i} - x_{i}^{'} β | < λ_{γ}

, we set

γ_{i}

to zero, which indicates the ith observation is not an outlier.

When there are massive outliers, diminishing the effect of outliers is often not sufficient. In contrast, the suggested methods have an effect of shifting the detected outliers onto the fitted regression line. Thus, most of the undesirable effects can be removed effectively without deleting the observations, which is our major contribution to the literature. Another advantage could be the broad application of this idea to general modelling procedures beyond ordinary least squares (OLS).

In Section 2.1, the estimation procedure is cast into the convex function optimization (Rockafellar 1997). The estimation of β and $γ = (γ_{1}, \dots, γ_{n})^{'}$ which satisfies the aforementioned conditions is also conducted through an iterative maximum likelihood-type approach, and we show that the two approaches are equivalent. The choice of $λ_{γ}$ is discussed in Section 2.2. The form of penalty for $γ_{i} s$ appears under the convex minimization approach only. The case-specific parameters and the corresponding penalty are applied to the OLS and Huber's robust regression (Huber and Ronchetti 2009) in Section 3.1.

As lasso (Tibshirani 1996) can do variable selection, which is a crucial part of data analysis, we extend the proposed method to lasso and Huberized lasso (Rosset and Zhu 2007) in Section 3. In Section 4, simulation studies are conducted to illustrate the suggested methods and to show the improvement compared to the methods by Lee et al. (2012) and other robust methods. Further simulation studies are contained in the supplementary materials. The applications of the proposed methods to real data sets are given in Section 5, where the first data set is a well-known small data set and the second data set is a larger data set with 3,000 samples.

2 Robust regression by shifting outliers

2.1 Model estimation

We set X to be a design matrix of size n by $p$ with full column rank, and Y to be a response vector of size n. For the purpose of estimating (1.2), we minimize the following objective function:

L (β, γ) = \frac{1}{2} \sum_{i = 1}^{n} (y_{i} - x_{i}^{'} β - γ_{i})^{2} + \sum_{i = 1}^{n} γ_{i}^{2} I (| y_{i} - x_{i}^{'} β | < λ_{γ}) subject to \sum_{i = 1}^{n} γ_{i} ε_{i} = 0 .

(2.1)

The first term on the right side of (2.1) is the measure of goodness of fits, and the second term penalizes

γ_{i}

only when

| y_{i} - x_{i}^{'} β |

is smaller than

λ_{γ}

. Therefore, this new penalty is devised to capture the outliers by shrinking the case-specific parameters

γ_{i}

s towards zero, only for non-outlying observations. On the other hand,

γ_{i}

can be non-zero only if the ith observation is a potential outlier. Since the first and second terms of (2.1) are quadratic, it is clear that

L (β, γ)

is a convex function, which will provide convergency of the suggested iterative algorithm.

Note that the form of the penalty for $γ_{i}$ s is similar to hard thresholding rule. When the hard threshold is used for $β$ , it is identical to best subset selection under the orthonormal design matrix. The proposed penalty is somewhat similar to the best subset selection regarding its format. The difference is that we are selecting the observations instead of variables.

The aforementioned model identifiability constraint for model (1.2) implies that $γ_{i} ε_{i} = 0$ , and $γ_{i}$ and $ε_{i}$ cannot be zero simultaneously. For implementation, we include an addition term of $\sum_{i = 1}^{n} γ_{i} ε_{i} I (| y_{i} - x_{i}^{'} β | < λ_{γ})$ to (2.1). However, it does not affect the value of the original objective function in (2.1) because $γ_{i} ε_{i} = 0$ . Now, as $γ_{i}^{2} + γ_{i} ε_{i} = γ_{i} (y_{i} - x_{i}^{'} β)$ from the model (1.2), the objective function in (2.1) is equivalent to

L (β, γ) = \frac{1}{2} \sum_{i = 1}^{n} (y_{i} - x_{i}^{'} β - γ_{i})^{2} + \sum_{i = 1}^{n} γ_{i} (y_{i} - x_{i}^{'} β) I (| y_{i} - x_{i}^{'} β | < λ_{γ}) .

(2.2)

Now, the minimization of

L (β, γ)

becomes more apparent. We then adopt an iterative algorithm to estimate

β

and

γ

alternately. Due to the iterative fashion, it is noteworthy that the degrees of freedom for estimating regression coefficients in each step are the same as those of OLS without considering variable selection procedure. That is, we adjust the

y_{i} s

with the estimates of

γ

, then find the estimates of

β

with the adjusted

y_{i} s

. Therefore, incorporating the case-specific parameters

γ_{i} s

does not increase model complexity, but requires additional computational cost of an iterative algorithm. However, the number of iteration is minimal (

< 10

) in practice. For the details of algorithm and the adjustment of

y_{i} s

, we first initialize

{\hat{γ}}^{(0)} = 0

, and find

{\hat{β}}^{(0)} =

arg

min_{β} L (β, 0)

, which is just the OLS estimate. We denoted it as

{\hat{β}}^{(0)} = {\hat{β}}^{OLS}

Note that the degrees of freedom (df) for OLS model have been traditionally defined as the number of independent residuals, that is, when we use p variables with sample size n, df is $(n - p)$ without intercept. When we remove m outliers for better fit, the df will be decreased to $(n - m - p)$ for OLS model. In contrast, df of ${OLS}^{S}$ is maintained as $(n - p)$ with possibly removing major effect of m outliers. Of course, this statement holds only when we do not consider variable selection procedure (or, only when we employ the same p covariates for both $OLS$ and ${OLS}^{S}$ ). When we adopt variable selection as a part of modelling, then the selected variables for OLS and ${OLS}^{S}$ can be different, therefore, their df may not be equal. Since the variable selection is embedded when incorporating a lasso penalty, the df of lasso regression can be different from that of the suggested methods with lasso regression.

Next, $\hat{γ}$ can be updated conditional on ${\hat{β}}^{(0)}$ via minimizing

L ({\hat{β}}^{(0)}, γ) = \frac{1}{2} \sum_{i = 1}^{n} (r_{i} - γ_{i})^{2} + \sum_{i = 1}^{n} γ_{i} r_{i} I (| r_{i} | < λ_{γ}),

where,

r_{i} = y_{i} - x_{i}^{'} {\hat{β}}^{(0)}

. By solving

\frac{\partial}{\partial γ_{i}} L ({\hat{β}}^{(0)}, γ) = 0

for each i, we have

{\hat{γ}}_{i}^{(1)} = \{\begin{matrix} r_{i} & if & | r_{i} | \geq λ_{γ}, \\ 0 & if & | r_{i} | < λ_{γ} . \end{matrix}

(2.3)

Now, given

{\hat{γ}}^{(1)}

, we have

L (β, {\hat{γ}}^{(1)}) = \frac{1}{2} \sum_{i = 1}^{n} (y_{i} - x_{i}^{'} β - {\hat{γ}}_{i}^{(1)})^{2} + \sum_{i = 1}^{n} {\hat{γ}}_{i}^{(1)} (y_{i} - x_{i}^{'} β) I (| y_{i} - x_{i}^{'} β | < λ_{γ}) .

(2.4)

Note that the second term of right hand side in (2.4) will be zero when

| y_{i} - x_{i}^{'} β | \geq λ_{γ}

. When

| y_{i} - x_{i}^{'} β | < λ_{γ}

, we have

{\hat{γ}}_{i}^{(1)} = 0

by (2.3). Thus,

L (β, {\hat{γ}}^{(1)})

in (2.4) is rewritten as

\frac{1}{2} \sum_{i = 1}^{n} (y_{i} - x_{i}^{'} β - {\hat{γ}}_{i}^{(1)})^{2}

. We can easily get its minimizer

{\hat{β}}^{(1)} = (X^{'} X)^{- 1} X^{'} (Y - {\hat{γ}}^{(1)})

. Again, note that this is a convex minimization with given estimates of

γ

. We call this estimating method ‘outlier shifting least squares’. The tuning parameter

λ_{γ}

needs to be pre-specified, and we defer its estimation to Section 2.2.

The convergency of the iterative algorithm for convex function minimization is widely studied by Rockafellar (1997). Specifically, the convergency with case-specific parameters is discussed in Lee et al. (2012). Proposition 3 in Lee et al. (2007) can be directly applied to show the convergence of Algorithm 1 to unique minimizer ( $β^{*}, γ^{*}$ ) under the objective function in (2.1).

Next, we summarize the algorithm that can provide more general application beyond the discussed squared error loss function. Algorithm 1:

Initialize ${\hat{γ}}^{(0)} = 0$ , and obtain ${\hat{β}}^{(0)} = arg min_{β} L (β, 0)$ , and $Y^{(0)}$ .

Update ${\hat{γ}}^{(m + 1)}$ by finding the minimizer of $L ({\hat{β}}^{(m)}, γ)$ with $Y^{(m)}$ .

Update $Y^{(m + 1)}$ with $Y^{(m)} - {\hat{γ}}^{(m + 1)}$ .

Update ${\hat{β}}^{(m + 1)}$ by finding the minimizer of $L (β, γ^{(m + 1)})$ with $Y^{(m + 1)}$ .

Iterate between Step 2 and Step 4 until $| | β^{(m + 1)} - β^{(m)} | |^{2}$ is small.

But, there are certain restrictions of the application such as binomial loss where Y is not quantitative. In this case, it is hard to understand and/or interpret the adjustment in Step 3. We set $Y^{(0)}$ to be the original observations of the response variable. Note that the estimates of $γ$ are used to modify the original observations in Step 3. Thus, when we estimate $β$ in Step 4, the df for model is not affected by the number of non-zero estimate of $γ$ , but the outlying observations are modified only. Once $y_{i}$ is judged as an outlier, then the selected $y_{i} s$ are modified as $y_{i}^{'} = y_{i} - {\hat{γ}}_{i} = y_{i} - r_{i}$ . As we subtract the residual from the original observation, it shifts the outlier onto the fitted model. If it is not selected, that is, if ${\hat{γ}}_{i} = 0$ , then there is no adjustment in $y_{i}$ . As we employ an iterative algorithm which finds the minimizer of $β$ and $γ$ alternately, this adjustment can happen multiple times as shown in Figure 1. This reflects certain advantage of introducing the case-specific parameters to control outliers rather than merely removing them.

A toy example in Figure 1 illustrates how the fitted line and $y_{i}^{(m)}$ are changed as we iterate the estimation procedure described in Algorithm 1.

In Figure 1, only $y_{1}$ is detected as an outlier, and is adjusted as $y_{1} - {\hat{γ}}_{1}^{(1)} - {\hat{γ}}_{1}^{(2)}$ after two iterations of Algorithm 1. When the algorithm is converged after two iterations, the final fitted line (solid line) is very close to the underlying truth (red solid line). We also show the fitted line from OLS after removing $y_{1}$ (dashed line). In this simple example, removing an observation improves model fit at the sacrifice of one degrees of freedom. In contrast, the proposed method illustrates how the huge outlier can be adjusted effectively without being disregarded. In Figure 1, we observe that the obtained ${\hat{γ}}_{1}^{(1)}$ and ${\hat{γ}}_{1}^{(2)}$ have the effect of shifting the potential outlier $y_{1}$ towards the fitted line. Thus, we call the suggested penalty in (2.1) outlier-shifting penalty.

Figure 1

A toy example for the illustration of Algorithm 1 with 10 observations where the first observation $y_{1}$ is highly corrupted. ‘1st fit’ is the initial fit with unadjusted observations, and ‘2nd fit’ is the fit after one iteration. The ‘final fit’ is from the estimates when the Algorithm 1 is converged after two iterations. ‘Deleted’ stands for the fitted line by OLS after removing $y_{1}$ . The thick solid line (in red) represents the underlying truth.

Remark 1. When we replace the penalty for $γ_{i}$ in (2.1) with $λ_{γ} \sum_{i = 1}^{n} | γ_{i} |$ , it will result in the estimator in Lee et al. (2012). In this case, we have ${\hat{γ}}_{i} = sgn (r_{i}) (| r_{i} | - λ_{γ})_{+}$ . The idea of Lee et al. (2012) is good in the sense that outliers are detected and diminished in an automatic fashion. However, large residuals are (even if they are all large) in various sizes, and the universal shrinkage by $λ_{γ}$ does not seem to be the most effective way. Extremely large outliers often have high impact on the fitted model after being reduced by $λ_{γ}$ . Using extremely large value of $λ_{γ}$ could be a solution, but then, we cannot treat moderately large outliers. This is the main motivation of the suggested penalty in this article.

2.2 A default value of

λ_{γ}

Because we target detecting highly contaminated data by non-zero $γ_{i}$ s, it is reasonable to assume that $ε_{i}$ s are not contaminated and normally distributed. Thus, we propose to select $λ_{γ}$ such that

P (| ε_{i} | \geq λ_{γ}) \leq n_{o} / n,

(2.5)

where

n_{o}

is the size of highly contaminated

y_{i} s

. As we are interested in highly corrupted data, we assume the proportion of outliers do not disappear. That is,

n_{o} / n > 0

n \to \infty

Under the assumption of $ε_{i} \sim N (0, σ^{2})$ , the inequality (2.5) is equivalent to

λ_{γ} \geq σ Ψ^{- 1} (\frac{2 n - n_{o}}{2 n}),

where

Ψ^{- 1} (\cdot)

is the inverse cumulative distribution function (C.D.F) of the standard normal distribution. Herein,

σ

should be estimated and so does the number of outliers

n_{o}

. Since

λ_{γ}

is required to be specified prior to the model estimation procedure, it is desirable to use a robust estimate for

σ

, such as interquartile range, median absolute deviation (MAD), or robust scale estimators suggested by Rousseeuw and Croux (1993). Investigation of the residuals from OLS often suggests reasonable number of outliers

n_{o}

. Then we set

λ_{γ}

to be the following value,

{\hat{λ}}_{γ} = \hat{σ} Ψ^{- 1} (\frac{2 n - {\hat{n}}_{0}}{2 n}) .

(2.6)

The above default value, which is devised to catch highly contaminated observations, differs from a conventional rule of defining outliers as a fixed percentage (typically, 1% or 5%) of the observations. Instead,

{\hat{λ}}_{γ}

is adaptive to the data through the proportion of outliers

n_{o} / n

2.3 Asymptotic properties

We investigate the asymptotic properties of ${\hat{β}}^{(m)}$ under the outlier-shifting least squares and show its limiting distribution. Note that the classical linear model in (1.1) has iid errors, which is far from the reality when outliers exist. To incorporate the outliers in the error term, errors should be a mixture of iid errors and outliers. For this purpose, we may use an error of $ε_{i}^{*} = (1 - α) ε_{i} + α ξ_{i}$ where $ε_{i}$ is iid errors with mean 0 and finite variance, and $ξ_{i}$ s are independent and follow some heavy-tailed distribution with finite mean and variance. Finally, $α (> 0)$ is the proportion of outliers. Note that $ε_{i}^{*}$ s are not necessarily identical or symmetric. Then, under the model of

y_{i} = x_{i}^{'} β + ε_{i}^{*},

(2.7)

we have

E ({\hat{β}}^{OLS}) = β + (X^{'} X)^{- 1} E (X^{'} ε^{*})

. For

{\hat{β}}^{OLS}

to be consistent,

E (X^{'} ε^{*})

needs to approach to zero with probability 1. However, this is not the case when the distribution of

ξ_{i}

does not have mean zero. So, we impose a restriction on

ε_{i}^{*}

that the mean of

ε_{i}^{*}

tends to zero with probability 1. This implies that asymptotically equal amount of outliers exist in each tail. Then, under some mild conditions on X,

{\hat{β}}^{OLS}

is consistent by Linderberg–Feller type of central limit theorem under independent but not identical errors. Here is the list of conditions we consider. \\ CONDITIONS

$lim_{n \to \infty} \frac{1}{n} \sum_{i = 1}^{n} x_{i} x_{i}^{'} = Q$ , where $Q$ is a non-negative definite matrix.

$E (x_{i}) < \infty$ for all i.

$λ_{γ} = O (1)$

$E (ξ_{i}) \to 0$ w.p. 1.

Conditions 1 and 2 are the standard moment conditions for X. Condition 3 arises from the assumption that the proportion of outliers $n_{o} / n$ approaches to a non-zero constant discussed in Section 2.2. Condition 4 implies the amount of outliers in each tail is asymptotically equal. In practice, the proposed method works well when all the outliers are in one tail as they will be effectively shifted. For the asymptotic results, we need further clarification. When the number of outliers are ignorable ( $n_{o} / n \to 0$ ), then we do not need the Condition 4. Since $n_{o} = \sqrt{n}$ is also ignorable, we may not need Condition 4 in many cases except for the case of $n_{o} / n \to c$ , where $c > 0$ . Then, the proposed estimator is $\sqrt{n}$ -consistent as shown in Theorem 1.

Theorem Assuming Conditions 1–4 hold, we have, for any m,

\sqrt{n} ({\hat{β}}^{(m)} - β) \overset{d}{\to} N (0, σ^{2} Q^{- 1}) .

Proof. See Appendix.

3 Other applications of shifting outliers

3.1 Outlier-shifting in robust loss function

A well-known robust approach for improving the least squares method is Huber's loss function:

ρ^{H} (u) = \frac{u^{2}}{2} I (| u | \leq c) + c (| u | - \frac{c}{2}) I (| u | > c) .

(3.1)

Even for this robust loss function, we observe from our empirical investigations that it is still possible to improve the performance in the presence of massively contaminated observations. Therefore, we add the case-specific parameters with the outlier-shifting penalty to Huber's loss function, resulting in the following outlier-shifting Huber's loss function:

L (β, γ) = \sum_{i = 1}^{n} ρ^{H} (y_{i} - x_{i}^{'} β - γ_{i}) + \sum_{i = 1}^{n} γ_{i}^{2} I (| y_{i} - x_{i}^{'} β | < λ_{γ}) .

Algorithm 1 can be implemented to estimate the parameters. For the threshold value c in (3.1), we set it scale invariant as

c = c_{0} σ

, and also use a robust scale estimator for

σ

as discussed in Section 2.2. We conduct extensive numeric comparisons among the estimates from ordinary least squares, outlier-shifting least squares, Huber's loss and outlier-shifting Huber's loss in Section 4.

The Huber's approach shares the idea of M-quantile (Breckling and Chambers 1988) for robust regression. The two methods have the same form of loss function when M-quantile targets conditional median. Thus, we naturally see that M-quantile regression might get benefit by adopting the proposed outlier-shifting penalty for the quantile regression. We may need to adjust the form of the proposed penalty due to the asymmetric loss function in quantile regression except for the median. The details of adopting the proposed penalty in the light of quantile regression will not be considered here; here is the sketch of the idea. Under the linear model of $y_{i} = x_{i}^{'} β + ε$ , quantile regression estimator for $q$ th quantile is defined as the minimizer of the check loss function:

L_{q} (β) = \sum_{i = 1}^{n} (y_{i} - x_{i}^{'} β) {q - 1 + I (y_{i} - x_{i}^{'} β \geq 0)} .

Then, we can incorporate an asymmetric outlier shifting penalty with case-specific parameter

γ_{i}

\begin{matrix} L_{q} (β, γ) & = \sum_{i = 1}^{n} (y_{i} - x_{i}^{'} β - γ_{i}) {q - 1 + I (y_{i} - x_{i}^{'} β - γ_{i} \geq 0)} \\ + \sum_{i = 1}^{n} γ_{i} q I (0 < y_{i} - x_{i}^{'} β < λ_{γ}) + \sum_{i = 1}^{n} γ_{i} (1 - q) I (- λ_{γ} < y_{i} - x_{i}^{'} β < 0), \end{matrix}

which may result in robustified quantile regression.

3.2 Outlier-shifting in Lasso-type models

In this section, we switch to the often encountered high-dimensional problem with a large number of covariates. To address this problem, the lasso-type penalties are widely used. The idea of case-specific parameters and shifting outlier penalty is directly applicable to lasso (Tibshirani 1996) and Huberized lasso (Rosset and Zhu 2007) to improve robustness. Lasso is defined as finding a minimizer of

L (β) = \frac{1}{2} \sum_{i = 1}^{n} (y_{i} - x_{i}^{'} β)^{2} + λ_{β} \sum_{j = 1}^{p} | β_{j} | .

(3.2)

We denote

{\hat{β}}^{lasso}

as the minimizer of (3.2). By introducing case-specific parameters

γ

and the outlier-shifting penalty used in (2.1), the outlier-shifting lasso procedure is defined as finding the minimizer of

L (β, γ) = \frac{1}{2} \sum_{i = 1}^{n} (y_{i} - x_{i}^{'} β - γ_{i})^{2} + λ_{β} \sum_{j = 1}^{p} | β_{j} | + \sum_{i = 1}^{n} γ_{i}^{2} I (| y_{i} - x_{i}^{'} β | < λ_{γ}) .

(3.3)

To estimate $β$ and $γ$ , Algorithm 1 in Section 2.1 can be directly implemented. However, the L1-type penalty on $β$ is non-differentiable. Thus, we resort to coordinate the decent type of methods (Wu and Lange 2008).

With initial values of ${\hat{γ}}^{(0)} = 0$ and ${\hat{β}}^{(0)} = {\hat{β}}^{lasso}$ , we proceed to iterate updating $\hat{γ}$ and $\hat{β}$ alternately as in Algorithm 2. At the $(m + 1)$ th iteration, the outlier-shifting lasso estimates are

\begin{matrix} {\hat{γ}}_{i}^{(m + 1)} & = & (y_{i} - x_{i}^{'} {\hat{β}}^{(m)}) I (| y_{i} - x_{i}^{'} {\hat{β}}^{(m)} | \geq λ_{γ}), for i = 1, \dots, n, \\ {\hat{β}}^{(m + 1)} & = & arg min_{β \in^{p}} {\frac{1}{2} \sum_{i = 1}^{n} (y_{i} - x_{i}^{'} β - {\hat{γ}}_{i}^{(m + 1)})^{2} + λ_{β} \sum_{j = 1}^{p} | β_{j} |} . \end{matrix}

We use the same estimated value of

λ_{γ}

provided in Section 2.2.

Rosset and Zhu (2007) considered Huberized lasso where the squared error loss in (3.2) is replaced by Huber's loss function. In fact, Huberized lasso is coincidentally the same as what proposed by Lee et al. (2012) with case-specific parameters. Analogous to Section 3.1, we can add the outlier-shifting penalty to Huberized lasso to intensify the robustness. Then, the outlier-shifting Huberized lasso is defined as the minimizer of

L (β, γ) = \frac{1}{2} \sum_{i = 1}^{n} ρ^{H} (y_{i} - x_{i}^{'} β - γ_{i}) + λ_{β} \sum_{j = 1}^{p} | β_{j} | + \sum_{i = 1}^{n} γ_{i}^{2} I (| y_{i} - x_{i}^{'} β | < λ_{γ}),

(3.4)

where

ρ^{H} (\cdot)

is given in (3.1). The performance of lasso outlier-shifting lasso, Huberized lasso and outlier-shifting Huberized lasso is investigated through the simulation studies in the next section.\\ Remark. Although we choose lasso penalty for variable selection, the choice of penalty for

β

is not a crucial part. It is obvious that other penalties for

β

can be combined with outlier-shifting penalty. For example, smoothly clipped absolute deviation (SCAD) penalty proposed by Fan and Li (2001) can replace lasso penalty.

3.3 Connection to maximum likelihood approach

Now, treating the negative log likelihood as a loss function, we consider a specific form of Algorithm 1 with normal likelihood. Because $γ_{i}$ s are employed to account for highly corrupted outliers, we may assume normally distributed $ε_{i}$ s in Model (1.2). Since $γ_{i}$ s explain extreme outliers, the normality assumption is reasonable. And by the constraint given in Section 1, $ε_{i}$ s are zero when $γ_{i}$ s are non-zero. Then, the likelihood is simply

\prod_{i = 1}^{n} \frac{1}{\sqrt{2 π σ^{2}}} exp \{- \frac{(y_{i} - x_{i}^{'} β - γ_{i})^{2}}{2 σ^{2}}\} .

The maximum likelihood estimate of

β

σ^{2}

and

γ

can be written as

\begin{matrix} \hat{β} & = & (X^{'} X)^{- 1} X^{'} (Y - \hat{γ}), \\ {\hat{σ}}^{2} & = & | | Y - X \hat{β} - \hat{γ}) | |^{2} / n, \\ \hat{γ} & = & Y - X \hat{β} . \end{matrix}

Note that

\hat{γ} = (γ_{1}, \dots, γ_{n})^{'}

are in fact the residuals, which take non-zero values for both the outliers and non-outlying observations. To be coherent with the description of

γ_{i}

in Section 1,

{\hat{γ}}_{i}

is modified to be

{\hat{γ}}_{i} = (y_{i} - x_{i}^{'} \hat{β}) I (| y_{i} - x_{i}^{'} \hat{β} | \geq λ_{γ})

, for

i = 1, . . ., n

, where

λ_{γ}

is a pre-specified threshold. With this modified

\hat{γ}

, we define an updated response variable

y_{i}^{*} = y_{i} - {\hat{γ}}_{i}

. It implies that the observations of

y = (y_{1}, \dots, y_{n})^{'}

are treated as outliers and accordingly adjusted only if

| y_{i} - x_{i}^{'} \hat{β} |

is larger than

λ_{γ}

. Or equivalently,

y_{i}^{*} = x_{i}^{'} \hat{β}

| r_{i} | \geq λ_{γ}

and

y_{i}^{*} = y_{i}

, otherwise. This adjustment is to shift the outlying observations onto the fitted values while preserving the other observations. Treating the modified

y_{i}^{*}

as a new observation, we repeat the estimation procedure to update

\hat{β}

{\hat{σ}}^{2}

and

\hat{γ}

. Here is the detailed algorithm of the estimation procedure. Algorithm 2:

Initialize $β^{(0)}$ and $σ^{2 (0)}$ with their least square estimates, set ${\hat{γ}}^{(0)} = 0$ , and use the original observations of $y$ as $Y^{(0)}$ .

Update ${\hat{γ}}_{i}^{(m + 1)}$ with $(y_{i}^{(m)} - x_{i}^{'} {\hat{β}}^{(m)}) I (| y_{i}^{(m)} - x_{i}^{'} {\hat{β}}^{(m)} | \geq λ_{γ})$ , for $i = 1, \dots, n$ .

Update $Y^{(m + 1)}$ with $(Y^{(m)} - {\hat{γ}}^{(m + 1)})$ , where ${\hat{γ}}^{(m + 1)} = ({\hat{γ}}_{1}^{(m + 1)}, \dots, {\hat{γ}}_{n}^{(m + 1)})^{'}$ .

Update ${\hat{β}}^{(m + 1)}$ with $(X^{'} X)^{- 1} X^{'} Y^{(m + 1)}$ .

Update ${\hat{σ}}^{2 (m + 1)}$ with $| | Y^{(m + 1)} - X {\hat{β}}^{(m + 1)} | |^{2} / n$ .

Iterate Steps 2 to 5 until $| | β^{(m + 1)} - β^{(m)} | |^{2}$ is small.

Remark 2. Notice that the obtained $\hat{γ}$ is identical to (2.3), which is estimated with the outlier-shifting penalty described in the previous section. Therefore, this provides an intuitive interpretation of the outlier-shifting penalty as well as an alternative algorithm to optimize the objective function (2.1).

Interestingly, this outlier adjustment method also shares some similarity with the Expectation–Maximization (EM) algorithm (Dempster et al. 1977) in terms of shrinkage towards the fitted value in an iterative manner. A clear distinction between the two is that the outlier adjustment method automatically determines the unidentified outliers via the tuning parameter $λ_{γ}$ , while the missing observations are identified in advance in the EM algorithm.

4 Simulation studies

We conduct two sets of simulations, investigating various sizes and proportions of contaminated observations. The second set is in the supplementary materials. In the simulations, we compare the suggested estimators to existing methods. They are (i) ordinary least squares (OLS); (ii) outlier-shifting least squares ( ${OLS}^{S}$ ); (iii) Huber's estimator ( $H$ ); (iv) outlier-shifting Huber's estimator ( $H^{S}$ ); (v) median regression ( $med$ ); (vi) $lasso$ ; (vii) outlier-shifting lasso ( ${lasso}^{S}$ ); (viii) Huberized lasso ( $Hlasso$ ); and (ix) outlier-shifting Huberized lasso ( ${Hlasso}^{S}$ ). The proposed outlier-shifting penalty based approaches are notated with superscript ‘ $S$ ’.

For ${OLS}^{S}$ , we take the pre-specified value in (2.6) as $λ_{γ}$ , wherein MAD/0.6745 of residuals from fitting the median regression is used as the robust estimator of $σ$ . The MAD should be scaled to be a scale estimator. Thus, we divided it by 0.6745, which is known to be a robust estimator of $σ$ when the errors do not follow normal distribution (Huber and Ronchetti 2009). For Huber's function defined in (3.1), we report the result at $c = 1.5$ that is observed to have the best empirical performance among various examined values of c. To implement $H$ , we use rlm function in MASS package. For the implementation of $lasso$ , glmnet function in glmnet package is utilized. For the selection of penalty parameter in $lasso$ , we use 10-fold cross-validation with 100-candidate values of the penalty parameter, $λ_{β}$ .

In all the simulation studies, we generate $K = 500$ data sets with sample size $n = 100$ from the model,

Y = X β + σ ε,

where the components of X and

ε

follow standard normal distribution. We use

σ = 1

for the base model. We consider contaminating 10%, 20% and 30% of the response in the base model with

σ = 3, 6

and 10, respectively, which produces nine different combinations of the scenarios. Thus, the generated errors follow mixture normal distribution of

(1 - π) N (0, 1) + π N (0, σ)

, where

π = 0.1

, 0.2 and 0.3. For each case, we report the mean squared error (

MSE

) between the true and the fitted regression lines over

K

simulated data set to assess the performance of each estimator,

MSE = \sum_{k = 1}^{K} ({\hat{β}}^{k} - β)^{'} E (X^{'} X) ({\hat{β}}^{k} - β) / K

, where

{\hat{β}}^{k}

is the parameter estimate from the

k

th simulated data set.

Simulation set 1: We consider $p = 8$ covariates and set $β = (3, 1.5, 0, 0, 2, 0, 0, 0)^{'}$ and take the correlation between $x_{i}$ and $x_{j}$ to be $ϱ^{| i - j |}$ with $ϱ = 0.5$ , which are also used in Tibshirani (1996). (The results based on orthogonal covariates with $ϱ = 0$ is in supplementary material.) The $MSE$ values obtained by the nine estimators are provided in Table 1.

Table 1

Point estimate of MSE, and its standard error (in parentheses) multiplied by 1000 based on 500 simulated data sets obtained by OLS, ${OLS}^{S}$ , $H$ , $H^{S}$ , $med$ , $lasso$ , ${lasso}^{S}$ , $Hlasso$ and ${Hlasso}^{S}$ under $σ$ = 3, 6 and 10 with the contaminated proportions of 10%, 20% and 30%.

,	10%	20%	30%
		$σ = 3$
OLS	161.6 (4.3)	234.7 (6.1)	312.0 (8.1)
${OLS}^{S}$	122.8 (3.5)	165.1 (4.7)	220.8 (6.4)
$H$	121.3 (3.0)	159.5 (4.1)	212.2 (5.7)
$H^{S}$	126.3 (3.2)	156.0 (4.2)	197.5 (5.3)
$med$	168.1 (4.1)	200.3 (4.9)	243.2 (6.3)
$lasso$	131.1 (4.2)	191.9 (5.9)	250.3 (7.7)
${lasso}^{S}$	103.1 (3.2)	136.2 (4.2)	177.3 (5.7)
$Hlasso$	101.7 (3.2)	136.9 (4.3)	181.7 (5.8)
${Hlasso}^{S}$	102.3 (3.1)	129.3 (4.0)	167.8 (5.3)
		$σ = 6$
OLS	393.7 (11.6)	712.3 (19.2)	1,045 (27.0)
${OLS}^{S}$	198.5 (6.3)	341.7 (10.8)	536.8 (16.5)
$H$	134.4 (3.4)	205.9 (5.5)	340.5 (9.9)
$H^{S}$	128.3 (3.2)	159.4 (4.0)	223.3 (6.2)
$med$	175.0 (4.3)	221.5 (5.5)	296.4 (7.9)
$lasso$	320.5 (11.3)	585.8 (19.1)	845.7 (26.2)
${lasso}^{S}$	127.8 (4.1)	194.6 (6.8)	288.5 (10.4)
$Hlasso$	127.6 (4.1)	218.0 (7.1)	358.6 (12.2)
${Hlasso}^{S}$	107.6 (3.3)	143.7 (4.5)	212.6 (7.6)
		$σ = 10$
OLS	943.7 (29.2)	1,843 (51.0)	2,776 (71.4)
${OLS}^{S}$	394.2 (13.8)	776.7 (26.1)	1,286 (40.5)
$H$	141.3 (3.6)	234.7 (6.3)	450.7 (13.7)
$H^{S}$	128.1 (3.2)	162.3 (4.1)	248.3 (7.3)
$med$	178.2 (4.5)	230.6 (5.8)	321.6 (8.8)
$lasso$	774.8 (28.7)	1,532 (51.9)	2,276 (70.0)
${lasso}^{S}$	204.5 (8.8)	356.4 (15.3)	553.3 (21.7)
$Hlasso$	171.9 (5.8)	368.3 (13.7)	724.3 (26.3)
${Hlasso}^{S}$	111.1 (3.6)	161.5 (5.5)	271.4 (10.3)

As expected, we see that OLS is the most vulnerable to outliers. MSE values obtained by the methods of ${OLS}^{S}$ , $H$ and $H^{S}$ are similar when 10% of the data is contaminated with $σ = 3$ . As the coverage and the magnitude of contamination increases, the advantage of the proposed $H^{S}$ and ${Hlasso}^{S}$ over the other methods becomes much more prominent; it is not surprising that $H$ performs well. For the base model without outliers, OLS, ${OLS}^{S}$ , $H$ and $H^{S}$ yield MSE values at the similar levels of 0.0927, 0.0924, 0.0937 and 0.1012, respectively. In summary, $H^{S}$ and ${Hlasso}^{S}$ , which are doubly robustified by Huber's loss function and the proposed outlier-shifting penalty, manifest the most benefit for a severely contaminated data set. We also consider a non-parametric robust approach—median regression (Koenker 2005)—for comparison in this simulation study. It pursues the conditional median regression model rather than the mean models considered thus far. However, the mean and median are equivalent in the considered scenario of symmetric error distributions. The median regression approach is always inferior to $H^{S}$ and superior to ${lasso}^{S}$ , and fall between ${OLS}^{S}$ and $H$ most of the time. Because the median is known to be less efficient than the mean under normal distributions, we can expect that the median regression obtains larger variances of parameter estimates than the mean regression. The phenomenon is demonstrated by Figure 2, which shows the distributions of the regression coefficient estimates obtained by the investigated five methods. We see that the estimates of the median regression approach have slightly larger variance than that of $H^{S}$ . Overall, the difference in the variances of the five estimators enlarges as the level of contamination increases. And $H^{S}$ produces the smallest variance among the five methods as shown in Figure 2.

Figure 2

Distribution of estimates from 500 simulated samples under the simulation set 1. See the text for the details of simulation setting. Horizontal lines represent the true regression coefficient.

As the four regression parameter penalized methods ( $lasso$ , ${lasso}^{S}$ , $Hlasso$ and ${Hlasso}^{S}$ ) can produce sparse solutions, we examine the false negative rates and false positive rates. Note that there are three non-zero coefficients in true $β$ . We record the average probability of identifying the non-zero coefficient from 500 simulated data sets in Table 2. We also present the average false positive rate of the five truly zero coefficients in the parentheses in the same table. As the standard errors of the average false positive rates are all around 0.01, they are not presented.

From Table 2, the four methods select more variables than they should as the false positive rates are considerably high. As a consequence, the true positive rates are all close to 1. The false discovery rates from the suggested methods are higher than their competitors.

Table 2

Average false negative rate of three non-zero coefficients and average false positive rate of five zero coefficients (in parenthesis) based on 500 simulated data sets obtained by $lasso$ , ${lasso}^{S}$ , $Hlasso$ and ${Hlasso}^{S}$ with $σ$ = 6 and 10 with contaminated proportions of 10%, 20% and 30%

,		$σ = 6$			$σ = 10$
	10%	20%	30%	10%	20%	30%
$lasso$	0 (.445)	0 (.443)	0 (.466)	0.001 (.447)	0.005 (.435)	0.024 (.457)
${lasso}^{S}$	0 (.519)	0 (.548)	0 (.581)	0 (.570)	0 (.572)	0 (.583)
$Hlasso$	0 (.471)	0 (.501)	0 (.508)	0 (.517)	0 (.522)	0 (.531)
${Hlasso}^{S}$	0 (.492)	0 (.548)	0 (.579)	0 (.517)	0 (.568)	0 (.626)

Finally, we examine the performance of ${\hat{γ}}_{i}$ s. Note that we generate iid errors from $N (0, 1)$ first. Then, 10% (or 20% or 30%) of the errors (out of n) are multiplied by 3 (6 or 10). Thus, the generated errors follow mixture normal distribution of $(1 - π) N (0, 1) + π N (0, σ)$ , where $π = 0.1$ , 0.2 and 0.3, and $σ = 3$ , 6 and 10 are considered. It is very possible that some errors are still quite small after being multiplied by 10. Thus, we regard $ε_{i}$ s that satisfy $| ε_{i} | > 2.5$ as outliers. As the sample size n = 100 and $σ = 1$ in the data set without contamination, we would expect to observe about one $ε_{i}$ which satisfies the above criterion under the assumed normality. Then, we count the number of non-zero ${\hat{γ}}_{i}$ s that capture the outliers. As the contaminated proportion increases, we expect larger number of outliers are generated. Thus, we report the average percentage of correct detection of the outliers in Table 3. As standard errors are all small (around 0.01, or smaller), we omit them in the table. In general, as the contaminated proportion decreases, and as the contaminated magnitude increases, we see a higher rate of identifying the outliers.

Table 3

Averaged probability of identifying the outliers with $| ε_{i} | > 2.5$ based on 500 simulated data sets obtained by ${OLS}^{S}$ , $H^{S}$ , ${lasso}^{S}$ and ${Hlasso}^{S}$ with contaminated proportions of 10%, 20% and 30% and magnitude of $σ$ = 3, 6 and 10

,		$σ = 3$			$σ = 6$			$σ = 10$
	10%	20%	30%	10%	20%	30%	10%	20%	30%
${OLS}^{S}$	0.82	0.78	0.70	0.88	0.83	0.73	0.87	0.81	0.71
$H^{S}$	0.88	0.85	0.78	0.94	0.92	0.86	0.96	0.94	0.89
${lasso}^{S}$	0.87	0.85	0.77	0.92	0.88	0.81	0.91	0.87	0.79
${Hlasso}^{S}$	0.59	0.56	0.51	0.74	0.73	0.69	0.88	0.88	0.82

5 Applications

5.1 Analysis of stack loss data

The stack loss data set (Brownlee 1965) has been examined by many researchers (Daniel and Wood 1980; Chambers and Heathcote 1981; Carroll and Ruppert 1985). The details of pre-analysis can be found in Chapter 5.4 of Daniel and Wood (1971), and using non-linear covariates are included in Denby and Mallows (1977).

The data set has 21 observations with 3 explanatory variables. The data describe the operation of a plant. Researchers have sought to explain the percentage of unconverted ammonia that escapes from the plant (Y) by the flow of cooling air ( $X_{1}$ ), inlet temperature of cooling water ( $X_{2}$ ) and concentration of acid ( $X_{3}$ ). Thus, we consider the following standard linear model:

y_{i} = β_{0} + \sum_{j = 1}^{3} β_{j} x_{ij} + ε_{i}, i = 1, \dots, n,

where

ε_{i}

s are iid random variables with mean 0 and variance

σ^{2}

. Most people concluded that the Observations 1, 3, 4 and 21 are outliers. Chambers and Heathcote (1981) and Hoeting et al. (1996) suggested deleting the two observations of 4 and 21 can improve the fitted surface. Thus, we compose three different data sets of containing (a) all 21 observations, (b) 19 observations (Observations 4 and 21 are deleted), and (c) 17 observations (Observations 1, 3, 4 and 21 are deleted), which were examined by Chambers and Heathcote (1981). We fit

OLSs

model to these three data sets, then we record the change in the estimates of the regression coefficients. To estimate the scale parameter, the median absolute deviation of the residuals from fitting a median regression is utilized. For comparison, the estimates from OLS,

H

and Tukey's robust regression (

Tukey

) by (Beaton and Tukey 1974) are examined. The loss function for Tukey's robust regression is defined as,

ρ^{T} (u) = \frac{1}{6} [1 - {1 - {(u / (6 c))}^{2}}^{3}] I (| u | \leq 6 c) .

(5.1)

The values used for the bending constant c in

ρ^{H} (u)

and

ρ^{T} (u)

are 1.345 and 4.685, respectively. Huber and Ronchetti (2009) argued that c = 1.345 in

ρ^{H} (c)

is a good choice, and showed that asymptotically, it is 95% as efficient as least squares if the true distribution is normal. Similarly, c = 4.645 in

ρ^{T} (c)

provides 95% of asymptotic efficiency under the normal distribution. rlm function in R package MASS is used for implementing

H

and Tukey's robust regression. The estimates from the above four methods are represented in Table 5.1.

Table 4

Estimates of $β_{0}$ , $β_{1}$ , $β_{2}$ , $β_{3}$ for the stack loss data. Data 1, data 2 and data 3 stand for the data set with 21, 19 and 17 observations, respectively

,	${\hat{β}}_{0}$	${\hat{β}}_{1}$	${\hat{β}}_{2}$	${\hat{β}}_{3}$
		Data 1
OLS	17.524	6.561	4.094	$-$ 0.815
$H$	17.596	7.604	2.927	$-$ 0.685
$Tukey$	17.800	8.504	2.057	$-$ 0.602
${OLS}^{S}$	17.112	7.614	1.781	$-$ 0.387
		Data 2
OLS	17.688	8.770	1.756	$-$ 0.583
$H$	17.593	8.600	1.702	$-$ 0.581
$Tukey$	17.625	8.673	1.703	$-$ 0.603
${OLS}^{S}$	17.228	7.914	1.431	$-$ 0.517
		Data 3
OLS	16.943	7.313	1.825	$-$ 0.359
$H$	17.006	7.503	1.643	$-$ 0.397
$Tukey$	17.018	7.508	1.636	$-$ 0.390
${OLS}^{S}$	17.131	7.722	1.441	$-$ 0.484

To measure the sensitivity of the estimates as we remove two and four observations, the absolute differences between the minimum and the maximum estimates among the three data sets are calculated. For example, sensitivity of OLS for $β_{1}$ is $∣ 6.561 - 8.770 ∣$ = 2.209. Therefore, estimator with large value of sensitivity indicates that the estimator is sensitive to outliers. Three variables are standardized to have scale invariant sensitivity. The sensitivity of the estimates from the five methods are demonstrated in Table 5.

Table 5

Sensitivity of five estimators in Table 5.1 at $β_{1}$ , $β_{2}$ , $β_{3}$ with the stack loss data

,	${\hat{β}}_{1}$	${\hat{β}}_{2}$	${\hat{β}}_{3}$
OLS	2.2092	2.3381	0.4558
$H$	1.0967	1.2842	0.2878
$Tukey$	1.1642	0.4205	0.2125
${OLS}^{S}$	0.3005	0.3498	0.1293

Among the four aforementioned methods, the proposed method is the most stable; That is, the amount of change in the estimates is the smallest for the suggested method as we change the sample size by removing two and four outliers.

Figure 3 shows the residual plots from the four methods of OLS, $H$ , $Tukey$ and ${OLS}^{S}$ with all observations. In the lower right panel, the residuals from the Observations 1, 3, 4 and 21 (marked with *) are very close to zero due to shifting. This partly reflects that our method reduces the undesirable impact effectively. Overall, the residuals from ${OLS}^{S}$ are smaller than the other methods. Another interesting point is that the residual from Observation 2 is very close to zero with our method, even if it is not shifted. Some researchers (Andrews and Pregibon 1978; Dempster and Gasko-Green 1981; Rousseeuw and Leroy 2005) reported that Observation 2 is a moderate outlier, although it is not clearly judged by Figure 3. Our conjecture is that when the injurious effects of Observations 1, 3, 4 and 21 are significantly reduced by our method, Observation 2 is no more an outlier. Especially, the fit of Observation 21 is significantly improved by the suggested method as shown in Figure 3, whereas the other methods still suffer from the large residual from Observation 21.

Figure 3

Residual plots from the four methods of OLS, $H$ , $Tukey$ and ${OLS}^{S}$ with all observations. In the ${OLS}^{S}$ , the four observations with $*$ are detected as outliers (observations 1, 3, 4 and 21) and shifted by our method.

5.2 Analysis of mid-Atlantic wage data

For analysis of real data, we numerically compare the eight models (OLS, ${OLS}^{S}$ , $H$ , $H^{S}$ , $lasso$ , ${lasso}^{S}$ , $Hlasso$ and ${Hlasso}^{S}$ ) employed in simulation studies in Section 4. The mid-Atlantic wage data is manually assembled by Steve Miller of Open BI (www.openbi.com) from the March 2011 supplement to Current Population Survey and is available in R package ISLR. The data set is from 3,000 males in the mid-Atlantic region of the United States. There are eight explanatory variables of year, age, married, race, edu, class, health, and ins, and one response variable wage. Here are the details of each variable:

wage: Workers raw wage in 1,000 USD

year: 2003 to 2009

age: 18 or above

married: A factor with levels (a) Married and (b) The others (composed of never married, widowed, divorced, separated)

race: A factor with levels (a) White, (b) African American and (c) The others

edu: A factor with levels (a) Less than high school graduate (b) High school graduate, (c) Some college, (d) College graduate and (e) Advanced degree

class: A factor with levels (a) Industrial and (b) Information indicating type of job

health: A factor indicating health status of (a) Good or less than good and (b) Very good

ins: A factor for status of having health insurance. (a) Yes and (b) No

We perform pre-analysis and modify the original variables slightly and do log transformation of the response variable due to heteroscedasticity. For the variable married, we combine the four categories of never married, widowed, divorced and separated into one category since there seems little difference. Similarly, the others category in race includes Asian which was a separate category in the original data in ISLR package. We treat year and age as continuous variables and the others as categorical variables. Since ${age}^{2}$ is a significant variable (details can be found in Chapter 1 of Games et al. 2013), we add it into the model. To incorporate the categorical variables in the regression model, we make indicator variables. For example, we make two indicator variables for race where the first level is considered as a baseline. Thus, we have ${race}_{2 i}$ and ${race}_{3 i}$ , where ${race}_{2 i} = 1$ if ith worker is an African American and 0 otherwise. In each categorical variable, level of 1 is considered as a baseline. Finally, we take log transformation on the response variable. The model we consider is

\begin{matrix} log ({wage}_{i}) & = & β_{0} + β_{1} {year}_{i} + β_{2} {age}_{i} + β_{3} {age}_{i}^{2} + β_{4} {married}_{i} + β_{5} {race}_{2 i} + β_{6} {race}_{3 i} \\ + β_{7} {edu}_{2 i} + β_{8} {edu}_{3 i} + β_{9} {edu}_{4 i} + β_{10} {edu}_{5 i} + β_{11} {class}_{i} + β_{12} {health}_{i} \\ + β_{13} {ins}_{i} + ε_{i}, \end{matrix}

(5.2)

where

ε_{i}

s are independent error with mean 0 and finite variance. To compare the eight candidate models, we perform cross-validation studies. In detail, first, we randomly split the data into two parts where the first part with 2,000 samples is regarded as a training data set, and the other one with 1,000 samples is test dataset. Second, we fit the candidate models and predict the value of response variable in the test data. Then, the absolute deviation between predicted value and observed value is measured. We incorporate a robust measure of absolute deviation since the squared distance is very sensitive to outliers. Finally, cross-validated score (CVS) is calculated by

\sum_{i = 1}^{1000} | y_{i} - y_{i}^{pred} | / 1000

, where

y_{i}

is the observed value of the response variable (

log ({wage}_{i})

) in the test data and

y_{i}^{pred}

is the predicted value of

y_{i}

. To reduce variation raised from random split, we repeat the above procedure for 500 times.

The mean of 500 $CVS$ values is in Table 6. Standard errors of the means are all very close to 0.02, thus omitted in the table. The results of the predictive performance of the eight models illustrate that $lasso$ and ${OLS}^{S}$ perform better than the others. The value of c in Table 6 is the threshold in Huber's loss function in (3.1).

As the improvement by our methods is quite minor, the values of the estimates by the above methods are similar in general. When income is treated as a response variable in economic data, extreme outliers are often detected even if the response variable is log transformed. A robust model may improve the prediction as our analysis shows here.

Table 6

Mean of $CVS$ based on 500 random split of mid-Atlantic wage data by models of OLS, ${OLS}^{S}$ , $H$ , $H^{S}$ , $lasso$ , ${lasso}^{S}$ , $Hlasso$ and ${Hlasso}^{S}$ . All values are multiplied by 100.

,OLS	${OLS}^{S}$	$H (c = 2)$	$H (c = 2.5)$
20.13	19.95	20.08	20.08
$H^{S} (c = 2)$	$H^{S} (c = 2.5)$	$lasso$	${lasso}^{S}$
20.07	20.06	19.95	20.07
$Hlasso (c = 2)$	$Hlasso (c = 2.5)$	${Hlasso}^{S} (c = 2.5)$	${Hlasso}^{S} (c = 2.5)$
20.07	20.06	20.08	20.10

6 Discussion

In this work, we propose a method of automatically adjusting severely outlying observations by using case-specific parameters and shifting outlier penalty to reduce the undesirable impact. We provide a default value for the threshold of the outliers to capture huge outliers and to attain consistency. The outlier shifting trace is demonstrated with the example displayed in Figure 1. In terms of model estimation, we demonstrate the equivalence between the maximum likelihood-based and convex minimization approaches in the case of least squares estimation. The convex minimization approach enables applying the proposed method to a wide range of modelling procedures, such as Huber's regression, lasso and Huberized lasso. Our empirical investigation shows that the advantage of the proposed method in accurately capturing the true regression surface increases with the extent of data contamination. Despite of our focus on linear models in this work, the new method can be tailored to non-linear models involving convex minimization problems. We will extensively investigate more complicated cases in the future work.

Appendix

Proof of Theorem 1 First consider $m = 1$ . For all i, we have,

P (| {\hat{γ}}_{i}^{(1)} | > 0) = P (| y_{i} - x_{i}^{'} {\hat{β}}^{(0)} | > λ_{γ}) .

(6.1)

Since

| y_{i} - x_{i}^{'} {\hat{β}}^{(0)} | = O_{p} (1)

and

λ_{γ} = O (1)

P (| {\hat{γ}}_{i}^{(1)} | > 0)

does not approach to zero as n increases. Thus,

{\hat{γ}}_{i}^{(1)} = O_{p} (1)

. On the other hand,

{\hat{β}}^{(1)} = (X^{'} X)^{- 1} X^{'} (Y^{(0)} - {\hat{γ}}^{(1)}) = {\hat{β}}^{OLS} - (X^{'} X)^{- 1} X^{'} {\hat{γ}}^{(1)}

. Now, a slight modification shows that

(X^{'} X)^{- 1} X^{'} {\hat{γ}}^{(1)} = {(\frac{1}{n} \sum_{i = 1}^{n} x_{i} x_{i}^{'})}^{- 1} \frac{1}{n} \sum_{i = 1}^{n} x_{i} {\hat{γ}}_{i}^{(1)} .

Then, we have

(\frac{1}{n} \sum_{i = 1}^{n} x_{i} x_{i}^{'})^{- 1} \overset{p}{\to} Q^{- 1}

by Condition 1. Since

{\hat{γ}}_{i}^{(1)} = 0

r_{i}

from (2.3),

\frac{1}{n} \sum_{i = 1}^{n} {\hat{γ}}_{i}^{(1)} \overset{p}{\to} E (ξ_{i})

. When this result is combined with Condition 4, we have

\frac{1}{n} \sum_{i = 1}^{n} x_{i} {\hat{γ}}_{i}^{(1)} = o_{p} (1)

. Thus,

{\hat{β}}^{(1)} = {\hat{β}}^{OLS} + o_{p} (1)

. Now, we consider

{\hat{γ}}_{i}^{(2)}

. First, we replace

{\hat{γ}}_{i}^{(1)}

with

{\hat{γ}}_{i}^{(2)}

, and

{\hat{β}}^{(0)}

with

{\hat{β}}^{(1)}

in (6.1), which readily shows that

{\hat{γ}}_{i}^{(2)} = O_{p} (1)

. Next, we consider

{\hat{β}}^{(2)}

\begin{matrix} {\hat{β}}^{(2)} & = & (X^{'} X)^{- 1} X^{'} (Y^{(1)} - {\hat{γ}}^{(2)}) \\ = & (X^{'} X)^{- 1} X^{'} (Y^{(0)} - {\hat{γ}}^{(1)} - {\hat{γ}}^{(2)}) \\ = & {\hat{β}}^{OLS} - (X^{'} X)^{- 1} X^{'} ({\hat{γ}}^{(1)} + {\hat{γ}}^{(2)}) . \end{matrix}

Again,

(X^{'} X)^{- 1} X^{'} ({\hat{γ}}^{(1)} + {\hat{γ}}^{(2)})

o_{p} (1)

by the similar argument. Since

\sqrt{n} ({\hat{β}}^{OLS} - β) \overset{d}{\to} N (0, σ^{2} Q^{- 1})

from the central limit theorem, we have the same result for

{\hat{β}}^{(2)} + (X^{'} X)^{- 1} X^{'} ({\hat{γ}}^{(1)} + {\hat{γ}}^{(2)})

, and this holds for general m with

{\hat{β}}^{(m)}

, which completes the proof.

□

Footnotes

Acknowledgments

Jung's work was supported by University of Waikato Research Trust (#103406), Lee's work was supported by Hankuk University of Foreign Studies Research Fund of 2015 and Hu's work was partially supported by the National Institute of Health Grants R01GM080503 and R01CA158113.

References

Andrews

Pregibon

(1978) Finding the outliers that matter. Journal of the Royal Statistical Society, Series B, Methodological , 40, 85–93.

Beaton

Tukey

(1974) The fitting of power series, meaning polynomials, illustrated on band-spectroscopic data. Technometrics , 16, 147–185.

Bloomfield

Steiger

(1983) Least Absolute Deviations: Theory, Applications and Algorithms (1st ed.). Boston: Birkhauser Boston ISBN 0817631577, 9780817631574.

Breckling

Chambers

(1988) M-quantiles. Biometrika , 75, 761–771.

Brownlee

(1965) Statistical Theory and Methodolooy in Science and Engineering (2nd ed.). New York: Wiley.

Carroll

Ruppert

(1985) Transformations in regression: A robust analysis. Technometrics , 27, 1–12.

Chambers

Heathcote

(1981) On the estimation of slope and the identification of outliers in linear regression. Biometrika , 68, 21–33.

Daniel

Wood

(1971) Fitting Equations to Data: Computer Analysis of Multifactor Data for Scientists and Engineers . New York: Wiley.

Daniel

Wood

(1980) Fitting Equations to Data: Computer Analysis of Multifactor Data (2nd ed.). New York: John Wiley & Sons Inc. ISBN 0471053708.

10.

Dempster

Gasko-Green

(1981) New tools for residual analysis. The Annals of Statistics , 9, 945–959.

11.

Dempster

Laird

Rubin

(1977) Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society, Series B, Methodological , 39, 1–38.

12.

Denby

Mallows

(1977) Two diagnostic displays for robust regression analysis. Technometrics , 19, 1–13.

13.

Fan

(2001) Variable selection via nonconcave penalized likelihood and its oracle properties. Journal of the American Statistical Association , 96, 1348–1360.

14.

Games

Witten

Hastie

Tibshirani

(2013) An Introduction to Statistical Learning: With Applications in R . New York: Springer-Verlag. Retrieved 29 December 2015 from www.StatLearning.com\

15.

Hadi

Simonoff

(1993) Procedures for the identification of multiple outliers in linear models. Journal of the American Statistical Association , 88, 1264–1272.

16.

Hoeting

Raftery

Madigan

(1996) A method for simultaneous variable selection and outlier identification in linear regression. Computational Statistics & Data Analysis , 22, 251–270.

17.

Huber

Ronchetti

(2009) Robust Statistics (2nd ed.). Hoboken, NJ: John Wiley & Sons Inc. ISBN 978-0-470-12990-6.

18.

Huber

(1973) Robust regression: Asymptotics, conjectures, and monte carlo. The Annals of Statistics , 1, 799–821.

19.

Koenker

(2005) Quantile Regression . New York, NY: Cambridge University Press. ISBN 0521608279.

20.

Kosinski

(1999) A procedure for the detection of multivariate outliers. Computational Statistics & Data Analysis , 29, 145–161.

21.

Lange

Little

RJA

Taylor

JMG

(1989) Robust statistical modeling using the t distribution. Journal of the American Statistical Association , 84, 881–896.

22.

Lee

MacEachern

Jung

(2007) Regularization of case-specific parameters for robustness and efficiency. Technical Report No. 799, The Ohio State University, Columbus, OH.

23.

Lee

MacEachern

Jung

(2012) Regularization of case-specific parameters for robustness and efficiency. Statistical Science , 27, 350–372.

24.

Pena

Yohai

(1999) A fast procedure for outlier diagnostics in large regression problems. Journal of the American Statistical Association , 94, 434–445.

25.

Rockafellar

(1997) Convex Analysis (Princeton Landmarks in Mathematics and Physics) . Princeton, NJ: Princeton University Press.

26.

Rosset

Zhu

(2007) Piecewise linear regularized solution paths. The Annals of Statistics , 35, 1012–1030.

27.

Rousseeuw

Croux

(1993) Alternatives to the median absolute deviation. Journal of the American Statistical Association , 88, 1273–1283.

28.

Rousseeuw

Leroy

(2005) Robust Regression and Outlier Detection . Hoboken, NJ: John Wiley & Sons Inc. ISBN 9780471852339.

29.

Tibshirani

(1996) Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society, Series B, Methodological , 58, 267–288.

30.

Lange

(2008) Coordinate descent algorithms for lasso penalized regression. The Annals of Applied Statistics , 2, 224–244.

Robust regression for highly corrupted response by shifting outliers

Abstract

Keywords

1 Introduction

2.1 Model estimation

3.1 Outlier-shifting in robust loss function

4 Simulation studies

Table 1

Point estimate of MSE, and its standard error (in parentheses) multiplied by 1000 based on 500 simulated data sets obtained by OLS, OLS S , H , H S , med , lasso , lasso S , Hlasso and Hlasso S under σ = 3, 6 and 10 with the contaminated proportions of 10%, 20% and 30%.

Distribution of estimates from 500 simulated samples under the simulation set 1. See the text for the details of simulation setting. Horizontal lines represent the true regression coefficient.

Average false negative rate of three non-zero coefficients and average false positive rate of five zero coefficients (in parenthesis) based on 500 simulated data sets obtained by lasso , lasso S , Hlasso and Hlasso S with σ = 6 and 10 with contaminated proportions of 10%, 20% and 30%

Averaged probability of identifying the outliers with | ε i | > 2.5 based on 500 simulated data sets obtained by OLS S , H S , lasso S and Hlasso S with contaminated proportions of 10%, 20% and 30% and magnitude of σ = 3, 6 and 10

5.1 Analysis of stack loss data

Estimates of β 0 , β 1 , β 2 , β 3 for the stack loss data. Data 1, data 2 and data 3 stand for the data set with 21, 19 and 17 observations, respectively

Sensitivity of five estimators in Table 5.1 at β 1 , β 2 , β 3 with the stack loss data

Residual plots from the four methods of OLS, H , Tukey and OLS S with all observations. In the OLS S , the four observations with * are detected as outliers (observations 1, 3, 4 and 21) and shifted by our method.

Mean of CVS based on 500 random split of mid-Atlantic wage data by models of OLS, OLS S , H , H S , lasso , lasso S , Hlasso and Hlasso S . All values are multiplied by 100.

Appendix

Footnotes

Acknowledgments

References

Averaged probability of identifying the outliers with $| ε_{i} | > 2.5$ based on 500 simulated data sets obtained by ${OLS}^{S}$ , $H^{S}$ , ${lasso}^{S}$ and ${Hlasso}^{S}$ with contaminated proportions of 10%, 20% and 30% and magnitude of $σ$ = 3, 6 and 10

Estimates of $β_{0}$ , $β_{1}$ , $β_{2}$ , $β_{3}$ for the stack loss data. Data 1, data 2 and data 3 stand for the data set with 21, 19 and 17 observations, respectively

Sensitivity of five estimators in Table 5.1 at $β_{1}$ , $β_{2}$ , $β_{3}$ with the stack loss data

Residual plots from the four methods of OLS, $H$ , $Tukey$ and ${OLS}^{S}$ with all observations. In the ${OLS}^{S}$ , the four observations with $*$ are detected as outliers (observations 1, 3, 4 and 21) and shifted by our method.

Mean of $CVS$ based on 500 random split of mid-Atlantic wage data by models of OLS, ${OLS}^{S}$ , $H$ , $H^{S}$ , $lasso$ , ${lasso}^{S}$ , $Hlasso$ and ${Hlasso}^{S}$ . All values are multiplied by 100.