Outlier detection in gamma regression using Pearson residuals: Simulation and an application

1.   Introduction

Regression analysis is the main tool to study the relationship between and dependence of one or more variables on one or more independent variables. These relationships can be observed in every research area. Regression analysis has a wide variety of applications in different fields ^[1,2,3,4] and Sarhan et al. ^[5]. Regression results are reliable only if the quality of data is good and the selected regression model is correct. If the quality of data and the regression model are not appropriate, then the results of measuring such relationships are incorrect. The quality of data refers to the outlier free data set ^[2,6]. The correct regression model refers to identifying the distribution of the response variable.

An outlier is a point that is far from the rest of the data points ^[7,8,9]. In the regression context, Desgagné ^[10] stated that the observations with more distant regression errors conflict with most of the errors originating from the assumed normal distribution. An outlier may or may not affect the regression inferences. An outlier may be in one or multiple variables. Outlier detection in univariate analysis has been done by several researchers i.e. many studies in the literature have focused on univariate analysis ^{[7,8,9,11,12]}. Outlier diagnostics in the linear regression model (LRM) has also gained much attention from researchers ^[13,14]. Balasooriya et al. ^[15] compared some well-known outlier detection methods by using the LRM. They concluded that all methods do not agree with each other for the detection of outliers.

Regression analysis is used to determine the model for forecasting/prediction purposes. There is a variety of regression models, e.g., LRMs, generalized linear models (GLMs) and non-linear models, Gamma regression (GR) is employed, when the distribution of the dependent variable is gamma. GR has a variety of applications in the literature with examples in health sciences, industries and environment, for more details see ^{[16,17,18,19,20,21,22,23]}.

Outlier detection using univariate gamma response without considering any independent variable is also available in the literature ^{[24,25,26,27,28]}. Shayib and Young ^[29] first studied the extreme residuals in GR and proposed the Pearson and Anscombe residuals with modified forms; they concluded that the modified forms of these residuals are not good.

The detection of outliers in the GR model has not been addressed in the literature. This paper deals with outlier detection in the GR model by using a new approach i.e., an adjusted Pearson residuals (PRs) approach. Outlier detection with the adjusted form of residuals (other than PRs) was first studied by Tiao and Guttman ^[30]. This approach was also proposed for some of the GLM responses. Cordeiro ^[31] introduced the adjusted PRs for the Poisson regression model. In addition, adjusted Wald residuals and PRs for beta regression have been suggested by various authors ^[32,33]. In these studies, most of the researchers focused on the adjusted PRs probability distributions. They observed that the adjusted residuals performed better than others.

The main objectives of the current research were to propose some outlier diagnostics based upon Pearson (standardized and adjusted) residuals in GR, modify some available LRM-based outlier detection methods for GR and then make a comparison of these modified and proposed outlier detection methods with the help of simulations and a real data set.

2.   Materials and methods

The probability density function of the gamma response variable (y) is given by

The mean and variance of Eq (2.1) are respectively given as $E(y) = \nu\tau$ and $Var(y) = \nu\tau^{2}$. According to Hardin and Hilbe ^[34], Eq (2.1) can be transformed with parameters $\nu = \phi^{- 1}$ and $\tau = \mu\phi$. Given these parameters, the gamma density for y is given by

The mean and variance of Eq (2.2) are respectively given as $E(y) = \mu$ and $Var(y) = \phi\mu^{2}$. For the ith observation, let $x_{i1}, \ldots x_{\text{ip}}$ represent the p independent variables. Then, the GR for the mean of the response variable y is given by

where ${\bf{x}}_{\bf{i}}^{{\bf{T}}} = \left(1, x_{i1}, \ldots, x_{\text{ip}} \right), \ \beta = \left(\beta_{0}, \beta_{1}, \ldots, \beta_{p} \right)$ is the vector of regression coefficients including intercept and g(.) is the link function. This link function in the GR can be reciprocal or log.

Let $l_{i}$ be the log-likelihood function of the response variable of Eq (2.2) which is mathematically defined by

Let $\widehat{ \beta ,} \ \ \widehat{ \mu }$ and $\widehat{\phi}$ be the maximum likelihood estimates (MLEs) which are obtained by maximizing the log-likelihood of Eq (2.3) using the Netwon-Raphson iterative method. The MLE of $\beta$ is computed by solving the system of equations. For this purpose, we equate the first derivative of Eq (2.3) to zero; then, we have

