Robust estimation in beta regression via maximum L_q-likelihood

Random phenomena that produce data in the unit interval $(0,1)$ are common in many research areas, including medicine, finance, economics, environment, hydrology, and psychology. Data of this type are usually rates, proportions, percentages, or fractions. A few examples are unemployment rates, the proportion of exports in total sales, the fraction of surface covered by vegetation, Gini’s index, the loss given default, health-related quality of life scores, and neonatal mortality rates. Beta regression models (Ferrari and Cribari-Neto, 2004; Smithson and Verkuilen, 2006) are widely employed for modeling continuous data bounded within the unit interval. The response variable is assumed to follow a beta distribution indexed by the mean and a precision parameter, which in turn are modeled through a variety of regression structures.

Let the response variable $y_{i}$ have a beta distribution with mean $0<\mu_{i}<1$ and precision parameter $\phi_{i}>0$. The probability density function of $y_{i}$ is

where $B(.,.)$ is the beta function. We write $y_{i}\sim{\rm Beta}(\mu_{i},\phi_{i})$, $i=1,\ldots,n$, and assume that $y_{1},\ldots,y_{n}$ are independent. Since ${\rm Var}(y_{i})=\mu_{i}(1-\mu_{i})/(1+\phi_{i})$, the parameter $\phi_{i}$ is taken as a precision parameter because it is inversely related to the dispersion of $y_{i}$, for fixed mean $\mu_{i}$. The beta density (1) assumes different forms depending on the parameters. It has a single mode in $(0,1)$ if $\mu_{i}\phi_{i}>1$ and $(1-\mu_{i})\phi_{i}>1$, and its possible forms include J, inverted-J, tilde, inverted tilde, bathtub and uniform shapes. When $\mu_{i}\phi_{i}<1$ or $(1-\mu_{i})\phi_{i}<1$ the beta density is unbounded at one or both of the boundaries. We will assume that the regression structures for $\mu_{i}$ and $\phi_{i}$ are, respectively,

where ${\bm{\beta}}=(\beta_{1},\ldots,\beta_{p_{1}})^{\top}\in\mathbb{R}^{p_{1}}$ and $\mbox{\boldmath$\gamma$\unboldmath}=(\gamma_{1},\ldots,\gamma_{p_{2}})^{\top}\in\mathbb{R}^{p_{2}}$ are vectors of unknown regression parameters, $\mbox{\boldmath$X$\unboldmath}_{i}=(x_{i1},\ldots,x_{ip_{1}})^{\top}$ and $\mbox{\boldmath$Z$\unboldmath}_{i}=(z_{i1},\ldots,z_{ip_{2}})^{\top}$ are vectors of covariates of lengths $p_{1}$ and $p_{2}$ ($p_{1}+p_{2}=p<n$), respectively, which are assumed fixed and known. In addition, $g_{\mu}(\cdot):(0,1)\rightarrow\mathbb{R}$ and $g_{\phi}(\cdot):(0,\infty)\rightarrow\mathbb{R}$ are strictly monotonic and twice differentiable link functions. The model introduced by Ferrari and Cribari-Neto (2004) is a constant precision (or, equivalently, constant dispersion) beta regression model in the sense that $\phi_{i}=\phi$, for $i=1,\ldots,n$.

Inference for the beta regression model (1)-(3) is commonly based on the maximum likelihood estimator (MLE). However, it may be highly influenced by the presence of outlying observations in the data.

As an illustration of the sensitivity of the MLE to outliers, we consider a dataset on risk management practices of 73 firms introduced and analyzed by Schmit and Roth (1990); see also Gómez-Déniz et al. (2014). The response variable is Firm cost, defined as premiums plus uninsured losses as a percentage of the total assets. Smaller values of Firm cost are attributed to firms that have a good risk management performance. We take only two covariates, namely Ind cost, a measure of the firm’s industry risk, and Size log, the logarithm of its total asset. Figure 1 displays the scatter plot of Firm cost versus Ind cost along with the fitted lines based on the MLE for a beta regression with varying precision for the full data and the data without the most eye-catching outliers, labeled in the plots as observations 15, 16 and 72.²²2The scatter plots were produced by setting the value of the remaining covariate at its sample median value. These observations have disproportionate effects on the fitted regression curves.

