A robust approach for generalized linear models based on maximum L$q$-likelihood procedure

The maximum L$q$-likelihood estimator for $\beta$ in the class of generalized linear models is a solution of the estimating equation:

with

corresponds to the weight arising from the estimation function defined in (1.4), with $c_{q}(y_{i},\phi)=(1-q)c(y_{i},\phi)$, while $\mu_{i}=\operatorname{E}(y_{i})=\dot{b}(\theta_{i})$, $\operatorname{var}(y_{i})=\phi^{-1}V_{i}$, with $V_{i}=\ddot{b}(\theta_{i})$ being the variance function, and $W_{i}=V_{i}\{\dot{k}(\eta_{i})\}^{2}$, we shall use dots over functions to denote derivatives, $\dot{k}(u)={\operatorname{d}}k/{\operatorname{d}}u$ and $\ddot{k}(u)={\operatorname{d}}^{2}k/{\operatorname{d}}u^{2}$. In this paper, we consider $q<1$ in which case the weights $U_{i}(\mbox{\boldmath$\beta$})$ $(i=1\dots,n)$ provide a mechanism for downweighting outlying observations. The individual score function adopts the form $\mbox{\boldmath$s$}_{i}(\mbox{\boldmath$\beta$})=\phi W_{i}^{1/2}r_{i}\mbox{\boldmath$x$}_{i}$, where $r_{i}=(y_{i}-\mu_{i})/\sqrt{V_{i}}$ denotes the Pearson residual. It is straightforward to notice that the estimation function for $\beta$ can be written as:

where $\mbox{\boldmath$U$}=\operatorname{diag}(U_{1},\dots,U_{n})$, with $U_{i}=U_{i}(\mbox{\boldmath$\beta$})$ is given in (2.2), while $\mbox{\boldmath$W$}=\operatorname{diag}(W_{1},\dots,W_{n})$, $\mbox{\boldmath$V$}=\operatorname{diag}(V_{1},\dots,V_{n})$, $\mbox{\boldmath$Y$}=(y_{1},\dots,y_{n})^{\top}$ and $\mbox{\boldmath$\mu$}=(\mu_{1},\dots,\mu_{n})^{\top}$.

In Appendix A we show that the estimating function in (2.1) is unbiased for the surrogate parameter $\mbox{\boldmath$\beta$}_{*}$, that is $\operatorname{E}_{0}\{\mbox{\boldmath$\Psi$}_{n}(\mbox{\boldmath$\beta$}_{*})\}=\mbox{\boldmath$0$}$, where $\operatorname{E}_{0}(\cdot)$ denotes expectation with respect to the true distribution $f(\cdot;\theta_{0},\phi)$. Thereby, this allows introducing a suitable calibration function, say $\tau_{q}(\cdot)$, such that $\widehat{\mbox{\boldmath$\beta$}}_{q}=\tau_{q}(\widehat{\mbox{\boldmath$\beta$}}{}_{q}^{*})$ which is applied to rescale the parameter estimates in order to achieve Fisher-consistency for $\mbox{\boldmath$\beta$}_{0}$.

When partial derivatives and expectations exist, we define the matrices

where $\mbox{\boldmath$J$}=\operatorname{diag}(J_{1},\dots,J_{n})$ whose diagonal elements are given by $J_{i}=J_{q}^{-\phi}(\theta_{i})$, for $i=1,\dots,n$, with $J_{q}(\theta)=\exp\{-qb(\theta)+b(q\theta)\}$, $q\theta\in\Theta$ was defined by Menendez (2000), whereas $\mbox{\boldmath$G$}=\operatorname{diag}(\dot{g}(\theta_{1}^{*}),\dots,\dot{g}(\theta_{n}^{*}))$ and $\mbox{\boldmath$K$}=\operatorname{diag}(\dot{k}(\eta_{1}),\dots,\dot{k}(\eta_{n}))$ with ${\operatorname{d}}\eta_{*}/{\operatorname{d}}\eta=(1/q)\dot{g}(\theta_{*})\dot{k}(\eta)$ and $g(\theta_{*})=k^{-1}(\theta_{*})$ corresponding to the inverse function of the $\theta$-link evaluated on the surrogate parameter. We propose to use the Newton-scoring algorithm (Jørgensen and Knudsen, 2004) to address the parameter estimation associated with the estimating equation $\mbox{\boldmath$\Psi$}_{n}(\mbox{\boldmath$\beta$})=\mbox{\boldmath$0$}$, which assumes the form

