Robustness against outliers: A new variance inflated regression model for proportions

Abstract

This article addresses the issue of building regression models for bounded responses, which are robust in the presence of outliers. To this end, a new distribution on (0,1) and a regression model based on it are proposed and some properties are derived. The distribution is a mixture of two beta components. One of them, showing a higher variance (variance inflated) is expected to capture outliers. Within a Bayesian approach, an extensive robustness study is performed to compare the new model with three competing ones present in the literature. A broad range of inferential tools are considered, aimed at measuring the influence of various outlier patterns from diverse perspectives. It emerges that the new model displays a better performance in terms of stability of regression coefficients’ posterior distributions and of regression curves under all outlier patterns. Moreover, it exhibits an adequate behaviour under all considered settings, unlike the other models.

Keywords

Bayesian inference Beta regression bounded response Hamiltonian Monte Carlo mixture model outlier

1 Introduction

In many disciplines, the issue of modelling a bounded continuous response variable (e.g., a rate or proportion) as a function of covariates is a relevant one (Kieschnick and McCullough, 2003), some examples being geology, medicine, economics, psychology, environmetrics, psephology, etc. Clearly, linear regression models are unsuitable for these data. To properly address this issue, the beta regression (BReg) model has been proposed and extensively studied (Ferrari and Cribari-Neto, 2004; Paolino, 2001; Smithson and Verkuilen, 2006; Espinheira et al., 2008a,b; Ferrari et al., 2011; Pereira et al., 2013). Moreover, Ospina et al. (2006), Ferrari and Pinheiro (2011) and Bayer and Cribari-Neto (2013) have addressed small sample likelihood asymptotic issues, whereas Cribari-Neto and Souza (2012), Bayer and Cribari-Neto (2017), Simas et al. (2010) and Canterle and Bayer (2017) have addressed testing, model selection, and alternative regression functions. Inflated beta models which account for the presence of mass at 0 and 1 response values are dealt with in Ospina and Ferrari (2010, 2012) and Pereira and Cribari-Neto (2014). For a thorough review of the whole subject, see Liu and Eugenio (2016)

In the basic beta framework, the mean response is related to a linear predictor through a link function. A precision parameter is also present, which may be modelled as a function of covariates too. Despite its variety of shapes, the beta distribution fails to model a wide range of phenomena (Bayes et al., 2012; García et al., 2011) including the heavy tails phenomenon, which is strictly related to the presence of outliers, a critical problem in regression analysis (Rousseeuw and Leroy, 2005). Outliers are observations which are statistically unreasonable with respect to a given probability model (Barnett and Lewis, 1995) and can be interpreted as generated by an alternative random mechanism (Box and Tiao, 1968). Indeed, several authors highlighted the lack of robustness to outliers connected with standard inferential procedures for the BReg model. In particular, Espinheira et al. (2008a,b) adopting a classical approach, propose ways of identifying outliers or influential observations so as to remove them before applying maximum likelihood inference. Instead, Ghosh (2017), again under a classical framework, suggests a different estimation method aimed at producing more stable results.

Bayes et al. (2012) adopt a third and radically different approach by proposing, under a Bayesian setting, a new more sophisticated model specifically devised to account for outliers. Here the aim is not only to achieve point estimate robustness but also to model both the main data body and outliers. Bayes’ model is based on the beta rectangular (BR) distribution (Hahn, 2008), which is a mixture of a uniform and a beta distribution, where the uniform component is expected to capture outliers. A simulation study proved its greater robustness with respect to the BReg model.

To achieve a greater general flexibility, a different mixture regression model can be found in Migliorati et al. (2018), where the response is a flexible beta (FB), that is, a mixture of two betas with common variance and different means (mean-shift type model). This proposal is capable of taking heavy tails, asymmetry, and bimodality into account. Although not specifically meant to address robustness, the flexible beta regression (FBReg) is shown to behave satisfactorily in the presence of some outlier patterns.

In the spirit of Bayes’ work, the main objective of the present article is to introduce and study, within a Bayesian approach, a new mixture regression model for bounded responses capable of handling outliers. To this end, a new beta mixture distribution is proposed and employed for the response variable.

The motivation behind this new regression model comes from some limitations of the two existing mixture ones. In particular, the mean-shift structure of the FBReg model is consistent with the presence of only all upper or all lower outliers, which may not be suitable in some contexts. The beta rectangular regression (BRReg) model does not share this difficulty, but it models the outliers through a fixed component, the uniform, implying a given mean and variance, which may lead to a poor fit.

A natural and simple way to overcome these possible drawbacks is to consider for the outliers a component which has a variable centre, equal to the main data component mean, and a larger (and not fixed) variance. The resulting distribution, called variance inflated beta (VIB), is thus a mixture of two betas with the same mean but different variances.

To understand the potential of the proposed regression model, we perform a comprehensive and detailed robustness comparison of all regression models, aimed at providing an evaluation of the strength of various types of inferential effects produced under different outlier scenarios. This will also allow us to provide practical guidelines for selecting the most appropriate model under various contexts.

In particular, we contaminate a simulated dataset with three different types of outlier patterns, representative of most important cases. In each pattern, outliers are introduced with variable degrees of departure from the rest of the data. To this end, we implement a perturbation operation for the observations specifically designed for constrained data. We then assess outlier influence from a variety of different perspectives. Specifically, we consider their impact on marginal and joint parameter posteriors, fitting measures suitably built to cope with mixture models (widely applicable information criterion WAIC), residuals and predictive intervals. We also propose a measure of the outlier impact on the whole regression curve, as, in some contexts, this gives further complementary information with respect to the regression coefficients’ behaviour.

Differently from the other models, the variance inflated beta regression (VIBReg) will be shown to have a satisfactory performance under almost all the examined conditions.

The rest of the article is organized as follows. The VIB distribution and VIBReg model are introduced in Section 2, where some properties (e.g., identifiability and likelihood boundedness) are derived as well. Details concerning the Bayesian estimation procedure, including a highly efficient Hamiltonian Monte Carlo (HMC) algorithm, are provided in Section 3. Section 4 is devoted to the robustness analysis of the four regression models. An application to a well-known dataset is illustrated in Section 5, whereas Section 6 draws some final conclusions and a discussion of main results.

2 A new variance inflated beta model for bounded continuous variables

2.1 Introduction

In the regression models available in the literature (as well as in the new variance inflated one), the overall mean of the response variable $E [Y_{i}] = μ_{i}$ is directly modelled as a suitable function of covariates. More precisely, following a GLM methodology (McCullagh and Nelder, 1989), a regression model for the mean is defined as:

g (μ_{i}) = x_{i}^{T} β i = 1, \dots, n

(2.1)

where $x_{i}^{T} = (x_{i 0}, x_{i 1}, \dots, x_{ik})$ is a vector of covariates, $β^{T} = (β_{0}, β_{1}, \dots, β_{k})$ is a vector of regression parameters and $g (\cdot)$ is an adequate link function, the most popular one being the logit: $logit (μ_{i}) = \log (μ_{i} / (1 - μ_{i}))$ .

The basic model is the beta regression one BReg( $β, ϕ$ ) (Ferrari and Cribari-Neto, 2004), which assumes that the responses $Y^{T} = (Y_{1}, \dots, Y_{i}, \dots, Y_{n})$ are independent with a beta distribution. In a regression perspective, rather than working with the standard beta parameterization, the mean-precision parameterization is adopted. This contemplates the parameter $0 < μ_{i} < 1$ , which identifies the mean of $Y_{i}$ ( $i = 1, \dots, n$ ) and the precision parameter $ϕ > 0$ .

Estimation issues within the BReg model have been dealt with by maximum likelihood (Ferrari and Cribari-Neto, 2004) and through Bayesian inference (Branscum et al., 2007).

Then, Bayes et al. (2012) proposed the BRReg( $β, ϕ, α$ ) model specifically designed to handle outliers. This model is based on the BR distribution (Hahn, 2008), which is a mixture of a uniform distribution (aimed at capturing outliers) and the beta one. In the BRReg model, the overall regression curve $μ_{i} = g^{- 1} (x_{i}^{T} β)$ can be viewed as a weighted mean of the beta component (which depends on covariates), and a constant identically equal to 0.5 (associated with the uniform component). It is important to note that the weight given to the uniform component is not fixed, but depends on covariates through $μ_{i}$ . This implies, in particular, that outlying observations are deemed to occur more likely for central (i.e., close to 0.5) values of the response variable than for extreme values.

Recently, Migliorati et al. (2018) proposed the FBReg( $β, ϕ, w, p$ ) model where the ith response variable $Y_{i}$ is FB distributed (Ongaro and Migliorati, 2013), that is, it can be viewed as a special mixture of two beta distributions with arbitrary means and common variance.

