Weighted Maximum Likelihood Technique for Logistic Regression

Idriss Abdelmajid Idriss¹, Weihu Cheng¹, Yemane Hailu Fissuh²

¹Department of Statistics, Faculty of Science, Beijing University of Technology, Beijing, China.
²Department of Statistics, Aksum University, Tigray, Ethiopia.

DOI: 10.4236/ojs.2023.136041                

Abstract

In this paper, a weighted maximum likelihood technique (WMLT) for the logistic regression model is presented. This method depended on a weight function that is continuously adaptable using Mahalanobis distances for predictor variables. Under the model, the asymptotic consistency of the suggested estimator is demonstrated and properties of finite-sample are also investigated via simulation. In simulation studies and real data sets, it is observed that the newly proposed technique demonstrated the greatest performance among all estimators compared.

Keywords

Logistic Regression, Clean Model, Robust Estimation, Contaminated Model, Weighted Maximum Likelihood Technique

1. Introduction

Logistic regression is a method used in statistics for modeling the relationship between a binary dependent variable and one or more explanatory variables. It is usually utilized in various fields. Often, the estimation method for parameters of a statistical model uses the maximum likelihood estimator by optimizing the likelihood function. In logistic regression models, MLE is utilized to estimate the coefficients of the predictor variables that most effectively fit the data. Unfortunately, this approach is not robust to unusual data findings. Several robust estimators have been presented as an alternative to MLE to address this issue. The MLE in logistic regression was shown to be very sensitive to outlying data by [1] , who also devised a diagnostic assessment of outlying observations; see also [2] . For binary regression, [3] studied several M-estimators, their estimations depending on leverage down weight, a Mallows-type estimator. [4] created a robust estimate for the logistic regression model based on a modified median estimator, and they also investigated a Wald-type test statistic. [5] created incredibly robust projection estimators for the GLM, but their calculation is quite difficult. [6] proposed a quasi-likelihood estimator by substituting the least absolute deviation estimator (L1 norm) for the least squares estimator (L2 norm) in the definition of quasi-likelihood. [7] suggested resilient estimators and testing techniques for Poisson and binomial models depending on the idea of a quasi-likelihood estimator created by [8] . The breakdown of the MLE in the logistic model was investigated in [9] , whereas [10] offered a method for the logistic regression model. [11] created a new robust technique for Poisson regression. [12] provided reliable estimators for generalized linear models; the fundamental concept is to change the response via a variance stabilizing transformation before using an estimation. [13] indicated that was a very consistent and reliable estimator. [14] provides a reliable and efficient approach for computing the M-estimator described in [13] . [15] presented a quick technique for the generalized linear model based on breakdown points of the trimmed likelihood. Fisher-consistent estimators are another kind of robust estimator, as introduced in [13] . [16] investigated a resistant robust estimator and this estimation relied on the misclassification model. [17] presented a new family of robust methods for logistic regression. [18] compared minimal distance approaches to more robust methods discovered a unique weighted likelihood and used it in Poisson and binary regression. The ideally bounded scoring functions described by [19] for linear models were applied to the logistic model in [20] .

All of these estimators possess significant disparities in their resistance to outliers and effectiveness within the framework of the model. In this paper, we have conducted a comprehensive investigation into the behavior of some of these estimators, both in terms of their asymptotic properties and their behavior in finite samples. Our findings indicate that the Mallows-type estimator proposed by [3] is very robust to outlier contamination but inefficient under the model, while Schweppe-type estimators proposed by [2] are very efficient under the model but show a poor outlier resistance. In this paper, we propose an estimator that can be as robust as Mallows-estimators under contamination but is much more efficient under the model, this is achieved by an adaptive continuous weight. This continuous weighted maximum likelihood estimator depends on the annoyance parameter estimator as a result of Kolmogorov-Smirnov statistics. The maximum likelihood estimators and a logistic regression model are covered in Section 2. Robust Methods for Logistic Regression Model discussed in section 3. We propose a robust technique for logistic regression in section 4. Section 5 displays the findings of the Monte Carlo simulation research and real data analysis. Section 6 contains the conclusions.

2. Logistic Regression Model and ML Estimator

The logistic regression model is a popular technique used to examine the relationship between a categorical variable and one or more predictor variables. It is often used in binary classification issues when there are only two potential outcomes for the dependent variable (e.g., true or false, yes or no, etc.). The logistic function is the foundation of the logistic regression model, which transforms a continuous input variable into a probability value between 0 and 1. Assume a random sample of observations $\left(x_1\mathrm{,}{y}_1\right)\mathrm{,}\cdots \mathrm{,}\left(x_n\mathrm{,}{y}_n\right)$ , where $X_i$ represents a p predictor variables and $y_i\in \left(\mathrm{0,1}\right)$ is a binary variable and suppose the probability of positive response ${\mu }_i=p\left(y_i=1|x_i\right)$ is associated with the covariates through the relation $g\left(\mu \right)=X\beta$ , such that $g^{-1}\left(X\beta \right)$ is known as the logit link function converts the covariate values in the range $\left(\mathrm{0,1}\right)$ .

Using the logit link function The multiple logistic regression model may be created by:

$p\left(Y=1|X=x_i\right)=F\left(x_i^{\text{T}}\beta \right)=\frac{\exp \left(x_i^{\text{T}}\beta \right)}{1+\exp \left(x_i^{\text{T}}\beta \right)},i=1,\cdots ,n,$ (1)

where $X=\left(\mathrm{1,}{x}_1\mathrm{,}\cdots \mathrm{,}{x}_p\right)$ are the predictor variables values and ${\beta }^{\text{T}}=\left({\beta }_0\mathrm{,}{\beta }_1,\cdots \mathrm{,}{\beta }_p\right)$ represents an unknown parameter vector. We may characterize the binary regression model as follows:

${\eta }_i=x_i^{\text{T}}\beta \mathrm{,}$

where ${\eta }_i$ is a linear predictor that is sometimes referred to as a transformation function, where ${\eta }_i=\ln \left(\frac{{\mu }_i}{1-{\mu }_i}\right)$ . The MLE is a method for parameter estimation

of a statistical model by maximization of the likelihood function. The MLE assumes that the data are generated by a specific probability distribution (in this case, the logistic distribution) and finds the parameter values that make it most likely for that distribution to have generated those data. The MLE is often used in logistic regression because it provides unbiased estimates of model parameters and has good statistical properties. Assume that the dependent variables $y_i$ have a Bernoulli distribution and we can derive the probability distribution for the i^th observation as follows:

$f\left(y_i\right)={\mu }_i^{y_i}{\left(1-{\mu }_i\right)}^{1-y_i},i=1,2,\cdots ,n,$

as well as each observation $y_i$ takes the value 1 with probability ${\mu }_i$ or the value 0 with probability ( $1-{\mu }_i$ ). The probability function is defined as follows

$\ln \left(\beta ;y_1,y_2,\cdots ,y_n\right)={\displaystyle \underset{i=1}{\overset{n}{\prod }}}f\left(y_i\right)={\displaystyle \underset{i=1}{\overset{n}{\prod }}}{\mu }_i^{y_i}{\left(1-{\mu }_i\right)}^{1-y_i},$ (2)

Then we compute the log-likelihood of the preceding formula:

$\ln \left(\beta ;y_1,y_2,\cdots ,y_n\right)=\ln {\displaystyle \underset{i=1}{\overset{n}{\prod }}}f\left(y_i\right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}}\left[y_i\ln \left(\frac{{\mu }_i}{1-{\mu }_i}\right)\right]+{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\ln \left(1-{\mu }_i\right),$

where ${\eta }_i=\ln \left[{\mu }_i/\left(1-{\mu }_i\right)\right]=x_i^{\text{T}}\beta$ and $1-{\mu }_i={\left[1+\exp \left(x_i^{\text{T}}\beta \right)\right]}^{-1}$ . So, we can express the log-likelihood as:

$\begin{array}{c}\ln l\left(\beta ;Y\right)={\displaystyle \underset{i=1}{\overset{n}{\sum }}}\ln y_i{x}_i^{\text{T}}\beta -{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\ln \left[1+\exp \left(x_i^{\text{T}}\beta \right)\right]\\ ={\beta }^{\text{T}}XY-{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\ln \left[1+\exp \left(x_i^{\text{T}}\beta \right)\right].\end{array}$ (3)

In experimental design, we execute repeated observations at each level of independent variables (x). Let ${\eta }_i$ be the trials number at each predictor level and $y_i$ be the number of 1’s observed at the i^th observations with $n=n_1,n_2,\cdots ,n_m$ . Therefore, we may express the log-likelihood as:

$\ln l\left(\beta ;Y\right)={\beta }^{\text{T}}XY-{\displaystyle \underset{i=1}{\overset{m}{\sum }}}{n}_i\ln \left[1+\exp \left(x_i^{\text{T}}\beta \right)\right],$ (4)

nevertheless, it is possible to maximize the likelihood function by differentiating it with respect to $\beta$ :

$\frac{\partial \ln l\left(\beta ,Y\right)}{\partial \beta }=X^{\text{T}}Y-{\displaystyle \underset{i=1}{\overset{m}{\sum }}}\left[\frac{n_i}{1+{\text{e}}^{x_i^{\text{T}}\beta }}\right]{\text{e}}^{x_i^{\text{T}}\beta }{X}_i,$

where, ${\text{e}}^{x_i^{\text{T}}}/\left(1+{\text{e}}^{x_i^{\text{T}}\beta }\right)=1/\left(1+{\text{e}}^{-x_i^{\text{T}}\beta }\right)={\mu }_i$ , we have

$\frac{\partial \ln \left(\beta ;Y\right)}{\partial \beta }=X^{\text{T}}Y-{\displaystyle \underset{i}{\overset{m}{\sum }}}{n}_i{\mu }_i{X}_i,$

where $n_i{\mu }_i$ denotes the average of the binomial variable, the preceding formula may be written in matrix form as $X^{\text{T}}\left(Y-\mu \right)$ , where

$\mu =\left[\begin{array}{c}{\mu }_1\\ {\mu }_2\\ ⋮\\ {\mu }_m\end{array}\right]\mathrm{.}$

Hence, the MLE is normally computed by resolving the scoring equation:

$X^{\text{T}}\left(Y-\mu \right)=0.$ (5)

Since Equation (5) is a nonlinear function of beta, the iteratively weighted least squares (IWLS) technique may be used to find a solution.

3. Robust Methods for Logistic Regression Model

In logistic regression models, robust estimators are statistical methods for estimating parameters that are less sensitive to outliers and influential observations. These techniques are designed to provide reliable estimates of the regression coefficients even when the data contain extreme values or other anomalies that can distort the results. There are several robust estimators that can be used in logistic regression. [13] suggested robust estimators via the management of deviations to get theoretically unbiased estimators; however, an extra bias-correction component must be introduced, which makes the calculation of their estimator very complex and the estimator itself not straightforward. [2] developed Mallows-type estimators by separately manipulating the variables and residuals in the estimation equation. In the instance of logistic regression, they were categorized as solutions of

${\displaystyle \underset{i=1}{\overset{n}{\sum }}}w\left(x_i;\hat{\eta }\right){\pi }_b\left(y_i-\pi \left(x_i^{\text{T}}{\hat{\beta }}_i\right)-c\left(\pi \left(x_i^{\text{T}}\hat{\beta }\right),b\right)\right)x_i=0,$ (6)

where $\hat{\eta }$ represents the annoyance parameters (location and scatter estimate of covariates), and ${\Phi }_b$ is often considered as Huber’s function ${\pi }_b\left(t\right)=\left(-b\right)\vee \left(t\wedge b\right)$ , $c\left(t\mathrm{,}b\right)$ is the bias-correction function expressed as

$c(t,b)=(\begin{array}{ll}b\pi /1-\pi \left(t\right)-\pi \left(t\right)\hfill & \text{if}t<0,b<1-\pi \left(t\right)\hfill \\ 1-\pi \left(t\right)-b1-\pi \left(t\right)/\pi \left(t\right)\hfill & \text{if}t>0,b<\pi \left(t\right)\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}$

the weights $w\left(x_i\mathrm{<mo>\; ;\; </mo>}\hat{\eta }\right)$ often depend only on continuous variables. Suppose we write $x_i^{\text{T}}=\left(u_i^t\mathrm{,}{z}_i^{\text{T}}\right)$ , where $u_i\in R^{p-q}$ are qualitative variables and $z_i\in R^q$ are the continuous variables, thereafter, the weights usually take the form $w\left(x_i\mathrm{<mo>\; ;\; </mo>}\hat{\eta }\right)=w\left({\left(z_i-\hat{\mu }\right)}^{\text{T}}{\hat{\Sigma }}^{-1}\left(z_i-\hat{\mu }\right)/t\right)$ , with $w\mathrm{<mo>\; :\; </mo>}{R}_{+1}\to R_{+1}$ a function that does not increase, $\hat{\mu }$ is the robust estimator of the location of the $z_i$ , $\hat{\Sigma }$ is the robust estimator of scatter of the $z_i$ and t is the threshold value (usually $t={\chi }_{q\mathrm{,1}-\alpha }$ for some $\alpha =\left(\mathrm{0,1}\right)$ ). The initial robust estimator of location and scale for predictor variables $\hat{\mu }$ and $\hat{\Sigma }$ can be calculated utilizing the minimum covariance determinant methods. MCD was one of the first multivariate location and scatter estimators that were both affine equivariant and relatively robust. It discovers the $h\left(>n/2\right)$ observations $x_i$ whose matrix of traditional covariance

$v=\frac{1}{h}{\displaystyle \underset{i}{\sum }}\left(x_i-t\right){\left(x_i-t\right)}^{\text{T}}$

has the smallest possible determinant, where $t=x$ the average of those h points. As observed, the residual weight ${\varphi }_b$ and covariate weight w are independent; moreover, it will decrease the effectiveness of the resultant estimators due to the fact that the estimation equation will downweight well-fitted observations with extreme variables. [21] robustly presented a family of resilient adaptive weighted maximum likelihood techniques for logistic regression models. These adaptive weights are dependent on adaptive cut-off thresholds to regulate observations with extreme covariables. They demonstrated that estimators based on adaptive thresholds are more efficient than estimators based on non-adaptive thresholds for clean models and have equivalent resilience for polluted models.

The lack of dependence between the weights assigned to covariates and the weights assigned to deviances in Equation (6) is the underlying cause for the generally lower efficiency of Mallows-type estimators compared to Schweppetype estimators. This occurs due to the downweighting of observations with extreme covariates, even if they exhibit good fit. It is evident that enhancing the efficiency of Mallows-type estimators can be achieved by reducing the thresholding proportions, although this may compromise the robustness of the estimator. In order to simultaneously achieve both high efficiency and high robustness, it becomes necessary to employ adaptive thresholds, as detailed in the next section.

4. Weighted Maximum Likelihood Technique (WMLT)

In this section, we build a novel class of continuous weighted maximum likelihood estimators based on the annoyance parameter estimator as a function of the Kolmogorov-Smirnov statistics. We’ll refer to these estimators as WMLT (weighted maximum likelihood technique). First build two estimators ${\hat{\mu }}^{\left(0\right)}$ and ${\hat{\Sigma }}^{\left(0\right)}$ that are the initial estimates of location and scatter of the predictor variables ${z^{\prime }}_{i}s$ , thereafter, calculate the Mahalanobis distances squared of ${z^{\prime }}_{i}s$ that

is characterized by $m^2={\left(z_i-{\hat{\mu }}^{\left(0\right)}\right)}^{\text{T}}{\left({\hat{\Sigma }}^{\left(0\right)}\right)}^{-1}\left(z_i-{\hat{\mu }}^{\left(0\right)}\right)$ . Furthermore, the empirical distribution function of $m_i^2$ may be expressed as:

$F_n\left(t\right)=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}I\left(\left|m_i^2\right|\le t\right),$

when ${z^{\prime }}_{i}s$ have a normal distribution, $F_n$ converge to $F_{{\chi }_q^2}$ ( ${\chi }_q^2$ distribution function). Then we can estimate the outliers proportion in covariates by [21]

${\alpha }_n={\sup }_{t\ge F_{{\chi }_q^2}^{-1}\left(1-\delta \right)}{\left\{F_{{\chi }_q^2}\left(t\right)-F_n\left(t\right)\right\}}_{+}=\underset{i\ge i_0}{\max }{\left\{F_{{\chi }_q^2}\left(m_{\left(i\right)}^2\right)-\frac{i-1}{n}\right\}}_{+},$

where ${\left\{\mathrm{.}\right\}}_{+}$ represents the positive part, $\delta$ measures the length of the tail ( $\delta =0.25$ ) is an acceptable option, and $i_o=\min \left\{i:m_i^2\ge F_{{\chi }_q^2}^{-1}\left(1-\delta \right)\right\}$ . When $\left|F_{{\chi }_q^2}\left(t\right)-F_n\left(t\right)\right|$ is large for a large t, it indicates that the sample contains outliers. Hence, an adaptive threshold may be described as

$t_n=F_n^{-1}\left(1-{\alpha }_n\right)=m_{\left(n-\left[n{\alpha }_n\right]\right)}^2.$

[21] proposed the adaptive threshold estimators. Specifically, they developed an estimator of the Mallows-type estimator with weights $w\left(x_i\mathrm{<mo>\; ;\; </mo>}\hat{\eta }\right)=w\left(m_i^2/t_n\right)$ which are basically weighted maximum likelihood estimators. The proposed weight function may be defined as

$w\left(x_i\mathrm{<mo>\; ;\; </mo>}{\alpha }_n\right)=m\left({\alpha }_n{m}_i^2\right)\mathrm{,}$

we assume m is a completely non-increasing continuous mapping from ${ℝ}^{+}$ to $\left(\mathrm{0,1}\right)$ such that $m\left(0\right)=1$ , ${\sup }_{x>0}\left\{xm\left(x\right)\right\}<\infty$ , and the first derivative is bounded with $m^{\left(1\right)}\left(0\right)=0$ . Determine an objective function

${\psi }_{\beta ,\eta }\left(x,y\right)=w_{\eta }\left(x\right){\varphi }_{\beta }\left(x,y\right),x,\beta \in {ℝ}^p,y\in \left(0,1\right),$ (7)

where $\eta =\left(\alpha \mathrm{,}\mu \mathrm{,}\Sigma \right)$ is the location, scale, and goodness of fit of the explanatory variables respectively, an adaptive weight function $w_{\eta }\left(x\right)=m\left(\alpha {\left(z-\mu \right)}^{\text{T}}{\Sigma }^{-1}\left(z-\mu \right)\right)$ and ${\varphi }_{\beta }\left(x\mathrm{,}y\right)=\left(y-\left(x^{\text{T}}\beta \right)\right)x$ , with $\mu \in {ℝ}^a$ , $\Sigma$ is $q\times q$ real matrix and $x^{\text{T}}=\left(u^{\text{T}}\mathrm{,}{z}^{\text{T}}\right)$ . Finally, we define our adaptive estimator ${\hat{\beta }}_n$ of $\beta$ as the solution to the estimating equation

${\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\psi }_{\beta },\hat{\eta }\left(x_i,y_i\right)=0,$

where $\hat{\eta }$ is a consistent estimator of $\eta =\left(\alpha \mathrm{,}\mu \mathrm{,}\Sigma \right)$ .

If $z_i$ are normally distributed, then $t_n=\infty$ and the reweighted estimators are asymptotically equivalent to the sample mean and the sample covariance matrix, and therefore fully efficient. This efficient carries over to the adaptive Mallows-type estimators, as shown in the next subsection.

4.1. Asymptotic of Proposed Method

The estimating equation (6) can be written as ${\displaystyle {\sum }_{i=1}^n}\Psi \left(x_i,y_i;\hat{\beta };\hat{\eta }\right)=0$ , with $\Psi \left(x,y;\beta ;\eta \right)=w\left(x;\eta \right){\psi }_b\left(y-F\left(x^{\text{T}}\beta \right)-c\left(x^{\text{T}}\beta ,b\right)\right)x$ . Under appropriate regularity conditions, the classical asymptotic of M-estimators hold, see [22] . Let ${\beta }_0$ be the model parameter and $E_0$ denote expectation under the model; define $M_0\left(\beta \right)=E_0\left\{\Psi \left(x\mathrm{,}y\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{<mo>\; ;\; </mo>}{\eta }_0\right)\right\}$ , with ${\eta }_0=p\lim \hat{\eta }$ . If $C_0={D{M}_0\left(\beta \right)|}_{\beta ={\beta }_0}$ , where D denote the differential and $A_0=E_0\left\{\Psi \left(x\mathrm{,}y\mathrm{<mo>\; ;\; </mo>}{\beta }_0\mathrm{<mo>\; ;\; </mo>}{\eta }_0\right)\Psi {\left(x\mathrm{,}y\mathrm{<mo>\; ;\; </mo>}{\beta }_0\mathrm{<mo>\; ;\; </mo>}{\eta }_0\right)}^{\text{T}}\right\}$ , then

$\sqrt{n}\left(\hat{\beta }-{\beta }_0\right)\to N_p\left(\mathrm{0,}{C}_0^{-1}{A}_0{C}_0^{-1}\right)\mathrm{.}$

This result is valid for non-adaptive and adaptive weights alike, as long as $\hat{\eta }$ converges to ${\eta }_0$ in probability. For the weighted maximum likelihood estimator given by ${\psi }_b\left(u\right)=u$ , $A_0$ and $C_0$ have simple expressions

$C_0=-E\left\{w\left(x\mathrm{<mo>\; ;\; </mo>}{\eta }_0\right)F^{\prime }\left(x^{\text{T}}{\beta }_0\right)x{x}^{\text{T}}\right\}\mathrm{,}$ (8)

$A_0=E\left\{w^2\left(x\mathrm{<mo>\; ;\; </mo>}{\eta }_0\right)F\left(x^{\text{T}}{\beta }_0\right)\left(1-F\left(x^{\text{T}}{\beta }_0\right)\right)x{x}^{\text{T}}\right\}\mathrm{.}$ (9)

Using the asymptotic normality of $\hat{\beta }$ it is possible to construct confidence ellipsoids for ${\beta }_0$ . First, we estimate the matrices (9) and (6) with

$\hat{C}=-\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}w\left(x;\hat{\eta }\right)F^{\prime }\left(x_i^{\text{T}}\hat{\beta }\right)x_i{x}_i^{\text{T}},$

$\hat{A}=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}w\left(x_i;\hat{\beta }\right)F\left(x_i^{\text{T}}\hat{\beta }\right)\left(1-F\left(x_i^{\text{T}}\hat{\beta }\right)\right)x_i{x}_i^{\text{T}}.$

