# Information Complexity Criterion for Model Selection in Robust Regression
Using A New Robust Penalty Term

Esra Pamukçu –  (1) Fırat University, Faculty of Science, Department of
Statistics, Elazıǧ, Turkey<EMAIL_ADDRESS>
Mehmet Niyazi Çankaya –  Department of International Trading and Finance,
Faculty of Applied Sciences, Uşak University & Department of Statistics,
Faculty of Arts and Science, Uşak University, Uşak<EMAIL_ADDRESS>*Corresponding Author)

###### Abstract

Model selection is basically a process of finding the best model from the
subset of models in which the explanatory variables are effective on the
response variable. The log likelihood function for the lack of fit term and a
specified penalty term are used as two parts in a model selection criteria. In
this paper, we derive a new tool for the model selection in robust regression.
We introduce a new definition of relative entropy based on objective
functions. Due to the analytical simplicity, we use Huber’s objective function
$\rho_{H}$ and propose our specified penalty term $C_{0}^{\rho_{H}}$ to derive
new Information Complexity Criterion ($RICOMP_{C_{0}^{\rho_{H}}}$) as a robust
model selection tool. Additionally, by using the properties of
$C_{0}^{\rho_{H}}$, we propose a new value of tuning parameter called
$k_{C_{0}}$ for the Huber’s $\rho_{H}$. If a contamination to normal
distribution exists, $RICOMP_{C_{0}^{\rho_{H}}}$ chooses the true model better
than the rival ones. Monte Carlo Simulation studies are carried out to show
the utility both of $k_{C_{0}}$ and $RICOMP_{C_{0}^{\rho_{H}}}$. A real data
example is also given.

Key-Words:

  * •

Information theory, Model selection, Robust statistics.

##  1 Introduction

Since the early 1970s, it has been possible to come across many studies on the
model selection algorithm and information criteria. These include classical
methods and methods based on information criteria. Classical model selection
methods are generally performed by hypothesis testing. However, most
statisticians and some other scientists stress that the confidence interval
which is used in the classical model selection methods is baseless [23]. Some
scientists have stated that the hypothesis testing approach does not have a
theoretical accuracy and is generally insufficiently valid [6].

Model selection procedures based on information criteria (IC) provide an
alternative to classical approaches. The basic idea of such criteria is to
minimize the Kullback-Leibler distance between the true model distribution and
the distribution of response variable which belongs to most appropriate model.
[1] and [4] with the successive published articles has led to developments in
the field of statistical modeling and statistical model selection or
evaluation. Therefore, he is considered as a first researcher in this field.
Akaike Information Criterion (AIC) proposed by Akaike is given by

$AIC=-2\log L(v;\hat{\boldsymbol{\beta}},\hat{\sigma})+2p,$ (1.1)

where $\hat{\boldsymbol{\beta}}$ is a vector of parameters
$\beta_{1},\beta_{2},\cdots,\beta_{p}$ and $\sigma$ is scale parameter.
$\hat{\boldsymbol{\beta}}$ and $\hat{\sigma}$ are ML estimators of these
parameters of probability density (p.d.) function
$f(v;\boldsymbol{\beta},\sigma)$. $p$ is the number of estimable parameters.

In the following years, many model selection criteria based on AIC were
proposed. Some of these are Takeuchi’s Information Criteria also known as the
Generalized Akaike Information Criteria proposed by [33], Corrected
Information Criteria developed by [32], the Bayesian Information Criterion
also known as the Bayesian Information Criterion developed by [28], Hannan-
Quinn Information Criterion developed by [16], Consistent Akaike Information
Criterion developed by [7].

ICOMP (I; Information-COMP; Complexity) which are information criteria
developed by [8] is a tool for model selection in multivariate linear and
nonlinear models. Whereas AIC is only intended to strike a balance between the
lack of fit and the penalty term, ICOMP aims to establish this balance by
taking into account a complexity measure that measures how the parameters in
the model relate to each others. Therefore, instead of penalizing the number
of independent parameters directly, it penalizes the covariance complexity of
the model. In the most general form, ICOMP is given by

$ICOMP=-2\log L(v;\hat{\boldsymbol{\beta}},\hat{\sigma})+C(\Sigma).$ (1.2)

The lack of fit part $-2\log L$ remains unchanged while the penalty term
$C(\Sigma)$ can be chosen by user with different ways to increase the
performance of ICOMP.

ICOMP is defined as follow,

$ICOMP=-2\log
L(v;\hat{\boldsymbol{\beta}},\hat{\sigma})+2C_{1}(\Sigma(\hat{\boldsymbol{\beta}})),$
(1.3)

where $\Sigma$ is the estimated covariance matrix of the model and
$C_{1}(\Sigma)$ is the maximal entropic complexity defined by [8]. We will
provide clear explanations about the complexity $C_{1}$ (see section 2.2).

There are several forms and theoretical justifications of ICOMP defined as a
general form in Eq. (1.3). So far, ICOMP have been compared with other model
selection criteria on many different methods and it has been observed that
they give better results clearly [22, 2, 18, 21, 25].

In this paper, our main idea is to produce a new robust model selection
criterion when a contamination into normal distribution occurs. We follow
$C_{0}$ based on work of [12] and introduce our penalty term
$C_{0}^{\rho_{H}}$ by using Huber’s objective function ($\rho_{H}$) in the
framework of $C_{0}$. We choose Huber’s $\rho$ function because of its
analytical tractability for calculating the expectation $E(\cdot)$. The normal
distribution is chosen as the underlying distribution for calculating
$E(\cdot)$ in Shannon’s relative entropy or Shannon’s divergence.

The rest of paper is organized as follows. Section 2.1 introduces the robust
model selection criteria. Section 2.2 provides the theoretical baseline of the
information complexity criteria based on Shannon’s relative entropy. We
propose our penalty term $C_{0}^{\rho_{H}}$ in Section 2.3. Section 2.4
presents the proposition of the new value of Huber’s tuning parameter $k$. In
section 2.5, we propose the robust ICOMP based on Huber’s $\rho_{H}$ function.
The numerical examples are demonstrated in Section 3. Finally, Section 4 is
divided by conclusion and brief discussion.

##  2 Material and Method

###  2.1 Robust Model Selection

Over the past 50 years, many robust estimation techniques have been proposed
as alternative approaches to the sum of least squares technique for regression
modeling when there are small deviations from the model assumption about
underlying distribution and unusual observations in the data sets [5, 17, 15].

Whereas there are a lot of works about robust estimation methods for
regression modeling [19, 3], model selection in robust regression is neglected
in the literature. Model selection criteria which are based on Gaussian
assumptions discussed above are not robust since the log likelihood of normal
distribution is used as lack of fit term. First studies about the robustifed
version of AIC were proposed by [14] and [27]. For the penalty term, they use
a generalization of AIC proposed by [6] and given as the following form:

$AIC=-2\log L(v;\hat{\boldsymbol{\beta}},\hat{\sigma})+\kappa p.$ (2.1)

A choice for the parameter $\kappa p$ in Eq. (2.1) is given by [31] as follow:

$AIC=-2\log
L(v;\hat{\boldsymbol{\beta}},\hat{\sigma})+\text{trace}(L_{2}^{-1}(\hat{\boldsymbol{\beta}},\hat{\sigma})L_{1}(\hat{\boldsymbol{\beta}},\hat{\sigma})),$
(2.2)

where

$L_{2}(\boldsymbol{\theta})=E\left(\frac{\partial^{2}}{\partial\boldsymbol{\theta}\partial\boldsymbol{\theta}^{T}}\log
f(v;\boldsymbol{\theta})\right),$ (2.3)
$L_{1}(\boldsymbol{\theta})=E\left(\frac{\partial}{\partial\boldsymbol{\theta}}\log
f(v;\boldsymbol{\theta})\left(\frac{\partial}{\partial\boldsymbol{\theta}}\log
f(v;\boldsymbol{\theta})\right)^{T}\right).$ (2.4)

$L_{2}$ and $L_{1}$ are also known as the definitions of Fisher Information
Matrix (FIM). Let us introduce robust versions of Eqs. (2.3) and (2.4) using
the objective function $\rho$ instead of $\log(f)$, as considered by [14] and
[27].

$L_{2}(\boldsymbol{\theta})=E\left(\frac{\partial^{2}}{\partial\boldsymbol{\theta}\partial\boldsymbol{\theta}^{T}}\rho(v;\boldsymbol{\theta})\right),$
(2.5)
$L_{1}(\boldsymbol{\theta})=E\left(\frac{\partial}{\partial\boldsymbol{\theta}}\rho(v;\boldsymbol{\theta})\left(\frac{\partial}{\partial\boldsymbol{\theta}}\rho(v;\boldsymbol{\theta})\right)^{T}\right).$
(2.6)

Due to fact that
$\psi(v;\boldsymbol{\theta})=\frac{\partial}{\partial\boldsymbol{\theta}}\rho(v;\boldsymbol{\theta})$,
we can rewrite Eqs. (2.5) and (2.6) as follows:

$L_{2}(\boldsymbol{\theta})=E\left(\frac{\partial}{\partial\boldsymbol{\theta}}\psi(v;\boldsymbol{\theta})\right),$
(2.7)
$L_{1}(\boldsymbol{\theta})=E(\psi(v;\boldsymbol{\theta})\psi(v;\boldsymbol{\theta})^{T}).$
(2.8)

As a consequent, the robust penalty term for AIC is proposed by [14]:

$\footnotesize\alpha_{H}=\alpha_{R}+[E(\frac{\partial}{\partial\boldsymbol{\theta}}\psi(v;\boldsymbol{\theta}))]^{-1}E(\psi(v;\boldsymbol{\theta})\psi(v;\boldsymbol{\theta})^{T})[E(\frac{\partial}{\partial\boldsymbol{\theta}}\psi(v;\boldsymbol{\theta}))]^{-1}p.$
(2.9)

and alternative proposition, $\alpha_{R}$, is introduced by [27]

$\alpha_{R}=[E(\frac{\partial}{\partial\boldsymbol{\theta}}\psi(v;\boldsymbol{\theta}))]^{-1}E(\psi(v;\boldsymbol{\theta})\psi(v;\boldsymbol{\theta})^{T})p.$
(2.10)

For the lack of fit part, they use objective function which defines robust
estimators instead of using the likelihood function in AIC to reach the
robustness. Therefore, $AIC_{\text{Hampel}}$ (AICH) and
$AIC_{\text{Ronchetti}}$ (AICR) can be defined as the following forms:

$AIC_{\text{H}}=2\sum_{i=1}^{n}\rho(v;\boldsymbol{\theta})+\alpha_{H},$ (2.11)
$AIC_{\text{R}}=2\sum_{i=1}^{n}\rho(v;\boldsymbol{\theta})+\alpha_{R}.$ (2.12)

When the objective function $\rho$ is chosen as Huber’s $\rho_{H}$, the
derived forms of Eqs. (2.9) and (2.10) are given in Appendix A.

The robust version of ICOMP called as RICOMP(IFIM) defined by [24] is given as
follow

$RICOMP(IFIM)=-2\log
L(\hat{\boldsymbol{\theta}_{R}})+2C_{1}(\mathcal{F}_{R}^{-1}),$ (2.13)

where
$\hat{\boldsymbol{\theta}}=(\hat{\boldsymbol{\beta}}_{R},\hat{\sigma}_{R}^{2})$
represents the robust estimate of the parameter vector which can be obtained
from the M-estimation and $\mathcal{\hat{F}}_{R}^{-1}$ is Inverse of Fisher
Information Matrix (IFIM) is given by

$\mathcal{\hat{F}}_{R}^{-1}=\widehat{Cov}(\hat{\boldsymbol{\beta}_{R}},\hat{\sigma}_{R}^{2})=\begin{bmatrix}\widehat{Cov}(\hat{\boldsymbol{\beta}_{R}})&0\\\
0&2\hat{\sigma}_{R}^{2}\\\ \end{bmatrix}.$ (2.14)

So, $C_{1}(\mathcal{\hat{F}}_{R}^{-1})$ is computed by

$C_{1}(\mathcal{\hat{F}}_{R}^{-1})=\frac{s}{2}\log\left[s^{-1}tr(\hat{\mathcal{F}}_{R}^{-1})\right]-\frac{1}{2}\log|\hat{\mathcal{F}}_{R}^{-1}|.$
(2.15)

For the robustness of lack of fit part in ICOMP, [24] used the Least Favorable
distribution instead of normal distribution in $\log L$ and proposed
$RICOMP(IFIM)$ as follow:

$\displaystyle RICOMP(IFIM)$ $\displaystyle=$ $\displaystyle
n\log(2\pi)+2n\log(\sigma)-2\log w(x_{i})$ $\displaystyle+$ $\displaystyle
2\sum_{i=1}^{n}\rho\left(\frac{w(x_{i})}{\hat{\sigma}}(y_{i}-x_{i}^{T}\hat{\boldsymbol{\beta}})\right)+2C_{1}(\mathcal{\hat{F}}_{R}^{-1}).$

For the special case of the function $w(x_{i})=1$, Eq. (2.1) drops to Huber’s
$\rho$ function [17].

Another definition of robust ICOMP called as RICOMPM from M-estimation method
is defined by [13]

$RICOMP_{M}=2\sum_{i=1}^{n}\rho(v_{i};\hat{\boldsymbol{\beta}_{M}},\hat{\sigma_{M}})+2C(\Sigma_{M}),$
(2.17)

where $C$ is complexity part of ICOMP. In this paper, we use $\rho$ as Huber’s
$\rho_{H}$ function and $C$ as maximal complexity of Bozdogan’s $C_{1}$.

###  2.2 Information theoretic model selection criterion based on Shannon’s
relative entropy

Let us consider a continuous d-dimensional distribution with joint density
function $f(v_{1},v_{2},\cdots,v_{d})$ and marginal density functions
$f_{j}(v_{j})$, $j=1,2,\cdots,d$.

Shannon’s marginal entropy is given by

$H(X_{j})=-E[\log f_{j}(v_{j})]=-\int f_{j}(v_{j})\log
f(v_{j})\mathrm{d}v_{j}$ (2.18)

and joint or global entropy is given by [29]

$\displaystyle H(v_{1},v_{2},\cdots,v_{d})$ $\displaystyle=$
$\displaystyle-E[\log f(v_{1},v_{2},\cdots,v_{d})]$ $\displaystyle=$
$\displaystyle-\int\cdots\int f(v_{1},v_{2},\cdots,v_{d})\log
f(v_{1},v_{2},\cdots,v_{d})\mathrm{d}v_{1}\mathrm{d}v_{2}\cdots\mathrm{d}v_{d}.$

This is also called the Shannon’s complexity [26]. Kullback-Leibler is used
Eqs. (2.18) and (2.2) to define information divergence against independence.
It is given by [20]

$KL(f_{j}(v_{j})||f(v_{1},v_{2},\cdots,v_{d}))=\sum_{j=1}^{d}H(v_{j})-H(v_{1},v_{2},\cdots,v_{d}).$
(2.20)

This is also known as the mutual information. It was used to introduce
information theoretic measure of covariance complexity by using Gaussian
distribution by [12]. The initial definition of covariance complexity is

$C_{0}(\sigma_{jj},\Sigma)=\frac{1}{2}\sum_{j=1}^{p}\log(\sigma_{jj})-\frac{1}{2}\log|\Sigma|,$
(2.21)

where $\sigma_{jj}$= $\sigma_{j}^{2}$ is the j-th diagonal element of $\Sigma$
and $p$ is the dimension of $\Sigma$. Another definition of complexity which
is obtained by maximization of $C_{0}$ according to set of all orthonormal
transformations was proposed by [8]:

$C_{1}(\Sigma)=\frac{p}{2}\log\left[\frac{tr(\Sigma)}{p}\right]-\frac{1}{2}\log|\Sigma|.$
(2.22)

Since the complexities in Eqs. (2.21) and (2.22) are based on the multivariate
Gaussian distributional assumption, we can say that they are not robust. Thus,
we expect that they do not work properly in a Robust Model Selection Criteria
(RMSC). In this paper, we introduce a new robust complexity term based on the
structure of complexity in Eq. (2.21) by using Huber’s $\rho_{H}$ objective
function. Next section represents the theoretical details.

###  2.3 Huber’s $\rho_{H}$ complexity based on KL divergence for the normal
distribution

Any objective function $\rho$ can be used instead of using $-\log(f)$ in Eq.
(2.18) [17, 9]. We use this similarity in construction of our complexity term.
Let us redefine the Shannon’s entropy based on $\rho$.

$\displaystyle H(v)$ $\displaystyle=$
$\displaystyle\int-\log[f(v)]f(v)\mathrm{d}v$ (2.23) $\displaystyle
H^{\rho}(v)$ $\displaystyle=$ $\displaystyle\int\rho(v)f(v)\mathrm{d}v,$

where $f$ is an underlying distribution.

To obtain Huber’s $\rho_{H}$ complexity ( $C_{0}^{\rho_{H}}$) based on
Shannon’s entropy, $\rho$ and $f$ in Eq. (2.23) are chosen as $\rho_{H}$ and
normal distribution, respectively.

Firstly, we consider to get univariate case, that is,
$\sum_{j=1}^{p}E_{g}(\rho(u))=\int\rho(u)g(u)\mathrm{d}u$. It should be noted
that the underlying distribution $g$ for expectation $E(\cdot)$ is chosen as
standard normal distribution. Considering Huber’s $\rho_{H}$ given as,

$\rho_{H}(u)=\left\\{\begin{array}[]{ll}u^{2}\nobreak\leavevmode&,|u|\leq k\\\
2k|u|-k^{2}&,|u|>k,\end{array}\right.$ (2.24)

we have

$E(\rho_{H}(u))=\int_{-\infty}^{-k}\left(-ku-\frac{k^{2}}{2}\right)g(u)\mathrm{d}u+\int_{-k}^{k}u^{2}g(u)\mathrm{d}u+\int_{k}^{\infty}\left(ku-\frac{k^{2}}{2}\right)g(u)\mathrm{d}u.$
(2.25)

As is seen, there are three parts in Eq. (2.25). If we take,
$u=\frac{v-\mu}{\sigma}$. Then $\mathrm{d}u=\mathrm{d}v/\sigma$. Thus,

$\displaystyle E(\rho_{H}(v;\mu,\sigma))$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{-k}\left(-k(\frac{v-\mu}{\sigma})-\frac{k^{2}}{2}\right)f(v;\mu,\sigma)\mathrm{d}v$
$\displaystyle+$
$\displaystyle\int_{-k}^{k}\left(\frac{v-\mu}{\sigma}\right)^{2}f(v;\mu,\sigma)\mathrm{d}v$
$\displaystyle+$
$\displaystyle\int_{k}^{\infty}\left(k(\frac{v-\mu}{\sigma})-\frac{k^{2}}{2}\right)f(v;\mu,\sigma)\mathrm{d}v,$

where
$f(v;\mu,\sigma)=\frac{1}{\sigma\sqrt{2\pi}}\exp\\{-\left(\frac{v-\mu}{\sqrt{2}\sigma}\right)^{2}\\}$.

For calculating the integrals in Eq. (2.3), we use the integral kernels in
[10]. Then, we have

$\displaystyle E(\rho_{H}(v;\mu,\sigma))$ $\displaystyle=$
$\displaystyle\frac{1}{\sigma\sqrt{\pi}}\left[\gamma(\frac{3}{2},\frac{k^{2}}{2})-\frac{k^{2}}{2}\Gamma(\frac{1}{2},\frac{k^{2}}{2})\right],$
(2.27)

where $\gamma$ is lower and $\Gamma$ is upper incomplete gamma functions,
respectively. $k$ is the tuning parameter of $\rho_{H}$.

Let us consider the multivariate case, we can write

$E(\rho_{H}(v_{1},v_{2},\cdots,v_{n}))=\int\rho_{H}(\left((v-\boldsymbol{\mu})^{T}\Sigma^{-1}(v-\boldsymbol{\mu})\right)^{1/2})f(v;\boldsymbol{\mu},\Sigma)\mathrm{d}v,$
(2.28)

where the multivariate form
$\left((v-\boldsymbol{\mu})^{T}\Sigma^{-1}(v-\boldsymbol{\mu})\right)^{1/2}$
can be equivalent to $u=\frac{v-\mu}{\sigma}$ in the univariate case.

For the integral of middle part in Eq. (2.28), we have to calculate

$\int\frac{1}{2}(v-\boldsymbol{\mu})^{T}\Sigma^{-1}(v-\boldsymbol{\mu})f(v;\boldsymbol{\mu},\Sigma)\mathrm{d}v,$
(2.29)

where
$f(v;\boldsymbol{\mu},\Sigma)=(2\pi)^{-p/2}|\Sigma|^{-1/2}\exp\\{-\frac{1}{2}[(v-\boldsymbol{\mu})^{T}\Sigma^{-1}(v-\boldsymbol{\mu})]\\}$.

Let $u^{2}$ be $(v-\boldsymbol{\mu})^{T}\Sigma^{-1}(v-\boldsymbol{\mu})$.
After taking the derivatives,
$u\mathrm{d}u=(v-\boldsymbol{\mu})\Sigma^{-1}\mathrm{d}v$ and
$u^{-2}(v-\boldsymbol{\mu})^{T}\Sigma^{-1}u\mathrm{d}u=\mathrm{d}v$. With
little rearrangement, we have

$\Sigma(v-\boldsymbol{\mu})^{-1}u\mathrm{d}u=\mathrm{d}v.$ (2.30)

If we take $v-\boldsymbol{\mu}=z$, we can write $u^{2}=z^{T}\Sigma^{-1}z$.
Then $u\mathrm{d}u=z\Sigma^{-1}\mathrm{d}z$. Thus, we can rewrite Eq. (2.30)
as follow,

$\Sigma z^{-1}z\Sigma^{-1}\mathrm{d}z=\mathrm{d}v.$ (2.31)

By using the integral kernels in Ref. [10], we obtain following result for the
integral of middle part in Eq. (2.28) given by

$\frac{1}{|\Sigma|^{1/2}(2\pi)^{p/2}}\sqrt{2}\gamma(\frac{3}{2},\frac{k^{2}}{2}).$
(2.32)

For the tail part of $\rho_{H}$ function, we have to calculate following
integral.

$\int(k|(v-\boldsymbol{\mu})^{T}\Sigma^{-1}(v-\boldsymbol{\mu})|^{1/2}-\frac{k^{2}}{2})f(v;\boldsymbol{\mu},\Sigma)\mathrm{d}v.$
(2.33)

The absolute term
$|(v-\boldsymbol{\mu})^{T}\Sigma^{-1}(v-\boldsymbol{\mu})|^{1/2}$ will produce
the negative value if $v<\boldsymbol{\mu}$. Otherwise, it will be positive. So
the result of left of integral is

$\int
k|(v-\boldsymbol{\mu})^{T}\Sigma^{-1}(v-\boldsymbol{\mu})|^{1/2}g(v;\boldsymbol{\mu},\Sigma)\mathrm{d}v=0,$
(2.34)

and right part of integral is

$-\int\frac{k^{2}}{2}g(v;\boldsymbol{\mu},\Sigma)\mathrm{d}v=-\frac{1}{|\Sigma|^{1/2}(2\pi)^{p/2}}\frac{k^{2}}{\sqrt{2}}\Gamma(\frac{1}{2},\frac{k^{2}}{2}).$
(2.35)

Note that the variable transformation also includes to change the limit values
of integrals. So we omit to write the limits of integrals.

Finally, we define our new robust penalty term as follow:

$\displaystyle C_{0}^{\rho_{H}}(\Sigma(\hat{\boldsymbol{\beta}}),k)$
$\displaystyle=$
$\displaystyle\sum_{j=1}^{d}\frac{1}{\sigma_{jj}\sqrt{\pi}}\left[\gamma(\frac{3}{2},\frac{k^{2}}{2})-\frac{k^{2}}{2}\Gamma(\frac{1}{2},\frac{k^{2}}{2})\right]$
$\displaystyle-\frac{1}{|\Sigma|^{1/2}(2\pi)^{p/2}}\left[\sqrt{2}\gamma(\frac{3}{2},\frac{k^{2}}{2})-\frac{k^{2}}{\sqrt{2}}\Gamma(\frac{1}{2},\frac{k^{2}}{2})\right].$

The first part in Eq. (2.3) is sum of the marginal Shannon’s entropy based on
$\rho_{H}$ and the second part is joint Shannon’s entropy based on $\rho_{H}$.

###  2.4 A new proposition for the values of Huber’s tuning parameter using
$C_{0}^{\rho_{H}}$

It is important that $C_{0}^{\rho_{H}}$ depends on Huber’s tuning parameter
$k$. As is well-known, the value of tuning parameter $k$ is equal to $1.345$
in literature of robustness [5, 17]. However, for $\Sigma=I_{(p\times p)}$, if
we take $k=1.345$, we obtain $C_{0}^{\rho_{H}}\neq 0$ ($0.632784$). In this
case is against to the definition of complexity of a covariance matrix.
Additionally, since $C_{0}^{\rho_{H}}$ depends on the incomplete gamma
functions, it does not equal to zero exactly. So we have to find which value
of $k$ will give us $C_{0}^{\rho_{H}}(I_{(p\times p)})\approx 0$. By using the
simple grid search, we find $k=0.8875916$ Finally, $C_{0}^{\rho_{H}}$ has
following properties:

  * •

$C_{0}^{\rho_{H}}(I_{(p\times p)})\approx 0$, if $\Sigma=I_{(p\times p)}$ and
$k=0.8875916$.

  * •

$C_{0}^{\rho_{H}}(\Sigma)\rightarrow\infty$, if $|\Sigma|\rightarrow 0$ for
all values of $k$.

This is an important finding and we call this new value of $k$ as $k_{C_{0}}$.
It will be discussed how to obtain good estimates of parameters
$\boldsymbol{\beta}$ using $k_{C_{0}}$ in the simulation section.

###  2.5 A new robust information complexity criterion:
$RICOMP_{C_{0}^{\rho_{H}}}$

In this section, we define our new robust information complexity criterion
($RICOMP_{C_{0}^{\rho_{H}}}$). The robust penalty term, that is the second
component of $RICOMP_{C_{0}^{\rho_{H}}}$, is obtained in Section 2.3. For the
robustness of first component of $RICOMP_{C_{0}^{\rho_{H}}}$, that is lack of
fit term, we use the objective function instead of $-\log L$ in the lack of
fit part by inspiring from [14, 27, 24]. Therefore, our new criterion is
defined as follow

$RICOMP_{C_{0}^{\rho_{H}}}=n\log(2\pi)+n\log(\hat{\sigma}^{2})+2\sum_{i=1}^{n}\rho_{H}(v_{i};\hat{\boldsymbol{\beta}},\hat{\sigma},k)+2C_{0}^{\rho_{H}}(\Sigma).$
(2.37)

Eq. (2.37) has four components. These are $n\log(2\pi)$, $n\log(\sigma^{2})$,
$\rho_{H}(v_{i};\hat{\boldsymbol{\beta}},\hat{\sigma},k)$ and
$C_{0}^{\rho_{H}}(\Sigma)$. $n\log(2\pi)$, $n\log(\sigma^{2})$ are from
$RICOMP$ in Eq. (2.1) constructed based on the Least Favorable distribution to
add scale parameter into model selection. $n\log(\sigma^{2})$ is also a part
for lack of fit. Another part for lack of fit as a main part is represented by
$\rho_{H}(v_{i};\hat{\boldsymbol{\beta}},\hat{\sigma},k)$ to perform robust
modeling. The last component is robust penalty term based on the Huber’s
objective function. Note that since we use $\rho_{H}$ for lack of fit part,
the penalty term should be based on $\rho_{H}$.

##  3 Numerical examples

###  3.1 Monte Carlo simulation study for $k_{C_{0}}$

We discuss how $k_{C_{0}}$ will affect the robust estimation of parameters in
the robust regression.

We show the simulation protocol to generate the explanatory variables
$x=(x_{1},x_{2},\cdots,x_{5})$ which are dependent each other, as given by the
following design:

  * •

$x_{i}=\sqrt{1-\alpha^{2}}z_{i}+\alpha z_{6}$, $i=1,2,\cdots,5$

where the variables $z_{1},z_{2},\cdots,z_{6}$ are distributed from Burr III
distribution with $c=2$ and $k=20$ as shape parameters. Thus, the variables
$x_{1},x_{2},\cdots x_{5}$ have big values due to the heavy-tailedness
property of BIII distribution [11]. In our case, $\alpha=0.9$.

Let us introduce how to generate random variables including outliers for the
dependent variable $y$ in the regression model.

  * •

We have two populations with the sample sizes $n_{1}$ and $n_{2}$ from
$N_{1}(\mu_{1}=0,\sigma_{1}=1)$ and $N_{2}(\mu_{2}=5,\sigma_{2}=10)$. We also
use other contamination as $N_{2}(0,50)$. We combine these artificial data
sets, i.e., $n=(n_{1},n_{2})$ with different percent of the Level of
Contamination (LC) such as $\%10,\%20,\%30$ and $\%50$, i.e.,
$(1-\tau)N_{1}(0,1)+\tau N_{2}(5,10)$ and $(1-\tau)N_{1}(0,1)+\tau
N_{2}(0,50)$. $\tau\in[0,1]$ is LC. The population from $N_{2}$ is regarded as
contamination into the population from $N_{1}$. These two contaminations
$N_{2}(\mu_{2}=5,\sigma_{2}=10)$ and $N_{2}(\mu_{2}=0,\sigma_{2}=50)$ are used
to generate outliers which are distant from underlying distribution
$N_{1}(\mu_{1}=0,\sigma_{1}=1)$. Thus, we test the robustness property of
$RICOMP_{C_{0}^{\rho_{H}}}$ in Eq. (2.37). Thus, we have the random numbers
for $\varepsilon$ with $n\times 1$.

  * •

We generate dependent random variable as a true model:
$y=x_{1}+x_{2}+x_{3}+\sigma\varepsilon$, where $\sigma=1$.

  * •

’robustfit’ module in MATLAB2015a$\copyright$ is used to estimate regression
and scale parameters robustly.

    * –

$\hat{Y_{1}}=X\hat{\boldsymbol{\beta}}_{R}$ for $k=1.345$

    * –

$\hat{Y_{2}}=X\hat{\boldsymbol{\beta}}_{R}$ for $k_{C_{0}}=0.8875916$.

    * –

In order to compare the prediction performance of $\hat{Y_{1}}$ and
$\hat{Y_{2}}$, we use mean absolute error (MAE), given by

$\text{MAE}_{j}=\frac{1}{r}\sum_{t=1}^{r}\left(\frac{1}{n}\sum_{i=1}^{n}|Y_{i}-\hat{Y}_{j_{i}}|\right)\leavevmode\nobreak\
\leavevmode\nobreak\ \text{for}\leavevmode\nobreak\ \leavevmode\nobreak\
j=1,2.$ (3.1)

where $r$ and $n$ are the numbers of replication and sample size,
respectively. Using MAE for testing the performance of prediction is better
than mean squared error, i.e.,
$\frac{1}{r}\sum_{t=1}^{r}\left(\frac{1}{n}\sum_{i=1}^{n}(Y_{i}-\hat{Y}_{j_{i}})^{2}\right)$
[30].

Table 1: Performance comparison of $k$ and $k_{C_{0}}$ in robust regression

|  |  |  | $(1-\tau)N_{1}(0,1)+\tau N_{2}(5,10)$ | $(1-\tau)N_{1}(0,1)+\tau N_{2}(0,50)$
---|---|---|---|---|---
n |  LC* ($\%$) | $n_{1}$ | $n_{2}$ |  MAE1 |  MAE2 |  MAE1 |  MAE2
n=30 | 10 | 27 | 3 | 11.020 | 10.861 | 13.707 | 13.521
20 | 24 | 6 | 11.334 | 11.186 | 16.802 | 16.572
30 | 21 | 9 | 11.673 | 11.518 | 19.794 | 19.504
50 | 15 | 15 | 12.315 | 12.184 | 25.595 | 25.277
n=50 | 10 | 45 | 5 | 11.312 | 11.165 | 14.111 | 13.950
20 | 40 | 10 | 11.822 | 11.674 | 17.227 | 17.038
30 | 35 | 15 | 12.113 | 11.977 | 20.379 | 20.166
50 | 25 | 25 | 12.912 | 12.793 | 26.507 | 26.270
n=100 | 10 | 90 | 10 | 11.653 | 11.512 | 14.359 | 14.207
20 | 80 | 20 | 12.034 | 11.396 | 17.675 | 17.513
30 | 70 | 30 | 12.474 | 12.340 | 20.799 | 20.627
50 | 50 | 50 | 13.274 | 13.163 | 27.131 | 26.961

*LC represents Level of Contamination. MAE1 and MAE2 represent MAE of $k=1.345$ and $k_{C_{0}}=0.8875916$, respectively. $\tau\in[0,1]$ is LC.
---

All computations were done on a MacBook Pro with 8 GB Memory, 3.8 GHz,
IntelCore 8th generation i5. The results of MAEj are given in Table 1 The
number of replication is $r=10000$. According to Table 1, it is clear that
robust regression models with $k_{C_{0}}$ has superior performance by yielding
less MAE. We propose that $k_{C_{0}}$ can be used as a new tuning parameter in
robust modeling. Therefore, we will use $k_{C_{0}}$ in the computation of
robust regression models in the next simulation in order to be able to compare
the model selection criteria.

###  3.2 Monte Carlo simulation study for $RICOMP_{C_{0}^{\rho_{H}}}$

In this section, we present Monte Carlo simulation study to show the
performance of $RICOMP_{C_{0}^{\rho_{H}}}$ in choosing true model by comparing
other robust information criteria. The simulation protocol for generating the
explanatory and dependent variables is same with previous simulation. We only
take $N_{2}(\mu_{2}=5,\sigma_{2}=10)$ for the population with sample size
$n_{2}$ as outliers.

After having the artificial data sets, we use following steps in order to
select the best model in competing models using RMSC.

  * •

All subset models constructed from the explanatory variables
$x_{1},x_{2},\cdots,x_{5}$ are used. There are $2^{5}-1=31$ models to compare
in competing subset models and we know that only subset of
$\\{x_{1},x_{2},x_{3}\\}$ is true model.

  * •

Robust regression are performed to these competing models.

  * •

For each regression models, all RMSC are computed and listed.

  * •

The smallest value in the list is chosen.

  * •

The models which correspond to the smallest model selection values are checked
whether or not they are equivalent to true model.

  * •

If they are the true model, then the criteria which select true model have
earned one point and this score is stored.

This process for all scenarios with different LC levels is repeated 10000
times. The results are given in Table 2.

Table 2: The number of choosing true model for the model selection criteria in
10000 runs

|  LC* ($\%$) | $n_{1}$ | $n_{2}$ | $RICOMP_{C_{0}^{\rho_{H}}}$ | $RICOMP(IFIM)$ | $RICOMP_{M}$ | $AIC_{\text{H}}$ | $AIC_{\text{R}}$
---|---|---|---|---|---|---|---|---
n=30 | 10 | 27 | 3 | 2798 | 2801 | 273 | 1262 | 195
20 | 24 | 6 | 2637 | 2117 | 214 | 1101 | 146
30 | 21 | 9 | 2885 | 1946 | 147 | 1395 | 79
50 | 15 | 15 | 2281 | 1604 | 107 | 1388 | 24
n=50 | 10 | 45 | 5 | 3082 | 3505 | 158 | 564 | 115
20 | 40 | 10 | 3108 | 2501 | 100 | 397 | 77
30 | 35 | 15 | 3076 | 2050 | 62 | 392 | 51
50 | 25 | 25 | 3665 | 2362 | 31 | 971 | 17
n=100 | 10 | 90 | 10 | 3243 | 4347 | 42 | 134 | 37
20 | 80 | 20 | 3242 | 2966 | 6 | 44 | 7
30 | 70 | 30 | 3198 | 2146 | 8 | 45 | 4
50 | 50 | 50 | 4535 | 3088 | 1 | 89 | 2

*LC represents Level of Contamination. $AIC_{\text{H}}:AIC_{\text{Hampel}}$, $AIC_{\text{R}}$: $AIC_{\text{Ronchetti}}$
---

  * •

The bold number represents the model selection criterion which has the highest
number of choosing the true model.

  * •

For all scenario, $RICOMP_{C_{0}^{\rho_{H}}}$ has the highest number of
choosing true model for the different LC except for $\%10$. It is observed
that $RICOMP_{C_{0}^{\rho_{H}}}$ and $RICOMP(IFIM)$ have close performance in
choosing the true model for LC$=\%10$. Due to fact that LC is low, it can be
noted that the empirical distribution of data set is near to normal. Since the
penalty term of $RICOMP(IFIM)$ is based on the normal distribution, this
result is expected.

  * •

Evaluating the result of $AIC_{\text{H}}$ in itself, it works well for $n=30$.
However, its performance decreases noticeably when the number of observations
increases.

  * •

It is generally observed that $RICOMP_{C_{0}^{\rho_{H}}}$ and $RICOMP(IFIM)$
show better performance when compared other criteria.

  * •

$RICOMP_{C_{0}^{\rho_{H}}}$ yields more stable result even if LC increases. It
can be said that $C_{0}^{\rho_{H}}$ has important role in robust model
selection when the outliers exist in data sets.

  * •

For all scenario, $RICOMP_{M}$, $AIC_{\text{H}}$ and $AIC_{\text{R}}$ have
poor performance when compared others.

###  3.3 Real data example

We present the performance of $RICOMP_{C_{0}^{\rho_{H}}}$ and the other model
selection criteria to select the explanatory variables which can affect to the
response variable. For this purpose, we use the bridge construction data in
[34]. The variables in data are as follows:

  * •

Time: Time of construction ($y$)

  * •

CCost: Cost of construction ($x_{1}$)

  * •

Dwgs: Structural drawings number ($x_{2}$)

  * •

Length: Length of bridge ($x_{3}$)

  * •

Spans: Spans number ($x_{4}$)

  * •

DArea: Deck of area of bridge ($x_{5}$)

The investigation about relationship among variables of data set was performed
by [34] and illustrated the scatter plot of data. [34] noted that there exists
the striking nonlinear pattern among the variables, some of variables seem to
be skewed and there are multicollinearity among some explanatory variables.
After they applied log transformation on data set, they estimated regression
parameters by using the Modified Least Trimmed Square Counter
Multicollinearity (MLTSCM) as a robust regression method. From this point, the
regression techniques used by us and them for estimating regression parameters
are similar. These techniques can be comparable in the framework of efficient
and robust modeling on the dependent variable $y$. However, applying the
logarithmic transformation on the variables of a real data set may change the
information. So it may not a good way to regularize the variables for
modeling. If it is possible, modeling by using raw data without transformation
on the variables may yield better result because type of transformation is not
always one-to-one. In our case, we do not perform transformation and we
analyze the data for all subset competing models by using robust regression
with Huber M-estimator where $k=k_{C_{0}}=0.8875916$. For all competing
models, RMSC are computed and assigned the best model according to RMSC’s
smallest values. The results are in Table 3.

Table 3: The best subset of predictors of regression model using RMSC for the
bridge construction data set.

Variables | RMSC | Regression Model | MAE
---|---|---|---
$\\{x_{1}\\}$ |  $AIC_{H}:$ 3125.4 | $y=65.5715+0.2658x_{1}$ | 48.9173
$\\{x_{4}\\}$ |  $AIC_{R}:$ 31.205 | $y=37.5115+45.8683x_{4}$ | 46.164
$\\{x_{4}\\}$ |  $RICOMP_{M}:$23.1171
$\\{x_{1},x_{4}\\}$ |  $RICOMP(IFIM):$ 439.691 | $y=45.0571+0.1664x_{1}+22.7530x_{4}$ | 44.9872
$\\{x_{2},x_{5}\\}$ |  $RICOMP_{C_{0}^{\rho_{H}}}:$ 415.2527 | $y=-13.0449+17.7851x_{2}+2.7736x_{5}$ | 36.5115
*$\\{x_{1},x_{2},x_{4}\\}$ | - | $\log(y)=2.0304+0.3056\log(x_{1})+$ | 36.146
$0.6210\log(x_{2})+0.0657\log(x_{4})$

RMSC: Robust Model Selection Criteria. *Ref. [34]
---

According to Table 3, it is observed the best subsets of predictors are
$S_{1}=\\{x_{2},x_{5}\\}$ and $S_{2}=\\{x_{1},x_{2},x_{4}\\}$. Note that the
subset of $S_{1}$ is obtained by using $RICOMP_{C_{0}^{\rho_{H}}}$ and the
subset of $S_{2}$ is provided by Ref. [34].

As it is seen, MAE values are very close for $S_{1}$ and $S_{2}$. However, it
is important finding to obtain similar performance with fewer variables. This
shows utility and effectiveness of our proposed method with
$RICOMP_{C_{0}^{\rho_{H}}}$.

Let us give a few comment about performances of other RMSC. $RICOMP(IFIM)$ can
be considered to have better performance than others. The worst performance is
observed in $AIC_{H}$. We consider that it is due to penalty term in which has
a part $x^{T}x$ (see Appendix A). As is well-known, when the multicollinearity
problem occurs in the data, $x^{T}x$ could be singular.

(a) The true and predicted values for $y$

(b) Probability plot of residuals from Huber’s robust regression

(c) Probability plot of residuals obtained by Ref. [34]

Figure 1: (Color online) True and predicted values ($\hat{Y}$) from regression
models and probability plots of residuals

Fig. 1 (a) shows three lines to observe the nearness of prediction values
obtained from regression models in last two lines in Table 3 to true values of
$y$. Note that the values of MAE in Table 3 also show the prediction
performance of regression models. Fig. 1 (b) and (c) show the goodness of fit
by means of the probability plot. The dashed line represents the probability
values of normal distribution as the theoretical values. The blue pluses in
Fig 1 (b) and (c) represent the empirical probabilities. When we compare the
pluses in Fig. 1 (b) and (c), it is observed that the pluses in Fig. 1 (b) can
fit the theoretical values better than that in Fig. 1 (c).

##  4 Conclusion and discussion

In this paper, we presented several forms of the robust model selection
criteria. These are $AIC_{H}$, $AIC_{R}$, $RICOMP(IFIM)$ and $RICOMP_{M}$. We
have extended information theoretic measure of complexity ICOMP class of
criteria for robust regression in order to be resistant to high amount of
outliers in data set.

We used Huber’s $\rho_{H}$ to produce our robust penalty term for ICOMP and we
introduced ${C_{0}^{\rho_{H}}}$. Then, we used as a specific form of the lack
of fit part of $RICOMP(IFIM)$ by using Huber’s objective function. Finally, we
proposed a new robust information complexity criterion,
$RICOMP_{C_{0}^{\rho_{H}}}$. Additionaly, by using the properties of
${C_{0}^{\rho_{H}}}$, we proposed a new value of tuning parameter $k$ of
Huber’s $\rho_{H}$ as $0.8875916$. We performed a Monte Carlo simulation to
show important role of the value of $k$ in estimating regression parameters
robustly. The simulation results demonstrated clearly that our new value of
$k$ called $k_{C_{0}}$ can be used as a tuning parameter in robust regression
modeling. Moreover, $RICOMP_{C_{0}^{\rho_{H}}}$ was demonstrated on both Monte
Carlo simulation study and real data to show the utility and effectiveness.
According to numerical examples, the proposed criteria
$RICOMP_{C_{0}^{\rho_{H}}}$ is successful to select the best model in robust
regression for the data set which has high amount of outliers. In the future,
$RICOMP_{C_{1}^{\rho_{H}}}$ will be derived from [8].

## ACKNOWLEDGMENTS

We would like to thank referees, Associative Editor(s) and Editor. We would
like to say our special thanks to Prof. Dr. Hamparsum Bozdogan for his
constructive comments.

## References

  * [1] Akaike H (1973) Maximum likelihood identification of Gaussian autoregressive moving average models. Biometrika, 60(2), 255-265
  * [2] Akbilgic O, Bozdogan H (2011) Predictive subset selection using regression trees and RBF neural networks hybridized with the genetic algorithm. European Journal of Pure and Applied Mathematics, 4(4), 467-485
  * [3] Alfò M, Salvati N, Ranallli MG (2017) Finite mixtures of quantile and M-quantile regression models. Statistics and Computing, 27(2), 547-570
  * [4] Akaike H (1974) A new look at the statistical model identification. In Selected Papers of Hirotugu Akaike (pp. 215-222). Springer, New York, NY
  * [5] Andrews DF, Bickel PJ, Hampel FR, Huber PJ, Rogers WH, Tukey JW (1972) Robust Estimates of Location: Survey and Advances, Princeton University, Princeton, N.J.
  * [6] Bhansali RJ, Downham DY (1977) ”Some properties of the order of an autoregressive model selected by a generalization of Akaike? s EPF criterion.” Biometrika 64.3 : 547-551
  * [7] Bozdogan H (1987) Model selection and Akaike’s information criterion (AIC): The general theory and its analytical extensions. Psychometrika 52.3 : 345-370
  * [8] Bozdogan H (1988) ICOMP: A new model selection criterion. In: Hans H. Bock (Ed.), Classification and Related Methods of Data Analysis, North-Holland, Amsterdam, April, 599-608
  * [9] Çankaya, MN (2015) ”M-Estimators with asymmetric influence function: Properties and their applications.” PhD diss., University of Ankara
  * [10] Çankaya M (2018) Asymmetric bimodal exponential power distribution on the real line. Entropy, 20(1), 23
  * [11] Çankaya M N, Yalçınkaya A, Altındaǧ Ö, Arslan O (2019) On the robustness of an epsilon skew extension for Burr III distribution on the real line. Computational Statistics, 1-27
  * [12] Emden V (1971) An analysis of complexity. Mathematiscxh Centrum
  * [13] Güney Y (2018) Robust Model Selection Using Information Complexity (ICOMP) Crtiterion for Symmetric and Skew Distributions. PhD diss., Ankara University
  * [14] Hampel FR (1983) Some aspects of model choice in robust statistics. In Proceedings of the 44th Session of the ISI (Vol. 2, pp. 767-771
  * [15] Hampel FR, Ronchetti EM, Rousseeuw PJ, Stahel WA (1986) Robust Statistics: The Approach Based on Influence Functions, Wiley Series in Probability and Statistics, New York, USA
  * [16] Hannan EJ, Quinn BG (1979) The determination of the order of an autoregression. Journal of the Royal Statistical Society: Series B (Methodological) 41(2) 190-195
  * [17] Huber PJ (1981) Robust Statistics, Wiley, New York
  * [18] Howe ED, Bozdogan H, Kiroglu G (2012) Performance of information complexity criteria in structural equation models with applications. European Journal of Pure and Applied Mathematics, 5(3), 282-301
  * [19] Jurečková J, Picek J, Schindler M (2019) Robust statistical methods with R. Chapman and Hall/CRC
  * [20] Kullback, S. (1968) Information Theory and Statistics. Dover, New York
  * [21] Koc EK, Bozdogan H. (2015) Model selection in multivariate adaptive regression splines (MARS) using information complexity as the fitness function. Machine Learning, 101(1-3), 35-58
  * [22] Kwon Y, Bozdogan H, Bensmail H (2008) Performance of model selection criteria in Bayesian threshold VAR (TVAR) models. Econometric Reviews, 28(1-3), 83-101
  * [23] Linhart H, Zucchini W (1986) Model selection. John Wiley and Sons
  * [24] Liu Y (2007) Robust and Misspecification Resistant Model Selection in Regression Models with Information Complexity and Genetic Algorithms
  * [25] Pamukcu E, Bozdogan H, Çalık S (2015) A novel hybrid dimension reduction technique for undersized high dimensional gene expression data sets using information complexity criterion for cancer classification. Computational and mathematical methods in medicine, 2015
  * [26] Rissanen J (1989) Stochastic complexity in statistical inquiry
  * [27] Ronchetti E (1985) Robust model selection in regression. Statistics & Probability Letters 3 21-23
  * [28] Schwarz G (1978) Estimating the dimension of a model. The annals of statistics, 6(2), 461-464
  * [29] Shannon CE (1948) A mathematical theory of communication. Bell Systems Technology Journal, 27, 379-423
  * [30] Shao J (2003) Mathematical Statistics. Springer Texts in Statistics
  * [31] Stone M (1977) An asymptotic equivalence of choice of model by cross-validation and Akaike’s criterion. Journal of the Royal Statistical Society: Series B (Methodological), 39(1), 44-47
  * [32] Sugiura N (1978) Further analysts of the data by akaike’s information criterion and the finite corrections: Further analysts of the data by akaike’s. Communications in Statistics-Theory and Methods 7.1 13-26
  * [33] Takeuchi K (1976) The distribution of information statistics and the criterion of goodness of fit of models. Mathematical Science 153 12-18
  * [34] Roozbeh M, Babaie-Kafaki S, Sadigh, AN (2018) A heuristic approach to combat multicollinearity in least trimmed squares regression analysis. Applied Mathematical Modelling, 57, 105-120

##  Appendix A Appendix: $F(\boldsymbol{\theta})$ and $R(\boldsymbol{\theta})$
of Huber’s $\rho_{H}$ function

Eq. (2.2) includes the penalty term which is
$\text{trace}(F^{-1}(\hat{\boldsymbol{\beta}},\hat{\sigma})R(\hat{\boldsymbol{\beta}},\hat{\sigma}))$.
We give the definition of $F$ and $R$ given by

$F(\boldsymbol{\theta})=E\left(\frac{\partial^{2}\rho(v;\boldsymbol{\theta})}{\partial\boldsymbol{\theta}\partial\boldsymbol{\theta}^{T}}\right)=E\left(\frac{\partial}{\partial\boldsymbol{\theta}}\psi(v;\boldsymbol{\theta})\right)$
(1.1)

and

$R(\boldsymbol{\theta})=E\left(\left(\frac{\partial\rho(v;\boldsymbol{\theta})}{\partial\boldsymbol{\theta}}\right)\left(\frac{\partial\rho(v;\boldsymbol{\theta})}{\partial\boldsymbol{\theta}}\right)^{T}\right)=E\left(\psi(v;\boldsymbol{\theta})\psi(v;\boldsymbol{\theta})^{T}\right),$
(1.2)

where $\boldsymbol{\theta}=(\boldsymbol{\beta},\sigma)$ [31, 27]. In our case,
we choose Huber’s $\rho_{H}$ function for two cases.

Firstly, let us give the derivatives of middle of Huber’s $\rho_{H}$ with
respect to (w.r.t.) parameters.

$\frac{\partial^{2}\rho_{H}}{\partial\boldsymbol{\beta}\partial\boldsymbol{\beta}^{T}}=\frac{x^{T}x}{\sigma^{2}},$
(1.3)
$\frac{\partial}{\partial\sigma}\frac{\partial\rho_{H}}{\partial\boldsymbol{\beta}}=\frac{\partial}{\partial\boldsymbol{\beta}}\frac{\partial\rho_{H}}{\partial\sigma}=\frac{2x}{\sigma^{2}}\left(\frac{y-x^{T}\boldsymbol{\beta}}{\sigma}\right),$
(1.4)
$\frac{\partial^{2}\rho_{H}}{\partial\sigma^{2}}=\frac{3}{\sigma^{4}}(y-x^{T}\boldsymbol{\beta})^{2}.$
(1.5)

Let us give the derivatives of tail part of Huber’s $\rho_{H}$ w.r.t.
parameters.

$\frac{\partial^{2}\rho_{H}}{\partial\boldsymbol{\beta}\partial\boldsymbol{\beta}^{T}}=0x^{T}x=\boldsymbol{0}_{p\times
p},$ (1.6)
$\frac{\partial}{\partial\sigma}\frac{\partial\rho_{H}}{\partial\boldsymbol{\beta}}=\frac{kx}{\sigma^{2}}\text{sign}(y-x^{T}\boldsymbol{\beta}),$
(1.7)
$\frac{\partial^{2}\rho_{H}}{\partial\sigma^{2}}=\frac{2k}{\sigma^{3}}|y-x^{T}\boldsymbol{\beta}|.$
(1.8)

Now we provide the $R$ case. For middle part, we have

$\frac{\partial\rho_{H}}{\partial\boldsymbol{\beta}}=\frac{-x}{\sigma^{2}}(y-x^{T}\boldsymbol{\beta}),$
(1.9)
$\frac{\partial\rho_{H}}{\partial\sigma}=\frac{-1}{\sigma^{3}}(y-x^{T}\boldsymbol{\beta})^{2}.$
(1.10)

We have products of the derivatives of $\rho_{H}$ w.r.t. parameters:

$\frac{\partial\rho_{H}}{\partial\boldsymbol{\beta}}\frac{\partial\rho_{H}}{\partial\boldsymbol{\beta}}=\frac{x^{T}x}{\sigma^{4}}(y-x^{T}\boldsymbol{\beta})^{2},$
(1.11)
$\frac{\partial\rho_{H}}{\partial\sigma}\frac{\partial\rho_{H}}{\partial\sigma}=\frac{1}{\sigma^{6}}(y-x^{T}\boldsymbol{\beta})^{4},$
(1.12)
$\frac{\partial\rho_{H}}{\partial\boldsymbol{\beta}}\frac{\partial\rho_{H}}{\partial\sigma}=\frac{x}{\sigma^{5}}(y-x^{T}\boldsymbol{\beta})^{3}.$
(1.13)

Similarly, one can obtain the derivatives of tail part of Huber’s $\rho_{H}$
function.

Secondly, we give expectation for the elements of matrices. We use the
prepared kernels from Ref. [10] to calculate the integrals of $F$ and $R$. The
underlying distribution of expectation is chosen as normal distribution.

$\displaystyle
E\left(\frac{\partial^{2}\rho_{H}}{\partial\boldsymbol{\beta}\partial\boldsymbol{\beta}^{T}}\right)$
$\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\left(\frac{\partial^{2}\rho_{H}}{\partial\boldsymbol{\beta}\partial\boldsymbol{\beta}^{T}}\right)f(v;\boldsymbol{\beta},\sigma)$
$\displaystyle=$
$\displaystyle\frac{x^{T}x}{\sigma^{2}\sqrt{\pi}}\gamma(1/2,k^{2}/2),$

where
$f(v;\boldsymbol{\beta},\sigma)=\frac{1}{\sigma\sqrt{2\pi}}\exp\\{-\left(\frac{y-x^{T}\boldsymbol{\beta}}{\sqrt{2}\sigma}\right)^{2}\\}$.
The undiagonal elements of matrix $F$ and $R$, that is,

$E\left(\frac{\partial}{\partial\boldsymbol{\beta}}\frac{\partial\rho_{H}}{\partial\sigma}\right)=E\left(\frac{\partial}{\partial\sigma}\frac{\partial\rho_{H}}{\partial\boldsymbol{\beta}}\right)=E\left(\frac{\partial\rho_{H}}{\partial\boldsymbol{\beta}}\frac{\partial\rho_{H}}{\partial\sigma}\right)=0.$
(1.15)

Note that regression and scale parameters are independent for the regression
and scale case in Eq. (1.15). One can show these results easily.

$E\left(\frac{\partial^{2}\rho_{H}}{\partial\sigma^{2}}\right)=\frac{6}{\sigma^{2}\sqrt{\pi}}\gamma(3/2,k^{2}/2)+\frac{2\sqrt{2}k}{\sigma^{2}\sqrt{\pi}}\Gamma(1,k^{2}/2).$
(1.16)

Eqs. (A)-(1.16) are for $F$ case. Now let us give the integral results for $R$
case.

$E\left(\frac{\partial\rho_{H}}{\partial\boldsymbol{\beta}}\frac{\partial\rho_{H}}{\partial\boldsymbol{\beta}}\right)=\frac{x^{T}x}{\sigma^{2}\sqrt{\pi}}\left[2\gamma(3/2,k^{2}/2)+k^{2}\Gamma(1/2,k^{2}/2)\right],$
(1.17)
$E\left(\frac{\partial\rho_{H}}{\partial\sigma}\frac{\partial\rho_{H}}{\partial\sigma}\right)=\frac{2}{\sigma^{2}\sqrt{\pi}}\left[2\gamma(5/2,k^{2}/2)+k^{2}\Gamma(3/2,k^{2}/2)\right].$
(1.18)

Eqs. (A)-(1.18) show the calculated form of Eq. (2.9). Let us give matrix
elements of Eq. (2.9):

$\displaystyle A_{11}$ $\displaystyle=$ $\displaystyle\boldsymbol{1}_{p\times
p}\left[2\gamma(3/2,k^{2}/2)+k^{2}\Gamma(1/2,k^{2}/2)\right]/\gamma(1/2,k^{2}/2)$
$\displaystyle+$
$\displaystyle\frac{2\gamma(3/2,k^{2}/2)+k^{2}\Gamma(1/2,k^{2}/2)}{\gamma^{2}(1/2,k^{2}/2)}(x^{T}x)^{-1}\sigma^{2}\sqrt{\pi},$
$A_{12}=A_{21}=0,$ (1.20) $\displaystyle A_{22}$ $\displaystyle=$
$\displaystyle\frac{2\gamma(5/2,k^{2}/2)+k^{2}\Gamma(3/2,k^{2}/2)}{3\gamma(3/2,k^{2}/2)+\sqrt{2}k\Gamma(1,k^{2}/2)}$
$\displaystyle\left(1+\frac{\sigma^{2}\sqrt{\pi}}{2\left[3\gamma(3/2,k^{2}/2)+\sqrt{2}k\Gamma(1,k^{2}/2)\right]}\right).$

Eq. (A) has $(x^{T}x)^{-1}$.

$A=\begin{bmatrix}A_{11}&&A_{12}\\\ A_{21}&&A_{22}\end{bmatrix}.$ (1.22)