To deal with atypical observations, models based on distributions that are more flexible than the beta distribution may be employed. For instance, Bayes et al. (2012), Migliorati et al. (2018), and Di Brisco et al. (2020) proposed different kinds of mixtures of beta distributions. However, more flexibility comes at the cost of model complexity and additional parameters, making model specification, inference, and interpretation more difficult. Our approach keeps the model simple, i.e. the beta distribution remains as the postulated model, but the inference procedure is replaced by a robust method.

Ghosh (2019) proposed a robust estimator in beta regression models based on the minimization of an empirical density power divergence introduced by Basu et al. (1998). The estimator, called the minimum density power divergence estimator (MDPDE), is similar to the approach we will propose here. The main difference between the two estimation procedures is that they are based on two different families of divergences that have the Kullback-Leibler divergence as a special case. In other words, both have the MLE as a special case. Ghosh (2019) presented some simulations of the MDPDE, but only for constant precision beta regression models with fixed values for the tuning constant. Here, we propose a data-driven method to select the tuning constant of the MDPDE which ensures full efficiency in the absence of contamination. Additionally, we evaluate the performance of the MDPDE under constant and varying precision using the proposed selection algorithm.

The chief goal of this paper is to propose a new robust estimator under the beta regression model (1)-(3) with a scheme determined by the data to choose its optimum tuning constant. In Section 2, we introduce a robust estimator based on the maximization of a reparameterized L_q-likelihood, and derive some properties. In Section 3, we introduce a data-driven method to the selection of the tuning constant. In Section 4, we provide Monte Carlo evidence on the proposed robust inference approach. Section 5 contains applications of the robust estimation for three datasets. Finally, Section 6 closes the paper with some concluding remarks. Some technical details are left for an appendix. The on-line Supplementary Material contains further technical information and additional numerical results. The codes for the simulations and applications and the datasets are available at https://github.com/terezinharibeiro/RobustBetaRegression.

Let $y_{1},\ldots,y_{n}$ be $n$ independent observations, where $y_{i}$ comes from a postulated model indexed by a common unknown parameter $\mbox{\boldmath$\theta$\unboldmath}\in\Theta\subset\mathbb{R}_{p}$. For the beta regression model, the postulated density for $y_{i}$ is given by (1)-(3), and we will denote it, from now on, by $f_{\bm{\theta}}(\cdot;\mu_{i},\phi_{i})$. Let $\ell_{q}(\mbox{\boldmath$\theta$\unboldmath})$ denote the L_q-likelihood (Ferrari and Yang, 2010),

where $0<q\leq 1$ is a tuning constant, and

is the Box-Cox transformation (Box and Cox, 1964), also called distorted logarithm. The population version of $-\ell_{q}(\mbox{\boldmath$\theta$\unboldmath})/n$ with respect to the density $g(\cdot)$ is

which is the $q$-entropy of $f_{\bm{\theta}}$ with respect to $g$ (Tsallis, 1988) and reduces to the Shannon entropy when $q=1$. For $q=1$, $\ell_{q}(\mbox{\boldmath$\theta$\unboldmath})$ is the usual log-likelihood function and maximizing $\ell_{q}(\mbox{\boldmath$\theta$\unboldmath})$ leads to the maximum likelihood estimator (MLE).

The estimator obtained by maximizing (4), will be denoted by $\widetilde{{\bm{\theta}}}_{q}$. It comes from the estimating equation

where $U(y_{i},{\bm{\theta}})=\nabla_{\bm{\theta}}\log[f_{\bm{\theta}}(y_{i};\mu_{i},\phi_{i})]$, with $\nabla_{{\bm{\theta}}}$ denoting the gradient with respect to $\theta$. The score vector for $\theta$ corresponding to the $i$-th observation is given by

where $y_{i}^{\star}=\log(y_{i}/(1-y_{i}))$, $y_{i}^{\dagger}=\log(1-y_{i})$, $\mu_{i}^{\star}=\mbox{E}(y_{i}^{\star})=\psi(\mu_{i}\phi_{i})-\psi((1-\mu_{i})\phi_{i})$, and $\mu_{i}^{\dagger}=\mbox{E}(y_{i}^{\dagger})=\psi((1-\mu_{i})\phi_{i})-\psi(\phi_{i})$, with $\psi(\cdot)$ denoting the digamma function.