Here the parameters $0 < w < 1$ and $0 < p < 1$ represent the distance between the means of two components and, respectively, the mixing proportion. Also in this model, the overall regression curve $μ_{i} = g^{- 1} (x_{i}^{T} β)$ can be viewed as a weighted mean of two curves both depending on covariates. One curve dominates the other for any given value of covariates and their difference can be controlled through a specific parameter.

2.2 The variance inflated beta distribution VIB (

μ, ϕ, k, p

)

Our motivation to devise a new regression model, which produces robust inferences in presence of outliers, has originated from the attempt to remove some possible limitations of the BRReg and the FBReg models. The very same structure of the FBReg model favours the presence of only upper or lower outliers, which may not be suitable in all contexts. On the other hand, the BRReg forces the mixture component which captures the outliers to have a fixed mean and variance ( $\frac{1}{2}$ and $\frac{1}{12}$ , respectively), which seems quite a restrictive requirement.

A natural and simple way to remove these limitations is to consider for the outliers a component which has the same centre as the main data component, but which displays a larger (and variable) variance. The resulting distribution (hereafter, VIB $[μ, ϕ, k, p]$ ) has therefore probability density function (pdf) given by:

f_{VIB} (y; μ, ϕ, k, p) = p f_{B}^{*} (y; μ, ϕ k) + (1 - p) f_{B}^{*} (y; μ, ϕ) 0 < y < 1,

(2.2)

where

f_{B}^{*} (y; μ, ϕ) = \frac{Γ (ϕ)}{Γ (μ ϕ) Γ ((1 - μ) ϕ)} y^{μ ϕ - 1} (1 - y)^{(1 - μ)) ϕ - 1}

(2.3)

is the beta density in the mean-precision parametrization.

The parameter $0 < μ < 1$ represents the overall (as well as mixture component) mean, $ϕ > 0$ can be interpreted as a precision parameter since the variance is decreasing with respect to it (see formula 2.5), while $0 < p < 1$ is the mixing proportion. Finally, $0 < k < 1$ determines the extent of the variance inflation.

VIB raw and central moments of any order (as well as the moment generating function) can be easily and explicitly derived. This is because, thanks to the special beta mixture, moments and moment generating function can be written as mixtures of the corresponding quantities of the two beta components, which are available in explicit form. For example, the moment generating function is:

E (e^{tY}) = p_{1} F_{1} (μ ϕ k; ϕ k; t) + (1 - p)_{1} F_{1} (μ ϕ; ϕ; t)

where $_{1} F_{1}$ is the confluent hypergeometric function of the first kind (see e.g., Johnson et al., 1995), and the raw moments are given by:

E (Y^{r}) = p \frac{{[μ ϕ k]}^{(r)}}{{[ϕ k]}^{(r)}} + (1 - p) \frac{{[μ ϕ]}^{(r)}}{ϕ^{(r)}}

(2.4)

where $x^{(r)}$ is the rising factorial defined as $x^{(r)} = \prod_{j = 0}^{r - 1} (x + j)$ .

In particular, the variance of the VIB distribution is equal to:

Var (Y) = μ (1 - μ) \frac{1 + ϕ (k (1 - p) + p)}{(ϕ + 1) (ϕ k + 1)} = p σ_{1}^{2} + (1 - p) σ_{2}^{2},

(2.5)

where $σ_{1}^{2} = μ (1 - μ) / (1 + ϕ k)$ and $σ_{2}^{2} = μ (1 - μ) / (1 + ϕ)$ are the mixture comp- onent variances. The skewness coefficient of the VIB distribution takes the form:

\frac{2 (1 - 2 μ)}{{[1 + ϕ (k (1 - p) + p)]}^{3 / 2}} \sqrt{\frac{(1 + ϕ) (1 + ϕ k)}{μ (1 - μ)}} (\frac{p (1 + ϕ)}{2 + ϕ k} + \frac{(1 - p) (1 + ϕ k)}{2 + ϕ}) .

Note that the sign of the skewness is determined by μ, $μ > 1 / 2$ ( $μ < 1 / 2$ ) denoting positive (negative) asymmetry. Moreover, the skewness coefficient is null if and only if $μ = 1 / 2$ , which is easily seen to be the only case of symmetric VIB distribution.

The kurtosis coefficient is equal to:

\frac{p H (ϕ k) (1 + ϕ)^{2} + (1 - p) H (ϕ) (1 + ϕ k)^{2}}{{[1 + ϕ (k (1 - p) + p)]}^{2}}

where

H (x) = \frac{6}{3 + x} (\frac{(1 - 2 μ)^{2} (1 + x)}{μ (1 - μ) (2 + x)} - 1) + 3 .

(2.6)

The VIB kurtosis formula is involved and for space constraints we do not perform a full study here. Rather, we limit ourselves to show VIB kurtosis greater flexibility than the beta one, in particular as long as high values are concerned. The beta kurtosis, which is given by $H (ϕ)$ (see formula 2.6), can tend to infinity only if $μ \to 0$ or $μ \to 1$ . For any given value of the mean μ, it cannot reach arbitrary high values. For example, if $μ \in (1 / 2 - \sqrt{5} / 10, 1 / 2 + \sqrt{5} / 10)$ , one can show that the supremum of the beta kurtosis is equal to 3. On the contrary, the VIB kurtosis is unbounded for any given μ. For example, it is simple to check that when $ϕ \to \infty$ , $ϕ k \to q$ (where $q$ is a non-negative constant) and $p \to 0$ the kurtosis diverges for any given μ. Analogously, divergence is obtained when $ϕ \to \infty$ , $ϕ k \to \infty$ , $p = k^{τ} \to 0$ , where $0 < τ < 2$ .

Indeed, the VIB distribution extends the shapes of the beta and BR ones mainly in terms of heavy tails. More precisely, any unimodal beta distribution necessarily displays zero limits at zero and at one. Positive and finite limits can be obtained, within the beta distribution, only with monotone (or uniform) densities. Furthermore, any unimodal BR has both limits positive and equal.

On the contrary, the VIB distribution allows to have unimodal densities with one or both positive limits (see Figures 1a and 1c). Moreover, the VIB can also model the presence of extremely low (close to zero) and/or extremely high (i.e., close to one) values, maintaining a unimodal shape for the central data body (see Figure 1b), unlike both the beta and the BR.

Figure 1:

Some examples of VIB distributions

2.3 The variance inflated beta regression model VIBReg (

β, ϕ, k, p

)

The VIBReg is obtained by letting the ith response variable to be VIB distributed as $Y_{i} \sim VIB (μ_{i}, ϕ, k, p)$ ( $i = 1, \dots, n$ ), where, as in the previous models, the mean is given by

g (μ_{i}) = x_{i}^{T} β i = 1, \dots, n,

(2.7)

g being a suitable link function.

It is worthwhile to highlight the variance behaviour entailed by the model. By expression (2.5), $Var (Y_{i})$ is proportional to $μ_{i} (1 - μ_{i})$ , so that the VIBReg model induces a form of heteroscedasticity which is a very natural one for bounded random variables on $(0, 1)$ . Indeed, for such random variables, the variance is strictly dependent on the mean $μ_{i}$ , going to zero when $μ_{i}$ tends to zero or one. More precisely, for any given mean, the maximum possible variance is exactly equal to $μ_{i} (1 - μ_{i})$ , which is therefore a natural scaling factor.

The VIBReg model shows other relevant properties such as strong identifiability. The latter means that two elements of a parametric family are equal if and only if the corresponding parameters are identical. Typically, general mixture models (including mixtures of arbitrary beta distributions) do not possess this property, often displaying invariance under permutation of the components instead. This may have negative consequences both on the interpretation of parameters (and therefore on the elicitation of priors), and on computational grounds, where, for example, labelling problems may arise in the estimation process. Under completely general conditions, the VIBReg does not share these difficulties, as shown by the following proposition (see 7 for the proof).

Proposition 1 (Identifiability of the VIBReg model)

Consider a vector of independent response variables $Y_{i} \sim VIB (θ_{i})$ , $(i = 1, \dots, n)$ with $θ_{i} = (μ_{i}, ϕ, k, p)$ . Suppose that

g (μ_{i}) = x_{i}^{T} β

where $x_{i}$ $(i = 1, \dots, n)$ is a vector of covariates of dimension $q$ (with $q \leq n$ ), $β$ is the vector of corresponding regression parameters, and g is a strictly monotone function from $(0, 1)$ to $ℝ$ . Let us denote the corresponding regression model by $VIBReg (η)$ , where $η = (β, ϕ, k, p)$ . Let $Y \sim VIBReg (η)$ and $Y^{'} \sim V I B Re g (η^{'})$ . Then, if the design matrix $X$ is of full rank, $Y \sim Y^{'}$ if and only if $η = η^{'}$ .