Therefore, the estimation procedure adopts an iteratively reweighted least squares structure, as

where $\mbox{\boldmath$Z$}_{q}=\mbox{\boldmath$\eta$}+\mbox{\boldmath$J$}^{-1}\mbox{\boldmath$G$}^{-1}\mbox{\boldmath$K$}^{-1}\mbox{\boldmath$W$}^{-1/2}\mbox{\boldmath$UV$}^{-1/2}(\mbox{\boldmath$Y$}-\mbox{\boldmath$\mu$})$ denotes the working response which must be evaluated at $\mbox{\boldmath$\beta$}=\mbox{\boldmath$\beta$}^{(t)}$. Let $\widehat{\mbox{\boldmath$\beta$}}{}_{q}^{*}$ be the value obtained at the convergence of the algorithm in (2.5). Thus, in order to obtain a corrected version of the maximum L$q$-likelihood estimator, we must take

Evidently, for the canonical $\theta$-link, $k(u)=u$, it is sufficient to consider $\widehat{\mbox{\boldmath$\eta$}}_{q}=q\widehat{\mbox{\boldmath$\eta$}}{}_{q}^{*}$, or equivalently $\widehat{\mbox{\boldmath$\beta$}}_{q}=q\widehat{\mbox{\boldmath$\beta$}}{}_{q}^{*}$. In addition, for this case we have $\mbox{\boldmath$GK$}=\mbox{\boldmath$I$}_{n}$, with $\mbox{\boldmath$I$}_{n}$ being the identity matrix of order $n$, which leads to a considerably simpler version of the algorithm proposed in (2.5). Qin and Priebe (2017) suggested another alternative to achieve Fisher-consistency which is based on re-centering the estimation equation, considering a bias-correction term, and who termed it as the bias-corrected maximum L$q$-likelihood estimator. We emphasize the simplicity of the calibration function in (2.6), which introduces a reparameterization that allows us to characterize the asymptotic distribution of $\widehat{\mbox{\boldmath$\beta$}}_{q}$.

The following result describes the asymptotic normality of the $\widehat{\mbox{\boldmath$\beta$}}_{q}$ estimator in the context of generalized linear models.

From the general theory of robustness, we can characterize a slight misspecification of the proposed model by the influence function of a statistical functional $\mbox{\boldmath$T$}(F)$, which is defined as (see Hampel et al., 1986, Sec. 4.2),

where $F_{\epsilon}=(1-\epsilon)F+\epsilon\delta_{Y}$, with $\delta_{Y}$ a point mass distribution in $Y$. That is, $\operatorname{IF}(Y;\mbox{\boldmath$T$})$ reflects the effect on $T$ of a contamination by adding an observation in $Y$. It should be noted that the maximum L$q$-likelihood estimation procedure for $q$ fixed corresponds to an $M$-estimation method, where $q$ acts as a tuning constant. It is well known that the influence functions of the estimators obtained by the maximum likelihood and L$q$-likelihood methods adopt the form

respectively, where $\mbox{\boldmath$\mathcal{F}$}_{n}=\operatorname{E}\{-\partial\mbox{\boldmath$s$}_{n}(\mbox{\boldmath$\beta$})/\partial\mbox{\boldmath$\beta$}^{\top}\}=\phi^{-1}\mbox{\boldmath$X$}^{\top}\mbox{\boldmath$WX$}$ is the Fisher information matrix, with $\mbox{\boldmath$s$}_{n}(\mbox{\boldmath$\beta$})=\sum_{i=1}^{n}\mbox{\boldmath$s$}_{i}(\mbox{\boldmath$\beta$})$ the score function, whereas $\{f(y;k(\mbox{\boldmath$x$}^{\top}\mbox{\boldmath$\beta$}),\phi)\}^{1-q}$ corresponds to the weighting defined in (2.2) and $\mbox{\boldmath$s$}(\mbox{\boldmath$\beta$})=\phi W^{1/2}(y-\mu)\mbox{\boldmath$x$}/\sqrt{V}$, with $\mbox{\boldmath$x$}=(x_{1},\dots,x_{p})^{\top}$. It is possible to appreciate that when $q<1$ the weights $U_{i}(\mbox{\boldmath$\beta$})$ provide a mechanism to limit the influence of extreme observations. This allows downweighting those observations that are strongly discrepant from the proposed model. The properties of robust procedures based on the weighting of observations using the working statistical model have been discussed, for example in Windham (1995) and Basu et al. (1998), who highlighted the connection between these procedures with robust methods (see also, Ferrari and Yang, 2010). Indeed, we can write the asymptotic covariance of $\widehat{\mbox{\boldmath$\beta$}}_{q}$ equivalently as,