Then the estimated asymptotic variance of $\sqrt{n}\left(\hat{\beta }-{\beta }_0\right)$ is $\hat{V}={\hat{C}}^{-1}\hat{A}{\hat{C}}^{-1}$ . The asymptotic confidence ellipsoid of level $\left(1-\alpha \right)$ for ${\beta }_0$ is given $\xi \left({\beta }_0\right)=\left\{\beta \in R^p\mathrm{<mo>\; :\; </mo>}n{\left(\hat{\beta }-\beta \right)}^{\text{T}}{\hat{V}}^{-1}\left(\hat{\beta }-\beta \right)\le {\chi }_{p\mathrm{,1}-\alpha }^2\right\}$ . This can be generalized to linear transformations of ${\beta }_0$ .

4.2. Asymptotic Properties of WMLT

This subsection focuses on the asymptotic features of the suggested estimator ${\hat{\beta }}_n$ described in the preceding section. We will show that the estimator is asymptotically consistent based on some general assumptions about the moments of predictor variables. Suppose that ${\beta }_0$ , ${\mu }_0$ and ${\Sigma }_0$ are the actual values of $\beta$ , $\mu$ and $\Sigma$ respectively, and the independent sample $\left(x_1\mathrm{,}{y}_1\right)\mathrm{,}\cdots \mathrm{,}\left(x_n\mathrm{,}{y}_n\right)$ follows the logistic model $P\left(y_i=1\right)={\mu }_i$ , $i=1,\cdots ,n$ . We define functions

${\mathbb{P}}_n{\psi }_{\beta }={\Psi }_n\left(\beta \right)=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\psi }_{\beta },\hat{\eta }\left(x_i,y_i\right)$

and

$\Psi \left(\beta \right)=P_{\psi \beta ,0}=P_{\varphi \beta },$

where ${\psi }_{\beta \mathrm{,0}}$ is the calculated by Equation (7) with $\eta$ replaced with ${\eta }_0=\left(\mathrm{0,}{\mu }_0\mathrm{,}{\Sigma }_0\right)$ and p indicates the joint probability distribution of the (x,y)’s. Theorem 4.2 establishes that ${\hat{\beta }}_n$ is consistent, which makes use of the conclusions of Lemma 1, 2, and 3 mentioned below. To prove the Lemmas and the Theorems, the following assumptions must be made:

B1: $\hat{\mu }{\to }^p{\mu }_0$ and $\hat{\Sigma }{\to }^p{\Sigma }_0$ .

B2: $E{\left(xx\right)}^{\text{T}}$ is nonsingular.

B3: $E_{G_0}\left({\Vert x\Vert }^4\right)<\infty$ .

B4: The $m\left(x\right)$ is a continuous weight function and has a first derivative that is bounded, $m\left(0\right)=1$ and $m^1\left(0\right)=0$ .

B1 is met for the vast majority of well-known initial robust estimators, such as the Minimum Covariance Determinant used in the simulation experiments. In the following lemmas and theorem, the asymptotic characteristics are assumed to be as $n\to \infty$ . The proof for Lemma 1 is given in [21] .

Lemma 1. In case B1 is holds, then ${\alpha }_n=0_p\left(1\right)$ .

Lemma 2. In case B2 is holds, then $\Vert \Psi \left({\beta }_n\right)\Vert \to 0$ implies $\Vert {\beta }_n-{\beta }_0\Vert \to 0$ for any sequence $\left\{{\beta }_n\right\}\in \Theta$ .

Proof of Lemma 2. Consider the following

$\begin{array}{c}\Vert \Psi \left({\beta }_n\right)\Vert =\Vert P{\psi }_{{\beta }_n},0\Vert =\Vert P{\varphi }_{{\theta }_n}\Vert =\Vert P\left(y-\left(x^{\text{T}}{\beta }_n\right)\right)x\Vert \\ =\Vert P\left(y-\left(x^{\text{T}}{\beta }_0\right)+\left(x^{\text{T}}{\beta }_0\right)-\left(x^{\text{T}}{\beta }_n\right)\right)x\Vert \\ =\Vert P\left(\left(x^{\text{T}}{\beta }_0\right)\right)x+P\left(\left(x^{\text{T}}{\beta }_0\right)-\left(x^{\text{T}}{\beta }_n\right)\right)x\Vert \end{array}$

If $\Vert \Psi \left({\beta }_n\right)\Vert \to 0$ , then $\Psi \left({\beta }_n\right)\to 0$ . Since $P\left(y-\left(x^{\text{T}}{\beta }_n\right)\right)x=0$ , we may see from the equality above that $P\left(\left(x^{\text{T}}{\beta }_0\right)-\left(x^{\text{T}}{\beta }_n\right)\right)x\to 0$ . Note

$\begin{array}{c}P\left(\left(x^{\text{T}}{\beta }_0\right)-\left(x^{\text{T}}{\beta }_n\right)\right)x=P\left({}^{\left(1\right)}\left(c\right)\left(x^{\text{T}}{\beta }_n-x^{\text{T}}{\beta }_0\right)x\right)\\ \le \frac{1}{4}P\left(x^{\text{T}}\left({\beta }_n-{\beta }_0\right)x\right),\end{array}$