The score vector of the $i$-th observation is weighted by $f_{\bm{\theta}}(y_{i};\mu_{i},\phi_{i})^{1-q}$, which depends on the postulated model and the tuning constant $q$. A choice of $q$ not too close to 1 produces a robust estimation procedure, because observations that are inconsistent with the assumed model receive smaller weights. Furthermore, (5) defines an M-estimation procedure, first introduced by Huber (1964). The class of M-estimators is widely used for obtaining robust estimators, and has interesting properties such as the asymptotic normality (Huber, 1981).

The estimating function $U(\mbox{\boldmath$y$\unboldmath};\mbox{\boldmath$\theta$\unboldmath})=\sum_{i=1}^{n}U(y_{i},{\bm{\theta}})f_{\bm{\theta}}(y_{i};\mu_{i},\phi_{i})^{1-q}$ is biased unless $q=1$. Hence, the estimator $\widetilde{{\bm{\theta}}}_{q}$ is not Fisher-consistent. Ferrari and La Vecchia (2012) showed that a suitable calibration function may be applied to rescale the parameter estimates and to achieve Fisher-consistency if the family of postulated distributions is closed under the power transformation. Given a density $h$ and a constant $\alpha\in(0,\infty)$, the power transformation is defined as

provided that $\int h(y)^{\alpha}dy<\infty$, for $0<\alpha<\infty$. For a family of postulated densities $\{h_{\bm{\theta}}(\cdot),{\bm{\theta}}\in\Theta\}$ that is closed under (8) for all $0<\alpha<\infty$, let $\tau_{\alpha}(\mbox{\boldmath$\theta$\unboldmath}):\Theta\mapsto\Theta$ be an invertible continuous function satisfying $h_{\tau_{\alpha}(\bm{\theta})}(y)=h_{\bm{\theta}}^{(\alpha)}(y)$ for all $y$ in the support, which is assumed not to depend on $\theta$. Note that the function $\tau_{\alpha}(\cdot)$ maps each density $h_{\bm{\theta}}(\cdot)$ with a unique power-transformed density $h_{\tau_{\alpha}(\bm{\theta})}(\cdot)$.

The $q$-entropy $H_{q}(f_{\bm{\theta}},g)$ is related with the family of divergences of $f_{\bm{\theta}}$ with respect to $g$

Ferrari and La Vecchia (2012, Lemma 1) showed that $D_{q}(f_{\bm{\theta}},g^{(1/q)})\ \propto\ H_{q}(f_{\bm{\theta}},g)-H_{q}(g^{(1/q)},g)$. Consequently, the minimizer $\widetilde{{\bm{\theta}}}_{q}$ of the empirical $q$-entropy with respect to $\theta$ equals the minimizer of the empirical version of $D_{q}(f_{\bm{\theta}},g^{(1/q)})$. In the minimization of $D_{q}(f_{\bm{\theta}},g^{(1/q)})$, the target density is the transformed density $g^{(1/q)}$ instead of $g$. Ferrari and La Vecchia (2012, Proposition 1) show that $\widetilde{{\bm{\theta}}}_{q}$ is not Fisher-consistent for $\theta$ but it is Fisher-consistent for $\tau_{1/q}(\mbox{\boldmath$\theta$\unboldmath})$. Hence, $\widehat{{\bm{\theta}}}_{q}=\tau_{q}(\widetilde{{\bm{\theta}}}_{q})$ is Fisher-consistent for $\theta$. Alternatively, $\widehat{{\bm{\theta}}}_{q}$ is obtained by maximizing the L_q-likelihood in the parameterization $\tau_{1/q}(\bm{\theta})$ (La Vecchia et al., 2015). Another proof of the Fisher-consistency of $\widehat{{\bm{\theta}}}_{q}$ is given in the Supplementary Material.

For the beta regression model (1)-(3), $\widehat{{\bm{\theta}}}_{q}$ is the maximizer of

where $f^{(1/q)}_{\bm{\theta}}(y_{i};\mu_{i},\phi_{i})=f_{\bm{\theta}}(y_{i};\mu_{i,q^{-1}},\phi_{i,q^{-1}})$, for $0<q<1$, with