Additionally, it is known that $\mbox{\boldmath$B$}_{n}^{-1}\mbox{\boldmath$A$}_{n}\mbox{\boldmath$B$}_{n}^{-1}-\mbox{\boldmath$\mathcal{F}$}_{n}^{-1}$, is a positive semidefinite matrix. That is, the $\widehat{\mbox{\boldmath$\beta$}}_{q}$ estimator has asymptotic covariance matrix always larger than the covariance of $\widehat{\mbox{\boldmath$\beta$}}_{\mathrm{ML}}$ for the working model defined in (1.1) and (1.2). To quantify the loss in efficiency of the maximum L$q$-likelihood estimator, is to consider the asymptotic relative efficiency of $\widehat{\mbox{\boldmath$\beta$}}_{q}$ with respect to $\widehat{\mbox{\boldmath$\beta$}}_{\mathrm{ML}}$, which is defined as,

Hence, $\operatorname{ARE}(\widehat{\mbox{\boldmath$\beta$}}_{q},\widehat{\mbox{\boldmath$\beta$}}_{\mathrm{ML}})\geq 1$. This leads to a procedure for the selection of the distortion parameter $q$. Following the same reasoning as Windham (1995) an alternative is to choose $q$ from the available data as that $q_{\mathrm{opt}}$ value that minimizes $\operatorname{tr}(\mbox{\boldmath$B$}_{n}^{-1}\mbox{\boldmath$A$}_{n}\mbox{\boldmath$B$}_{n}^{-1})$. In other words, the interest is to select the model for which the loss in efficiency is the smallest. This procedure has also been used in Qin and Priebe (2017) for adaptive selection of the distortion parameter.

In this paper we follow the guidelines given by La Vecchia et al. (2015) who proposed a data-driven method for the selection of $q$ (see also Ribeiro and Ferrari, 2023). In this procedure the authors introduced the concept of stability condition as a mechanism for determining the distortion parameter by considering an equally spaced grid, say $1\geq q_{1}>q_{2}>\cdots>q_{m}$ and for each value in the grid $\widehat{\mbox{\boldmath$\beta$}}_{q_{j}}$, for $j=1,\dots,m$, is computed. To achieve robustness we can select an optimal value $q_{\mathrm{opt}}$, as the greatest value $q_{j}$ satisfying $QV_{j}<\rho$, where $QV_{j}=\|\widehat{\mbox{\boldmath$\beta$}}_{q_{j}}-\widehat{\mbox{\boldmath$\beta$}}_{q_{j+1}}\|$ with $\rho$ a threshold value. Based on experimental results, La Vecchia et al. (2015) proposed to construct a grid between $q_{1}=1$, and some minimum value, $q_{m}$ with increments of 0.01, whereas $\rho=0.05\,\|\widehat{\mbox{\boldmath$\beta$}}_{q_{m}}\|$. One of the main advantages of their procedure is that, in the absence of contamination, it allows to select the distortion parameter as close to 1 as possible, thas is, this leads to the maximum likelihood estimator. Yet another alternative that has shown remarkable performance for $q$-selection based on parametric bootstrap is reported by Giménez et al. (2022a, b).

The goodness of fit in generalized linear models, considering the estimation method based on the L$q$-likelihood function, requires appropriate modifications to obtain robust extensions of the deviance function. Consider the parameter vector to be partitioned as $\mbox{\boldmath$\beta$}=(\mbox{\boldmath$\beta$}_{1}^{\top},\mbox{\boldmath$\beta$}_{2}^{\top})^{\top}$ and suppose that our interest consists of testing the hypothesis $H_{0}:\mbox{\boldmath$\beta$}_{1}=\mbox{\boldmath$0$}$, where $\mbox{\boldmath$\beta$}_{1}\in\mathbb{R}^{r}$. Thus, we can evaluate the discrepancy between the model defined by the null hypothesis versus the model under $H_{1}:\mbox{\boldmath$\beta$}_{1}\neq\mbox{\boldmath$0$}$ using the proposal of Qin and Priebe (2017), as

