A residual for outlier identification in zero adjusted regression models

TL;DR

This paper introduces a new residual for zero adjusted regression models that improves outlier detection beyond the existing randomized quantile residual, enhancing diagnostic accuracy.

Contribution

A novel residual specifically designed for outlier detection in zero adjusted regression models, outperforming existing methods in identifying anomalies.

Findings

The new residual effectively detects outliers missed by the randomized quantile residual.

Simulation studies show the residual has good statistical properties.

Application results confirm improved outlier identification.

Abstract

Zero adjusted regression models are used to fit variables that are discrete at zero and continuous at some interval of the positive real numbers. Diagnostic analysis in these models is usually performed using the randomized quantile residual, which is useful for checking the overall adequacy of a zero adjusted regression model. However, it may fail to identify some outliers. In this work, we introduce a residual for outlier identification in zero adjusted regression models. Monte Carlo simulation studies and an application suggest that the residual introduced here has good properties and detects outliers that are not identified by the randomized quantile residual.

Tables

Table 1. Table 1: Descriptive statistics for the percentage of residuals in each interval - ZABE regression model - Scenario 1

	Residual	Theoretical	Simulation results
Residual	interval	Value	Min	Q1	Median	Mean	Q3	Max
$r_i^{\star q}$		0.13	0.04	0.10	0.11	0.11	0.12	0.20
		2.28	1.84	2.15	2.23	2.24	2.31	2.53
		15.87	15.41	15.90	16.05	16.04	16.16	16.73
	$>1$	15.87	15.19	15.82	15.99	16.00	16.21	16.70
	$>2$	2.28	1.89	2.21	2.27	2.28	2.37	2.56
	$>3$	0.13	0.06	0.10	0.12	0.11	0.13	0.19
$r_i^{\star asb}$		0.13	0.13	0.20	0.24	0.24	0.27	0.37
		2.28	2.30	2.55	2.62	2.63	2.73	2.96
		15.87	14.80	15.33	15.59	15.58	15.80	16.35
	$>1$	15.87	15.33	16.03	16.28	16.27	16.44	17.09
	$>2$	2.28	1.48	1.75	1.81	1.84	1.93	2.26
	$>3$	0.13	0.00	0.03	0.04	0.04	0.06	0.12
$r_i^{\star efc1}$		0.13	0.12	0.20	0.25	0.24	0.28	0.37
		2.28	2.24	2.54	2.64	2.63	2.75	3.03
		15.87	14.77	15.34	15.59	15.58	15.81	16.34
	$>1$	15.87	15.41	16.10	16.28	16.26	16.44	17.06
	$>2$	2.28	1.39	1.73	1.84	1.83	1.96	2.22
	$>3$	0.13	0.00	0.02	0.04	0.04	0.06	0.10
$r_i^{\star efc2}$		0.13	0.16	0.23	0.26	0.27	0.30	0.40
		2.28	2.43	2.67	2.76	2.77	2.87	3.14
		15.87	14.96	15.56	15.79	15.79	16.04	16.64
	$>1$	15.87	15.59	16.30	16.50	16.50	16.65	17.27
	$>2$	2.28	1.56	1.86	1.97	1.97	2.09	2.39
	$>3$	0.13	0.00	0.04	0.05	0.05	0.07	0.12
$r_i^{of}$		0.13	0.36	0.56	0.75	0.74	0.91	1.31
		2.28	3.19	3.63	3.92	3.87	4.10	4.51
		15.87	11.86	13.24	13.90	13.87	14.61	15.57
	$>1$	15.87	12.29	13.73	14.39	14.32	14.92	15.74
	$>2$	2.28	2.68	2.99	3.09	3.09	3.19	3.46
	$>3$	0.13	0.17	0.24	0.29	0.30	0.34	0.47

Table 2. Table 2: Descriptive statistics for the percentage of residuals in each interval - ZAIG regression model - Scenario 1

	Residual	Theoretical	Simulation results
Residual	interval	Value	Min	Q1	Median	Mean	Q3	Max
		0.13	0.06	0.10	0.12	0.12	0.13	0.18
		2.28	1.54	2.22	2.33	2.30	2.41	2.64