provided that $0<\mu_{i,q^{-1}}<1$ and $\phi_{i,q^{-1}}>0$, or equivalently, $\mu_{i}\phi_{i}>1-q$ and $(1-\mu_{i})\phi_{i}>1-q$. Then, $f_{\bm{\theta}}(y_{i};\mu_{i},\phi_{i})$ satisfies (8) for all $\alpha=1/q>0$, if $\mu_{i}\phi_{i}\geq 1$ and $(1-\mu_{i})\phi_{i}\geq 1$, i.e. if $f_{\bm{\theta}}(y_{i};\mu_{i},\phi_{i})$ is bounded. From (2) and (3), $f^{(1/q)}_{\bm{\theta}}(y_{i};\mu_{i},\phi_{i})$ corresponds to a modified beta regression model defined by the beta density (1) with mean and precision submodels given respectively by

and we denote it by $f_{\bm{\theta}}^{*}(y_{i};\mu_{i},\phi_{i})$.

Therefore, we propose the estimator that is obtained through the maximization of

which leads to the following estimating equation

where $U^{*}(y_{i},{\bm{\theta}})=\nabla_{\bm{\theta}}\log[f^{*}_{\bm{\theta}}(y_{i};\mu_{i},\phi_{i})]$ being a modified score vector for $\theta$ corresponding to the $i$-th observation given by

The final estimator, $\widehat{{\bm{\theta}}}_{q}$, will be simply called the surrogate maximum likelihood estimator (SMLE).

From the theory of M-estimation, $\widehat{{\bm{\theta}}}_{q}\mbox{ }{\overset{a}{\sim}}\mbox{ }\mbox{N}(\mbox{\boldmath$\theta$\unboldmath},V_{q}(\mbox{\boldmath$\theta$\unboldmath})),$ where $\overset{a}{\sim}$ denotes asymptotic distribution and

with

see Ferrari and La Vecchia (2012, Section 3.2). For the beta regression model (1)-(3), $J_{q}(\mbox{\boldmath$\theta$\unboldmath})$ and $K_{q}(\mbox{\boldmath$\theta$\unboldmath})$ are given in the Appendix. The matrix $V_{q}(\mbox{\boldmath$\theta$\unboldmath})$ is well-defined provided that $0<\mu_{i,2-q}<1$ and $\phi_{i,2-q}>0$, or equivalently $\mu_{i}\phi_{i}>2(1-q)/(2-q)$ and $(1-\mu_{i})\phi_{i}>2(1-q)/(2-q)$. These conditions hold for all $0<q\leq 1$, if $\mu_{i}\phi_{i}\geq 1$ and $(1-\mu_{i})\phi_{i}\geq 1$.

Likelihood-based tests are usually constructed from the three well-known test statistics, namely the likelihood ratio, Wald, and score statistics, and they rely on maximum likelihood estimation. Following Heritier and Ronchetti (1994) (see also Heritier et al. (2009, Section 2.5.3)), robust versions of these statistics may be derived by replacing the MLE, the log-likelihood function, the score vector, and the asymptotic covariance matrix respectively by the SMLE, the reparameterized $L_{q}$-likelihood in (10), the estimating function in (11), and the matrix in (14). In this paper, we focus on the robust Wald statistic, referred here as the Wald-type statistic, of the null hypothesis $\beta_{k}=\beta_{k}^{(0)}$ against a two-sided alternative defined as $\{(\widehat{\beta}_{kq}-\beta_{k}^{(0)})/{\rm se}(\widehat{\beta}_{kq})\}^{2}$, which has a $\chi^{2}_{1}$ distribution under the null hypothesis. Here, $\widehat{\beta}_{kq}$ is the SMLE of $\beta_{k}$ and ${\rm se}(\widehat{\beta}_{kq})$ is its asymptotic standard error that is computed from (14).