Another typical issue within mixture models is likelihood unboundedness. Possible unbounded likelihood points may generate computational difficulties and are generally of no inferential relevance, being often caused by an overfitting problem which can occur, for example, when a component of the mixture is entirely dedicated to a single observation. The VIBReg likelihood is not bounded, but very simple and mild restrictions on its parameter space can be adopted to avoid unbounded points.

Proposition 2 (Boundedness of the VIBReg model likelihood) Consider a vector of independent response variables $Y_{i} \sim VIB (θ_{i})$ , $(i = 1, \dots, n)$ with $θ_{i} = (μ_{i}, ϕ, k, p)$ and $g (μ_{i}) = x_{i}^{T} β$ , with $β$ of dimension $k + 1$ . Suppose that $n > k + 1$ . Then, the likelihood $L (β, ϕ, k, p | y)$ is a.s. bounded under each of the following constraints on the parametric space: (i) $ϕ \leq M$ for $M > 0$ , (ii) $k \geq ε$ for $0 < ε < 1$ , (iii) $k \geq ε ϕ^{- α}$ when $ϕ \geq M$ , where $0 < ε < M^{α}$ and $0 \leq α < 1$ .

In Appendix B, we give the complete proof of Proposition 2. Note that the constrained parametric spaces implied by (i) and (ii) can be viewed as restrictions of the one corresponding to case (iii), which is the more general. These constrained spaces have the following implications: case (i) implies a lower bound on the variance of the two mixture components, and a bound on their ratio, case (ii) implies a bound on the variance ratio only, case (iii) implies a bound on neither of them. Indeed, it is easy to see that, by choosing sufficiently low values of $ε$ and high values of $α$ and $M$ in case (iii), virtually the whole original parametric space can be spanned.

The two aforementioned propositions apply, in particular, to an i.i.d. sample from the VIB distribution, being this a special case of the VIBReg model, namely the case with a single (fictitious) constant covariate.

For comparison, recall that the FBReg model is both identifiable and with bounded likelihood, whereas the BRReg possesses neither property (see (Migliorati et al., 2018)).

3 Estimation issues

The likelihood of the VIBReg model for independent observations $y^{T} = {y_{1}, \dots, y_{n}}$ is given by:

L (η | y) = \prod_{i = 1}^{n} f_{VIB} (y_{i} | μ_{i}, ϕ, k, p),

(3.1)

where $η = {β, ϕ, k, p}$ , the overall mean $μ_{i} = g^{- 1} (x_{i}^{T} β)$ is a function of covariates and $f_{VIB} (y_{i} | μ_{i}, ϕ, k, p)$ is given by (2.2). Note that no explicit solution to the (incomplete data) estimation problem exists, since the allocation of each ith observation to one of the two mixture components is unknown.

We decided to adopt a Bayesian approach, which makes it easy to cope with complex models with many parameters and, moreover, with small sample problems affecting maximum likelihood estimates. In order to apply Markov chain Monte Carlo (MCMC) techniques, a recent solution is the adaptive HMC algorithm (Duane et al., 1987). The latter is a generalization of the Metropolis algorithm, which combines MCMC and deterministic simulation methods. This is achieved through augmentation of the posterior density by an independent distribution of an auxiliary variable, introduced to enable faster movements of the algorithm through the parameter space. Typically, this leads to greater efficiency with respect to Gibbs sampler or Metropolis techniques, by strongly reducing autocorrelation among draws (Brooks et al., 2011). One of the main drawbacks of HMC samplers, when applied to mixture models, is that they are only defined on continuous distributions. Indeed, in any mixture model it is necessary to introduce n-dimensional vectors of latent indicators for mixture components. Therefore, a partition of the parameter space must be introduced, and HMC updates for continuous parameters and Metropolis updates for discrete ones must be alternated. Finally, the discrete parameters can be summed out of the model.

Clearly, a complete specification of the Bayesian model requires the elicitation of the parameters’ prior distributions. We assumed a priori independence (a usual choice when no prior information is available) and adopted flat priors. More precisely, for the regression parameters β, we selected a multivariate normal prior with zero mean vector and a diagonal covariance matrix with ‘large’ values for the variances. For the remaining parameters, we chose a gamma distribution $Ga (g, g)$ for ϕ and uniform priors for the remaining parameters $k \sim U (0, 1)$ and $p \sim U (0, 1)$ . A numerical investigation showed that there is robustness (limited impact on inferential conclusions) with respect to different choices of the hyperparameters (i.e., prior variances of the normal distribution and g).

We implemented the HMC algorithm resorting to the Stan program (Stan Development Team, 2016; Gelman et al., 2014) and to its default No-U-Turn Sampler (NUTS) (Hoffman and Gelman, 2014).

4 Robustness evaluation

In this Section, we perform an extensive simulation study to assess the robustness of the proposed VIBReg model in presence of various configurations and various degrees of intensity of outliers and to compare it with competing models. Robustness will be measured through a broad range of tools taking different perspectives into account. In particular, the impact of outliers shall be assessed by considering marginal and global posteriors of regression parameters, measures of distance between regression curves, fitting measures, as well as predictive intervals and residuals.

To contaminate a subset of response values with the purpose of generating artificial outliers, we apply the perturbation operation on the simplex, which can be viewed as the analogous of addition onto real space (Pawlowsky-Glahn et al., 2015). Indeed, simply shifting the observed responses by a given range of quantities may generate values outside the interval $(0, 1)$ .

In general terms, the perturbation of vector $z$ by perturbing vector $t$ , both defined on a $D -$ dimensional simplex $𝒮^{D} = \{w \in ℝ^{D} : w_{i} > 0, i = 1, \dots, D,$ $\sum_{i = 1}^{D} w_{i} = 1\}$ , is equal to:

z \oplus t = C {z_{1} t_{1}, \dots, z_{D} t_{D}},

(4.1)

where 𝒞 is the closure operation being:

C {q} = \{\frac{q_{1}}{\sum_{i = 1}^{D} q_{i}}, \dots, \frac{q_{D}}{\sum_{i = 1}^{D} q_{i}}\},

where $q \in R_{+}^{D}$ , that is, $q_{i} > 0, \forall i = 1, \dots, D$ . The vector resulting from perturbation operation in (4.1) lies on $𝒮^{D}$ .

In a 2-dimensional simplex, the perturbing vector can be set equal to $(δ, 1 - δ)$ for $0 \leq δ \leq 1$ . In this way, as $δ$ increases from zero to one, any given response value $y_{i}$ will span the whole interval $[0, 1]$ . The neutral element of the perturbation, that is, the one which leaves the perturbed value unchanged, is $δ = 0.5$ , whereas a value greater (lower) than 0.5 determines an increase (decrease) of the perturbed value.

4.1 Outlier scenarios

The initial simulation framework considers a BReg model with a quantitative covariate x. Specifically, we simulate $n = 100$ i.i.d. observations from $X \sim Unif (0, 1)$ and we define $logit (μ_{i}) = β_{0} + β_{1} x_{i}$ , for $i = 1, \dots, n$ having fixed $β_{0} = 0.5$ and $β_{1} = 3$ . Then, we simulate a baseline dataset of n independent values from $Y_{i} \sim Beta (μ_{i}, ϕ)$ for several values of ϕ and replicate it $N = 500$ times. For space constraints, we report results only for $ϕ = 30$ , which is a representative case. We also performed our analysis with a different sample size, namely $n = 50$ . As we obtained fairly similar results, in the following we present only the $n = 100$ case.

Each baseline dataset has been further contaminated according to three patterns of perturbation as follows:

Increase by a perturbing factor $δ$ 4 response values associated with low x values;

Increase by a perturbing factor $δ$ 4 response values associated with spread x values;

Increase by a perturbing factor $δ$ 2 values associated with low x values and decrease by $(1 - δ)$ 2 values associated with high x values.

For each pattern of perturbation, different degrees of intensity have been evaluated by letting $δ$ assume values 0.65,0.75, 0.85, and 0.95.

Therefore, all the graphs reported further, which help to illustrate the various effects of outliers, have $δ$ on the horizontal axis, where $δ = 0.5$ corresponds to the uncontaminated (baseline) case.

The scatter plots for $δ = 0.75$ and $δ = 0.85$ for the first replication in each scenario are reported in Figure 2.

Figure 2:

Note: Filled-in circles for outliers.

The simulated scenarios have been defined following a long-established scheme (beginning with Peña and Guttman [1993], amongst other) of contamination of a baseline case. Here, the scenarios we devised are comparable to the ones in Bayes et al. (2012), the main differences being the introduction of the perturbation operation which replaces the addition and the removal of their Scenario II (i.e., decrease responses for high x values), which is symmetric to our Scenario I and gives similar results.

Under each of the four regression models, the parameters have been estimated by HMC through the Stan software, as described in Section 3, and the convergence to the equilibrium distribution has been checked by means of diagnostic tests for stationarity (Geweke and Heidel diagnostics) and autocorrelation (Raftery diagnostic) (Mengersen et al., 1999).