$r_i^{\star q}$		15.87	14.99	15.76	16.00	15.94	16.17	16.62
	$>1$	15.87	14.46	15.84	16.03	16.01	16.20	16.86
	$>2$	2.28	1.58	2.16	2.26	2.23	2.36	2.56
	$>3$	0.13	0.03	0.08	0.10	0.09	0.11	0.16
		0.13	0.11	0.16	0.18	0.18	0.21	0.28
		2.28	2.24	2.78	2.92	2.92	3.09	3.54
$r_i^{\ast dev}$		15.87	16.09	17.28	17.79	17.87	18.51	20.33
	$>1$	15.87	12.21	13.98	14.54	14.34	14.76	15.71
	$>2$	2.28	1.35	1.78	1.88	1.87	1.97	2.14
	$>3$	0.13	0.03	0.06	0.07	0.07	0.08	0.12
		0.13	0.00	0.00	0.00	0.01	0.01	0.04
		2.28	0.09	0.92	1.40	1.24	1.57	1.91
$r_i^{\ast pea}$		15.87	15.98	16.84	17.09	17.04	17.32	17.83
	$>1$	15.87	14.08	14.71	14.95	14.93	15.14	15.98
	$>2$	2.28	2.20	2.90	3.11	3.10	3.30	3.92
	$>3$	0.13	0.11	0.34	0.44	0.47	0.58	0.96
		0.13	0.06	0.12	0.13	0.13	0.15	0.20
		2.28	1.74	2.60	2.75	2.72	2.88	3.16
$r_i^{\ast ans}$		15.87	15.80	17.12	17.69	17.70	18.39	20.06
	$>1$	15.87	11.94	13.79	14.36	14.20	14.63	15.54
	$>2$	2.28	1.18	1.60	1.74	1.72	1.85	2.02
	$>3$	0.13	0.01	0.04	0.06	0.06	0.07	0.10
		0.13	0.10	0.15	0.17	0.17	0.20	0.26
		2.28	2.18	2.72	2.89	2.88	3.06	3.46
$r_i^{\ast wil}$		15.87	16.09	17.27	17.78	17.85	18.50	20.27
	$>1$	15.87	12.27	14.02	14.55	14.36	14.77	15.72
	$>2$	2.28	1.39	1.82	1.91	1.91	2.01	2.19
	$>3$	0.13	0.04	0.07	0.08	0.08	0.09	0.13
		0.13	0.44	0.60	0.72	0.71	0.81	0.96
		2.28	3.25	3.66	3.81	3.78	3.93	4.20
$r_i^{\ast \ast q}$		15.87	9.86	12.09	12.84	12.69	13.68	14.76
	$>1$	15.87	9.66	12.16	12.92	12.73	13.42	14.57
	$>2$	2.28	3.10	3.55	3.74	3.72	3.87	4.21
	$>3$	0.13	0.32	0.54	0.65	0.65	0.75	0.96

Table 3. Table 3: The final ZABE regression model for the writing essay score

Equation	Variable	Estimate	Standard	P-value
			error
$\mu$	Intercept	-3.3924	0.1511
	Natural sciences score	0.0015	0.0004
	Human sciences score	0.0019	0.0004
	Language score	0.0025	0.0004	0.0011
	Mathematics score	0.0009	0.0003
	Male student	-0.1968	0.0351	0.0004
	Age greather than 25	-0.1417	0.0396
$\varphi$	Intercept	-1.4090	0.1751
	Human sciences score	-0.0009	0.0003	0.0084
	Language score	0.0016	0.0004
$\alpha$	Intercept	6.0367	1.1238
	Human sciences score	-0.0094	0.0025	0.0002
	Language score	-0.0094	0.0028	0.0006
	Age greater than 25	0.7444	0.3166	0.0189

Equations