where $U\left( \beta \right)\ $is the score vector of the order$\ (p + 1) \times 1$. Since Eq (2.4) is nonlinear in $\beta$, Newton-Raphson methods can be employed for the estimation of $\beta$ ^[34]. Suppose $\beta ^{m}$ is the approximated MLE of $\beta$ at the $m$th iteration; then, the iterative reweighted method ^[35] gives the following expression

where $I\left( \beta ^{m} \right)$ is the $(p + 1) \times (p + 1)$ fisher information matrix at the $m$th iteration. Applying convergence in deviance to Eq (2.5), the unknown parameters can be computed as

where ${\widehat{z}}_{i} = {\widehat{\eta}}_{i} + \frac{y_{i} - {\widehat{\mu}}_{i}}{{\widehat{\mu}}_{i}^{2}}$ is the adjusted response variable, $\widehat{{\bf{W}}} = diag\left({\widehat{\mu}}_{1}^{2}, {\ldots, \widehat{\mu}}_{n}^{2} \right)$ and ${\widehat{\mu}}_{i} = \frac{1}{x_{i}^{T}\widehat{ \beta }}.$

Several types of the GLM residuals are available in the literature ^[33] but we consider the most popular PR for the detection of an outlier.

The PRs for GR are defined by

The standardized PRs are characterized as

where ${h}_{\text{ii}} = diag\left({\bf{H}} = {\widehat{\bf{w}}}^{\frac{\bf{1}}{\bf{2}}}{\bf{X(}}{\bf{X}}^{{\bf{T}}}\widehat{{\bf{W}}}{\bf{X}}\bf{)}^{\bf{- 1}}{\bf{X}}^{{\bf{T}}}{\widehat{{\bf{W}}}}^{\frac{\bf{1}}{\bf{2}}}\right)$.

Generally, $E\left(\chi_{i} \right)$ to order $O\left(n^{- 1} \right)\ $ does not converge to zero and $\text{Var}\left(\chi_{i} \right)$ does not tend to one. Here n is the sample size. To handle such a situation, we require some adjustments to these residuals. To do this, Cox and Snell ^[36] obtained some matrix formulae for the adjusted residuals.

Various criteria are available in the literature ^[2,12] for testing the quality of regression models. These criteria include mean quadratic error prediction (MEP), the Akaike information criterion (AIC), standard errors and coefficients of determination. As our study is concerned with the GR model, we consider some different criteria for testing the goodness of the GR model after diagnosing the outlying points. These include the Pearson chi-square statistic $\left(\chi^{2} \right)$, MEP, the AIC, Efron's pseudo r-squared criterion $\left(R^{2}_{Efron} \right)$ and the dispersion parameter $\left(\hat{\phi} \right)$; these criteria are computed according to the following relations:

$\chi^{2} = \sum_{i = 1}^{n} \chi^{2}_{i}$, where $\chi_{i}$ is the ith PR defined in Eq (2.7).

$MEP = \sum_{i = 1}^{n}\dfrac{V^{1/2}\left(\mu_{i}\right)\chi^{2}_{i}}{\left(1-h_{ii}\right)n}$, where $V\left(\mu_{i}\right) = \mu_{i}^{2}$ is the variance function of the GR model.

$AIC = \dfrac {-l_{i}\left(M_{K}\right)+2p}{n}$, where $l_{i}$ is the log likelihood function of the GR model defined in Eq (2.3).

$R^{2}_{Efron} = 1- \dfrac{\sum_{i = 1}^{n}\left(y_{i}-\hat{\mu}_{i}\right)^{2}}{\sum_{i = 1}^{n}\left(y_{i}-\bar{y}\right)^{2}}$,

$\hat{\phi} = \dfrac{\chi^{2}}{n-p-1}$.

3.   Outlier detection method

This section comprises two subsections. In the first subsection, the proposed outlier detection methods based on the PRs of the GR model are presented. The second one comprises a review of some existing outlier detection methods.

3.1.   Proposed outlier detection methods

The PRs have a significant contribution in regression diagnostics. Here we propose an outlier detection method based upon PRs.

3.1.1.   Standardized PRs

