Robust Relative Error Estimation

In regression analysis, many analysts use the (penalized) least squares estimation, which aims at minimizing the mean squared prediction error [1]. On the other hand, the relative (percentage) error is often more useful and/or adequate than the mean squared error. For example, in econometrics, the comparison of prediction performance between different stock prices with different units should be made by relative error; we refer to [2,3] among others. Additionally, the prediction error of photovoltaic power production or electricity consumption is evaluated by not only mean squared error but also relative error (see, e.g., [4]). We refer to [5] regarding the usefulness and importance of the relative error.

In relative error estimation, we minimize a loss function based on the relative error. An advantage of using such a loss function is that it is scale free or unit free. Recently, several researchers have proposed various loss functions based on relative error [2,3,6,7,8,9]. Some of these procedures have been extended to the nonparameteric model [10] and random effect model [11]. The relative error estimation via the $L_1$ regularization, including the least absolute shrinkage and operator (lasso; [12]), and the group lasso [13], have also been proposed by several authors [14,15,16], to allow for the analysis of high-dimensional data.

In practice, a response variable $y(>0)$ can turn out to be extremely large or close to zero. For example, the electricity consumption of a company may be low during holidays and high on exceptionally hot days. These responses may often be considered to be outliers, to which the relative error estimator is sensitive because the loss function diverges when $y\to \infty$ or $y\to 0$. Therefore, a relative error estimation that is robust against outliers must be considered. Recently, Chen et al. [8] discussed the robustness of various relative error estimation procedures by investigating the corresponding distributions, and concluded that the distribution of least product relative error estimation (LPRE) proposed by [8] has heavier tails than others, implying that the LPRE might be more robust than others in practical applications. However, our numerical experiments show that the LPRE is not as robust as expected, so that the robustification of the LPRE is yet to be investigated from the both theoretical and practical viewpoints.

To achieve a relative error estimation that is robust against outliers, this paper employs the $\gamma$-likelihood function for regression analysis by Kawashima and Fujisawa [17], which is constructed by the $\gamma$-cross entropy [18]. The estimating equation is shown to have a redescending property, a desirable property in robust statistics literature [19]. To find a minimizer of the negative $\gamma$-likelihood function, we construct a majorize-minimization (MM) algorithm. The loss function of our algorithm at each iteration is shown to be convex, although the original negative $\gamma$-likelihood function is nonconvex. Our algorithm is guaranteed to decrease the objective function at each iteration. Moreover, we derive the asymptotic normality of the corresponding estimator together with a simple consistent estimator of the asymptotic covariance matrix, which enables us to straightforwardly create approximate confidence sets. Monte Carlo simulation is conducted to investigate the performance of our proposed procedure. An analysis of electricity consumption data is presented to illustrate the usefulness of our procedure. Supplemental material includes our R package rree (robust relative error estimation), which implements our algorithm, along with a sample program of the rree function.

The reminder of this paper is organized as follows: Section 2 reviews several relative error estimation procedures. In Section 3, we propose a relative error estimation that is robust against outliers via the $\gamma$-likelihood function. Section 4 presents theoretical properties: the redescending property of our method and the asymptotic distribution of the estimator, the proof of the latter being deferred to Appendix A. In Section 5, the MM algorithm is constructed to find the minimizer of the negative $\gamma$-likelihood function. Section 6 investigates the effectiveness of our proposed procedure via Monte Carlo simulations. Section 7 presents the analysis on electricity consumption data. Finally, concluding remarks are given in Section 8.

Suppose that $x_i$ = ${(x_{i1},\cdots ,x_{ip})}^T$$(i=1,...,n)$ are predictors and $\mathbf{y}={(y_1,...,y_n)}^T$ is a vector of positive responses. Consider the multiplicative regression model where $\mathbf{\beta }={({\beta }_1,\cdots ,{\beta }_p)}^T$ is a p-dimensional coefficient vector, and ${\epsilon }_i$ are positive random variables. Predictors $x_i\in {\mathbb{R}}^p$ may be random and serially dependent, while we often set $x_{i1}=1$, that is, incorporate the intercept in the exponent. The parameter space $\mathcal{B}\subset {\mathbb{R}}^p$ of $\mathbf{\beta }$ is a bounded convex domain such that ${\mathbf{\beta }}_0\in \mathcal{B}$. We implicitly assume that the model is correctly specified, so that there exists a true parameter ${\mathbf{\beta }}_0=({\beta }_{1,0},\cdots ,{\beta }_{p,0})\in \mathcal{B}$. We want to estimate ${\mathbf{\beta }}_0$ from a sample $\{(x_i,y_i)$, $i=1,\cdots ,n\}$.