where $c\in \left(x^{\text{T}}{\beta }_0\mathrm{,}{x}^{\text{T}}{\beta }_n\right)$ and $h^{\prime }\left(\mu \right)$ is the initial derivation for logistic link, $h^{\prime }\left(\mu \right)=\frac{{\text{e}}^x}{{\left({\text{e}}^x=1\right)}^2}\in \left(0,1/4\right)$ . We also have $X={\left(x_1\mathrm{,}{x}_2\mathrm{,}\cdots \mathrm{,}{x}_p\right)}^{\text{T}}$ , ${\beta }_0-{\beta }_n={\beta }_d={\left({\beta }_{d1},{\beta }_{d2},\cdots ,{\beta }_{dp}\right)}^{\text{T}}$ .

$P\left(x^{\text{T}}\left({\beta }_0-{\beta }_n\right)x\right)=P\left({\displaystyle \underset{i=1}{\overset{p}{\sum }}}{\beta }_{di}\left[\begin{array}{c}{x}_1{x}_i\\ x_2{x}_i\\ ⋮\\ x_p{x}_i\end{array}\right]\right)={\displaystyle \underset{i=1}{\overset{p}{\sum }}}{\beta }_{di}P\left(\left[\begin{array}{c}{x}_1{x}_i\\ x_2{x}_i\\ ⋮\\ x_p{x}_i\end{array}\right]\right)\to 0.$

Since $P\left(x{x}^{\text{T}}\right)=E\left(x{x}^{\text{T}}\right)$ is nonsingular, we have ${\beta }_0-{\beta }_1\to 0$ , is shown by the fact that $\Vert {\beta }_n-{\beta }_0\Vert \to 0$ .

Lemma 3. Suppose that B3 and B4 are true. Then the class $\left\{{\psi }_{\eta ,\beta }:\beta \in \Theta ,\Vert \eta -{\eta }_0\Vert <\delta \right\}$ is P-Glivenko-Cantelli for some $\delta >0$ , where ${\eta }_0=\left(\mathrm{0,}{\mu }_0\mathrm{,}{\Sigma }_0\right)$ .

Proof of Lemma 3. To demonstrate that a class $Ϝ$ of vector-valued functions $\psi \mathrm{<mo>\; :\; </mo>}\left(x\mathrm{,}y\right)\to {ℝ}^p$ to be Glivenko-Cantelli, we must display each of the coordinate classes ${\phi }^i\mathrm{<mo>\; :\; </mo>}\left(x\mathrm{,}y\right)\to ℝ$ with $\psi ={\left({\psi }^1\mathrm{,}\cdots \mathrm{,}{\psi }^p\right)}^{\text{T}}$ ranging over $Ϝ\left(i=\mathrm{1,2,}\cdots \mathrm{,}p\right)$ is Glivenko-Cantelli.

The class $Ϝ=\left\{{\psi }_{\gamma }\mathrm{<mo>\; :\; </mo>}\gamma =\left(\beta \mathrm{,}\eta \right)=\left(\beta \mathrm{,}\alpha \mathrm{,}\mu \mathrm{,}\Sigma \right)\mathrm{,}\alpha \in \mathrm{<mo\; stretchy="false">\; [\; </mo>0,1\; <mo\; stretchy="false">\; ]\; </mo>,}\mu \in {ℝ}^q\mathrm{,}\Sigma \in {\mathbb{S}}_{+}^q\mathrm{,}\beta \in \Theta \mathrm{,}\Vert \eta -{\eta }_0\Vert <\delta \right\}$ is a set of measurable functions that are indexed by a bounded subset in $\Gamma \subset \Theta \times {ℝ}^q\times ℝ\times {ℝ}^{q\times q}$ , and $S_{+}^q$ represents a collection of positive semidefinite matrices described in ${ℝ}^{q\times q}$ . This is because $\Sigma$ is basically a variance-covariance matrix of continuous predictor variables, hence it is positive semidefinite and symmetric. For the norm, we use

$\Vert \eta -{\eta }_n\Vert ={\left({\alpha }^2+{\Vert \mu -{\mu }_0\Vert }^2+\Vert \Sigma -{\Sigma }_0\Vert +{\Vert \beta -{\beta }_0\Vert }^2\right)}^{\frac{1}{2}}$ , with $\Vert .\Vert$ indicates the Euclidean norm of vectors $\mu \mathrm{,}\beta$ and for matrix $\Sigma$ , $\Vert .\Vert$ represents the induced norm in general. For two values ${\gamma }_i=\left({\beta }_i\mathrm{,}{\eta }_i\right)$ , $\left(i=1,2\right)$ of $\gamma$ , we have

$\begin{array}{l}\left|{\psi }_{{\gamma }_1}^i\left(x\mathrm{,}y\right)-{\psi }_{{\gamma }_2}^i\left(x\mathrm{,}y\right)\right|\\ =\left|w_{{\eta }_1}\left(x\right){\varphi }_{{\beta }_1}^i\left(x\mathrm{,}y\right)-w_{{\eta }_2}\left(x\right){\varphi }_{{\beta }_2}^i\left(x\mathrm{,}y\right)\right|\\ =\left|w_{{\eta }_1}\left(x\right){\varphi }_{{\beta }_1}^i\left(x\mathrm{,}y\right)-w_{{\eta }_1}\left(x\right){\varphi }_{{\beta }_1}^i\left(x\mathrm{,}y\right)+w_{{\eta }_1}\left(x\mathrm{,}y\right){\varphi }_{{\beta }_2}^i\left(x\mathrm{,}y\right)-w_{{\eta }_2}\left(x\right){\varphi }_{{\beta }_2}^i\left(x\mathrm{,}y\right)\right|\\ \le \left|w_{{\eta }_1}\left(x\right){\varphi }_{{\beta }_1}^i\left(x\mathrm{,}y\right)-w_{{\eta }_1}\left(x\right){\varphi }_{{\beta }_2}^i\left(x\mathrm{,}y\right)\right|+\left|w_{{\eta }_1}\left(x\right){\varphi }_{{\beta }_2}^i\left(x\mathrm{,}y\right)-w_{{\eta }_2}\left(x\right){\varphi }_{{\beta }_2}^i\left(x\mathrm{,}y\right)\right|\mathrm{.}\end{array}$ (10)

But

$\begin{array}{l}\left|w_{{\eta }_1}\left(x\right){\varphi }_{{\beta }_1}^i\left(x,y\right)-w_{{\eta }_1}\left(x\right){\varphi }_{{\beta }_2}^i\left(x,y\right)\right|\\ =w_{{\eta }_1}\left(x\right)\left|\left(y-\left(x^{\text{T}}{\beta }_1\right)\right)x_i-\left(y-\left(x^{\text{T}}{\beta }_2\right)\right)x_i\right|\le \left|\left(x^{\text{T}}{\beta }_1\right)-\left(x^{\text{T}}{\beta }_2\right)\right|\left|x_i\right|\\ \le K_0\left|x_i\right|\left|x^{\text{T}}{\beta }_2-x^{\text{T}}{\beta }_1\right|\le K_0\left|x_i\right|\Vert c\Vert \Vert {\beta }_2-{\beta }_1\Vert ,\end{array}$ (11)

where $K_0$ represents the upper limit of the first derivation of the link function $h^{\prime }\left(\mu \right)$ .

Making use of the Mean Value Theorem, for each $x_1$ and $x_2$ , there exists $c\in \left(x_1,x_2\right)$ such that

$\Vert \mu \left(x_1\right)-\mu \left(x_2\right)\Vert =h^{\prime }\left(c\right)\left|x_1-x_2\right|<K_0\left|x_1-x_2\right|.$

In addition, we have

$\begin{array}{l}\left|w_{{\eta }_1}\left(x\right){\varphi }_{{\beta }_2}^i\left(x,y\right)-w_{{\eta }_2}\left(x\right){\varphi }_{{\beta }_2}\left(x,y\right)\right|=\left|w_{{\eta }_1}\left(x\right)-w_{{\eta }_2}\left(x\right)\right|\left|y-\left(x^{\text{T}}{\beta }_2\right)\right|\left|x_i\right|\\ =\left|w_{{\alpha }_1,{\mu }_1,{\Sigma }_1}\left(x\right)-w_{{\alpha }_1,{\mu }_2,{\Sigma }_2}\left(x\right)+w_{{\alpha }_1,{\mu }_2,{\Sigma }_2}\left(x\right)-w_{{\alpha }_2,{\mu }_2,{\Sigma }_2}\left(x\right)\right|\left|y-\left(x^{\text{T}}{\beta }_2\right)\right|\left|x_i\right|\\ \le \left(\left|w_{{\alpha }_1,{\mu }_1,{\Sigma }_1}\left(x\right)-w_{{\alpha }_1,{\mu }_2,{\Sigma }_2}\left(x\right)\right|+\left|w_{{\alpha }_1,{\mu }_2,{\Sigma }_2}\left(x\right)-w_{{\alpha }_2,{\mu }_2,{\Sigma }_2}\left(x\right)\right|\right)\left|x_i\right|,\end{array}$ (12)

then, we get

