Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2700846

Formats

Article sections

- Summary
- 1. Introduction
- 2. Weighted Wilcoxon-type Smoothly Clipped Absolute Deviation Method
- 3. Asymptotic Properties
- 4. Numerical Examples
- 5. Summary
- Supplementary Material
- References

Authors

Related links

Biometrics. Author manuscript; available in PMC 2010 June 1.

Published in final edited form as:

Published online 2008 July 18. doi: 10.1111/j.1541-0420.2008.01099.x

PMCID: PMC2700846

NIHMSID: NIHMS103887

Lan Wang, School of Statistics, University of Minnesota, 224 Church Street SE, Minneapolis, MN 55455, U.S.A., email: ude.nmu.tats@nal;

The publisher's final edited version of this article is available at Biometrics

See other articles in PMC that cite the published article.

Shrinkage-type variable selection procedures have recently seen increasing applications in biomedical research. However, their performance can be adversely influenced by outliers in either the response or the covariate space. This paper proposes a weighted Wilcoxon-type smoothly clipped absolute deviation (WW-SCAD) method, which deals with robust variable selection and robust estimation simultaneously. The new procedure can be conveniently implemented with the statistical software R. We establish that the WW-SCAD correctly identifies the set of zero coefficients with probability approaching one and estimates the nonzero coefficients with the rate *n*^{−1/2}. Moreover, with appropriately chosen weights the WW-SCAD is robust with respect to outliers in both the **x** and *y* directions. The important special case with constant weights yields an oracle-type estimator with high efficiency at the presence of heavier-tailed random errors. The robustness of the WW-SCAD is partly justified by its asymptotic performance under local shrinking contamination. We propose a BIC-type tuning parameter selector for the WW-SCAD. The performance of the WW-SCAD is demonstrated via simulations and by an application to a study that investigates the effects of personal characteristics and dietary factors on plasma beta-carotene level.

In biomedical research, statisticians often need to analyze data sets with a non-normally distributed response variable and/or many covariates that potentially contain multiple high leverage points. This often imposes serious problems for variable selection and the subsequent inference. Existing work on robust variable selection are mostly robust best-subset procedures, such as robust AIC or BIC, see Ronchetti (1985), Hurvich and Tsai (1990), Burman and Nolan (1995), Ronchetti and Staudte (1994), Ronchetti, Field and Blanchard (1997), Wisnowski et al. (2003) and Müller and Welsh (2005), among others. The best-subset type procedures are computationally intensive even for moderately large number of covariates; and are known to have inherent instability (Brieman, 1996) due to their discrete nature. Moreover, these approaches in general are only robust against outliers in the response space but are still sensitive to high-leverage points. This paper introduces a novel unified framework called the weighted Wilcoxon-type smoothly clipped absolute deviation method (WW-SCAD, for short) for automatic robust variable selection and robust estimation that can effectively handle the above concerns.

In Section 4, we analyzed a data set from a study investigating the effects of personal characteristics and dietary factors on plasma beta-carotene level. It has been observed that low plasma concentrations of beta-carotene might be associated with increased risk of developing certain types of cancer. Due to the nature of the study, many patients have rather low plasma beta-carotene levels. This results in a long-tailed and highly skewed distribution for the response variable (plasma beta-carotene level, ng/ml), see the histogram depicted in Figure 1(a). Also, two of the ten covariates: *x*_{8} (number of alcoholic drinks consumed per week) and *x*_{9} (cholesterol consumed per day) clearly contain multiple outliers, some are even quite extreme, as revealed by their boxplots in Figure 1(b). This leads us to propose a procedure that is robust on both the covariate and response spaces to analyze this data set. Some covariates may not have effects on the plasma beta-carotene level. Thus, it is of great interest to further develop variable selection procedure in robust statistical modeling. From the analysis in section 4.2, the newly proposed robust variable selection procedure reduces the median prediction error on the validation data to about 68% of that given by its nonrobust alternative.

Plasma beta-carotene level data: (a) histogram of *y*, (b) boxplots of *x*_{8} (number of alcoholic drinks consumed per week) and *x*_{9} (cholesterol consumed per day)

The WW-SCAD procedure is motivated by recent developments in shrinkage-type variable selection procedures such as LASSO (Tibshirani, 1996) and SCAD (Fan and Li, 2001). Distinguished from the robust subset-type procedures, the WW-SCAD simultaneously selects covariates and estimates parameters by minimizing an objective function which is the sum of the weighted Wilcoxon-type dispersion function and the smoothly clipped absolute deviation (SCAD) penalty function, see Section 2.2. The penalty term shrinks the estimated small coefficients to zero, thus serves the purpose of variable selection.