We first remark that the condition $x_{i1}=1$ ensures that the model (1) is scale-free regarding variables ${\epsilon }_i$, which is an essentially different nature from the linear regression model $y_i=x_i^T\mathbf{\beta }+{\epsilon }_i$. Specifically, multiplying a positive constant $\sigma$ to ${\epsilon }_i$ results in the translation of the intercept in the exponent: so that the change from ${\epsilon }_i$ to $\sigma {\epsilon }_i$ is equivalent to that from ${\beta }_1$ to ${\beta }_1+log\sigma$. See Remark 1 on the distribution of ${\epsilon }_1$.

To provide a simple expression of the loss functions based on the relative error, we write

Chen et al. [6,8] pointed out that the loss criterion for relative error may depend on $|(y_i-t_i)/y_i|$ and / or $|(y_i-t_i)/t_i|$. These authors also proposed general relative error (GRE) criteria, defined as where $g:[0,\infty )\times [0,\infty )\to [0,\infty )$. Most of the loss functions based on the relative error are included in the GRE. Park and Stefanski [2] considered a loss function $g(a,b)=a^2$. It may highly depend on a small $y_i$ because it includes $1/y_i^2$ terms, and then the estimator can be numerically unstable. Consistency and asymptotic normality may not be established under general regularity conditions [8]. The loss functions based on $g(a,b)=max\{a,b\}$ [3] and $g(a,b)=a+b$ (least absolute relative error estimation, [6]) can have desirable asymptotic properties [3,6]. However, the minimization of the loss function can be challenging, in particular for high-dimensional data, when the function is nonsmooth or nonconvex.

In practice, the following two criteria would be useful:

The loss functions of LPRE and LSRE are smooth and convex, and also possess desirable asymptotic properties [8]. The above-described GRE criteria and their properties are summarized in Table 1. Particularly, the “convexity” in the case of $g(a,b)=a+b$ holds when ${\epsilon }_i>0$, ${\epsilon }_i\ne 1$, and ${\sum }_{i=1}^n{x}_i{x}_i^T$ is positive definite, since the Hessian matrix of the corresponding $G\left(\mathbf{\beta }\right)$ is ${\sum }_{i=1}^n|{\epsilon }_i-{\epsilon }_i^{-1}|x_i{x}_i^T$ a.s.

Although not essential, we assume that the variables ${\epsilon }_i$ in Equation (1) are i.i.d. with common density function h. As in Chen et al. [8], we consider the following class of h associated with g: where and $C\left(g\right)$ is a normalizing constant ($\int h\left(\epsilon \right)d\epsilon =1$) and $I_{+}$ denotes the indicator function of set $(0,\infty )$. Furthermore, we assume the symmetry property $g(a,b)=g(b,a)$, $a,b\ge 0$, from which it follows that ${\epsilon }_1\sim {\epsilon }_1^{-1}$. The latter property is necessary for a score function to be associated with the gradient of a GRE loss function, hence being a martingale with respect to a suitable filtration, which often entails estimation efficiency. Indeed, the asymmetry of $g(a,b)$ (i.e., $g(a,b)\ne g(b,a)$) may produce a substantial bias in the estimation [3]. The entire set of our regularity conditions will be shown in Section 4.3. The conditions therein concerning g are easily verified for both LPRE and LSRE.

In this paper, we implicitly suppose that $i=1,\cdots ,n$ denote “time” indices. As usual, in order to deal with cases of non-random and random predictors in a unified manner, we employ the partial-likelihood framework. Specifically, in the expression of the joint density (with the obvious notation for the densities) we ignore the first product $\{\cdots \}$ and only look at the second one $\{\cdots \}$, which is defined as the partial likelihood. We further assume that the ith-stage noise ${\epsilon }_i$ is independent of $(x_1,\cdots ,x_i,y_1,\cdots ,y_{i-1})$, so that, in view of Equation (1), we have

The density function of response y given $x_i$ is

From Equation (3), we see that the maximum likelihood estimator (MLE) based on the error distribution in Equation (5) is obtained by the minimization of Equation (2). For example, the density functions of LPRE and LSRE are where $K_{\nu }\left(z\right)$ denotes the modified Bessel function of third kind with index $\nu \in \mathbb{R}$: and $C_{LSRE}$ is a constant term. Constant terms are numerically computed as $K_0\left(2\right)\approx 0.1139$ and $C_{LSRE}\approx 0.911411$. Density (6) is a special case of the generalized inverse Gaussian distribution (see, e.g., [20]).