$\begin{array}{l}\left|w_{{\alpha }_1,{\mu }_2,{\Sigma }_2}\left(x\right)-w_{{\alpha }_2,{\mu }_2,{\Sigma }_2}\left(x\right)\right|\\ =\left|m\left({\alpha }_1{\left(Z-{\mu }_2\right)}^{\text{T}}{\Sigma }_2^{-1}\left(Z-{\mu }_2\right)\right)-m\left({\alpha }_2{\left(Z-{\mu }_2\right)}^{\text{T}}{\Sigma }_2^{-1}\left(Z-{\mu }_2\right)\right)\right|\\ \le K_1\left|{\alpha }_2-{\alpha }_1\right|{\left(Z-{\mu }_2\right)}^{\text{T}}{\Sigma }_2^{-1}\left(Z-{\mu }_2\right),\end{array}$ (13)

and

$\begin{array}{l}{w}_{{\alpha }_1,{\mu }_1,{\Sigma }_1}\left(x\right)-w_{{\alpha }_1,{\mu }_2,{\Sigma }_2}\left(x\right)\\ =\left|m\left({\alpha }_1{\left(Z-{\mu }_1\right)}^{\text{T}}{\Sigma }_1^{-1}\left(Z-{\mu }_1\right)\right)-m\left({\alpha }_1{\left(Z-{\mu }_2\right)}^{\text{T}}{\Sigma }_2^{-1}\left(Z-{\mu }_2\right)\right)\right|\\ \le K_1{\alpha }_1\left|{\left(Z-{\mu }_1\right)}^{\text{T}}{\Sigma }_1^{-1}\left(Z-{\mu }_1\right)-{\left(Z-{\mu }_2\right)}^{\text{T}}{\Sigma }_2^{-1}\left(Z-{\mu }_2\right)\right|\\ \le K_1\left|{\left(Z-{\mu }_1\right)}^{\text{T}}{\Sigma }_1^{-1}\left(Z-{\mu }_1\right)-{\left(Z-{\mu }_1\right)}^{\text{T}}{\Sigma }_2^{-1}\left(Z-{\mu }_1\right)\right|\\ +K_1\left|{\left(Z-{\mu }_1\right)}^{\text{T}}{\Sigma }_2^{-1}\left(Z-{\mu }_1\right)-{\left(Z-{\mu }_2\right)}^{\text{T}}{\Sigma }_2^{-1}\left(Z-{\mu }_2\right)\right|.\end{array}$ (14)

Given that $\Sigma$ is a positive semidefinite matrix, as well ${\Sigma }^{-1}$ . Indicate what the eigenvalues of ${\Sigma }_2^{-1}$ are as ${\lambda }_1\ge {\lambda }_2\ge \cdots \ge {\lambda }_p\ge 0$ . Consequently, ${\Sigma }_2^{-1}$ has a set of orthogonal eigenvectors, say $p_1\mathrm{,}\cdots \mathrm{,}{p}_q$ , s.t. ${\Sigma }_2^{-1}{p}_i={\lambda }_i{p}_i$ . There exists an orthogonal matrix in matrix format Q s.t.

$Q^{-1}{\Sigma }_2^{-1}Q=Q^{\text{T}}{\Sigma }_2^{-1}Q=\Lambda .$

Then we get

$\begin{array}{l}\left|{\left(Z-{\mu }_1\right)}^{\text{T}}{\Sigma }_2^{-1}\left(Z-{\mu }_1\right)-{\left(Z-{\mu }_2\right)}^{\text{T}}{\Sigma }_2^{-1}\left(-{\mu }_2\right)\right|\\ =\left|{\displaystyle \underset{i=1}{\overset{q}{\sum }}}{\lambda }_i{\left(p_i^{\text{T}}\left(-{\mu }_1\right)\right)}^2-{\displaystyle \underset{i=1}{\overset{q}{\sum }}}{\lambda }_i{\left(p_i^{\text{T}}\left(Z-{\mu }_2\right)\right)}^2\right|\\ =\left|{\displaystyle \underset{i=1}{\overset{q}{\sum }}}{\lambda }_i{p}_i^{\text{T}}\left(2Z-{\mu }_1-{\mu }_2\right)p_i^{\text{T}}\left({\mu }_2-{\mu }_1\right)\right|\\ \le {\displaystyle \underset{i=1}{\overset{q}{\sum }}}{\lambda }_i{\Vert p_i\Vert }^2\Vert 2Z-{\mu }_1-{\mu }_2\Vert \Vert {\mu }_2-{\mu }_1\Vert \end{array}$

$\begin{array}{l}\le {\lambda }_1\Vert {\mu }_2-{\mu }_1\Vert \Vert 2Z-{\mu }_1-{\mu }_2\Vert {\displaystyle \underset{i=1}{\overset{q}{\sum }}}{\Vert p_i\Vert }^2\\ =\rho \left({\Sigma }_2^{-1}\right)\Vert {\mu }_2-{\mu }_1\Vert \Vert 2Z-{\mu }_1-{\mu }_2\Vert {\displaystyle \underset{i=1}{\overset{q}{\sum }}}{\Vert p_i\Vert }^2\\ \le \Vert {\Sigma }_2^{-1}\Vert \Vert {\mu }_2-{\mu }_1\Vert \Vert 2Z-{\mu }_1-{\mu }_2\Vert ,\end{array}$ (15)

where $\rho \left({\Sigma }_2^{-1}\right)$ is the spectral radius of ${\Sigma }_2^{-1}$ and $\rho \left(A\right)\le \Vert A\Vert$ applicable to any induced norm. Similarly ${\Sigma }_2^{-1}-{\Sigma }_1^{-1}$ is also symmetric, thus indicate the eigenvalues of ${\Sigma }_2^{-1}-{\Sigma }_1^{-1}$ as $p_1^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}\cdots \mathrm{,}{p}_q^{\mathrm{<mo>\; *\; </mo>}}$ s.t. $\left({\Sigma }_2^{-1}-{\Sigma }_1^{-1}\right)p_i^{\mathrm{<mo>\; *\; </mo>}}={\lambda }_i^{\mathrm{<mo>\; *\; </mo>}}{p}_i^{\mathrm{<mo>\; *\; </mo>}}$ . There exists an orthogonal matrix in matrix format $Q^{\mathrm{<mo>\; *\; </mo>}}$ s.t.

${\left(Q^{\mathrm{<mo>\; *\; </mo>}}\right)}^{-1}\left({\Sigma }_2^{-1}-{\Sigma }_1^{-1}\right)Q^{\mathrm{<mo>\; *\; </mo>}}={\left(Q^{\mathrm{<mo>\; *\; </mo>}}\right)}^{\text{T}}\left({\Sigma }_2^{-1}-{\Sigma }_1^{-1}\right)Q^{\mathrm{<mo>\; *\; </mo>}}={\Lambda }^{\mathrm{<mo>\; *\; </mo>}}\mathrm{.}$

Next, we have

$\begin{array}{l}\left|{\left(Z-{\mu }_1\right)}^{\text{T}}{\Sigma }_2^{-1}\left(Z-{\mu }_1\right)-{\left(Z-{\mu }_1\right)}^{\text{T}}{\Sigma }_1^{-1}\left(Z-{\mu }_1\right)\right|\\ =\left|{\left(Z-{\mu }_1\right)}^{\text{T}}\left({\Sigma }_2^{-1}-{\Sigma }_1^{-1}\right)\left(Z-{\mu }_1\right)\right|\\ =\left|{\left(Z-{\mu }_1\right)}^{\text{T}}{Q}^{*}{\Lambda }^{*}{\left(Q^{*}\right)}^{\text{T}}\left(Z-{\mu }_1\right)\right|\\ =\left|{\left[{\left(Q^{*}\right)}^{\text{T}}\left(Z-{\mu }_1\right)\right]}^{\text{T}}{\Lambda }^{*}\left[{\left(Q^{*}\right)}^{\text{T}}\left(Z-{\mu }_1\right)\right]\right|\end{array}$

$\begin{array}{l}=\left|{\displaystyle \underset{i=1}{\overset{q}{\sum }}}{\lambda }_i{\left({\left(p_i^{*}\right)}^{\text{T}}\left(Z-{\mu }_1\right)\right)}^2\right|\\ \le \underset{i=1,\cdots ,q}{\max }\left|{\lambda }_i^{*}\right|{\displaystyle \underset{i=1}{\overset{q}{\sum }}}{\left({\left(p_i^{*}\right)}^{\text{T}}\left(Z-{\mu }_1\right)\right)}^2\\ \le \rho \left({\Sigma }_2^{-1}-{\Sigma }_1^{-1}\right){\Vert Z-{\mu }_1\Vert }^2{\displaystyle \underset{i=1}{\overset{q}{\sum }}}\Vert {\left(p_i^{*}\right)}^2\Vert \\ \le {\Vert Z-{\mu }_1\Vert }^2\Vert {\Sigma }_2^{-1}-{\Sigma }_1^{-1}\Vert .\end{array}$ (16)

[23] presented a simple and popular M-estimator that minimizes a “bounded” version of the sum of residuals. The estimating equation is

${\hat{\theta }}_M={\displaystyle \underset{i=1}{\overset{n}{\sum }}}\psi \left(\frac{r_i\left(\hat{\theta }\right)}{\sigma }\right)x_i=0$ (17)

Then, using (12), (13), (14), (15) and (17) we get