$$f_Y(y;\alpha ,\mu ,\varphi )=\{\begin{array}{ll}\alpha & \text{if}y\in \{0\},\\ (1-\alpha )f_W(y;\mu ,\varphi )& \text{if}y\in \mathcal{Y},\end{array}$$

$$g_j({\theta }_{ij})=\{\begin{array}{l}{g}_1({\mu }_i)={\eta }_{i1}\\ g_2({\varphi }_i)={\eta }_{i2}\\ g_3({\alpha }_i)={\eta }_{i3}\end{array}\}={\eta }_{ij},i=1,2,\dots ,n;j=1,2,3,$$

$$r_{i}^{q}=\left\{\begin{array}[]{ccc}\Phi^{-1}(u_{i}),&\text{if}&y_{i}=0,\\ \Phi^{-1}(F(y_{i};\hat{\theta}_{i})),&\text{if}&y_{i}>0,\end{array}\right.$$

$$r_i^{\star }=\{\begin{array}{ll}{\Phi }^{-1}[(1-\Phi (|r_i|))(1-{\hat{\alpha }}_i)],& \text{if ri<0},\\ {\Phi }^{-1}[1-(1-\Phi (|r_i|))(1-{\hat{\alpha }}_i)],& \text{if ri>0}.\end{array}$$

$$r_i^{\star }=\{\begin{array}{ll}{\Phi }^{-1}[\Phi (r_i)(1-{\hat{\alpha }}_i)],& \text{if ri<0},\\ {\Phi }^{-1}[{\hat{\alpha }}_i+\Phi (r_i)(1-{\hat{\alpha }}_i)],& \text{if ri>0}.\end{array}$$

$$\text{Pr}(r_i^{\star }<-k)=\text{Pr}(r_i^{\star }>k)=1-\Phi (k).$$

$$\{\begin{array}{ll}\log \left(\frac{{\mu }_i}{1-{\mu }_i}\right)& ={\beta }_{11}+{\beta }_{12}{x}_{i12}+{\beta }_{13}{x}_{i13},\\ \log ({\varphi }_i)& ={\beta }_{21},\\ \log \left(\frac{{\alpha }_i}{1-{\alpha }_i}\right)& ={\beta }_{31}+{\beta }_{32}{x}_{i32}+{\beta }_{33}{x}_{i33},\end{array}$$

$$r_i^{\star q}$$

$$\text{Pr}(r_i^{\ast }<-k)$$

Full text

1Department of Statistics, Federal University of São Carlos

2Department of Applied Mathematics and Statistics, University of São Paulo

3Department of Statistics, Federal University of Campina Grande

4Department of Statistics, University of São Paulo

Key words: diagnostic analysis, outliers, randomized quantile residual, zero adjusted regression models.

1 Introduction

Zero adjusted regression models are used to fit variables that are discrete at zero and continuous at some interval of the positive real numbers. They are useful in many areas such as insurance (Bortoluzzo et al.,, 2011), microbiology (Rocha et al.,, 2017), credit risk (Tong et al.,, 2016), hydrology (Serinaldi & Kilsby,, 2014), biodiversity (Rubec et al.,, 2016) and geology (DeCarlo,, 2015). When the support of the response variable is the interval $[0;\infty[$, some possible models that can be used are the zero adjusted gamma regression model (Tong et al.,, 2013), the zero adjusted gaussian inverse regression model (Heller et al.,, 2006) and the zero adjusted Birnbaum-Saunders regression model (Tomazella et al.,, 2018). The zero adjusted beta regression model (Cook et al.,, 2008; Ospina & Ferrari,, 2012), also known as the zero inflated beta regression model, is usually used when the response variable assumes values in the interval $[0;1[$.

Residuals play an important role in checking model adequacy and in the identification of outliers and influential observations. For each observation of the sample data, a residual is a measure of disagreement between the observed and the fitted value by a regression model (Cook & Weisberg,, 1982). When the distribution of the residual is well approximated by a known probability distribution, we can classify an observation as a possible outlier if the absolute value of its residual exceeds a threshold. Otherwise, an observation is classified as a possible outlier if the absolute value of its residual is considerably greater than that of almost all observations. Since outliers are observations poorly fitted by the model, it is important to detect them for many reasons. First, outliers can be a result of an operational error that should be solved before using data and the regression model to make inference. Second, especially if many outliers are identified, they can suggest that the fitted regression model is not adequate for these data. Finally, outliers are usually also influential, which are observations that cause major changes in inferential analysis results (Kutner et al.,, 2004).

In zero adjusted regression models, there are two fitted values for each observation that are the main interest. The first is the estimate of the probability of the observation assuming the zero value. The other is the estimate of the expected value of the response variable given that the observation assumed a non-zero value. For observations that assumed the zero value, only the former is useful, but both are important for the other observations. As the number of fitted values of interest is different according to the assumed value, it is difficult to define a single residual for all observations.

The randomized quantile residual (Dunn & Smyth,, 1996) is useful for checking the overall adequacy of a zero adjusted regression model (Ospina & Ferrari,, 2012) and can be obtained for all observations. However, it may fail to identify some outliers when the probability of the response variable assuming positive values close to zero is very small (see Section 2). Therefore, for outlier identification, it is convenient to define a different residual for the discrete and for the continuous component of the model. For the discrete component, a residual defined for regression models for binary dependent variables (Hosmer & Lemeshow,, 2000, Chapter 5) can be used. Similarly, for the continuous component, we can use a residual commonly used in regression models for continuous response variables (McCullagh & Nelder,, 1989). However, for observations that assumed a non-zero value, this strategy is not appropriate because the residual is not a function of the estimate of the probability of the observation assuming the zero value.

For the zero adjusted beta regression model, Ospina & Ferrari, (2012) proposed a residual that is a function of a residual proposed by Espinheira et al., (2008) for beta regression and the estimate of the probability of the observation assuming the zero value. Nonetheless, simulation studies (see Section 3) suggest that this residual is not good for outlier identification in the zero adjusted beta regression model.

In this work, we introduce a class of residuals for outlier identification in zero adjusted regression models. This new class is a function of a residual defined for the continuous component of the model and the maximum likelihood estimate of the probability of the observation assuming the zero value. Simulation studies (Section 3) suggest that one of the residuals of this class is adequate for outlier identification in this class of regression models. Moreover, this residual detects outliers that are not identified by the randomized quantile residual (see Section 4).

The remainder of this paper is organized as follows. Section 2 presents a new class of residuals and its properties. In Section 3, Monte Carlo simulation studies are performed to examine the properties of this class of residuals. In the following section, an application to real data is shown considering the zero adjusted beta regression model. Concluding remarks are provided in Section 5.

2 A residual for outlier identification

In this section, we define the zero adjusted regression models and introduce a class of residuals for this class of models.

2.1 Zero adjusted regression models

A zero adjusted distributed random variable has probability density function (p.d.f.) given by

[TABLE]

where $\alpha=\textrm{Pr}(Y=0)$ is a mixture parameter, $f_{W}(y;\mu,\phi)$ is the p.d.f. of a continuous distribution with support on $\mathcal{Y}\subset\mathbb{R}^{+}$, and $\mu$ and $\phi$ are, respectively, location and dispersion (or precision) parameters of $f_{W}(y;\mu,\phi)$. Using this p.d.f., the class of regression models considered in this paper is defined as following.

Definition 2.1**.**

Let $Y_{1},Y_{2},\ldots,Y_{n}$ be independent random variables, where $Y_{i}$ is zero adjusted distributed with parameter vector $\theta_{i}=(\mu_{i},\phi_{i},\alpha_{i})^{\top}$. The zero adjusted regression (ZAR) models are defined by (2.1) and the following systematic components

[TABLE]

where $\eta_{ij}=x_{ij}^{\top}\beta_{j}$ are linear predictors, $\beta_{j}$ is a $p_{j}\times 1$ vector of unknown parameters, $x_{ij}^{\top}=(x_{ij1},x_{ij2},\ldots,x_{ijp_{j}})$ represents the values of $p_{j}$ predictor variables and the link functions $g_{j}(\cdot)$ are strictly monotone and at least twice differentiable.

Some particular cases of the model (2.2) are the zero adjusted beta (ZABE) regression model (Cook et al.,, 2008; Ospina & Ferrari,, 2012), the zero adjusted gamma (ZAGA) regression model (Tong et al.,, 2013), the zero adjusted gaussian inverse (ZAIG) regression model (Heller et al.,, 2006) and the zero adjusted reparameterized Birnbaum-Saunders (ZARBS) regression model (Tomazella et al.,, 2018). The model (2.2) is also a particular case of the generalized additive models for location, scale and shape (Rigby & Stasinopoulos,, 2005). Parameters of the ZAR model can be estimated by maximum likelihood using a numerical nonlinear optimization algorithm (Nocedal & Wright,, 2006, Section 6.1). Confidence intervals regarding the parameters can be obtained using the asymptotic properties of the maximum likelihood estimator and hypothesis testing can be conducted using Wald, likelihood ratio or score statistic (Sen et al.,, 2010, Section 8.6).

2.2 The randomized quantile residual

The randomized quantile residual was introduced by Dunn & Smyth, (1996). For the ZAR models, it is defined as

[TABLE]

where $\hat{\theta}_{i}$ is the maximum likelihood estimators of $\theta_{i}$, $\Phi(\cdot)$ is the cumulative distribution function of the standard normal distribution, $F(\cdot)$ is the cumulative distribution function of the response variable and $u_{i}$ is a uniformly distributed random variable on the interval $(0,\hat{\alpha}_{i})$. The residual $r_{i}^{q}$ is interesting for checking the overall adequacy of a ZAR model when sample size is not small, because it is asymptotically standard normally distributed.

Note from (2.3) that $r_{i}^{q}\geq\Phi^{-1}(\hat{\alpha}_{i})$ if $y_{i}>0$, because, in this case, $F(y_{i};\hat{\theta}_{i})\geq F(0;\hat{\theta}_{i})=\hat{\alpha}_{i}$. Therefore, if $y_{i}$ is close to zero and $\hat{\alpha}_{i}$ is not very small and below 0.5, $r_{i}^{q}$ will be negative and its absolute value will not be large. For this reason, $r_{i}^{q}$ may fail to identify some outliers in ZAR models. Suppose for example that a commercial bank want to model the amount of investment that customers have in the bank twelve month after opening a checking account. An important predictor variable in this case is the customer income. If a customer have a very high income but an amount of investment close to zero, this observation should be classified as an outlier. However, using $r_{i}^{q}$, this observation will not be classified as an outlier if $\hat{\alpha}_{i}$ is not very small. Rich customers may have $\hat{\alpha}_{i}$ not so close to zero, because many of them may prefer to invest their money on investment banks.

2.3 A residual for outlier identification in ZAR models

Considering the limitations of $r_{i}^{q}$, for outlier identification, it is interesting to use a different residual for the discrete and for the continuous component of the model. As mentioned earlier, for the discrete component, a residual defined for regression models for binary dependent variables (Hosmer & Lemeshow,, 2000, Chapter 5) can be used. Similarly, for the continuous component, we can use a residual commonly used in generalized linear models if the fitted model is the ZAGA or the ZAIG regression model or a residual usually used in beta regression if the fitted model is the ZABE regression model. However, for observations that assumed a non-zero value, the residual will not be a function of the estimate of the probability of the observation assuming the zero value and therefore it may not be adequate in some situations. Suppose for example, that we are fitting a ZAGA regression model and we are using the deviance residual for the continuous component of the model. Assume also that, for observations $s$ and $t$, $y_{s}$ and $y_{t}$ are positive and $r_{t}^{d}=r_{s}^{d}=3.5$, where $r_{i}^{d}$ is the deviance residual for observation $i$ considering the continuous component of the ZAGA regression model. In these case, both observations would be considered equally discrepant. Nonetheless, if $\hat{\alpha}_{s}$ is considerably greater than $\hat{\alpha}_{t}$ then observation $s$ should be considered more discrepant than observation $t$.

According to the previous considerations, it seems adequate to define a residual that is a function of $\hat{\alpha}_{i}$ and $r_{i}$, where $r_{i}$ is any residual used in regression models for continuous response variables. The estimate $\hat{\alpha}_{i}$ is greater than zero and lower than one and $r_{i}$ can assume any real number. To make them comparable and inspired by the randomized quantile residual, we propose to apply the standard normal distribution function on the absolute value of $r_{i}$. Therefore, we propose the following residual for outlier identification in ZAR models:

[TABLE]

The residual $r_{i}^{\star}$ is easily rewritten as

[TABLE]

Note that (2.4) defines a class of residuals for ZAR models, since $r_{i}$ is any residual used in regression models for continuous response variables (e.g. deviance residual).

A natural choice for $r_{i}$ is the quantile residual. We will use the term quantile residual and not randomized quantile residual to make clear that it refers to the continuous component of the ZAR model, in which it is not necessary to generate uniformly distributed random variables. We denote $r_{i}^{\star}$ as $r_{i}^{\star q}$ when $r_{i}$ is the quantile residual and we named it as zero adjusted quantile residual (ZAQR). For $r_{i}^{\star q}$, the following theorem holds. The proof is in appendix.

Theorem 2.1**.**

If $r_{i}^{q}>0$, then $r_{i}^{\star q}=r_{i}^{q}$.

Theorem 2.1 suggests that $r_{i}^{\star q}$ can be interpreted as a correction for outlier identification of the randomized quantile residual when $r_{i}^{q}<0$.

Another result about $r_{i}^{\star}$ also proved in appendix is stated in the following theorem.

Theorem 2.2**.**

If $r_{i}$ is standard normally distributed and if $\alpha_{i}$ is known $\forall i$ then $\forall k>\Phi^{-1}(1-0.5(1-\alpha_{i}))=\Phi^{-1}(0.5+0.5\alpha_{i})$,

[TABLE]

According to Theorem 2.2, if $r_{i}\sim N(0,1)$ and $\alpha_{i}$ is known, the residual $r_{i}^{\star}$ is quite interesting because it has a behavior similar to a normally distributed random variable. A standard normally distributed residual is desirable because it facilitates the definition of a threshold to the residuals. However, in general, the distribution of $r_{i}$ is not normal (not even the quantile residual is standard normally distributed in small samples (Scudilio & Pereira,, 2017)) and the parameter $\alpha_{i}$ is unknown. Thus, to study if $\textrm{Pr}(r_{i}^{\star}>k)$ is close to $1-\Phi(k)$ in practical situations, it is necessary to use Monte Carlo simulation studies.

3 Simulation studies

We performed Monte Carlo simulation studies considering the ZABE regression model and using the Ox language (Doornik,, 2009). For each case, we used 25,000 replications with $n=100$.

In the first scenario, we considered the following ZABE regression model:

[TABLE]

where $\beta_{11}=-1.5$, $\beta_{12}=-1.0$, $\beta_{13}=1.0$, $\beta_{21}=4.0$, $\beta_{31}=-0.5$, $\beta_{32}=0.5$ and $\beta_{33}=-1.0$, which results in $\mu\in(0.076,0.378)$, $\phi=54.6$ and $\alpha\in(0.182,0.500)$. The explanatory variables $x_{i12}$ and $x_{i13}$ were generated as independent draws from the the standard uniform distribution and we assumed $x_{i32}=x_{i12}$ and $x_{i33}=x_{i13}$. They remained constant throughout the simulations.

In the simulation studies, we compare $r_{i}^{\star q}$ with four other residuals. The first three belongs to the class introduced in this paper, but using other residuals as $r_{i}$ instead of the quantile residual. We used the standardized weighted residual 1 (Espinheira et al.,, 2008), the standardized weighted residual 2 (Espinheira et al.,, 2008) and the adjusted standardized weighted residual 1 (Anholeto et al.,, 2014) and they are denoted, respectively, by $r_{i}^{\star efc1}$, $r_{i}^{\star efc2}$ and $r_{i}^{\star asb}$. The last one, $r_{i}^{of}$, is the residual proposed by Ospina & Ferrari, (2012) mentioned earlier.

For each of the observations and residuals, we calculated the percentage of the residuals smaller than $-3,-2,-1$ and greater than $1,2,3$ among the $25,000$ replications. When $y_{i}=0$, the residuals are not defined and we considered that they are neither smaller than -1 and nor greater than 1. The goal of the simulation studies is to verify if the behavior of the residuals considered in the analysis are similar to the standard normal distribution. For this reason, we calculated descriptive statistics for the $n$ observations of each of the percentages mentioned above and compared them with the theoretical values of the standard normal distribution.

Table 1 presents the simulation results for Scenario 1. The mean and median of the calculated percentages for $r_{i}^{\star q}$ are close to the theoretical values of the standard normal distribution for all intervals. If we consider sampling variability, even the minimum and the maximum of the calculated percentages for $r_{i}^{\star q}$ are not so far from the theoretical values of the standard normal distribution. On the other hand, none of the other residuals has mean and median of the calculated percentages close to the theoretical values of the standard normal distribution for the 6 considered intervals.

Many other scenarios were evaluated changing the values of the parameters, the distribution of the explanatory variables, sample size and including explanatory variables in the submodel for the parameter $\phi$. Tables were omitted for the sake of brevity. In all scenarios, the distribution of $r_{i}^{\star q}$ is closer to the standard normal distribution than that of the other residuals. Additionally, in general, the results for $r_{i}^{\star q}$ are similar to Scenario 1. However, when sample sizes reduces to 50 or when we include explanatory variables in the submodel for the parameter $\phi$, the calculated percentages for $r_{i}^{\star q}$ are slightly far from the theoretical values of the standard normal distribution than in Scenario 1.

We also performed Monte Carlo simulation studies considering the ZAIG and the ZAGA regression model. Here, we considered as $r_{i}$ the quantile residual and four others: deviance residual (Davison & Gigli,, 1989), the Pearson residual (McCullagh & Nelder,, 1989), the Williams residual (Williams,, 1987) and the Anscombe residual (Pierce & Schafer,, 1986). They are denoted as $r_{i}^{*dev}$, $r_{i}^{*pea}$, $r_{i}^{*wil}$and $r_{i}^{*ans}$, respectively. We compared these residuals with $r_{i}^{**q}$, a natural extension to ZAGA and ZAIG regression model of the residual proposed by Ospina & Ferrari, (2012) for the ZABE regression model, but here using the quantile residual.

Many scenarios were also considered, but for the sake of brevity we only present the results of Scenario 1 for the ZAIG regression model (Table 2). In this scenario, the explanatory variables were generated as independent draws from the the standard uniform distribution and the values of the parameters result in $\mu\in(20.9,403.4)$, $\phi=0.02$ and $\alpha\in(0.27,0.62)$. As in the ZABE regression model, the calculated percentages for $r_{i}^{\star q}$ are close to the theoretical values of the standard normal distribution for all intervals. On the other hand, the other residuals has calculated percentages far from standard normal values for some intervals. In the other scenarios and in the scenarios for the ZAGA regression model, the distribution of $r_{i}^{\star q}$ is also closer to the standard normal distribution than that of the other residuals. In addition, in all of the scenarios for ZAIG or ZAGA regression model, the calculated percentages are not far from the theoretical values of the standard normal distribution.

4 Application

In this section we use data from 1000 applicants of the National Exam of Upper Secondary Education (ENEM) from 2014. The brazilian exam was created in 1998 as an alternative to higher education entrance examinations and focus on a comprehensive and structural concept of human intelligence. Data were obtained from http://inep.gov.br/web/guest/microdados.

The exam has 5 different tests: natural sciences, human sciences, language, mathematics and writing essay. The score of an applicant on the writing essay test is limited on the interval $[0;1000]$. Here we considered the score on the writing essay test divided by 1000 and our goal is to study the relation of this response variable and the scores of the other tests, the gender and the age of the applicant. In our sample, none of the applicants reached the maximum score and 51 obtained a zero score. Based on this information, the ZABE regression model is a natural choice to fit the response variable. Considering the results of exploratory analysis, the age of the applicant was categorized in two levels: greater than 25 and others.

All variables were initially included in the three submodels of the ZABE regression model. Variables were selected for the final model in a backward way using Wald tests. Figure 1 presents a half-normal residual plot with simulated envelope (Atkinson,, 1985, Section 4.2) for the final model using $r_{i}^{q}$. The plot does not suggest model misspecification.

To perform outlier identification, Figure 2 presents plots of residuals against $\mu_{i}$ considering $r_{i}^{q}$ (left) and $r_{i}^{\star q}$ (right) and using the 949 applicants that obtained a score in the writing essay test greater than zero. None of the observations has absolute value of the residuals considerable greater than the others. However, using $r_{i}^{\star q}$, four observations stand out. Three of them also emerge when $r_{i}^{q}$ is used, but observation #242 does not have a high absolute value of $r_{i}^{q}$. This applicant obtained score 0.12 on the writing essay test, considerably lower than the fitted value $\hat{\mu}_{242}=0.437$. However, as $\hat{\alpha}_{242}$ is not very low $(0.037)$, the absolute value of $r_{i}^{q}$ is not high for this applicant and this residual fails to identify one of the worst fitted observations.

We fitted the final ZABE regression model excluding, one at a time, the four observations pointed out in Figure 2 and compared the estimates of the parameters of each of these models with that obtained using all observations. The changes in the submodel for the parameter $\mu$ are lower than $6\%$. In addition, only when we exclude observation #662, the change in one of the estimates in the submodel for the parameter $\phi$ is greater than $8\%$, and it is not so high $(13.3\%)$. As we considered in Figure 2 only the 949 applicants that obtained a score in the writing essay test greater than zero, the changes in the estimates in the submodel for the parameter $\alpha$ are very small. We also analyzed the observation with zero value in the response variable and none of these observations have absolute value of the deviance residual greater than 3. Therefore, there are not observations that affect considerably the results of the final ZABE regression model.

Assume now that the value of the response variable for observation #242 was incorrectly recorded as 0.001, the minimum possible score greater than zero. In this case, after fitting the ZABE regression model, we obtain $r_{i}^{q}=-1.81$ and $r_{i}^{\star q}=-7.86$. Note that $r_{i}^{\star q}$ clearly identify the observation as outlier, while $r_{i}^{q}$ fails to do it. Considering the ZABE regression model, the estimate of the probability of an observation assuming a value greater than zero and lower than 0.01 if it has the covariates of the observation #242 is very small $(1.7\times 10^{-9})$. Excluding this observation, four estimates of parameters change more than $14\%$, two of them are in submodel for $\mu$, which is the main one in a ZABE regression model. In addition, the change in the estimate of the parameter associated with human sciences score in the submodel for $\phi$ is very high $(52\%)$. Based on these information, observation #242 should be identified as an outlier, which could enable to an applied statistician identifying that this observations was incorrectly recorded. This example illustrates that $r_{i}^{\star q}$ is better than $r_{i}^{q}$ for outlier identification in zero adjusted regression models.

Table 3 presents the results of the final ZABE regression model for the writing essay score. As expected, the mean of the continuous component of the distribution of the writing essay score increases as the score in the other tests increase. However, only for two out of the four other tests, there are evidences that an increase in the score of these tests, reduces the probability of the applicant obtaining a zero score in the writing essay test. In addition, older applicants has a lower mean and a higher probability of obtaining a zero score, while female applicants has a better mean performance.

5 Concluding remarks

In this work, we introduced a class of residuals for zero adjusted regression models. Monte Carlo simulation studies and an application were performed to compare these class of residuals with other residuals that can be used to perform diagnostic analysis in these regression models. The simulation studies suggested that one of the residual of this new class, named ZAQR, has some similar properties to a standard normally distributed variable. In addition, the application indicated that ZAQR can identify outliers in situations in which randomized quantile residual fails to do it. For these reasons, ZAQR seems to be the best option for outlier identification in zero adjusted regression models.

It is also interesting to mention that the randomized quantile residual may also fail to identify some outliers in regression models in which the response variable is discrete in more than one point, as the zero-and-one adjusted beta regression model (Ospina,, 2008; Stasinopoulos et al.,, 2007) and the three-point adjusted beta regression model (Pereira et al.,, 2013; Gray & Alava,, 2018). Therefore, future works can be developed extending ZAQR for outlier identification in these class of models.

Acknowledgments

This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brazil (CAPES) - Finance Code 001 and by São Paulo Research Foundation (FAPESP) - Grant #2013/17876-6.

Appendix A Appendix: Proof of theorems of Section 2.3

In this appendix, we prove the theorems presented in Section 2.3.

A.1 Proof of Theorem 2.1

If $r_{i}^{q}>0$ then

[TABLE]

A.2 Proof of Theorem 2.2

If $r_{i}$ is normally distributed with mean 0 and variance 1 and if $\alpha_{i}$ is known $\forall i$ then $\forall k>\Phi^{-1}(1-0.5(1-\alpha_{i}))=\Phi^{-1}(0.5+0.5\alpha_{i})$, then

[TABLE]

In a similar way, it is shown that $\textrm{Pr}(r_{i}^{\star}>k)=1-\Phi(k)$.