In practice, outliers can often be observed. For example, the electricity consumption data can have the outliers on extremely hot days. The estimation methods via GRE criteria, including LPRE and LSRE, are not robust against outliers, because the corresponding density functions are not generally heavy-tailed. Therefore, a relative error estimation method that is robust against the outliers is needed. To achieve this, we consider minimizing the negative $\gamma$-(partial-)likelihood function based on the $\gamma$-cross entropy [17].

We now define the negative $\gamma$-(partial-)likelihood function by where $\gamma >0$ is a parameter that controls the degrees of robustness; $\gamma \to 0$ corresponds to the negative log-likelihood function, and robustness is enhanced as $\gamma$ increases. On the other hand, a too large $\gamma$ can decrease the efficiency of the estimator [18]. In practice, the value of $\gamma$ may be selected by a cross-validation based on $\gamma$-cross entropy (see, e.g., [18,21]). We refer to Kawashima and Fujisawa [22] for more recent observations on comparison of the $\gamma$-divergences between Fujisawa and Eguchi [18] and Kawashima and Fujisawa [17].

There are several likelihood functions which yield robust estimation. Examples include the $L_q$-likelihood [23], and the likelihood based on the density power divergence [24], referred to as $\beta$-likelihood. It is shown that the $\gamma$-likelihood, the $L_q$-likelihood, and the $\beta$-likelihood are closely related. The negative $\beta$-likelihood function ${\ell }_{\alpha ,n}\left(\mathbf{\beta }\right)$ and the negative $L_q$-likelihood function ${\ell }_{q,n}\left(\mathbf{\beta }\right)$ are, respectively, expressed as

The difference between $\gamma$-likelihood and $\beta$-likelihood is just the existence of the logarithm on ${\ell }_{\gamma ,n}\left(\mathbf{\beta }\right)$. Furthermore, substituting $q=1-\alpha$ into Equation (9) gives us

Therefore, the minimization of the negative $L_q$-likelihood function is equivalent to minimization of the negative $\beta$-likelihood function without second term in the right side of Equation (8). Note that the $\gamma$-likelihood has the redescending property, a desirable property in robust statistics literature, as shown in Section 4.2. Moreover, it is known that the $\gamma$-likelihood is the essentially unique divergence that is robust against heavy contamination (see [18] for details). On the other hand, we have not shown whether the $L_q$-likelihood and/or the $\beta$-likelihood have the redescending property or not.

The integration $\int f{\left(y\right|x_i;\mathbf{\beta })}^{1+\gamma }dy$ in the second term on the right-hand side of Equation (7) is where is a constant term, which is assumed to be finite. Then, Equation (7) is expressed as where $C_0(\gamma ,h)$ is a constant term free from $\mathbf{\beta }$. We define the maximum $\gamma$-likelihood estimator to be any element such that

Even if the GRE criterion in Equation (2) is a convex function, the negative $\gamma$-likelihood function is nonconvex. Therefore, it is difficult to find a global minimum. Here, we derive the MM (majorize-minimization) algorithm to obtain a local minimum. The MM algorithm monotonically decreases the objective function at each iteration. We refer to Hunter and Lange [27] for a concise account of the MM algorithm.

Let ${\mathbf{\beta }}^{\left(t\right)}$ be the value of the parameter at the tth iteration. The negative $\gamma$-likelihood function in Equation (11) consists of two nonconvex functions, ${\ell }_1\left(\mathbf{\beta }\right)$ and ${\ell }_2\left(\mathbf{\beta }\right)$. The majorization functions of ${\ell }_j\left(\mathbf{\beta }\right)$, say ${\tilde{\ell }}_j\left(\mathbf{\beta }\right|{\mathbf{\beta }}^{\left(t\right)})$ ($j=1,2$), are constructed so that the optimization of ${min}_{\mathbf{\beta }}{\tilde{\ell }}_j\left(\mathbf{\beta }\right|{\mathbf{\beta }}^{\left(t\right)})$ is much easier than that of ${min}_{\mathbf{\beta }}{\ell }_j\left(\mathbf{\beta }\right)$. The majorization functions must satisfy the following inequalities:

Here, we construct majorization functions ${\tilde{\ell }}_j\left(\mathbf{\beta }\right|{\mathbf{\beta }}^{\left(t\right)})$ for $j=1,2$.