Robust Weighted LAD Regression Avi Giloni, Jeffrey S. Simonoff, and Bhaskar Senguptaasteriskmath February 25, 2005 Abstract The least squares linear regression estimator is well-known to be highly sensitive to unusual observations in the data, and as a result many more robust estimators have been proposed as alternatives. One of the earliest proposals was least-sum of absolute deviations (LAD) regression, where the regression coefficients are estimated through minimization of the sum of the absolute values of the residuals. LAD regression has been largely ignored as a robust alternative to least squares, since it can be strongly affected by a single observation (that is, it has a breakdown point of 1/n, where n is the sample size). In this paper we show that judicious choice of weights can result in a weighted LAD estimator with much higher breakdown point. We discuss the properties of the weighted LAD estimator, and show via simulation that its performance is competitive with that of high breakdown regression estimators, particularly in the presence of outliers located at leverage points. We also apply the estimator to several real data sets. Keywords: Breakdown point; Leverage points; Outliers; Robust regression. asteriskmathAvi Giloni is Assistant Professor, Sy Syms School of Business, Yeshiva University, 500 West 185th St, New York, NY 10033 (e-mail: agiloni@ymail.yu.edu). Jeffrey S. Simonoff is Professor, Leonard N. Stern School of Business, New York University, 44 West 4th Street, New York, NY, 10012 (e-mail: jsimonof@stern.nyu.edu). Bhaskar Sengupta is Section Head, Complex Systems Mod- eling, ExxonMobil Research and Engineering, 1545 Route 22 East, Annandale, NJ 08801 (e-mail: Bhaskar.Sengupta@exxonmobil.com). This work was done while he was at Yeshiva University. The au- thors would like to thank Steve Portnoy for helpful discussion of this material. 1 1 Introduction The linear regression problem is certainly one of the most important data analysis situa- tions (if not the most important). In this situation the data analyst is presented with n observations of a response variable y and some number p of predicting variables x1,...,xp satisfying y = Xbetaf + epsilon, (1) where epsilon is a vector of errors and y = parenlefttp parenleftexparenleftex parenleftexparenleftex parenleftbt y1 ? ? ? yn parenrighttp parenrightexparenrightex parenrightexparenrightex parenrightbt , X = parenlefttp parenleftexparenleftex parenleftexparenleftex parenleftbt x11 ???x1p ?? ?? ?? xn1 ???xnp parenrighttp parenrightexparenrightex parenrightexparenrightex parenrightbt = parenlefttp parenleftexparenleftex parenleftexparenleftex parenleftbt x1 ? ? ? xn parenrighttp parenrightexparenrightex parenrightexparenrightex parenrightbt =(x1,...,xp), epsilon = parenlefttp parenleftexparenleftex parenleftexparenleftex parenleftbt epsilon1 ? ? ? epsilonn parenrighttp parenrightexparenrightex parenrightexparenrightex parenrightbt (note that usually x1 is a column of ones, but this is not required). We assume that X is of full rank (that is, r(X)=p). It is well-known that if the errors epsilon are normally distributed with constant variance the optimal estimator of betaf is the ordinary least squares (OLS) estimator, based on minimizing the lscript2-norm bardbly - Xhatwidebardbl2 = summationtextni=1(yi - xihatwide)2 of the residuals. Unfortunately, it is also well- known that the least squares estimator is very nonrobust, being highly sensitive to unusual observations in the y space (outliers) and X space (leverage points). A way of quantifying this sensitivity is through the notion of the breakdown point of a regression estimator (Hampel, 1968). Suppose we estimate the regression parameters betaf by some technique tau from data (X,y), yielding the estimate betaftau. If we contaminate m (1 <=< m 0orr - i > 0 but not both, then |wi parenleftbigr+i - r-i parenrightbig | = wi parenleftbigr+i + r-i parenrightbig. Therefore, transforming the data by setting parenleftbig?i, ?iparenrightbig = wi parenleftbigxi,yiparenrightbig, implies parenleftbig?i - ?ibetafparenrightbig = wi parenleftbigr+i -r-i parenrightbig. This means that the linear program (4) can be reformulated as min eTn r+ + eTn r- such that wixibetaf + r+ -r- = wiyi for i =1,...,n betaf free, r+ greaterequal 0, r- greaterequal 0. That is, weighted LAD regression can be treated as LAD regression with suitably trans- formed data, and determining the breakdown of weighted LAD regression with known weights corresponds to determining the breakdown of LAD regression with data ( ? , ?). 5 The problem of choosing weights wto maximize the breakdown of the resultant weighted LAD estimator is a nonlinear mixed integer program. Giloni, Sengupta, and Simonoff (2004) showed that this problem is equivalent to a problem related to the knapsack problem, and solved a specififc form of the problem for the simple regression case with a uniform design, resulting in a WLAD breakdown over 30% (LAD regression has a 25% breakdown point in this situation). They also proposed a simple weighting scheme for multiple regression that yields breakdown points over 20% for various two-predictor problems. We propose here an improved weighting scheme that is fast computationally and leads to higher breakdown values. As was noted by Ellis and Morgenthaler (1992), the goal is to downweight observa- tions that are outlying in the predictor space (that is, are leverage points). Unfortunately, usual measures of leverage suffer from masking, in that multiple leverage points near each other can cause them to ?hide? each other. High-breakdown measures of leverage can be calculated, but these are computationally intensive. We will instead use conventional measures of leverage, but will measure outlyingness relative to what is (hopefully) a clean subset of the data. We defifne this subset to be the set of lscript observations that are closest to the coordinatewise medians, where distance is the sum of coordinatewise distances, and the variables have all been scaled to be in the range [0,1]. The size lscript should be large enough to include much of the data, but small enough so that it doesn?t include outlying observations; we use lscript = .6n here. The notion of identifying outlying observations by defifning a clean subset of the data and then measuring the distance of observations relative to that subset was introduced by Rosner (1975) for univariate Gaussian data. The fifrst application of this idea to multivariate and regression data appears to have been in Simonoff (1991). Other applications of this idea for multivariate and regression data can be found in Hadi (1992, 1994), Hadi and Simonoff (1993), Billor, Hadi, and Velleman (2000), and Atkinson and Riani (2000). The defifnition of the clean subset used here is that used by Billor, Hadi, and Velleman (2000) as version 2 of their algorithm for fifnding a clean subset of a multivariate data set (page 285), and as they note this is a robust method of choosing the clean subset. Let XS be this clean subset. The set of leverage values for an observation xi relative to the clean 6 subset is hi = xi(XprimeSXS)-1xiprime, using the usual (least squares) notion of leverage. Since the leverage is proportional to the squared Mahalanobis distance (Chatterjee and Hadi, 1988, page 102), the weight is taken to be wi = radicalbigminj(hj)/hi; that is, it is inversely proportional to the distance from the clean subset. WLAD regression is a particular example of generalized M (GM)-estimation (Mallows, 1975). It is known that the maximum breakdown point of all GM-estimators decreases as a function of p (Maronna, Bustos, and Yohai, 1979), implying that the breakdown point of WLAD regression cannot be arbitrarily high for models with many predictors. In the next section we use Monte Carlo simulations to study the properties of WLAD regression, and show that despite this, WLAD estimation can be competitive with high breakdown estimation even for reasonably large values of p. 3 The Properties of WLAD Regression Estimation We begin this section with discussion of the asymptotic properties of the WLAD estimator. Consider again the regression model (1), and assume that the errors are independent and identically distributed with cumulative distribution function F. Assume that max bardbl xi bardbl = o(n1/4), F is twice differentiable at 0, and f(0) = Fprime(0) > 0. Let ?w be the WLAD estimator of betaf. Let W= diag(w1,...,wn), and assume that the weights are known positive values that satisfy maxwi = O(1) and maxw-1i = O(1). Theorem 1 As n arrowrightinfinity, radicaln( ?w -betaf) is asymptotically p-variate normal with mean 0 and covariance matrix Q-1(XprimeW2X)Q-1omega2, where Q = limnarrowrightinfinity XprimeWX/n and omega =[2f(0)]-1. The proof of this theorem is given in the Appendix. The implication of this theorem is that confifdence regions forbetaf can be constructed based on ?w using the estimated asymptotic covariance, (XprimeWX)-1(XprimeW2X)(XprimeWX)-1 ?2, where omega is estimated in some reasonable way. In the simulations that follow we use a kernel estimator to estimate f(0), and hence omega. 7 We explore the fifnite-sample properties of WLAD regression using Monte Carlo simu- lations, performed using the R package (R Development Core Team, 2004). In addition to the WLAD estimator, we also report results for least squares, (unweighted) LAD, and MM estimators. The MM estimator (Yohai, 1987) uses an inefficient high-breakdown method as an initial estimate, but then uses M-estimation to improve efficiency while still maintain- ing a high breakdown point. We also included an M-estimator (Huber, 1973) and the least trimmed squares (LTS) high breakdown estimator (Rousseeuw, 1985) in the simulations, but the MM estimator consistently outperformed both, so we do not report those results here. The (W)LAD estimators were constructed using the quantreg package (Koenker, 2004), while the LTS and MM estimators were constructed using the MASS package (Venables and Ripley, 2002). We examine various values of sample size n and number of predictors k (note that p = k+1, as the models include an intercept term), and different outlier/leverage point proportion and position combinations. Predictors were generated multivariate normal (with each variable having mean 7.5 and standard deviation 4), with certain observations being modififed to be leverage points in some situations, and 500 simulations replications were generated for each setting. All regression functions had intercept equal to 0 and all slopes equal to 5, with variance of the Gaussian errors equal to 1. Figure 1 summarizes the results of simulations where n = 40 and k = 2. Each bar?s height represents n ? MSE, where MSE is the mean squared error of the slope estimate, separated by predictor, with the shaded portion corresponding to squared bias and the unshaded portion corresponding to variance. Although both predictors were generated to have variance equal to 16, by random chance the second predictor had much lower variability (sample variance 10.1), resulting in higher values of n ? MSE for that predictor?s slope estimate for all methods. As can be seen in the fifgure, the difference in variability of the predictors affects the relative performance of the methods. When there are no outliers, all of the estimators are (virtually) unbiased, and relative efficiency drives the results. As expected, OLS is most efficient, with the MM estimator close. The LAD estimators are less efficient, with WLAD having highest variability (as would be expected, since observations with predictor values farthest from the center are 8 downweighted, and it is these observations that increase efficiency). Relative performance changes markedly in the presence of outliers however, with the nonrobust OLS estimator no longer effective. The fifrst fifve plots refer to outliers with mean 3 standard deviations from the true expected value, while the last (?large outliers?) refers to outliers with mean 7 standard deviations from the true value. The last three plots refer to situations where the observations with outliers fifrst had their predictor values adjusted to make them leverage points (the predictor values were centered roughly 4 standard deviations from the predictor mean). Given that the weighting scheme of the WLAD estimator is designed to downweight leverage points, it is not surprising to see much better performance for WLAD in this situation compared to LAD (in fact, the breakdown point of LAD is 17.5% while that of WLAD is 22.5%, so deteriorating performance for LAD in the 20% outliers case would be expected). That breakdown is not the entire story, however, is clear from the much better performance of WLAD compared to MM, especially when the outliers are at leverage points. This is being driven by much lower squared bias, although the variance of WLAD is also lower when there are 20% large outliers at the leverage points. Figure 2 gives corresponding results where n = 100 and k = 6. Although it is not feasible to calculate exact breakdown values for the (W)LAD estimators, upper bounds on those values can be determined (by running the mixed integer program for many iterations, but not to convergence), and they are 16% (for LAD) and 20% (for WLAD), respectively, in this case for the designs with 20% leverage points. Despite the fact that the breakdown point for WLAD is not greater than the observed percentage of generated outlier obser- vations (implying that breakdown can occur), it is still an effective estimator, once again outperforming the MM estimator in the presence of leverage points. Once again, the es- timator exhibits good bias properties, although that is outweighed by high variance when the data do not contain leverage points. Figure 3 gives results for n = 1000 and k = 20. Even in this case, where the large number of predictors implies a maximum breakdown point of WLAD less than 2% for the design with 20% leverage points, the performance of WLAD is similar to that seen earlier, with the estimator outperforming the MM estimator for large outliers in the presence of leverage points. Thus, it is apparent that the goal of 9 increasing the breakdown point of the estimator is a reasonable one, even when there are more outliers than the breakdown value. This will also be evident in several of the data examples discusses in the next section. We also examined the usefulness of the asymptotic distribution derived in Theorem 1 as a tool for inference, by examining the properties of hypothesis tests based on the approximate normal distribution. Construction of such tests requires estimating omega =[2f(0)]-1, which is done here using a kernel density estimate (see, e.g., section 3.1 of Simonoff, 1996), with the amount of smoothing chosen using the bandwidth selector of Sheather and Jones (1991). The adequacy of the approximation was evaluated by determining the average rejection proportions (empirical sizes) of Wald (Gaussian-based) tests for each coefficient of the actual null value at a nominal .05 level and then averaging over all slopes, so the goal would be tests with size close to .05. Results for the simulation situations previously examined are summarized in Table 1. We give results for the actual tests, and also for tests where the coefficient estimates are recentered at the null value, so that the bias of the estimators does not affect performance. For the n = 1000 case with leverage points, separate average empirical sizes are given for the predictors with and without leverage points for the uncorrected tests, since their performances are very different. It can be seen that while the bias seen earlier (especially in the leverage point case) results in very anticonservative tests (which would presumably also be true for the other estimators, since they are even more biased), when this bias is corrected, the tests have size reasonably close to .05. Thus, the evidence suggests that the assumption of a Gaussian distribution for the coefficient estimator, and the implied estimates of standard errors of the coefficient estimates, are reasonable even for small samples. 4 Application to Real Datasets In this section we apply the WLAD estimator to several well-known datasets from the ro- bustness and outlier identififcation literature. The Hertzsprung-Russell stars data (Rousseeuw and Leroy, 1987, page 27) consist of measurements of the logarithm of light intensity versus logarithm of effective surface temperature for 47 stars. The data are given in the leftmost 10 plot in Figure 4. Although there is a direct relationship between the two variables for most of the observations, four stars have low temperature with high light intensity (these are so- called ?red giant? stars). The OLS fift (dotted line) is drawn to these outliers, as is the LAD fift (dashed line), but the WLAD fift (solid line) follows the general pattern of the points well (note that the LAD estimate actually passes through one of the outliers). It might be thought that the improved performance of WLAD over LAD is because of its higher breakdown point, but this is not, in fact, the case. The breakdown point of LAD here is 10.6% (5/47), so the four observed outliers are not enough to break down the estimator. This can be seen in the middle plot of Figure 4, where the four outliers have had their responses adjusted upwards by 10. The OLS line continues to follow the points, but the LAD line is virtually unchanged (and no longer passes through any of the outliers), because it has not broken down (not surprisingly, the WLAD line is unchanged). Thus, WLAD provides a better fift than LAD even when LAD has not broken down. On the other hand, the improved breakdown point of WLAD (which is 14/47 = 29.8%) becomes apparent if two more values are perturbed upwards (the rightmost plot of Figure 4). With six outliers LAD has broken down, and tracks the unusual points in a similar way to OLS (once again passing through one of them), while WLAD still follows the bulk of the points. Hawkins, Bradu, and Kass (1984) constructed an artififcial three-predictor data set with 75 observations, where outliers were placed at cases 1?10. The LAD estimator has break- down point 10.7% (8/75) and breaks down, with the fiftted regression hyperplane going through one of the outlier points (it goes through cases 5, 20, and 32). In contrast the WLAD estimator, with breakdown point 20% (15/75) works well, going through cases 18, 25, and 30, with each of the outlier observations having absolute residual at least 2.3 times that of any of the clean points. The fiftted WLAD regression is Y = -0.446 + 0.159X1 +0.090X2 -0.032X3, with z-statistics for the three slopes being 1.15, 0.78, and -0.37, respectively. That is, there is little evidence for any relationship here, which is consistent with the way the data were constructed, as none of the t-tests for the three predictors in an OLS fift on observations 11 11-75 are statistically signififcant (in contrast to the fift on all of the observations, where the outliers result in variables X2 and X3 being signififcant predictors). The fifnal data set examined here is the modififed wood gravity of Rousseeuw (1984). These data are based on a real data set (with n = 20 and k = 5), but were modififed to have outliers at cases 4, 6, 8, and 19. Both LAD and WLAD have the same breakdown point (15% = 3/20), but despite this, while LAD performs poorly (passing through the outlier case 8), WLAD performs well, with each of the four outlier cases having absolute residual more than 7.7 times that of any of the clean points. Thus, the WLAD weighting is benefifcial even when breakdown is not improved. The fiftted WLAD regression is Y =0.387 + 0.321X1 - 0.422X2 -0.541X3 - 0.336X4 +0.523X5, with z-statistics for the fifve slopes being 8.50,-2.64, -15.18, -6.32, and 7.79, respectively. An OLS fift on the clean data also identififes predictors 1, 3, 4, and 5 as being most important, although in that case the coefficient for variable 2 is not statistically signififcant (in contrast to an OLS fift on the entire data set, where variables 4 and 5 are insignififcant, because of the effect of the outliers). 5 Conclusion In this paper we have proposed a weighted version of LAD regression designed to increase the breakdown of the estimator that is easy to compute and has performance competi- tive with high breakdown estimators, particularly in the presence of leverage points. These weighting ideas also apply to other estimators. Given the good bias properties of the WLAD estimator, it is reasonable to wonder if a similar weighting scheme used for a more efficient GM-estimator, such as that of Krasker and Welsch (1982), would reduce the variance while preserving robustness, and be even more effective than WLAD. The LAD estimator is a special case of regression quantile estimators (Koenker and Bassett, 1978; Koenker, 2000), which have been shown to be useful in highlighting interesting structure in regression prob- lems, including nonnormality and heteroscedasticity in the error distribution; it would be interesting to see if weighted versions of such estimators would be more resistant to the 12 effects of unusual observations. Appendix Consider the model y = Xbetaf + deltanu, where deltan = diag(delta1,...,deltan), betaf elementRfracturp, and u is a vector of independent and identically distributed errors with cumulative distribution function F. Assume that the deltais are known values satisfying maxdeltai = O(1) and maxdelta-1i = O(1), and that maxbardbl xi bardbl = o(n1/4). Assume that F is twice differentiable at 0, and f(0) = Fprime(0) > 0. Let Qn = Xprimedelta-1n X/n = Q + O(n-1/4 logn). Let ? be the least absolute deviation (LAD) estimator of betaf that minimizes nsummationdisplay i=1 |yi - xibetaf|. Lemma 1 As n arrowrightinfinity, ?betaf - betaf = n-1Q-1n f(0) nsummationdisplay i=1 xiprimePsi(ui)+O((logn/n)3/4), where Psi(x)=I(x<0) -.5. Proof: This result follows from an adaptation of the proof of Theorem 2.1 of Zhou and Portnoy (1998) (hereafter ZP). In that theorem the multipliers deltai are estimated based on a linear function of a parameter vector gamma, so the results quoted here are based on taking ? = gamma in that proof. Let Wn(t)=summationtextni=1 xiprimePsi(ui). By Lemma A.1 of ZP, Wn(? -betaf)=O(n-3/4). (5) Let Mn = c0n-1/2(logn)1/2. Then, by Lemma A.2 of ZP, sup bardbltbardbl<=