4.2 Marginal posteriors

A first evaluation of the robustness of the four models can be obtained by analysing the variations induced by the outliers in the posterior distribution of each regression parameter. To this end, we plot the MC posterior mean and 95% credible interval (CI) as a function of $δ$ , separately for $β_{0}$ and $β_{1}$ , for each regression model, see Figures 3, 4 and 5.

Figure 3:

Scenario I: Posterior mean (solid line) and 95% credible interval (dashed line) of $β_{0}$ (top panel) and $β_{1}$ (bottom panel) for each regression model as $δ$ varies

In all scenarios, the regression parameters of the BReg seem to be highly affected by increasing values of $δ$ revealing a remarkable lack of robustness. Indeed, in Scenarios I and III, the 95% CIs associated with $β_{1}$ do not contain the true value for $δ = 0.95$ . A similar behaviour, albeit milder, is observed for the FBReg which, by dedicating one of its components to model the group of outlying observations, shifts the overall mean towards them. In these scenarios (namely I and III), the VIBReg model provides the most robust behaviour, and its superiority with respect to the BRReg one is mostly evident with respect to $β_{0}$ .

Furthermore, it should be noted that under the BRReg model, the intercept tends to move downwards (below the true value) even when only upper outliers are present.

This apparently odd (and in contrast with the other models) performance can be explained by recalling that the BRReg, by modelling the outliers with a uniform, assigns a 0.5 mean to them, which has the effect of always shrinking the overall mean towards such a value regardless the data pattern.

In Scenario II, the VIBReg and the BRReg models behave very well even for high degrees of perturbation.

Figure 4:

Scenario II: Posterior mean (solid line) and 95% credible interval (dashed line) of $β_{0}$ (top panel) and $β_{1}$ (bottom panel) for each regression model as $δ$ varies

Figure 5:

Scenario III: Posterior mean (solid line) and 95% credible interval (dashed line) of $β_{0}$ (top panel) and $β_{1}$ (bottom panel) for each regression model as $δ$ varies

4.3 Global regression impact

A second assessment of the outliers impact is measured in terms of modification of the joint, instead of the marginal, parameter posteriors. In particular, we focus on the regression coefficients, being typically the most relevant parameters of the model. We, therefore, compute the Kullback–Leibler (KL) divergence between the joint posterior referred to the unperturbed dataset, $π (β | y)$ , and to the perturbed dataset depending on $δ$ , $π^{δ} (β | y)$ . From the MCMC values, we estimate the posterior densities via kernel approach (Simonoff, 2012; Duong et al., 2007). Then the computation of the KL divergence is straightforward, having checked that the area under each estimated density is (nearly) equal to 1.

Figure 6:

KL divergence measure for VIBReg (solid line), FBReg (dashed line), BRReg (dotted line) and BReg (dotdash line)

The inspection of the KL divergence over the three simulated scenarios (Figure 6) confirms the better general robustness of the VIBReg model, whose distance curves either lie below or are comparable to the others. The VIBReg is the only model with non-(rapidly) increasing curves in all three scenarios. Interestingly, the proposed model is remarkably stable for extreme levels of perturbation, in all scenarios, unlike the other models.

A careful investigation of the plots in Section 4.2 reveals that in some cases, $β_{0}$ and $β_{1}$ move towards opposite directions when outliers are introduced. This phenomenon may produce a balancing effect on the resulting regression curve, which may result only slightly altered. A clear example of this possibility is shown in the application of Section 5.1.

Therefore, it seems relevant to directly evaluate the global impact of the outliers on regression curves. To this end, we consider the $L_{1}$ distance between unperturbed and perturbed fitted regression curves, given by:

\int |g^{- 1}({\bar{β}}_{0}^{δ} + {\bar{β}}_{1}^{δ} x)- g^{- 1}({\bar{β}}_{0} + {\bar{β}}_{1} x)| d x

(4.2)

where the integral is over the range of observed x-values, $g (\cdot)$ is the logit function, while $({\bar{β}}_{0}, {\bar{β}}_{1})$ and $({\bar{β}}_{0}^{δ}, {\bar{β}}_{1}^{δ})$ are the posterior means referred to the unperturbed dataset and, respectively, to the dataset perturbed by $δ$ .

Figure 7:

Curve distances for VIBReg (solid line), FBReg (dashed line), BRReg (dotted line) and BReg (dotdash line)

The results reported in Figure 7 are largely consistent with the ones relative to the joint posterior KL distance, with an even more evident better VIB performance.

The main difference lies in a general worse behaviour of the BRReg model, especially in Scenario III where the BRReg is outperformed even by the BReg for moderate perturbation. This may be due to the balancing effect of the estimates of the regression parameters which is weaker for the BRReg.

4.4 Fitting measures

The aim of this Section is to consider the effects of outliers on the fitting ability of the models. A well-established way of comparing models is through the value of proper fitting criteria, which typically mediate between goodness of fit and complexity. Different measures have been proposed within a Bayesian framework, the most popular being the deviance information criterion (DIC) (Spiegelhalter et al., 2002) and the Bayesian counterparts of the well-known Akaike information criterion (AIC) (Akaike, 1998) and Bayesian information criterion (BIC) (Schwarz, 1978) criteria, namely expected AIC (EAIC) and expected BIC (EBIC) (Brooks, 2002). DIC is defined as:

DIC = \bar{D} + p_{D} = 2 \bar{D} - D (\bar{η})

where $D (η) = - 2 \log [L (η | y)]$ is the model deviance, $\bar{D}$ is the posterior mean of the deviance, $\bar{η}$ is the vector of posterior means of the parameters, and $p_{D} = \bar{D} - D (\bar{η})$ can be interpreted as an estimate of the ‘true’ number of parameters. Indeed, the latter suggests how many parameters one loses by adopting a Bayesian approach (Brooks, 2002) due to the prior-induced parameter reduction.

The other two criteria take the form:

EAIC = \bar{D} + 2 p, EBIC = \bar{D} + p \log (n)

where $p$ is the number of parameters and n is the sample size.

These criteria are known to have some problems, especially in mixture models contexts. In particular, DIC is not defined for non-regular statistical models and can produce negative estimates of $p_{D}$ (typically for non-identifiable models). EAIC and EBIC use a penalization based on $p$ , ignoring the ‘true’ parameter number issue.

Recently, a new criterion called WAIC has been proposed (Watanabe, 2010). This is a fully Bayesian measure in that it uses the entire posterior distribution (instead of point estimates only), the main consequences being invariance to reparameterizations and possibility of computation also for non-regular models (Gelman et al., 2014; Vehtari et al., 2017).

WAIC evaluates the goodness of fit of a model through the log pointwise predictive density ( $lppd$ ). Such quantity is evaluated using draws $η^{s}$ ( $s = 1, \dots, S$ ) from the simulated posterior for a value of $S$ large enough to fully capture the entire posterior, namely:

\hat{lppd} = \sum_{i = 1}^{n} \log (\frac{1}{S} \sum_{s = 1}^{S} f (y_{i} | η^{s}))

where $f (y_{i} | η)$ is the density function of the ith observation.

The model complexity adjustment is essential here, as the $lppd$ , referred to observed data, overestimates the corresponding quantity expected for future data. Therefore, a correction for effective number of parameters allows to adjust for overfitting. This can be done in two ways, giving rise to two variants of the criterion.

The first penalization follows a logic similar to that used for DIC and it can be estimated by:

{\hat{p}}_{WAIC 1} = 2 \sum_{i = 1}^{n} \{\log (\frac{1}{S} \sum_{s = 1}^{S} f (y_{i} | η^{s})) - \frac{1}{S} \sum_{s = 1}^{S} \log (f (y_{i} | η^{s}))\} .

(4.3)

The second is based on the posterior variance of the log predictive density, which can be estimated as:

{\hat{p}}_{WAIC 2} = 2 \sum_{i = 1}^{n} V_{s = 1}^{S} (f (y_{i} | η^{s}))

where $V_{s = 1}^{S} a_{s} = \frac{1}{S - 1} \sum_{s = 1}^{S} (a_{s} - \bar{a})^{2}$ is the sample variance.

In both cases WAIC is given by −2 times the difference between $\hat{lppd}$ and the chosen measure of model complexity. Clearly, the smaller the values of the comparison criteria, the better the model.

In our simulation study, both WAIC criteria have been computed, but only the first is reported, say WAIC₁ (see penalization 4.3), since they provide completely consistent conclusions. In particular, Figure 8 shows the WAIC₁ values as δ varies.

Figure 8:

Mean values of WAIC₁ criterion for VIBReg (solid line), FBReg (dashed line), BRReg (dotted line) and BReg (dotdash line) models as $δ$ varies in the three scenarios

As a general pattern, for low values of $δ$ , corresponding to weakly perturbed datasets, the VIBReg, FBReg and BReg models show similar goodness of fit, while the BRReg one shows a slightly worse fit.