In regression analysis, the standardized residuals are generally used for the detection of outliers. So, here we consider standardized PRs, which have been defined in Section 2 by Eq (2.8) as

On the basis of standardized PRs, the ith point is considered to be an outlier if $\left|\chi_{i}^{'}\right|>3$.

3.1.2.   Jackknife PRs

There are some analytical methods that are used for the detection of outliers. One of these analytical methods is the use of jackknife residuals. Cook and Weisberg ^[37] suggested that the outliers can be detected with the help of jackknife residuals. They defined the jackknife residuals for LRMs as

where $r_{i}$ is the standardized residual of the LRMs. The decision rule for the detection outlier is that if $\left|e_{Ji} \right|>t_{\frac{\alpha}{2n}}$ (with $n-p-1$ degrees of freedom), then there is an indication for the existence of outliers. The application of these residuals for the detection of outliers in chemometrics with reference to an LRM has been studied by Meloun and Militky ^[2]. Now, we modify Eq (3.2) for the GR by following Amin et al. ^[1] and obtaining

To identify outliers of the GR, we propose the cut-off point for the jackknife PR to be if $\left| \chi_{J, i} \right| > t_{(1 - \alpha)}(n - p - 1)$, then the ith observation is declared as an outlier t is the student t-distribution with $(n-p-1)$ degrees of freedom.

3.1.3.   Adjusted PRs

Amin et al. ^[1] proposed the adjusted PRs in the inverse Gaussian regression for the detection of single influential points in chemometrics. These residuals are proposed for the GR as follows:

where $r_{i} = \left(E\left(R_{i} \right) \right)^{T} = - \frac{\sqrt{\phi}}{2}\left({\bf{I - H}} \right)J{\bf{z}}$, where ${\bf{H}} = {\bf{W}}^{\frac{\bf{1}}{\bf{2}}}{\bf{X(}}{\bf{X}}^{{\bf{T}}}{\bf{WX}}{\bf{)}}^{\bf{- 1}}{\bf{X}}^{{\bf{T}}}{\bf{W}}^{\frac{\bf{1}}{\bf{2}}},$ ${\bf{J}} = diag\left(2\mu^{2} \right)$ and ${\bf{z}} = (z_{11}, \ldots, z_{\text{nn}})^{T}$is a vector; also, $v_{i} = \left(\text{Var}\left(R_{i} \right) \right)^{T} = 1 + \frac{\phi}{2}\left({\bf{QHJ}} - {\bf{T}} \right){\bf{z}}\ $, where ${\bf{Q}} = diag(2)$ and ${\bf{T}} = diag\left(\left(2\phi^{- 1} + 6 \right)\mu^{2} \right)$.

Note that ${\widehat{r}}_{i}\ and\ {\widehat{v}}_{i}$ are computed by using ${\widehat{\mu}}_{i}$ instead of $\mu_{i}$.

Amin et al. ^[1] stated that the adjusted residuals can be used for the detection of an outlier. The decision rule to declare the ith observation is an outlier is that $\left| \chi_{i}^{A} \right| > 2.$

3.2.   Available outlier detection methods

In the literature, numerous methods have been recommended for the detection of outliers. We consider a few of them for comparison with the adjusted PR for the GR model.

3.2.1.   Z-method

For the identification of outliers in univariate cases, the Z-score is introduced based on the median and inter-quartile range (IQR) as below:

The Z-score method declares that the ith observation is an outlier if $|Z| > 3$.

3.2.2.   Modified Z-method

The modified Z-statistic (MZS) has been proposed for the detection of outliers based on the median ^[38]; MZS is defined as

where $\text{Median}|.|$ represents the median absolute deviation from the median. One can conclude that the ith observation is an outlier if $Z_{i}^{*} > 3.50$.

3.2.3.   Grubb's test

This test was introduced by Grubbs ^[13] for the detection of a single outlier in the univariate response variable. The Grubb's (G) test only detects single outliers. Therefore, it suspects that most of the observations are outliers. The G statistic can be defined as $G = \frac{\text{Max}\left| y_{i} - \overline{y} \right|}{S}$, where $s$ is the sample standard deviation. For our assumed model, the G statistic is modified as

The decision rule of the G statistic for detecting outliers is given as G > $\sqrt{\frac{t_{\left(\frac{\mathfrak{a}}{2n} \right)}^{2}, n - 2}{n - 2 + (t_{\left(\frac{\mathfrak{a}}{2n} \right)}^{2}, n - 2)}}$, where $\mathfrak{a}$ is the level of significance and n represents the sample size. The above decision rule may be unable to find the appropriate outlier. So, we propose another decision rule for Grubb's method to diagnose an outlier a GR model. For the ith data point, if $G \geq 2$, then this data point is declared as an outlier.

4.   Results and discussion

A comparison of the proposed methods of outlier detection with already available methods through the use of simulations and a real-life data set is presented in this section.

4.1.   Simulation study

This section explains the simulation experiment conducted to study the performance of two types of GR residuals, i.e., the standardized and adjusted residuals to detect outliers. We generated the dependent variable of the GR model with reciprocal link function to be $y_{i} \sim \text{Gamma}\left(\mu_{i}, \phi \right); i = 1, ..., n$, where $\mu_{i} = \left(\beta_{0} + \beta_{1}x_{i1} +... + \beta_{4}x_{i4} \right)^{- 1}$, where $x_{i}$'s are generated from two probability distributions, i.e., uniform (0, 1) and standard normal distributions. The true values of the regression parameters vector $\beta$ were selected as the normalized eigenvector corresponding to the largest eigenvalue of the ${\bf{X}}^{T}\bf{\hat{W}}{\bf{X}}$ matrix such that $\beta^{T}\beta$ = 1 ^[39]. The sample sizes were generated as $n = 25, 50, 75,100,125,150,175 \ and \ 200$ and we took on the values of dispersion, i.e., $\phi = 0.33, 0.67 \ and \ 2$. A single outlier was generated for the dependent variable, i.e., the $8^{\text{th}}$ observation was replaced as $y_{8} = y_{8} + a_{0}$, where $a_{0} = \overline{y_{i}} + 3\left(V\left(y_{i} \right) \right)$. The multiple outliers in the dependent variable are generated as $y_{ii} = y_{ii} + a_{0}$, where ii = 8, 15, 20 three outliers. The performance of these diagnostics was assessed by using gamma—produced samples for the identification of a generated outlier(s). The simulated results were computed with the help of the R-statistical language. The simulation study was replicated 1000 times to find the outliers in percentages.

In Table 1, the performance of our proposed methods is gauged for the detection of outliers based on the standard normal generated $x$'s, and by using PRs with dispersion and sample sizes. It can be observed that for $\phi < 1$, the performance of the $\chi_{i}^{A}$ method was better than that of the other methods in diagnosing the generated single outlier. For this dispersion, when sample sizes were increased to 100, the performance of the $\chi_{i}^{'}$, $\chi_{i}^{A}$ and $\chi_{J}$ outlier detection methods was improved. The performance of the $Z$, $Z^{*}$ and $G$ outlier diagnostic methods was not affected by the increase in sample size. It can also be observed that the performance of all outlier diagnostic methods increased with increasing in dispersion. Moreover, when $\phi = 2$ and $n > 50$, the performance of all the diagnostic methods seems to have been identical. This indicates that sample size and dispersion have some significant effect on our proposed method $\chi_{i}^{A}$ in the detection of a single outlying point. The comparison of the outlier diagnostic methods revealed that $\chi_{i}^{A}$ is better than other methods. It can be seen that the detection was better with the proposed method than with all of the other diagnostic methods, when the values of $X$ were generated from $U(0, \ 1)$. For further details, see Table 2. It can be seen that the performances of the available methods were not good as compared to the performance of our proposed method.