Since the MLE ($\widehat{\mbox{\boldmath$\theta$\unboldmath}}$) and the SMLE are M-estimators, their influence functions for the beta regression model (1)-(3) are $\mbox{IF}(y;\widehat{\mbox{\boldmath$\theta$\unboldmath}})=-K_{1}({\mbox{\boldmath$\theta$\unboldmath}})^{-1}U(y;{\mbox{\boldmath$\theta$\unboldmath}})$ and $\mbox{IF}(y;\widehat{{\bm{\theta}}}_{q})=-J_{q}(\mbox{\boldmath$\theta$\unboldmath})^{-1}U^{*}(y,{\bm{\theta}})f^{*}_{\bm{\theta}}(y;\mu,\phi)^{1-q}$, where $U(y;{\mbox{\boldmath$\theta$\unboldmath}})$ and $U^{*}(y,{\bm{\theta}})$ are given in (7) and (13), respectively, $K_{1}({\mbox{\boldmath$\theta$\unboldmath}})$ is the Fisher information matrix (see the Appendix), and $J_{q}(\mbox{\boldmath$\theta$\unboldmath})$ is given in (14). The influence function measures the asymptotic bias of the estimator caused by an infinitesimal contamination introduced in the data point $y$ (Hampel et al., 2011, Section 2.1.b). The individual score vector $U(y;{\mbox{\boldmath$\theta$\unboldmath}})$ is unbounded when $y\rightarrow 0$ or $y\rightarrow 1$, i.e., for $y$ tending to one of the boundaries of the support set. Hence, the influence function of the MLE is also unbounded. On the other hand, for bounded beta densities the weighted modified score vector $U^{*}(y,{\bm{\theta}})f^{*}_{\bm{\theta}}(y;\mu,\phi)^{1-q}$ and the influence function of the SMLE are bounded, which implies that the SMLE is B-robust (Hampel et al., 2011, Section 2.1). Figure 2 illustrates the behavior of the score vector and of the weighted modified score vector for the i.i.d. case. Note that the components of the score vector grow or decrease unboundedly when $y$ approaches zero or one, while the components of the weighted modified score vector has a redescending shape, approaching zero at the boundaries. This means that observations close to zero or one may be influential in the MLE fit but not in the SMLE fit.

We also studied the behavior of the change-of-variance function, that measures the bias on the covariance matrix caused by an infinitesimal contamination in the data point $y$. For the MLE, the change-of-variance function is unbounded. However, as presented in Ferrari and La Vecchia (2012, Proposition 2), the change-of-variance function of SMLE is bounded if the weighted modified score vector and its first derivative are bounded. We prove the boundedness of theses quantities for bounded beta densities, which implies that the SMLE is V-robust (Hampel et al., 2011, Section 2.5); see the Supplementary Material for details.

Ghosh (2019) proposed a robust estimator in beta regression models that minimizes an empirical version of the density power divergence introduced by Basu et al. (1998). The minimum density power divergence estimator (MDPDE) for the beta regression model (1)-(3) comes from the estimating equations

where $E_{i,q}(\bm{\theta})={\rm E}(U(y_{i},{\bm{\theta}})f_{\bm{\theta}}(y_{i};\mu_{i},\phi_{i})^{1-q})$. The components of $E_{i,q}(\bm{\theta})$ may be found in Ghosh (2019), and the asymptotic covariance matrix of the MDPDE is given in the Supplementary Material. The expectation $E_{i,q}(\bm{\theta})$ requires that $\mu_{i}\phi_{i}>(1-q)/(2-q)$ and $(1-\mu_{i})\phi_{i}>(1-q)/(2-q)$. Besides, the asymptotic covariance matrix of the MDPDE requires that $\mu_{i}\phi_{i}>2(1-q)/(3-2q)$ and $(1-\mu_{i})\phi_{i}>2(1-q)/(3-2q)$. Hence, these quantities are well-defined, for all $0<q\leq 1$, only for bounded beta densities.

Note that the estimating functions in (15) are unbiased. They agree with those derived from the $L_{q}$-likelihood in (5) except for the term $E_{i,q}(\bm{\theta})$ that plays the role of correcting the score bias. Recall that the score bias correction in the $L_{q}$-likelihood approach adopted here is achieved by rescaling the parameter estimates or, equivalently by reparameterizing the $L_{q}$-likelihood.

Ghosh (2019) presented simulations of the MDPDE only for constant precision beta regression models with fixed values for the tuning constant. However, in practical applications the choice of such a constant becomes a crucial issue. In Section 3, we will propose a data-driven method to select the tuning constant of the SMLE and the MDPDE which ensures full efficiency in the absence of contamination.

An important aspect in real data applications is the choice of the optimal value of the tuning constant $q$. Smaller values of $q$ increase robustness at the expense of reduced efficiency. La Vecchia et al. (2015) suggest to select the tuning constant closest to 1 (i.e., closest to MLE) such that the estimates of the parameters are sufficiently stable, ensuring full efficiency in the absence of contamination. This is the approach we follow here, but we standardize the estimates by ‘$\sqrt{n}\times\mbox{std. error}$’ to simultaneously remove the effect of the sample size and of the magnitude of estimates of different parameters. We then propose a data-driven algorithm that starts with $q_{0}=1$ (MLE) and follows the steps below.