where $\widehat{\mbox{\boldmath$\beta$}}_{q}$ and $\widetilde{\mbox{\boldmath$\beta$}}_{q}$ are the corrected maximum L$q$-likelihood estimates for $\beta$ obtained under $H_{0}$ and $H_{1}$, respectively, and using these estimates we have that $\widetilde{\mu}_{i}=\mu_{i}(\widetilde{\mbox{\boldmath$\beta$}}_{q})$ and $\widehat{\mu}_{i}=\mu_{i}(\widehat{\mbox{\boldmath$\beta$}}_{q})$, for $i=1,\dots,n$. Qin and Priebe (2017) based on results available in Cantoni and Ronchetti (2001), proved that the $D_{q}(\mbox{\boldmath$Y$},\widehat{\mbox{\boldmath$\mu$}})$ statistic has an asymptotic distribution following a discrete mixture of chi-squared random variables with one degree of freedom. We should note that the correction term proposed in Qin and Priebe (2017) in our case zeroed for the surrogate parameter.

Following Ronchetti (1997), we can introduce the Akaike information criterion based on the L$q$-estimation procedure, as follows

where $\widehat{\mbox{\boldmath$\beta$}}_{q}$ corresponds to the maximum L$q$-likelihood estimator obtained from the algorithm defined in (2.5). Using the definitions of $\mbox{\boldmath$A$}_{n}$ and $\mbox{\boldmath$B$}_{n}$, given in (2.3) and (2.4) leads to,

Evidently, when the canonical $\theta$-link is considered, we obtain that the penalty term assumes the form $2\operatorname{tr}(\widehat{\mbox{\boldmath$B$}}_{n}^{-1}\widehat{\mbox{\boldmath$A$}}_{n})=2p/(2-q)$ and in the case that $q\to 1$ we recover the usual Akaike information criterion, $\operatorname{AIC}$ for model selection in the context of maximum likelihood estimation. Some authors (see for instance, Ghosh and Basu, 2016), have suggested using the Akaike information criterion $\operatorname{AIC}_{q}$ as an alternative mechanism for the selection of the tuning parameter, i.e., the distortion parameter $q$.

A robust alternative for conducting hypothesis testing in the context of maximum L$q$-likelihood estimation has been proposed by Qin and Priebe (2017), who studied the properties of L$q$-likelihood ratio type statistics for simple hypotheses. As suggested in Equation (3.1) this type of development allows, for example, the evaluation of the fitted model. In this section we focus on testing linear hypotheses of the type

where $H$ es a known matrix of order $d\times p$, with $\operatorname{rank}(\mbox{\boldmath$H$})=d$ $(d\leq p)$ and $\mbox{\boldmath$h$}\in\mathbb{R}^{d}$. Wald, score-type and bilinear form (Crudu and Osorio, 2020) statistics for testing the hypothesis in (3.2) are given by the following result.

Robustness of the statistics defined in (3.3)-(3.5) can be characterized using the tools available in Heritier and Ronchetti (1994). In addition, it is relatively simple to extend these test statistics to manipulate nonlinear hypotheses following the results developed by Crudu and Osorio (2020).

Suppose that our aim is to evaluate the inclusion of a new variable in the regression model. Following Wang (1985), we can consider the added-variable model, which is defined through the linear predictor $\eta_{i}=\mbox{\boldmath$x$}_{i}^{\top}\mbox{\boldmath$\beta$}+z_{i}\gamma=\mbox{\boldmath$f$}_{i}^{\top}\mbox{\boldmath$\delta$}$, where $\mbox{\boldmath$f$}_{i}^{\top}=(\mbox{\boldmath$x$}_{i}^{\top},z_{i})$ and $\mbox{\boldmath$\delta$}=(\mbox{\boldmath$\beta$}^{\top},\gamma)^{\top}$, with $\theta_{i}=k(\eta_{i})$, for $i=1,\dots,n$. Based on the above results we present the score-type statistic for testing the hypothesis $H_{0}:\gamma=0$. Let $\mbox{\boldmath$F$}=(\mbox{\boldmath$X$},\mbox{\boldmath$z$})$ be the model matrix for the added-variable model, with $\mbox{\boldmath$X$}=(\mbox{\boldmath$x$}_{1},\dots,\mbox{\boldmath$x$}_{n})^{\top}$ and $\mbox{\boldmath$z$}=(z_{1},\dots,z_{n})^{\top}$. That notation, enables us to write the estimation function associated with the maximum L$q$-likelihood estimation problem as:

whereas,

It is easy to notice that the estimate for $\beta$ under the null hypothesis, $H_{0}:\gamma=0$ can be obtained by the Newton-scoring algorithm given in (2.5). Details of the derivation of the score-type statistic for testing $H_{0}:\gamma=0$ are deferred to Appendix E, where it is obtained

which for the context of the added variable model, takes the form:

where $\mbox{\boldmath$P$}=\mbox{\boldmath$X$}(\mbox{\boldmath$X$}^{\top}\mbox{\boldmath$WJGKX$})^{-1}\mbox{\boldmath$X$}^{\top}$. Thus, we reject the null hypothesis by comparing the value obtained for $R_{n}$ with some percentile value $(1-\alpha)$ of the chi-squared distribution with one degree of freedom. This type of development has been used, for example, in Wei and Shih (1994) to propose tools for outlier detection considering the mean-shift outlier model in both generalized linear and nonlinear regression models for the maximum likelihood framework. It should be noted that the mean-shift outlier model can be obtained by considering $\mbox{\boldmath$z$}=(0,\dots,1,\dots,0)^{\top}$ as the vector of zeros except for an 1 at the $i$th position, which leads to the standardized residual.

where $\widehat{U}_{i}=U_{i}(\widehat{\mbox{\boldmath$\beta$}})$, $\widehat{J}_{i}=J_{q}^{-\phi}(\theta_{i}(\widehat{\mbox{\boldmath$\beta$}}))$, $\widehat{g}_{i}=\dot{g}(\widehat{\theta}_{i}^{*})$, $\widehat{k}_{i}=\dot{k}(\widehat{\eta}_{i})$, and $\widehat{m}_{ii}$, $\widehat{m}_{ii}^{*}$ are, respectively, the diagonal elements of matrices $\widehat{\mbox{\boldmath$M$}}=\widehat{\mbox{\boldmath$P$}}\widehat{\mbox{\boldmath$W$}}\widehat{\mbox{\boldmath$J$}}$ and $\widehat{\mbox{\boldmath$M$}}{}^{2}$, for $i=1,\dots,n$. It is straightforward to note that $t_{i}$ has an approximately standard normal distribution and can be used for the identification of outliers and possibly influential observations. Evidently, we have that for $q=1$ we recover the standardized residual for generalized linear models. We have developed the standardized residual, $t_{i}$ based on the type-score statistic defined based on (3.4), further alternatives for standardized residuals can be constructed by using the Wald and bilinear-form statistics defined in Equations (3.3) and (3.5), respectively. Additionally, we can consider other types of residuals. Indeed, based on Equation (3.1) we can define the deviance residual as,

with $d(y_{i},\widehat{\mu}_{i})=2\{l_{q}(y_{i},\widehat{\mu}_{i})-l_{q}(y_{i},y_{i})\}$, the deviance associated with the $i$th observation. Moreover, it is also possible to consider the quantile residual (Dunn and Smyth, 1996), which is defined as

where $F(\cdot;\mu,\phi)$ is the cumulative distribution function associated with the random variable $y$, while $\Phi(\cdot)$ is the cumulative distribution function of the standard normal.

Throughout this work we have assumed $\phi$ as a fixed niusance parameter. Following Giménez et al. (2022a, b), we can obtain the maximum L$q$-likelihood estimator $\widehat{\phi}_{q}$ by solving the problem,

which corresponds to a profiled L$q$-likelihood function. It should be noted that to obtain $\widehat{\phi}_{q}$ it is sufficient to consider an one-dimensional optimization procedure.

