پیاده سازی برآوردگر بیانکو و Yohai برای رگرسیون لجستیک

A fast and stable algorithm to compute a highly robust estimator for the logistic regression model is proposed. A criterium for the existence of this estimator at finite samples is derived and the problem of the selection of an appropriate loss function is discussed. It is shown that the loss function can be chosen such that the robust estimator exists if and only if the maximum likelihood estimator exists. The advantages of using a weighted version of this estimator are also considered. Simulations and an example give further support for the good performance of the implemented estimators.

Let Yi, 1⩽i⩽n, be independent Bernoulli variables whose success probabilities depend on the values of p-dimensional explanatory variables X1,…,Xn through the relation where F is a strictly increasing cumulative distribution function. Taking F(u)=1/(1+exp(−u)) results in the logit model, which is the model we will consider in this paper. To simplify notation, we will use γ=(α,βt)t and zi=(1,xit)t for all 1⩽i⩽n. An estimator for γ computed from the sample Xn={(x1,y1),…,(xn,yn)} is denoted . The maximum likelihood (ML) estimator is defined as equation(1.1) where is the conditional log-likelihood function calculated in γ and d(zitγ;yi) is the deviance component given by The ML estimator is the most efficient estimator (asymptotically), but it may behave very poorly in presence of outliers. Therefore, robust alternatives need to be constructed. In this paper, focus is on a generalization of (1.1) which consists of replacing the function d by another one. The estimator of interest is therefore defined by equation(1.2) where ϕ is a positive and almost everywhere differentiable function. It needs to satisfy ϕ(s;0)=ϕ(−s;1) for any score s, where a score value si=zitγ is obtained as a linear combination of a given parameter vector γ. Instead of working with ϕ, we will most of the time use the univariate function φ(s)=ϕ(s;0). The value φ(s) corresponding to an observation with y=0 gives the impact of a particular score s on the value of the objective function in (1.2). This function is assumed to be nondecreasing since a large positive score s should not be attributed to observations having a null y value and should therefore receive a larger weight in the function to minimize. We further require that , implying that a large negative score, what we expect for a value y=0, is not contributing to the objective function in (1.2). As (1.2) shows, the estimator belongs to the class of M-type estimators and the associated first-order condition is given by equation(1.3) where Ψ(s;0)=∂ϕ(s;0)/∂s and Ψ(s;1)=−Ψ(−s;0). Notice that these first-order conditions may have multiple solutions and we need to use the value of the objective function in (1.2) to select the final estimator. Throughout the paper we will use the notation ψ(s)=Ψ(s;0)=φ′(s). Particular cases of (1.2) have already been considered in the literature. First, the ML estimator obviously belongs to this class of M-estimators with φML(s)=−ln(1−F(s)). A more robust proposal is due to Pregibon (1982) who suggested an estimator defined by where λ is a strictly increasing Huber's type function. This estimator was designed to give less weight to observations poorly accounted for by the model but did not downweight influential observations in the design space. Moreover, it was not consistent. Then, Bianco and Yohai (1996) constructed a consistent and more robust version of Pregibon's estimator by working with a bounded function ρ, and defining equation(1.4) with C(γtzi) a bias correction term given by C(s)=G(F(s))+G(1−F(s))−G(1) where Substracting the term G(1) to the bias correction originally proposed by Bianco and Yohai (1996) is necessary to consider this estimator as a particular case of (1.2). Indeed, the estimator of Bianco and Yohai, simply denoted as BY from now on, corresponds to the univariate function φ given by equation(1.5) which satisfies the requirement to go to 0 at −∞. The construction of the BY estimator shows that the resulting φBY function only depends on the choice of the ρ function. In their 1996 paper, BY suggested using the following function: equation(1.6) where c is a tuning parameter, but stressed that other choices are possible. For this particular case, the right panel of Fig. 1 represents the corresponding φBY function. The left panel is for the classical ML. While φBY yields large but bounded values for large positive scores (corresponding to misclassified observations), the φML function returns extremely large values for these. Full-size image (3 K) Fig. 1. φ functions for the ML (left) and BY (right) estimators. Figure options This paper investigates some properties of estimators defined by an optimization problem of type (1.2), where we focus on the BY estimator. The latter estimator has been shown to be consistent, asymptotically normal, and to have good bias and robustness properties, while still being sufficiently efficient (see Bianco and Yohai, 1996). Numerical experiments we carried out confirmed these findings, but we did find one major inconveniency: when working with the loss function (1.6) it occurred frequently that the BY estimator was not existing. Even for samples without any outliers it might happen that the BY-estimator “explodes to infinity,” meaning that the minimum in (1.2) is attained at the edge of the real vector space. The contribution of the paper is twofold: (i) we propose a ρ function guaranteeing existence of the BY estimator as soon as the ML estimator exists; (ii) we propose a stable and quite fast algorithm to compute the BY estimator. The paper is organized as follows: Section 2 derives a criterium for the existence of the estimator at finite samples. This leads to the overlap condition of Albert and Anderson (1984) and Santner and Duffy (1986) for the existence of the ML estimator and to a short discussion in Section 3 of influence functions of M-type estimators in logistic regression. Weighted versions of M-type estimators are introduced in Section 4 while Section 5 outlines a fast and stable algorithm for the computation of the estimator. In Section 6, a simulation experiment and an example compare the BY-estimator with other robust estimators in logistic regression. Finally, Section 7 draws some conclusions.

In this paper, the estimator of BY for logistic regression is revisited. Bianco and Yohai (1996) proved consistency and asymptotical normality of their estimator. Furthermore, the bias behavior of the estimator was shown to be very attractive with respect to other existing estimators. However, due to the lack of efficient algorithms to compute it, this estimator has not been used much in practice. In this paper, a fast and stable algorithm for the BY estimator was proposed. It has been shown that, using the implemented algorithm, this estimator can without any doubt compete with other existing procedures. The good performance of this estimator at finite samples requires a slight modification of the original proposal of Bianco and Yohai (1996). For the newly introduced loss function, the estimator exists for any data set with overlap. This condition is equivalent to the usual requirement for the existence of the classical ML procedure, and is the best one can get. The robustness of the estimator of BY was already observed by Bianco and Yohai (1996). It has been assessed here by means of a study of the influence function of the estimator and its weighted version. Computing other measures of robustness turns out to be less evident in the logistic regression model. For example, computation of the breakdown point, as was done for the ML estimator by Croux et al. (2002), appears to be much more difficult. Furthermore, focus was only put on the estimation part for the logistic regression model, leaving aside the testing problem. Cantoni and Ronchetti (2001) obtained results for the testing part for M-estimators, which could be applicable to the BY estimator as well.