When we applied these methods for the detection of multiple outlying points with standard normally generated independent variables, the performance of our proposed method seemed to be better and more consistent than those of the other methods. The detection performance of the other diagnostic methods was reduced to 50%. On the other hand, the outlying diagnostic performance of all methods was reduced to some extent when the independent variables were generated from the uniform distribution. The performance of all outlier diagnostic methods increased rapidly as the dispersion crossed 1.0 (see Tables 1–4). We found that our proposed method performed better than the other methods including the case when the $X$'s were generated from the standard normal distribution instead of the $x$'s being generated from the uniform distribution.

Based on our findings (Tables 1–4), we can rank the methods as the adjusted PR being the first and the Grubb method being the second best method for outlier detection. Moreover, the results show that the performance of outlier detection increases with an increase in sample size. In studying the effect of dispersions on the outlier detection methods, we have found that outlier detection methods are affected directly by the dispersion. This means that upon increasing the dispersions, the outlier detection efficiencies of different methods are increased. The maximum numbers outliers were detected for larger dispersions.

4.2.   Application: ARDENNES dataset

Now, we will evaluate the performance of the proposed methods with the help of a real application. For this purpose, we applied the ARDENNES data taken from Barnard et al. ^[40]. The main use for this data set was to determine the first etch biopsy, i.e., in the beginning of a layer of extracted incisor enamel (y) based on two explanatory variables for the data on 55 children. These explanatory variables included the etched depth (x₁), which was estimated from the amount of calcium removed for the duration of the etch biopsy as the first explanatory variable. The age of the child (x₂) which had been transformed to the decimal system from years and months was considered as the second explanatory variable. Mallet et al. ^[41] applied linear regression and other models to this data set using log (y) as a response variable. After fitting the LRM, they explored the outlier detection analysis of this data and found that the 23^rd, 48^th and 52^nd points were the outliers. However, this data set is not well fitted to the normal distribution since the trend of the dependent variable is positively skewed. From the distribution of the fitting test, we observed that the GR model is well fitted to this data set, the results are reported in Table 5.

