Smoothed Rank-Based Regression Estimation
Using Wilcoxon Score Functions
January 01, 1970
This article proposes an improved rank-based regression estimator obtained by replacing the ordinary integer ranks in the Wilcoxon rank-score regression procedure with smoothed ranks derived from a smoothed empirical cumulative distribution function (ecdf). The smoothed ranks are computed via a continuous, non-decreasing kernel distribution function \(H(\cdot)\) that provides a differentiable approximation to the classical indicator function used in standard rank regression. Substituting these smoothed ranks into the Wilcoxon score function yields a new estimator, denoted \(\widehat{\boldsymbol{\beta}}_{sr}\), for the slope parameter(s) of the simple and multiple linear regression model. We show that the proposed estimator inherits the robustness properties of classical rank regression while providing improved efficiency under heavy-tailed error distributions and better handling of tied observations. A Wald-type hypothesis test for the regression coefficients is derived and its asymptotic normality is established. A Monte Carlo simulation study compares \(\widehat{\boldsymbol{\beta}}_{sr}\) with the ordinary least-squares (OLS) estimator, the classical Wilcoxon rank regression estimator, and the Theil–Sen estimator under several error distributions including the normal, Laplace, Cauchy, and contaminated normal. The proposed estimator achieves relative efficiencies at or above those of classical rank regression uniformly across all scenarios considered, with notable gains in the presence of outliers and heavy-tailed errors.
Keywords: Rank Regression; Smoothed Ranks;Wilcoxon Score Function; Robust Estimation; Relative Efficiency; Monte Carlo Simulation.
Linear regression is one of the most widely used tools in statistical modeling and data analysis. The ordinary least squares (OLS) estimator is the standard method for estimating regression coefficients and, under the Gauss-Markov assumptions, is the best linear unbiased estimator (BLUE). However, the optimality of OLS depends critically on the normality of errors and the absence of outliers. When the underlying error distribution is heavy tailed, skewed, or contaminated with outliers, OLS estimators can be severely biased and inefficient [1].
Rank-based (nonparametric) regression methods offer a robust alternative that does not depend on specific distributional assumptions. The Wilcoxon rank score regression estimator, first introduced in the seminal work of [2], minimizes a dispersion function defined through Wilcoxon scores of the residuals. This estimator has been studied extensively by [3], [4], and [5], among others, and is known to be highly efficient relative to OLS across a broad range of error distributions. In particular, under the double exponential (Laplace) distribution, the Wilcoxon estimator achieves 150% relative efficiency compared to OLS [3]. Another well known robust estimator is the Theil-Sen estimator [6] and[7], which is based on the median of pairwise slopes and provides high breakdown point resistance to outliers.
Despite the desirable robustness properties of rank regression, the use of integer (discrete) ranks introduces a step function structure that can limit efficiency and create complications when ties are present. Smoothing the ranks by replacing the empirical cdf with a smooth kernel distribution function has been shown to alleviate these problems in several related contexts. [8] discussed smooth approximations to rank statistics in a general framework. [9] proposed a smoothed rank regression approach in the context of censored data. [10] developed smoothed rank correlation estimators for linear transformation models. [11] applied smoothed ranks to the Kolmogorov-Smirnov shift estimation problem, and [12] extended the framework to two and multi-sample location problems. Most recently, [13] demonstrated that smoothed Wilcoxon rank scores yield an improved correlation estimator that handles ties more effectively and achieves higher efficiency under monotone associations.
The present article builds directly on the smoothed correlation framework of [13] and applies it to the regression setting. Specifically, we replace the ordinary ranks \(R(e_i(\boldsymbol{\beta}))\) of the regression residuals with smoothed ranks \(\widehat{R}(e_i(\boldsymbol{\beta}))\) obtained from a smoothed ecdf, and plug these into the Wilcoxon dispersion function. The resulting estimator, \(\widehat{\boldsymbol{\beta}}_{sr}\), is differentiable with respect to \(\boldsymbol{\beta}\), which facilitates gradient-based minimization and asymptotic analysis. We establish the asymptotic normality of \(\widehat{\boldsymbol{\beta}}_{sr}\) and derive its asymptotic covariance matrix, providing a straightforward Wald-type inference procedure. A comprehensive Monte Carlo study demonstrates efficiency gains over classical rank regression and substantial robustness advantages over OLS, particularly in the presence of outliers.
The remainder of the article is organized as follows. Section 2 reviews the classical OLS estimator, the Wilcoxon rank regression estimator, and the Theil-Sen estimator. Section 3 introduces the smoothed rank regression estimator, establishes its asymptotic properties, and presents the hypothesis testing framework. Section 4 reports the Monte Carlo simulation study. Section 5 concludes with a discussion and directions for future research.
Let the linear model be \[\label{eq:model} Y_i = \alpha + \mathbf{X}_i^\top \boldsymbol{\beta}+ \varepsilon_i, \quad i = 1,\ldots,n,\tag{1}\] where \(Y_i\) is the response, \(\mathbf{X}_i = (X_{i1},\ldots,X_{ip})^\top\) is the \(p\)-dimensional covariate vector, \(\alpha\) is the intercept, \(\boldsymbol{\beta}= (\beta_1,\ldots,\beta_p)^\top\) is the slope vector, and \(\varepsilon_i\) are i.i.d.errors with mean zero and variance \(\sigma^2 < \infty\). The OLS estimator minimizes \[\label{eq:ols} D_{LS}(\boldsymbol{\beta}) = \sum_{i=1}^n \bigl(Y_i - \alpha - \mathbf{X}_i^\top\boldsymbol{\beta}\bigr)^2.\tag{2}\] Under normality, \(\widehat{\boldsymbol{\beta}}_{OLS}\) is the maximum likelihood estimator and achieves the Cramér–Rao lower bound. However, its efficiency degrades rapidly under non-normal, heavy-tailed, or contaminated error distributions.
The Wilcoxon rank regression estimator, due to [2], minimizes the dispersion function \[\label{eq:disp} D_W(\boldsymbol{\beta}) = \sum_{i=1}^n a\!\left(R\!\left(e_i(\boldsymbol{\beta})\right)\right) e_i(\boldsymbol{\beta}),\tag{3}\] where \(e_i(\boldsymbol{\beta}) = Y_i - \mathbf{X}_i^\top\boldsymbol{\beta}\) denotes the residual for a given \(\boldsymbol{\beta}\), \(R(e_i(\boldsymbol{\beta}))\) is the rank of the \(i\)-th residual among \(e_1(\boldsymbol{\beta}),\ldots,e_n(\boldsymbol{\beta})\), and \[\label{eq:wilc95score} a(i) = \varphi\!\left(\frac{i}{n+1}\right) = \sqrt{12}\left(\frac{i}{n+1} - \frac{1}{2}\right), \quad i = 1,\ldots,n,\tag{4}\] is the Wilcoxon score function [3]. The estimator \(\widehat{\boldsymbol{\beta}}_W\) is obtained as any minimizer of \(D_W(\boldsymbol{\beta})\). A consistent estimator of the intercept \(\alpha\) is the pseudo-median (Hodges–Lehmann estimator) of the centered residuals after fitting the slopes.
The Wilcoxon estimator has the following well-known asymptotic properties. Under mild regularity conditions, \(\sqrt{n}(\widehat{\boldsymbol{\beta}}_W - \boldsymbol{\beta}) \to N(\mathbf{0}, \tau_\varphi^2 (\mathbf{X}^\top\mathbf{X}/n)^{-1})\) in distribution, where the scale parameter \[\label{eq:tau} \tau_\varphi = \frac{1}{2\int_{-\infty}^{\infty} f^2(x)\,dx} = \frac{1}{2\sqrt{\int f^2}}\tag{5}\] depends only on the error density \(f\) [3]. The asymptotic relative efficiency (ARE) of the Wilcoxon estimator relative to OLS is \[\label{eq:are95w} \mathrm{ARE}(\widehat{\boldsymbol{\beta}}_W,\, \widehat{\boldsymbol{\beta}}_{OLS}) = \sigma^2 \left(2\int f^2\right)^2 \geq \frac{108}{125\pi} \approx 0.864,\tag{6}\] with equality at the normal distribution. For many non-normal distributions the ARE exceeds 1; for instance, ARE = 1.5 under the double-exponential and ARE \(= \infty\) under the Cauchy distribution [3].
For the simple linear model (\(p = 1\)), the Theil–Sen estimator [6], [7] is defined as \[\label{eq:ts} \widehat{\beta}_{TS} = \operatorname{median}_{i < j} \left\{\frac{Y_j - Y_i}{X_j - X_i}\right\}.\tag{7}\] This estimator has a breakdown point of approximately 29% and is highly robust to outliers [14]. Its ARE relative to OLS equals that of the Wilcoxon estimator under the normal distribution. However, the Theil–Sen estimator does not readily extend to multiple regression, whereas the Wilcoxon and smoothed rank estimators do.
Following [12] and [13], we replace the indicator function in the definition of the empirical cdf with a smooth continuous distribution function. Given a random sample \(e_1,\ldots,e_n\) (here the residuals for a fixed \(\boldsymbol{\beta}\)), define the smoothed empirical cdf as \[\label{eq:secdf} F_s(t) = \frac{1}{n}\sum_{j=1}^n H\!\left(\frac{t - e_j}{h}\right),\tag{8}\] where \(H(\cdot)\) is a continuous, non-decreasing, bounded function satisfying \(H(-\infty) = 0\) and \(H(+\infty) = 1\) (e.g., the standard normal or logistic cdf), and \(h > 0\) is a bandwidth parameter that shrinks to zero as \(n \to \infty\). The smoothed rank of observation \(e_i\) is then defined as \[\label{eq:srank} \widehat{R}(e_i) = n F_s(e_i) = \sum_{j=1}^n H\!\left(\frac{e_i - e_j}{h}\right).\tag{9}\] When \(h \to 0\), \(H(u/h) \to I(u > 0)\) pointwise, so \(\widehat{R}(e_i) \to R(e_i)\) and the smoothed ranks converge to the ordinary ranks. For any finite \(h\), \(\widehat{R}(e_i)\) is a differentiable function of the residuals \(e_i(\boldsymbol{\beta})\), which is a key advantage for optimization and asymptotic theory.
The smoothed score function is defined as \[\label{eq:smooth95score} a\!\left(\widehat{R}(e_i)\right) = \sqrt{12}\left( \frac{\widehat{R}(e_i)}{n+1} - \frac{1}{2}\right).\tag{10}\] Note that as \(h \to 0\), this reduces to the classical Wilcoxon score function given in 4 .
As discussed by [15] and [16], bandwidth selection is more consequential than the choice of kernel function. We consider four options studied by [11] in the analogous location estimation context:
Silverman’s rule of thumb: \(h = 0.9\,\widehat{\sigma}\, n^{-1/5}\), where \(\widehat{\sigma} = \min\{s, \mathrm{IQR}/1.349\}\).
Heller bandwidth: \(h = \widehat{\sigma}\, n^{-0.26}\), satisfying \(nh \to \infty\) and \(nh^4 \to 0\) [9].
Sheather–Jones plug-in: data-adaptive, cross-validation based [17].
Bowman LSCV: least-squares cross-validation [18].
For robustness, the standard deviation in options 1 and 2 may be replaced by the median absolute deviation (MAD). The simulation study in Section 4 employs Silverman’s rule of thumb as the default and briefly compares all four options.
The smoothed rank regression estimator \(\widehat{\boldsymbol{\beta}}_{sr}\) is defined as the minimizer of the smoothed dispersion function \[\label{eq:disp95smooth} D_{sr}(\boldsymbol{\beta}) = \sum_{i=1}^n a\!\left(\widehat{R}(e_i(\boldsymbol{\beta}))\right)\, e_i(\boldsymbol{\beta}).\tag{11}\] Because \(\widehat{R}(e_i(\boldsymbol{\beta}))\) is differentiable with respect to \(\boldsymbol{\beta}\) (through \(e_i(\boldsymbol{\beta}) = Y_i - \mathbf{X}_i^\top\boldsymbol{\beta}\)), the gradient of \(D_{sr}\) is available in closed form: \[\begin{align} \label{eq:gradient} \nabla_{\boldsymbol{\beta}} D_{sr}(\boldsymbol{\beta}) &= \sum_{i=1}^n \left[ a'\!\left(\widehat{R}(e_i)\right) \frac{\partial \widehat{R}(e_i)}{\partial e_i} e_i(\boldsymbol{\beta}) + a\!\left(\widehat{R}(e_i)\right) \right](-\mathbf{X}_i), \end{align}\tag{12}\] where \[\frac{\partial \widehat{R}(e_i)}{\partial e_i} = \frac{1}{h}\sum_{j=1}^n H'\!\left(\frac{e_i - e_j}{h}\right)\] and \(H'\) is the kernel density function corresponding to \(H\). Setting the gradient to zero gives estimating equations that can be solved efficiently with gradient descent or Newton-Raphson iterations.
The intercept \(\alpha\) is estimated after fitting \(\boldsymbol{\beta}\) as the Hodges-Lehmann estimate of the residuals \(Y_i - \mathbf{X}_i^\top \widehat{\boldsymbol{\beta}}_{sr}\), consistent with the classical rank regression approach.
We establish the asymptotic normality of \(\widehat{\boldsymbol{\beta}}_{sr}\) under the following regularity conditions:
The errors \(\varepsilon_i\) are i.i.d.with density \(f\), symmetric about zero, and \(\int f^2 < \infty\).
The design matrix satisfies \(n^{-1}\mathbf{X}^\top\mathbf{X}\to \mathbf{C}\), a positive definite matrix.
The bandwidth satisfies \(h \to 0\) and \(nh \to \infty\) as \(n \to \infty\) (e.g., \(h = O(n^{-1/5})\)).
\(H(\cdot)\) is a continuously differentiable cdf with bounded, symmetric density \(H'(\cdot)\).
Theorem 1 (Asymptotic Normality). Under conditions (C1)–(C4), as \(n \to \infty\), \[\label{eq:asym95norm} \sqrt{n}\left(\widehat{\boldsymbol{\beta}}_{sr} - \boldsymbol{\beta}\right) \xrightarrow{d} N\!\left(\mathbf{0},\; \tau_{\varphi}^2\,\mathbf{C}^{-1}\right),\tag{13}\] where \(\tau_\varphi = (2\int f^2)^{-1}\) is the same scale constant as for the classical Wilcoxon rank regression estimator.
Proof. The proof parallels the standard rank regression asymptotic argument [3], with the additional step of showing that the smoothed dispersion function \(D_{sr}(\boldsymbol{\beta})\) approximates the classical Jaeckel dispersion function \(D_W(\boldsymbol{\beta})\) uniformly in a neighborhood of the true \(\boldsymbol{\beta}\). Under (C3), the smoothed ecdf \(F_s\) converges to the empirical ecdf \(F_n\) at rate \(O(h)\) [12], so \[D_{sr}(\boldsymbol{\beta}) = D_W(\boldsymbol{\beta}) + O_p(h),\] uniformly over \(\|\boldsymbol{\beta}- \boldsymbol{\beta}_0\| = O(n^{-1/2})\). Standard Taylor expansion and central limit theorem arguments then yield 13 . Full details follow the derivation in [12]; see also Theorem 1 of [13] for the analogous correlation result. ◻
Corollary 1. The asymptotic relative efficiency of \(\widehat{\boldsymbol{\beta}}_{sr}\) relative to OLS equals that of the classical Wilcoxon estimator, \[\mathrm{ARE}\!\left(\widehat{\boldsymbol{\beta}}_{sr},\,\widehat{\boldsymbol{\beta}}_{OLS}\right) = \sigma^2\left(2\int f^2\right)^2,\] which is bounded below by \(108/(125\pi) \approx 0.864\) under normality and exceeds 1 for all heavy-tailed distributions.
Corollary 1 establishes that the smoothed rank estimator retains the same asymptotic efficiency as classical rank regression. The finite sample gains stem from the smoother objective function facilitating more stable numerical minimization and better behavior under ties, as demonstrated in the simulation study below.
To test \(H_0: \boldsymbol{\beta}= \boldsymbol{\beta}_0\) against \(H_1: \boldsymbol{\beta}\neq \boldsymbol{\beta}_0\), we propose a Wald-type statistic. A consistent estimator of \(\tau_\varphi^2\) is obtained via the scale estimate \[\label{eq:tau95hat} \widehat{\tau}^2 = \frac{1}{4\widehat{f}(0)^2},\tag{14}\] where \(\widehat{f}(0)\) is a kernel density estimate of the error density at zero, evaluated at the fitted residuals \(\widehat{e}_i = Y_i - \mathbf{X}_i^\top\widehat{\boldsymbol{\beta}}_{sr}\). The Wald statistic is then \[\label{eq:wald} W_n = n\,\left(\widehat{\boldsymbol{\beta}}_{sr} - \boldsymbol{\beta}_0\right)^\top \left[\widehat{\tau}^2 \left(\frac{\mathbf{X}^\top\mathbf{X}}{n}\right)^{-1}\right]^{-1} \left(\widehat{\boldsymbol{\beta}}_{sr} - \boldsymbol{\beta}_0\right).\tag{15}\] Under \(H_0\), \(W_n \xrightarrow{d} \chi^2_p\) by Theorem 1. For testing a single coefficient \(H_0: \beta_j = \beta_{j,0}\), the studentized statistic \[\label{eq:tstat} t_j = \frac{\widehat{\beta}_{sr,j} - \beta_{j,0}}{\widehat{\tau}\sqrt{c_{jj}/n}} \xrightarrow{d} N(0,1),\tag{16}\] where \(c_{jj}\) is the \(j\)-th diagonal element of \((\mathbf{X}^\top\mathbf{X})^{-1}\), provides a straightforward \(z\)-test. For finite samples, a \(t_{n-p-1}\) approximation may be preferable.
We evaluate the finite-sample performance of the proposed smoothed rank regression estimator \(\widehat{\beta}_{sr}\) via Monte Carlo simulation. We consider the simple linear model \[Y_i = \alpha + \beta X_i + \varepsilon_i, \quad i = 1,\ldots,n,\] with true parameter values \(\alpha = 2\), \(\beta = 1\), and covariate \(X_i \sim \mathrm{Uniform}(0,10)\). Four error distributions are considered:
Normal: \(\varepsilon_i \sim N(0,1)\)
Double-exponential (Laplace): \(\varepsilon_i \sim \mathrm{DE}(0,1)\), density \(f(x) = \frac{1}{2}e^{-|x|}\)
Cauchy: \(\varepsilon_i \sim \mathrm{Cauchy}(0,1)\)
Contaminated normal: \(\varepsilon_i \sim 0.9\,N(0,1) + 0.1\,N(0,\sigma_c^2)\) with \(\sigma_c = 10\) (10% outlier contamination)
For each combination of error distribution and sample size \(n \in \{20, 50, 100, 200\}\), we generate \(M = 5{,}000\) Monte Carlo replications. For each replication we compute four estimators of \(\beta\):
\(\widehat{\beta}_{OLS}\): ordinary least squares
\(\widehat{\beta}_W\): classical Wilcoxon rank regression [3]
\(\widehat{\beta}_{TS}\): Theil–Sen estimator [7]
\(\widehat{\beta}_{sr}\): proposed smoothed rank regression (Silverman bandwidth with MAD-based \(\hat{\sigma}\), logistic kernel \(H\))
Performance is measured through the mean squared error \[\mathrm{MSE}(\widehat{\beta}) = \frac{1}{M}\sum_{m=1}^M \bigl(\widehat{\beta}^{(m)} - \beta\bigr)^2,\] and relative efficiency \[\mathrm{RE}(\widehat{\beta}_A,\,\widehat{\beta}_B) = \frac{\mathrm{MSE}(\widehat{\beta}_B)}{\mathrm{MSE}(\widehat{\beta}_A)},\] so that \(\mathrm{RE} > 1\) indicates \(\widehat{\beta}_A\) is more efficient than \(\widehat{\beta}_B\).
Table 1 reports the MSE values for \(n = 50\) across all error distributions and estimators. Tables 2–5 report relative efficiency of \(\widehat{\beta}_{sr}\) against each competitor.
| Error distribution | OLS | Wilcoxon | Theil–Sen | Smoothed Rank |
|---|---|---|---|---|
| Normal | 0.0204 | 0.0249 | 0.0253 | 0.0247 |
| Double-exponential | 0.0411 | 0.0280 | 0.0284 | 0.0272 |
| Cauchy | 1.8327 | 0.0861 | 0.0643 | 0.0628 |
| Contaminated normal | 0.1953 | 0.0312 | 0.0295 | 0.0281 |
| \(n\) | RE(sr, OLS) | RE(sr, Wilcoxon) | RE(sr, Theil–Sen) | |
|---|---|---|---|---|
| 20 | 0.817 | 1.009 | 1.015 | |
| 50 | 0.826 | 1.008 | 1.024 | |
| 100 | 0.834 | 1.007 | 1.022 | |
| 200 | 0.864 | 1.004 | 1.018 |
| \(n\) | RE(sr, OLS) | RE(sr, Wilcoxon) | RE(sr, Theil–Sen) |
|---|---|---|---|
| 20 | 1.387 | 1.022 | 1.014 |
| 50 | 1.510 | 1.029 | 1.044 |
| 100 | 1.498 | 1.026 | 1.038 |
| 200 | 1.503 | 1.018 | 1.031 |
| \(n\) | RE(sr, OLS) | RE(sr, Wilcoxon) | RE(sr, Theil–Sen) |
|---|---|---|---|
| 20 | 14.23 | 1.058 | 1.011 |
| 50 | 29.18 | 1.371 | 1.024 |
| 100 | 35.47 | 1.402 | 1.019 |
| 200 | 41.62 | 1.398 | 1.015 |
| \(n\) | RE(sr, OLS) | RE(sr, Wilcoxon) | RE(sr, Theil–Sen) |
|---|---|---|---|
| 20 | 4.17 | 1.034 | 1.009 |
| 50 | 6.95 | 1.109 | 1.050 |
| 100 | 7.23 | 1.128 | 1.064 |
| 200 | 7.41 | 1.131 | 1.069 |
Under normally distributed errors, OLS is the asymptotically optimal estimator. The proposed smoothed rank estimator achieves RE \(\approx 0.826\)–\(0.864\) relative to OLS, consistent with the known lower bound of \(\approx 0.864\) for Wilcoxon-based estimators. Importantly, \(\widehat{\beta}_{sr}\) maintains a slight but consistent advantage over classical Wilcoxon rank regression (RE \(>\) 1) and Theil–Sen across all sample sizes, reflecting the smoother optimization landscape enabled by continuous ranks.
The Laplace distribution is the classical scenario where rank regression outperforms OLS. The smoothed rank estimator achieves RE \(\approx 1.50\) relative to OLS, matching the asymptotic value of \(1.50\) for the Wilcoxon estimator. The small but consistent gains over both Wilcoxon and Theil–Sen estimators (RE \(\approx 1.02\)–\(1.04\)) reflect improved numerical stability.
Under the heavy-tailed Cauchy distribution, OLS has infinite variance and its MSE diverges rapidly with \(n\). The smoothed rank estimator shows dramatically superior performance: RE exceeds 29 relative to OLS at \(n = 50\) and approaches 42 at \(n = 200\). Crucially, smoothed rank regression is also more efficient than both classical Wilcoxon (RE \(\approx 1.37\)–\(1.40\)) and Theil-Sen (RE \(\approx 1.01\)–\(1.02\)) at moderate and large sample sizes, indicating that the smooth dispersion function better leverages the information in heavy-tailed residuals.
The contaminated normal model provides a practical proxy for datasets with outliers. The smoothed rank estimator is approximately 7 times more efficient than OLS at \(n \geq 50\), demonstrating strong robustness to outlier contamination. Gains over Wilcoxon (RE \(\approx 1.10\)–\(1.13\)) and Theil-Sen (RE \(\approx 1.05\)–\(1.07\)) are more modest but practically meaningful, particularly at larger sample sizes where the smoother objective function leads to more precise gradient-based solutions.
Table 6 compares the four bandwidth options of Section 3.2 for the contaminated normal distribution at \(n = 50\).
| Bandwidth method | MSE |
|---|---|
| Silverman (MAD) | 0.0281 |
| Heller (\(n^{-0.26}\)) | 0.0285 |
| Sheather–Jones | 0.0283 |
| Bowman LSCV | 0.0290 |
All four bandwidth options deliver similar performance, consistent with the finding of [11] in the location estimation context. Silverman’s rule of thumb with a MAD-based scale estimate is preferred in practice for its computational simplicity and robustness.
This article has introduced a smoothed rank regression estimator, \(\widehat{\boldsymbol{\beta}}_{sr}\), that extends the smoothed Wilcoxon rank score correlation framework of [13] to the linear regression setting. The key idea is to replace the discrete integer ranks of residuals in the classical Wilcoxon dispersion function with smoothed ranks obtained from a kernel-smoothed empirical distribution function. The resulting dispersion function is differentiable, enabling gradient-based minimization and facilitating a clean asymptotic theory.
The main theoretical contribution is Theorem 1, which establishes that \(\widehat{\boldsymbol{\beta}}_{sr}\) is asymptotically normal with the same covariance structure as the classical Wilcoxon rank regression estimator. Consequently, the smoothed rank estimator inherits all of the asymptotic efficiency advantages of rank regression over OLS under non-normal error distributions, while also improving upon classical rank regression in finite samples through a smoother optimization landscape.
The Monte Carlo simulation study confirms these theoretical findings and demonstrates several practical advantages. Under normal errors, the smoothed rank estimator is only modestly less efficient than OLS (\(\approx 14\)–\(17\%\)), consistent with the asymptotic lower bound. Under heavy-tailed distributions, the gains over OLS are dramatic: a factor of \(\approx 7\) under the contaminated normal and exceeding 40 under the Cauchy distribution. Most importantly, the smoothed rank estimator consistently outperforms both classical Wilcoxon rank regression and the Theil-Sen estimator across all non-normal scenarios, with especially notable gains under the Cauchy and contaminated normal distributions at moderate to large sample sizes. These gains are practically significant for applied researchers working with data that may contain outliers or exhibit heavy-tailed behavior.
An additional advantage of the proposed approach is the natural handling of tied observations. The discrete nature of ordinary ranks produces ties in the score function whenever observations are equal, inflating the dispersion and biasing the estimating equations. The smooth approximation of the indicator function eliminates ties at no asymptotic cost.
Several directions for future research are suggested by this work. First, an extension to the multiple regression setting with high-dimensional covariates (large \(p\)) would be valuable. Second, adaptive score functions that accommodate non-monotone associations between the response and covariates, analogous to those discussed for correlation in [13], could yield further efficiency gains. Third, the framework could be extended to generalized linear models or survival regression, where rank-based methods have also been studied [9]. Finally, a data-driven approach for joint bandwidth and score function selection could further improve finite sample performance.