We revisit the dataset introduced by Finney (1947) which aims to assess the effect of the volume and rate of inspired air on transient vasoconstriction in the skin of the digits. The response is binary, indicating the occurrence or non-occurrence of vasoconstriction. Following Pregibon (1981), a logistic regression model using the linear predictor $\eta=\beta_{1}+\beta_{2}\log(\text{volume})+\beta_{3}\log(\text{rate})$, was considered. Observations 4 and 18 have been identified as highly influential on the estimated coefficients obtained by the maximum likelihood method (Pregibon, 1981). In fact, several authors have highlighted the extreme characteristic of this dataset since the removal of these observations brings the model close to indeterminacy (see, for instance Künsch et al., 1989, Sec. 4). Table 4.1, reports the parameter estimates and their standard errors obtained from three robust procedures. Here CR denotes the fit using the method proposed by Cantoni and Ronchetti (2001), BY indicates the estimation method for logistic regression developed by Bianco and Yohai (1996) both available in the robustbase package (Todorov and Filzmoser, 2009) from the R software (R Development Core Team, 2022), while ML$q$, denotes the procedure proposed in this work where we have selected $q$ using the mechanism described in Section 2.2, in which $q=0.79$ was obtained. Furthermore, the results of the maximum likelihood fit are also reported and estimation was carried out using the weighted Bianco and Yohai estimator (Croux and Haesbroeck, 2003). However, since the explanatory variables do not contain any outliers the estimates obtained using BY and weighted BY methods are identical and therefore they are not reported here.

From Table 4.1 we must highlight the drastic change in the estimates of the regression coefficients by maximum likelihood when observations 4 and 18 are deleted (see also Pregibon, 1981; Wang, 1985). It is interesting to note that the results for this case are very close to those offered by the estimation procedure based on the $L_{1}$-norm introduced by Morgenthaler (1992). Moreover, as discussed in Sec. 6.17 of Atkinson and Riani (2000) the removal of observations 4, 18 and 24 leads to a perfect fit. Indeed, the elimination of observations 4 and 18 leads to a huge increment in the standard errors.

Figure 4.1 (a) shows a strong overlap in the groups with zero and unit response. The poor fit is also evident from the QQ-plot of Studentized residuals with simulated envelopes (see, Figure 4.1 (b)). The estimated weights $\widehat{U}_{i}=U_{i}(\widehat{\mbox{\boldmath$\beta$}})$, $(i=1,\dots,39)$ are displayed in Figure 4.2, where it is possible to appreciate the ability of the algorithm to downweight observations 4, 18 and 24. It should be stressed that the estimation of the regression coefficients by the maximum L$q$-likelihood method with $q=0.79$ is very similar to that obtained by the resistant procedure discussed by Pregibon (1982). Figure 4.3, reveals that the maximum L$q$-likelihood procedure offers a better fit and the QQ-plot with simulated envelope correctly identifies observations 4 and 18 as outliers.

In Appendix F the fit for the proposed model using values ranging from $q=1$ to 0.74, is reported. Table F.1 shows that as the value of $q$ decreases, there is an increase in the standard error of the estimates. Additionally, Figures F.1 and F.2 from Appendix F show how as $q$ decreases the weights of observations 4 and 18 decrease rapidly to zero. Indeed, when $q\leq 0.76$ the problem of indeterminacy underlying these data becomes evident.

To evaluate the performance of the proposed estimation procedure we developed a small Monte Carlo simulation study. We consider a Poisson regression model with logarithmic link function, based on the linear predictor

where $x_{ij}\sim\mathsf{U}(0,1)$ for $j=1,2,3$ and $\mbox{\boldmath$\beta$}=(\beta_{1},\beta_{2},\beta_{3})^{\top}=(1,1,1)^{\top}$. We construct $M=1000$ data sets with sample size $n=25,100$ and $400$ from the model (4.1). Following Cantoni and Ronchetti (2006) the observations are contaminated by multiplying by a factor of $\nu$ a percentage $\varepsilon$ of randomly selected responses. We have considered $\nu=2,5$ and $\varepsilon=0.00,0.05,0.10$ and $0.25$ as contamination percentages. For contaminated as well as non-contaminated data we have carried out the estimation by maximum likelihood, using the method of Cantoni and Ronchetti (2001) and the procedure proposed in this paper.