The WW-SCAD is robust against outliers in both the **x** and *y* directions with appropriately chosen weights. This is different from the LAD-LASSO procedure based on the least absolute deviation regression (Wang, Li and Jiang, 2007) and the penalized composite quantile regression (Zou and Yuan, 2007), which provide a certain degree of protection against outliers in the response space but are vulnerable to high leverage points. We provide theoretical justification for the robustness of the WW-SCAD by studying its performance under shrinking local contamination. Under the local contamination, we reveal that the WW-SCAD still identifies zero coefficients with probability approaching one and estimates nonzero coefficients with a bias bounded in (**x**, *y*) when the weights are appropriately chosen.

The WW-SCAD with constant weights leads to an important special case that is closely related to the classical Wilcoxon inference based on Jaeckel’s (1972) dispersion function with Wilcoxon scores. In this case, with a proper tuning parameter the resulted estimator possesses the oracle property (Fan and Li, 2001) and often significantly improves the efficiency of the LS-SCAD (least-squared procedure with SCAD penalty) in the presence of heavy-tailed errors. The tuning parameter in the WW-SCAD controls the model complexity and plays an important role in the variable selection procedure. In practice, it is desirable to select the tuning parameter using a data-driven method. We propose a BIC-type tuning parameter selector and show that with probability tending to one, the WW-SCAD with the BIC-selector can identify the most parsimonious correct model.

Rank-based statistical procedures have wide applications in biomedical research due to their robustness and high efficiency; see Jin et al. (2003), Jung and Ying (2003), Mahfoud and Randles (2005), Rosner, Glynn and Lee (2006a, 2006b), Heller (2007), Datta and Satten (2008), Wang and Zhao (2008) and the references therein. However, the aforementioned work mainly focuses on estimation and hypothesis testing. Our proposal therefore extends rank-based nonparametric analysis to the important area of variable selection.

The rest of the paper is structured as follows. In the next section, we introduce the WW-SCAD procedure and discuss its implementation via the software package R. In Section 3, we establish the asymptotic normality and consistency of selection, and provide justification for robustness by considering the asymptotic distribution under local contamination. Furthermore, we introduce a BIC-type procedure for selecting the tuning parameter. In Section 4, we demonstrate the performance of the WW-SCAD by Monte Carlo studies and apply it to analyze the plasma beta-carotene level data set. Section 5 summarizes the paper.

Consider a linear regression model

$$\mathbf{Y}=\alpha {\mathbf{1}}_{n}+\mathbf{X}\mathit{\beta}+\mathit{\epsilon},$$

where **Y** = (*Y*_{1},…, *Y _{n}*)′ is an

The penalized weighted Wilcoxon method estimates ** β** by minimizing

$${n}^{-1}\sum _{i<j}{b}_{ij}\mid {e}_{i}-{e}_{j}\mid +n\sum _{j=1}^{d}{p}_{\lambda}(\mid {\beta}_{j}\mid ),$$

where the *b _{ij}*’s are positive and symmetric weights,
${e}_{i}={Y}_{i}-{\mathbf{x}}_{i}^{\prime}\mathit{\beta}$ with

Directly minimizing *n*^{−2} Σ_{i}_{<}* _{j} b_{ij}*|

Fan and Li (2001) provided deep insights into the principles of choosing an appropriate penalty function. They proposed the smoothly clipped absolute deviation (SCAD) penalty function, which satisfies *p _{λ}*(0) = 0 and has the first-order derivative

$${p}_{\lambda}^{\prime}(\theta )=\lambda \left\{I(\theta \le \lambda )+\frac{{(a\lambda -\theta )}_{+}}{(a-1)\lambda}I(\theta >\lambda )\right\}$$

(1)

for some *a* > 2 and *θ* > 0. Following Fan and Li (2001), we will use *a* = 3.7 throughout this paper. Recently, Zou and Li (2007) proposed a local linear approximation to the SCAD penalty, which retains the same asymptotic properties and at the same time significantly improves the computational efficiency of Fan and Li’s LS-SCAD. Adopting this idea, we propose a WW-SCAD procedure for robust simultaneous variable selection and estimation. Formally, the WW-SCAD method estimates ** β** by minimizing

$${n}^{-1}\sum _{i<j}{b}_{ij}\mid {e}_{i}-{e}_{j}\mid +n\sum _{j=1}^{d}{p}_{\lambda}^{\prime}(\mid {\beta}_{j}^{0}\mid )\mid {\beta}_{j}\mid ,$$

(2)

where
${p}_{\lambda}^{\prime}(\xb7)$ is defined in (1),
${\sum}_{j=1}^{d}{p}_{\lambda}^{\prime}(\mid {\beta}_{j}^{0}\mid )\mid {\beta}_{j}\mid $ is the linearized SCAD penalty (Zou and Li, 2007) and *β*^{0} is an initial estimator, which we set to be the unpenalized weighted Wilcoxon estimator. Unlike the objective function for the LS-SCAD, the objective function defined in (2) is convex in ** β**.