For higher $δ$ values, BReg is the worst model in all scenarios. In Scenarios I and II, the FBReg model provides the best fit. On the contrary, in the last scenario, VIBReg and BRReg display the best performance. VIBReg always displays a better (or similar) behaviour than the BRReg, behaviour which seems in general sufficiently satisfactory even in the first two scenarios.

4.5 Residuals and predictive intervals

To obtain further insight into the performances of the four models in the various outlier scenarios, we analyse residuals and predictive intervals based on a given sample with $n = 100$ . Specifically, for each simulated scenario, we take one sample and fix a perturbing value of $δ$ so as to guarantee at least three visible outlying observations among the four perturbed responses. In particular, $δ = 0.75$ has been chosen for Scenarios I and II, while $δ = 0.85$ for Scenario III (see the scatter plots in Figure 2).

We examine Bayesian residuals (Gelman et al., 2014), defined as:

r_{i} = \frac{y_{i} - {\hat{μ}}_{i}}{\sqrt{\hat{Var} (y_{i})}}, i = 1, \dots, n

(4.4)

where $y_{i}$ is the observed response variable, ${\hat{μ}}_{i} = g^{- 1} ({\hat{β}}_{0} + {\hat{β}}_{1} x_{i 1})$ is the predicted mean and $\hat{Var} (y_{i})$ is the variance of $y_{i}$ , both evaluated at the posterior means of the parameters.

Figure 9:

Note: Filled-in circles for outliers.

BReg, FBReg and VIBReg exhibit a similar outlier detection capability. Roughly speaking, they recognize third out of the fourth introduced outliers in the first two scenarios, and all of them in the third one (see Figure 9). In addition, all these models identify as an anomalous response the non-perturbed value with a high negative residual. This seems a reasonable outcome in the light of Figure 2. Note that the anomalous non-perturbed value is the rightmost among the lowest ones, occurring where the predictive mean is close to one and, consequently, the variance is very low. On the contrary, BRReg only detects the two most extreme lower outliers introduced in Scenario III. This is because the uniform component of the BR, which models outliers, induces a particularly high response variability. Indeed, the estimated response variance under the BRReg is always (in each scenario and for any covariate value) at least three times as large as the variance under the other three models. See Figure 10, where the ratios of the BRReg, VIBReg and BReg variances to the FBReg one are reported. From this figure, it also emerges that the other models display a similar variance.

Figure 10:

Estimated response variance ratio with respect to FBReg in the three scenarios. VIBReg (solid line), BRReg (dotted line) and BReg (dotdash line)

Finally, we simulate draws from the posterior predictive distribution for each model (see Figures 16, 17 and 18 in 16). The FBReg provides the smallest and best centred intervals in the first two scenarios, whereas the BReg always has the poorest performance. The BRReg almost always displays very large intervals for the outliers (often nearly covering the whole (0,1) interval), thus being poorly informative. Again, this is caused by its uniform component. The VIBReg shows a slightly worse performance than the FBReg in the first two scenarios and a comparable behaviour in the third one.

4.6 Overall simulation findings analysis

It is useful to perform an overall comparison of the regression models under the various scenarios and different perspectives considered in the simulation in order to obtain useful indications on which model is preferable in different contexts.

Let us first compare the models in terms of residuals, as here results are quite clear-cut and somewhat special. Indeed, the BRReg is the only model unable to identify nearly all the outliers in all scenarios. This is because the uniform component of such model greatly increases the response variability, thereby producing a masking effect of the outliers. On the contrary, the VIBReg and FBReg manage to accurately identify anomalous observations, and exhibit a considerable robustness at the same time.

Let us now compare the models excluding the effects on residuals.

In all scenarios, the VIBReg has better or comparable performances than all other models in terms of robustness of marginal and joint posteriors as well as regression curves. Very importantly, it is the only model always showing a great stability when extreme outliers are present. The FBReg prevails in terms of predictive intervals and fitting ability in the first two scenarios. The VIBReg is the best performing model in Scenario III with respect to all considered effects. In particular, it displays a substantial larger stability than the other models in terms of regression coefficients and curve.

In no occasion, the BRReg dominates the other models (nor the BReg) resulting particularly unsatisfactory with respect to stability of the regression curve, residuals and predictive intervals.

5 Application to Australian Institute of Sport data

With the purpose of comparing the four competing models, we provide an application to the Australian Institute of sport dataset, which is included in the R library sn. This dataset has been analysed by Bayes et al. (2012) to show BRReg robustness. It consists of 37 rowing athletes for whom the body fat percentage (Bfat, response variable), the lean body mass (lbm, a quantitative covariate) and gender (encoded as 1 if female and 0 if male) have been collected.

5.1 Regression model with one covariate

In all four models, a regression model for the mean is defined as follows:

logit (μ_{i}) = β_{0} + β_{1} X_{i}, for i = 1, \dots, n,

(5.1)

where $X_{i}$ is the covariate lbm.

Here, the presence of two lower outliers (filled-in circles corresponding to subjects 16 $th$ and 30 $th$ in Figure 11) is well-known and evident, which is a case a priori favourable to the FBReg, as confirmed by simulations.

The estimation of vector of parameters, both on the entire dataset and after removing the outliers, has been dealt with through HMC algorithm (see Section 3). Table 1 reports the parameter posterior means as well as the fitting measure WAIC₁ (see Section 4.4).

Table 1:

Posterior parameter means and fitting measure WAIC₁ for the four models under the regression model with one covariate

		WAIC₁	$β_{0}$	$β_{1}$	ϕ	$p$	$w$	$α$	$k$
All observations	VIBReg	−144.8	0.631	−0.035	231	0.219			0.157
	FBReg	−153.2	0.738	−0.037	223	0.917	0.498
	BRReg	−147.2	0.902	−0.036	228			0.193
	BReg	−135.5	0.098	−0.027	91
No outliers	VIBReg	−161.3	0.839	−0.038	303	0.483			0.655
	FBReg	−170.7	0.895	−0.039	860	0.449	0.122
	BRReg	−160.2	0.861	−0.037	229			0.074
	BReg	−161.9	0.839	−0.038	229

The best fit is provided by the FBReg model both on the entire dataset and after removing the outliers. The regression coefficient estimates are similar among models in the case with no outliers. In terms of single regression coefficients, the BReg results the most affected by outliers, whereas the BRReg is the most stable. However, the overall impact of outliers on the regression curves leads to different conclusions. To see this, in Figure 11, we report the fitted curves for the four models. As they are virtually overlapping when outliers are absent, in this case, we plot only one curve to preserve readability.

Figure 11:

Note: The thick solid line represents the common fitted regression curve after removing outliers.

Comparing the curves with and without outliers, one can see that the VIBReg and the FBReg exhibit quite a smaller variation than the other two models. Indeed, this is confirmed by the curve distances (see formula 4.2), which take the values 0.264 (VIBReg), 0.406 (FBReg), 0.709 (BRReg) and 0.906 (BReg). Interestingly, the regression curve of the BRReg is shifted upwards, that is, in the direction opposite to the outliers, unlike the other models. This behaviour, as already noticed, is due to the uniform component which drives the overall regression curve towards 0.5.

An analysis of the residuals, of the estimated response variances, and of the predictive intervals (see Figures 12, 13 and 14) confirms the conclusions drawn in Section 4.5.

Figure 12:

Standardized residuals for each subject and for each model on the entire dataset

Figure 13:

Estimated response variance ratio with respect to FBReg. VIBReg (solid line), BRReg (dotted line) and BReg (dotdash line)

Figure 14:

Posterior predictive means (cross symbols) and 90% predictive intervals for each subject and for each model (observed values in circles, filled-in for outliers) on the entire dataset

Specifically, only the BRReg model residuals do not identify the two outliers. This is because its response variance is at least 20 times larger than the variances of the other models (see Figure 13). In terms of predictive intervals, only the BReg fails to account for the outliers. Moreover, the FBReg shows the smallest length of the intervals corresponding to outliers, while this length is larger for the VIBReg (approximately equal to 0.5) and especially for the BRReg (approximately equal to 1).

5.2 Regression model with three covariates

The analysis on the Sport dataset has been further enhanced by evaluating the impact of gender on the response variable. To this purpose, the regression function for the mean, in all four competing models, is set as follows:

logit (μ_{i}) = β_{0} + β_{1} X_{1 i} + β_{2} X_{2 i} + β_{3} X_{3 i} i = 1, \dots, n

where $X_{1 i}$ is the covariate $lbm$ , $X_{2 i}$ is the covariate about gender encoded as 1 if female and 0 if male, and $X_{3 i}$ is the interaction between these two covariates.

To better explore the flexibility provided by our new model, we additionally regressed onto covariates the parameter $k$ of the VIBReg model as follows:

logit (k_{i}) = η_{0} + η_{1} X_{2 i} .

(5.2)

This is consistent with data pattern which displays a different shape for male and female athletes (see top left panel of Figure 15). When accounting for gender, a third outlier clearly shows up belonging to the female group (lowest filled squared point).

Furthermore, one of the two previously identified outliers (squared empty point in Figure 15) in the new setting becomes a very extreme one. For this reason, we decided to conduct our analysis by considering the following three cases: (1) entire dataset, (2) dataset without the extreme male outlier, (3) dataset without all three outliers.

Figure 15:

Note: Filled points refer to female, squared points refer to outliers.

As in the one covariate example, we report the posterior means of the parameters as well as the fitting measure WAIC₁ of the various models (see Table 2). It is worth noting that the BRReg model always provides the worst fit, whereas the VIBReg outperforms all the other models when outliers are present.

Table 2:

Parameter estimates and WAIC₁ under the regression models with three covariates

Case		WAIC₁	$β_{0}$	$β_{1}$	$β_{2}$	$β_{3}$	ϕ	$p$	$w$	$α$	$η_{0}$	$η_{1}$
(1)	VIBReg	−162.1	−1.604	−0.008	0.565	0.002	382.570	0.503			0.167	−0.468
	FBReg	−160.9	−1.556	−0.009	0.549	0.002	205.772	0.698	0.222
	BRReg	−156.9	−1.441	−0.009	0.492	0.003	185.866			0.085
	BReg	−160.7	−1.594	−0.008	0.537	0.002	198.606
(2)	VIBReg	−159.5	0.540	−0.035	−1.572	0.029	476.527	0.543			0.102	−0.665
	FBReg	−157.7	0.473	−0.034	−1.467	0.027	212.989	0.684	0.216
	BRReg	−155.2	0.588	−0.034	−1.540	0.028	191.870			0.086
	BReg	−157.1	0.550	−0.035	−1.585	0.029	203.860
(3)	VIBReg	−165.4	0.559	−0.035	−1.070	0.021	467.640	0.493			0.469	−0.093
	FBReg	−165.4	0.532	−0.035	−1.020	0.020	306.838	0.594	0.200
	BRReg	−163.7	0.573	−0.034	−1.013	0.020	265.703			0.075
	BReg	−165.7	0.489	−0.034	−0.982	0.019	315.159

A graphical representation of the fitted regression curves for the four models can be found in Figure 15. The VIBReg, FBReg and BReg regression curves are basically overlapping in all scenarios. Furthermore, no model seems to be stable when the extreme male outlier is introduced (case 1). Indeed, this extreme point causes a substantial flattening effect for the male group and is never recognized as an outlier. However, fairly satisfactory results are obtained when such a point is excluded. It should also be noted that the regression curve of the BRReg model is still affected by the shifting upwards phenomenon, already observed in one covariate example. This may explain its poor fit.

Thus, although the VIBReg is not more stable than the other models in this particular example, its special structure does allow for a better fitting in presence of outliers.

It appears also of interest to examine how the additional flexibility, introduced by letting the parameter $k$ vary with gender, is exploited by the VIBReg model. In Table 3, we computed the posterior estimates of inflation parameter $k$ separately for males and females together with the estimates of precision parameters. When no outliers are included, case (3), the effect of shrinkage on parameter ϕ is similar across groups. Conversely, a relevant difference in terms of variance inflation is achieved when outliers are present, specially in case (2) where the extreme male outlier is excluded. This suggests that it may be convenient, in some more complicated contexts, to model the inflation parameter $k$ as a function of appropriate covariates to achieve a better data description.

Table 3:

Posterior estimates of inflation parameter $k$ and precision parameters separately for males and females

Case	ϕ	$k_{M}$	$k_{F}$	$k_{M} ϕ$	$k_{F} ϕ$
(1)	382.570	0.542	0.425	207.353	162.592
(2)	476.527	0.525	0.363	250.177	172.979
(3)	467.640	0.615	0.593	287.599	277.311

Even in this case, we have computed residuals and posterior predictive intervals. We decided not to report them, because they display the same characteristics already highlighted in the one covariate example.

6 Discussion and final remarks

In this article, we introduced a new regression model for proportions aimed at handling outlier observations, and we compared it with the BReg model and two other mixture regression models for bounded responses. All three mixture models share a common beta mixture structure with two components: an arbitrary beta component is intended to account for the main body of data and the other one is devoted to capture outlying observations. This second component is expected to better account for the increased variability and to make the parameter estimates and regression curves less affected by outliers. In particular, the second component is taken to be a uniform distribution in the BRReg, a beta distribution with the same variance as the main component but a shifted mean in the FBReg, and a beta distribution with the same mean but an inflated variance in the VIBReg model.

Indeed, the results of our robustness study show that, in many respects, mixture models are by far more stable than the (BReg) one.

More specifically, in all scenarios, the VIBReg has better or comparable performances than all other models in terms of robustness of estimates as well as regression curves, with specially outstanding results when extreme outliers or both upper and lower outliers are present. The FBReg prevails in terms of predictive intervals and fitting ability in the first two scenarios, with the VIBReg being just slightly worse. It follows that the VIBReg, unlike the other models, appears a satisfactory choice under all examined conditions.

From a computational perspective, the use of Stan software, instead of the more accessible to practitioners WinBUGS, made the application of the estimation algorithms quite efficient in terms both of runtimes and of convergence characteristics. This is particularly relevant, since, within mixture models, low convergence rates or even convergence failures often occur.

Although in the article, we have considered response variables with values in $(0, 1)$ , the VIBReg model (as well as all the other discussed models) can be straightforwardly extended to deal with variables defined on an arbitrary bounded interval. It is sufficient to apply a suitable linear transformation, which does not alter properties and interpretation of the model.

Finally, another relevant extension of the VIBReg can be achieved by letting some other model parameters (besides the mean) to depend on covariates. An example in this direction has been given in the application of Section 5.2, where the variance inflation parameter is considered. In addition, although the model implies a natural form of heteroscedasticity, as discussed in Section 2.3, it may be useful in some settings to model the variability (i.e., the parameter ϕ) as a function of covariates as shown, for example, by Smithson and Verkuilen (2006). This can be achieved within the VIBReg without special difficulties.

Supplementary materials

Supplementary materials for this article are available at http://www.statmod.org/smij/archive.html.

Acknowledgements

We are grateful to the referees and to the editor for their constructive comments, which helped to improve the article.

Appendices

Appendix A Proof of Proposition 1 (Identifiability of the VIBReg model)

Proof. First we show that the VIB distribution is identifiable. Let $Y \sim VIB (θ)$ and $Y^{'} \sim VIB (θ^{'})$ . Then, we have to show that $Y \sim Y^{'}$ if and only if $θ = θ^{'}$ . It is sufficient to prove that $Y \sim Y^{'}$ implies $θ = θ^{'}$ , the converse being obvious. Clearly $E (Y) = E (Y^{'})$ and therefore $μ = μ^{'}$ . $Y$ has density

\begin{matrix} f (y; θ) & = x^{μ ϕ k - 1} (1 - x)^{(1 - μ) ϕ k - 1} [p c (ϕ k, μ) \\ + (1 - p) c (ϕ, μ) x^{μ ϕ (1 - k)} (1 - x)^{(1 - μ) ϕ (1 - k)}] \end{matrix}

where

c (λ, δ) = \frac{Γ (λ)}{Γ (δ λ) Γ ((1 - δ) λ)},

and it is a continuous function on $(0, 1)$ . As $f (y; θ) = f (y; θ^{'})$ a.s., continuity implies that the two densities must be identically equal for any $y \in (0, 1)$ . Consider now the behaviour of $f$ at the origin. It is not difficult to check that, when $y \to 0$ , $f (y; θ) y^{α}$ converges to a strictly positive finite constant if and only if $α = 1 - ϕ μ k$ . In this case, the limit as $y$ approaches zero is $p c (ϕ k, μ)$ .

It follows that $f (y; θ^{'}) y^{1 - ϕ μ k}$ must converge to a strictly positive finite constant when $y \to 0$ , which is possible only if $ϕ μ k = ϕ^{'} μ^{'} k^{'}$ . Hence, $ϕ k = ϕ^{'} k^{'}$ . Furthermore, the two limits must coincide, that is, $p c (ϕ k, μ) = p^{'} c (ϕ^{'} k^{'}, μ^{'})$ This, given the previous results, implies that $p = p^{'}$ . By using the equality $Var (Y) = Var (Y^{'})$ , where

Var (Y) = μ (1 - μ) (\frac{p}{ϕ k + 1} + \frac{1 - p}{ϕ + 1}),

we obtain $ϕ = ϕ^{'}$ . As we have shown that $ϕ k = ϕ^{'} k^{'}$ , this implies that $k = k^{'}$ , thereby completing the proof of the VIB identifiability.

Let us now focus on the VIBReg identifiability. Obviously, $η = η^{'}$ implies that $Y \sim Y^{'}$ . Suppose now that $Y \sim Y^{'}$ . Then $Y \sim Y^{'}$ , $i = 1, \dots, n$ . By the just shown identifiability of the VIB distribution, we have that $p = p^{'}$ , $ϕ = ϕ^{'}$ , $k = k^{'}$ , and $g^{- 1} (x_{i}^{T} β) = g^{- 1} (x_{i}^{T} β^{'})$ , $i = 1, \dots, n$ . It follows that $X^{T} X β {= X}^{T} X β^{'}$ and, therefore, $β = β^{'}$ , being $X$ of full rank $q$ . $□$

Appendix B Proof of Proposition 2 (Boundedness of the VIBReg model likelihood)

Proof. Let us first study some asymptotic properties of the log-likelihood $l (μ, ϕ | y) = \log f_{B}^{*} (y; μ, ϕ))$ of a single beta distributed observation $y$ , written in the mean-precision parametrization (2.3). This log-likelihood can be written as