$\begin{array}{l}\Vert w_{{\eta }_1}\left(x\right){\varphi }_{{\beta }_2}^i\left(x,y\right)-w_{{\eta }_2}\left(x\right){\varphi }_{{\beta }_2}\left(x,y\right)\Vert \\ <K_1\left|{\alpha }_2-{\alpha }_1\right|{\left(Z-{\mu }_2\right)}^{\text{T}}{\Sigma }_2^{-1}\left(Z-{\mu }_2\right)\left|x_i\right|\\ +K_1\Vert {\Sigma }_2^{-1}\Vert \Vert {\mu }_2-{\mu }_1\Vert \Vert 2Z-{\mu }_1-{\mu }_2\Vert \Vert x_i\Vert \\ +K_1{\Vert Z-{\mu }_1\Vert }^2\Vert {\Sigma }_2^{-1}-{\Sigma }_1^{-1}\Vert \left|x_i\right|.\end{array}$ (18)

Since $\left|{\alpha }_2-{\alpha }_1\right|,\Vert {\mu }_2-{\mu }_1\Vert ,\Vert {\beta }_2-{\beta }_1\Vert <\Vert {\gamma }_2-{\gamma }_1\Vert$ , from (11), (10) and (16) now, we can provide a bound for $\psi$ :

$\begin{array}{l}\left|{\psi }_{{\gamma }_1}^i\left(x,y\right)-{\psi }_{{\gamma }_2}\left(x,y\right)\right|\\ <K_0\Vert {\beta }_2-{\beta }_1\Vert \left|x_i\right|\Vert x\Vert +K_1\left|{\alpha }_2-{\alpha }_1\right|{\left(Z-{\mu }_2\right)}^{\text{T}}{\Sigma }_2^{-1}\left(Z-{\mu }_2\right)\left|x_i\right|\\ +K_1\Vert {\Sigma }_2^{-1}\Vert \Vert {\mu }_2-\mu +1\Vert \Vert 2Z-{\mu }_1-{\mu }_2\Vert \left|x_i\right|+K_1{\Vert Z-{\mu }_1\Vert }^2\Vert {\Sigma }_2^{-1}-{\Sigma }_1^{-1}\Vert \left|x_i\right|\\ \le (K_0\left|x_i\right|\Vert x\Vert +K_1{\left(Z-{\mu }_2\right)}^{\text{T}}{\Sigma }_2^{-1}\left(Z-{\mu }_2\right)\left|x_i\right|\\ +K_1\Vert {\Sigma }_2^{-1}\Vert \Vert 2Z-{\mu }_1-{\mu }_2\Vert \left|x_i\right|+K_1{\Vert Z-{\mu }_1\Vert }^2\left|x_i\right|)\Vert {\gamma }_2-{\gamma }_1\Vert \\ =L^i\left(x\right)\Vert {\gamma }_2-{\gamma }_1\Vert ,\end{array}$

for every ${\gamma }_1\mathrm{,}{\gamma }_2$ . We have constructed a Lipschitz criterion for each ${\psi }^i,i=1,\cdots ,p$ , so for $\psi$ we have $\left|{\psi }_{{\gamma }_1}\left(x,y\right)-{\psi }_{{\gamma }_2}\left(x,y\right)\right|<L\left(x\right)\Vert {\gamma }_2-{\gamma }_1\Vert$ , every ${\gamma }_1\mathrm{,}{\gamma }_2$ ,

where

$\begin{array}{c}L\left(x\right)=K_{0}x\Vert x\Vert +K_1{\left(Z-{\mu }_2\right)}^{\text{T}}{\Sigma }_2^{-1}\left(Z-{\mu }_2\right)x\\ +K_1\Vert {\Sigma }_2^{-1}\Vert \Vert 2Z-{\mu }_1-{\mu }_2\Vert \left|x_i\right|+K_1{\Vert Z-{\mu }_1\Vert }^{2}x\mathrm{.}\end{array}$

Now investigate the bracketing entropy in relative to the $L_r\left(P\right)$ -norm

${\Vert {\psi }_i\Vert }_{Pr}={\left(P{\left|\psi \right|}^r\right)}^{\frac{1}{r}}\mathrm{.}$

Use brackets of the type $\left[{\psi }_{\gamma }-\epsilon L\mathrm{,}{\psi }_{\gamma }+\epsilon L\right]$ for $\gamma$ spanning over a sufficient chosen subset of $\Gamma$ and these brackets have $L_r\left(P\right)-size2\epsilon {\Vert L\Vert }_{pr}$ . If $\gamma$ ranges over a grid of mesh width $\epsilon$ over $\Gamma$ then the brackets $\left[{\psi }_{\gamma }-\epsilon L\mathrm{,}{\psi }_{\gamma }+\epsilon L\right]$ range over $Ϝ$ . Based on the Lipschitz condition we get

${\psi }_{\gamma 1}-\epsilon L\le {\psi }_{\gamma 2}\le {\psi }_{\gamma 1}+\epsilon L\mathrm{,}\text{if}\Vert {\gamma }_2-{\gamma }_1\Vert \le \epsilon \mathrm{,}$

Hence, many brackets are required as in radial balls $\frac{\epsilon }{2}$ to cover $\Gamma$ , alternatively we require less than ${\left(diam\Gamma /\epsilon \right)}^{2p+2}$ cubes with size $\epsilon$ to cover parameter space $\Gamma$ . If $P{\left|L_i\right|}^r<\infty$ , then a constant J exists, depends only on $\Gamma$ and P, in such a way that the bracketing numbers fulfill

$N_{\mathrm{<mo\; stretchy="false">\; [\; </mo>\; <mo\; stretchy="false">\; ]\; </mo>}}\left(\epsilon \mathrm{,}Ϝ\mathrm{,}{L}_R\left(P\right)\right)\le J{\left(\frac{diam\Gamma }{\epsilon }\right)}^p\mathrm{,}\text{every}0<\epsilon <diam\Gamma \mathrm{.}$

Given that all $\psi \in Ϝ$ are continuous functions, they can be measured. If B3 is fulfilled, then $P\left|L\right|<\infty$ , and as a result, class $Ϝ$ is P-Glivenko-Cantelli from the Theorem (19.4) (Glivenko-Cantelli) in [24] .¨

If B1, B2, B3 and B4 are hold, then estimators ${\hat{\beta }}_n$ are used as the solving to the equation $\Psi \left({\hat{\beta }}_n\right)=0$ converges to ${\beta }_0$ in probability.

Proof of Theorem 1. denote

$\Psi \left(\beta \right)=P{\psi }_{\beta ,0}=P{\varphi }_{\beta }\Psi \left(\beta ,\eta \right)=P{\psi }_{\beta ,\eta }=P{\varphi }_{\beta }{W}_{\eta }$

$\Psi \left(\beta \right)={\Psi }_n\left(\beta ,\eta \right)=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\psi }_{\beta ,\eta }{\Psi }_n\left(\beta ,\hat{\eta }\right)=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\psi }_{\beta ,\hat{\eta }}.$

Note that $\Psi \left({\hat{\beta }}_n\right)=\left(\Psi \left({\hat{\beta }}_n\right)-{\Psi }_n\left({\hat{\beta }}_n\right)\right)+{\Psi }_n\left({\hat{\beta }}_n\right)$ . Then we can establish that

${\sup }_{\beta \in \Theta }\Vert {\Psi }_n\left(\beta \right)-\Psi \left(\beta \right)\Vert =0_p\left(1\right)$ , then $\Vert \Psi \left({\hat{\beta }}_n\right)\Vert =0_p\left(1\right)$ since ${\Psi }_n\left({\hat{\beta }}_n\right)=0_p\left(1\right)$ . Thus, based upon Lemma 2, we have $\Vert {\hat{\beta }}_n-{\beta }_0\Vert =0_p\left(1\right)$ . To show ${\sup }_{\beta \in \Theta }\Vert {\Psi }_n\left(\beta \right)-\Psi \left(\beta \right)\Vert =0_p\left(1\right)$ , we consider it as follows:

$\begin{array}{l}\underset{\beta \in \Theta }{\sup }\Vert {\Psi }_n\left(\beta \mathrm{,}\hat{\eta }\right)-\Psi \left(\beta \right)\Vert \\ =\underset{\beta \in \Theta }{\sup }\Vert {\Psi }_n\left(\beta \mathrm{,}\hat{\eta }\right)+\Psi \left(\beta \mathrm{,}\eta \right)+{\Psi }_n\left(\beta \mathrm{,}\eta \right)-\Psi \left(\beta \mathrm{,}\eta \right)+\Psi \left(\beta \mathrm{,}\eta \right)-\Psi \left(\beta \mathrm{,0}\right)\Vert \\ \le J_1+J_2+J_3\mathrm{,}\end{array}$

where

$J_1=\underset{\beta \in \Theta }{\sup }\Vert {\Psi }_n\left(\beta ,\hat{\eta }\right)-{\Psi }_n\left(\beta ,\eta \right)\Vert$

$J_2=\underset{\beta \in \Theta }{\sup }\Vert {\Psi }_n\left(\beta ,\eta \right)-\Psi \left(\beta ,\eta \right)\Vert$

$J_3=\underset{\beta \in \Theta }{\sup }\Vert \Psi \left(\beta ,\eta \right)-\Psi \left(\beta ,0\right)\Vert .$

Lemma 3 informs us that $Ϝ$ is P-Glivenko-Cantelli class, so $J_2\to 0$ . For $J_1$ ,