So the appropriate regression model to determine the etch biopsy based on these two explanatory variables is the GR model.

The fitted GR is given by

where the square brackets contain the standard errors of the estimated parameters. The letter $N$ represents the non-significance and $S$ represents the significance of the regression coefficients. The fitted GR model achieved the following results MEP = 267.08, AIC = 847.19, and Pearson chi-square statistic = 11.61. After fitting the GR model, we next computed the outlier statistics, which are plotted in Figure 1.

From Figure 1, we can observe that the proposed outlier detection methods detected the 9^th, 23^rd, 48^th and 52^nd points as outliers. The 9^th point was not detected as an outlier in the original work due to misidentification of the regression model but this did not affect the GR estimates. Hossain and Naik ^[42] also indicated that an outlier may or may not affect the regression estimates. So, in this case this outlying point had minimal effect on the GR estimates while $R_{Efron}^{2}$ also decreased because outliers are related to the response variable. $R_{Efron}^{2}$ only improved due to elimination of the influential points which are related to the explanatory variables for more details see ^[5,43,44]. All detected outliers indicated that these children had high lead levels in their enamel, but there were no extreme points in the explanatory variables. Moreover, these outliers were also identified successfully by using the Jackknife residuals and adjusted PRs.

After deleting the identified outliers, the refitted GR model was given by

with $MEP = 169.85$, $AIC = 756.15$ and $Pearson\ chi-square\ statistic = 6.89$. These results indicate how much the variation decreased after deleting the identified outliers, e.g., the Pearson chi-square statistic (variation measure) was reduced to 59% (see Table 6). Another noticeable point is that the two independent variables were found to be significant after deleting the identified outliers. From Table 6, it can be observed that an individually high outlier is the 48^th point, which affected the fitted GR results. Collectively, all detected outlying points affected the GR estimates of $\beta_{0}$ and $\beta_{2}$ respectively.

5.   Conclusions

Outlier detection in regression models is an important step to getting reliable and valid results. These detections methods are based on some diagnostic test statistics, which can be calculated based on the regression results. To detect the outliers, the first and most important step is the choice of an appropriate regression model, because, sometimes, outliers may arise due to an inappropriate regression model. For the selection of an appropriate regression model, one should test the distribution of the response variable. If the probability distribution of the dependent variable is gamma, the appropriate choice of model is the GR model. In this paper, we proposed outlier diagnostics based on the use of the Pearson (standardized and adjusted) residuals in GR model. Some available LRM-based outlier detection methods were modified for the GR model. These modified methods were compared with our proposed outlier detection methods with the help of simulations and a real data set. These results indicate that our proposed methods for the detection of outliers are better than available methods in terms of improving the results of the selected model to enable better decisions in statistics and other disciplines.

Acknowledgments

The authors would like to thank the Deanship of Scientific Research at Umm Al-Qura University for supporting this work by (Grant Code 22UQU4310063DSR05).

Conflict of interest

The authors declare that they have no conflicts of interest regarding the publication of this article.