l (μ, ϕ | y) = \log Γ (ϕ) - \log Γ (ϕ μ) - \log Γ (ϕ (1 - μ)) + ϕ (μ \log y + (1 - μ) \log (1 - y)) .

(B.1)

By Lemma 1 in the Appendix of Migliorati et al. (2017), a necessary condition for $l (μ, ϕ | y)$ to diverge to $+ \infty$ is that both $ϕ μ$ and $ϕ (1 - μ)$ tend to $+ \infty$ . In such a case, an expansion of the relevant terms in (B.1) produces the following asymptotic approximation:

l (μ, ϕ | y) ≃ - ϕ K (y, μ) + (1 / 2) \log ϕ,

(B.2)

where

K (y, μ) = μ \log \frac{μ}{y} + (1 - μ) \log \frac{1 - μ}{1 - y}

is the KL divergence between the vector of probabilities $(y, 1 - y)$ and $(μ, 1 - μ)$ . Note that $K (y, μ) \geq 0$ , equality holding if and only if $y = μ$ . Approximation (B.2) is based on the following asymptotic expansion valid for $x \to + \infty$ :

\log Γ (x) = (x - 1 / 2) \log x - x + (1 / 2) \log π + O (\frac{1}{x}) .

Let us now consider the log-likelihood $l (β, ϕ, k, p)$ for the VIBReg model:

\begin{matrix} l (β, ϕ, k, p) & = \sum_{i = 1}^{n} \log [p f_{B}^{*} (y_{i}; μ_{i}, ϕ k) + (1 - p) f_{B}^{*} (y_{i}; μ_{i}, ϕ)] \\ = \sum_{i = 1}^{n} \log f_{VIB} (y_{i}; μ_{i}, ϕ, k, p), \end{matrix}

(B.3)

where $μ_{i} = g^{- 1} (x_{i}^{T} β)$ . It is easy to see that the parametric subspaces corresponding to cases (i) and (ii) are contained in the one defined by (iii) for appropriate choices of the constants. In particular, case (ii) is obtained by choosing $α = 0$ .

Therefore, it is enough to prove likelihood boundedness for case (iii). A necessary condition for $l (β, ϕ, k, p)$ to diverge is that at least one of the n terms in (B.3), say the $i th$ , tends to $+ \infty$ . By previous results on beta likelihood, this can happen only if both $ϕ μ_{i} \to + \infty$ and $ϕ (1 - μ_{i}) \to + \infty$ , which implies that $ϕ \to + \infty$ . In this case, under condition (iii), $ϕ k$ diverges as well. If either $ϕ k μ_{i}$ or $ϕ k (1 - μ_{i})$ does not diverge (implying that $μ_{i}$ tends to zero or one) then $f_{B}^{*} (y_{i}; μ_{i}, ϕ k)$ is bounded. Furthermore, by (B.2), $f_{B}^{*} (y_{i}; μ_{i}, ϕ)$ behaves asymptotically as

p e^{- ϕ K (y_{i}, μ_{i})} \sqrt{ϕ},

thus converging to zero, as $μ_{i}$ tends to zero or one and $y_{i}$ can be assumed different from zero or one. We can therefore assume that both $ϕ μ_{i} k$ and $ϕ (1 - μ_{i}) k$ diverge. In this case, by (B.2), $f_{VIB} (y_{i}; μ_{i}, ϕ, k, p)$ admits the following asymptotic approximation

f_{VIB} (y_{i}; μ_{i}, ϕ, k, p) ≃ p e^{- ϕ k K (y_{i}, μ_{i})} \sqrt{ϕ k} + (1 - p) e^{- ϕ K (y_{i}, μ_{i})} \sqrt{ϕ} .

(B.4)

It follows that $f_{VIB} (y_{i}; μ_{i}, ϕ, k, p)$ diverges only if $μ_{i} = y_{i}$ (or $μ_{i} \to y_{i}$ ) and the rate of the divergence is at most $\max \{p \sqrt{ϕ k}, (1 - p) \sqrt{ϕ}\}$ .

On the other hand, if there exists an observation $y_{j}$ which is different from $μ_{j}$ , then $f_{VIB} (y_{j}; μ_{j}, ϕ, k, p)$ (see 8.4) converges to zero with rate

\max \{p e^{- ϕ k K (y_{j}, μ_{j})} \sqrt{ϕ k}, (1 - p) e^{- ϕ K (y_{j}, μ_{j})} \sqrt{ϕ}\} .

As, under condition (iii), $ϕ k \geq ε ϕ^{1 - α}$ for large ϕ, the exponential rate of convergence to zero of $f_{VIB} (y_{j}; μ_{j}, ϕ, k, p)$ dominates the divergence rate of $f_{VIB} (y_{i}; μ_{i}, ϕ, k, p)$ . It follows that the likelihood diverges only if $y_{i}$ is equal to $μ_{i}$ $\forall i$ . This can happen with probability one only if $n \leq k + 1$ . $□$

Appendix C Predictive intervals

Figure A1:

Scenario I: Posterior predictive means (cross symbols) and 90% predictive intervals for each subject (outlying observations are filled-in circles)

Figure A2:

Scenario II: Posterior predictive means (cross symbols) and 90% predictive intervals for each subject (outlying observations are filled-in circles)

Figure A3:

Scenario III: Posterior predictive means (cross symbols) and 90% predictive intervals for each subject (outlying observations are filled-in circles)

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship and/or publication of this article.

Funding

Research partially financially supported by the Italian Ministry of University and Research, Grants F.A. from the University of Milano-Bicocca.

References

Akaike

(1998) Information theory and an extension of the maximum likelihood principle, pp. 199–213. New York, NY: Springer. ISBN 978-1- 4612-1694-0. URL https://doi.org/10.1007/978-1-[4612-1694-0_15] (last accessed 23 November 2018).

Barnett

Lewis

(1995) Outliers in statistical data . Chichester: J. Wiley & Sons.

Bayer

Cribari-Neto

(2013) Bartlett corrections in beta regression models. Journal of Statistical Planning and Inference , 143, 531–47.

Bayer

Cribari-Neto

(2017) Model selection criteria in beta regression with varying dispersion. Communications in Statistics ’ Simulation and Computation , 46, 729–46.

Bayes

Bazn

Garca

(2012) A new robust regression model for proportions. Bayesian Analysis , 7, 841–66.

Box

Tiao

(1968) A bayesian approach to some outlier problems. Biometrika , 55, 119–29.

Branscum

Johnson

Thurmond

(2007) Bayesian beta regression: Applications to household expenditure data and genetic distance between foot-and-mouth disease viruses. Australian & New Zealand Journal of Statistics , 49, 287–301.

Brooks

(2002) Discussion on the article by spiegelhalter, best, carlin and van der linde.

Brooks

Gelman

Jones

Meng

X-L

(2011) Handbook of Markov chain Monte Carlo . Boca Raton, FL: CRC press.

10.

Canterle

Bayer

(2017) Variable dispersiong beta regressions with parametric link functions. Statistical Papers . doi: 10.1007/s00362-017-0885-9

11.

Cribari-Neto

Souza

(2012) Testing inference in variable dispersion beta regressions. Journal of Statistical Computation and Simulation , 82, 1827–43.

12.

Duane

Kennedy

Pendleton

Roweth

(1987) Hybrid Monte Carlo. Physics letters B , 195, 216–22.

13.

Duong

(2007) ks: Kernel density estimation and kernel discriminant analysis for multivariate data in r. Journal of Statistical Software , 21, 1–16.

14.

Espinheira

Ferrari

Cribari-Neto

(2008a) Influence diagnostics in beta regression. Computational Statistica and Data Analysis , 52, 4417–31.

15.

Espinheira

Ferrari

Cribari-Neto