1. 1.

   Define an equally spaced grid for $q$: $q_{0}>q_{1}>q_{2}>\cdots>q_{m},$ with $q_{m}\geq q_{\rm min}$;

2. 2.

   compute the estimates $\widehat{\theta}^{1}_{q_{k}},\ldots,\widehat{\theta}^{p}_{q_{k}}$ with the corresponding asymptotic standard errors, $\mbox{se}(\widehat{\theta}^{1}_{q_{k}}),\ldots,$ $\mbox{se}(\widehat{\theta}^{p}_{q_{k}})$;

3. 3.

   compute the vectors of standardized estimates,

   $$\mbox{\boldmath$z$\unboldmath}_{q_{k}}=\left(\frac{\widehat{\theta}^{1}_{q_{k}}}{\sqrt{n}\ \mbox{se}(\widehat{\theta}^{1}_{q_{k}})},\cdots,\frac{\widehat{\theta}^{p}_{q_{k}}}{\sqrt{n}\ \mbox{se}(\widehat{\theta}^{p}_{q_{k}})}\right)^{\top},$$

   and the standardized quadratic variations (SQV) defined by

   $$\mbox{SQV}_{q_{k}}=p^{-1}||\mbox{\boldmath$z$\unboldmath}_{q_{k}}-\mbox{\boldmath$z$\unboldmath}_{q_{k+1}}||;$$

4. 4.

   for each $q_{k}$, check whether the stability condition defined by $\mbox{SQV}_{q_{k}}<L$ is satisfied, where $L>0$ is a pre-defined threshold;

5. 5.

   if the stability condition is satisfied for all $q_{k}$, set the optimal value of $q$ at $q^{*}=\max q_{k}$;

6. 6.

   if the stability condition is not satisfied for some $q_{k}$, set a new grid starting from the smallest $q_{k}$ for which the condition does not hold;

7. 7.

   repeat steps 1-6 until achieving stability of the estimates for the current grid or reaching the minimum value $q_{\rm min}$; if stability is achieved, set the optimal value of $q$ at $q^{*}=\max q_{k}$, otherwise at $q^{*}=1$.

This procedure chooses the optimum value for the tuning constant that is closest to 1 (MLE) and still guarantees stability of the estimates for $m$ consecutive equally spaced values of $q$. It may occur that such an optimal $q$ is not reached. In this case, the algorithm sets the optimal $q$ at 1 i.e. the MLE is chosen. This choice may seem unreasonable but recall that the SMLE is not expected to be robust for unbounded beta densities, which might lead to unstable estimates even for reasonably small $q$. In such a situation, it is not advisable to replace the MLE by the SMLE. Our experience with simulated samples suggests to set the grid spacing at 0.02, the size of the grids at $m=3$, the minimum value of $q$ at $q_{\rm min}=0.5$ and the threshold at $L=0.02$. These values were used in the simulations and applications reported below and worked very well.

To evaluate the performance and compare the MLE, SMLE and MDPDE, Monte Carlo simulations in the presence and absence of contamination in the data were carried out. The sample sizes considered are $n=40,80,160,320$. The covariate values were set for the sample size $n=40$ and replicated twice, four times and eight times to obtain the covariate values corresponding to the sizes $n=80,$ 160, 320, respectively. This scheme guarantees that the heteroskedasticity intensity, $\mbox{max}(\phi_{i})/\mbox{min}(\phi_{i})$, be constant for all the sample sizes (Espinheira et al., 2017, 2019). For all the scenarios we employ logit and logarithmic links in the mean and precision submodels, respectively (2) and (3). The models include intercepts, i.e. $x_{1i}=z_{1i}=1$, for all $i=1,2,\ldots,n$, and the other covariates for the mean submodel are taken as random draws from a standard uniform distribution and kept constant throughout the simulated samples. For non-constant precision scenarios, the same covariates are used in the mean and precision submodels. The percentage of contamination in the sample for all the scenarios was set at 5%. The selection of the tuning constant $q$ for SMLE and MDPDE was conducted following the proposed algorithm described in Section 3. All the results are based on 1,000 Monte Carlo replications and were carried out using the R software environment (R Core Team, 2020).

Different configurations of parameter values and patterns of data contamination were considered. We now describe the results for two scenarios; results for other scenarios are collected in the Supplementary Material.