---
abstract: 'Generalized linear models are widely used in many areas of knowledge. As in other classes of regression models, it is desirable to perform diagnostic analysis in generalized linear models using residuals that are approximately standard normally distributed. Diagnostic analysis in this class of models are usually performed using the standardized Pearson residual or the standardized deviance residual. The former has skewed distribution and the latter has negative mean, specially when the variance of the response variable is high. In this work, we introduce the adjusted quantile residual for generalized linear models. Using Monte Carlo simulation techniques and two applications, we compare this residual with the standardized Pearson residual, the standardized deviance residual and two other residuals. Overall, the results suggest that the adjusted quantile residual is a better tool for diagnostic analysis in generalized linear models.'
address:
- 'Maths Science Institute and Computing, University of São Paulo, São Carlos, SP, Brazil.'
- ' Department of Statistics, Federal University of São Carlos, São Carlos, SP, Brazil.'
author:
- Juliana Scudilio
- 'Gustavo H. A. Pereira'
bibliography:
- 'bibliografia.bib'
title: '[ Adjusted quantile residual for generalized linear models ]{}'
---

Diagnostic analysis; generalized linear models; quantile residual.

Introduction
============

Generalized linear models (GLMs) were introduced by [@nelder72] and are widely used to study the relationship between a response variable and a set of predictor variables. It assumes that the distribution of the response variable given the set of predictors belongs to the exponential family and can be used to model many types of variables, such as binary variables, count variables and continuous skewed variables. The routines to fit GLMs are implemented in many statistical softwares and these models are described in details in [@mccullagh] and [@hardin2007generalized].

In regression models, it is important to check model adequacy and to identify outliers and influential observations. To reach these goals, it is interesting to use residuals whose distribution is well approximated by the standard normal distribution. When the residuals do not have this property, it is not unusual that they are also non-identically distributed (see, for example, Table 1 of [@anholeto2014]). In these cases, it is hard to interpret residuals plots and well-fitted models can be discarded (see Application 3 of [@pereira2017quantile]).

Many residuals are used to perform diagnostic analysis in GLMs, such as the deviance residual [@davison], the Pearson residual [@mccullagh], the Williams residual [@williams] and the Anscombe residual [@pierce1986residuals]. Usually they are standardized so that their variance become close to one. However, the distributions of these residuals are not approximated standard normal when the variance of the response variable is high (see Section \[simu\]). Some works introduced residuals for GLMs that are approximated standard normal distributed [@cordeiro2004pearson; @urbano2012wald], but they are complex and are not implemented in the statistical softwares.

The quantile residual [@dunn] is simple and general, but is not often used in GLMs. It is usually used to perform diagnostic analysis in complex regression models, such as the generalized additive models for location, scale and shape [@rigby2005]. The quantile residual is asymptotically standard normally distributed if the parameters are consistently estimated [@dunn], but its distribution is unknown in small sample sizes. [@feng2017randomized] compared the quantile residual with the deviance residual and the Pearson residual in generalized linear models. Using simulation studies, the authors concluded that the distribution of the quantile residual is better approximated by the standard normal distribution than that of the others residuals. In addition, [@feng2017randomized] showed that the quantile residual is the best for detecting lack of fit. However, the authors considered sample size equal to $1000$ and they did not note that the variance of the quantile residual is far from one in small sample sizes (see Section \[qu\]).

In this work, we introduce the adjusted quantile residual for GLMs. Using Monte Carlo simulation studies and two applications, we compare the behavior of the adjusted quantile residual and other residuals usually used in GLMs.

This paper is organized as follows. Section \[method\] presents the models used in this work, introduces the adjusted quantile residual and defines other residuals considered in this paper. In Section \[simu\], Monte Carlo simulation studies are performed to compare the properties of the defined residuals in several scenarios. In the following section, we study the ability of the residuals to identify model misspecification and outliers. Finally, in Section \[conc\], we present the conclusions.