(2008b) On beta regression residuals. Journal of Applied Statistics , 35, 407–19.

16.

Ferrari

Cribari-Neto

(2004) Beta regression for modelling rates and proportions. Journal of Applied Statistics , 31, 799–815.

17.

Ferrari

Pinheiro

(2011) Improved likelihood inference in beta regression. Journal of Statistical Computation and Simulation , 81, 431–43.

18.

Ferrari

Espinheira

Cribari-Neto

(2011) Diagnostic tools in beta regression with varying dispersion. Statistica Neerlandica , 65, 337–51.

19.

García

Pérez

van Dorp

(2011) Modeling heavy-tailed, skewed and peaked uncertainty phenomena with bounded support. Statistical Methods & Applications , 20, 463–86.

20.

Gelman

Carlin

Stern

Rubin

(2014) Bayesian data analysis, Vol. 2. Milton Park: Taylor & Francis.

21.

Ghosh

(2017) Robust inference under the beta regression model with application to health care studies. Statistical Methods in Medical Research . doi: 10.1177/096228 0217738142

22.

Hahn

(2008) Mixture densities for project management activity times: A robust approach to (PERT). European Journal of Operational Research , 188, 450–59.

23.

Hoffman

Gelman

(2014) The No-U-Turn Sampler: Adaptively setting path lengths in Hamiltonian Monte Carlo. Journal of Machine Learning Research , 15, 1593–1623.

24.

Johnson

Kotz

Balakrishnan

(1995) Continuous Univariate Distributions, Vol. 2. Hoboken, NJ: Wiley.

25.

Kieschnick

McCullough

(2003) Regression analysis of variates observed on (0, 1): Percentages, proportions and fractions. Statistical Modelling , 3, 193–213.

26.

Liu

Eugenio

(2016) A review and comparison of bayesian and likelihood-based inferences in beta regression and zero-or- one-inflated beta regression. Statistical Methods in Medical Research , 27, 1024–44. doi: 10.1177/0962280216650699

27.

McCullagh

Nelder

(1989) Generalized Linear Models , Vol. 37. Boca Raton, FL: CRC press.

28.

Mengersen

Robert

Guihenneuc- Jouyaux

(1999) Mcmc convergence diagnostics: A review. Bayesian Statistics , 6, 415–40.

29.

Migliorati

Ongaro

Monti

(2017) A structured Dirichlet mixture model for compositional data: Inferential and applicative issues. Statistics and Computing , 27, 963–83. doi: 10.1007/s11222-016-9665-y

30.

Migliorati

Di Brisco

Ongaro

(2018) A new regression model for bounded responses. Bayesian Anal ., 13, 845–72. doi: 10.1214/17-BA1079

31.

Ongaro

Migliorati

(2013) A generalization of the Dirichlet distribution. Journal of Multivariate Analysis , 114, 412–26.

32.

Ospina

Ferrari

(2010) Inflated beta distributions. Statistical Papers , 51, 111–26.

33.

Ospina

Ferrari

(2012) A general class of zero-or-one inflated beta regression models. Computational Statistics & Data Analysis , 56, 1609–23. URL https://doi.org/10.[1016/j.csda.2011.10.005]

34.

Ospina

Cribari-Neto

Vasconcellanos

(2006) Improved point and interval estimation for a beta regression model. Computational Statistica and Data Analysis , 51, 960–81.

35.

Paolino

(2001) Maximum likelihood estimation of models with beta-distributed dependent variables. Political Analysis , 9, 325–46.

36.

Pawlowsky-Glahn

Egozcue

Tolosana-Delgado

(2015) Modeling and Analysis of Compositional Data . Hoboken, NJ: John Wiley & Sons.

37.

Peña

Guttman

(1993) Comparing probabilistic methods for outlier detection in linear models. Biometrika , 80, 603–10.

38.

Pereira

Botter

Sandoval

(2013) A regression model for special proportions. Statistical Modelling 13, 125–51. doi: 10.1177/1471082X13478274

39.

Pereira

Cribari-Neto

(2014) Detecting model misspecification in inflated beta regressions. Communications in Statistics–Simulation and Computation , 43, 631–56.

40.

Rousseeuw

Leroy

(2005) Robust regression and outlier detection , Vol 589. Hoboken, NJ: John Wiley & Sons.

41.

Schwarz

(1978) Estimating the dimension of a model. The Annals of Statistics , 6, 461–64.

42.

Simas

Barreto-Souza

Rocha

(2010) Improved estimators for a general class of beta regression models. Computational Statistics & Data Analysis , 54, 348–66.

43.

Simonoff

(2012) Smoothing methods in statistics . Berlin: Springer Science & Business Media.

44.

Smithson

Verkuilen

(2006) A better lemon squeezer? maximum-likelihood regression with beta-distributed dependent variables. Psychological Methods , 11, 54.

45.

Spiegelhalter

Best

Carlin

Van Der Linde

(2002) Bayesian measures of model complexity and fit (with discussion). Journal of the Royal Statistical Society: Series B , 64, 583–639.

46.

Stan Development Team (2016) Stan Modeling Language Users Guide and Reference Manual. URL http://mc-stan.org/

47.

Vehtari

Gelman

Gabry

(2017) Practical bayesian model evaluation using leave-one-out cross-validation and waic. Statistics and Computing , 27, 1413–32.

48.

Watanabe

(2010, December) Asymptotic equivalence of bayes cross validation and widely applicable information criterion in singular learning theory. Journal of Machine Learning Research , 11, 3571–94.

Robustness against outliers: A new variance inflated regression model for proportions

Abstract

Keywords

1 Introduction

2 A new variance inflated beta model for bounded continuous variables

2.1 Introduction

Some examples of VIB distributions

Figure 3:

Scenario I: Posterior mean (solid line) and 95% credible interval (dashed line) of β 0 (top panel) and β 1 (bottom panel) for each regression model as δ varies

Scenario II: Posterior mean (solid line) and 95% credible interval (dashed line) of β 0 (top panel) and β 1 (bottom panel) for each regression model as δ varies

Scenario III: Posterior mean (solid line) and 95% credible interval (dashed line) of β 0 (top panel) and β 1 (bottom panel) for each regression model as δ varies

Figure 6:

KL divergence measure for VIBReg (solid line), FBReg (dashed line), BRReg (dotted line) and BReg (dotdash line)

Curve distances for VIBReg (solid line), FBReg (dashed line), BRReg (dotted line) and BReg (dotdash line)

Mean values of WAIC1 criterion for VIBReg (solid line), FBReg (dashed line), BRReg (dotted line) and BReg (dotdash line) models as δ varies in the three scenarios

Estimated response variance ratio with respect to FBReg in the three scenarios. VIBReg (solid line), BRReg (dotted line) and BReg (dotdash line)

5 Application to Australian Institute of Sport data

5.1 Regression model with one covariate

Posterior parameter means and fitting measure WAIC1 for the four models under the regression model with one covariate

Standardized residuals for each subject and for each model on the entire dataset

Estimated response variance ratio with respect to FBReg. VIBReg (solid line), BRReg (dotted line) and BReg (dotdash line)

Posterior predictive means (cross symbols) and 90% predictive intervals for each subject and for each model (observed values in circles, filled-in for outliers) on the entire dataset

Posterior estimates of inflation parameter k and precision parameters separately for males and females

Supplementary materials

Acknowledgements

Appendices

Appendix A Proof of Proposition 1 (Identifiability of the VIBReg model)

Appendix B Proof of Proposition 2 (Boundedness of the VIBReg model likelihood)

Appendix C Predictive intervals

Scenario I: Posterior predictive means (cross symbols) and 90% predictive intervals for each subject (outlying observations are filled-in circles)

Scenario II: Posterior predictive means (cross symbols) and 90% predictive intervals for each subject (outlying observations are filled-in circles)

Scenario III: Posterior predictive means (cross symbols) and 90% predictive intervals for each subject (outlying observations are filled-in circles)

Declaration of conflicting interests

Funding

References

Scenario I: Posterior mean (solid line) and 95% credible interval (dashed line) of $β_{0}$ (top panel) and $β_{1}$ (bottom panel) for each regression model as $δ$ varies

Scenario II: Posterior mean (solid line) and 95% credible interval (dashed line) of $β_{0}$ (top panel) and $β_{1}$ (bottom panel) for each regression model as $δ$ varies

Scenario III: Posterior mean (solid line) and 95% credible interval (dashed line) of $β_{0}$ (top panel) and $β_{1}$ (bottom panel) for each regression model as $δ$ varies

Mean values of WAIC₁ criterion for VIBReg (solid line), FBReg (dashed line), BRReg (dotted line) and BReg (dotdash line) models as $δ$ varies in the three scenarios

Posterior parameter means and fitting measure WAIC₁ for the four models under the regression model with one covariate

Posterior estimates of inflation parameter $k$ and precision parameters separately for males and females