$\begin{array}{c}{J}_1=\underset{\beta \in \Theta }{\sup }\Vert \frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\varphi }_{\beta }\left(x_i,y_i\right)\left(w_{\hat{\eta }}\left(x_i\right)-w_{\eta }\left(x_i\right)\right)\Vert \\ \le \frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\underset{\beta \in \Theta }{\sup }\Vert {\varphi }_{\beta }\left(x_i,y_i\right)\left(w_{\hat{\eta }}\left(x_i\right)-w_{\eta }\left(x_i\right)\right)\Vert \\ =\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}\left|w_{\hat{\eta }}\left(x_i\right)-w_{\eta }\left(x_i\right)\right|\underset{\beta \in \Theta }{\sup }\Vert {\varphi }_{\beta }\left(x_i,y_i\right)\Vert .\end{array}$

Using (16), we get

$\begin{array}{c}\left|w_{\hat{\eta }}\left(x_i\right)-w_{\eta }\left(x_i\right)\right|<K_1\left|\hat{\alpha }-\alpha \right|{\left(Z-\mu \right)}^{\text{T}}{\Sigma }^{-1}\left(Z-\mu \right)\\ +K_1\Vert {\Sigma }^{-1}\Vert \Vert \hat{\mu }-\mu \Vert \Vert 2Z-\hat{\mu }-\mu \Vert \\ +K_1{\Vert Z-\hat{\mu }\Vert }^2\Vert {\hat{\Sigma }}^{-1}-{\Sigma }^{-1}\Vert \mathrm{.}\end{array}$ (19)

Since $\hat{\eta }-\eta =\left(\hat{\alpha }-\alpha ,\hat{\mu }-\mu ,\hat{\Sigma }-\Sigma \right)=0_p\left(1\right)$ , we have $\hat{\alpha }-\alpha =0_p\left(1\right)$ , $\hat{\mu }-\mu =0_p\left(1\right)$ and $\hat{\Sigma }-\Sigma =0_p\left(1\right)$ . Also, ${\hat{\Sigma }}^{-1}-\Sigma ={\hat{\Sigma }}^{-1}\left(\Sigma -\hat{\Sigma }\right){\Sigma }^{-1}$ , as a result ${\hat{\Sigma }}^{-1}-{\Sigma }^{-1}=0_p\left(1\right)$ . Therefore, (19) follows that $\left|w_{\hat{\eta }}\left(X_i\right)\right|=0_p\left(1\right)$ . Because B2 and B3 are fulfilled, we also have ${\sup }_{\beta \in \Theta }{\Vert {\varphi }_{\beta }\left(x,y\right)\Vert }^2<\infty$ , then we get $J_1{\to }^{p}0$ . For $J_3$

$\begin{array}{c}{J}_3=\underset{\beta \in \Theta }{\sup }\Vert P{\varphi }_{\beta }\left(x,y\right)\left(w_{\eta }\left(x\right)-w_{{\eta }_0}\left(x\right)\right)\Vert \\ \le \underset{\beta \in \Theta }{\sup }P\Vert {\varphi }_{\beta }\left(x,y\right)\left(w_{\eta }\left(x\right)-w_{{\eta }_0}\left(x\right)\right)\Vert \\ \le \underset{\beta \in \Theta }{\sup }{P}^{\frac{1}{2}}{\Vert {\varphi }_{\beta }\left(x,y\right)\Vert }^2{P}^{\frac{1}{2}}{\left|w_{\eta }\left(x\right)-w_{{\eta }_0}\left(x\right)\right|}^2\\ =P^{\frac{1}{2}}{\left|w_{\eta }\left(x\right)-w_{{\eta }_0}\left(x\right)\right|}^2\underset{\beta \in \Theta }{\sup }{P}^{\frac{1}{2}}{\Vert {\varphi }_{\beta }\left(x,y\right)\Vert }^2,\end{array}$

and the Cauchy-Schwarz inequality was employed in the final inequality. We already know

$w_{\eta }\left(x\right)-w_{{\eta }_0}\left(x\right)=1-m\left(\eta {\left(Z-\mu \right)}^{\text{T}}{\Sigma }^{-1}\left(Z-\mu \right)\right)\to 0.$

Then $\hat{\eta }=0_p\left(1\right)m\left(0\right)=1$ and ${\left|w_{\hat{\eta }}\left(x\right)-w_{{\eta }_0}\left(x\right)\right|}^2$ is bounded by 1, using the dominated convergence theorem we get

$\underset{n\to \infty }{lim}{P}^{\frac{1}{2}}{\left|w_{\eta }\left(x\right)-w_{{\eta }_0}\left(x\right)\right|}^2\to 0.$

We also have

$\underset{\beta \in \Theta }{\sup }P{\Vert {\varphi }_{\beta }\left(x,y\right)\Vert }^2=\underset{\beta \in \Theta }{\sup }P{\Vert \left(y-x^{\text{T}}\beta \right)x\Vert }^2\le P{\Vert x\Vert }^2.$

and from assumption B3, we obtain ${\sup }_{\beta \in \Theta }P{\Vert {\varphi }_{\beta }\left(x\mathrm{,}y\right)\Vert }^2<\infty$ . It then follows that $J_3\to 0$ .¨

5. Assess Robustness of the Estimators

To evaluate the robustness of the methods, two methods have been employed. The first includes simulated models for contrasting the novel technique with the traditional MLE, Mallows-type estimator (Mallows) and (CUBI). In the second, we used the actual leukemia data set and The Erythrocyte Sedimentation Rate Data.

5.1. Simulation study

A Monte Carlo simulation analysis was performed for this subsection to assess the efficiency and robustness of the suggested estimator ${\hat{\beta }}_n$ . For the initial robust estimators of the scatter and location $\hat{\Sigma }$ and $\hat{\mu }$ utilized the minimum covariance determinant (MCD). We calculated the following estimators for comparison: MLE, conditionally unbiased bounded influence (CUBI) of [2] , and the Mallows-type estimator (Mallows) of [3] . In the simulation, the following weight function $m\left(x\right)$ was applied:

$m\left(x\right)=a\left(1-{\left(x-1\right)}^6\right)+b\mathrm{,}$

where $a=0.8$ and $b=0.2$ .

Three models are involved in the simulation study: a clean logistic regression model, a contaminated model with a 10% contamination rate, and a contaminated model with a 20% contamination rate. First, clean model, the standard normal distribution was used to generate two predictor variables, $x_1~N\left(\mathrm{0,1}\right)$ and $x_2~N\left(\mathrm{0,1}\right)$ . Three sample sizes were used: ( $n=100,300,500$ ) and $p=2$ . We generate the response variable according to the Bernoulli distribution with a parameter equal to ${\pi }_i=\exp \left({\beta }_0+{\beta }_1{x}_1+{\beta }_2{x}_2\right)/\left(1+\exp \left({\beta }_0+{\beta }_1{x}_1+{\beta }_2{x}_2\right)\right)$ . The values of the true parameters $\beta$ are taken $\left(\mathrm{0,1.6,1.2}\right)$ for three models. Second, the percentage of contamination in data equal 10%, and their predicted variables are generated from a normal distribution with ( $\mu =0$ ) and ( $\sigma =1$ ). Third, the percentage of contamination in data equals 20% in a similar manner to the above.

The performance of these estimators is evaluated using Bias and mean squared error for various models. Nevertheless, the estimator with the smallest Bias and MSE is the best. Each scenario simulation included over 1000 repetitions. Consequently, for each parameter, the following are the calculations for bias and mean squared error:

$\text{Bias}=\left|\frac{1}{1000}{\displaystyle \underset{i=1}{\overset{1000}{\sum }}{\beta }_i-\beta }\right|\mathrm{,}$

and

$\text{MSE}=\frac{1}{1000}{\displaystyle \underset{i=1}{\overset{n}{\sum }}}{\left|{\hat{\beta }}_i-\beta \right|}^2\mathrm{.}$

5.2. Numerical Results

The numerical results, displayed in this paper, are based on simulation studies, and two real data applications. This numerical result is expected to evaluate, the performance of a proposed model. Table 1 shows the bias and mean squared errors of the four estimation techniques for the clean model. The findings show that the bias and MSE of the MLE, Mallows, and CUBIF estimators are relatively similar, while the WMLT estimator performs worse than the other estimators. When the sample size increases, the bias and mean squared errors are observed to decrease. As shown in Table 2, under 10% of the data were contaminated, so the new robust approach WMLT has the greatest overall performance among all comparable estimators for varied sample sizes. Table 3 demonstrates that when 20% of the data are contaminated, our proposed technique (WMLT) outperforms other estimators in terms of bias and mean squared errors. Due to the sensitivity of anomalies, conventional maximum likelihood estimates perform inadequately in the contaminated model. In conclusion, the proposed method outperforms all other methods compared with contaminated data. Furthermore, the new estimator performs reasonably well with clear data.

Table 1. Bias and mean squared errors of estimators for clean model.

Table 2. Bias and (MSE) of estimators for second model (10% of data are contaminated).

Table 3. Bias and (MSE) of estimators for second model (20% of data are contaminated).

5.3. Leukemia Data

This study uses data from [25] which includes information from 33 people who died from acute myeloid leukemia. Each patient was measured for three variables: AG, WBC, and time. The response variable represents the survival time in weeks of the patient; it was converted into a binary variable with (Y = 1) indicating patients whose survival time exceeded 52 weeks and (Y = 0) indicating those who did not. WBC represents the white blood cell concentration of the patient. Whereas AG (present = 1, absent = 0) measured the presence or absence of a morphologic characteristic of white blood cells. The observation number 17 appears to be atypical. Using AG and WBC as predictor variables and binary survival time y as the response variable, a logistic regression model was constructed. The estimators analyzed here are weighted maximum likelihood technique (WMLT), MLE, MLE17 (MLE17 is the maximum likelihood estimator for clean data when observation 17 is excluded.), Mallows (estimator of the Mallows type), and CUBIF (conditionally unbiased bounded-influence function estimator).

