Robust Three-Step Regression Based on Comedian and Its Performance in Cell-Wise and Case-Wise Outliers

Both cell-wise and case-wise outliers may appear in a real data set at the same time. Few methods have been developed in order to deal with both types of outliers when formulating a regression model. In this work, a robust estimator is proposed based on a three-step method named 3S-regression, which uses the comedian as a highly robust scatter estimate. An intensive simulation study is conducted in order to evaluate the performance of the proposed comedian 3S-regression estimator in the presence of cell-wise and case-wise outliers. In addition, a comparison of this estimator with recently developed robust methods is carried out. The proposed method is also extended to the model with continuous and dummy covariates. Finally, a real data set is analyzed for illustration in order to show potential applications.


Introduction
Regression models are one of most used statistical tools for diverse practitioners [1]. A well-known assumption to infer the parameters of these models is that the error term follows a normal distribution. However, as it was discussed in [2], the presence of outliers could affect the estimation of parameters under normality and lead to inaccurate results [3]. Because outliers frequently are present in real data sets, robust estimation methods are called to be considered in practice to avoid this inaccuracy [4][5][6][7].
A type of model introduced in [8] is known as three-step (3S) regression, which proceed as follows. In the first step, a univariate filter is applied to remove cell-wise outliers. In the second step, a generalized S-estimator (GSE) is used to down-weight the effect of case-wise outliers. In the third step, the regression coefficients are estimated. Nevertheless, some limitations of the 3S-regression method are mentioned in [9]. For instance, the GSE employed in the second step loses robustness against case-wise outliers, when the dimension is greater than ten. In addition, the extended minimum volume ellipsoid (EMVE) estimator, also utilized by the GSE as an initial value, is computationally slow and it does not scale well to higher dimensions. Hence, a new robust estimator, called generalized Rocke S-estimator (GRE), is proposed in [9] to replace the GSE in the second step. Furthermore, a cluster-based algorithm introduced in [9] for faster and more reliable sampling when computing the EMVE estimator, named EMVE-C, is used as an initial value for the GRE algorithm. To our best knowledge, the comedian [10,11] has not been employed as an initial estimate in 3S-regression methods.
In robust linear regression, the case-wise contamination model is known as Tuckey-Huber contamination (THC). This model has been broadly studied in the literature, but its use in practice is not frequent. In the THC model, a small proportion of cases are contaminated. When the contamination is carried out using cell-wise, the independent contamination (IC) model arises. Note that the IC model has a small proportion of individual cells in the covariates that are independently contaminated [12]. However, the literature about the IC model is limited. To our best knowledge, there is few works that can deal with both types of outliers (case-wise and cell-wise) at the same time.
A classical robust regression for case-wise contamination is the least median of square (LMS) method. The LMS regression was proposed in order to optimize the median of the squares of residuals and it allows for a breakdown point of 50%, but it is computationally inefficient [13,14]. Two alternative methods that use iterative strategies are the regressions via the estimation of: (i) the minimal covariance determinant (MCD) [15,16]; and (ii) the iterative re-weighted least square (IRLS) [17,18]. The MCD regression minimizes the covariance matrix determinant of the central points. The IRLS regression includes additional information regarding error variance and covariance by incorporating a weight matrix into the model estimation, whose diagonal elements depend on a loss function. High breakdown affine equivariant estimators provide down-weighting to outlying cases, such as the least trimmed square (LTS) regression [14], S-regression [19], and MM-regression [20]; see also [21,22] for M-estimators in regression. All of these methods work well in practice under the THC model.
A number of authors [8,[23][24][25][26][27][28] have proposed robust regression models that are resilient to case-wise and cell-wise outliers by robustifying the components of the covariance matrix in the solution of the least square (LS) optimization problem. Additionally, the multivariate S-estimator is incorporated instead of the empirical covariance and mean [24,25]. It has been also showed that, under mild assumptions (including symmetry and independence in the residuals), the 2S-regression estimator [8] is Fisher consistent and asymptotically normal, even if the multivariate S-estimators are not. Based on incomplete data, two kinds of estimators were constructed [26]: (i) the GSE and (ii) the extended S-estimator (ESE), which match with the multivariate S-estimator under complete data. Note that the GSE needs a robust initial estimate. Furthermore, the extended EMVE estimator is introduced in [26] as a particular case of the ESE. The EMVE estimator can be considered to be an initial value and a generalization of the minimum volume ellipsoid (MVE) estimator proposed in [15]. Moreover, the shooting S-estimator that is derived in [28] assigns individual weights to each cell in the data table, combining the shooting algorithm [29] and the simple S-regression [19]. Observe that the data may be snipped replacing cell-wise outliers by missing values NA [27]. Moreover, the Gervini-Yohai univariate filter [30] can be used followed by the GSE [26,31]. Notice that the 3S-regression [8] considers an estimator which is analogous to one defined in [26], but with the filter that is consistent for a broader range of distributions.
Based on this bibliographical review, the objective of this study is to propose a comedian-three-step (C3S) regression estimator that considers the comedian as the initial robust scatter value for the GRE algorithm. The proposed estimation method: (i) utilizes an adaptive consistent univariate filter to control the effect of extreme cell-wise outlier propagation; (ii) applies the GRE algorithm, but modified using the sample comedian matrix and wise-median as an initial robust scatter estimate for the filtered data; and (iii) estimates the regression coefficients using the GRE algorithm in the previous step.
This paper is organized, as follows. In Section 2, the general context and notations are provided. Then, the consistent factor in median absolute deviation (MAD), the comedian function, and the empirical comedian covariance matrix are defined. Section 3 introduces the models with continuous and dummy covariates and proposes an estimator based on the C3S-regression, as well as its asymptotic properties. In Section 4, an extensive simulation study is conducted in order to compare the performance of proposed estimator with recently developed robust methods. Additionally, in this section, a real data example is used for illustration and for showing potential applications. Section 5 describes some conclusions and ideas for possible future works.

Comedian Covariance Matrix and Comedian Matrix
In this section, the general context and notations used in this work are presented. The consistent factor in MAD, the comedian function, and the empirical comedian covariance matrix are also defined here.

General Context and Notations
Let X and Y be two continuous random variables. Subsequently, the MAD of X, the comedian between X and Y, and the correlation median between X and Y (δ) are, respectively, given as [10] MAD(X) = median(|X − median(X)|), Note that the MAD(X) defined in Equation (1) is a robust measure of dispersion (or scatter) of X, COM(X, Y) is a robust measure of the covariance between X and Y, while δ is a robust measure of the correlation between X and Y. When X = Y, (COM(X, X)) 1/2 can be used as a robust measure of variability, such as it occurs with the standard deviation (SD) and the covariance of X, that is, SD(X) = (COV(X, X)) 1/2 . Then, a robust measure of the covariance and correlation matrices for any random vector may be obtained by utilizing the robust measures defined in (1). Notice that the usual covariance between X and Y can be obtained as where g is the comedian function stated as g( ) = COM(X, Y) (see Lemma 2.1 in [10]) and is the correlation coefficient between X and Y. Observe that may be represented in terms of the correlation median as = COR(X, Y) = g −1 (g(1)δ).
The comedian function under a non-degenerate bivariate normal distribution was studied in [10], obtaining g(1) = (Φ −1 (0.75)) 2 , where Φ is the standard normal cumulative distribution function and Φ −1 is the inverse of Φ or normal quantile function. In this work, we also extend g(1) to other non-normal distributions. Note that Equation (2) can be written as and the consistent factors b X and b Y depend upon the marginal distributions of X and Y, respectively. The consistent factors for some distributions have been obtained and are presented in Table 1. Detailed calculations of these factors are available upon request from the authors. Notice that the comedian can be considered as a robust initial scatter estimate for the GRE algorithm. Let (X 1 , Y 1 ), . . . , (X n , Y n ) be independent random vectors following a bivariate distribution. Subsequently, the empirical comedian is established by COM n (X, Y) = median i∈{1,...,n} (X i − median n (X))(Y i − median n (Y)), (5) where median n (X) and median n (Y) denote the sample medians of X 1 , . . . , X n and Y 1 , . . . , Y n , respectively. Thus, the covariance matrix that is defined in Equation (2) is a sophisticated scatter estimate based on the comedian matrix. Table 1. Consistent factor b X of the indicated distribution.

Distribution of X Notation
* m median is the solution of m to the non-linear equations for t(ν) and Wei(α, β) given, respectively, by

The Consistent Factor in MAD
The MAD is a very robust scatter estimate, which has 50% breakdown point (the best possible). According to [32], if we want to estimate the SD consistently, the MAD must be multiplied by a correction factor. Thereby, an alternative robust estimate of the SD of X is given by S X = b X MAD n (X), where MAD n (X) = median{|X i − median n (X)|, i ∈ {1, . . . , n}}. The consistent factor b X depends exclusively on the distribution of the random variable X. If the marginal distribution of X is unknown, the consistent factor can be estimated via the non-parametric bootstrapping method [33]. Let F n be the empirical cumulative distribution function of the random variable X. Subsequently, the bootstrapping process used to obtain the consistent factor b X is summarized in Algorithm 1.

Algorithm 1
Bootstrapping process used in order to obtain the consistent factor b X . 1: Generate X * 1 , . . . , X * n ∼ F n (x) randomly.

The Comedian and Empirical Comedian Covariance
The comedian function g( ) is needed in order to estimate the correlation median defined in Equation (3). The empirical correlation median δ n = COM n (X, Y)( MAD n (X) MAD n (Y)) −1 stated in Equation (5) can be seen as a robust estimate of the correlation coefficient by The function g was analyzed in [10] when (X, Y) has a bivariate normal distribution, but an explicit form was not obtained. However, it may be approximated through Monte Carlo simulations. We conduct an extensive Monte Carlo simulation study for g via the R software [34]. This simulation was carried out by using an R package named MASS and its mvrnorm function. The empirical medians of 10 000 000 random numbers from a bivariate normal distribution were also used, with varying from −1 to 1 by 0.01 when / ∈ [−0.1, 0.1], and by 0.001 when ∈ [−0.1, 0.1]. The number of replicates is N = 10 and a visualization of the approximation for g is shown in Figure 1. In general, the exact image value of the inverse comedian function g −1 cannot be obtained, but it may be estimated through an approximation. A discrete approximation of the comedian function is obtained at the aforementioned simulation study. The expression that is given in Equation (6) may be approximated as n = g −1 (g(1) δ n ), where g −1 is an estimate of g −1 by interpolating the approximated points while using a cubic spline. To carry this method, it is necessary to know all of the values corresponding to the support of the inverse comedian function. If the consistent factors of the marginal distributions defined in Equation (4) are estimated properly, then n for other bivariate distributions can be obtained.
By using the empirical marginal distributions, the consistent factors may be estimated via bootstrapping, as described in Section 2. Thus, from Equation (6) Let X = (X 1 , . . . , X p ) be a set of p covariates. Then, a robust version of the empirical covariance between any pair of covariates (X i , X j ), for i, j ∈ {1, . . . , p}, can be stated as where S c is an element of the robust version of the empirical covariance matrix S c XX . We also propose to use the expression given in Equation (7) as an initial scatter estimate for the GRE algorithm instead of the EMVE estimator. This proposal is called here the full version of the C3S-regression estimator.

Comedian Three-Step Regression
In this section, the model with continuous and dummy covariates is introduced. Moreover, the proposed estimator that is based on the 3S-regression, as well as its asymptotic properties, are developed.

The Proposed Estimator
A multiple regression is used to model the linear relationship between a dependent (response) variable Y and p independent (covariates) variables X = (X 1 , . . . , X p ) with observed values for the case i denoted by x i = (x i1 , . . . , x ip ) . Subsequently, the multiple regression model can be written as where the error terms ε i , for i ∈ {1, . . . , n}, are independent and identically distributed random variables, which are also independent of the values of the covariates The LS estimates of the parameters (β 0 , β) are defined as the solution to an optimization problem in order to minimize the sum squares of residuals as The solution of Equation (9) can be explicitly given by where Σ XX and Σ XY are the components of the empirical covariance matrix, and µ Y and µ X are empirical means of Y and X, respectively. As suggested by a number of authors [8,[23][24][25][26][27][28], the components of the solution stated in Equation (10) can be robustified to immunize the estimator against case-wise and cell-wise outliers. Inspired by [9], we use a modified version of the GRE algorithm to obtain robust estimates of means and covariances needed in the solution presented in Equation (10). The modification proposed by us is that the GRE algorithm considers the comedian as an initial value instead of the EMVE-C estimate. The robust method basically utilizes the empirical median and the robust version of the covariance, introduced in Section 2, as the initial location and scatter estimates for the GRE algorithm. The proposed estimator uses the univariate filter given in [8] and the GRE algorithm for incomplete data developed in [9]. Our proposal works similarly to that used in the 3S-regression, but employing in the second step a different initial robust estimate for the GRE algorithm. In the present work, the initial estimates of location and scatter, the empirical median, and a robust estimate of the covariance, are computed after snipping the data. Therefore, the proposed robust regression estimator (C3S-regression) is established as where both m and S come from the modified GRE algorithm proposed in this work, and they are computed as in Algorithm 2.
Algorithm 2 Computation of m and S from the modified GRE algorithm.
1: Filter extreme cell-wise outliers using a univariate filter to prevent cell-wise contaminated cases. 2: Compute the wise-median and the robust version of the covariance matrix (or comedian matrix) as initial robust location and scatter estimates. 3: Down-weight the effect of case-wise outliers by applying the GRE algorithm for computing robust location and scatter estimates with the filtered data from Step 1.
Now, consider a set of n data with observed covariates {x 1 , . . . , x n } and the corresponding response variable {Y 1 , . . . , Y n }. Let {z 1 , . . . , z n } be the joint data with z i = (Y i , x i ) . In the first step, a univariate filter, as described in [8], is applied to each observed covariate x j , for j ∈ {1, . . . , p}. Let Z = (z 1 , . . . , z n ) and U denote the resulting auxiliary matrices of zeros and ones, with zeros indicating the filtered (missing) entries. Subsequently, based on the GRE algorithm, we obtain where m GRE and S GRE are robust location and scatter based on the GRE algorithm for the incomplete data (Z, U). Computation of the 3S-regression and C3S-regression estimates is summarized Algorithm 3.

Algorithm 3 Computation of 3S-regression and C3S-regression estimates.
1: Snip data. 2: Apply the GRE algorithm with robust initial location and scatter estimates. 3: Estimate the regression coefficients as in Equation (10).
Note that the C3S-regression uses a robust estimated covariance matrix as an initial scatter value instead of the EMVE estimate. In addition, the biflat ρ function [35] is employed instead of the Tukey bisquare ρ function for the GSE. More details of the definitions and algorithms of the EMVE estimate and GSE can be found in [26]. Furthermore, the GRE algorithm was studied in [9], showing that, in large dimension, the Rocke biflat function is more robust than the Tuckey bisquare function [35,36].

Models with Continuous and Dummy Covariates
Notice that the M-regression and 3S-regression were used in [8] in order to deal with continuous and dummy covariates. There, a 3S-regression was employed to estimate the coefficients of the continuous covariates, whereas an M-regression with the Huber ρ function given by ρ H (t) = min(1, t 2 /2) [37] was considered to estimate the coefficients of the dummy covariates. This is a modification of the M-regression and 3S-regression proposed in [38]. We act similarly by using the C3S-regression to estimate the coefficients of the continuous covariates and the M-regression for the dummy coefficients. Consider the model with continuous and dummy covariates defined as where the columns in X and D are linearly independent. More precisely, our method that is based on the M-regression and C3S-regression works as where h denotes the operator of a C3S-regression in each iteration for (X, Y); while M denotes the operator of an M-regression with no intercept for (D, Y), as stated in (11). To control the effect of propagation of cell-wise outliers, let X be the imputed X with the filtered entries by the linear predictor using ( m (r) , S (r) ) as defined in Equation (12) at the rth iteration of the GRE algorithm. The method that is presented in Equation (14) needs initial estimates ( β 2 ) to start the algoritm until a maximum of R = 20 iterations [8]. Then, we first remove the effect of d i from the continuous covariates and response. Let Y = Y − D t and X = X − DT, where t = M(D, Y) and T is a p 1 × p 2 matrix with the jth column as T j = M(D, (x ij , . . . , x nj ) ). Subsequently, the initial estimates are defined by

Asymptotic Properties of the Comedian Three-Step Regression
The strong consistency of the empirical comedian is proved in [10], as well as its asymptotic normality. The strong consistency, asymptotic normality, and regularity assumptions for the GSE were established in [8]. Because the respective estimates under the 3S-regression and C3S-regression are based on the same GSE, independent of the differences in the initial estimates and weight functions, asymptotic properties of the corresponding estimators from the C3S-regression and 3S-regression are guaranteed. Note that the 3S-regression and C3S-regression become a 2S-regression for an enough large n. Therefore, the estimators obtained from the C3S-regression inherit the asymptotic properties of the estimators obtained from the 2S-regression such as the 3S-regression does. The properties of the corresponding asymptotic covariance matrix are also found in [8].
Let H be the distribution of (X , Y), ( m, S) be the GSE, and ( β 0 C3S , β C3S ) be the estimated 3S-regression coefficients. Subsequently, , where x i is the best linear prediction of x i . Thus, the asymptotic covariance matrix is estimated through the asymptotic S-estimator variance (ASV) matrix stated as ASV(H) = C3S , and w(d n ( z i )) = ρ R (d n ( z i )), with ρ R being the Rocke biflat function and d n defined in Equation (13).

Numerical Studies
In this section, the computational framework and simulation scenarios are described. Subsequently, we report the results of an intensive simulation study, which is conducted to evaluate the statistical performance of the C3S-regression coefficient estimators and to compare the proposed estimator and other existing estimators. In addition, the illustration with real data is provided.

Computational Framework and Simulation Scenarios
Our simulation study is performed by utilizing the R software with a Hewlett-Packard HP Compaq computer, Pro 6300 SFF with 8 cores processor GenuineIntel Intel(R) Core(TM) i7-3770 CPU @ 3.40 GHz. The simulation is similar to the one carried out in [8] using the same criteria for evaluating the performance of the C3S-regression. The method proposed in the present paper is also compared with the performance of the LS regression and following two robust alternatives: • The first one is the 2S-regression [25] that uses an MVE estimate as an initial value. The MVE estimate is computed by means of an iterative subsampling with a concentration step. The MVE estimate is implemented in an R package named rrcov, using the function CovSest with the option method = "bisquare" [39].

•
The second one is the 3S-regression [8] that reduces the high computational burden of uniform subsampling for the EMVE estimate. The GSE with bisquare ρ function is computed by an iterative algorithm that employs the EMVE-C estimate as an initial value. The 3S-regression without modifications is implemented in an R package, named robreg3S, using the function robreg3S as the default option [40]. Nevertheless, the GSE with the EMVE-C estimate as an initial value is implemented in the GSE package, using the function GSE with the option init = "emve_c" [41].
The univariate filter that is needed by the C3S-regression is implemented in the robreg3S package, while the GRE is computed using the GSE function with the option method = "rocke". Two versions of the C3S-regression to be considered are: (i) the full version using the comedian covariance matrix; and, (ii) the light version using the raw comedian matrix as an initial scatter estimate. The last one is called light version, because it uses less operations to compute than the full version. From now on, the C3S-regression is referred to both versions, unless that an indication is done.
Next, the regression model presented in Equation (8) with p = 15 and n ∈ {150, 300, 500, 1000, 5000} is considered. The values of covariates x i , for i ∈ {1, . . . , n}, are generated from a multivariate normal distribution N p (µ, Σ). We set µ = 0 and Σ jj = 1 [8], for j ∈ {1, . . . , p}, without loss of generality, because the GSE used by the C3S-regression is location and scale equivariant. (Note that from the location equivariant of the GSE, β 0 = 0 can be set.) To address the fact that the C3S-regression and 3S-regression are neither affine-equivariant nor regression equivariant, the correlation structure Σ may be used. Observe that this correlation structure is described in [27], with the condition number fixed at 100 and random generation of β as β = Rb. Let R = 10 and b follow a uniform distribution on the unit spherical surface. The response variable Y i is given by Y i = x i β + ε i , where ε i ∼ N(0, σ = 0.5), for i ∈ {1, . . . , n}, are independent. The scenarios assumed in the simulation study are: S1 Clean data: the data generation is not altered. S2 Cell-wise contamination: randomly replace a proportion q of the cells in the covariates by outliers x cont ij = E(X ij ) + k SD(X ij ) and of the responses by outliers Y cont Case-wise contamination: randomly replace a proportion q of the cases by leverage outliers Here, v is the eigenvector corresponding to the smallest eigenvalue of Σ with length such that (v − µ) Σ −1 (v − µ) = 1. To compute the value of c ∈ {1, . . . , 100}, we follow the same process introduced in [8,27]; that is, a Monte Carlo study with the same number of replicates N = 500. We observe that c = 22 is the value that produces the worst performance of the scatter estimator.
Let q ∈ {0.01, 0.05, 0.09} for the cell-wise contamination. From the fact that case-wise outliers are unusual in practice, we consider q = 0.03 for the case-wise contamination. The number of replicates for each setting is N = 1000. In addition, the simulation study is also carried out in order to consider the regression model presented by Equation (13) with p 1 = 12 continuous covariates, p 2 = 3 dummy covariates, and n ∈ {500, 1000}. Then, the performance of the M-regression and C3S-regression is evaluated. The values of covariates (x i , d i ), for i ∈ {1, . . . , n}, are first generated from a multivariate normal distribution N p 1 +p 2 (0, Σ), where Σ is the randomly generated correlation matrix with a fixed condition number of 100. Subsequently, d ij is dichotomized at Φ −1 (π j ), with π j ∈ {1/4, 1/3, 1/2}, for j ∈ {1, 2, 3}, respectively. The generation of the model with continuous and dummy covariates follows the scenarios S1-S2 and for the case-wise contamination follows the scenario S3.
Let Σ 1 be a sub-matrix of Σ, which quantifies the covariance of the continuous covariates. In this new scenario, randomly replace a proportion q of the cases in X by leverage outliers (x cont Here, v is now the eigenvector corresponding to the smallest eigenvalue of Σ 1 , with length such that (v − µ) Σ −1 1 (v − µ) = 1, and the corresponding least favorable case-wise contamination size for the twelve continuous variables is c = 18.
Once again, we consider q ∈ {0.01, 0.05, 0.09} for the cell-wise contamination and q = 0.03 for the case-wise contamination. The number of replicates for each setting is N = 1000. Furthermore, the simulation study is conducted for non-normal covariates to compare the performance of the C3S-regression, 3S-regression, 2S-regression and LS estimators. For the C3S-regression, the full and light versions of the proposed estimator are considered. The same regression model with p = 15 and n = 500 is used, but the covariates are generated from a non-normal distribution [8]. The covariates X i , for i ∈ {1, . . . , n}, are first generated from a multivariate normal distribution with zero mean and covariance matrix Σ, which is, X i ∼ N p (0, Σ), where, again, Σ is a randomly generated correlation matrix with a fix condition number of 100. Subsequently, the covariates are transformed by means of (X i1 , . . . , ). We consider a distribution for G j as: N(0, 1), with j ∈ {1, 2, 3}; χ 2 (20), with j ∈ {4, 5, 6}; F(90, 10), with j ∈ {7, 8, 9}; χ 2 (1), with j ∈ {10, 11, 12}; and Pareto(1, 3), with j ∈ {13, 14, 15}. The scenarios that are evaluated in this simulation study are as S1. For the cell-wise contamination, we replace q = 0.05 by the proportion of cells in the covariates with outliers x cont ij = kG j (0.999), and by the proportion of responses with outliers Y cont ij = E(Y ij ) + k SD(ε i ).

Simulation Results
The statistical performance in the estimation of regression coefficients due to the effect of cell-wise and case-wise outliers can be evaluated using the empirical mean squared error (MSE), defined as where β (i) j is the estimate of β (i) j at the ith Monte Carlo replicate. Tables 2 and 3 report the MSE defined in Equation (15) for k = 1 in all the settings with n ∈ {500, 1000}. The results for k = {5, 10} are omitted, because they are similar to k = 1. Figures 2 and 3 show curves of MSE for cell-wise and case-wise contamination in models with p = 15 continuous covariates and n = 1000.     Figure 6 shows curves of MSE for cell-wise and case-wise contamination in models with continuous covariates and n = 500. Note that models with continuous covariates in the M-regression and C3S-regression outperform in all of the assumed scenarios for both cell-wise and case-wise contaminations. In addition, in the four panels of Figure 6, both versions of the C3S-regression have almost the same behavior for all settings assumed. The full version of the C3S-regression is little more robust than its light version, but the estimates of both are almost equal for all contamination settings. The results for n = 1000 are similar to the cell-wise contamination settings. In the cell-wise contamination setting for small and moderate contamination proportions (q ≤ 0.05), the C3S-regression is highly robust against moderate and large cell-wise outliers (k ≥ 3), but less robust against inliers (k ≤ 2). The 3S-regression and C3S-regression perform similarly for moderate and large outliers, but in the presence of inliers (k ≤ 3), the 3S-regression is less robust; see first two panels of Figure 6.   The 2S-regression and 3S-regression perform similarly in the presence of inliers, as expected from the simulation studies carried out in [8]. However, the 2S-regression breaks down in cases when the proportion of contaminated cells is q > 0.5; that is, when the propagation of large cell-wise outliers is expected to affect more than 50% of the cases.
For a large contamination proportion (q = 0.09), the C3S-regression, 3S-regression, and 2S-regression perform similarly in the presence of inliers (k ≤ 3), but the 3S-regression breaks down for moderate and large cell-wise outliers (k ≥ 4). However, the C3S-regression is highly robust against large cell-wise outliers (k ≥ 5) although less robust against moderate outliers. In the case-wise contamination setting, the C3S-regression, 3S-regression and 2S-regression perform fairly well and similarly. Nevertheless, the 2S-regression has the best performance, followed by the 3S-regression, which is followed in performance by the C3S-regression.
We also study the performance of the estimator with moderate and large case-wise contamination levels of 10% and 20%, in which at a size of leverage outliers of 22, the C3S-regression and 3S-regression break down as k increases. In this settings, the C3S-regression outperforms the 3S-regression, but, as expected, the 2S-regression maintains its robustness for any contamination level.
Note that, in practice, it is unusual to find case-wise outliers and even more at moderate or large levels. Thus, the loss of robustness for the C3S-regression and 3S-regression does not present a disadvantage. We detect that models with continuous and dummy covariates in the M-regression and C3S-regression outperform in all assumed scenarios. Table 4 reports a summary of the performance of the estimators evaluated by MSE. The performance of the 3S-regression considering non-normal covariates is comparable to all the other estimators for clean data. However, both versions of the C3S-regression outperform all other estimators for any contamination size k in the cell-wise contamination setting. In the cases of non-normal covariates, the C3S-regression maintains its competitive performance, followed by the 3S-regression, while the 2S-regression, as expected, breaks down in the presence of moderate and large cell-wise outliers proportion. Next, the statistical performance of confidence intervals (CIs) for the regression coefficients based on the asymptotic covariance matrix, as described in Subsection 3.3, is evaluated. The asymptotic 100(1 − τ)% CIs for the coefficients of C3S-regression can be established as The performance of CIs defined in Equation (16) may be evaluated using the empirical mean coverage rate (CR) given by and the empirical mean CI length (CIL) defined as Table 5 reports the average CIL defined in Equation (18) obtained from the C3S-regression and 3S-regression in the case of clean data and contaminated data with 1% cell-wise (k = 9), 5% cell-wise (k = 6), 9% cell-wise (k = 3) and 3% case-wise (k = 3), for n ∈ {150, 300, 500, 1000, 5000}. The results of the LS and 2S-regression estimates are not included here, because we are interested in comparing the CIL between the 3S-regression and C3S-regression. The CIL that is obtained from the C3S-regression is comparable to that of the 3S-regression for all considered scenarios. The CIL reached from the 3S-regression are shorter than that for the C3S-regression with clean data and data with small and moderate cell-wise contamination levels. For data with large cell-wise contamination levels or case-wise contamination, the CILs of the C3S-regression are shorter than the CILs of the 3S-regression. Moreover, for any assumed scenario, CILs of the 3S-regression and C3S-regression decrease as the sample size n increases. Figure 7 shows the CR defined in Equation (17) in the case of clean data and contaminated data with 5% cell-wise contamination (k = 5), and 3% case-wise contamination (k = 3), and for different sample sizes n ∈ {150, 300, 500, 1000}. Although the results for the sample size n = 5000 are not shown here for visualization, it can be noticed that, for the C3S-regression and 3S-regression, the evaluations of CR under n = 5000 are better than those when n = 1000. For contamination settings, the 3S-regression yields the best CR, which is the closest to the nominal level. In general, the CR for the C3S-regression is similar to that of the 3S-regression, and it tends to be equal as the sample size n increases.

Analysis of Real Data
The airfoil self-noise data set is used for the illustration purpose. These data were obtained from a series of aerodynamic and acoustic tests of two and three-dimensional airfoil blade sections conducted in an anechoic wind tunnel by the NASA. The data set comprises airfoils of different sizes at various wind tunnel speeds and angles of attack with n = 1503 observations (cases). For this data set, Table 6 shows five covariates and one response variable along with their statistical summaries. This data set is available at the UCI repository [42]. The aim of this empirical study is to predict the noise generated by an airfoil, from dimensions, speed and angle of attack. Specifically, the objective is to explain the scaled sound pressure level. The data set is fitted with the model given by where the log function is used for X 1 due to its wide range and high skewness, while the log function is employed for Y in order to improve the R2-adjusted. The corresponding parameters with C3S-regression (in both versions and full version computed by bootstrap estimation), 2S-regression, 3S-regression, and LS estimates are obtained. The regression coefficient estimates and the corresponding p-values are reported in Table 7. Note that the regression coefficients are similar for all the estimates, except for the covariate X 5 , (that is, the suction side displacement thickness). The coefficient of X 5 estimated by 3S-regression and 2S-regression are similar, but are very different from the C3S-regression and LS estimates. For the C3S-regression, X 5 is highly not significant, while for the 2S-regression and 3S-regression, it is only not significant. However, the LS method indicates that X 5 is significant. The squared norm distance, defined as SND = n ∑ p j=1 ( β j,A − β j,B ) 2 MAD(X ij , . . . , X nj ) 2 , is used to compare the four estimators. Table 8 reports the corresponding SND, which shows that these distances from each two pairs are not large. Therefore, it suggests that the data are not contaminated or the contamination level is very small (inliers).

Conclusions and Future Works
We have provided a new form for robustifying the estimation of parameters of a linear regression model in order to immunize these estimators against case-wise and cell-wise outliers. The main idea here was to modify the generalized Rocke S-estimator in order to obtain robust estimators of the corresponding means and covariances. The difference in our proposal was changing, in the generalized Rocke S-estimator, the initial scatter estimate from the extended minimum volume ellipsoid estimate by the empirical median. The proposed estimator used a univariate filter introduced in the literature and the generalized Rocke S-estimator modified for incomplete data. Our method worked well and similar to that used in the 3S-regression, but in the second step with a different initial robust estimate for the generalized Rocke S-estimator. The initial estimates of location and scatter, the empirical median, and the robust version of the covariance, were computed after snipping the data. Therefore, we have obtained the following findings: • A new method, called comedian-three-step regression, was proposed, which showed an overall outperformance over the recent developed robust methods. • An exact correction factor (b X ) was calculated in order to estimate consistently the standard deviation by using the median absolute deviation for the exponential, logistic and uniform distributions. In addition, a numerical solution for this correction factor was introduced in the Student-t and Weibull distributions.
• In continuous covariates, for small contamination proportion and large cell-wise outliers, the 3S-regression performed similarly to the C3S-regression. However, in general, the C3S-regression outperformed the 3S-regression when the cell-wise contamination proportion increase.

•
In continuous and dummy covariates, the C3S-regression outperformed both the 3S-regression and 2S-regression, for different contamination proportions. However, for the case-wise outliers, the performance of the three estimators was quite similar.

•
The performance of the full version of the C3S-regression estimator proposed in this work was better than its light version. However, the latter one is computationally faster and it can also be used without significant loss of robustness.
Therefore, we have contributed to the robust statistic literature modifying the original three-step regression model by introducing a new family of initial estimates based on the comedian. Our method and the original one are useful to deal with both cell-wise and case-wise outliers. Nevertheless, the numerical results reported that the method proposed in the present study showed an overall outperformance over the recent developed robust methods and a better performance for models with continuous and dummy covariates.
The following aspects derived of this paper may be considered for future work: • The C3S-regression and 3S-regression estimators work well for cell-wise contamination. However, the performance of these estimators with moderate and large case-wise contamination levels (for example, between 10% and 20%) do not work well when the contamination level increases. Some new kind of shrinkage estimator for the initial scatter estimate should be investigated.

•
A bivariate filter can be considered in the first step in order to snip deviation of cells, which could improve the performance of the estimator.

•
A numerical procedure must be studied to calculate the correction factor for any distribution.