# Wilcoxon-type Multivariate Cluster Elastic Net

Mayu Hiraishi
Graduate School of Culture and Information Science
Doshisha University
Kyoto, Japan
<EMAIL_ADDRESS>
& Kensuke Tanioka
Department of Biomedical Sciences and Informatics
Doshisha University
Kyoto, Japan
<EMAIL_ADDRESS>
& Hiroshi Yadohisa
Department of Culture and Information Science
Doshisha University
Kyoto, Japan
<EMAIL_ADDRESS>
Use footnote for providing further information about author (webpage,
alternative address)— _not_ for acknowledging funding agencies.

###### Abstract

We propose a method for high dimensional multivariate regression that is
robust to random error distributions that are heavy-tailed or contain
outliers, while preserving estimation accuracy in normal random error
distributions. We extend the Wilcoxon-type regression to a multivariate
regression model as a tuning-free approach to robustness. Furthermore, the
proposed method regularizes the $L1$ and $L2$ terms of the clustering based on
$k$-means, which is extended from the multivariate cluster elastic net. The
estimation of the regression coefficient and variable selection are produced
simultaneously. Moreover, considering the relationship among the correlation
of response variables through the clustering is expected to improve the
estimation performance. Numerical simulation demonstrates that our proposed
method overperformed the multivariate cluster method and other methods of
multiple regression in the case of heavy-tailed error distribution and
outliers. It also showed stability in normal error distribution. Finally, we
confirm the efficacy of our proposed method using a data example for the gene
associated with breast cancer.

_Keywords_ MM algorithm $\cdot$ multivariate regression $\cdot$ robust
statistics

##  1 Introduction

Numerous studies on multivariate regression analysis have proposed methods to
adequately demonstrate the efficiency of the estimation. As part of this
framework, it is expected that the accuracy of the estimation is enhanced by
considering correlations among the objective variables. Breiman and Friedman
(1997) proposed the Curds and Whey method, which predicts multivariate
response variables with an optimal linear combination of least squared
predictors. Rothman et al. (2010) proposed a multivariate regression with
covariance estimation, which estimates regression coefficients and inverse
error covariance simultaneously. To improve estimation, multivariate
regression has been extended to include a sparse regularization term for high-
dimensional data, such as genetic data (Peng et al., 2010; Kim and Xing, 2012;
Chen et al., 2017). These methods consider the correlations of response
variables based on a priori data. Methods for reducing the dimension of the
response have also been introduced by Cook et al. (2010), Cook and Zhang
(2015), and Sun et al. (2015). The multivariate cluster elastic net, in
contrast, is proposed by Price and Sherwood (2018) as a multivariate
regression that can take into account the correlation of the response
variables without a prior information. This method includes an $L_{1}$ penalty
(Tibshirani, 1996) for the variable selection and $L_{2}$ penalty for
identifying clusters based on $k$-means (Forgy, 1965) that are associated with
the response variables. This is different from the previously mentioned
methods in that the estimation of the regression coefficients and grouping of
the response variables is conducted simultaneously. For data with a small
number of samples and a large number of explanatory variables, such as genetic
data, the lasso regularization term helps to estimate a sparse coefficient
matrix. The multivariate cluster elastic net can be also applied to high
dimensional data with multiple response variables. Another characteristic of
the multivariate cluster elastic net is the grouping of the fitted values of
the response variables by the $k$-means clustering term. With this term, the
multivariate cluster elastic net can improve estimation accuracy by
considering correlations between the response variables.

Meanwhile, in a regression model, the error distributions are assumed to be
Gaussian. However, when the errors follow a heavy-tailed distribution or
contain outliers, this may affect the estimation results [e.g., Huber, 2011].
Various robust regression methods have been proposed to solve this problem
(Fan et al., 2017; Loh, 2017; Sun et al., 2019; Lozano et al., 2016; Wang et
al., 2012a; Avella-Medina and Ronchetti, 2018; Prasad et al., 2020; Belloni et
al., 2011; Bradic et al., 2011; Wang et al., 2012b; Wang, 2013; Fan et al.,
2014; Sun and Zhang, 2012). Huber regression reduces the impact of outliers by
placing an adjustment parameter (Huber, 1964). Instead of least squares, the
lad-lasso ( (Wang et al., 2007) sets the least absolute deviation as its loss
function, and Zou and Yuan (2008) introduced the composite quantile
regression. The rank-based estimate with Wilcoxon scores (Jureckova, 1971;
Jaeckel, 1972; Hettmansperger and McKean, 1977) is also known to be more
efficient than least squared procedures when the data deviates from a normal
distribution. This is because it is robust against a response containing
outliers (Hettmansperger and McKean, 1998). This Wilcoxon estimate does not
need a tuning parameter for robust estimation. Wang et al. (2020) and Wang and
Li (2009) use a Wilcoxon estimate as their loss function with a regularization
penalty for robust estimation in a high-dimensional setting. However, these
methods remain in the framework of multiple regression.

We propose a Wilcoxon-type regression that is extended to a multivariate
cluster elastic net. It is a high-dimensional regression with robustness to
heavy-tailed error distribution or outliers containing random errors. We call
this new method the Wilcoxon-type multivariate cluster elastic net (WMCEN).
The Wilcoxon-type loss function is expected to be robust to heavy-tailed error
distribution or outliers. The proposed method is also suitable for a high-
dimensional situation because the variable selection uses the $L_{1}$
regularization penalty. Moreover, the proposed method obtains more efficiency
because the correlation of the response variables is considered, following a
similar framework to that of the multivariate cluster elastic net. When
updating the parameters of the proposed method, we derive the majorizing
function of the updated formula based on the majorize-minimization (MM)
algorithm (Hunter and Lange, 2004). With the majorizing function, an updated
formula based on the squared error criterion can be calculated, which ensures
the updated formula can be solved easily.

In Section $2$ of this paper, we show the objective function of our proposed
method, and we explain the derivation of the majorizing function. After
describing the algorithm and the updated formula of the proposed method, we
show the efficiency of the proposed method through numerical simulation in
Section $3$. In Section $4$, we apply the proposed method to real genetic data
related to breast cancer. After discussing based on the results, Section $5$
concludes the article with a discussion of further future issues.

##  2 Wilcoxon-type Multivariate Cluster Elastic Net (WMCEN)

In this section, we explain the optimization problem of our proposed method
and introduce the majorizing function corresponding to the new method. Before
that, we first explain the objective function of Wilcoxon-type regression,
which is the basis of the extension of our proposed method. Let $y_{i}\
(i=1,2,\cdots,n)$ be response variables, where $n$ is the number of subjects.
Let ${\bm{x}_{i}}\in\mathbb{R}^{p}\ (i=1,2,\cdots,n)$ be covariates in $p$
dimension, and ${\bm{\beta}}\in\mathbb{R}^{p}$ be unknown regression
coefficients in $p$ dimension. The Wilcoxon-type regression method for a
single outcome is estimated by minimizing $\bm{\beta}$:

$\displaystyle\sum_{i<j}|e_{i}-e_{j}|$ (1)

where
$e_{i}=y_{i}-\bm{x}_{i}^{T}\bm{\beta},e_{j}=y_{j}-\bm{x}_{j}^{T}\bm{\beta}\;(i,j=1,2,\cdots,n)$.
$e_{i}$ and $e_{j}$ are expressed in the form of a linear regression model,
and Eq. (1) represents the absolute sum of the difference of the residuals
between $i$ and $j$. This function in Eq. (1) is equivalent to Jaeckel’s
Wilcoxon-type dispersion function (Jaeckel, 1972):

$\displaystyle\sqrt{12}\sum_{i=1}^{n}\left[\frac{\bm{R}(y_{i}-\bm{x}_{i}^{T}\bm{\beta})}{n+1}-\frac{1}{2}\right](y_{i}-\bm{x}_{i}^{T}\bm{\beta})$
(2)

where $\bm{R}(y_{i}-\bm{x}_{i}^{T}\bm{\beta})$ is the rank of
$(y_{i}-\bm{x}_{i}^{T}\bm{\beta})\;(i=1,2,\cdots,n)$ (Hettmansperger and
McKean, 1978). That is to say, Eq. (2) can be regarded as a weighted
regression.

###  2.1 Optimization problem of Wilcoxon-type multivariate cluster elastic
net

Before describing the optimization problem of the proposed method, we
introduce the multivariate regression model in this subsection. Let
$\bm{Y}=(\bm{y_{1}},\bm{y_{2}},\cdots,\bm{y_{q}})^{T}\in\mathbb{R}^{n\times
q}$ be matrices of responses,
${\bm{X}}=(\bm{x_{1}},\bm{x_{2}},\cdots,\bm{x_{p}})^{T}\in\mathbb{R}^{n\times
p}$ be regression covariates,
${\bm{B}}=(\bm{\beta_{1}},\bm{\beta_{2}},\cdots,\bm{\beta_{q}})^{T}\in\mathbb{R}^{p\times
q}$ be unknown regression coefficients, and
$\bm{E}=(\bm{\epsilon_{1}},\bm{\epsilon_{2}},\cdots,\bm{\epsilon_{n}})^{T}$ be
i.i.d. random errors. The linear regression model is then represented as
follows:

$\displaystyle\bm{Y}=\bm{XB}+\bm{E}.$

Here, we explain the optimization problem of the WMCEN. To achieve this, we
extend the Wilcoxon-type regression method to a multivariate regression method
as follows:

$\displaystyle
L^{*}(\\{{\bm{\beta}}_{s}\\}_{s=1}^{q})=\sum_{s=1}^{q}\sum_{i<j}|e_{is}-e_{js}|\rightarrow{\rm
Min}$ (3)

where
$e_{is}=y_{is}-\bm{x}_{i}^{T}\bm{\beta}_{s},e_{js}=y_{js}-\bm{x}_{j}^{T}\bm{\beta}_{s}\;(i,j=1,2,\cdots,n;s=1,2,\cdots,q)$.
Eq.(3) represents the sum of the objective function of the Wilcoxon-type
regression method for each response variable. We extend the Wilcoxon-type
regression to the multivariate cluster elastic net (MCEN), and formulate the
optimization problem of WMCEN based on Eq.(3) as follows:

$\displaystyle
L(\\{{\bm{\beta}}_{s}\\}_{s=1}^{q},{\bm{U}},\\{{\bm{v}}_{\ell}\\}_{\ell=1}^{k})=\sum_{s=1}^{q}\sum_{i<j}|e_{is}-e_{js}|+\lambda\sum^{q}_{s=1}\|{\bm{\beta}}_{s}\|_{1}+\frac{\gamma}{2}\sum_{s=1}^{q}\sum_{\ell=1}^{k}u_{s\ell}\|{\bm{X}}{\bm{\beta}}_{s}-{\bm{X}}{\bm{v}}_{\ell}\|_{F}^{2}\rightarrow{\rm
Min}$ (4)

where $\lambda\;(\lambda>0)$ and $\gamma\;(\gamma>0)$ are the tuning
parameters. $\|\cdot\|_{F}$ is the Frobenius norm, $\|\cdot\|_{1}$ is the
$L_{1}$ norm, and $\|\cdot\|_{2}$ is the $L_{2}$ norm. The second and third
term of Eq. (4) are the same regularization terms as the multivariate cluster
elastic net (Price and Sherwood, 2018). ${\bm{U}}=(u_{s\ell})\
(s=1,2,\cdots,q;\ell=1,2,\cdots,k)$ is the indicator matrix, which represents
the degree of belonging to cluster $\ell$ on the response variable $s$. When
$u_{s\ell}=1$, ${\bm{X}}{\bm{\beta}}_{s}$ belongs to the $\ell$th cluster,
otherwise $u_{s\ell}=0$. $v_{\ell}\;(\ell=1,2,\cdots,k$) is the partial
regression coefficient for cluster centroid. The second formula term is the
lasso penalty (Tibshirani, 1996), and the third term of Eq.(4) is the
$k$-means clustering function.

In this study, we derive the updated formula using the MM algorithm (Hunter
and Lange, 2004). To solve the problem of applying the MM algorithm to the
lasso penalty (Hunter and Li, 2005), we substitute the perturbed version of
the lasso penalty term for the traditional lasso penalty term (Yu et al.,
2015).

$\displaystyle
L^{{\dagger}}(\\{{\bm{\beta}}_{s}\\}_{s=1}^{q},{\bm{U}},\\{{\bm{v}}_{\ell}\\}_{\ell=1}^{k})=$
$\displaystyle\sum_{s=1}^{q}\sum_{i<j}|e_{is}-e_{js}|+\lambda\sum_{s=1}^{q}\sum_{h=1}^{p}\left(|\beta_{sh}|-\epsilon\log\left(1+\frac{|\beta_{sh}|}{\epsilon}\right)\right)$
$\displaystyle+\frac{\gamma}{2}\sum_{s=1}^{q}\sum_{\ell=1}^{k}u_{s\ell}\|{\bm{X}}{\bm{\beta}}_{s}-{\bm{X}}{\bm{v}}_{\ell}\|_{F}^{2}$
(5)

where $\epsilon\;(\epsilon>0)$ is the hyper parameter of the lasso penalty.
When $\epsilon$ is approaching $0$, the second term of Eq. (5) converges to
the second term of Eq. (4).

###  2.2 Deriving the majorization function of WMCEN

In this subsection, we derive the majorizing function from Eq. (3), as shown
in Lemma 1, and the majorizing function of the perturbed lasso penalty is
derived in Lemma 2. With Lemma 1 and Lemma 2, we provide the majorizing
function of the proposed method in Proposition 1.

We first derive the majorizing function of the first term of Eq. (5). The
objective function of WMCEN in Eq. (3), that is the first term of Eq. (5), can
be transformed as follows:

$\displaystyle\sum_{s=1}^{q}\sum_{i<j}|e_{is}-e_{js}|$
$\displaystyle=\sum_{s=1}^{q}\sum_{i<j}|(y_{is}-\bm{x}_{i}^{T}\bm{\beta}_{s})-(y_{js}-\bm{x}_{j}^{T}\bm{\beta}_{s})|$
$\displaystyle=\sum_{s=1}^{q}\sum_{i<j}|(y_{is}-y_{js})-(\bm{x}_{i}-\bm{x}_{j})^{T}\bm{\beta}_{s}|$
$\displaystyle=\sum_{s=1}^{q}\sum_{i<j}|g_{ijs}-\bm{r}_{ij}^{T}\bm{\beta}_{s}|$
(6)

where $g_{ijs}=y_{is}-y_{js}$, and
$\bm{r}_{ij}=\bm{x}_{i}-\bm{x}_{j}\;(i,j=1,2,\cdots,n;\ i<j)$. $g_{ijs}$
represents the difference between the $i$ and $j$ of the $s$th response
variable, and $\bm{r}_{ij}$ is the difference between the $i$ and $j$ of the
covariate vectors. By forming them in Eq. (6), the difference of the residuals
can be expressed as a regression.

Here, we explain the majorizing function. $\theta$ is set as the parameter of
the interested real-valued objective function $f(\theta)$, and $\theta^{(t)}$
denotes the fixed value of parameter $\theta$. In the algorithm, for example,
$\theta^{(t)}$ represents the estimated $\theta$ at the $t$th step.
$g(\theta|\theta^{(t)})$ also denotes the function, and the
$g(\theta|\theta^{(t)})$ can easily derive the updated formula when given
$\theta^{(t)}$. The function $g(\theta|\theta^{(t)})$ is defined as the
majorizing function of $f(\theta)$ at $\theta^{(t)}$ when the following two
conditions are met:

$\displaystyle g(\theta|\theta^{(t)})$ $\displaystyle\geq f(\theta)\ {\rm for\
all}\ \theta,$ $\displaystyle g(\theta^{(t)}|\theta^{(t)})$
$\displaystyle=f(\theta^{(t)}).$

Next, we obtain the first term of the majorizing function from Eq. (6).

######  Lemma 1.

Given $g_{ijs},{\bm{r}}_{ij}$, and ${\bm{\beta}}_{s}^{(t)}$, for any
${\bm{\beta}}_{s}$, the following inequality holds:

$\displaystyle\sum_{s=1}^{q}\sum_{i<j}|g_{ijs}-{\bm{r}}_{ij}^{T}{\bm{\beta}}_{s}|\leq\sum_{s=1}^{q}\sum_{i<j}\frac{1}{2}\frac{|g_{ijs}-{\bm{r}}_{ij}^{T}{\bm{\beta}}_{s}|^{2}}{|g_{ijs}-{\bm{r}}_{ij}^{T}{\bm{\beta}}_{s}^{(t)}|}+\sum_{s=1}^{q}\sum_{i<j}\frac{1}{2}|g_{ijs}-{\bm{r}}_{ij}^{T}{\bm{\beta}}_{s}^{(t)}|.$
(7)

where ${\bm{\beta}}_{s}^{(t)}$ is the $t$th fixed ${\bm{\beta}}_{s}$ in the
algorithm. If ${\bm{\beta}}_{s}^{(t)}={\bm{\beta}}_{s}$ for all
$s=1,2,\cdots,q$, the equality of Eq. (7) holds.

###### Proof.

First, the inequality formula used to derive the majorizing function of the
objective function is explained. Eq. (8) is called the arithmetic-geometric
mean inequality (Hunter and Lange, 2004).

For real numbers $z_{1},z_{2},\cdots,z_{N}\geq 0$,

$\displaystyle\sqrt[N]{\prod_{l=1}^{N}z_{l}}\leq\frac{1}{N}\sum_{l=1}^{N}z_{l},$
(8)

only when $z_{1}=z_{2}=\cdots=z_{N}$, the equality of Eq. (8) holds.

Here, when $N=2$, $z_{1}=|g_{ijs}-{\bm{r}}_{ij}^{T}{\bm{\beta}}_{s}|$, and
$z_{2}=|g_{ijs}-{\bm{r}}_{ij}^{T}{\bm{\beta}}_{s}^{(t)}|$, it can be described
for any $i,j$, and $s$ by using Eq. (8) as follows:

$\displaystyle\sqrt{|g_{ijs}-{\bm{r}}_{ij}^{T}{\bm{\beta}}_{s}||g_{ijs}-{\bm{r}}_{ij}^{T}{\bm{\beta}}_{s}^{(t)}|}$
$\displaystyle\leq\frac{1}{2}(|g_{ijs}-{\bm{r}}_{ij}^{T}{\bm{\beta}}_{s}|+|g_{ijs}-{\bm{r}}_{ij}^{T}{\bm{\beta}}_{s}^{(t)}|)$
$\displaystyle\Longleftrightarrow\ \
|g_{ijs}-{\bm{r}}_{ij}^{T}{\bm{\beta}}_{s}||g_{ijs}-{\bm{r}}_{ij}^{T}{\bm{\beta}}_{s}^{(t)}|$
$\displaystyle\leq\frac{1}{4}(|g_{ijs}-{\bm{r}}_{ij}^{T}{\bm{\beta}}_{s}|+|g_{ijs}-{\bm{r}}_{ij}^{T}{\bm{\beta}}_{s}^{(t)}|)^{2}$
$\displaystyle\Longleftrightarrow\ \ \
\quad\quad\quad\qquad\quad|g_{ijs}-{\bm{r}}_{ij}^{T}{\bm{\beta}}_{s}|$
$\displaystyle\leq\frac{1}{2}\frac{|g_{ijs}-{\bm{r}}_{ij}^{T}{\bm{\beta}}_{s}|^{2}}{|g_{ijs}-{\bm{r}}_{ij}^{T}{\bm{\beta}}_{s}^{(t)}|}+\frac{1}{2}|g_{ijs}-{\bm{r}}_{ij}^{T}{\bm{\beta}}_{s}^{(t)}|.$
(9)

Because Eq. (9) holds for any $i,j$, and $s$, we have the following
inequality:

$\displaystyle\sum_{s=1}^{q}\sum_{i<j}|g_{ijs}-{\bm{r}}_{ij}^{T}{\bm{\beta}}_{s}|\leq\sum_{s=1}^{q}\sum_{i<j}\frac{1}{2}\frac{|g_{ijs}-{\bm{r}}_{ij}^{T}{\bm{\beta}}_{s}|^{2}}{|g_{ijs}-{\bm{r}}_{ij}^{T}{\bm{\beta}}_{s}^{(t)}|}+\sum_{s=1}^{q}\sum_{i<j}\frac{1}{2}|g_{ijs}-{\bm{r}}_{ij}^{T}{\bm{\beta}}_{s}^{(t)}|.$

∎

From Lemma 1, the first term of right side in Eq. (5) can be derived as a
quadratic function. Then, for any ${\bm{\beta}}_{s}$, this majorizing function
can be described:

$\displaystyle\sum_{s=1}^{q}\sum_{o=1}^{n(n-1)/2}w_{os}{|g_{os}-{\bm{r}}_{o}^{T}{\bm{\beta}}_{s}|^{2}}+w_{os}$
(10)

where $w_{os}=\frac{1}{2|g_{os}-{\bm{r}}_{o}^{T}{\bm{\beta}}_{s}^{(t)}|}$. Let
${\bm{\beta}}_{s}^{(t)}$ be a fixed $p$-dimensional vector, and the equality
of Eq. (10) is established when ${\bm{\beta}}_{s}={\bm{\beta}}_{s}^{(t)}$.
Because Eq. (10) is a weighted multivariate regression, it is possible to
obtain predicted values utilizing the conventional squared error criterion.

The second term of the majorizing function is then provided from the second
term of Eq. (5).

######  Lemma 2.

Given $\lambda\ (\lambda>0)$, $\epsilon\ (\epsilon>0)$, and fixed
$\beta_{sh}^{(t)}$, for any ${\bm{\beta}}_{s}$, the following inequality
holds:

$\displaystyle\lambda\sum_{s=1}^{q}\sum_{h=1}^{p}\left(|\beta_{sh}|-\epsilon\log\left(1+\frac{|\beta_{sh}|}{\epsilon}\right)\right)\leq\lambda\sum_{s=1}^{q}\sum_{h=1}^{p}\left(|{\beta}_{sh}^{(t)}|-\epsilon\log\left(1+\frac{|{\beta}_{sh}^{(t)}|}{\epsilon}\right)+\frac{({\beta}_{sh})^{2}-({\beta}_{sh}^{(t)})^{2}}{2(|{\beta}_{sh}^{(t)}|+\epsilon)}\right)$
(11)

where the equality of Eq. (11) holds only when $\beta_{sh}=\beta_{sh}^{(t)}$
for all $\beta_{sh}$.

###### Proof.

Eq. (11) can be proved in the same manner as Yu et al. (2015), and it can be
extended to the case of multivariate regression. For any $\beta_{sh}$, the two
functions are defined as below:

$\displaystyle Q_{\epsilon}(\beta_{sh})$
$\displaystyle=\left(|\beta_{sh}|-\epsilon\log\left(1+\frac{|\beta_{sh}|}{\epsilon}\right)\right),$
$\displaystyle S_{\epsilon}(\beta_{sh}|\beta_{sh}^{(t)})$
$\displaystyle=|{\beta}_{sh}^{(t)}|-\epsilon\log\left(1+\frac{|{\beta}_{sh}^{(t)}|}{\epsilon}\right)+\frac{({\beta}_{sh})^{2}-({\beta}_{sh}^{(t)})^{2}}{2(|{\beta}_{sh}^{(t)}|+\epsilon)}.$

Here, we define the following function:

$\displaystyle
H_{\epsilon}(\beta_{sh}|\beta_{sh}^{(t)})=S_{\epsilon}(\beta_{sh}|\beta_{sh}^{(t)})-Q_{\epsilon}(\beta_{sh}).$
(12)

We prove $H_{\epsilon}(\beta_{sh}|\beta_{sh}^{(t)})\geq 0$. First, when
$\beta_{sh}=\beta_{sh}^{(t)}$,

$\displaystyle H_{\epsilon}(\beta_{sh}^{(t)}|\beta_{sh}^{(t)})$
$\displaystyle=S_{\epsilon}(\beta_{sh}^{(t)}|\beta_{sh}^{(t)})-Q_{\epsilon}(\beta_{sh}^{(t)})$
$\displaystyle=\left(|{\beta}_{sh}^{(t)}|-\epsilon\log\left(1+\frac{|{\beta}_{sh}^{(t)}|}{\epsilon}\right)+\frac{({\beta}_{sh}^{(t)})^{2}-({\beta}_{sh}^{(t)})^{2}}{2(|{\beta}_{sh}^{(t)}|+\epsilon)}\right)-\left(|{\beta}_{sh}^{(t)}|-\epsilon\log\left(1+\frac{|{\beta}_{sh}^{(t)}|}{\epsilon}\right)\right)$
$\displaystyle=|{\beta}_{sh}^{(t)}|-\epsilon\log\left(1+\frac{|{\beta}_{sh}^{(t)}|}{\epsilon}\right)+\frac{({\beta}_{sh}^{(t)})^{2}-({\beta}_{sh}^{(t)})^{2}}{2(|{\beta}_{sh}^{(t)}|+\epsilon)}-|{\beta}_{sh}^{(t)}|+\epsilon\log\left(1+\frac{|{\beta}_{sh}^{(t)}|}{\epsilon}\right)$
$\displaystyle=0.$ (13)

Eq. (13) holds for any $\beta_{sh}=\beta_{sh}^{(t)}$; therefore, the equality
of Eq. (11) is satisfied.

Next, when $\beta_{sh}\neq\beta_{sh}^{(t)}$, we prove
$H_{\epsilon}(\beta_{sh}|\beta_{sh}^{(t)})>0$ in the cases below:

Case 1

If $0\leq\beta_{sh}<\beta_{sh}^{(t)}$,
$H_{\epsilon}(\beta_{sh}|\beta_{sh}^{(t)})$ monotonically decreases.

Eq. (12) is differentiated with respect to $\beta_{sh}$,

$\displaystyle\frac{\partial
H_{\epsilon}(\beta_{sh}|\beta_{sh}^{(t)})}{\partial{\beta}_{sh}}$
$\displaystyle=2\times\frac{1}{2(|{\beta}_{sh}^{(t)}|+\epsilon)}\times{\beta}_{sh}-1+\frac{\epsilon}{\epsilon+|{\beta}_{sh}|}$
$\displaystyle=\frac{{\beta}_{sh}}{|{\beta}_{sh}^{(t)}|+\epsilon}+\frac{\epsilon}{|{\beta}_{sh}|+\epsilon}-1.$
(14)

For Eq. (14), because ${\beta}_{sh}<{\beta}_{sh}^{(t)}$, the following
inequality holds:

$\displaystyle\frac{{\beta}_{sh}}{|{\beta}_{sh}^{(t)}|+\epsilon}+\frac{\epsilon}{|{\beta}_{sh}|+\epsilon}-1<\frac{{\beta}_{sh}}{|{\beta}_{sh}|+\epsilon}+\frac{\epsilon}{|{\beta}_{sh}|+\epsilon}-1=\frac{{\beta}_{sh}+\epsilon}{{\beta}_{sh}+\epsilon}-1=0.$

Therefore,

$\displaystyle\frac{\partial
H_{\epsilon}(\beta_{sh}|\beta_{sh}^{(t)})}{\partial{\beta}_{sh}}=\frac{{\beta}_{sh}}{|{\beta}_{sh}^{(t)}|+\epsilon}+\frac{\epsilon}{|{\beta}_{sh}|+\epsilon}-1<0.$
(15)

Eq. (15) shows that $H_{\epsilon}(\beta_{sh}|\beta_{sh}^{(t)})$ monotonically
decreases for $0\leq\beta_{sh}<\beta_{sh}^{(t)}$.

Case 2

If $0\leq\beta_{sh}^{(t)}<\beta_{sh}$,
$H_{\epsilon}(\beta_{sh}|\beta_{sh}^{(t)})$ monotonically increases.

Because $\beta_{sh}^{(t)}<\beta_{sh}$, the following inequality is satisfied:

$\displaystyle\frac{\partial
H_{\epsilon}(\beta_{sh}|\beta_{sh}^{(t)})}{\partial{\beta}_{sh}}$
$\displaystyle=\frac{{\beta}_{sh}}{|{\beta}_{sh}^{(t)}|+\epsilon}+\frac{\epsilon}{|{\beta}_{sh}|+\epsilon}-1>\frac{{\beta}_{sh}}{|{\beta}_{sh}|+\epsilon}+\frac{\epsilon}{|{\beta}_{sh}|+\epsilon}-1=\frac{{\beta}_{sh}+\epsilon}{{\beta}_{sh}+\epsilon}-1=0.$

Therefore,

$\displaystyle\frac{\partial
H_{\epsilon}(\beta_{sh}|\beta_{sh}^{(t)})}{\partial{\beta}_{sh}}=\frac{{\beta}_{sh}}{|{\beta}_{sh}^{(t)}|+\epsilon}+\frac{\epsilon}{|{\beta}_{sh}|+\epsilon}-1>0.$
(16)

Eq. (16) shows that $H_{\epsilon}(\beta_{sh}|\beta_{sh}^{(t)})$ monotonically
increases when $0\leq\beta_{sh}^{(t)}<\beta_{sh}$.

Case 3

If $\beta_{sh}^{(t)}<\beta_{sh}\leq 0$,
$H_{\epsilon}(\beta_{sh}|\beta_{sh}^{(t)})$ monotonically increases.

Eq. (12) is differentiated with respect to $\beta_{sh}$.

$\displaystyle\frac{\partial
H_{\epsilon}(\beta_{sh}|\beta_{sh}^{(t)})}{\partial{\beta}_{sh}}$
$\displaystyle=2\times\frac{1}{2(|\beta_{sh}^{(t)}|+\epsilon)}\times\beta_{sh}-(-1)-\frac{\epsilon}{\epsilon+|\beta_{sh}|}$
$\displaystyle=\frac{\beta_{sh}}{|\beta_{sh}^{(t)}|+\epsilon}-\frac{\epsilon}{|\beta_{sh}|+\epsilon}+1.$

For $\beta_{sh}^{(t)}<\beta_{sh}$, the following inequality holds:

$\displaystyle\frac{\partial
H_{\epsilon}(\beta_{sh}|\beta_{sh}^{(t)})}{\partial{\beta}_{sh}}=\frac{{\beta}_{sh}}{|{\beta}_{sh}^{(t)}|+\epsilon}-\frac{\epsilon}{|{\beta}_{sh}|+\epsilon}+1>\frac{{\beta}_{sh}}{|{\beta}_{sh}|+\epsilon}-\frac{\epsilon}{|{\beta}_{sh}|+\epsilon}+1.$

Because $\beta_{sh}\leq 0$,

$\displaystyle\frac{-|{\beta}_{sh}|-\epsilon}{|{\beta}_{sh}|+\epsilon}+1=\frac{-(|{\beta}_{sh}|+\epsilon)}{|{\beta}_{sh}|+\epsilon}+1=0.$

Therefore,

$\displaystyle\frac{\partial
H_{\epsilon}(\beta_{sh}|\beta_{sh}^{(t)})}{\partial{\beta}_{sh}}=\frac{{\beta}_{sh}}{|{\beta}_{sh}^{(t)}|+\epsilon}-\frac{\epsilon}{|{\beta}_{sh}|+\epsilon}+1>0.$
(17)

Eq. (17) indicates that $H_{\epsilon}({\beta_{sh}|\beta_{sh}^{(t)}})$
monotonically increases for $0>\beta_{sh}>-\beta_{sh}^{(t)}$.

Case 4

If $\beta_{sh}<\beta_{sh}^{(t)}\leq 0$,
$H_{\epsilon}(\beta_{sh}|\beta_{sh}^{(t)})$ monotonically decreases.

Because $\beta_{sh}<-\beta_{sh}^{(t)}$, the following inequality holds:

$\displaystyle\frac{\partial
H_{\epsilon}(\beta_{sh}|\beta_{sh}^{(t)})}{\partial{\beta}_{sh}}=\frac{{\beta}_{sh}}{|{\beta}_{sh}^{(t)}|+\epsilon}-\frac{\epsilon}{|{\beta}_{sh}|+\epsilon}+1<\frac{{\beta}_{sh}}{|{\beta}_{sh}|+\epsilon}-\frac{\epsilon}{|{\beta}_{sh}|+\epsilon}+1.$

Because $\beta_{sh}\leq 0$,

$\displaystyle\frac{-|{\beta}_{sh}|-\epsilon}{|{\beta}_{sh}|+\epsilon}+1=\frac{-(|{\beta}_{sh}|+\epsilon)}{|{\beta}_{sh}|+\epsilon}+1=0.$

Therefore,

$\displaystyle\frac{\partial
H_{\epsilon}(\beta_{sh}|\beta_{sh}^{(t)})}{\partial{\beta}_{sh}}=\frac{{\beta}_{sh}}{|{\beta}_{sh}^{(t)}|+\epsilon}-\frac{\epsilon}{|{\beta}_{sh}|+\epsilon}+1<0.$
(18)

Eq. (18) shows that $H_{\epsilon}({\beta_{sh}|\beta_{sh}^{(t)}})$
monotonically increases for $\beta_{sh}<\beta_{sh}^{(t)}\leq 0$.

From Eq. (13) and ${\rm Cases1,2,3}$, and $4$, for any $\beta_{sh}$,
$H_{\epsilon}(\beta_{sh}|\beta_{sh}^{(t)})\geq 0$. Therefore, Eq. (11) is
satisfied.

∎

With ${\rm Lemma2}$, we show the majorizing function can be derived from the
perturbed version of the lasso penalty as the squared form. The right-hand
side of Eq. (11) can be expressed as

$\displaystyle\lambda\sum_{s=1}^{q}\sum_{h=1}^{p}\left(|{\beta}_{sh}^{(t)}|-\epsilon\log\left(1+\frac{|{\beta}_{sh}^{(t)}|}{\epsilon}\right)+\frac{({\beta}_{sh})^{2}-({\beta}_{sh}^{(t)})^{2}}{2(|{\beta}_{sh}^{(t)}|+\epsilon)}\right)=\|{\bm{\Psi}}_{s}{\bm{\beta}}_{s}\|_{F}^{2}+C$
(19)

where

$\displaystyle{\bm{\Psi}}_{s}={\rm
diag}\left(\frac{1}{\sqrt{2({\beta}_{s1}^{(t)}+\epsilon)}},\frac{1}{\sqrt{2({\beta}_{s2}^{(t)}+\epsilon)}},\cdots,\frac{1}{\sqrt{2({\beta}_{sp}^{(t)}+\epsilon)}}\right),$

and $C$ represents constant values that are not relevant to $\beta_{sh}$. From
Lemma 1 and Lemma 2, the following proposition holds.

######  Proposition 1.

Given ${\bm{Y}},{\bm{X}},\lambda\;(\lambda>0),\gamma\;(\gamma>0),$ and
$\epsilon\;(\epsilon>0)$, the following inequality is satisfied:

$\displaystyle
L^{{\dagger}}(\\{{\bm{\beta}}_{s}\\}_{s=1}^{q},{\bm{U}},\\{{\bm{v}}_{\ell}\\}_{\ell=1}^{k})$
$\displaystyle=$
$\displaystyle\sum_{s=1}^{q}\sum_{i<j}|e_{is}-e_{js}|+\lambda\sum_{s=1}^{q}\sum_{h=1}^{p}\left(|\beta_{sh}|-\epsilon\log\left(1+\frac{|\beta_{sh}|}{\epsilon}\right)\right)+\frac{\gamma}{2}\sum_{s=1}^{q}\sum_{\ell=1}^{k}u_{s\ell}\|{\bm{X}}{\bm{\beta}}_{s}-{\bm{X}}{\bm{v}}_{\ell}\|_{F}^{2}$
$\displaystyle\leq$
$\displaystyle\sum_{s=1}^{q}\sum_{o=1}^{n(n-1)/2}w_{os}{|g_{os}-{\bm{r}}_{o}^{T}{\bm{\beta}}_{s}|^{2}}+\lambda\sum_{s=1}^{q}\|{\bm{\Psi}}_{s}{\bm{\beta}}_{s}\|_{F}^{2}+C+\frac{\gamma}{2}\sum_{\ell=1}^{k}\sum_{s=1}^{q}u_{s\ell}\|{\bm{X}}{\bm{\beta}}_{s}-{\bm{X}}{\bm{v}}_{\ell}\|_{F}^{2}$
$\displaystyle=$ $\displaystyle
M(\\{{\bm{\beta}}_{s}\\}_{s=1}^{q},{\bm{U}},\\{{\bm{v}}_{\ell}\\}_{\ell=1}^{k}|\\{{\bm{\beta}}_{s}^{(t)}\\}_{s=1}^{q})$
(20)

where the equality holds if ${\bm{\beta}}_{s}={\bm{\beta}}_{s}^{(t)}$ in all
$(s,h)\;(s=1,2,\cdots,q;\ h=1,2,\cdots,p)$.

From Eq. (20), the updated formula of the proposed method can be calculated
based on the least squared criteria. Therefore, it can be simply expressed.

###  2.3 Algorithm and updated formula

This subsection explains the algorithm and updated formula of the proposed
method. The proposed method updates ${\bm{\beta}}_{s}$, ${\bm{U}}$, and
${\bm{v}}_{\ell}$ based on the alternate least squares criterion (Young, 1981)
by using the majorizing function of the proposed method.

In Proposition 2, we show the updated formula of ${\bm{\beta}}_{s}$ based on
the majorizing function. We then show the procedure for the estimation of
${\bm{U}}$ and ${\bm{v}}_{\ell}$ based on $k$-means (Forgy, 1965). After that,
we present the algorithm of the proposed method.

######  Proposition 2.

Given
${\bm{X}},w_{os},g_{os},{\bm{r}}_{o}\;(o=1,2,\cdots,n(n-1)/2;s=1,2,\cdots,q),\lambda\;(\lambda>0),\gamma\;(\gamma>0),{\bm{\Psi}}_{s},$
and $u_{s\ell}\;(s=1,2,\cdots,q;\ell=1,2,\cdots,k),$ the updated formula of
${\bm{\beta}}_{d}\ (d=1,2,\cdots,q)$ is derived as follows:

$\displaystyle{\bm{\beta}}_{d}^{(t+1)}\leftarrow$
$\displaystyle\left(2\left({\sum_{o=1}^{n(n-1)/2}w_{os}^{(t+1)}{\bm{r}}_{o}{\bm{r}}_{o}^{T}}\right)+2\left(\lambda{\bm{\Psi}}_{d}^{(t+1)T}{\bm{\Psi}}_{d}^{(t+1)}\right)+\gamma\sum_{\ell=1}^{k}u_{d\ell}^{(t^{\dagger})}{\bm{X}}^{T}{\bm{X}}-\gamma\sum_{\ell=1}^{k}\frac{1}{q_{\ell}}u_{d\ell}^{(t^{\dagger})}{\bm{X}}^{T}{\bm{X}}\right)^{-1}$
$\displaystyle\left(2\left(\sum_{o=1}^{n(n-1)/2}w_{os}^{(t+1)}g_{os}{\bm{r}}_{o}\right)+\gamma\sum_{\ell=1}^{k}\frac{1}{q_{\ell}}u_{d\ell}^{(t^{\dagger})}\sum_{m\neq
d}^{q}u_{m\ell}^{(t^{\dagger})}{\bm{X}}^{T}{\bm{X}}{\bm{\beta}}_{m}^{(f(m,s))}\right)$
(21)

where $q_{\ell}$ indicates the number of ${\bm{X}}{\bm{\beta}}_{s}$ belonging
to cluster $\ell$. ${\bm{\beta}}_{s}$ is assumed to be updated in the order
$s=1,2,\cdots,d,\cdots,q$ as Algorithm 1. ${\bm{\beta}}_{s}^{(t)}$ and
${\bm{\beta}}_{s}^{(t+1)}$ are coefficients at the $t$th and $(t+1)$th steps,
respectively. The function $f(m,s)$ is defined as $t+1$ when $m<s$, whereas it
is defined as $t$ when $m>s$.

###### Proof.

We provide the results of the derivatives of each term on the right-hand side
of Eq. (20) and an updated formula based on these results.

For the first term on the right side of Eq. (20),

$\displaystyle M_{1}(\\{{\bm{\beta}}_{s}\\}_{s=1}^{q})$
$\displaystyle=\sum_{s=1}^{q}\sum_{o=1}^{n(n-1)/2}w_{os}{|g_{os}-{\bm{r}}_{o}^{T}{\bm{\beta}}_{s}|^{2}}$
$\displaystyle=\sum_{s=1}^{q}\sum_{o=1}^{n(n-1)/2}w_{os}(g_{os}^{2}-2g_{os}{\bm{r}}_{o}^{T}{\bm{\beta}}_{s}+{\bm{\beta}}_{s}^{T}{\bm{r}}_{o}{\bm{r}}_{o}^{T}{\bm{\beta}}_{s})$
$\displaystyle=\sum_{s=1}^{q}\sum_{o=1}^{n(n-1)/2}w_{os}g_{os}^{2}-2\sum_{s=1}^{q}{\sum_{o=1}^{n(n-1)/2}w_{os}g_{os}{\bm{r}}_{o}^{T}}{\bm{\beta}}_{s}+\sum_{s=1}^{q}{\bm{\beta}}_{s}^{T}\left({\sum_{o=1}^{n(n-1)/2}w_{os}{\bm{r}}_{o}{\bm{r}}_{o}^{T}}\right){\bm{\beta}}_{s}$
(22)

where the second and third terms in Eq. (22) are related to
${\bm{\beta}}_{s}$.

Here, Eq. (22) is differentiated by ${\bm{\beta}}_{d}\;(d=1,2,\cdots,q)$, and
we have

$\displaystyle\frac{\partial
M_{1}(\\{{\bm{\beta}}_{s}\\}_{s=1}^{q})}{\partial{\bm{\beta}}_{d}}=-2\left(\sum_{o=1}^{n(n-1)/2}w_{od}g_{od}{\bm{r}}_{o}\right)+2\left({\sum_{o=1}^{n(n-1)/2}w_{od}{\bm{r}}_{o}{\bm{r}}_{o}^{T}}\right){\bm{\beta}}_{d}.$
(23)

Next, the second term of Eq. (20) can be described as follows:

$\displaystyle
M_{2}(\\{{\bm{\beta}}_{s}\\}_{s=1}^{q})=\lambda\sum_{s=1}^{q}\|{\bm{\Psi}}_{s}{\bm{\beta}}_{s}\|_{F}^{2}=\lambda\sum_{s=1}^{q}({\bm{\Psi}}_{s}{\bm{\beta}}_{s})^{T}({\bm{\Psi}}_{s}{\bm{\beta}}_{s})=\lambda\sum_{s=1}^{q}({\bm{\beta}}_{s}^{T}{\bm{\Psi}}_{s}^{T}{\bm{\Psi}}_{s}{\bm{\beta}}_{s})=\sum_{s=1}^{q}{\bm{\beta}}_{s}^{T}(\lambda{\bm{\Psi}}_{s}^{T}{\bm{\Psi}}_{s}){\bm{\beta}}_{s}.$
(24)

Here, Eq. (24) is differentiated by ${\bm{\beta}}_{d}$, and we have

$\displaystyle\frac{\partial
M_{2}(\\{{\bm{\beta}}_{s}\\}_{s=1}^{q})}{\partial{\bm{\beta}}_{d}}=2\lambda{\bm{\Psi}}_{d}^{T}{\bm{\Psi}}_{d}{\bm{\beta}}_{d}.$
(25)

Then, the third term of the right side of Eq. (20) can be described as
follows:

$\displaystyle M_{3}(\\{{\bm{\beta}}_{s}\\}_{s=1}^{q})=$
$\displaystyle\frac{\gamma}{2}\sum_{\ell=1}^{k}\sum_{s=1}^{q}u_{s\ell}\|{\bm{X}}{\bm{\beta}}_{s}-{\bm{X}}{\bm{v}}_{\ell}\|^{2}$
$\displaystyle=$
$\displaystyle\frac{\gamma}{2}\sum_{\ell=1}^{k}\sum_{s=1}^{q}u_{s\ell}({\bm{X}}{\bm{\beta}}_{s}-{\bm{X}}{\bm{v}}_{\ell})^{T}({\bm{X}}{\bm{\beta}}_{s}-{\bm{X}}{\bm{v}}_{\ell})$
$\displaystyle=$
$\displaystyle\frac{\gamma}{2}\sum_{\ell=1}^{k}\sum_{s=1}^{q}u_{s\ell}({\bm{\beta}}_{s}^{T}{\bm{X}}^{T}{\bm{X}}{\bm{\beta}}_{s}-2{\bm{\beta}}_{s}^{T}{\bm{X}}^{T}{\bm{X}}{\bm{v}}_{\ell}+{\bm{v}}_{\ell}^{T}{\bm{X}}^{T}{\bm{X}}{\bm{v}}_{\ell})$
$\displaystyle=$
$\displaystyle\frac{\gamma}{2}\sum_{\ell=1}^{k}\sum_{s=1}^{q}u_{s\ell}{\bm{\beta}}_{s}^{T}{\bm{X}}^{T}{\bm{X}}{\bm{\beta}}_{s}-\gamma\sum_{\ell=1}^{k}\sum_{s=1}^{q}u_{s\ell}{\bm{\beta}}_{s}^{T}{\bm{X}}^{T}{\bm{X}}{\bm{v}}_{\ell}+\frac{\gamma}{2}\sum_{\ell=1}^{k}\sum_{s=1}^{q}u_{s\ell}{\bm{v}}_{\ell}^{T}{\bm{X}}^{T}{\bm{X}}{\bm{v}}_{\ell}.$
(26)

Before deriving the updated formula, to find $v_{\ell}$ for the optimization
of the right side of Eq. (20), we have :

$\displaystyle\frac{\partial M_{3}(\\{{\bm{\beta}}_{s}\\}_{s=1}^{q})}{\partial
v_{\ell}}=-\gamma\sum_{s=1}^{q}u_{s\ell}{\bm{X}}^{T}{\bm{X}}{\bm{\beta}}_{s}+\gamma\sum_{s=1}^{q}u_{s\ell}{\bm{X}}^{T}{\bm{X}}{\bm{v}}_{\ell}={\bm{0}}_{p}$
$\displaystyle\Longleftrightarrow-\gamma\sum_{s=1}^{q}u_{s\ell}{\bm{X}}^{T}{\bm{X}}{\bm{\beta}}_{s}+\gamma\sum_{s=1}^{q}u_{s\ell}{\bm{X}}^{T}{\bm{X}}{\bm{v}}_{\ell}={\bm{0}}_{p}$
$\displaystyle\Longleftrightarrow\sum_{s=1}^{q}u_{s\ell}{\bm{X}}^{T}{\bm{X}}{\bm{v}}_{\ell}=\sum_{s=1}^{q}u_{s\ell}{\bm{X}}^{T}{\bm{X}}{\bm{\beta}}_{s}$
$\displaystyle\Longleftrightarrow{\bm{v}}_{\ell}=\frac{1}{q_{\ell}}\sum_{s=1}^{q}u_{s\ell}{\bm{\beta}}_{s}.$
(27)

where ${\bm{0}}_{p}$ is the element is $0$ with the length of $p$.

Here, when substituting Eq. (27) into Eq. (26), we have

$\displaystyle
M_{3}^{\dagger}(\\{{\bm{\beta}}_{s}\\}_{s=1}^{q})=\frac{\gamma}{2}\sum_{\ell=1}^{k}\sum_{s=1}^{q}u_{s\ell}{\bm{\beta}}_{s}^{T}{\bm{X}}^{T}{\bm{X}}{\bm{\beta}}_{s}$
$\displaystyle-\gamma\sum_{\ell=1}^{k}\sum_{s=1}^{q}u_{s\ell}{\bm{\beta}}_{s}^{T}{\bm{X}}^{T}{\bm{X}}\left(\frac{1}{q_{\ell}}\sum_{m=1}^{q}u_{m\ell}{\bm{\beta}}_{m}\right)$
$\displaystyle+\frac{\gamma}{2}\sum_{\ell=1}^{k}\sum_{s=1}^{q}u_{s\ell}\left(\frac{1}{q_{\ell}}\sum_{m=1}^{q}u_{m\ell}{\bm{\beta}}_{m}\right)^{T}{\bm{X}}^{T}{\bm{X}}\left(\frac{1}{q_{\ell}}\sum_{m=1}^{q}u_{m\ell}{\bm{\beta}}_{m}\right)$
(28)

To derive the updated formula of ${\bm{\beta}}_{d}$, the second term of Eq.
(28) can be divided into a term relevant to ${\bm{\beta}}_{d}$ and a term not
relevant to ${\bm{\beta}}_{d}$:

$\displaystyle M^{{\dagger}}_{3T2}(\\{{\bm{\beta}}_{s}\\}_{s=1}^{q})=$
$\displaystyle-\gamma\sum_{\ell=1}^{k}\sum_{s=1}^{q}u_{s\ell}{\bm{\beta}}_{s}^{T}{\bm{X}}^{T}{\bm{X}}\left(\frac{1}{q_{\ell}}\sum_{m=1}^{q}u_{m\ell}{\bm{\beta}}_{m}\right)$
$\displaystyle=$
$\displaystyle-\gamma\sum_{\ell=1}^{k}u_{d\ell}{\bm{\beta}}_{d}^{T}{\bm{X}}^{T}{\bm{X}}\left(\frac{1}{q_{\ell}}u_{d\ell}{\bm{\beta}}_{d}+\frac{1}{q_{\ell}}\sum_{m\neq
d}^{q}u_{m\ell}{\bm{\beta}}_{m}\right)$
$\displaystyle-\gamma\sum_{\ell=1}^{k}\sum_{s\neq
d}^{q}u_{s\ell}{\bm{\beta}}_{s}^{T}{\bm{X}}^{T}{\bm{X}}\left(\frac{1}{q_{\ell}}u_{d\ell}{\bm{\beta}}_{d}+\frac{1}{q_{\ell}}\sum_{m\neq
d}^{q}u_{m\ell}{\bm{\beta}}_{m}\right)$ $\displaystyle=$
$\displaystyle-\gamma\sum_{\ell=1}^{k}u_{d\ell}{\bm{\beta}}_{d}^{T}{\bm{X}}^{T}{\bm{X}}\frac{1}{q_{\ell}}u_{d\ell}{\bm{\beta}}_{d}-\gamma\sum_{\ell=1}^{k}u_{d\ell}{\bm{\beta}}_{d}^{T}{\bm{X}}^{T}{\bm{X}}\frac{1}{q_{\ell}}\sum_{m\neq
d}^{q}u_{m\ell}{\bm{\beta}}_{m}$
$\displaystyle-\gamma\sum_{\ell=1}^{k}\sum_{s\neq
d}^{q}u_{s\ell}{\bm{\beta}}_{s}^{T}{\bm{X}}^{T}{\bm{X}}\frac{1}{q_{\ell}}u_{d\ell}{\bm{\beta}}_{d}-\gamma\sum_{\ell=1}^{k}\sum_{s\neq
d}^{q}u_{s\ell}{\bm{\beta}}_{s}^{T}{\bm{X}}^{T}{\bm{X}}\frac{1}{q_{\ell}}\sum_{m\neq
d}^{q}u_{m\ell}{\bm{\beta}}_{m}$ $\displaystyle=$
$\displaystyle-{\bm{\beta}}_{d}^{T}\gamma\sum_{\ell=1}^{k}{\bm{X}}^{T}{\bm{X}}\frac{1}{q_{\ell}}u_{d\ell}{\bm{\beta}}_{d}-{\bm{\beta}}_{d}^{T}\gamma\sum_{\ell=1}^{k}u_{d\ell}{\bm{X}}^{T}{\bm{X}}\frac{1}{q_{\ell}}\sum_{m\neq
d}^{q}u_{m\ell}{\bm{\beta}}_{m}$
$\displaystyle-{\bm{\beta}}_{d}^{T}\gamma\sum_{\ell=1}^{k}u_{d\ell}{\bm{X}}^{T}{\bm{X}}\frac{1}{q_{\ell}}\sum_{m\neq
d}^{q}u_{m\ell}{\bm{\beta}}_{m}-\gamma\sum_{\ell=1}^{k}\sum_{s\neq
d}^{q}u_{s\ell}{\bm{\beta}}_{s}^{T}{\bm{X}}^{T}{\bm{X}}\frac{1}{q_{\ell}}\sum_{m\neq
d}^{q}u_{m\ell}{\bm{\beta}}_{m}$ $\displaystyle=$
$\displaystyle-{\bm{\beta}}_{d}^{T}\left(\gamma\sum_{\ell=1}^{k}{\bm{X}}^{T}{\bm{X}}\frac{1}{q_{\ell}}u_{d\ell}\right){\bm{\beta}}_{d}-{\bm{\beta}}_{d}^{T}\left(2\gamma\sum_{\ell=1}^{k}u_{d\ell}{\bm{X}}^{T}{\bm{X}}\frac{1}{q_{\ell}}\sum_{m\neq
d}^{q}u_{m\ell}{\bm{\beta}}_{m}\right)$
$\displaystyle-\gamma\sum_{\ell=1}^{k}\sum_{s\neq
d}^{q}u_{s\ell}{\bm{\beta}}_{s}^{T}{\bm{X}}^{T}{\bm{X}}\frac{1}{q_{\ell}}\sum_{m\neq
d}^{q}u_{m\ell}{\bm{\beta}}_{m}.$ (29)

Similarly, the third term of Eq. (28) can be divided into a term relevant to
${\bm{\beta}}_{d}$ and one not relevant to ${\bm{\beta}}_{d}$.

$\displaystyle M^{{\dagger}}_{3T3}(\\{{\bm{\beta}}_{s}\\}_{s=1}^{q})=$
$\displaystyle\frac{\gamma}{2}\sum_{\ell=1}^{k}q_{\ell}\left(\frac{1}{q_{\ell}}\sum_{m=1}^{q}u_{m\ell}{\bm{\beta}}_{m}\right)^{T}{\bm{X}}^{T}{\bm{X}}\left(\frac{1}{q_{\ell}}\sum_{m=1}^{q}u_{m\ell}{\bm{\beta}}_{m}\right)$
$\displaystyle=$
$\displaystyle\frac{\gamma}{2}\sum_{\ell=1}^{k}q_{\ell}\left(\frac{1}{q_{\ell}}u_{d\ell}{\bm{\beta}}_{d}+\frac{1}{q_{\ell}}\sum_{m\neq
d}^{q}u_{m\ell}{\bm{\beta}}_{m}\right)^{T}{\bm{X}}^{T}{\bm{X}}\left(\frac{1}{q_{\ell}}u_{d\ell}{\bm{\beta}}_{d}+\frac{1}{q_{\ell}}\sum_{m\neq
d}^{q}u_{m\ell}{\bm{\beta}}_{m}\right)$ $\displaystyle=$
$\displaystyle\frac{\gamma}{2}\sum_{\ell=1}^{k}q_{\ell}\left(\frac{1}{q_{\ell}}u_{d\ell}{\bm{\beta}}_{d}+\frac{1}{q_{\ell}}\sum_{m\neq
d}^{q}u_{m\ell}{\bm{\beta}}_{m}\right)^{T}{\bm{X}}^{T}{\bm{X}}\left(\frac{1}{q_{\ell}}u_{d\ell}{\bm{\beta}}_{d}\right)$
$\displaystyle+\frac{\gamma}{2}\sum_{\ell=1}^{k}q_{\ell}\left(\frac{1}{q_{\ell}}u_{d\ell}{\bm{\beta}}_{d}+\frac{1}{q_{\ell}}\sum_{m\neq
d}^{q}u_{m\ell}{\bm{\beta}}_{m}\right)^{T}{\bm{X}}^{T}{\bm{X}}\left(\frac{1}{q_{\ell}}\sum_{m\neq
d}^{q}u_{m\ell}{\bm{\beta}}_{m}\right)$ $\displaystyle=$
$\displaystyle\frac{\gamma}{2}\sum_{\ell=1}^{k}u_{d\ell}{\bm{\beta}}_{d}^{T}{\bm{X}}^{T}{\bm{X}}\frac{1}{q_{\ell}}u_{d\ell}{\bm{\beta}}_{d}+\frac{\gamma}{2}\sum_{\ell=1}^{k}\sum_{m\neq
d}^{q}u_{m\ell}{\bm{\beta}_{m}^{T}}{\bm{X}}^{T}{\bm{X}}\frac{1}{q_{\ell}}u_{d\ell}{\bm{\beta}}_{d}$
$\displaystyle+\frac{\gamma}{2}\sum_{\ell=1}^{k}u_{d\ell}{\bm{X}}^{T}{\bm{X}}\frac{1}{q_{\ell}}\sum_{m\neq
d}^{q}u_{m\ell}{\bm{\beta}}_{m}+\frac{\gamma}{2}\sum_{\ell=1}^{k}\left(\sum_{m\neq
d}^{q}u_{m\ell}{\bm{\beta}}_{m}^{T}\right){\bm{X}}^{T}{\bm{X}}\left(\frac{1}{q_{\ell}}\sum_{m\neq
d}^{q}u_{m\ell}{\bm{\beta}}_{m}\right)$ $\displaystyle=$
$\displaystyle{\bm{\beta}}_{d}^{T}\left(\frac{\gamma}{2}\sum_{\ell=1}^{k}u_{d\ell}{\bm{X}}^{T}{\bm{X}}\frac{1}{q_{\ell}}u_{d\ell}\right){\bm{\beta}}_{d}+{\bm{\beta}}_{d}^{T}\left(\gamma\sum_{\ell=1}^{k}u_{d\ell}{\bm{X}}^{T}{\bm{X}}\frac{1}{q_{\ell}}\sum_{m\neq
d}^{q}u_{m\ell}{\bm{\beta}}_{m}\right)$
$\displaystyle+\frac{\gamma}{2}\sum_{\ell=1}^{k}\sum_{m\neq
d}u_{m\ell}{\bm{\beta}}_{m}^{T}{\bm{X}}^{T}{\bm{X}}\frac{1}{q_{\ell}}\sum_{m\neq
d}^{q}u_{m\ell}{\bm{\beta}}_{m}$ (30)

where the first and second terms of Eq. (30) are relevant to
${\bm{\beta}}_{d}$. Using Eq. (29) and Eq. (30), Eq.(28) can be described as
follows:

$\displaystyle M_{3}^{*}(\\{{\bm{\beta}}_{s}\\}_{s=1}^{q})=$
$\displaystyle\frac{\gamma}{2}\sum_{\ell=1}^{k}\sum_{s=1}^{q}u_{s\ell}{\bm{\beta}}_{s}^{T}{\bm{X}}^{T}{\bm{X}}{\bm{\beta}}_{s}-{\bm{\beta}}_{d}^{T}\left(\gamma\sum_{\ell=1}^{k}{\bm{X}}^{T}{\bm{X}}\frac{1}{q_{\ell}}u_{d\ell}\right){\bm{\beta}}_{d}-{\bm{\beta}}_{d}^{T}\left(2\gamma\sum_{\ell=1}^{k}u_{d\ell}{\bm{X}}^{T}{\bm{X}}\frac{1}{q_{\ell}}\sum_{m\neq
d}^{q}u_{m\ell}{\bm{\beta}}_{m}\right)$
$\displaystyle-\gamma\sum_{\ell=1}^{k}\sum_{s\neq
d}^{q}u_{s\ell}{\bm{\beta}}_{s}^{T}{\bm{X}}^{T}{\bm{X}}\frac{1}{q_{\ell}}\sum_{m\neq
d}^{q}u_{m\ell}{\bm{\beta}}_{m}$
$\displaystyle+{\bm{\beta}}_{d}^{T}\left(\frac{\gamma}{2}\sum_{\ell=1}^{k}u_{d\ell}{\bm{X}}^{T}{\bm{X}}\frac{1}{q_{\ell}}u_{d\ell}\right){\bm{\beta}}_{d}+{\bm{\beta}}_{d}^{T}\left(\gamma\sum_{\ell=1}^{k}u_{d\ell}{\bm{X}}^{T}{\bm{X}}\frac{1}{q_{\ell}}\sum_{m\neq
d}^{q}u_{m\ell}{\bm{\beta}}_{m}\right)$
$\displaystyle+\frac{\gamma}{2}\sum_{\ell=1}^{k}\sum_{m\neq
d}u_{m\ell}{\bm{\beta}}_{m}^{T}{\bm{X}}^{T}{\bm{X}}\frac{1}{q_{\ell}}\sum_{m\neq
d}u_{m\ell}{\bm{\beta}}_{m}$ $\displaystyle=$
$\displaystyle\frac{\gamma}{2}\sum_{\ell=1}^{k}u_{d\ell}{\bm{\beta}}_{d}^{T}{\bm{X}}^{T}{\bm{X}}{\bm{\beta}}_{d}-{\bm{\beta}}_{d}^{T}\gamma\sum_{\ell=1}^{k}u_{d\ell}{\bm{X}}^{T}{\bm{X}}\frac{1}{q_{\ell}}u_{d\ell}{\bm{\beta}}_{d}-{\bm{\beta}}_{d}^{T}2\gamma\sum_{\ell=1}^{k}u_{d\ell}{\bm{X}}^{T}{\bm{X}}\frac{1}{q_{\ell}}\sum_{m\neq
d}^{q}u_{m\ell}{\bm{\beta}}_{m}$
$\displaystyle+{\bm{\beta}}_{d}^{T}\frac{\gamma}{2}\sum_{\ell=1}^{k}u_{d\ell}{\bm{X}}^{T}{\bm{X}}\frac{1}{q_{\ell}}u_{d\ell}{\bm{\beta}}_{d}+{\bm{\beta}}_{d}^{T}\gamma\sum_{\ell=1}^{k}u_{d\ell}{\bm{X}}^{T}{\bm{X}}\frac{1}{q_{\ell}}\sum_{m\neq
d}^{q}u_{m\ell}{\bm{\beta}}_{m}+C_{2}$ $\displaystyle=$
$\displaystyle{\bm{\beta}}_{d}^{T}\left(\frac{\gamma}{2}\sum_{\ell=1}^{k}u_{d\ell}{\bm{X}}^{T}{\bm{X}}\right){\bm{\beta}}_{d}-{\bm{\beta}}_{d}^{T}\left(\frac{\gamma}{2}\sum_{\ell=1}^{k}u_{d\ell}{\bm{X}}^{T}{\bm{X}}\frac{1}{q_{\ell}}u_{d\ell}\right){\bm{\beta}}_{d}$
$\displaystyle-{\bm{\beta}}_{d}^{T}\left(\gamma\sum_{\ell=1}^{k}u_{d\ell}{\bm{X}}^{T}{\bm{X}}\frac{1}{q_{\ell}}\sum_{m\neq
d}^{q}u_{m\ell}{\bm{\beta}}_{m}\right)+C_{2}$ (31)

where $C_{2}$ is the constant value that is not relevant to
${\bm{\beta}}_{d}$.

Here, Eq. (31) is differentiated by ${\bm{\beta}}_{d}\;(d=1,2,\cdots,q)$, and
we have

$\displaystyle\frac{\partial
M_{3}^{*}(\\{{\bm{\beta}}_{s}\\}_{s=1}^{q})}{\partial{\bm{\beta}}_{d}}=\gamma\sum_{\ell=1}^{k}u_{d\ell}{\bm{X}}^{T}{\bm{X}}{\bm{\beta}}_{d}-\gamma\sum_{\ell=1}^{k}\frac{1}{q_{\ell}}u_{d\ell}{\bm{X}}^{T}{\bm{X}}{\bm{\beta}_{d}}-\gamma\sum_{\ell=1}^{k}u_{d\ell}\frac{1}{q_{\ell}}\sum_{m\neq
d}^{q}u_{m\ell}{\bm{X}}^{T}{\bm{X}}{\bm{\beta}}_{d}$ (32)

From the above, with Eq. (23), Eq. (25), and Eq. (32), we have the updated
formula of ${\bm{\beta}}_{d}$ as follows:

$\displaystyle\frac{\partial
M^{\dagger}({\bm{\beta}}_{s})}{\partial{\bm{\beta}}_{d}}=$
$\displaystyle-2\left(\sum_{o=1}^{n(n-1)/2}w_{od}g_{od}{\bm{r}}_{o}\right)+2\left({\sum_{o=1}^{n(n-1)/2}w_{od}{\bm{r}}_{o}{\bm{r}}_{o}^{T}}\right){\bm{\beta}}_{d}+2\left(\lambda{\bm{\Psi}}_{d}^{T}{\bm{\Psi}}_{d}\right){\bm{\beta}}_{d}$
$\displaystyle+\gamma\sum_{\ell=1}^{k}u_{d\ell}{\bm{X}}^{T}{\bm{X}}{\bm{\beta}}_{d}-\gamma\sum_{\ell=1}^{k}\frac{1}{q_{\ell}}u_{d\ell}{\bm{X}}^{T}{\bm{X}}{\bm{\beta}_{d}}-\gamma\sum_{\ell=1}^{k}u_{d\ell}\frac{1}{q_{\ell}}u_{d\ell}\sum_{m\neq
d}^{q}u_{m\ell}{\bm{X}}^{T}{\bm{X}}{\bm{\beta}}_{d}={\bm{0}}_{p}$
$\displaystyle\Longleftrightarrow$
$\displaystyle\left(2\left({\sum_{o=1}^{n(n-1)/2}w_{os}{\bm{r}}_{o}{\bm{r}}_{o}^{T}}\right)+2\left(\lambda{\bm{\Psi}}_{d}^{T}{\bm{\Psi}}_{d}\right)+\gamma\sum_{\ell=1}^{k}u_{d\ell}{\bm{X}}^{T}{\bm{X}}-\gamma\sum_{\ell=1}^{k}\frac{1}{q_{\ell}}u_{d\ell}{\bm{X}}^{T}{\bm{X}}\right){\bm{\beta}}_{d}$
$\displaystyle=2\left(\sum_{o=1}^{n(n-1)/2}w_{od}g_{od}{\bm{r}}_{o}\right)+\gamma\sum_{\ell=1}^{k}\frac{1}{q_{\ell}}u_{d\ell}\sum_{m\neq
d}^{q}u_{m\ell}{\bm{X}}^{T}{\bm{X}}{\bm{\beta}}_{m}$
$\displaystyle\Longleftrightarrow{\bm{\beta}}_{d}$
$\displaystyle=\left(2\left({\sum_{o=1}^{n(n-1)/2}w_{os}{\bm{r}}_{o}{\bm{r}}_{o}^{T}}\right)+2\left(\lambda{\bm{\Psi}}_{d}^{T}{\bm{\Psi}}_{d}\right)+\gamma\sum_{\ell=1}^{k}u_{d\ell}{\bm{X}}^{T}{\bm{X}}-\gamma\sum_{\ell=1}^{k}\frac{1}{q_{\ell}}u_{d\ell}{\bm{X}}^{T}{\bm{X}}\right)^{-1}$
$\displaystyle\left(2\left(\sum_{o=1}^{n(n-1)/2}w_{os}g_{os}{\bm{r}}_{o}\right)+\gamma\sum_{\ell=1}^{k}\frac{1}{q_{\ell}}u_{d\ell}\sum_{m\neq
d}^{q}u_{m\ell}{\bm{X}}^{T}{\bm{X}}{\bm{\beta}}_{m}\right)$ (33)

As Proposition ${\rm{2}}$, because ${\bm{\beta}}_{s}$ is updated in the order,
function $f(m,s)$ is defined as $t+1$ when $m<s$, whereas it is defined as $t$
when $m>s$.

In Algorithm 1, when updating ${\bm{\beta}}_{s}$, $u_{s\ell}$ is fixed as
$u_{s\ell}^{t^{{\dagger}}}$. Therefore, the updated formula of
${\bm{\beta}}_{s}$ is obtained as Eq. (21). ∎

Next, we explain the updated rule of ${\bm{U}}$ and the updated formula of
${\bm{v}}_{\ell}\ (\ell=1,2,\cdots,k)$. ${\bm{U}}$ and ${\bm{v}}_{\ell}$ are
updated in the same manner as the $k$-means.

Update ${\bm{U}}$

${\bm{U}}$ is updated for each $s$. $u_{s\ell^{\dagger}}$ is updated as

$\displaystyle
u_{s\ell^{\dagger}}^{(t^{\dagger}+1)}\leftarrow\begin{cases}1&\big{(}\|{\bm{X}}{\bm{\beta}}_{s}^{(t)}-{\bm{X}}{\bm{v}}_{\ell^{\dagger}}^{(t^{\dagger})}\|_{F}^{2}\leq\|{\bm{X}}{\bm{\beta}}_{s}^{(t)}-{\bm{X}}{\bm{v}}_{\ell}^{(t^{\dagger})}\|_{F}^{2}\
\rm{for}\ \rm{any}\ \ell\ (\ell=1,2,\cdots,k)\big{)}\\\
0&({\rm{otherwise}})\end{cases}$ (34)

which is applied to all $s\;(s=1,2,\cdots,q)$ and $\ell^{\dagger}\
(\ell^{\dagger}=1,2,\cdots,k$).

Update ${\bm{v}}_{\ell}\ (\ell=1,2,\cdots,k)$

${\bm{v}}_{\ell}$ is updated as

$\displaystyle{\bm{v}}_{\ell}^{(t^{\dagger}+1)}\leftarrow\frac{1}{q_{\ell}}\sum_{s=1}^{q}u_{s\ell}^{(t^{\dagger}+1)}{\bm{\beta}}_{s}^{(t)}\quad(\ell=1,2,\cdots,k).$
(35)

Our proposed method updates ${\bm{\beta}}_{s}$, auxiliary variables
${\bm{w}}_{os}$ and ${\bm{\Psi}}_{s}$, indicator matrix ${\bm{U}}$, and
cluster centroid coefficient $\bm{v}_{\ell}$ alternatively based on the MM
algorithm. The details of the MM algorithm are described in Algorithm 1.

Algorithm 1 Wilcoxon-type MCEN based on MM algorithm

0: ${\bm{X}},{\bm{Y}},k,\;\lambda>0,\gamma>0,\;\epsilon>0;$ threshold for this
algorithm $\epsilon^{*}>0$

0:
${\bm{\beta}}_{s}(s=1,2,\cdots,q),{\bm{U}},{\bm{v}}_{\ell}(\ell=1,2,\cdots,k)$

1: Set $t\leftarrow 1$ and $t^{{\dagger}}\leftarrow 1$

2: for $s=1$ to $q$ do

3: Set initial values ${\bm{\beta}_{s}^{(t)}}$

4: end for

5: Set initial ${\bm{U}}^{(t^{{\dagger}})}$ and
${\bm{v}}_{\ell}^{(t^{{\dagger}})}(\ell=1,2,\cdots,k)$ by $k$-means

6: while
$L^{{\dagger}}(\\{{\bm{\beta}}_{s}^{(t)}\\}_{s=1}^{q},{\bm{U}}^{(t^{{\dagger}})},\\{{\bm{v}}_{\ell}^{(t^{{\dagger}})}\\}_{\ell=1}^{k})-L^{{\dagger}}(\\{{\bm{\beta}}_{s}^{(t+1)}\\}_{s=1}^{q},{\bm{U}}^{(t^{{\dagger}}+1)},\\{{\bm{v}}_{\ell}^{(t^{{\dagger}}+1)}\\}_{\ell=1}^{k})\geqq\epsilon^{*}$
do

7: while
$L^{{\dagger}}(\\{{\bm{\beta}}_{s}^{(t)}\\}_{s=1}^{q},{\bm{U}}^{(t^{{\dagger}})},\\{{\bm{v}}_{\ell}^{(t^{{\dagger}})}\\}_{\ell=1}^{k})-L^{{\dagger}}(\\{{\bm{\beta}}_{s}^{(t+1)}\\}_{s=1}^{q},{\bm{U}}^{(t^{{\dagger}})},\\{{\bm{v}}_{\ell}^{(t^{{\dagger}})}\\}_{\ell=1}^{k})\geqq\epsilon^{*}$
do

8: for $s=1$ to $q$ do

9: for $o=1$ to $(n(n-1))/2$ do

10: $w_{os}^{(t+1)}\leftarrow
1/(2|g_{os}-{\bm{r}}_{o}^{T}{\bm{\beta}}_{s}^{(t)}|)$

11: end for

12: ${\bm{\Psi}}_{s}^{(t+1)}\leftarrow{\rm
diag}\left(\frac{1}{\sqrt{2(\beta_{s1}^{(t)}+\epsilon)}},\frac{1}{\sqrt{2(\beta_{s2}^{(t)}+\epsilon)}},\cdots,\frac{1}{\sqrt{2(\beta_{sp}^{(t)}+\epsilon)}}\right)$

13:
${\bm{\beta}}_{s}^{(t+1)}\leftarrow\left(2\left({\sum_{o=1}^{n(n-1)/2}w_{os}^{(t+1)}{\bm{r}}_{o}{\bm{r}}_{o}^{T}}\right)+2\left(\lambda{\bm{\Psi}}_{s}^{(t+1)T}{\bm{\Psi}}_{s}^{(t+1)}\right)+\gamma\sum_{\ell=1}^{k}u_{s\ell}^{(t^{{\dagger}})}{\bm{X}}^{T}{\bm{X}}-\gamma\sum_{\ell=1}^{k}\frac{1}{q_{\ell}}u_{s\ell}^{(t^{{\dagger}})}{\bm{X}}^{T}{\bm{X}}\right)^{-1}$
$\times\left(2\left(\sum_{o=1}^{n(n-1)/2}w_{os}^{(t+1)}g_{os}{\bm{r}}_{o}\right)+\gamma\sum_{\ell=1}^{k}\frac{1}{q_{\ell}}u_{s\ell}^{(t^{{\dagger}})}\sum_{m\neq
s}^{q}u_{m\ell}^{(t^{{\dagger}})}{\bm{X}}^{T}{\bm{X}}{\bm{\beta}}_{m}^{(f(m,s))}\right)$

14: end for

15: $t\leftarrow t+1$

16: end while

17: while
$L^{{\dagger}}(\\{{\bm{\beta}}_{s}^{(t)}\\}_{s=1}^{q},{\bm{U}}^{(t^{{\dagger}})},\\{{\bm{v}}_{\ell}^{(t^{{\dagger}})}\\}_{\ell=1}^{k})-L^{{\dagger}}(\\{{\bm{\beta}}_{s}^{(t)}\\}_{s=1}^{q},{\bm{U}}^{(t^{{\dagger}}+1)},\\{{\bm{v}}_{\ell}^{(t^{{\dagger}}+1)}\\}_{\ell=1}^{k})\geqq\epsilon^{*}$
do

18: for $s=1$ to $q$ do

19: Update ${\bm{U}}^{(t^{{\dagger}})}$ to ${\bm{U}}^{(t^{{\dagger}}+1)}$
based on Eq. (34)

20: end for

21: for $\ell=1$ to $k$ do

22: Update ${\bm{v}}_{\ell}^{(t^{{\dagger}})}$ to
${\bm{v}}_{\ell}^{(t^{{\dagger}}+1)}$ based on Eq. (35)

23: end for

24: $t^{{\dagger}}\leftarrow t^{{\dagger}}+1$

25: end while

26: end while

##  3 Numerical Simulation

In this section, we demonstrate the efficiency of the proposed method through
a numerical simulation.

###  3.1 Simulation design

This simulation generates covariate matrix ${\bm{X}}$ and response variables
${\bm{Y}}$ with a true coefficient matrix. We show the simulation setting
patterns in Table 1. The simulation evaluates the prediction accuracy and mean
squared error (MSE) of the coefficient matrix. In this simulation, we use
RStudio Version 1.4.1103.

We generate the simulation data, adjusted based on Price and Sherwood (2018).
The sample size is $50$, and the number of response variables is $9$. We set
the number of covariate matrix $p$ as $12$ and $100$, which will be discussed
in Factor 2. First, let

$\displaystyle\tilde{\bm{\Sigma}}=(\tilde{\sigma}_{j^{\prime}j}),\tilde{\sigma}_{jj}=1,\tilde{\sigma}_{j^{\prime}j}=0.7,(j^{\prime},j=1,2,\cdots,12).$
(36)

The covariate vectors are generated by
${\bm{X}}_{i}\sim\mathcal{N}({\bm{0}}_{12},\tilde{\bm{\Sigma}})$ for $p=12$.
In the case of $p=100$,
${\bm{X}}_{i}\sim\mathcal{N}({\bm{0}}_{100},{\bm{\Sigma}}^{\dagger})$, where

$\displaystyle{\bm{\Sigma}}^{\dagger}=\left(\begin{array}[]{ll}\tilde{\bm{\Sigma}}&{\bm{O}_{12,88}}\\\
{\bm{O}_{88,12}}&{\bm{I}}_{88}\end{array}\right).$ (39)

${\bm{O}_{12,88}}$ is a zero matrix with $12$ rows and $88$ columns,
${\bm{O}_{88,12}}$ is a zero matrix with $88$ rows and $12$ columns, and
${\bm{I}}_{88}$ is the identity matrix.

To make the matrix of ${\bm{B}}$, in the case of $p=12$, we set the
${\bm{b}}_{4}(\eta,\xi)=({\bm{\eta}}_{4}-\xi,{\bm{\eta}}_{4},{\bm{\eta}}_{4}+\xi)\in\mathbb{R}^{4\times
3}$. ${\bm{\eta}}_{K}$ is a vector with the length of $K$, whose elements are
all $\eta$. Details on how the values of $\eta$ and $\xi$ are determined will
be explained in Factor $3$ and Factor $4$, respectively.

Here, the true ${\bm{B}}^{*}_{\eta,\xi}$ in $p=12$ is set as follows:

$\displaystyle{\bm{B}}^{*}_{\eta,\xi}=\left(\begin{array}[]{lll}{\bm{b}}_{4}(\eta,\xi)&{\bm{O}}_{4,3}&{\bm{O}}_{4,3}\\\
{\bm{O}}_{4,3}&{\bm{b}}_{4}(\eta,\xi)&{\bm{O}}_{4,3}\\\
{\bm{O}}_{4,3}&{\bm{O}}_{4,3}&{\bm{b}}_{4}(\eta,\xi)\end{array}\right).$ (43)

For the coefficient groups in $p=100$, as is the case of $p=12$,
${\bm{b}}_{10}(\eta,\xi)=({\bm{\eta}}_{10}-\xi,{\bm{\eta}}_{10},{\bm{\eta}}_{10}+\xi)\in\mathbb{R}^{10\times
3}$.

$\displaystyle{\bm{B}}^{*}_{\eta,\xi}=\left(\begin{array}[]{lll}{\bm{b}}_{10}(\eta,\xi)&{\bm{O}}_{10,3}&{\bm{O}}_{10,3}\\\
{\bm{O}}_{10,3}&{\bm{b}}_{10}(\eta,\xi)&{\bm{O}}_{10,3}\\\
{\bm{O}}_{10,3}&{\bm{O}}_{10,3}&{\bm{b}}_{10}(\eta,\xi)\\\
{\bm{O}}_{70,3}&{\bm{O}}_{70,3}&{\bm{O}}_{70,3}\\\ \end{array}\right).$ (48)

Next, we explain how to generate the response variable, with ${\bm{X}}$,
${\bm{B}}$. The response variable is generated as follows:

$\displaystyle{\bm{y}}_{i}={{\bm{B}}^{*}_{\eta,\xi}}^{T}{\bm{X}}_{i}+{\bm{\varepsilon}}_{i}.$
(49)

Here, ${\bm{\varepsilon}}_{i}$ is the random error. Several different settings
are employed, which are described in detail in Factor 5. The number of total
patterns in this simulation setting is $4$ (Factor 1) $\times 2$ (Factor 2)
$\times 4$ (Factor 3) $\times 3$ (Factor 4) $\times 4$ (Factor 5) $=384$. For
each pattern combination, we generate $50$ learning samples and $1000$ test
samples from Eq. (49), and we repeat the calculation $100$ times.

Next, we explain the evaluation index. In this simulation, we evaluate the
prediction variable and coefficient matrix. The median of the absolute
prediction error (APE) is

$\displaystyle{\rm median}\big{\\{}|y^{*}_{is}-\hat{y}_{is}|,\
i=1,2,\cdots,1000;s=1,2,\cdots,9\big{\\}}$

where $y_{is}^{*}$ represents the testing samples and $\hat{y}_{is}$
represents the prediction variables. We also compare the mean squared error of
estimator

$\displaystyle\frac{1}{p\times
9}\sum_{s=1}^{9}\|{\hat{\bm{\beta}}}_{s}-{\bm{\beta}}_{s}^{*}\|^{2}_{2}$

where ${\bm{\beta}}_{s}^{*}$ is the true ${\bm{\beta}}_{s}$, and
${\hat{\bm{\beta}}_{s}}$ are the predicted coefficient vectors.

Table 1: Factors of the simulation design Factor No. | Factor name | level
---|---|---
Factor 1 | Method | 4
Factor 2 | Number of covariate variable | 2
Factor 3 | Number of parameter $\eta$ for true coefficient matrix | 4
Factor 4 | Number of parameter $\xi$ for true coefficient matrix | 3
Factor 5 | Number of error distribution | 4

Now, we explain each simulation setting factor denoted in Table 1.

Factor 1: Method

We apply the $4$ methods to compare their performance.

In addition to the proposed method, we apply three methods for comparison: the
multivariate cluster elastic net (MCEN) (Price and Sherwood, 2018), separate
elastic net (Zou and Hastie, 2005), and Wilcoxon-type lasso (Wang et al.,
2020). The optimization problem of MCEN is shown as follows:

$\displaystyle\frac{1}{2n}\|{\bm{Y}}-{\bm{X}}{\bm{B}}\|_{F}^{2}+\lambda\sum^{q}_{s=1}\|{\bm{\beta}}_{s}\|_{1}+\frac{\gamma}{2}\sum_{s=1}^{q}\sum_{\ell=1}^{k}u_{s\ell}\|{\bm{X}}{\bm{\beta}}_{s}-{\bm{X}}{\bm{v}}_{\ell}\|_{2}^{2}$

where $k$ is the number of clusters, and $\lambda\ (\lambda>0)$ and $\gamma\
(\gamma>0)$ are the tuning parameters. ${\bm{U}}=(u_{s\ell})\
(s=1,2,\cdots,q;\ell=1,2,\cdots,k)$ is the indicator matrix representing the
degree of belonging to cluster $\ell$ on the response variable $s$.
$v_{\ell}\;(\ell=1,2,\cdots,k$) is the partial regression coefficient for the
cluster centroid. Elastic net is applied to each response variable separately,
denoted as SEN. The objective function of elastic net is as follows:

$\displaystyle\|{\bm{y}}_{s}-{\bm{X}}{\bm{\beta}}_{s}\|^{2}_{2}+\delta_{s}\|{\bm{\beta}}_{s}\|_{1}+\gamma_{s}\|{\bm{\beta}}_{s}\|^{2}_{2}\quad(s=1,2,\cdots,q)$

where $\delta_{s}\;(\delta_{s}>0)$ and $\gamma_{s}\;(\gamma_{s}>0)$ are tuning
parameters. The third compared method is lasso regression, for which the loss
function is a Wilcoxon-type regression function. It is denoted by WLASSO. The
objective function of WLASSO is as follows:

$\displaystyle\sum_{i<j}|e_{is}-e_{js}|+\gamma_{s}\|{\bm{\beta}}_{s}\|_{1}\quad(s=1,2,\cdots,q)$

where
$e_{is}=y_{is}-\bm{x}_{i}^{T}\bm{\beta}_{s},e_{js}=y_{js}-\bm{x}_{j}^{T}\bm{\beta}_{s}$,
and $\gamma_{s}\;(\gamma_{s}>0)$ is the tuning parameter.

MCEN uses the mcen package in R software (Sherwood and Price, 2020), and SEN
are fitted by the glmnet package in R software (Friedman et al., 2008). For
WLASSO, we use LADlasso (Wang et al., 2007) in the package called MTE (Li and
Qin, 2021), as mentioned in Wang et al. (2020). The tuning parameters of all
methods are chosen by five-fold cross validation, and the cluster sizes of $2$
and $3$ are candidates for the proposed method and MCEN.

Factor 2: Covariate variable

The number of the covariate vectors is set as $p=12$ and $100$. The covariate
matrix for $p=12$ is described in Eq. (36), and it is described for $p=100$ in
Eq. (39).

Factor 3: Parameter $\eta$ for true coefficient matrix
${\bm{B}}^{*}_{\eta,\xi}$

$\eta$ is a factor to control each value of the non-zero part of a true
coefficient matrix. The candidates for $\eta$, one of the parameters for the
true coefficient matrix, are set as $0.25,0.5,0.75$, and $1$.

Factor 4: Parameter $\xi$ for true coefficient matrix
${\bm{B}}^{*}_{\eta,\xi}$

$\xi$ is a factor to control the difference between each non-zero true
coefficient vector in ${\bm{B}}^{*}_{\eta,\xi}$. The candidates for $\xi$, the
other parameter for the true coefficient matrix, are set as $0.02,0.05$, and
$0.10$. With $\eta$ and $\xi$, the true coefficient matrix generates
${\bm{B}}^{*}_{\eta,\xi}$ in Eq. (43) and Eq. (48).

Factor 5: Error distribution

We set four different patterns for ${\bm{\varepsilon}}_{i}=(\varepsilon_{is})$
based on Wang et al. (2020).

Error 1 is set as

$\displaystyle\varepsilon_{is}\sim N(0,1)\ (i=1,2,\cdots,n;s=1,2,\cdots,q).$

Error 2 is set as

$\displaystyle\varepsilon_{is}\sim 0.95N(0,1)+0.05N(0,100)\
(i=1,2,\cdots,n;s=1,2,\cdots,q).$

Error 3 is set as

$\displaystyle\varepsilon_{is}\sim\sqrt{2}t(4)\
(i=1,2,\cdots,n;s=1,2,\cdots,q)$

where $t(4)$ represents the $t$ distribution with $4$ degree of freedom.

Error 4 is set as

$\displaystyle\varepsilon_{is}\sim Cauchy(0,1)\
(i=1,2,\cdots,n;s=1,2,\cdots,q)$

where $Cauchy(0,1)$ represents a Cauchy distribution with a location parameter
of $0$ and a scale parameter of $1$. Error 1 assumes a normal distribution,
Error 2 assumes a distribution containing outliers, and Error 3 and Error 4
assume heavy-tailed distribution.

###  3.2 Simulation result

The results of the simulation are displayed by error distribution in Tables 2
to 5 for $p=12$, and Tables 6 to 9 for $p=100$. First, as seen by each outcome
for Factor 1, the results of the proposed method were better than the those of
the other compared methods in terms of both evaluation indices in all patterns
of ($\eta$, $\xi$) in Error 2 for $p=100$ (Table 7) and Error 4 for both
$p=12$ (Table 5) and $100$ (Table 9). Focusing on the median APE, the proposed
method was better in all patterns of ($\eta$, $\xi$) in Error 3 for $p=12$
(Table 4); Error 1 for $p=12$ (Table 2), except $\eta=0.25$; and Error 2 of
$p=12$ (Table 3), except $(\eta,\xi)=(0.75,0.05)$. In terms of the MSE of
${\bm{\beta}}_{s}$ (the second evaluation index), the proposed method had the
smallest values in all patterns except $(\eta,\xi)=(0.75,0.05)$ in Error 2 for
$p=12$ (Table 3); and Error 3 for $p=12$ except $\eta=0.25$. For the compared
methods, the MCEN had smaller values in almost all patterns of Error 3 for
$p=100$ (Table 8) in both evaluation indices. For median APE, the MCEN was
better in several patterns of Error 1 for $p=100$ (Table 6) as well as in
Error 1 for $p=12$ (Table 2) and $100$ (Table 6) for the MSE of
${\bm{\beta}}_{s}$. SEN was better for several patterns of $\eta=0.25$ of
Error 1 for $p=12$ (Table 2) for both evaluations indices as well as in Error
1 for $p=100$ (Table 6) for the median APE. WLASSO was only better than the
other methods in $\eta=0.75$ and $\xi=0.05$ in Error 2 for $p=12$ (Table 3).

Next, observing in the covariate variables of Factor 2 on the median of APE,
the proposed method was superior to SEN in $\eta=0.50,0.75,$ and $1.00$ in
$p=12$ for Error 1. However, the values of MCEN and SEN were smaller than the
proposed method in the same pattern for $p=100$. In Error 3 for $p=12$, the
proposed method outperformed MCEN for most patterns of both evaluation
measures, while MCEN was better for $p=100$. In Error 1 and Error 3, the
values of the proposed method and MCEN were close each other, regardless of
the value of the covariate variables. The results on Error 2 and Error 4 did
not vary depending on the covariate variables. As for Factor 3, parameter
$\eta$, the values of all methods increased in proportion to the augment of
$\eta$. Meanwhile, the values did not tend to change depending on $\xi$,
Factor $4$.

Finally, we see the results of the error distribution, Factor 5. In Error 1,
the proposed method was almost stable for $p=12$ in the median APE, while SEN
and MCEN were better than the proposed method for $p=100$. MCEN had the
smallest MSE of ${\bm{\beta}}_{s}$ for both $p=12$ and $p=100$. However, the
value differences among the methods were slight in both the median APE and the
MSE of ${\bm{\beta}}_{s}$. In the case of Error 2, the situation where the
error contains outliers, for both $p=12$ and $p=100$, the proposed method was
better than the other compared methods in both the median APE and the MSE of
${\bm{\beta}}_{s}$. Focusing on the differences, the value differences between
the proposed method and WLASSO were closer than those between the proposed
method and MCEN and SEN. In Error 3, whose error distribution follows $t$
distribution, the median APE of the proposed method was better in
$(\eta,\xi)=(0.25,0.05),(0.25,0.10)$, and $(0.50,0.10)$, and MCEN had better
results in the other pattern of $\eta$ and $\xi$. In the MSE of
${\bm{\beta}}_{s}$, the proposed method was the smallest for all for $p=12$,
except $\eta=0.25$. For $p=100$, MCEN was smaller than the proposed method;
however, the differences were minor. Finally, in Error 4, the other heavy-
tailed error distribution, the proposed method was better than all the other
compared methods. For $p=12$ in Error 4, WLASSO was also stable in terms of
both the median of APE and MSE of ${\bm{\beta}}_{s}$ compared to MCEN and SEN.

Table 2: Result of the simulation $p=12$ for Error $1$

|  | the median of APE (sd of APE) | MSE of ${\bm{\beta}}_{s}$ (sd of MSE)
---|---|---|---
| $\xi$ | 0.02 | 0.05 | 0.10 | 0.02 | 0.05 | 0.10
$\eta$ |  |  |  |  |  |  |
0.25 | WMCEN | 0.706 (0.010) | 0.706 (0.010) | 0.705 (0.011) | 0.023 (0.004) | 0.023 (0.004) | 0.021 (0.004)
| MCEN | 0.712 (0.005) | 0.711 (0.014) | 0.715 (0.014) | 0.019 (0.004) | $\bm{0.017(0.004)}$ | 0.019 (0.003)
| SEN | $\bm{0.703(0.010)}$ | $\bm{0.705(0.010)}$ | $\bm{0.705(0.010)}$ | $\bm{0.018(0.003)}$ | 0.019 (0.004) | $\bm{0.018(0.004)}$
| WLASSO | 0.720 (0.012) | 0.720 (0.012) | 0.714 (0.011) | 0.035 (0.007) | 0.034 (0.007) | 0.029 (0.005)
0.50 | WMCEN | $\bm{0.720(0.013)}$ | $\bm{0.724(0.013)}$ | $\bm{0.721(0.012)}$ | 0.035 (0.007) | 0.037 (0.007) | 0.035 (0.006)
| MCEN | 0.741 (0.034) | 0.741 (0.033) | 0.746 (0.033) | $\bm{0.029(0.008)}$ | $\bm{0.029(0.008)}$ | $\bm{0.032(0.007)}$
| SEN | 0.725 (0.012) | $\bm{0.724(0.012)}$ | 0.724 (0.012) | 0.035 (0.006) | 0.036 (0.006) | 0.035 (0.007)
| WLASSO | 0.730 (0.013) | 0.735 (0.014) | 0.734 (0.014) | 0.043 (0.008) | 0.047 (0.008) | 0.046 (0.008)
0.75 | WMCEN | $\bm{0.724(0.014)}$ | $\bm{0.724(0.013)}$ | $\bm{0.736(0.014)}$ | 0.038 (0.008) | 0.038 (0.008) | 0.039 (0.008)
| MCEN | 0.764 (0.062) | 0.765 (0.062) | 0.766 (0.062) | $\bm{0.028(0.009)}$ | $\bm{0.028(0.009)}$ | $\bm{0.041(0.009)}$
| SEN | 0.735 (0.013) | 0.734 (0.013) | 0.729 (0.013) | 0.046 (0.009) | 0.043 (0.009) | 0.041 (0.009)
| WLASSO | 0.739 (0.015) | 0.739 (0.015) | 0.739 (0.015) | 0.051 (0.0010) | 0.051 (0.010) | 0.051 (0.010)
1.00 | WMCEN | $\bm{0.725(0.014)}$ | $\bm{0.725(0.015)}$ | $\bm{0.726(0.014)}$ | 0.039 (0.009) | 0.040(0.009) | 0.052 (0.009)
| MCEN | 0.805 (0.100) | 0.803 (0.098) | 0.803 (0.099) | $\bm{0.034(0.013)}$ | $\bm{0.034(0.013)}$ | $\bm{0.034(0.013)}$
| SEN | 0.735 (0.014) | 0.735 (0.015) | 0.737 (0.014) | 0.048 (0.010) | 0.045 (0.010) | 0.049 (0.010)
| WLASSO | 0.740 (0.015) | 0.740 (0.015) | 0.740 (0.015) | 0.052 (0.010) | 0.052 (0.010) | 0.052 (0.010)

Table 3: Result of the simulation $p=12$ for Error $2$

|  | the median of APE (sd of APE) | MSE of ${\bm{\beta}}_{s}$ (sd of MSE)
---|---|---|---
| $\xi$ | 0.02 | 0.05 | 0.10 | 0.02 | 0.05 | 0.10
$\eta$ |  |  |  |  |  |  |
0.25 | WMCEN | $\bm{0.756(0.014)}$ | $\bm{0.757(0.014)}$ | $\bm{0.757(0.014)}$ | $\bm{0.025(0.006)}$ | $\bm{0.026(0.006)}$ | $\bm{0.026(0.006)}$
| MCEN | 2.712 (1.039) | 2.622 (1.020) | 2.667 (1.024) | 1.026 (1.285) | 0.567 (0.655) | 0.772 (1.004)
| SEN | 2.298 (1.175) | 3.074 (1.751) | 2.884 (1.618) | 11.18 (11.689) | 17.530 (16.730) | 15.666 (15.252)
| WLASSO | 0.778 (0.017) | 0.778 (0.019) | 0.778 (0.019) | 0.046 (0.013) | 0.045 (0.013) | 0.045 (0.013)
0.50 | WMCEN | $\bm{0.778(0.017)}$ | $\bm{0.779(0.017)}$ | $\bm{0.779(0.018)}$ | $\bm{0.043(<0.001)}$ | $\bm{0.044(0.010)}$ | $\bm{0.044(0.010)}$
| MCEN | 2.742 (1.061) | 2.748 (1.065) | 2.761 (1.074) | 0.822 (0.699) | 0.858 (0.734) | 0.953 (0.822)
| SEN | 3.403 (1.790) | 3.662 (2.013) | 3.286 (1.739) | 18.138 (16.549) | 20.243 (18.213) | 17.781 (16.601)
| WLASSO | 0.808 (0.024) | 0.808 (0.024) | 0.808 (0.024) | 0.070 (0.020) | 0.070 (0.020) | 0.070 (0.020)
0.75 | WMCEN | $\bm{0.794(0.024)}$ | 0.833 (0.031) | $\bm{0.791(0.021)}$ | $\bm{0.057(0.017)}$ | 0.093 (0.028) | $\bm{0.054(0.015)}$
| MCEN | 2.910 (1.093) | 2.958 (1.138) | 2.962 (1.142) | 1.730 (1.634) | 2.086 (1.900) | 2.135 (1.946)
| SEN | 3.498 (1.736) | 3.662 (1.890) | 3.538 (1.783) | 19.120 (17.141) | 19.014 (16.954) | 18.692 (16.883)
| WLASSO | 0.814 (0.026) | $\bm{0.814(0.026)}$ | 0.814 (0.026) | 0.075 (0.022) | $\bm{0.075(0.022)}$ | 0.075 (0.022)
1.00 | WMCEN | $\bm{0.791(0.026)}$ | $\bm{0.796(0.026)}$ | $\bm{0.793(0.026)}$ | $\bm{0.055(0.020)}$ | $\bm{0.059(0.021)}$ | $\bm{0.056(0.019)}$
| MCEN | 3.011 (1.118) | 3.01 2(1.118) | 3.021 (1.120) | 2.135 (1.717) | 2.063 (1.736) | 2.096 (1.757)
| SEN | 3.705 (1.910) | 3.868 (2.000) | 3.551 (1.773) | 19.082 (16.757) | 20.503 (17.810) | 17.927 (16.542)
| WLASSO | 0.816 (0.027) | 0.816 (0.028) | 0.816 (0.028) | 0.078 (0.024) | 0.078 (0.024) | 0.078 (0.024)

Table 4: Result of the simulation $p=12$ for Error $3$

|  | the median of APE (sd of APE) | MSE of ${\bm{\beta}}_{s}$ (sd of MSE)
---|---|---|---
| $\xi$ | 0.02 | 0.05 | 0.10 | 0.02 | 0.05 | 0.10
$\eta$ |  |  |  |  |  |  |
0.25 | WMCEN | $\bm{1.097(0.018)}$ | $\bm{1.101(0.019)}$ | $\bm{1.100(0.018)}$ | 0.039 (0.008) | 0.042 (0.009) | 0.040 (0.009)
| MCEN | 1.126 (0.022) | 1.122 (0.022) | 1.122 (0.021) | $\bm{0.035(0.010)}$ | $\bm{0.031(0.008)}$ | $\bm{0.030(0.008)}$
| SEN | 1.117 (0.023) | 1.107 (0.022) | 1.121 (0.024) | 0.043 (0.014) | 0.037 (0.011) | 0.052 (0.016)
| WLASSO | 1.121 (0.021) | 1.122 (0.020) | 1.121 (0.021) | 0.061 (0.015) | 0.062 (0.015) | 0.060 (0.015)
0.50 | WMCEN | $\bm{1.128(0.019)}$ | $\bm{1.128(0.019)}$ | $\bm{1.126(0.019)}$ | $\bm{0.069(0.011)}$ | $\bm{0.069(0.011)}$ | $\bm{0.067(0.010)}$
| MCEN | 1.171 (0.029) | 1.169 (0.029) | 1.172 (0.030) | 0.070 (0.016) | 0.071 (0.017) | 0.070 (0.016)
| SEN | 1.150 (0.026) | 1.155 (0.027) | 1.146 (0.024) | 0.083 (0.023) | 0.086 (0.024) | 0.077 (0.017)
| WLASSO | 1.159 (0.021) | 1.159 (0.022) | 1.158 (0.022) | 0.104 (0.0214) | 0.104 (0.021) | 0.103 (0.021)
0.75 | WMCEN | $\bm{1.150(0.022)}$ | $\bm{1.150(0.022)}$ | $\bm{1.156(0.021)}$ | $\bm{0.093(0.015)}$ | $\bm{0.093(0.015)}$ | $\bm{0.098(0.016)}$
| MCEN | 1.213 (0.048) | 1.213 (0.048) | 1.214 (0.048) | 0.099 (0.022) | 0.099 (0.022) | 0.100 (0.023)
| SEN | 1.189 (0.029) | 1.176 (0.025) | 1.179 (0.028) | 0.121 (0.028) | 0.121 (0.023) | 0.110 (0.026)
| WLASSO | 1.177 (0.023) | 1.178 (0.023) | 1.177 (0.023) | 0.124 (0.024) | 0.124 (0.023) | 0.124 (0.023)
1.00 | WMCEN | $\bm{1.159(0.024)}$ | $\bm{1.163(0.024)}$ | $\bm{1.159(0.022)}$ | $\bm{0.104(0.020)}$ | $\bm{0.108(0.021)}$ | $\bm{0.105(0.019)}$
| MCEN | 1.240 (0.071) | 1.252 (0.071) | 1.252 (0.072) | 0.105 (0.033) | 0.115 (0.030) | 0.116 (0.029)
| SEN | 1.200 (0.029) | 1.200 (0.028) | 1.198 (0.029) | 0.138 (0.028) | 0.139 (0.027) | 0.141 (0.030)
| WLASSO | 1.198 (0.027) | 1.198 (0.027) | 1.198 (0.027) | 0.151 (0.029) | 0.151 (0.029) | 0.151 (0.028)

Table 5: Result of the simulation $p=12$ for Error $4$

|  | the median of APE (sd of APE) | MSE of ${\bm{\beta}}_{s}$ (sd of MSE)
---|---|---|---
| $\xi$ | 0.02 | 0.05 | 0.10 | 0.02 | 0.05 | 0.10
$\eta$ |  |  |  |  |  |  |
0.25 | WMCEN | $\bm{1.134(0.033)}$ | $\bm{1.102(0.027)}$ | $\bm{1.105(0.026)}$ | $\bm{0.072(0.023)}$ | $\bm{0.047(0.013)}$ | $\bm{0.048(0.013)}$
| MCEN | 2.263 (0.629) | 2.264 (0.627) | 2.276 (0.631) | 696.075 (5826.668) | 654.725 (5501.132) | 668.054 (5430.069)
| SEN | 2.365 (0.686) | 2.389 (0.758) | 2.167 (0.636) | 53547.438 (505912.597) | 52569.256 (506213.536) | 53526.572 (505817.297)
| WLASSO | 1.205 (0.054) | 1.198 (0.053) | 1.186 (0.050) | 0.135 (0.062) | 0.130 (0.059) | 0.120 (0.058)
0.50 | WMCEN | $\bm{1.157(0.031)}$ | $\bm{1.157(0.031)}$ | $\bm{1.159(0.031)}$ | $\bm{0.084(0.018)}$ | $\bm{0.084(0.018)}$ | $\bm{0.084(0.018)}$
| MCEN | 2.385 (0.602) | 2.387 (0.607) | 2.386 (0.612) | 809.840 (6028.106) | 798.554 (5957.003) | 843.446 (6447.461)
| SEN | 2.533 (0.658) | 2.526 (0.632) | 2.573 (0.731) | 53541.222 (505771.267) | 53537.523 (505972.399) | 53556.395 (505970.137)
| WLASSO | 1.248 (0.055) | 1.249 (0.056) | 1.253 (0.057) | 0.166 (0.064) | 0.166 (0.064) | 0.169 (0.066)
0.75 | WMCEN | $\bm{1.221(0.041)}$ | $\bm{1.203(0.037)}$ | $\bm{1.204(0.037)}$ | $\bm{0.137(0.034)}$ | $\bm{0.118(0.025)}$ | $\bm{0.118(0.025)}$
| MCEN | 2.490 (0.598) | 2.489 (0.601) | 2.491 (0.604) | 798.148 (6304.464) | 794.190 (6289.904) | 776.378 (6257.876)
| SEN | 2.811 (0.728) | 2.824 (0.749) | 2.748 (0.718) | 53569.628 (506048.462) | 53576.791 (506079.745) | 53492.565 (505240.466)
| WLASSO | 1.296 (0.06) | 1.294 (0.059) | 1.295 (0.058) | 0.208 (0.070) | 0.206 (0.069) | 0.205 (0.068)
1.00 | WMCEN | $\bm{1.249(0.049)}$ | $\bm{1.249(0.048)}$ | $\bm{1.233(0.043)}$ | $\bm{0.160(0.041)}$ | $\bm{0.161(0.042)}$ | $\bm{0.142(0.031)}$
| MCEN | 2.607 (0.600) | 2.602 (0.595) | 2.600 (0.587) | 163.618 (1110.611) | 1439.451 (9552.887) | 1072.825 (8082.635)
| SEN | 2.899 (0.726) | 2.921 (0.719) | 2.964 (0.709) | 53575.183 (506125.472) | 53581.492 (506154.804) | 53582.836 (506128.047)
| WLASSO | 1.300 (0.058) | 1.296 (0.059) | 1.303 (0.058) | 0.201 (0.0504) | 0.197 (0.053) | 0.204 (0.056)

Table 6: Result of the simulation $p=100$ for Error $1$

|  | the median of APE (sd of APE) | MSE of ${\bm{\beta}}_{s}$ (sd of MSE)
---|---|---|---
| $\xi$ | 0.02 | 0.05 | 0.10 | 0.02 | 0.05 | 0.10
$\eta$ |  |  |  |  |  |  |
0.25 | WMCEN | 0.848 (0.016) | 0.844 (0.017) | 0.844 (0.017) | 0.007 (<0.001) | 0.007 (<0.001) | 0.007 (<0.001)
| MCEN | $\bm{0.846(0.027)}$ | 0.848 (0.027) | $\bm{0.808(0.023)}$ | $\bm{0.005(<0.001)}$ | $\bm{0.005(<0.001)}$ | $\bm{0.005(<0.001)}$
| SEN | $\bm{0.846(0.016)}$ | $\bm{0.817(0.014)}$ | 0.826 (0.015) | 0.007 (<0.001) | $\bm{0.005(<0.001)}$ | 0.006 (<0.001)
| WLASSO | 1.111 (0.034) | 1.111 (0.035) | 1.100 (0.035) | 0.020 (0.002) | 0.020 (0.002) | 0.020 (0.002)
0.50 | WMCEN | 0.968 (0.033) | 0.999 (0.032) | 1.002 (0.33) | 0.014 (0.002) | 0.015 (0.002) | 0.015 (0.002)
| MCEN | $\bm{0.926(0.054)}$ | $\bm{0.920(0.054)}$ | $\bm{0.910(0.053)}$ | $\bm{0.009(0.002)}$ | $\bm{0.009(0.002)}$ | $\bm{0.009(0.002)}$
| SEN | 0.979 (0.031) | 0.964 (0.028) | 1.004 (0.029) | 0.014 (0.002) | 0.013 (0.002) | 0.015 (0.002)
| WLASSO | 1.203 (0.040) | 1.197 (0.041) | 1.195 (0.041) | 0.027 (0.003) | 0.027 (0.003) | 0.027 (0.003)
0.75 | WMCEN | 1.078 (0.051) | 1.078 (0.051) | 1.081 (0.050) | 0.021 (0.003) | 0.021 (0.003) | 0.021 (0.003)
| MCEN | 1.068 (0.093) | 1.068 (0.092) | $\bm{1.068(0.093)}$ | $\bm{0.016(0.003)}$ | $\bm{0.016(0.003)}$ | $\bm{0.016(0.003}$
| SEN | $\bm{1.045(0.050)}$ | $\bm{1.044(0.044)}$ | 1.070 (0.046) | 0.019 (0.003) | 0.019 (0.003) | 0.021 (0.003)
| WLASSO | 1.260 (0.057) | 1.260 (0.056) | 1.253 (0.056) | 0.032 (0.004) | 0.032 (0.004) | 0.032 (0.004)
1.00 | WMCEN | 1.180 (0.064) | 1.183 (0.064) | 1.115 (0.067) | 0.028 (0.005) | 0.028 (0.005) | 0.024 (0.005)
| MCEN | 1.165 (0.132) | 1.166 (0.132) | 1.165 (0.131) | $\bm{0.021(0.004)}$ | $\bm{0.021(0.004)}$ | $\bm{0.021(0.004)}$
| SEN | $\bm{1.101(0.066)}$ | $\bm{1.111(0.065)}$ | $\bm{1.088(0.063)}$ | 0.024 (0.005) | 0.026 (0.005) | 0.023 (0.005)
| WLASSO | 1.285 (0.068) | 1.272 (0.068) | 1.283 (0.068) | 0.035 (0.006) | 0.034 (0.005) | 0.035 (0.006)

Table 7: Result of the simulation $p=100$ for Error $2$

|  | the median of APE (sd of APE) | MSE of ${\bm{\beta}}_{s}$ (sd of MSE)
---|---|---|---
| $\xi$ | 0.02 | 0.05 | 0.10 | 0.02 | 0.05 | 0.10
$\eta$ |  |  |  |  |  |  |
0.25 | WMCEN | $\bm{0.871(0.024)}$ | $\bm{0.891(0.024)}$ | $\bm{0.874(0.022)}$ | $\bm{0.005(<0.001)}$ | $\bm{0.006(0.001)}$ | $\bm{0.005(0.001)}$
| MCEN | 5.893 (2.714) | 6.308 (2.863) | 5.290 (2.382) | 1.141 (0.681) | 1.1357 (0.846) | 0.888 (0.524)
| SEN | 10.536 (6.366) | 11.010 (6.620) | 10.019 (6.162) | 5.240 (4.030) | 5.576 (4.226) | 4.762 (3.708)
| WLASSO | 1.718 (0.940) | 1.716 (0.937) | 1.718 (0.942) | 0.170 (0.437) | 0.169 (0.437) | 0.170 (0.437)
0.50 | WMCEN | $\bm{1.087(0.049)}$ | $\bm{1.090(0.052)}$ | $\bm{1.088(0.048)}$ | $\bm{0.016(0.002)}$ | $\bm{0.016(0.002)}$ | $\bm{0.016(0.002)}$
| MCEN | 7.287 (3.411) | 7.295 (3.416) | 7.301 (3.422) | 1.768 (1.085) | 1.772 (1.087) | 1.792 (1.095)
| SEN | 11.232 (6.546) | 10.761 (6.334) | 11.375 (6.601) | 5.545 (4.160) | 5.190 (3.945) | 5.673 (4.2560)
| WLASSO | 2.028 (1.115) | 2.166 (1.365) | 2.23 (1.463) | 0.270 (0.599) | 0.371 (0.793) | 0.389 (0.814)
0.75 | WMCEN | $\bm{1.285(0.112)}$ | $\bm{1.245(0.095)}$ | $\bm{1.224(0.086)}$ | $\bm{0.028(0.006)}$ | $\bm{0.026(0.005)}$ | $\bm{0.026(0.005)}$
| MCEN | 9.525 (4.610) | 9.520 (4.596) | 9.528 (4.585) | 3.145 (2.011) | 3.143 (2.011) | 3.143 (2.011)
| SEN | 11.895 (6.747) | 11.519 (6.682) | 11.421 (6.522) | 6.049 (4.485) | 5.802 (4.364) | 5.802 (4.175)
| WLASSO | 2.434 (1.405) | 2.433 (1.405) | 2.398 (1.203) | 0.423 (0.774) | 0.422 (0.773) | 0.422 (0.583)
1.00 | WMCEN | $\bm{1.596(0.323)}$ | $\bm{1.609(0.300)}$ | $\bm{1.606(0.305)}$ | $\bm{0.055(0.029)}$ | $\bm{0.057(0.028)}$ | $\bm{0.057(0.028)}$
| MCEN | 10.07 (4.807) | 10.069 (4.805) | 10.066 (4.809) | 3.443 (2.236) | 3.442 (2.236) | 3.442 (2.237)
| SEN | 11.785 (6.666) | 11.590 (6.533) | 11.962 (6.800) | 5.933 (4.358) | 5.698 (4.193) | 6.154 (4.564)
| WLASSO | 2.595 (1.405) | 2.598(1.414) | 2.600 (1.420) | 0.469 (0.808) | 0.470 (0.814) | 0.472 (0.820)

Table 8: Result of the simulation $p=100$ for Error $3$

|  | the median of APE (sd of APE) | MSE of ${\bm{\beta}}_{s}$ (sd of MSE)
---|---|---|---
| $\xi$ | 0.02 | 0.05 | 0.10 | 0.02 | 0.05 | 0.10
$\eta$ |  |  |  |  |  |  |
0.25 | WMCEN | 1.293 (0.025) | $\bm{1.205(0.024)}$ | $\bm{1.240(0.022)}$ | 0.011 (0.001) | $\bm{0.006(0.001)}$ | 0.009 (0.001)
| MCEN | $\bm{1.256(0.033)}$ | 1.247 (0.033) | 1.246 (0.036) | $\bm{0.008(0.001)}$ | 0.007 (0.001) | $\bm{0.007(0.001)}$
| SEN | 1.280 (0.038) | 1.262 (0.036) | 1.271 (0.036) | 0.010 (0.002) | 0.010 (0.003) | 0.010 (0.003)
| WLASSO | 1.956 (0.090) | 1.961 (0.093) | 1.941 (0.086) | 0.064 (0.013) | 0.064 (0.013) | 0.062 (0.012)
0.50 | WMCEN | 1.420 (0.036) | 1.438 (0.039) | $\bm{1.405(0.035)}$ | 0.020 (0.002) | 0.022 (0.002) | 0.019 (0.002)
| MCEN | $\bm{1.382(0.062)}$ | $\bm{1.396(0.061)}$ | 1.422 (0.060) | $\bm{0.016(0.003)}$ | $\bm{0.017(0.003)}$ | $\bm{0.018(0.003)}$
| SEN | 1.485 (0.050) | 1.466 (0.046) | 1.469 (0.060) | 0.024 (0.004) | 0.022 (0.003) | 0.024 (0.005)
| WLASSO | 2.049 (0.087) | 2.074 (0.092) | 2.034 (0.082) | 0.075 (0.013) | 0.077 (0.013) | 0.073 (0.012)
0.75 | WMCEN | 1.584 (0.053) | 1.588 (0.058) | 1.588 (0.062) | 0.033 (0.003) | 0.034 (0.004) | 0.033 (0.004)
| MCEN | $\bm{1.526(0.097)}$ | $\bm{1.509(0.098)}$ | $\bm{1.513(0.097)}$ | $\bm{0.025(0.005)}$ | $\bm{0.024(0.005)}$ | $\bm{0.024(0.005)}$
| SEN | 1.631 (0.066) | 1.642 (0.070) | 1.664 (0.065) | 0.037 (0.006) | 0.039 (0.006) | 0.039 (0.005)
| WLASSO | 2.195 (0.098) | 2.176 (0.095) | 2.199 (0.105) | 0.093 (0.015) | 0.091 (0.015) | 0.095 (0.017)
1.00 | WMCEN | 1.754 (0.085) | 1.754 (0.086) | 1.752 (0.082) | 0.049 (0.007) | 0.049 (0.007) | 0.049 (0.007)
| MCEN | $\bm{1.638(0.134)}$ | $\bm{1.638(0.134)}$ | $\bm{1.656(0.132)}$ | $\bm{0.033(0.007)}$ | $\bm{0.033(0.007)}$ | $\bm{0.034(0.007)}$
| SEN | 1.761 (0.083) | 1.744 (0.084) | 1.808 (0.083) | 0.050 (0.007) | 0.049 (0.008) | 0.054 (0.008)
| WLASSO | 2.288 (0.120) | 2.273 (0.115) | 2.288 (0.019) | 0.108 (0.019) | 0.105 (0.017) | 0.108 (0.019)

Table 9: Result of the simulation $p=100$ for Error $4$

|  | the median of APE (sd of APE) | MSE of ${\bm{\beta}}_{s}$ (sd of MSE)
---|---|---|---
| $\xi$ | 0.02 | 0.05 | 0.10 | 0.02 | 0.05 | 0.10
$\eta$ |  |  |  |  |  |  |
0.25 | WMCEN | $\bm{1.377(0.043)}$ | $\bm{1.378(0.043)}$ | $\bm{1.401(0.045)}$ | $\bm{0.011(0.001)}$ | $\bm{0.011(0.001)}$ | $\bm{0.012(0.002)}$
| MCEN | 2.473 (0.542) | 2.473 (0.540) | 2.472 (0.535) | 1.479 (9.455) | 1.512 (9.630) | 1.461 (9.968)
| SEN | 2.393 (0.643) | 2.638 (0.789) | 2.425 (0.669) | 582.884 (5031.615) | 670.449 (5871.294) | 535.364 (4371.789)
| WLASSO | 3.202 (0.401) | 3.202 (0.400) | 3.206 (0.399) | 0.218 (0.144) | 0.218 (0.144) | 0.219 (0.144)
0.50 | WMCEN | $\bm{1.654(0.057)}$ | $\bm{1.640(0.058)}$ | $\bm{1.639(0.062)}$ | $\bm{0.026(0.003)}$ | $\bm{0.025(0.003)}$ | $\bm{0.024(0.003)}$
| MCEN | 3.514 (0.890) | 2.895 (0.545) | 2.901 (0.549) | 39.493 (243.653) | 1.099 (6.496) | 0.833 (5.079)
| SEN | 3.247 (0.812) | 3.355 (0.893) | 3.029 (0.726) | 730.828 (6363.460) | 608.572 (5044.895) | 713.699 (6253.216)
| WLASSO | 3.507 (0.422) | 3.383 (0.397) | 3.382 (0.394) | 0.303 (0.173) | 0.248 (0.143) | 0.249 (0.143)
0.75 | WMCEN | $\bm{1.916(0.111)}$ | $\bm{1.971(0.123)}$ | $\bm{1.930(0.110)}$ | $\bm{0.041(0.005)}$ | $\bm{0.045(0.005)}$ | $\bm{0.042(0.004)}$
| MCEN | 3.565 (0.798) | 3.433 (0.722) | 3.556 (0.791) | 25.040 (139.366) | 18.090 (97.046) | 24.118 (133.008)
| SEN | 3.776 (0.838) | 3.756 (0.894) | 3.769 (0.889) | 666.696 (5723.088) | 631.024 (5330.103) | 723.491 (6093.913)
| WLASSO | 3.978 (0.484) | 3.979 (0.483) | 3.980 (0.484) | 61.319 (510.994) | 61.317 (510.980) | 61.315 (510.968)
1.00 | WMCEN | $\bm{2.217(0.104)}$ | $\bm{2.261(0.207)}$ | $\bm{2.253(0.202)}$ | $\bm{0.070(0.008)}$ | $\bm{0.066(0.009)}$ | $\bm{0.067(0.009)}$
| MCEN | 3.807 (0.754) | 3.789 (0.741) | 3.777 (0.763) | 23.917 (123.683) | 22.480 (114.956) | 20.724 (105.953)
| SEN | 4.413 (0.924) | 4.428 (1.017) | 4.351 (0.947) | 612.964 (5383.766) | 646.733 (5487.569) | 688.426 (6132.423)
| WLASSO | 4.248 (0.506) | 4.169 (0.482) | 4.063 (0.463) | 61.386 (511.104) | 0.686 (0.931) | 0.405 (0.175)

##  4 Real Example

In this section, we applied the proposed method to a genetic dataset named
"chin07" from package lol in R software (Chin et al., 2007) to verify the
method’s usefulness. The dataset consists of copy number patterns and mRNA
expression levels in genomic regions with candidate oncogenes for breast
cancer. In recent years, with the development of computational technologies,
there has been a lot of research on identifying genomic regions with candidate
oncogenes through genome-wide profiling, with the goal of suppressing the
expression of cancer. Altered DNA copy numbers are considered to be one of the
causes of genetic abnormalities leading to disease. Identifying which DNA copy
number alternation patterns affect the mRNA expression level from this dataset
is key. Therefore, it is necessary to improve the accuracy of estimation along
with variable selection. This dataset consists of a matrix of seven types of
mRNA expressions highly relevant to breast cancer with 106 samples and a
matrix of DNA copy number data for 339 regions with 106 samples. We set the
mRNA expression levels as response variables and the DNA copy numbers as
explanatory variables. Figure 1 draws a scatter plot, histogram, and
correlations of the response variables. GI$\\_$17318566-A, Hs.500472-S, and
GI$\\_$4758297-S appear to have right-tailed distribution, while
GI$\\_$38505204-S has a left-tailed distribution. GI$\\_$17318560-A contains
outliers. We compare the proposed method with MCEN, SEN, and WLASSO just as we
did with the numerical simulation.

Figure 1 draws a scatter plot, histogram, and correlations of the response
variables. GI$\\_$17318566-A, Hs.500472-S, and GI$\\_$4758297-S appear to have
right-tailed distribution, while GI$\\_$38505204-S has a left-tailed
distribution. GI$\\_$17318560-A contains outliers. We compare the proposed
method with MCEN, SEN, and WLASSO just as we did with the numerical
simulation.

Figure 1: Scatter plot, histogram, and correlations of the response variables

Next, we describe the evaluation procedures. The subjects are randomly split
into 80 samples for the training data and 26 samples for the test data. Each
method is then used to estimate ${\bm{\beta}}_{s}(s=1,2,\cdots,7)$ with the
training data, and an evaluation index is calculated by applying the estimated
${\bm{\beta}}_{s}(s=1,2,\cdots,7)$ to the test data. For the evaluation index,
the median APE is $|y_{is}-\hat{y}_{is}|(i=1,2,\cdots,26,s=1,2,\cdots,7)$,
where $\hat{y}_{is}$ is the predicted value, which is set in the same manner
as Wang and Li (2009). In this evaluation, both the response variables and
explanatory variables were standardized in the training and test data in the
same manner as Mukherjee and Zhu (2011). A five-fold cross-validation is
applied to the training data to decide the tuning parameter in all methods.
Candidate $k$ is set as 2 and 3 for the proposed method and MCEN. MCEN, SEN,
and WLASSO use the same R software packages as in the numerical simulation. We
compare the four methods by each response variable.

The results for the response variables are plotted in Figure 2. The vertical
axis describes the median APE. We compare the result of the proposed method to
each compared method. First, compared to MCEN, the proposed method had good
results in the response variables GI$\\_$17318566-A, GI$\\_$31543215-S,
Hs.500472-S, GI$\\_$38505204-S, and GI$\\_$16950654-S. The proposed method
also performed better than WLASSO for all response variables. In terms of the
comparison with SEN, the proposed method had smaller evaluation values in the
response variables GI$\\_$17318566-A, Hs.500472-S, GI$\\_$38505204-S, and
GI$\\_$17318560-S.

The results differed for each response variable, so we discuss the results
from the distribution of the response variables in Figure 1. The proposed
method was superior to the other methods in the response variables that have a
heavy tail such as GI$\\_$17318566-A, Hs.500472-S, GI$\\_$38505204-S, and
GI$\\_$4758297-S. We focus on the shapes of the distribution of the response
variables, for which the proposed method performed well. In the case of the
comparison to MCEN, the proposed method was better than MCEN, regardless of
the shape of the distribution of the response variables. Next, in the results
compared to SEN, the proposed method had better outcomes for the variables
GI$\\_$17318566-A, Hs.500472-S and GI$\\_$38505204-S, which follow heavy-
tailed distributions. In the variables GI$\\_$31543215-S, GI$\\_$4758297-S,
and GI$\\_$17318560-A, the proposed method was minimally inferior to SEN,
which might be due to the selection of the tuning parameters on each response
variable in SEN.

Figure 2: Results of the real data for each response variable

##  5 Discussion and Conclusion

In the numerical simulation, the proposed method was superior when the errors
followed heavy-tailed distributions or the data contained outliers, like Error
$2$, Error $3$, and Error $4$. Especially in Error $2$, which contained
outliers, and Error $4$, which had a strongly heavy-tailed setting, the
proposed method performed much better than the other compared methods in both
the median APE and MSE of ${\bm{\beta}}_{s}$. This indicates that the proposed
method is stable to heavy-tailed error distributions and outliers in the
error. From these results, our proposed method obtained the similar results as
Wang et al. (2020). In addition, among the compared methods, WLASSO had closer
results to the proposed method than MCEN and SEN in Error 2. This shows that a
Wilcoxon-type regression function is effective when there are outliers.
Moreover, our proposed method extended the Wilcoxon-type regression function
in the framework of multivariate regression, which enabled it to take into
account the correlation of response variables. Therefore, the results of the
proposed method were better than WLASSO. In Error $1$, where the error follows
normal distribution, the proposed method was superior to the other methods in
$p=12$. This result is consistent with the trend of results in the error with
normal distribution in Wang et al. (2020). However, SEN and MCEN are had
better results for $p=100$ due to the assumption of a normal distribution for
the residuals. When focusing on value differences, the values of the proposed
method were close to those of the compared methods in both $p=12$ and $p=100$.
Therefore, the proposed method can retain comparable estimation accuracy to
the compared methods in situations where the error distribution follows a
normal distribution, regardless of the number of variables.

In the real data application, the proposed method was better than the other
methods in the responses GI$\\_$17318566-A, Hs.500472-S, GI$\\_$38505204-S,
and GI$\\_$4758297-S, which followed heavy-tailed distribution. The proposed
method showed the stability of heavy-tailed situations in real data as well.
Meanwhile, in the variables GI$\\_$4758297-S and GI$\\_$17318560-A, the
proposed method was minimally inferior to SEN. The difference of the median
APE between the proposed method and SEN was not very broad. This is considered
to be due to the fact that SEN selects the tuning parameters for each response
variable, while the proposed method selected single tuning parameters for all
response variables in this study.

We showed the efficacy of the proposed method through the numerical simulation
and real data application. From these results, this method will be expected to
provide high estimation accuracy when applied to actual multivariate data
which includes heavy-tailed distributions and contains outliers. To improve
our proposed method, four points need to be considered. First, in this study,
we set common tuning parameters for all response variables. We expect that the
estimation accuracy of the method would enhance more if different tuning
parameters are set depending on each response variable. Secondly, we showed
that the Wilcoxon-type regression is more efficient with respect to the tailed
distribution or outliers of the response. However, Wilcoxon-type regression is
not very robust against the outliers in explanatory variables. To overcome
this, the weight of the Wilcoxon-type regression can be set, which has been
proposals in several studies (Naranjo and Hettmansperger, 1994; Sievers, 1983;
Wang and Li, 2009). This will need further consideration in situations where
the explanatory variables contain outliers. Third, our extension of the
Wilcoxon-type regression to multivariate regression is formulated differently
than the method proposed by Weihua (2010). Thus, we need to compare the two
methods from various perspectives. Finally, we extended the Wilcoxon-type
regression in the framework of multivariate regression. However, several
robust reduced-rank regressions for multivariate regression have also been
proposed (Chao et al., 2021; Ding et al., 2021; Wang and Karunamuni, 2022; Tan
et al., 2022; Mishra and MÃ¼ller, 2022; Zhao et al., 2017; She and Chen, 2017)
to deal with correlation of the response variables. The difference between our
proposed method and these robust reduced-rank methods is that the proposed
method has a clustering term to group the fitted values of the response
variables, which allows the correlation among the response variables to be
considered. Further consideration of this will be needed.

## References

  * Breiman and Friedman [1997] L. Breiman and J. H. Friedman.  Predicting multivariate responses in multiple linear regression.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 59(1):3–54, 1997.
  * Rothman et al. [2010] A. J. Rothman, E. Levina, and Z. Ji.  Sparse multivariate regression with covariance estimation.  _Journal of Computational and Graphical Statistics_ , 19(4):947–962, 2010.
  * Peng et al. [2010] J. Peng, J. Zhu, A. Bergamaschi, W. Han, D.Y. Noh, J.R. Pollack, and P. Wang.  Regularized multivariate regression for identifying master predictors with application to integrative genomics study of breast cancer.  _The annals of applied statistics_ , 4(1):53, 2010\.
  * Kim and Xing [2012] S. Kim and E.P. Xing.  Tree-guided group lasso for multi-response regression with structured sparsity, with an application to eqtl mapping.  _The Annals of Applied Statistics_ , 6(3):1095–1117, 2012.
  * Chen et al. [2017] Y. Chen, R. Iyengar, and G. Iyengar.  Modeling multimodal continuous heterogeneity in conjoint analysis—a sparse learning approach.  _Marketing Science_ , 36(1):140–156, 2017.
  * Cook et al. [2010] R. D. Cook, B. Li, and F. Chiaromonte.  Envelope models for parsimonious and efficient multivariate linear regression.  _Statistica Sinica_ , pages 927–960, 2010.
  * Cook and Zhang [2015] R.D. Cook and X. Zhang.  Foundations for envelope models and methods.  _Journal of the American Statistical Association_ , 110(510):599–611, 2015.
  * Sun et al. [2015] Q. Sun, H. Zhu, Y. Liu, and J.G. Ibrahim.  Sprem: sparse projection regression model for high-dimensional linear regression.  _Journal of the American Statistical Association_ , 110(509):289–302, 2015.
  * Price and Sherwood [2018] B. S. Price and B. Sherwood.  A cluster elastic net for multivariate regression.  _Journal of Machine Learning Research_ , pages 1–30, 2018.
  * Tibshirani [1996] R. Tibshirani.  Regression shrinkage and selection via the lasso.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 58(1):267–288, 1996.
  * Forgy [1965] E. W. Forgy.  Cluster analysis of multivariate data: efficiency versus interpretability of classifications.  _biometrics_ , 21:768–769, 1965.
  * Huber [2011] P.J. Huber.  Robust statistics.  In _International encyclopedia of statistical science_ , pages 1248–1251. Springer, 2011.
  * Fan et al. [2017] J. Fan, Q. Li, and Y. Wang.  Estimation of high dimensional mean regression in the absence of symmetry and light tail assumptions.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 79(1):247–265, 2017.
  * Loh [2017] P.L. Loh.  Statistical consistency and asymptotic normality for high-dimensional robust $m$-estimators.  _The Annals of Statistics_ , 45(2):866–896, 2017\.
  * Sun et al. [2019] Q. Sun, W.X. Zhou, and J. Fan.  Adaptive huber regression.  _Journal of the American Statistical Association_ , 2019.
  * Lozano et al. [2016] A.C. Lozano, N. Meinshausen, and E. Yang.  Minimum Distance Lasso for robust high-dimensional regression.  _Electronic Journal of Statistics_ , 10(1):1296 – 1340, 2016.
  * Wang et al. [2012a] L. Wang, Y. Wu, and R. Li.  Robust Variable Selection With Exponential Squared Loss.  _Journal of American Statistical Association_ , 108(1):632 – 643, 2012a.
  * Avella-Medina and Ronchetti [2018] M. Avella-Medina and E. Ronchetti.  Robust and consistent variable selection in high-dimensional generalized linear models.  _Biometrika_ , 105:31–44, 2018.
  * Prasad et al. [2020] A. Prasad, A.S. Sai, S. Balakrishnan, and P. Ravikumar.  Robust estimation via robust gradient estimation.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 82, 2020.
  * Belloni et al. [2011] A. Belloni, V. Chernozhukov, and L. Wang.  Square-root lasso: pivotal recovery of sparse signals via conic programming.  _Biometrika_ , 98(4):791–806, 2011.
  * Bradic et al. [2011] J. Bradic, J. Fan, and W. Wang.  Penalized composite quasi-likelihood for ultrahigh dimensional variable selection.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 73(3):325–349, 2011.
  * Wang et al. [2012b] L. Wang, Y. Wu, and R. Li.  Quantile regression for analyzing heterogeneity in ultra-high dimension.  _Journal of the American Statistical Association_ , 107(497):214–222, 2012b.
  * Wang [2013] L. Wang.  The l1 penalized lad estimator for high dimensional linear regression.  _Journal of Multivariate Analysis_ , 120:135–151, 2013\.
  * Fan et al. [2014] J. Fan, Y. Fan, and E. Barut.  Adaptive robust variable selection.  _Annals of statistics_ , 42 1:324–351, 2014.
  * Sun and Zhang [2012] T. Sun and C.H. Zhang.  Scaled sparse linear regression.  _Biometrika_ , 99(4):879–898, 2012.
  * Huber [1964] P.J. Huber.  Robust estimation of a location parameter.  _Annals of Mathematical Statistics_ , 35:492–518, 1964\.
  * Wang et al. [2007] H. Wang, G. Li, and G. Jiang.  Robust regression shrinkage and consistent variable selection through the lad-lasso.  _Journal of Business & Economic Statistics_, 25(3):347–355, 2007.
  * Zou and Yuan [2008] H. Zou and M. Yuan.  Composite quantile regression and the oracle model selection theory.  _The Annals of Statistics_ , 36(3):1108–1126, 2008.
  * Jureckova [1971] J. Jureckova.  Nonparametric estimate of regression coefficients.  _The Annals of Mathematical Statistics_ , pages 1328–1338, 1971.
  * Jaeckel [1972] L. A. Jaeckel.  Estimating regression coefficients by minimizing the dispersion of the residuals.  _The Annals of Mathematical Statistics_ , pages 1449–1458, 1972.
  * Hettmansperger and McKean [1977] T.P. Hettmansperger and J.W. McKean.  A robust alternative based on ranks to least squares in analyzing linear models.  _Technometrics_ , 19(3):275–284, 1977.
  * Hettmansperger and McKean [1998] T. P. Hettmansperger and J. W. McKean.  _Robust nonparametric statistical methods_.  Arnold, 1998.
  * Wang et al. [2020] L. Wang, B. Peng, J. Bradic, R. Li, and Y. Wu.  A tuning-free robust and efficient approach to high-dimensional regression.  _Journal of the American Statistical Association_ , 115(532):1700–1714, 2020.
  * Wang and Li [2009] L. Wang and R. Li.  Weighted wilcoxon-type smoothly clipped absolute deviation method.  _Biometrics_ , 65(2):564–571, 2009.
  * Hunter and Lange [2004] D. R. Hunter and K. Lange.  A tutorial on mm algorithms.  _The American Statistician_ , 58(1):30–37, 2004\.
  * Hettmansperger and McKean [1978] T. P. Hettmansperger and J. W. McKean.  Statistical inference based on ranks.  _Psychometrika_ , 43:69–79, 1978.
  * Hunter and Li [2005] D. R. Hunter and R. Li.  Variable selection using mm algorithm.  _Annals of statistics_ , 33(4):1617, 2005.
  * Yu et al. [2015] D. Yu, J. H. Won, T. Lee, J. Lim, and S. Yoon.  High-dimensional fused lasso regression using majorization-minimization and parallel processing.  _Journal of Computational and Graphical Statistics_ , 24(1):121–153, 2015.
  * Young [1981] F.W. Young.  Quantitative analysis of qualitative data.  _Psychometrika_ , 46:357–388, 1981.
  * Zou and Hastie [2005] H. Zou and T. Hastie.  Regularization and variable selection via the elastic net.  _Journal of the Royal Statistical Society Series B (statistical methodology)_ , 67(2):301–320, 2005.
  * Sherwood and Price [2020] B. Sherwood and B. Price.  _mcen: Multivariate Cluster Elastic Net_ , 2020.  URL https://CRAN.R-project.org/package=mcen.  R package version 1.2.
  * Friedman et al. [2008] J. Friedman, T. Hastie, and R. Tibshirani.  Regularization paths for generalized linear models via coordinate descent.  _Journal of statistical software_ , 33(1):1, 2008\.
  * Li and Qin [2021] S. Li and Y. Qin.  _MTE: Maximum Tangent Likelihood Estimation for Linear Regression_ , 2021.  URL https://CRAN.R-project.org/package=MTE.  R package version 1.0.1.
  * Chin et al. [2007] S. F. Chin, A. E. Teschendorff, J. C. Marioni, Y. Wang, N. L. Barbosa-Morais, N. P. Thorne, J. L. Costa, S. E. Pinder, M. A. Van de Wiel, A. R Green, et al.  High-resolution acgh and expression profiling identifies a novel genomic subtype of er negative breast cancer.  _Genome biology_ , 8(10):1–17, 2007.
  * Mukherjee and Zhu [2011] A. Mukherjee and J. Zhu.  Reduced rank ridge regression and its kernel extensions.  _Statistical analysis and data mining: the ASA data science journal_ , 4(6):612–622, 2011.
  * Naranjo and Hettmansperger [1994] J.D. Naranjo and T.P. Hettmansperger.  Bounded influence rank regression.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 56(1):209–220, 1994.
  * Sievers [1983] G.L. Sievers.  A weighted dispersion function for estimation in linear models.  _Communications in Statistics-Theory and Methods_ , 12(10):1161–1179, 1983.
  * Weihua [2010] Z. Weihua.  A multivariate wilcoxon regression estimate.  _Journal of Nonparametric Statistics_ , 22(7):859–877, 2010.
  * Chao et al. [2021] S. K. Chao, W. K. HÃ¤rdle, and M. Yuan.  Factorisable multitask quantile regression.  _Econometric Theory_ , 37(4):794–816, 2021.
  * Ding et al. [2021] H. Ding, S. Qin, Y. Wu, and Y. Wu.  Asymptotic properties on high-dimensional multivariate regression m-estimation.  _Journal of Multivariate Analysis_ , 183:104730, 2021.
  * Wang and Karunamuni [2022] Y. Wang and R. J. Karunamuni.  High-dimensional robust regression with lq-loss functions.  _Computational Statistics and Data Analysis_ , page 107567, 2022.
  * Tan et al. [2022] K. M. Tan, Q. Sun, and D. Witten.  Sparse reduced rank huber regression in high dimensions.  _Journal of the American Statistical Association_ , pages 1–11, 2022\.
  * Mishra and MÃ¼ller [2022] A. Mishra and C. L. MÃ¼ller.  Robust regression with compositional covariates.  _Computational Statistics and Data Analysis_ , 165:107315, 2022.
  * Zhao et al. [2017] W. Zhao, H. Lian, and S. Ma.  Robust reduced-rank modeling via rank regression.  _Journal of Statistical Planning and Inference_ , 180:1–12, 2017.
  * She and Chen [2017] Y. She and K. Chen.  Robust reduced-rank regression.  _Biometrika_ , 104(3):633–647, 2017.