Table 4 demonstrates that the MLE is extremely sensitive to influential observations. Furthermore, eliminating observation 17 lowered the impact of WBC to near nil. For the leukemia data, the new estimator (WMLT) demonstrated the greatest performance among all other estimators. However, Mallow’s estimators are reasonably close to the MLE17.

5.4. The Erythrocyte Sedimentation Rate Data

The Erythrocyte Sedimentation Rate (ESR) data. In this data, the primary objective was to determine if the levels of two plasma proteins (Fibrinogen and γ.Globulin) were responsible for the increase in ESR in healthy individuals. The research was conducted by the Institute of Medical Research in Kuala Lumpur, Malaysia, on 32 patients, and the original data was collected by [26] . The response of zero indicates a healthy person, whereas the response of one indicates an unwell person. Here, the continuous variables (FIB and γ.GLO) are compared to the binary response (ESR). In the original ESR data, [27] identified two outliers (cases 13 and 29) in X-space. Cases 14 and 15 are influential observations. Therefore, removing instances 14 and 15 would result in cases with no overlap.

Table 4. The estimated parameters and standard errors for the leukemia data.

Case 13 was modified so that $\left(y,x_1,x_2\right)=\left(3,6,37\right)$ in order to execute uncontaminated data analysis. From the uncontaminated data, the ESR data were contaminated where the occurrences (Y = 1) and non-occurrences (Y = 0) were replaced with each other for cases 14 and 15, and this might only leave one out of the three overlapping cases for the ESR data.

Under contaminated ESR data, β₀ and se(β₀) of all estimators are impacted by outliers as compared to the other parameters (see Table 5). The results shown in Table 5 also indicate that the MLE is primarily influenced by outliers. Following the modification of the tainted data, only one overlapping observation, case 13, remains. This is the reason why the coefficients and standard errors of the WMLT that downweights this observation are so large. Even though the WMLT has the smallest χ² value, the WMLT estimator should also be taken into consideration. The results shown in Table 5 indicate that the WMLT is a suitable estimator for the ESR data because its estimates are relatively closer to the MLE for uncontaminated data.

6. Conclusions

In this paper, we develop a new robust technique for logistic regression, also known as the weighted maximum likelihood technique (WMLT). The asymptotic consistency of the proposed estimator was demonstrated.

In order to evaluate the robustness of a new technique, we conducted simulation experiments using a variety of scenarios and data sets. Classical maximum likelihood estimates show the lack of robustness in the presence of outliers. Our simulation study for the clean model illustrated that the MLE, Mallows, and CUBIF estimators perform similarly, while the new weighted technique performs less effectively than the other estimators. The simulation study also shows that the WMLT technique outperforms other estimators when dealing with contaminated data and demonstrates the greatest performance among all estimators compared to various scenarios and real data sets. The proposed method (WMLT) can be applied to other generalized linear models (GLMs) and is expected to be superior to existing methods in practical applications. The findings

Table 5. Estimated coefficients, standard errors, and χ² for ESR.

of this chapter provide researchers and practitioners with a new approach to developing robust estimators for logistic regression and potentially other generalized linear models (GLMs).

Acknowledgments

We thank the Editor and the referee for their comments. Research of M. Baron is funded by the National Science Foundation grant DMS 1322353. This support is greatly appreciated.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Pregibon, D. (1981) Logistic Regression Diagnostics. Annals of Statistics, 4, 705-724. https://doi.org/10.1214/aos/1176345513
[2]	Kunsch, H.R., Stefonski, L.A. and Carroll, R.J. (1989) Conditionally Unbiased Bounded-Influence Estimation in General Regression Models, with Applications to Generalized Linear Models. Journal of the American Statistical Association, 84, 460-466. https://doi.org/10.1080/01621459.1989.10478791
[3]	Carroll, R.J. and Pederson, S. (1993) On Robustness in the Logistic Regression Model. Journal of the Royal Statistical Society, 55, 693-706. https://doi.org/10.1111/j.2517-6161.1993.tb01934.x
[4]	Hobza, T., Mart, N. and Pardo, L. (2017) A Wald-Type Test Statistic Based on Robust Modified Median Estimator in Logistic Regression Models. Journal of Statistical Computation and Simulation, 87, 2309-2333. https://doi.org/10.1080/00949655.2017.1330414
[5]	Bergesio, A. and Yohai, V.J. (2011) Projection Estimators for Generalized Linear Models. Journal of the American Statistical Association, 106, 661-671. https://doi.org/10.1198/jasa.2011.tm09774
[6]	Myers, M.D.C.R.H. and Vining, G.G. (1992) Least-Absolute-Deviations Fits for Generalized Linear Models. Biometrika, 79, 747-754. https://doi.org/10.1093/biomet/79.4.747
[7]	Cantoni, E. and Ronchetti, E. (2001) Robust Inference for Generalized Linear Models. Journal of the American Statistical Association, 96, 1022-1030. https://doi.org/10.1198/016214501753209004
[8]	Wedderburn, R.W.M. (1974) Quasi-Likelihood Functions, Generalized Linear Models, and the Gaussian Newton Method. Biometrika, 61, 439-447. https://doi.org/10.1093/biomet/61.3.439
[9]	Croux, C., Flandre, C. and Haesbroeck, G. (2002) The Breakdown Behavior of the Maximum Likelihood Estimator in the Logistic Regression Model. Statistics and Probability Letters, 60, 377-386. https://doi.org/10.1016/S0167-7152(02)00292-4
[10]	Rousseeuw, P.J. and Christmann, A. (2005) Robustness against Separation and Outliers in Logistic Regression. Quality Control and Applied Statistics, 50, 451-452.
[11]	Idriss, I.A. and Cheng, W. (2023) Robust Estimators for Poisson Regression. Open Journal of Statistics, 13, 112-118. https://doi.org/10.4236/ojs.2023.131007
[12]	Valdora, M. and Yohai, V.J. (2014) Robust Estimators for Generalized Linear Models. Journal of Statistical Planning and Inference, 146, 31-48. https://doi.org/10.1016/j.jspi.2013.09.016
[13]	Bianco, A.M. and Yohai, V.J. (1996) Robust Estimation in the Logistic Regression Model. Springer, New York, 17-34. https://doi.org/10.1007/978-1-4612-2380-1_2
[14]	Croux, C. and Haesbroeck, G. (2014) Implementing the Bianco and Yohai Estimator for Logistic Regression. Computational Statistics and Data Analysis Journal, 146, 31-48.
[15]	Muller, C.H. and Neykov, N. (2003) Breakdown Points of Trimmed Likelihood Estimators and Related Estimators in Generalized Linear Models. Journal of Statistical Planning and Inference, 44, 273-295. https://doi.org/10.1007/978-3-642-57338-5_24
[16]	Copas, J.B. (1988) Binary Regression Models for Contaminated Data. Journal of the Royal Statistical Society, 50, 225-265. https://doi.org/10.1111/j.2517-6161.1988.tb01723.x
[17]	Ahmed, I.A.I. and Cheng, W.H. (2020) The Performance of Robust Methods in Logistic Regression Model. Open Journal of Statistics, 10, 127-138. https://doi.org/10.4236/ojs.2020.101010
[18]	Ruckstuhl, A.F. and Welsh, A.H. (2001) Robust Fitting of the Binomial Model. The Annals of Statistics, 29, 1117-1136. https://doi.org/10.1214/aos/1013699996
[19]	Krasker, W.S. and Welsch, R.E. (1982) Efficient Bounded-Influence Regression Estimation. Journal of the American Statistical Association, 77, 595-604. https://doi.org/10.1080/01621459.1982.10477855
[20]	Stefanski, L.A., Carroll, R.J. and Ruppert, D. (1986) Optimally Bounded Score Functions for Generalized Linear Models with Applications to Logistic Regression. Biometrika, 73, 413-424. https://doi.org/10.1093/biomet/73.2.413
[21]	Gervini, D. (2005) Robust Adaptive Estimators for Binary Regression Models. Journal of Statistical Planning and Inference, 131, 297-311. https://doi.org/10.1016/j.jspi.2004.02.006
[22]	Vandev, D.L. and Neykov, N.M. (1998) About Regression Estimators with High Breakdown Point. Statistics, 32, 111-129. https://doi.org/10.1080/02331889808802657
[23]	Huber, P.J. (1964) Robust Estimation of a Location Parameter. The Annals of Mathematical Statistics, 35, 73-101. https://doi.org/10.1214/aoms/1177703732
[24]	Vaart, A.W. (2000) Asymptotic Statistics. Cambridge University Press, Cambridge.
[25]	Cox, D. and Oakes, D. (1984) Analysis of Survival Data. Chapman and Hall, New York.
[26]	Collett, D. and Jemain, A. (1985) Residuals, Outliers and Influential Observations in Regression Analysis. Sains Malaysiana, 14, 493-511.
[27]	Syaiba, B. and Habshah, M. (2010) Robust Logistic Diagnostic for the Identification of High Leverage Points in Logistic Regression Model. Journal of Applied Sciences, 10, 3042-3050. https://doi.org/10.3923/jas.2010.3042.3050