Methodology {#method}
===========

Generalized linear models {#mod}
-------------------------

Let $y_1, y_2, \ldots, y_n$ be independent random variables, where each $y_i$ has a density functions in the linear exponential family defined as $$\pi(y;\theta_i,\phi)= \exp\left[\phi \left\{ y\theta_i-b(\theta_i)+c(y)\right\}+d(y;\phi)\right],
  \label{glm}$$

where $b(\cdot)$ and $c(\cdot)$ are know functions and $\phi_i$ is a precision parameter. It can be proved that $E(y_i)=\mu_i = b'(\theta_i)$ and $Var(y_i)=\phi^{-1}V(\mu_i)$, where $V_i=V(\mu_i)=d\mu_i/d\theta_i$ is the variance function.

Generalized linear models (GLMs) are defined by the family of distributions specified in equation \[glm\] and the systematic component $$g(\mu_i) = \eta_i = x_i^\top\boldsymbol{\beta},$$ where $x_{i}=(x_{i1},x_{i2},\ldots,x_{ik})^\top$ is a vector of known covariates ($k < n$), $\boldsymbol{\beta}=(\beta_1,\beta_2,\ldots,\beta_k)^\top$ is a vector of unknown parameters ($\boldsymbol{\beta} \in \mathbb{R}^k$) and $g(.)$ is a strictly monotonic and twice differentiable link function.

This paper focus on generalized linear models with continuous and asymmetric response variables. We considered only continuous variable as response variable, because in the discrete case the residual distribution cannot be continuous without randomization. We also use asymmetric distributions because in the normal case a residual with good properties is already known.

The gamma regression models and the inverse Gaussian regression models are the most important regression models in the class of GLMs with asymmetric and continuous response variable. They are considered in this paper and we define these models using the parameterization used by [@stasinopoulos08].

The gamma distribution with parameters $\mu_i$ and $\sigma$ has a density function defined as

$$\begin{aligned}
\label{gama}
f_1(y;\mu_i,\sigma)=\frac{exp \left(-\frac{ y}{\sigma^2\mu_i}\right)}{\Gamma(\frac{1}{\sigma^2})} \frac{y^{(\frac{1}{\sigma^2}-1)}}{\left(\mu_i\sigma^2\right)^{\frac{1}{\sigma^2}}},\end{aligned}$$

where $y>0$, $\mu_i>0$ and $\sigma>0$. In this parameterization of the gamma distribution, $\mbox{E}(y_i)=\mu_i$ and $\mbox{Var}(y_i)=\sigma^2\mu_i^2$.

The inverse Gaussian distribution with parameters $\mu_i$ and $\sigma$ is defined as $$f_2(y;\mu_i,\sigma)=\frac{1}{\sqrt{(2\pi y^3 \sigma^2)}}exp\left\{-\frac{(y-\mu_i)^2}{2\mu_i^2\sigma^2y}\right\},
 \label{invga}$$ where $y>0$, $\mu_i>0$ and $\sigma>0$. In this parameterization of the inverse Gaussian, $\mbox{E}(y_i)=\mu_i$ and $\mbox{Var}(y_i)=\sigma^2\mu_i^3$.

The parameters of the GLMs can be estimated by maximum likelihood using numerical optimization algorithms and hypothesis tests can be conducted using Wald, score or likelihood ratio statistic. The *glm* function and the *gamlss* package, which are available for the R statistical software, can be used for fitting gamma and inverse Gaussian regression models.

Standardized quantile residual {#qu}
------------------------------

The quantile residual can be used in a wide class of regression models. When the response variable is discrete, it uses uniform distributed random variables and it is called randomized quantile residual. For the subclass of GLMs considered in this paper, it is defined as

$$r_i^{qu}=\Phi^{-1}\{F(y_i;\hat{\mu}_i,\hat{\sigma})\},
 \label{quanti}$$

where $\Phi(.)$ is the cumulative distribution function of the standard normal distribution, $F(.)$ is the cumulative distribution function of the gamma distribution or inverse Gaussian distribution and $\hat{\mu}_i$ and $\hat{\sigma}$ are the maximum likelihood estimates of the parameters $\mu_i$ and $\sigma$, respectively.

According to [@dunn], [$r_i^{qu}$]{} is asymptotically standard normally distributed, but it is important to study its properties when sample size is not large. We used Monte Carlo simulation studies to study the properties of the quantile residual in small sample sizes. Gamma and inverse Gaussian regression models were considered and we assumed the model $log(\mu_i)=\beta_0+\beta_{1}x_{i1}+\beta_{2}x_{i2}$. The covariate values were generated as independent draws from the standard uniform distribution and remained constant throughout the simulations. We considered $\beta_{0}=3$, $\beta_1=2$ and $\beta_2=1$, which resulted in $\mu \in (20{.}085 ; 403{.}429)$ and used $\sigma=0{.}1$ for gamma regression model and $\sigma=0{.}02$ for inverse Gaussian regression model.

The left side of Tables \[tab-quant\] and \[tab-quant1\] present the sample mean, variance, skewness and kurtosis coefficients for [$r_i^{qu}$]{} for gamma and inverse Gaussian regression models, respectively. The tables also contain the value of the Anderson-Darling (AD) statistic [@anderson], used to test whether each residual is standard normally distributed. The value of the Anderson-Darling statistic is used as a closeness measure between each residual distribution and standard normal distribution. Results are based on $5000$ Monte Carlo replications. Note that for both models the variance of [$r_i^{qu}$]{} is far from 1.

=0.07cm

\[tab-quant\]

In GLMs, Pearson and deviance residuals are standardized so that their variance becomes close to one. Our proposal is to use the same standardization term used for these residual in the quantile residual, ie, divide the residual [$r_i^{qu}$]{} by $\sqrt{(1-\hat{h}_{ii})}$, where $\hat{h}_{ii}$ is the $i$-th diagonal element of matrix $\hat{\boldmath{H}}=\hat{W}^{1/2}X(X^T\hat{W}X)^{-1}X^T\hat{W}^{1/2}$ and $\hat{W}=diag\{\hat{w_1},....,\hat{w_n}\}$ is the diagonal matrix of weights. In particular, for the gamma distribution and for the inverse Gaussian distribution, $w_i=({d\mu_i}/{d\eta_i})^2/\mu_i^2$ and $w_i=({d\mu_i}/{d\eta_i})^2/\mu_i^3$, respectively. Then, the adjusted quantile residual is defined as $$r_i^{*qu} = \frac{r_i^{qu}}{\sqrt{(1-\hat{h}_{ii})}}.$$ We used the term adjusted instead of standardized because [@klar2012specification] introduced a different residual called standardized quantile residual to derive a goodness-of-fit test for generalized linear models. We did not consider the standardized quantile residual in this work because it assumes that the variance of the quantile residual is constant across observations. Tables \[tab-quant\] and \[tab-quant1\] suggest that this assumption is not reasonable.

=0.07cm

\[tab-quant1\]

The results of the adjusted quantile residual are shown on the right side of Tables \[tab-quant\] and \[tab-quant1\]. Note that dividing [$r_i^{qu}$]{} by $\sqrt{(1-\hat{h}_{ii})}$ results in a residual with mean, skewness and kurtosis similar to [$r_i^{qu}$]{}, a variance close to $1$ and smaller AD than the [$r_i^{qu}$]{}. In other words, the distribution of [$r_i^{*qu}$]{} is closer to the standard normal distribution than that of the [$r_i^{qu}$]{}. We did the same analysis in several scenarios, which were omitted in this paper because they resulted on similar conclusions.

Other residuals in GLMs {#others}
-----------------------

The other residuals considered in this paper are the deviance residual, the Pearson residual, the Anscombe residual and the William residual. The standardized deviance residual is defined as $$r_i^{dev}=\frac{d(y_i;\hat{\mu}_i)}{\hat{\sigma}\sqrt{1-\hat{h}_{ii}}},$$ where $d(y_i;\hat{\mu}_i)=\left\{-log(\frac{y_i}{\hat{\mu_i}})+\frac{(y_i-\hat{\mu}_i)^2}{y_i\hat{\mu}_i^2}\right\}$ for the gamma distribution and $d(y_i;\hat{\mu}_i)=\left\{\frac{(y_i-\hat{\mu}_i)}{y_i\hat{\mu}_i^2}\right\}$ for the inverse Gaussian distribution.

The standardized Pearson residual is given by

$$r_i^{pea}=\frac{(y_i-\hat{\mu}_i)}{\hat{\sigma}\sqrt{\hat{\mu}_i^2(1-\hat{h}_{ii})}}.
    \label{pea}$$

The William residual is also used in GLMs and is defined as $$r_{i}^{wil}=sign\{y_i-\hat{\mu}_i\}[(1-\hat{h}_{ii})(r_i^{dev})^2+\hat{h}_{ii}(r_i^{pea})^2]^{1/2}.$$

The last residual used in this paper is the standardized Anscombe residual given by $$r_i^{ans}=\frac{\psi(y_i)-\psi(\hat{\mu}_i)}{\hat{\sigma}\hat{V}^{1/2}(\hat{\mu})\psi'(\hat{\mu}_i)\sqrt{(1-\hat{h}_{ii})}},$$ where $\psi(\mu_i)=3\mu_i^{1/3}$ and $V(\mu_i)=\mu^2$ for the gamma distribution and $\psi(\mu_i)=log(\mu_i)$ and $V(\mu_i)=\mu^3$ for the inverse Gaussian distribution.

Simulation studies {#simu}
==================

In this section, we use Monte Carlo simulation studies to compare the properties of the residuals in the GLMs. Here, we use the gamma and the inverse Gaussian regression models. We considered seven scenarios for each distribution to study the residuals’ properties. In the first six scenarios, we assumed the model as $$log(\mu_i)=\beta_{0}+\beta_{1}x_{i1}+\beta_{2}x_{i2}.$$

In the first four scenarios and in the last scenario for gamma and inverse Gaussian regression models, the covariates values were generated as independent draws from the standard uniform distribution. For the others scenarios, $x_{i1}$ values were generated from the normal distribution, ($x_{i1}\sim$ NO $(0.5,0.25^2)$), $x_{i2}$ values were generated from the inverse Gaussian distribution in scenarios V-a and V-b ($x_{i2}\sim $ IG $ (0.4,2)$) and from the gamma distribution in scenarios VI-a and VI-b ($x_{i2}\sim $ GA $ (0.4,1)$). The covariates values remained constant throughout the simulations.

In order to check if the results change considerably when other link function is used, in the last scenario we used the canonical link function. In Scenario VII-a, we used the gamma regression model with the inverse link function, in which $$\begin{aligned}
\frac{1}{\mu_i} = \beta_{0}+\beta_{1}x_{i1}+\beta_{2}x_{i2}.\end{aligned}$$

In the last scenario (Scenario VII-b), we used the inverse Gaussian regression model with the link function $1/\mu^2$, in which $$\begin{aligned}
    \frac{1}{\mu_i^2} = \beta_{0}+\beta_{1}x_{i1}+\beta_{2}x_{i2}.\end{aligned}$$

Tables \[tab:1\] and \[tab:2\] present the description of the scenarios for gamma and inverse Gaussian regression models. In the first scenario (Scenario I-a e I-b), we considered $\beta_{0}=3$, $\beta_1=2$ and $\beta_2=1$, which resulted in $\mu \in (20{.}085 ; 403{.}429)$ and used $\sigma=0{.}1$ for gamma regression model and $\sigma=0{.}02$ for inverse Gaussian regression model. In the second and third scenarios, we used the same $\beta_j$, $j=0,1,2$, of the Scenario I, but with different $\sigma$ values. In scenarios II-a and II-b we decreased the variance and in scenarios III-a and III-b we increased the variance. In the following scenarios (scenarios IV-a and IV-b), we changed the values of the coefficients to result in mean response values close to zero. In scenarios V and VI, we changed the generating distribution of $x_{i1}$ and $x_{i2}$. In the last scenarios (scenarios VII-a and VII-b), we changed the link function. For all scenarios, two different sample sizes were considered: $n=15$ and $n=50$. All results are based on $5000$ Monte Carlo replications. All simulations were performed using the software R. Results for [$r_i^{qu}$]{} were omitted from the tables, because the simulation studies showed that their distributions are worse approximated by standard normal distribution than that of [$r_i^{*qu}$]{} in all scenarios.

   Scenario                                                                                                                                                                  
  ---------- ------------------- ------------------------------------------------------------ --------------------------------------------------- -------------------------- ----------
     I-a        $log(\mu_i)$                      $\beta_0 = 3$, $\beta_1 = 2$, $\beta_2 = 1$              $x_{i1}\sim U(0,1); x_{i2}\sim U(0,1)$   $(20{.}085 ; 403{.}429)$   $0{.}10$
     II-a       $log(\mu_i)$                      $\beta_0 = 3$, $\beta_1 = 2$, $\beta_2 = 1$              $x_{i1}\sim U(0,1); x_{i2}\sim U(0,1)$   $(20{.}085 ; 403{.}429)$   $0{.}05$
    III-a       $log(\mu_i)$                      $\beta_0 = 3$, $\beta_1 = 2$, $\beta_2 = 1$              $x_{i1}\sim U(0,1); x_{i2}\sim U(0,1)$   $(20{.}085 ; 403{.}429)$   $0{.}50$
     IV-a       $log(\mu_i)$                 $\beta_0 = -3$, $\beta_1 = 1{.}5$, $\beta_2 = 1$              $x_{i1}\sim U(0,1); x_{i2}\sim U(0,1)$       $(0{.}049; 0{.}606)$   $0{.}10$
     V-a        $log(\mu_i)$                      $\beta_0 = 3$, $\beta_1 = 2$, $\beta_2 = 1$   $x_{i1}\sim NO(0.5,0.25^2); x_{i2}\sim IG(0.4,2)$   $(20{.}085 ; 403{.}429)$   $0{.}10$
     VI-a       $log(\mu_i)$                      $\beta_0 = 3$, $\beta_1 = 2$, $\beta_2 = 1$   $x_{i1}\sim NO(0.5,0.25^2); x_{i2}\sim GA(0.4,1)$   $(7{.}400 ; 1892{.}000)$   $0{.}10$
    VII-a     $\frac{1}{\mu_i}$    $\beta_0 =0{.}0025$ $\beta_1 = 0{.}04$, $\beta_2 = 0{.}01$              $x_{i1}\sim U(0,1); x_{i2}\sim U(0,1)$         $(19{.}04 ; 4000)$   $0{.}10$

  : Description of the scenarios for gamma regression model.

\[tab:1\]

   Scenario                                                                                                                                                                    
  ---------- --------------------- ------------------------------------------------------------ --------------------------------------------------- -------------------------- ----------
     I-b         $log(\mu_i)$                       $\beta_0 = 3$, $\beta_1 = 2$, $\beta_2 = 1$              $x_{i1}\sim U(0,1); x_{i2}\sim U(0,1)$   $(20{.}085 ; 403{.}429)$   $0{.}02$
     II-b        $log(\mu_i)$                       $\beta_0 = 3$, $\beta_1 = 2$, $\beta_2 = 1$              $x_{i1}\sim U(0,1); x_{i2}\sim U(0,1)$   $(20{.}085 ; 403{.}429)$   $0{.}01$
    III-b        $log(\mu_i)$                       $\beta_0 = 3$, $\beta_1 = 2$, $\beta_2 = 1$              $x_{i1}\sim U(0,1); x_{i2}\sim U(0,1)$   $(20{.}085 ; 403{.}429)$   $0{.}03$
     IV-b        $log(\mu_i)$                  $\beta_0 = -3$, $\beta_1 = 1{.}5$, $\beta_2 = 1$              $x_{i1}\sim U(0,1); x_{i2}\sim U(0,1)$       $(0{.}049; 0{.}606)$   $0{.}50$
     V-b         $log(\mu_i)$                       $\beta_0 = 3$, $\beta_1 = 2$, $\beta_2 = 1$   $x_{i1}\sim NO(0.5,0.25^2); x_{i2}\sim IG(0.4,2)$   $(20{.}085 ; 403{.}429)$   $0{.}02$
     VI-b        $log(\mu_i)$                       $\beta_0 = 3$, $\beta_1 = 2$, $\beta_2 = 1$   $x_{i1}\sim NO(0.5,0.25^2); x_{i2}\sim GA(0.4,1)$   $(7{.}400 ; 1892{.}000)$   $0{.}02$
    VII-b     $\frac{1}{\mu_i^2}$    $\beta_{0}=0{.}000006$, $\beta_1=0{.}002$, $\beta_2=0.001$              $x_{i1}\sim U(0,1); x_{i2}\sim U(0,1)$     $(18{.}24 ; 408{.}25)$   $0{.}02$

  : Description of the scenarios for inverse Gaussian regression model.

\[tab:2\]

Tables \[cen1ga\] and \[tab:cen1\] present the simulation results considering $n=15$ in the Scenario I-a for gamma regression model and in the Scenario I-b for inverse Gaussian regression model, respectively. The means of [$r_i^{*qu}$]{} and [$r_i^{pea}$]{} are very close to zero, but the other residuals have means slightly lower than zero. For all residuals, the variance is close to one, but none of the residuals have a variance close to one for all $15$ observations. The skewness coefficient is considerably greater than zero for [$r_i^{pea}$]{} specially for inverse Gaussian regression model, very close to zero for [$r_i^{wil}$]{} and slightly lower than zero for the other residuals. The excess kurtosis coefficients are similar and not too close to zero for the five residuals. The residual [$r_i^{*qu}$]{} has smaller AD statistic than the other residuals in 10 out of 15 observations for gamma regression model and 14 out of 15 observations for inverse Gaussian regression model. The AD statistic is smaller for the adjusted quantile residual, than that for the other residuals, because the mean of [$r_i^{*qu}$]{} is closer to zero than that of [$r_i^{dev}$]{}, [$r_i^{wil}$]{} and [$r_i^{ans}$]{} and the skewness is closer to zero than that of [$r_i^{pea}$]{}.

=0.07cm

\[cen1ga\]

=0.07cm

\[tab:cen1\]

Table \[tab:gama\] presents the average values for the mean, variance, skewness and excess kurtosis for all residuals and both sample sizes for gamma regression model. Changes in the range of $\mu$, in the distribution of the covariates and in the link function does not seem to affect considerably the distribution of the residuals, since the results of Scenarios IV-a to VII-a are very similar to that of Scenario I-a. However, the distributions of the five residuals are worse approximated by the standard normal distribution when variance increases. The skewness coefficients are far from zero in Scenario III-a than in Scenario I-a and, except for [$r_i^{pea}$]{}, the mean and variance are also far from zero and one, respectively. As a consequence, the AD statistics increase substantially from Scenario I-a to Scenario III-a, but the rise is smaller for [$r_i^{*qu}$]{}.

In general, when sample size increases from $n=15$ to $n=50$, all residuals’s distributions are better approximated by the standard normal distribution. In all scenarios, the absolute value of the average skewness coefficient reduces considerably, except for the residual [$r_i^{pea}$]{}, in which the average skewness coefficient rises slightly. Additionally, the absolute average excess kurtosis coefficient also decreases in all scenarios for all residuals, except for [$r_i^{pea}$]{} in Scenario III-a. Finally, in six out of seven scenarios, the AD statistic reduces for all residuals, but the reduction is greater for [$r_i^{*qu}$]{}.

Table \[adga\] summarizes the results of the AD statistic in all scenarios for $n=15$ and $n=50$ considering the gamma regression model. When $n=50$ the mean, the standard deviation and the maximum of the AD statistic are lower for [$r_i^{*qu}$]{} than for other residuals. When sample size is 15, the sample mean of the AD statistic is also smaller for [$r_i^{*qu}$]{} than for other residuals, except in Scenario IV-a in which [$r_i^{pea}$]{} and [$r_i^{*qu}$]{} have similar sample mean of the AD statistic. On the other hand, the standard deviation and the maximum are the lowest for [$r_i^{*qu}$]{} in some scenarios and for [$r_i^{pea}$]{} in others.

Tables \[tab:ig\] and \[adig\] present the average of the sample moments for the residuals and a descriptive summary of the results of the AD statistic for the inverse Gaussian regression model. For all residuals, in the most of the scenarios, the average values of the AD statistic are higher for the inverse Gaussian regression model than for the gamma regression model. A plausible explanation for this result is that, for the same mean and variance, the inverse Gaussian distribution presents a greater skewness than the gamma distribution. In general, the standard deviation of AD statistic are also higher for the inverse Gaussian regression model. Residual [$r_i^{*qu}$]{} presents lower mean of the AD statistic than the other residuals for all scenarios and both sample sizes. Additionally, the differences between the mean AD statistic for [$r_i^{*qu}$]{} and for the other residuals are greater for the inverse Gaussian regression model. Finally, except for Scenario V-a when $n=15$, the standard deviation of the AD statistic is lower for [$r_i^{*qu}$]{} than for the other residuals.

=0.1cm

\[tab:gama\]

=0.2cm

=0.1cm

\[tab:ig\]

Application {#apps}
===========

In Section \[simu\], we studied whether the distribution of the residuals considered in this paper are well approximated by the standard normal distribution. Other essential properties for residuals are their ability to identify model misspecification and outliers. We used two real datasets to study the residuals considered in this paper according to these aspects. In this section, we compare the [$r_i^{*qu}$]{} with the two most commonly used residual in the GLMs, the [$r_i^{dev}$]{} and the [$r_i^{pea}$]{}. Applications were performed using the *gamlss* package and the *glm* function of the R software.

Oil dataset
-----------

The first application uses a database with 1000 observations about the daily prices of WTI (West Texas Intermediate) oil price traded by NYMEX(New York Mercantile Exchange). We used this dataset to investigate whether the [$r_i^{dev}$]{}, [$r_i^{pea}$]{}and [$r_i^{*qu}$]{} are able to identify model misspecification. The response variable is the price of oil. Here, we use two explanatory variables: the lagged version of the response variable and the log price of front month heating oil contract traded by NYMEX, as proposed by [@rigby2015]. We fitted inverse Gaussian regression model with logit link function for the oil data and obtained the three sets of residuals. [@rigby2015] proposed the Sinh-Arcsinh (SHASH) distribution and splines to obtain a reasonable fit for this response variable. Therefore, lack of fit is expected when we fit the inverse Gaussian regression model. Figure \[segreg\] presents a plot of residuals versus linear predictor and a normal residual plot with simulated envelope [@atkinson] for [$r_i^{dev}$]{}, [$r_i^{pea}$]{}and [$r_i^{*qu}$]{}. Clearly, the plots for the three residuals suggest lack of fit for the inverse Gaussian regression model, indicating that these residuals are able to identify this kind of problem. A similar analysis was performed using the gamma regression model. The conclusions were the same and the plots were omitted for the sake of brevity.

![image](dev_oil1.pdf){width="0.35\linewidth"} ![image](env_dev_oil.pdf){width="0.35\linewidth"}\

![Residuals against linear predictor (left) and normal residual plots with simulated envelope (right) for oil dataset.[]{data-label="segreg"}](pea_oil.pdf "fig:"){width="0.35\linewidth"} ![Residuals against linear predictor (left) and normal residual plots with simulated envelope (right) for oil dataset.[]{data-label="segreg"}](env_pea_oil.pdf "fig:"){width="0.35\linewidth"}\
![Residuals against linear predictor (left) and normal residual plots with simulated envelope (right) for oil dataset.[]{data-label="segreg"}](qu_oil.pdf "fig:"){width="0.35\linewidth"} ![Residuals against linear predictor (left) and normal residual plots with simulated envelope (right) for oil dataset.[]{data-label="segreg"}](env_qu_oil.pdf "fig:"){width="0.35\linewidth"}\

Turbine dataset
---------------

The second real dataset consists of the time to evaluate the performance of five types of high-speed turbines for plane engines [@lawless]. The response variable is the time (in million of cycles units) until the loss of velocity and the covariate is the type of turbine. Our focus with this application is to check if [$r_i^{dev}$]{}, [$r_i^{pea}$]{}, [$r_i^{*qu}$]{} are good residuals to perform outlier identification. We fitted gamma regression model with identity link function to this dataset and obtained the three set of residuals. Figure \[turbine\] presents a plot of residuals against linear predictors and a normal residual plot with simulated envelope for [$r_i^{dev}$]{}, [$r_i^{pea}$]{}and [$r_i^{*qu}$]{}. None of the residuals suggest model misspecification.

![Residuals against linear predictor (left) and half-normal residual plots with simulated envelope (right) for turbine data.[]{data-label="turbine"}](dev-turbine.pdf "fig:"){width="0.35\linewidth"} ![Residuals against linear predictor (left) and half-normal residual plots with simulated envelope (right) for turbine data.[]{data-label="turbine"}](env-dev-turbine.pdf "fig:"){width="0.35\linewidth"}\
![Residuals against linear predictor (left) and half-normal residual plots with simulated envelope (right) for turbine data.[]{data-label="turbine"}](pea-turbine.pdf "fig:"){width="0.35\linewidth"} ![Residuals against linear predictor (left) and half-normal residual plots with simulated envelope (right) for turbine data.[]{data-label="turbine"}](env-pea-turbine.pdf "fig:"){width="0.35\linewidth"}\
![Residuals against linear predictor (left) and half-normal residual plots with simulated envelope (right) for turbine data.[]{data-label="turbine"}](qu-turbine.pdf "fig:"){width="0.35\linewidth"} ![Residuals against linear predictor (left) and half-normal residual plots with simulated envelope (right) for turbine data.[]{data-label="turbine"}](env-qu-turbine.pdf "fig:"){width="0.35\linewidth"}

For the three residuals, observation $40$ has the highest value and observation $1$ has the smallest. However, the value of the residual for observation $40$ is considerably greater for [$r_i^{pea}$]{} than for the other residuals and the absolute value of the residual for observation $1$ is substantially smaller for [$r_i^{pea}$]{} than for [$r_i^{dev}$]{} and [$r_i^{*qu}$]{}. In the simulation studies, we noted that the distribution of [$r_i^{pea}$]{} is right skewed and this explains the mentioned result. As a consequence, [$r_i^{pea}$]{} may incorrectly identify as an outlier an observation with positive value for the residual and may not identify an outlier with value of [$r_i^{pea}$]{} negative. On the other hand, the values of [$r_i^{dev}$]{} and [$r_i^{*qu}$]{} are close for all observations. However, the values of the positive residual are slightly greater for [$r_i^{*qu}$]{} than for [$r_i^{dev}$]{} and the absolute values of the negative residuals are slightly smaller for [$r_i^{*qu}$]{} than for [$r_i^{dev}$]{}. As [$r_i^{dev}$]{} has mean far from zero than [$r_i^{*qu}$]{} and similar variance, skewness and excess kurtosis, the latter seems to be a better measure of the discrepancy of an observation and, consequently, better to perform outlier identification.

Conclusions {#conc}
===========

In this work, we introduced the adjusted quantile residual to perform diagnostic analysis in generalized linear models. We compared the distribution of the adjusted quantile residual with four other residuals using Monte Carlos simulation studies. Additionally, we used two applications to investigate if the residuals are able to identify model misspecification and outliers.

It is very desirable to find a residual whose distribution is well approximated by the standard normal distribution. Our simulation studies suggest that the adjusted quantile residual’s distribution better approximates to the standard normal distribution in all scenarios than the other residual’s distribution, specially when sample size is not small. Simulation results also suggest that, when variance increases, all residuals worsen regarding normal approximation. However, when variance is high, the distribution of the quantile residual is much better approximated by the standard normal distribution than that of the other residuals, even in small samples.

The applications investigated other properties of the residuals. The first application suggested that the three residuals considered in the analysis can detect lack of fit in generalized linear models. Based on the second application, the adjusted quantile residual seems to be the best to perform outlier identification.

The standardized Pearson residual and the standardized deviance residual are calculated by many packages and statistical softwares in their generalized linear models routine, but commonly the quantile residual is not implemented. However, adjusted quantile residual is simple and easy to calculate using any statistical package. Considering the results of the simulations studies, the applications and its simplicity, the adjusted quantile residual is a better choice of to perform diagnostic analysis in generalized linear models than the competing residuals.

Acknowledgements
================

The authors thank “Coordenação de Aperfeiçoamento de Pessoal de Nível Superior” (CAPES) for the financial support received for this project.

References
==========