We complete this subsection with a brief discussion of estimating the intercept parameter *α.* Since (2) is invariant to a location change, *α* cannot be estimated simultaneously with ** β**. Instead,

An appealing feature of the WW-SCAD is that its computation can be conveniently carried out using the statistical software R. Our algorithm is similar to that of the LAD-LASSO (Wang, Li and Jiang, 2007). The key observation is that minimizing (2) can be achieved by fitting an *L*_{1} regression model based on the pseudo observations
$({\mathbf{x}}_{k}^{\ast},{Y}_{k}^{\ast}),k=1,\dots ,{\scriptstyle \frac{n(n-1)}{2}}+d$. The first *n*(*n* − 1)/2 pseudo observations correspond to (*b _{ij}*(

The unpenalized weighted Wilcoxon estimator
${\beta}_{j}^{0}$ can obtained using the function *wwfit* in the R software developed by Terpstra and McKean (2005) for rank regression (downloadable from http://www.stat.wmich.edu/mckean/HMC/Rcode). And the *L*_{1} regression model can be fitted using the R package *quantreg* by Roger Koenker for quantile regression.

It is worth emphasizing that in order to achieve practical robustness it is not sufficient to merely have a robust objective function. The algorithm itself is also of critical importance. For shrinkage-type procedures, special algorithms are often used for the implementation and most of these algorithms are sensitive to outlier contamination. As an example, if we approximate the WW-SCAD objective function quadratically and then apply the LARS algorithm (Efron et al., 2004), the resulting procedure is likely to be still sensitive to outliers.

We consider Mallows-type weights *b _{ij}* that possibly depend on the covariates in the form

$$h({\mathbf{x}}_{i})=min\left\{1,\frac{b}{({\mathbf{x}}_{i}-\widehat{\mu}{)}^{\prime}{S}^{-1}({\mathbf{x}}_{i}-\widehat{\mu})}\right\},$$

with (*, S*) being the robust minimum volume ellipsoid (MVE) estimator of the location and scatter (Rousseeuw and van Zomeren, 1990), and *b* being the 95th percentile of *χ*^{2}(*d*).

Following the notations in Naranjo and Hettmansperger (1994), let **W** be an *n × n* matrix of elements *w _{ij}*, where

$${w}_{ij}=\{\begin{array}{ll}-{n}^{-1}{b}_{ij},\hfill & i\ne j\hfill \\ {n}^{-1}{\sum}_{k\ne i}{b}_{ik},\hfill & i=j.\hfill \end{array}$$

Assume
${n}^{-1}{\mathbf{X}}^{\prime}\mathbf{WX}\stackrel{p}{\to}\mathbf{C},{n}^{-1}{\mathbf{X}}^{\prime}{\mathbf{W}}^{2}\mathbf{X}\stackrel{p}{\to}\mathbf{V}$ and
${n}^{-1}{\mathbf{X}}^{\prime}\mathbf{X}\stackrel{p}{\to}\mathbf{\sum}$, where **C**, **V** and **Σ** are positive definite matrices:

$$\begin{array}{l}\mathbf{C}=\frac{1}{2}\int \int ({\mathbf{x}}_{2}-{\mathbf{x}}_{1})({\mathbf{x}}_{2}-{\mathbf{x}}_{1}{)}^{\prime}b({\mathbf{x}}_{1},{\mathbf{x}}_{2})dM({\mathbf{x}}_{2})dM({\mathbf{x}}_{1}),\\ \mathbf{V}=\int \{({\mathbf{x}}_{2}-{\mathbf{x}}_{1})b({\mathbf{x}}_{2},{\mathbf{x}}_{1})dM({\mathbf{x}}_{2})\}\{({\mathbf{x}}_{2}-{\mathbf{x}}_{1})b({\mathbf{x}}_{2},{\mathbf{x}}_{1})dM({\mathbf{x}}_{2}){\}}^{\prime}dM({\mathbf{x}}_{1}),\\ \mathbf{\sum}=\frac{1}{2}\int \int ({\mathbf{x}}_{2}-{\mathbf{x}}_{1})({\mathbf{x}}_{2}-{\mathbf{x}}_{1}{)}^{\prime}dM({\mathbf{x}}_{2})dM({\mathbf{x}}_{1}),\end{array}$$

and *M*(**x**) denotes the cumulative distribution function of **x**.

We denote the true value of ** β** by
${\mathit{\beta}}_{0}=({\beta}_{10},\dots ,{\beta}_{d0}{)}^{\prime}=({\mathit{\beta}}_{10}^{\prime},{\mathit{\beta}}_{20}^{\prime}{)}^{\prime}$. Without loss of generality, we assume that

$$\mathbf{C}=\left[\begin{array}{cc}{\mathbf{C}}_{11}& {\mathbf{C}}_{12}\\ {\mathbf{C}}_{21}& {\mathbf{C}}_{22}\end{array}\right],\phantom{\rule{0.38889em}{0ex}}\mathbf{V}=\left[\begin{array}{cc}{\mathbf{V}}_{11}& {\mathbf{V}}_{12}\\ {\mathbf{V}}_{21}& {\mathbf{V}}_{22}\end{array}\right]\phantom{\rule{0.38889em}{0ex}}\text{and}\phantom{\rule{0.38889em}{0ex}}\mathbf{\sum}=\left[\begin{array}{cc}{\mathbf{\sum}}_{11}& {\mathbf{\sum}}_{12}\\ {\mathbf{\sum}}_{21}& {\mathbf{\sum}}_{22}\end{array}\right].$$

In addition to the above, we assume that the error density function *f*(·) is absolutely continuous with finite Fisher information, i.e., *∫*{*f*(*x*)}^{−1} *f*′(*x*)^{2}*dx* < ∞. And **X** and **WX** both satisfy Huber’s condition, a sufficient and necessary condition for the least-squares estimator to have an asymptotic normal distribution; see condition (D.2) of Hettmansperger and McKean (1998). Under these conditions, the unpenalized WW estimator is
$\sqrt{n}$-consistent for *β*_{0} and asymptotically normal.

Theorem 1 below presents the asymptotic property of the WW-SCAD as a simultaneous model selection and parameter estimation tool, and its proof is given in the Web Appendix.

Assume the regularity conditions in Section 3.1. If λ_{n} → 0 and
$\sqrt{n}{\lambda}_{n}\to \infty $ as n → ∞, then the WW-SCAD estimator
$\widehat{\mathit{\beta}}=({\widehat{\mathit{\beta}}}_{10}^{\prime},{\widehat{\mathit{\beta}}}_{20}^{\prime}{)}^{\prime}$ must satisfy that P(_{20} = **0**) = 1, and

$$\sqrt{n}({\widehat{\mathit{\beta}}}_{10}-{\mathit{\beta}}_{10})\to {N}_{s}(\mathbf{0},{\tau}^{2}{\mathbf{C}}_{11}^{-1}{\mathbf{V}}_{11}{\mathbf{C}}_{11}^{-1}),$$

where $\tau ={[\sqrt{12}\int {f}^{2}(u)du]}^{-1}$.

The case with constant weight *b _{ij}* 1 is particularly important due to its simplicity and its close connection with the familiar Wilcoxon inference. In this case, we have
${\mathbf{C}}_{11}={\mathbf{V}}_{11}={\mathbf{\sum}}_{11}={\mathbf{X}}_{1}^{\prime}{\mathbf{X}}_{1}$, thus we have the following corollary.

Assume the conditions in Theorem 3.1, then when *b*_{ij} 1, satisfies: P(_{20} = **0**) = 1 and
$\sqrt{n}({\widehat{\mathit{\beta}}}_{10}-{\mathit{\beta}}_{10})\to {N}_{s}(\mathbf{0},{\tau}^{2}{\mathbf{\sum}}_{11}^{-1})$.

Corollary 1 suggests that the Wilcoxon-SCAD, with *b _{ij}* 1 and a properly chosen tuning parameter, possesses the oracle property (Fan and Li, 2001). That is, with probability approach one, the WW-SCAD can correctly identify the nonzero coefficients, and estimate them as efficiently as the unpenalized WW rank regression does as if the true model were known in advance. Moreover, the WW-SCAD can be more efficient than the LS-SCAD for estimating

This asymptotic relative efficiency is the same as that of the one-sample Wilcoxon test with respect to the *t*-test. It is well known in the literature of rank analysis that the *ARE* is as high as 0.955 for normal error distribution, and can be significantly higher than 1 for many heavier-tailed distributions. For instance, *ARE* = 1:5 for the double exponential distribution, and *ARE* = 1:9 for the *t* distribution with 3 degrees of freedom. For symmetric error distributions with finite Fisher information, this asymptotic relative efficiency is known to have a lower bound equal to 0.864.

The asymptotic covariance matrix of _{10} in Corollary 1 can be shown to be equivalent to that in Theorem 2.1 of Zou and Yuan (2007) for composite quantile regression when *K*, the number of quantiles, goes to infinity. The composite quantile regression, however, is more computationally involved. Its objective function involves a mixture of quantile objective functions at different quantiles (the suggested value of *K* for practical use is 19). As a result, besides the regression parameters one also needs to estimate *K* additional parameters corresponding to *K* different quantiles of the error distribution.

With nonconstant weights *b _{ij}*, the WW-SCAD still consistently selects the correct model; however the asymptotic covariance matrix of
$\sqrt{n}({\widehat{\mathit{\beta}}}_{10}-{\mathit{\beta}}_{10})$ is slightly different from what one would obtain if the true model were known in advance. This can be seen by observing that in
${\mathbf{C}}_{11}={\mathbf{X}}_{1}^{\prime}\mathbf{W}{\mathbf{X}}_{1}$ (and similarly

Now we study the robustness property of the WW-SCAD. For robust estimation and hypothesis testing, the influence function approach offers a convenient and essential way to investigate the local robustness (Hampel, 1974, Hampel et al., 1986). This approach, however, is not adequate in the current setting as we perform variable selection and parameter estimation simultaneously.

Following the spirit of influence function approach, we directly study the performance of the WW-SCAD under infinitesimal contamination. Specifically, we consider the following local shrinking contamination as in Heritier and Ronchetti (1994):

$${H}_{n}^{\ast}(\mathbf{x},y)=\left(1-\frac{\delta}{\sqrt{n}}\right)H(\mathbf{x},y)+\frac{\delta}{\sqrt{n}}{\mathrm{\Delta}}_{({\mathbf{x}}^{\ast},{y}^{\ast})},$$

(3)

where *H*(**x***, y*) is the joint cumulative distribution function of the underlying distribution without contamination, Δ_{(}_{x}_{*,}_{y}_{*)} represents the point mass at (**x**^{*},*y*^{*}) and δ is a constant.

Assume the regularity conditions in Section 3.1. If λ_{n} → 0 and
$\sqrt{n}{\lambda}_{n}\to \infty $ as n → ∞, then under local shrinking contamination (3), the WW-SCAD estimator
$\widehat{\mathit{\beta}}=({\widehat{\mathit{\beta}}}_{10}^{\prime},{\widehat{\mathit{\beta}}}_{20}^{\prime}{)}^{\prime}$ must satisfy that P(_{20} = **0**) = 1, and

$$\sqrt{n}({\widehat{\mathit{\beta}}}_{10}-{\mathit{\beta}}_{10})\to {N}_{s}(\eta ,{\tau}^{2}{\mathbf{C}}_{11}^{-1}{\mathbf{V}}_{11}{\mathbf{C}}_{11}^{-1}),$$

where η = δ[2F (y^{*} − **x**^{*}**β**_{0}) − 1] ∫ b(**x**^{*}, **x**)(**x**^{*} − **x**)dM(**x**).

The proof of Theorem 2 is given in the Web Appendix. Theorem 2 indicates that under the local contamination (3), the WW-SCAD can still correctly identifies the set of zero coefficients with probability tending to one; but the contamination introduces a bias *η* in estimating the nonzero coefficients. Note that [2*F*(*y*^{*} − **x**^{*}*β*_{0}) − 1] ∫ *b*(**x**^{*}, **x**)(**x**^{*} − **x**)*dM*(**x**) is also the core part of the influence function of the unpenalized weighted Wilcoxon estimator. The bias *η* is bounded in *y**. With proper choice of weights *b _{ij}*, such as the GR weights introduced in Section 3.1,

The tuning parameter *λ* controls the model complexity and plays a critical role in the WW-SCAD procedure. It is desirable to select *λ* automatically by a data-driven method. Here we propose to select *λ* for the WW-SCAD by minimizing

$$BI{C}_{\lambda}=log\left({n}^{-2}\sum _{i<j}{b}_{ij}\mid ({Y}_{i}-{\mathbf{x}}_{i}^{\prime}{\widehat{\mathit{\beta}}}_{\lambda})-({Y}_{j}-{\mathbf{x}}_{i}^{\prime}{\widehat{\mathit{\beta}}}_{\lambda})\mid \right)+{df}_{\lambda}log(n)/n$$

(4)

over an interval [0, *λ _{max}*], where

To introduce the property of the BIC tuning parameter selector, we next define some notations. We use *S* = {*j*_{1},…,*j _{d}*

For a given candidate model *S*, let *β** _{S}* be the vector of parameters. The

- for any
*S**S*, ${L}_{n}^{S}\stackrel{p}{\to}{L}^{S}$ for some_{F}*L*> 0, where $\stackrel{p}{\to}$ means converges in probability;^{S} - for any
*S**S*, we have_{T}*L*>^{S}*L*.^{ST}

Note that
${L}_{n}^{S}$ is the objective function to obtain the weighted Wilcoxon estimator when model *S* is used. Conditions (1) and (2) are standard for investigating parameter estimation under model misspecification, see White (1981). Let
$R(\mathit{\beta})=0.5\int \int b({\mathbf{x}}_{1},{\mathbf{x}}_{2})\mid ({y}_{1}-{\mathbf{x}}_{1}^{\prime}\mathit{\beta})-({y}_{2}-{\mathbf{x}}_{2}^{\prime}\mathit{\beta})\mid dH({\mathbf{x}}_{1},{y}_{1})dH({\mathbf{x}}_{2},{y}_{2})$. Then for the true model *S _{T}*,

Assume the conditions above and the regularity conditions in Section 3.1, then P(S_{BIC}= S_{T}) → 1:

The proof of Theorem 3 is given in the Web Appendix. Theorem 3 indicates that * _{BIC}* leads to a WW-SCAD estimator which consistently yields the true model. The verification of this theorem is similar to that in Wang, Li and Tsai (2007), in which the authors proposed a novel BIC-selector for the SCAD penalized least squares procedures.

In the literature, the LS-SCAD has been compared with the nonnegative garrote (Breiman, 1995), the LASSO, and the best subset variable selection procedures such as AIC or BIC, see for example, Fan and Li (2001) and Zou and Li (2007). Our simulations are designed to demonstrate the robustness and the efficiency of the WW-SCAD, compared with the LS-SCAD which is computed with the BIC tuning parameter selector of Wang, Li and Tsai (2007). We also compare with the benchmark oracle procedure, which sets the estimate of zero coefficients to be zero and estimates the nonzero coefficients by excluding the covariates of zero coefficients.

We focus on examining the performance of the WW-SCAD in terms of model complexity and model errors (ME) defined by

$$\text{ME}(\widehat{\mathit{\beta}})=(\widehat{\mathit{\beta}}-{\mathit{\beta}}_{0}{)}^{\prime}E({\mathbf{x}}_{1}{\mathbf{x}}_{1}^{\prime})(\widehat{\mathit{\beta}}-{\mathit{\beta}}_{0}).$$

(5)

As in Tibshirani (1996) and Fan and Li (2001), data are generated from

$${Y}_{i}={\mathbf{x}}_{i}^{\prime}\mathit{\beta}+{\mathit{\epsilon}}_{i},\phantom{\rule{0.38889em}{0ex}}i=1,\dots ,100,$$

(6)

where ** β** = (3, 1.5, 0, 0, 2, 0, 0, 0) and

Simulation results are summarized in Table 1, in which we report the average number of correct 0’s (the average number of the five true zero coefficients that are correctly estimated to be zero) and the average number of incorrect 0’s (the average number of the three true zero coefficients that are incorrectly estimated to be zero). We also report the proportion of correctly fitted models. To evaluate the lack-of-fit of the selected model, we report the relative model error RME=ME/ME_{WilFull}, where ME_{WilFull} is the ME for fitting the full model with unpenalized Wilcoxon rank regression.

The simulation results are based on 500 runs. *C* is the average number of correct zeros; *IC* is the average number of incorrect zeros; Correct Fit (%) is the proportion of times the correct model is selected; and MRME is the median of relative model error. **...**

From Table 1, we can see that the median of the RME of the WW-SCAD is close to that of the *WW*_{oracle}, the weighted Wilcoxon estimator from the oracle procedure. In terms of model error, the performance of the WW-SCAD is similar to the LS-SCAD for normal error, but much better than the LS-SCAD for both *t*_{3} error and contaminated normal error in terms of model error. And the WW-SCAD gives significantly higher percentage of correctly fitted 0’s compared to the LS-SCAD.

We now investigate the effect of outliers in the **x** direction on model selection. For this purpose, we consider the same regression model (6) with the standard normal random errors. We consider a contamination of the covariate **x** by replacing a random 5% of **x** with **x** + **e**, where **e** = (*e*_{1},…,*e*_{8})′, with *e*_{3} having a *N* (0, 5) distribution and all the other *e _{i}*’s having independent

In this example, the relative model error is defined as RME=ME/ME_{GRFull}, where ME_{GRFull} is the model error obtained by fitting the full model with the unpenalized weighted Wilcoxon procedure and the GR weights. We use the weighted Wilcoxon procedure with the GR weights under the true mode as the benchmark here. The simulation results are summarized in Table 2, from which we observe that the GR weights lead to model selection procedures robust to outliers in the **x** direction; in contrast the performance of the LS-SCAD is adversely affected. We also note that the WW-SCAD with the Wilcoxon weights is not as seriously affected as the LS-SCAD but does not perform as well as the WW-SCAD with the GR weights.

Observational studies have suggested that low plasma concentration of beta-carotene might be associated with increased risk of developing certain types of cancer. We consider a data set from a cross-sectional study that consists of 273 female patients who had an elective surgical procedure during a three-year period to biopsy or remove a lesion of the lung, colon, breast, skin, ovary or uterus that was found to be non-cancerous (Nierenberg et al., 1989). The response variable *y* is the plasma beta-carotene level (ng/ml) and there are ten covariates: *x*_{1} is age, *x*_{2} is smoking status (1=never, 2=former smoker, 3=current smoker), *x*_{3} is quetelet (weight/height^{2}), *x*_{4} denotes vitamin use (1=yes, fairly often, 2=yes, not often, 3=no), *x*_{5} is the number of calories consumed per day, *x*_{6} is grams of fat consumed per day, *x*_{7} is grams of fiber consumed per day, *x*_{8} is number of alcoholic drinks consumed per week, *x*_{9} is cholesterol consumed (mg per day) and *x*_{10} is dietary beta-carotene consumed (mcg per day).

As revealed by Figure 1, the distribution of *y* is highly skewed, while *x*_{8} and *x*_{9} contain some obvious outliers. One may suggest log transform the response variable. However, our preliminary analysis indicates that the log transformed *y* is still nonnormal. And it becomes even harder to find an appropriate transformation for *x*_{8}, which is on ordinal scale. Since the transformation may not remove the outliers and often brings additional issues for interpretability, we choose to analyze the variables on their original scale.

We use the first 200 observations as a training data set to select and fit the model, and use the rest as a validation data set to evaluate the prediction ability (measured by the median absolute prediction error) of the selected model. The *λ* values selected by the BIC criterion are 1.249, 2.834 and 3.028 for the LS-SCAD, the WW-SCAD with the Wilcoxon score, and the WW-SCAD with the GR score, respectively. The resulting estimated models are displayed in Table 3. The LS-SCAD does not exclude any of the ten candidate covariates from the selected model. The WW-SCAD with either the Wilcoxon score or the GR score fits a much more succinct model that includes *x*_{5} and *x*_{10}, and suggests that increased plasma beta-carotene level is associated with increased dietary intake of beta-carotene and reduced number of calories consumed per day. The WW-SCAD with the Wilcoxon score also includes *x*_{9}.

Analysis of plasma beta-carotene level data. Note: The median absolute prediction error is calculated from the validation data set: median APE=median {|*y*_{i} − *ŷ*_{i}|,*i* = 1,…,73}, where *y*_{i} is the *i*th response in the validation data set **...**

In terms of the prediction performance on the validation data, the WW-SCAD with either the Wilcoxon score or the GR score yields a much smaller median absolute prediction error. The median absolute prediction error of the WW-SCAD with the GR score is only 68% of that given by the LS-SCAD.

We propose a novel robust framework called WW-SCAD for simultaneous variable selection and parameter estimation. This new procedure can be conveniently implemented using the statistical software R. It is much less computationally intensive compared with the best subset type procedures. With appropriately chosen weights, the WW-SCAD procedure can effectively handle outliers in both the **x** and *y* directions. Moreover, it loses very little efficiency with normal data and can be much more efficient than the LS-SCAD at the presence of heavier-tailed random errors. Although we have focused on studying the SCAD penalty, without any difficulty our method can be extended with other penalty functions.

The authors would like to thank the Co-Editor for constructive comments that substantially improved an earlier draft, and thank Dr. John Dziak for his help on presentation. Wang’s research is supported by National Science Foundation grant DMS-0706842; and Li’s research is supported by National Institute on Drug Abuse, NIH, 1 R21 DA024260, and National Science Foundation grant DMS-0348869.

Wed Appendix referenced in Section 3 and the R codes for implementing the WW-SCAD procedure are available under the Paper Information link at the the Biometrics website http://www.tibs.org/biometrics.

Lan Wang, School of Statistics, University of Minnesota, 224 Church Street SE, Minneapolis, MN 55455, U.S.A., email: ude.nmu.tats@nal.

Runze Li, Department of Statistics and The Methodology Center, Pennsylvania State University, University Park, PA 16802-2111, U.S.A., email: ude.usp.tats@ilr.

- Breiman L. Better subset regression using the nonnegative garrote. Technometrics. 1995;37:373–384.
- Breiman L. Heuristics of instability and stabilization in model selection. The Annals of Statistics. 1996;24:2350–2383.
- Breiman L, Friedman JH. Estimating optimal transformations for multiple regression and correlation. Journal of the American Statistical Association. 1985;80:580–598.
- Burman P, Nolan D. A general Akaike-type criterion for model selection in robust regression. Biometrika. 1995;82:877–886.
- Chang WH, McKean JW, Naranjo JD, Sheather SJ. High-breakdown rank regression. Journal of the American Statistical Association. 1999;94:205219.
- Datta S, Satten GA. A signed-rank test for clustered data. Biometrics. 2008 to appear. [PubMed]
- Efron B, Hastie T, Johnstone I, Tibshirani R. Least angle regression. The Annals of Statistics. 2004;32:407–499.
- Fan J, Li R. Variable selection via nonconcave penalized likelihood and its oracle properties. Journal of the American Statistical Association. 2001;96:1348–1360.
- Hampel FR. The influence curve and its role in robust estimation. Journal of the American Statistical Association. 1974;69:383–393.
- Hampel FR, Ronchetti EM, Rousseeuw PJ, Stahel WA. Robust Statistics: The Approach Based on Influence Functions. New York: John Wiley; 1986.
- Heller G. Smoothed rank regression with censored data. Journal of the American Statistical Association. 2007;102:552–559.
- Heritier S, Ronchetti E. Robust bounded-influence tests in general parametric models. Journal of the American Statistical Association. 1994;89:897–904.
- Hettmansperger TP, McKean JW. Robust Nonparametric Statistical Methods. London: Arnold; 1998.
- Hurvich CM, Tsai CL. Model selection for least absolute deviations regression in small samples. Statistics & Probability Letters. 1990;9:259–265.
- Jaeckel LA. Estimating regression coefficients by minimizing the dispersion of the residuals. The Annals of Mathematical Statistics. 1972;43:1449–1458.
- Jin Z, Lin DY, Wei LJ, Ying Z. Rank-based inference for the accelerated failure time model. Biometrika. 2003;90:341–353.
- Jung SH, Ying ZL. Rank-based regression with repeated measurements data. Biometrika. 2003;90:732–740.
- Mahfouda ZR, Randles RH. Practical tests for randomized complete block designs. Journal of Multivariate Analysis. 2005;96:73–92.
- Müller S, Welsh AH. Outlier robust model selection in linear regression. Journal of the American Statistical Association. 2005;100:1297–1310.
- Naranjo JD, Hettmansperger TP. Bounded influence rank regression. Journal of the Royal Statistical Society, Series B. 1994;56:209–220.
- Nierenberg DW, Stukel TA, Baron JA, Dain BJ, Greenberg ER. Determinants of plasma levels of beta-carotene and retinol. American Journal of Epidemiology. 1989;130:511–521. [PubMed]
- Ronchetti E. Robust model selection in regression. Statistics & Probability Letters. 1985;3:21–23.
- Ronchetti E, Staudte RG. A robust version of Mallows’
*C*. Journal of the American Statistical Association. 1994;89:550–559._{p} - Ronchetti E, Field C, Blanchard W. Robust linear model selection by cross-validation. Journal of the American Statistical Association. 1997;92:1017–1023.
- Rosner B, Glynn RJ, Lee MLT. The Wilcoxon signed rank test for paired comparisons of clustered data. Biometrics. 2006a;62:185192. [PubMed]
- Rosner B, Glynn RJ, Lee MLT. Extension of the rank sum test for clustered data: two-group comparisons with group membership defined at the subunit level. Biometrics. 2006b;62:12511259. [PubMed]
- Rousseeuw PJ, van Zomeren BC. Unmasking multivariate outliers and leverage points. Journal of the American Statistical Association. 1990;85:633–639.
- Sievers GL. A weighted dispersion function for estimation in linear models. Communications in Statistics, Theory and Methods. 1983;12:11611179.
- Terpstra J, McKean J. Rank-Based Analyses of Linear Models using R. Journal of Statistical Software. 2005;14(7)
- Tibshirani R. Regression shrinkage and selection via the LASSO. Journal of the Royal Statistical Society, Series B. 1996;58:267–288.
- Wang H, Li G, Jiang G. Robust regression shrinkage and consistent variable selection via the LAD-LASSO. Journal of Business & Economics Statistics. 2007;25:347–355.
- Wang H, Li R, Tsai CL. Tuning parameter selector for SCAD. Biometrika. 2007;94:553–568. [PMC free article] [PubMed]
- Wang L, Li R. Technical report # 665. School of Statistics, University of Minnesota; 2007. Weighted Wilcoxon-type smoothly clipped absolute deviation method.
- Wang YG, Zhao YD. Weighted rank regression for clustered data analysis. Biometrics. 2008;64:3945. [PubMed]
- White H. Consequences and detection of misspecified nonlinear regression models. Journal of the American Statistical Association. 1981;76:419–433.
- Wisnowski JW, Simpson JR, Montgomery DC, Runger GC. Resampling methods for variable selection in robust regression. Computational Statistics and Data Analysis. 2003;43:341–355.
- Zou H, Li R. One-step sparse estimates in nonconcave penalized likelihood models. The Annals of Statistics. 2007 to appear. [PMC free article] [PubMed]
- Zou H, Yuan M. Composite quantile regression and the oracle model selection theory. The Annals of Statistics. 2007 to appear.

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |