%\VignetteIndexEntry{clere}
\documentclass[a4paper]{report}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{RJournal}
\usepackage{hyperref}
\usepackage{url}
\usepackage{amsmath,amssymb,array}
\usepackage[english]{babel}
\usepackage{bm}
\usepackage{amssymb,amsmath}
\usepackage{booktabs}


\def\1{\mbox{1\hspace{-.22em}I}}

\newcommand{\X}{\mathbb{X}}
\newcommand{\Rd}{\mathbb{R}^d}
\newcommand{\Xd}{\mathbb{X}^d}

\newcommand{\argmax}{\mathop{\mathrm{argmax}}}
\newcommand{\trace}{\mathop{\mathrm{Tr}}}
\newcommand{\diag}{\mathop{\mathrm{diag}}}
\newcommand{\II}{1 \! \! 1}
\newcommand{\IR}{\mathbb{R}}
\newcommand{\IZ}{\mathbb{Z}}
\newcommand{\IN}{\mathbb{N}}
\newcommand{\IE}{\mathbb{E}}
\newcommand{\mat}[4]{\begin{array}{cc}#1 & #2 \\#3 & #4 \end{array}}
\newcommand{\matb}[4]{\begin{array}{cc}{\bf #1} & {\bf #2} \\{\bf #3} & {\bf #4} \end{array}}
\newcommand{\med}{\mathrm{med}}
\newcommand{\tr}{\mbox{trace}}
\newcommand{\tra}[1]{\mbox{tr}{\bf #1}}
% \newcommand{\var}{\mbox{var}}
\newcommand{\ba}{\mathbf{a}}
\newcommand{\bc}{\mathbf{c}}
\newcommand{\bg}{\mathbf{g}}
\newcommand{\be}{\mathbf{e}}
\newcommand{\bp}{\mathbf{p}}
\newcommand{\bu}{\mathbf{u}}
\newcommand{\bx}{\mathbf{x}}
\newcommand{\bX}{\mathbf{X}}
\newcommand{\bZ}{\mathbf{Z}}
\newcommand{\bz}{\mathbf{z}}
\newcommand{\bt}{\mathbf{t}}
\newcommand{\by}{\mathbf{y}}
\newcommand{\x}{x}
\newcommand{\mx}{x}
\newcommand{\muu}{u}
\newcommand{\mv}{v}
\newcommand{\my}{y}
\newcommand{\mz}{z}
\newcommand{\z}{z}
\newcommand{\y}{y}
\newcommand{\Y}{Y}
\newcommand{\Z}{Z}
\newcommand{\balpha}{\boldsymbol{\alpha}}
\newcommand{\bDelta}{\boldsymbol{\Delta}}
\newcommand{\bepsilon}{\boldsymbol{\epsilon}}
\newcommand{\bGamma}{\boldsymbol{\Gamma}}
\newcommand{\blambda}{\boldsymbol{\lambda}}
\newcommand{\bmu}{\boldsymbol{\mu}}
\newcommand{\bpi}{\boldsymbol{\pi}}
\newcommand{\bphi}{\boldsymbol{\phi}}
\newcommand{\brho}{\boldsymbol{\rho}}
\newcommand{\btheta}{\boldsymbol{\theta}}
\newcommand{\bTheta}{\boldsymbol{\Theta}}
\newcommand{\bvarepsilon}{\boldsymbol{\varepsilon}}
\def\1{\mbox{1\hspace{-.22em}I}}
\newtheorem{defn}{Definition}[section]
\newtheorem{thm}{Theorem}[section]
\newtheorem{proposition}{Proposition}[section]
\newtheorem{corollary}{Corollary}[section]
\newtheorem{remark}{Remark}[section]
\newtheorem{lemma}{Lemma}[section]
\newcommand{\demi}{\mbox{$\frac{1}{2}$}}


\newcommand{\ra}[1]{\renewcommand{\arraystretch}{#1}}
\newcommand{\argmin}{\operatornamewithlimits{argmin}}

%% load any required packages here

\begin{document}

%% do not edit, for illustration only
\sectionhead{Contributed research article}
\volume{XX}
\volnumber{YY}
\year{2016}
\month{MARCH}

%% replace RJtemplate with your article
\begin{article}
\title{Variable Clustering in
        High-Dimensional Linear Regression: The \textbf{R} Package \textbf{clere}}
\author{by Loïc Yengo, Julien Jacques, Christophe Biernacki and Mickael Canouil}

\maketitle

\abstract{
Dimension reduction is one of the biggest challenge in high-dimensional regression models. We recently introduced
a new methodology based on variable clustering as a means to reduce dimensionality. We present here an R package that
implements this methodology. An overview of the package
functionalities as well as examples to run an analysis are
described. Numerical experiments on real data were performed to
illustrate the good predictive performance of our parsimonious method compared to standard dimension reduction approaches.}

\section[Introduction]{Introduction}

High dimensionality is increasingly ubiquitous in numerous scientific
fields including genetics, economics and physics. Reducing the
dimensionality is a challenge that most statistical methodologies must
meet not only to remain interpretable but also to achieve reliable
predictions. In linear regression models, dimension reduction
techniques often refer to variable selection. Approaches for variable
selection are implemented in publicly available software, that involve
the well-known \textbf{R} packages \textbf{glmnet} [\cite{glmnet-ref}]
and \textbf{spikeslab} [\cite{spikeslab-ref}]. The \textbf{R} package
\textbf{glmnet} implements the Elastic net methodology
[\cite{zou2005}], which is a generalization of both the LASSO
[\cite{tibshirani1996}] and the ridge regression (RR)
[\cite{hoerl1970}]. The \textbf{R} package \textbf{spikeslab} in turn,
implements the Spike and Slab methodology [\cite{ishwaran2005}], which
is a Bayesian approach for variable selection. 

Dimension reduction can not however be restricted to variable
selection. Indeed, the field can be extended to include approaches
which aim is to create surrogate covariates that summarize the
information carried in initial covariates. Since the emblematic
Principal Component Regression (PCR) [\cite{joliffe82}], many of the
other methods spread in the recent literature. As specific examples,
we may refer to the OSCAR methodology [\cite{bondell2008}], or the
PACS methodology [\cite{pacs2013}] which is a generalization of the
latter approach. Those methods mainly proposed variables clustering
within a regression model as a way to reduce the
dimensionality. Despite their theoretical and practical appeal,
implementations of those methods were often proposed only through
\textbf{Matlab} or \textbf{R} scripts, limiting thus the flexibility
and the computational efficiency of their use. The CLusterwise Effect
REgression (CLERE) methodology [\cite{yengo2013}], was recently
introduced as a novel methodology for simultaneous variables
clustering and regression. The CLERE methodology is based on the
assumption that each regression coefficient is an unobserved random
variable sampled from a mixture of Gaussian distributions with an
arbitrary number $g$ of components. In addition, all components in the
mixture are assumed to have different means ($b_1,\ldots,b_g$) and
equal variances $\gamma^2$. 

In this paper, we propose two new features for the CLERE model.
First, the stochastic EM (SEM) algorithm is proposed as a more
computationally efficicient alternative to the Monte Carlo EM (MCEM)
algorithm previously introduced in [\cite{yengo2013}].  Secondly,
the CLERE model is enhanced with the possibility of constraining the
first component to have its mean equal to 0, i.e. $b1 = 0$.  This
enhancement mainly aimed at facilitating the interpretation of the
model.  Indeed when $b_1$ is set to $0$,  variables assigned to the cluster
associated with $b_1$ might  be  considered  less  relevant  than  other
variables provided $\gamma^2$ to be small enough. Those  two  new
features  were implemented in the \textbf{R} package
\textbf{clere}. The core of the package is a \textbf{C++} program
interfaced with \textbf{R} using \textbf{R} packages \textbf{Rcpp}
[\cite{Rcpp-ref}] and \textbf{RcppEigen} [\cite{RcppEigen-ref}]. The
\textbf{R} package \textbf{clere} can be downloaded from the
Comprehensive \textbf{R} Archive Network (CRAN) at \url{http://cran.r-project.org/web/packages/clere/}. 

The outline of the present paper is the following. In the following
section the definition of the model is recalled and the strategy to
estimate the model parameter is presented. Afterwards, the main
functionalities of the \textbf{R} package \textbf{clere} are
presented. Real data analyses are then presented, aiming at
illustrating the good predictive performances of CLERE, with
noticeable parsimony ability, compared to
standard dimension reduction methods. Finally, perspectives and
further potential improvements of the package are discussed in the
last section. 

\section[Model definition and notation]{Model definition and notation}
\label{sec:modeldescription}
Our model is defined by the following hierarchical relationships: 
  \begin{equation}
    \label{eq:clere0}
    \left\{
    \begin{array}{l}
      y_i \sim\mathcal{N}\left(\beta_0+\sum_{j=1}^{p}{\beta_j x_{ij}},\sigma^2\right)\\
      \beta_j |\mathbf{z}_j  \sim \mathcal{N}\left(\sum^{g}_{k=1}{ b_k z_{jk} },\gamma^2  \right)\\
      \mathbf{z}_j = \left(z_{j1},\ldots,z_{jg} \right)\sim \mathcal{M}\left(1,\pi_1,\ldots,\pi_g\right)\text{,}
    \end{array}
    \right.
  \end{equation}
  {where $\mathcal{N}(\mu,\sigma^2)$ is the normal
  distribution of center $\mu$ and variance $\sigma^2$, and
  $\mathcal{M}\left(1,\pi_1,\ldots,\pi_g\right)$ the one-order
  multinomial distribution of parameters $\bm
  \pi=\left(\pi_1,\ldots,\pi_g \right)$ such as, $\forall$ $
  k=1,\ldots,g$ $\pi_k>0$ and $\sum^g_{k=1}{\pi_k}=1$, and $\beta_0$ is a constant term.}
  For an individual $i=1,\ldots,n$, $y_i$ is the response and $x_{ij}$ is an observed value for the $j$-th covariate. $\beta_j$ is the regression
  coefficient associated with the $j$-th covariate ($j=1,\ldots,p$),  {which is assumed to follow a mixture of $g$ Gaussians. The variable $\mathbf{z}_j$ indicates from which mixture component $\beta_j$ is drawn ($z_{jk}=1$ if $\beta_j$ comes from component $k$ of the mixture, $z_{jk}=0$ otherwise).}
  Let's note that model (\ref{eq:clere0}) can be considered as a variable selection-like model by constraining the model parameter $b_1$ to be equal to 0. 
  {Indeed, assuming that one of the components is centered in zero means that a cluster of regression coefficients have null expectation, and thus that the corresponding variables are not significant for explaining the response variable.}
  This functionality is available in the package.\\
  Let $\bm \beta  = \left(\beta_1,\ldots,\beta_p  \right)$, $\bm y = (y_1,\ldots ,y_n)'$, $\mathbf{X} = (x_{ij})$,
  $\bm Z = (z_{jk})$, $\bm b = (b_1\ldots b_g)'$ and $\bm \pi=
  (\pi_1,\ldots,  \pi_g)'$.
  Moreover, $\log p(\bm y|\bm X;\bm \theta)$ denotes the log-likelihood of model (\ref{eq:clere0}) assessed for
  the parameter  $\bm \theta = \left(\beta_0,\bm b,\bm\pi,\sigma^2,\gamma^2\right)$. 
  Model (\ref{eq:clere0}) can be interpreted as a Bayesian approach. However, to be fully Bayesian a prior
  distribution for parameter $\bm \theta$ would have been
  necessary. Instead, we proposed to estimate $\bm \theta$ by maximizing the (marginal)
  log-likelihood, $\log p(\mathbf{y}|\mathbf{X};\bm \theta)$. This partially Bayesian approach is referred to as \textit{Empirical
  Bayes} (EB) [\cite{casella1985}]. Let $\mathcal{Z}$ be the set of $p\times g$-matrices partitioning $p$ covariates into $g$ groups. Those matrices are defined as
  \[
  \bm Z = \left(z_{jk}\right)_{1\leq j \leq p, 1\leq k \leq g} \in
  \mathcal{Z} \Leftrightarrow \forall j=1,\ldots,p\text{ }
  \begin{cases}
    \exists!\text{ } k\text{ such as }z_{jk} = 1\\
    \text{For all }k'\neq k\text{ }z_{jk'}=0.
  \end{cases}
  \] 
  The log-likelihood  $\log p(\bm y|\bm X;\bm \theta)$ is defined as
\begin{equation*}   
  \label{eq:likelihood}
  \log
  p(\bm y|\bm X;\bm \theta)
  = \log\left[ \sum_{\bm Z\in\mathcal{Z}} {\int_{\mathbb{R}^p}{p(\bm
	y,\bm \beta,\bm Z|\bm X;\bm \theta)}\mathrm{d}\bm \beta }\right]\text{.}
  \end{equation*}
Since it requires integrating over $\mathcal{Z}$ with cardinality $g^p$, evaluating the likelihood becomes rapidly
computationally unaffordable.

Nonetheless, maximum likelihood estimation is still achievable using the
expectation maximization (EM) algorithm [\cite{dempster1977}]. The
latter algorithm is an iterative method which starts with an initial
estimate of the parameter and updates this estimate until
convergence. Each iteration of the algorithm consists of two steps,
denoted as the \textit{E} and the \textit{M} steps. At each iteration
$d$ of the algorithm, the \textit{E  step} consists in calculating the
expectation of the log-likelihood of the complete data (observed +
unobserved) with respect to $p(\bm \beta,\bm Z|\bm y,\bm X;\bm
\theta^{(d)})$, the conditional distribution of the
unobserved data 
%(the regression coefficients $\bm \beta$ and the matrix  of cluster membership indicators $\bm Z$) 
given the observed data, 
%(the vector of response $\bm y$)
and the value of the parameter at the current
iteration, $\bm \theta^{(d)}$. This expectation, often denoted as $Q(\bm \theta|\bm \theta^{(d)})$
is then maximized with respect to $\bm \theta$ at the \textit{M step}.  

In model (\ref{eq:clere0}), the \textit{E  step} is analytically
intractable. A broad literature devoted to intractable \textit{E  steps}
recommends the use of a stochastic approximation of $Q(\bm \theta|\bm \theta^{(d)})$ 
through Monte Carlo (MC) simulations [\cite{wei1990}, \cite{levine2001}]. This
approach is referred to as the MCEM algorithm. Besides, mean-field-type
approximations are also proposed [\cite{govaert2008}, \cite{mariadassou2010}]. Despite their
computational appeal, the latter approximations do not generally ensure
convergence to the maximum likelihood
[\cite{Gunawardana05convergencetheorems}]. Alternatively, the
SEM algorithm [\cite{celeux:inria-00074164}] was introduced as a
stochastic version of the EM algorithm. In this algorithm, the
\textit{E  step} is replaced with a simulation step (\textit{S  step}) that consists in
generating a complete sample by simulating the unobserved data using
$p(\bm \beta,\bm Z|\bm y,\bm X;\bm \theta^{(d)})$ providing thus a sample $(\bm \beta^{(d)},\bm Z^{(d)})$. 
Note that the Monte Carlo algorithm we use to perform this simulation
is the Gibbs sampler. After the
\textit{S  step} follows the \textit{M step} which consists in
maximizing $p(\bm \beta^{(d)},\bm Z^{(d)}|\bm y,\bm X;\bm \theta)$
with respect to $\bm \theta$. Alternating those two steps generate a
sequence $\left(\bm \theta^{(d)}\right)$, which is a Markov chain
whose stationary distribution (when it exists) concentrates around a
local maximum of the likelihood. 

%A comparison of both algorithms has been carried out in \cite{Yengo2014}. As long as $p>n$, it has been shown that the SEM algorithm has a lower complexity than the MCEM algorithm because a closed form of $p(\bm y,\bm Z|\bm X;\bm \theta)$ is available. However, the computational time to run our SEM algorithm is more variable compared to MCEM as its \text{M step} does no longer have a closed form: an internal EM algorithm is thus required. The use of the MCEM algorithm is thus advocate only when $p<n$. 
%Note finally that the estimation of the model for different number of clusters can easily be parallelized (this option is available in the package).


\section{ { Estimation and model selection}}
In this section, two algorithms for model inference are presented: the Monte-Carlo Expectation
Maximization (MCEM) algorithm and the Stochastic Expectation Maximization (SEM) algorithm.
The section starts with the initialization strategy common to both algorithms and continues with the
detailed description of each algorithm. Then, model selection (for choosing g) and variable selection
are discussed.
\label{sec:algos}
\subsection{ { Initialization} }
 { 
  The two algorithms presented in this section are initialized using a primary estimate ${\beta_j}^{(0)}$ of each
  $\beta_j$. The latter can be chosen either at random, or obtained from univariate regression
  coefficients or penalized approaches like LASSO and ridge regression. For large SEM or MCEM chains, initialization is not a critical issue. The choice of the initialization strategy is therefore made to speed up the convergence of the chains. A Gaussian mixture model with $g$ component(s) is then fitted using ${\bm \beta}^{(0)} = \left(\beta_1^{(0)},\ldots,{\beta}^{(0)}_p\right)$ as observed data to produce starting values $\mathbf{b}^{(0)}$, $\bm \pi^{(0)}$ and ${\gamma^2}^{(0)}$ respectively for parameters $\mathbf{b}$, $\boldsymbol{\pi}$ and
  ${\gamma^2}$. Using maximum a posteriori (MAP) clustering, an initial partition $\mathbf{Z}^{(0)} =  \left(z^{(0)}_{jk}\right)\in\mathcal{Z}$ is obtained as
  \[
  \forall j\in\{1,\ldots,p\},\text{ }z^{(0)}_{jk} = \begin{cases}
    1 & \text{ if }k = \argmin_{k'\in\{1,\ldots,g\}}{\left({\beta_j}^{(0)}-b^{(0)}_{k'}\right)^2}\\
    0 & \text{otherwise.}
  \end{cases}
  \]
  $\beta_0$ and $\sigma^2$ are initialized using ${\bm \beta}^{(0)}$ as follows:
  \small
  $$ \beta_0^{(0)}  =
  \frac{1}{n}\sum^{n}_{i=1}{\left(y_i-\sum^{p}_{j=1}{\beta^{(0)}_j x_{ij}}\right)}
  \text{ and }
  {\sigma^2}^{(0)}  =
  \frac{1}{n}\sum^{n}_{i=1}{\left(y_i-\beta_0^{(0)}
  -\sum^{p}_{j=1}{\beta^{(0)}_j x_{ij}}\right)^2} .$$
  \normalsize
}
\subsection{ { MCEM algorithm}}
\label{sec:mcem}
\subsubsection{ { The Stochastic Approximation of the E step}}
 { 
  Suppose at iteration $d$ of the algorithm that we have
  $\left\{\left(\bm \beta^{(1,d)},\bm Z^{(1,d)}\right),\ldots,
  \left(\bm \beta^{(M,d)},\bm Z^{(M,d)}\right)\right\}$, $M$ samples
  from $p\left(\bm \beta,\bm Z|\bm y,\bm X;\bm  \theta^{(d)}\right)$.
  Then the MC approximation of the \textit{E}-step can
  be written 
  \begin{equation*}
    \label{eq:approxEstep}
    Q\left(\bm \theta|\bm \theta^{(d)}\right) = \mathbb{E}\left[\log p(\mathbf{y},\bm \beta,\mathbf{Z}|\mathbf{X};\bm \theta^{(d)}) |\mathbf{y},\mathbf{X};\bm \theta^{(d)}\right] \approx \frac{1}{M}\sum^{M}_{m=1}{\log p(\mathbf{y},\bm \beta^{(m,d)},\mathbf{Z}^{(m,d)}|\mathbf{X};\bm \theta^{(d)})}\text{.}
  \end{equation*}  
  However, sampling from $p\left(\bm \beta,\bm Z|\bm y,\bm X;\bm
  \theta^{(d)}\right)$ is not straightforward. However, we can use a
  Gibbs sampling scheme to simulate unobserved data, taking advantage
  of $p\left(\bm \beta|\bm Z,\bm y,\bm X;\bm \theta^{(d)}\right)$ and
  $p\left(\bm Z|\bm \beta,\bm y,\bm X;\bm \theta^{(d)}\right)$ from
  which it is easy to simulate. Those
  distributions, respectively Gaussian and multinomial, are described below in Equations (\ref{eq:gibbs1}) and (\ref{eq:gibbs2}).
  \begin{equation}
    \label{eq:gibbs1}
    \left\{
    \begin{array}{l}
      \bm \beta|\bm{Z,y,X};\bm \theta^{(d)} \sim \mathcal{N}\left(\bm \mu^{(d)},\bm \Sigma^{(d)} \right)\\
      \bm \mu^{(d)} = \left[ \bm{X'X} +
      \frac{{\sigma^2}^{(d)}}{{\gamma^2}^{(d)}}\text{I}_{p} \right]^{-1}\bm{X}'\left(\mathbf{y}-\beta^{(d)}_0\bm 1_p\right) +
      \frac{{\sigma^2}^{(d)}}{{\gamma^2}^{(d)}}\left[ \bm{X'X} +  \frac{{\sigma^2}^{(d)}}{{\gamma^2}^{(d)}}\text{I}_{p} \right]^{-1}\bm{Zb}^{(d)}\\
      \bm \Sigma^{(d)} = {\sigma^2}^{(d)}\left[ \bm{X'X} + \frac{{\sigma^2}^{(d)}}{{\gamma^2}^{(d)}}\text{I}_{p} \right]^{-1}
    \end{array}
    \right.
  \end{equation}
  and (note that $p\left(\bm Z|\bm \beta,\bm y,\bm X;\bm \theta^{(d)}\right)$ does not depend on $\bm{X}$ nor $\bm{y}$)
  \begin{equation}
    \label{eq:gibbs2}
    p\left(z_{jk}=1|\bm \beta; \bm \theta^{(d)}\right) \propto \pi^{(d)}_k\exp\left(-\frac{\left(\beta_j - b^{(d)}_k\right)^2}{2{\gamma^2}^{(d)}} \right).
  \end{equation}
  In Equation (\ref{eq:gibbs1}), $\text{I}_{p}$ and $\bm 1_p$ respectively stands for the identity matrix in dimension $p$ and the vector of $\mathbb{R}^p$ which all coordinates equal 1. To efficiently sample from $p\left(\bm \beta|\bm Z,\bm y,\bm X;\bm \theta^{(d)}\right)$ a preliminary singular vector decomposition of matrix $\bm X$ is necessary. Once   this decomposition is performed the overall complexity of the approximated \textit{E step} is $\mathcal{O}\left[M(p^2+pg)\right]$.
}
  \subsubsection{ { The M step}}
   { Using the $M$ draws obtained by Gibbs sampling at iteration $d$, the
  \textit{M step} is straightforward as detailed in Equations (\ref{eq:mstep_pi}) to (\ref{eq:mstep_sigma}). The overall computational complexity of that
  step is  $\mathcal{O}\left( Mpg\right)$.
  \begin{equation}
    \label{eq:mstep_pi}
    \pi^{(d+1)}_k = \frac{1}{Mp}\sum^{M}_{m=1}{  \sum^{p}_{j=1}{z^{(m,d)}_{jk} } }\text{,}
  \end{equation}  
  \begin{equation}
    \label{eq:mstep_b}
    b^{(d+1)}_k = \frac{1}{Mp\pi^{(d+1)}_k}\sum^{M}_{m=1}{\sum^{p}_{j=1}{ z^{(m,d)}_{jk} \beta_{j}^{(m,d)} } }\text{,}
  \end{equation}
  \begin{equation}
    \label{eq:mstep_gamma}
    {\gamma^2}^{(d+1)} =
    \frac{1}{Mp}\sum^{M}_{m=1}{\sum^{p}_{j=1}{  \sum^{g}_{k=1}{z^{(m,d)}_{jk}    \left(\beta_{j}^{(m,d)}-b^{(d+1)}_k\right)^2} } }\text{,}
  \end{equation}
  \begin{equation}
    \label{eq:mstep_beta0}
    \beta^{(d+1)}_0 = \frac{1}{n}\sum^{n}_{i=1}{ \left[y_i-\sum^{p}_{j=1}{\left(\frac{1}{M}\sum^{M}_{m=1}{\beta^{(m,d)}_j} \right)x_{ij}} \right]}\text{,}
  \end{equation}  
  \begin{equation}
    \label{eq:mstep_sigma}
    {\sigma^2}^{(d+1)} = \frac{1}{nM}\sum^{M}_{m=1}{\sum^{n}_{i=1}{ \left(y_i - \beta^{(d+1)}_0-\sum^{p}_{j=1}{ \beta^{(m,d)}_j x_{ij} }\right)^2  } }\text{.}
  \end{equation}
}
\subsection{ { SEM algorithm}}
   { In most situations, the SEM algorithm can be considered as a special case
   of the MCEM algorithm [\cite{celeux:inria-00074164}], obtained by setting
  $M=1$. In model (\ref{eq:clere0}), such a direct derivation leads to an algorithm which computational complexity remains quadratic with respect to $p$. To reduce that complexity, we propose a SEM algorithm based on the integrated complete data likelihood $p(\bm y, \bm Z|\bm X;\bm \theta)$ rather than $p(\bm y,\bm \beta,\bm Z|\bm X;\bm \theta)$.  A closed form of $p(\bm y, \bm Z|\bm X;\bm \theta)$ is available and given subsequently.
} 

  \subsubsection{ {Closed form of the integrated complete data likelihood}}
       { Let the SVD decomposition of matrix $\bm X$ be $\bm {USV'}$, where
         $\bm U$ and $\bm V$ are respectively $n \times n$ and 
    $p \times p$ orthogonal matrices, and $\bm S$ is $n \times p$
    rectangular diagonal matrix which diagonal terms are the eigenvalues $\left(\lambda^2_1,\ldots,\lambda^2_n\right)$
    of matrix $\bm {XX}'$. We now define $\bm X^u = \bm {U'X}$
    and $\bm y^u = \bm {U'y}$. Let $\bm M$ be the $n\times (g+1)$ matrix
    which first column is made  of 1's and which additional columns are
    those  of matrix $\bm X^u \bm Z$.
    Let	also $\bm t=\left(\beta_0,\bm b\right)\in \mathbb{R}^{(g+1)}$ and $\bm R$ be a $n\times n$ diagonal matrix which $i$-th diagonal term equal
    $\sigma^2 + \gamma^2\lambda^2_i$. With these notations
    we can express the complete data likelihood integrated over $\bm
    \beta$ as 
    \begin{align} 
      \label{eq:ICDLL}
      \log p\left(\bm{y,Z}|\bm{X};\bm \theta\right)  &= -\frac{n}{2}\log\left( 2\pi\right)-\frac{1}{2}\sum^{n}_{i=1}{\log\left( \sigma^2 + \gamma^2\lambda^2_i\right)}-\frac{1}{2}\left(\bm{y}^u-\bm {Mt}\right)'\bm R^{-1}\left(\bm{y}^u-\bm {Mt}\right) \nonumber\\
       & + \sum^{p}_{j=1}{\sum^{g}_{k=1}{z_{jk}\log \pi_k}}\text{.}
    \end{align}}


 \subsubsection{ { Simulation step}}
     { To sample from $p\left(\bm{Z}|\bm{y,X};\bm \theta\right)$ we use a
    Gibbs sampling strategy based on the conditional distributions
    $p\left(\bm{z}_j|\bm y,\bm Z^{-j}, \bm X;\bm \theta\right)$, $\bm
    Z^{-j}$ denoting the set of cluster membership indicators for all
    covariates but the $j$-th. Let $\bm w^{-j} =
    \left(w^{-j}_1,\ldots,w^{-j}_n\right)'$, where $w_i^{-j} =
    y^u_i-\beta_0 - \sum_{l\neq j}{\sum^{g}_{k=1}{z_{lk} x^u_{il}
	b_k}}$. The conditional distribution  $p(z_{jk} =
    1|\bm{Z}^{-j},\bm{y},\bm{X};\bm \theta)$ can be written
    \begin{equation}
      \label{eq:Gibbs}
      p(z_{jk} = 1|\bm{Z}^{-j},\bm{y},\bm{X};\bm \theta) \propto
      \pi_k\exp\left[-\frac{b^2_k}{2} \left(\bm{x}^{u}_j\right)'\bm{R}^{-1} \bm{x}^{u}_j+b_k \left(\bm{w}^{-j}\right)'\bm{R}^{-1}\bm x^{u}_j\right]\text{,}
    \end{equation}
    where $\bm{x}^{u}_j$ is the $j$-th column of $\bm X^u$. In the classical
    SEM algorithm, convergence to $p\left(\bm{Z}|\bm{y,X};\bm
    \theta\right)$ should be reached before updating $\bm
    \theta$. However, a valid inference can still be ensured in settings when $\bm \theta$ is updated only after
    one or few Gibbs iterations. These approaches are referred to as SEM-Gibbs
    algorithm [\cite{biernacki2013}]. The overall computational complexity of the
    simulation step is $\mathcal{O}\left( npg\right)$, so linear with
     $p$ and no more quadratic contrarily to the previous MCEM.\\
    To improve the mixing of the generated Markov chain, we start the
    simulation step at each iteration by creating a random permutation of 
    $\left\{1,\ldots,p\right\}$. Then, according to the order defined by
    that permutation, we update each $z_{jk}$ using $p(z_{jk} =
    1|\bm{Z}^{-j},\bm{y},\bm{X};\bm \theta)$.
  \subsubsection{Maximization step}
    $\log p\left(\bm{y,Z}|\bm{X};\bm \theta\right)$ corresponds to the
    marginal log-likelihood of a linear mixed model
    [\cite{searle1992variance}] which can be written 
    \begin{equation}
      \label{mixedmodel0}
      \bm{y}^u = \bm M \bm t + \bm {\lambda v} + \bm \varepsilon
    \end{equation}
    where $\bm v$ is an unobserved random vector such as $\bm v \sim
    \mathcal{N}\left(0,\gamma^2\text{I}_{n}\right)$, $\bm \varepsilon \sim
    \mathcal{N}\left(0,\sigma^2\text{I}_{n}\right)$ and $\bm \lambda =
    \text{diag}\left(\lambda_1,\ldots,\lambda_n\right)$. The estimation of
    the parameters of model (\ref{mixedmodel0}) can be performed using the
    EM algorithm, as in [\cite{searle1992variance}]. We adapt below the EM equations defined
    in [\cite{searle1992variance}], using our notations.
    At iteration $s$ of the internal EM algorithm, we define $\bm R^{(s)} =
    {\sigma^2}^{(s)}\bm I_{n} + {\gamma^2}^{(s)}\bm \lambda'\bm
    \lambda$. The detailed \textit{internal E} and \textit{M steps} are
    given below: \\
    \underline{\textit{Internal E step}}:\\
    \small
    \begin{eqnarray*}
      v^{(s)}_\sigma & = & \mathbb{E}\left[\left(\bm{y}^u - \bm{Mt}^{(s)} - \bm {\lambda v}\right)'\left(\bm{y}^u - \bm{Mt}^{(s)} - \bm {\lambda v}\right)|\bm{y}^u\right] \\
      & = & {\sigma^4}^{(s)}\left(\bm{y}^u - \bm{Mt}^{(s)}\right)'\bm R^{(s)}\bm R^{(s)}\left(\bm{y}^u - \bm{Mt}^{(s)}\right)  + n\times {\sigma^2}^{(s)} - {\sigma^4}^{(s)} \sum^{n}_{i=1}{\frac{1}{{\sigma^2}^{(s)}+{\gamma^2}^{(s)}\lambda^2_i}}\text{.}\\
      v^{(s)}_\gamma & = & \mathbb{E}\left[\bm v'\bm v|\bm{y}^u\right] \\
      & = & {\gamma^4}^{(s)}\left(\bm{y}^u - \bm{Mt}^{(s)} \right)'\bm R^{(s)}\bm \lambda'\bm \lambda \bm R^{(s)}\left(\bm{y}^u - \bm{Mt}^{(s)}\right) + n \times{\gamma^2}^{(s)} - {\gamma^4}^{(s)} \sum^{n}_{i=1}{\frac{\lambda^2_i}{{\sigma^2}^{(s)}+{\gamma^2}^{(s)}\lambda^2_i}}\text{.}\\
      \bm h^{(s)} & = & \mathbb{E}\left[\bm{y}^u - \bm {\lambda v}  |\bm{y}^u\right]  = \bm{Mt}^{(s)} + {\sigma^2}^{(s)} \{R^{(s)}\}^{-1} \left(\bm{y}^u - \bm{Mt}^{(s)}\right)\text{.}
    \end{eqnarray*}
    \normalsize
    \underline{\textit{Internal M step}}:\\
    \begin{equation*}
      {\sigma^2}^{(s+1)} = v^{(s)}_\sigma / n\text{.}
    \end{equation*}
    \begin{equation*}
      {\gamma^2}^{(s+1)} = v^{(s)}_\gamma / n\text{.}
    \end{equation*}
    \begin{equation*}
      \bm{t}^{(s+1)} = \left[\bm M'\bm M\right]^{-1}\bm M'\bm h^{(s)}\text{.}
    \end{equation*}
    Given a non-negative user-specified threshold $\delta$ and a maximum number $N_{max}$ of iterations, \textit{Internal E} and \textit{M steps}
    are alternated until 
    \begin{equation*}
      |\log p\left(\bm{y,Z}|\bm{X};\bm \theta^{(s)}\right)- \log p\left(\bm{y,Z}|\bm{X};\bm \theta^{(s+1)}\right)|<\delta\text{ or } s = N_{max}\text{.}
    \end{equation*}
    The computational complexity of the
    \textit{M step} is $\mathcal{O}\left( g^3 + ngN_{max}\right)$, thus not involving $p$. 
\subsubsection{Attracting and absorbing states}
\begin{itemize}
\item \textit{Absorbing states}. The SEM algorithm described above defines a Markov chain which stationnary distribution is concentrated around values of $\bm \theta$ corresponding to local maxima of the likelihood function. This chain has absorbing states in values of $\bm \theta$ such as $\sigma^2=0$ or $\gamma^2=0$. In fact, the \textit{internal M step} reveals that updated values for $\sigma^2$ and $\gamma^2$ are proportional to previous values of those parameters. 
\item \textit{Attracting states}. We empirically observed that
     attraction around $\sigma^2=0$ was quite frequent when using the
     MCEM algorithm, especially when $p>n$ and when the number $M$ of
     draws was small. We therefore advocate to use at least 5 draws
     ($M \geq 5$ using option \code{nsamp=} in the function \code{fitClere}).
\end{itemize}
}
\subsection{ { Model selection}}
 { Once the MLE $\widehat{\bm \theta}$ is calculated (using one or the other
  algorithm), the maximum log-likelihood and the posterior clustering matrix $\mathbb{E}\left[\bm Z|\bm{y, X};\widehat{\bm \theta} \right]$ are  approximated using MC simulations based on Equations (\ref{eq:ICDLL}) and  (\ref{eq:Gibbs}). The approximated maximum log-likelihood $\widehat{l}$, is then utilized to calculate AIC [\cite{akaike1974}] and BIC [\cite{schwarz1978}] criteria for model selection. In model (\ref{eq:clere0}), those criteria can be written as
  \begin{equation}
    \text{BIC} = -2\widehat{l} + 2(g+1)\log (n)\text{ and }\text{AIC} = -2\widehat{l} + 4(g+1)\text{.}
  \end{equation}
 An additional criterion for model selection, namely the ICL criterion [\cite{biernacki2000}] is also implemented in the \textbf{R} package \textbf{clere}. The latter criterion can be written
  \begin{equation}
    \text{ICL} = \text{BIC} - \sum^{p}_{j=1}{\sum^{g}_{k=1}{\pi_{jk} \log ( \pi_{jk} ) }}\text{,}
  \end{equation}
  where $\pi_{jk} = \mathbb{E}\left[z_{jk}|\bm{y, X};\widehat{\bm \theta} \right]$.}

\subsection{ { Interpretation of the special group of variables associated with $b_1=0$}}
 { The constraint $b_1=0$ is mainly driven by an interpretation purpose. The meaning of this group depends on both the total number $g$ of groups and the estimated value of parameter $\gamma^2$. In fact, when $g>1$ and $\gamma^2$ is small, covariates assigned to that group are likely less relevant to explain the response. Determining whether $\gamma^2$ is small enough is not straightforward. However, when this property holds, we may expect the groups of covariates to be separated. This would for example translate in the posterior probabilities $\pi_{j1}$ being larger than 0.7.  In addition to the benefit in interpretation, the constraint $b_1=0$, reduces the number of parameters to be estimated and consequently the variance of the predictions performed using the model.}



\section[Package functionalities]{Package functionalities}
\label{sec:package}
The \textbf{R} package \textbf{clere} mainly implements a function for parameter
estimation and model selection: the function \code{fitClere()}. Four additional
functions for graphical representation \code{plot()}, summarizing the
results \code{summary()}, for getting the predicted clusters of variables \code{clusters()} and for making predictions from new design matrices \code{predict()} are
also implemented in the package. Examples of calls for the functions presented in this section are given in the next Section. 

\subsection[The function fitClere()]{The main function \code{fitClere()}}

{The main function \code{fitClere()} has only three mandatory arguments: the vector of response \code{y} (size $n$), the matrix of explanatory variable \code{x} (size $n\times p$) and the number $g$ of groups of regression coefficients which is expected. The optional parameter \code{analysis}, when it takes the value  \code{aic},  \code{bic} or \code{icl}, allows to test all the possible number of groups between $1$ and $g$. The choice between the two proposed algorithms is possible thanks to the parameter \code{algorithm}, but we encourage the users to use the default value, the SEM algorithm, which generally overperforms the MCEM algorithm (see the first experiment of the next section).\\
Several other parameters allow to tune the different numbers of iterations of the estimation algorithm. Generally, higher are these parameters values, better is the quality of the estimation but heavier is the computing time. What we advice is to use the default values, and to graphically check the quality of the estimation with plots as in Figure \ref{fig:plot}: if the values of the model parameters are quite stable for a sufficient large part of the iterations, it is ok. If the stability is not reached sufficiently early before the end of the iterations, higher numbers of iterations should be chosen.\\
Finally, among the remaining parameters (the complete list can be obtained with \code{help(fitClere)}), two are especially important: \code{parallel} allows to run parallel computations (if compatible with the user's computer) and \code{sparse} allows to impose that one of the regression parameter is equal to 0 and thus to introduce a cluster of not significant explanatory variables.}

\subsection[Companion functions summary(), plot(), clusters() and predict()]{Secondary methods \code{summary()}, \code{plot()}, \code{clusters()} and \code{predict()}}

The \code{summary()} function prints an overview of the estimated parameters and returns the estimated likelihood and information based model selection criteria (AIC, BIC and ICL). 

The call of function \code{plot()} is similar to the one of function
 \code{summary()}. The latter function produces graphs such as ones
 presented in Figure \ref{fig:plot}. 

%\begin{figure}[!ht]    
%    \centerline{
%      \includegraphics[scale=0.8]{images/Figure1.pdf}
%    }
%    \caption{\label{fig:plot} Values of the model parameters in view of SEM algorithm iterations. The vertical grey line in each of the four plots, represents the number %\code{nBurn} of iterations discarded before calculating maximum likelihood estimates.}
%\end{figure}  

The function \code{clusters()}, takes one argument of class \code{Clere} and a \code{threshold} argument. This function assigns each variable to the group which associated conditional probability of membership is larger than the given \code{threshold}. If conditional probabilities of membership are larger than the specified threshold for more than one group, then the group having the largest probability is returned and warning is printed. If moreover, there are ex-aequos on that largest probability then the group with the smallest index is returned. When \code{threshold = NULL}, the maximum a posteriori (MAP) strategy is used to infer the clusters.

The \code{predict()} function has two arguments, being a \code{clere} and a design matrix $\bm X_{new}$. Using that new design matrix, the \code{predict()} function returns an approximation of $\mathbb{E}\left[\bm X_{new}\bm \beta|\bm y, \bm X; \hat{\bm \theta}\right]$.

\section[Numerical experiments on real datasets]{Numerical experiments}
\label{sec:numerical}
{This section presents two sets of numerical experiments. The first set of experiments aims at comparing the MCEM and SEM algorithms in terms of computational time and estimation or prediction accuracy. The second set of experiments aimed at comparing CLERE to standard dimension reduction techniques. The latter comparison is performed on both simulated and real data.}
\subsection{SEM algorithm versus MCEM algorithm}
\label{sec:semvsmcem}
\subsubsection{Description of the simulation study}
{
In this section, a comparison between the SEM algorithm and the MCEM algorithm is performed.
This comparison is performed using the four following performance indicators:
\begin{enumerate}
\item Computational time (CT) to run a pre-defined number of SEM/MCEM iterations. This number was set to 2,000 in this simulation study.
\item Mean squared estimation error (MSEE) defined as 
\begin{equation*} 
 MSEE_a = \mathbb{E}\left[(\bm{\theta-\widehat{\theta}_a})'(\bm{\theta-\widehat{\theta}_a})\right]\text{,}
\end{equation*}
where $a\in\left\{\text{\code{"SEM","MCEM"}}\right\}$ and $\widehat{\theta}_a$ is an estimated value for parameter $\bm \theta$ obtained with algorithm $a$. Since $\bm \theta$ is only known up to a permutation of the group labels, we chose the permutation leading to the smallest MSEE value.
\item Mean squared prediction error (MSPE) defined as 
\begin{equation*} 
 MSPE_a = \mathbb{E}\left[(\bm{y^v-X^v\widehat{\theta}_a})'(\bm{y^v-X^v\widehat{\theta}_a})\right]\text{,}
\end{equation*}
where $\bm y^v$ and $\bm X^v$ are respectively a vector of responses and a design matrix from a validation dataset.
\item Maximum log-likelihood (ML) reached. This quantity was approximated using 1,000 samples from $p(\bm Z|\bm y;\widehat{\bm \theta})$.
\end{enumerate}
Three versions of the MCEM algorithm were proposed for comparison with the SEM algorithm, depending on the number $M$ (or \code{nsamp}) of Gibbs iterations used to approximate the \textit{E step}. That number was varied between 5, 25 and 125. We chose these iterations numbers so as to cover different
situations, from a situation in which the number of iterations is too small to a situation in which that
number seems sufficient to expect having reached the convergence of the simulated Markov chain. Those versions were respectively denoted MCEM$_5$, MCEM$_{25}$ and MCEM$_{125}$. The comparison was performed using 200 simulated datasets. In order to consider high-dimensional situations with sizes
allowing to reproduce multiple simulations in a reasonable time, each training dataset consisted of $n=25$ individuals and $p=50$ variables. Validation datasets used to calculate MSPE consisted of 1,000 individuals each. All covariates were simulated  independently according to the standard Gaussian distribution:
\begin{equation*}
 \forall(i,j)\text{ }x_{ij}\sim\mathcal{N}(0,1)\text{.}
\end{equation*}
Both training and validation datasets were simulated according to model (\ref{eq:clere0}) using $\beta_0 = 0$, $\bm b = (0,3,15)'$, $\bm \pi = (0.64,0.20,0.16)'$, $\sigma^2=1$ and $\gamma^2=0$. This is equivalent to simulate data according to the standard linear regression model defined by:
\begin{equation*}
 y_i \sim \mathcal{N}\left(\sum^{32}_{j=1}{0 \times x_{ij}} + \sum^{42}_{j=33}{3 \times x_{ij}} + \sum^{50}_{j=43}{15 \times x_{ij}},1\right)
\end{equation*}
All algorithms were run using 10 different random starting points. Estimates yielding the largest likelihood were then used for the comparison.
\subsubsection{Results of the comparison}
Table \ref{tab:simulations} summarizes the results of the comparison
between the algorithms. The MCEM$_{5}$ algorithm was 1.3-fold faster
than the SEM algorithm however the latter algorithm poorly performed
regarding all other performance criteria (estimation quality,
prediction error, likelihood maximization). This observation
illustrates the importance of setting a large number $M$ of draws to
improve the estimation. Indeed, when increasing this number to 25 or
125, we observed an improvement in the estimation accuracy but no
noticeable improvement in the likehood. In turn, the SEM algorithm
was quite efficient compared to  MCEM$_{25}$ and  MCEM$_{125}$
algorithms. This algorithm not only ran faster (between 3 and 13-fold
faster than  MCEM$_{25}$ and  MCEM$_{125}$ - see
Table \ref{tab:simulations}) but also reached predictive performances close
to the oracle (i.e. using the true parameter). Those good performances
were mainly explained by the fact that the SEM algorithm most of the
time (66.5\% of the time) reached a better likelihood than the other
algorithms.\\
The results of this simulation study were made available
  as an internal dataset named \code{algoComp} in the \textbf{R} package
  \textbf{clere}. More details can be obtained using the command
  \code{help(algoComp)}.
}


<<echo=FALSE,results=tex>>=
library(clere)
data(algoComp)
nsim <- 200
lnames <- names(algoComp)
for (i in 1:length(lnames)) {
  Tmp <- algoComp[[i]][, 1:4]
  colnames(Tmp) <- c("MCEM5", "MCEM25", "MCEM125", "SEM")
  if (lnames[i] == "Liks") {
    tmp <- format(100 * table(factor(apply(Tmp, 1, which.max), levels = 1:4)) / nrow(Tmp), digit = 2, nsmall = 1)
    names(tmp) <- colnames(Tmp)
  } else {
    tmp <- 100 * table(factor(apply(Tmp, 1, which.min), levels = 1:4)) / nrow(Tmp)
    names(tmp) <- colnames(Tmp)
  }
  assign(paste("isBest", lnames[i], sep = ""), tmp)
  ## Mean
  avTmp <- format(apply(Tmp, 2, median), digit = 2, nsmall = 1)
  sdTmp <- format(1.253 * apply(Tmp, 2, sd) / sqrt(nsim), digit = 2, nsmall = 1)
  assign(paste("av", lnames[i], sep = ""), avTmp)
  assign(paste("sd", lnames[i], sep = ""), sdTmp)
}
avTPpred <- format(median(algoComp$Pred[, "ORACLE"]), digit = 2, nsmall = 1)
sdTPpred <- format(1.253 * sd(algoComp$Pred[, "ORACLE"]) / sqrt(nsim), digit = 2, nsmall = 1)
avTPliks <- format(median(algoComp$Liks[, "ORACLE"]), digit = 2, nsmall = 1)
sdTPliks <- format(1.253 * sd(algoComp$Liks[, "ORACLE"]) / sqrt(nsim), digit = 2, nsmall = 1)
 
tab <- paste(
  "\\begin{center}\n",
  "\\begin{table}[h!]\n",
  "\\begin{tabular}{llrr}",
  "\\toprule\n",
  "              &                  & \\% of times                    & Median       \\\\\n",
  "Performance indicators & Algorithms & the algorithm was best       & (Std. Err.)\\\\\n",
  "\\midrule\n",
  "CT (seconds)  &  SEM             &  ", isBestTime["SEM"], "           &  ", avTime["SEM"], " ( ", sdTime["SEM"], " ) \\\\\n",
  "              &  MCEM$_5$        &  ", isBestTime["MCEM5"], "         &  ", avTime["MCEM5"], " ( ", sdTime["MCEM5"], " ) \\\\\n",
  "              &  MCEM$_{25}$     &  ", isBestTime["MCEM25"], "        &  ", avTime["MCEM25"], " ( ", sdTime["MCEM25"], " ) \\\\\n",
  "              &  MCEM$_{125}$    &  ", isBestTime["MCEM125"], "       &  ", avTime["MCEM125"], " ( ", sdTime["MCEM125"], " ) \\\\\n",
  "\\\\\n",
  "\\midrule\n",
  "MSEE          &  SEM             &  ", isBestBias["SEM"], "           &  ", avBias["SEM"], " ( ", sdBias["SEM"], " ) \\\\\n",
  "              &  MCEM$_5$        &  ", isBestBias["MCEM5"], "         &  ", avBias["MCEM5"], " ( ", sdBias["MCEM5"], " ) \\\\\n",
  "              &  MCEM$_{25}$     &  ", isBestBias["MCEM25"], "        &  ", avBias["MCEM25"], " ( ", sdBias["MCEM25"], " ) \\\\\n",
  "              &  MCEM$_{125}$    &  ", isBestBias["MCEM125"], "       &  ", avBias["MCEM125"], " ( ", sdBias["MCEM125"], " ) \\\\\n",
  "\\\\\n",
  "\\midrule\n",
  "MSPE          &  SEM             &  ", isBestPred["SEM"], "           &  ", avPred["SEM"], " ( ", sdPred["SEM"], " ) \\\\\n",
  "              &  MCEM$_5$        &  ", isBestPred["MCEM5"], "         &  ", avPred["MCEM5"], " ( ", sdPred["MCEM5"], " ) \\\\\n",
  "              &  MCEM$_{25}$     &  ", isBestPred["MCEM25"], "        &  ", avPred["MCEM25"], " ( ", sdPred["MCEM25"], " ) \\\\\n",
  "              &  MCEM$_{125}$    &  ", isBestPred["MCEM125"], "       &  ", avPred["MCEM125"], " ( ", sdPred["MCEM125"], " ) \\\\\n",
  "              &  True parameter  &  ---                             &  ", avTPpred, " (", sdTPpred, " )\\\\\n",
  "\\\\\n",
  "\\midrule\n",
  "ML            &  SEM             &  ", isBestLiks["SEM"], "           &  ", avLiks["SEM"], " ( ", sdLiks["SEM"], " ) \\\\\n",
  "              &  MCEM$_5$        &  ", isBestLiks["MCEM5"], "         &  ", avLiks["MCEM5"], " ( ", sdLiks["MCEM5"], " ) \\\\\n",
  "              &  MCEM$_{25}$     &  ", isBestLiks["MCEM25"], "        &  ", avLiks["MCEM25"], " ( ", sdLiks["MCEM25"], " ) \\\\\n",
  "              &  MCEM$_{125}$    &  ", isBestLiks["MCEM125"], "       &  ", avLiks["MCEM125"], " ( ", sdLiks["MCEM125"], " ) \\\\\n",
  "              &  True parameter  &  ---                             &  ", avTPliks, " (", sdTPliks, " )\\\\\n",
  "\\bottomrule\n",
  "\\end{tabular}\n",
  "\\caption{\\label{tab:simulations} Performance indicators used to compare SEM and MCEM algorithms. Computational Time (CT) was measured on a Intel(R) Xeon(R) CPU E7- 4870  @ 2.40GHz processor. The best algorithm is defined as the one that either reached the largest log-likelihood (ML) or the lowest CT, Mean Squared Prediction Error (MSPE) and Mean Squared Estimation Error (MSEE).}\n",
  "\\end{table}\n",
  "\\end{center}\n"
)
cat(tab)
@
\subsection{Comparison with other methods}
\label{sec:simulation}
\subsubsection{Description of the methods}
{
In this section, we compare our model to standard dimension reduction
  approaches in terms of MSPE. Although a more detailed comparison was proposed in [\cite{yengo2013}], we propose here a quick illustration of the relative predictive performance of our model. The comparison is achieved using data simulated according 
  to the scenario described above in Section \textbf{\textit{Description of the simulation study}}.%\ref{sec:semvsmcem}.
  The methods selected for comparison are %the variable   selection using LARS algorithm [\cite{efron2004}],
  the ridge  regression [\cite{hoerl1970}], the elastic net [\cite{zou2005}], the LASSO
  [\cite{tibshirani1996}], PACS [\cite{pacs2013}], the method of Park and
  colleagues [\cite{park2007}] (subsequently denoted AVG) and the spike
  and slab model [\cite{ishwaran2005}] (subsequently denoted SS). The
  first three methods are implemented in the freely available
  \textbf{R} package \textbf{glmnet}. With the latter package, the tuning
  parameter \texttt{lambda} was selected using the function
  \texttt{cv.glm} aiming at minimizing the mean squarred error (option
  \texttt{type="mse"}). In particular for the Elastic net, the
  second tuning parameter \texttt{alpha} (measuring of the amount of
  mixture between the $L^1$ and $L^2$ penalty) was jointly selected
  with \texttt{lambda} to minimize the mean squarred error. The
  \textbf{R} package \textbf{glmnet} proposes a procedure for
  automatically selecting values for \texttt{lambda}. We therefore
  used this default procedure while we selected \texttt{alpha} among  $\{0, 0.1,0.2,\ldots,0.9,1\}$.\\
   PACS methodology proposes to estimate the regression coefficients
  by solving a penalized least squares problem. It imposes a
  constraint on $\bm \beta$ that is a weighted combination of the
  $L^1$ norm and the pairwise $L^\infty$ norm. Upper-bounding the
  pairwise  $L^\infty$ norm enforces the covariates to have close
  coefficients. When the constraint is strong enough, closeness 
  translates into equality achieving thus a grouping property. For
  PACS, no software was available. Only an \textbf{R} script was
  released on Bondell's webpage\footnote{\url{http://www4.stat.ncsu.edu/~bondell/Software/PACS/PACS.R.r}}.\\ 
  Since this \textbf{R} script was running very slowly, we decided to
  reimplement it in \textbf{C++} and observed a 30-fold speed-up of computational time. Similarly to Bondell's
  \textbf{R} script, our implementation uses two parameters \texttt{lambda} and
  \texttt{betawt}. Our reimplementation of Bondell's script was
  included in the \textbf{R} package \textbf{clere} under the function
  \texttt{fitPacs()}. In [\cite{pacs2013}], the authors suggest assigning \texttt{betawt} with
  the coefficients obtained from a ridge regression model after the
  tuning parameter was selected using AIC. In this simulation study we
  used the same strategy; however the ridge parameter was selected via
  5-fold cross validation. 5-fold CV was preferred to AIC since
  selecting the ridge parameter using AIC always led to estimated
  coefficients equal to zero. Once \texttt{betawt} was selected,
  \texttt{lambda} was chosen via 5-fold cross validation among the
  following values: 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1, 2, 5, 10, 20,
  50, 100, 200 and 500. All other default parameters of their script  were unchanged. \\
  The AVG method is a two-step approach. The first step uses hierarchical
  clustering of covariates to create surrogate covariates by averaging the variables within each
  group. Those new predictors are afterwards included in a linear regression
  model, replacing the primary variables. A variable selection
  algorithm is then applied to select the most predictive groups of
  covariates. To implement this method, we followed the algorithm
  described in [\cite{park2007}] and programmed it in \textbf{R}.\\
  The spike and slab model is a Bayesian approach for variable
  selection. It is based on the assumption that the regression
  coefficients are distributed according to a mixture of two centered
  Gaussian distributions with different variances. One component of
  the mixture (the spike)  is chosen to have a small variance, while
  the other component (the slab) is allowed to have a large
  variance. Variables assigned to the spike are dropped from the
  model. We used the \textbf{R} package \textbf{spikeslab} to run the spike and
  slab models. Especially, we used the function \textbf{spikeslab}
  from that package to detect influential variables. The number of
  iterations used to run the function \textbf{spikeslab} was 2,000   (1,000 discarded). \\ 
  When running \code{fitClere()}, the number \code{nItEM} of SEM
  iterations was set to 2,000. The number \code{g} of groups for CLERE
  was chosen between 1 and 5 using AIC (option
  \code{analysis="aic"}). Two versions of CLERE were considered: the
  one with all parameters estimated and the one with $b_1$ set to
  0. The latter approach is subsequently denoted CLERE$_0$ (option \code{sparse=TRUE}).
} 
  
\subsubsection{Results of the comparison}
{
  Figure \ref{fig:sims}, summarizes the comparison between the
  methods. In this simulated scenario, CLERE outperformed the other
  methods in terms of prediction error. Those good performances were
  improved when parameter $b_1$ was set to 0. CLERE was also the most
  parsimonous approach with an averaged number of estimated parameters
  equal to 7.7 (6.9 when $b_1=0$). The second best approach was PACS
  which also led to parsimonous models. The superiority of such methods could be expected since
in the simulation model the regression coefficients are gathered in three groups. Variable selection approaches
  as whole yielded the largest prediction error in this
  simulation. CLERE, PACS and Spike and Slab had the largest computational times (CT).
  For CLERE and PACS  this loss in CT were compensated by a a strong improvement in
  prediction error as explained above, while Spike and Slab yielded
  the worst prediction error in addition of being the slowest approach
  in this example.\\
The results of this simulation study were made available
  as an internal dataset named \code{numExpSimData} in the \textbf{R} package
  \textbf{clere}. More details can be obtained using the command
  \code{help(numExpSimData)}.
\begin{figure}[!ht]
\centerline{
<<echo=FALSE>>=
library(clere)
data(numExpSimData)
meths <- c(
  "CLERE0", "CLERE", "PACS", "LASSO",
  "AVG", "Ridge", "Elastic net",
  "Spike and Slab"
)
o <- order(apply(numExpSimData[, 1:9], 2, median))
dfs <- round(apply(numExpSimData[, 10:18][, o[1:8]], 2, mean), 1)
sdf <- round(apply(numExpSimData[, 10:18][, o[1:8]], 2, sd), 1)
tts <- round(apply(numExpSimData[, 19:27][, o[1:8]], 2, mean), 1)
stt <- round(apply(numExpSimData[, 19:27][, o[1:8]], 2, sd), 2)
cols <- rainbow(9)
par(mar = c(5, 2, 4, 7) + 0.1)
boxplot(numExpSimData[, 1:9][, o[1:8]], horizontal = TRUE, log = "x", col = cols, axes = FALSE, pch = 18, xlab = "Mean Squared Prediction Error")
axis(1)
labs <- paste(meths, "\ndf: ", dfs, " (", sdf, ")", sep = "")
axis(4, at = 1:8, labels = labs, las = 2)
cts <- paste(tts, "s (", stt, ")", sep = "")
legend("topleft", legend = cts, box.lty = 0, lwd = 2, lty = 1, col = cols, title = "Computational time")
@
%\includegraphics[scale=0.8]{images/Simulations.pdf}
}
\caption{\label{fig:sims} Comparison between CLERE and some
standard dimension reduction approaches. The number of estimated parameters (df: +/- standard error) is given in the right along with the name of the method utilized. The average computational time with its corresponding standard error (given in parenthesis) is also provided in each situation.}
\end{figure}
}

\subsection{Real datasets analysis}
\label{sec:realdata}
\subsubsection{Description of the datasets}
  We used in this section the real datasets \code{Prostate} and
  \code{eyedata} from the \textbf{R} packages \textbf{lasso2} and
  \textbf{flare} respectively. The \code{Prostate} dataset comes from
  a study that examined the correlation between the level of prostate
  specific antigen and a number of clinical measures in $n=97$ men who
  were about to receive a radical prostatectomy. This dataset is a benchmark dataset used in multiple publications
about high-dimensional regression model, including [\cite{tibshirani1996}] and was chosen
here in order to illustrate the performances of CLERE in comparaison of the competitor methods.  We used the prostate specific antigen
  (variable \code{lpsa}) as response variable and the $p=8$ other
  measurements as covariates.    
  The \code{eyedata} dataset is extracted from the published study of
    [\cite{scheetz2006regulation}]. This dataset consists in gene
expression levels measured at $p=200$ probes in $n=120$ rats. The
response variable utilized was the expression of the TRIM32 gene which
is a biomarker of the Bardet-Bidel Syndrome (BBS). {We chose this
  dataset to illustrate the performances of CLERE on a (manageable)
  high-dimensional problem which is the actual context for which this method
  was developped [\cite{yengo2013}]}.

  Those two datasets was utilized to compare CLERE to the same methods used
  above in the Section presenting the simulation stydy. All 
  methods were compared in term of out-of-sample prediction 
  error estimated using 5-fold cross-validation (CV). Those CV
  statistics were then averaged and compared across the methods in
  Table \ref{tab:realdata}.
 
\subsubsection{Running the analysis}
\par{
Before presenting the results of the comparison between CLERE and its competitors, we illustrate the command lines to run the analysis of the \code{Prostate} dataset.    The dataset is loaded by typing:
%<<options,echo=FALSE>>=
% options(prompt = "R> ", continue = " ",width = 85)
%@
<<>>=
library(clere)
data(Prostate, package = "lasso2")
y <- Prostate[, "lpsa"]
x <- as.matrix(Prostate[, -which(colnames(Prostate) == "lpsa")]) 
@
Possible training (\code{xt} and \code{yt}) and validation (\code{xv} and \code{yv}) sets are generated as follows:
<<>>=
itraining <- 1:(0.8 * nrow(x))
xt <- x[ itraining, ]
yt <- y[ itraining]
xv <- x[-itraining, ]
yv <- y[-itraining]
@
The \code{fitClere()} function is run using AIC criterion to select the number of groups between 1 and 5. To lessen the impact of local minima in the model selection, 5 random starting points are used. This can be implemented as written below
<<>>=
Seed <- 1234
mod <- fitClere(
  y = yt, x = xt, g = 5, analysis = "aic", parallel = FALSE, nstart = 5,
  sparse = TRUE, nItEM = 2000, nBurn = 1000, nItMC = 10, dp = 5, nsamp = 1000,
  seed = Seed
)
summary(mod)
@
\begin{figure}[!ht]    
\centerline{
<<echo=FALSE>>=
plot(mod)
@
%\includegraphics[scale=0.8]{images/Figure1.pdf}
}
\caption{\label{fig:plot} Values of the model parameters in view of SEM algorithm iterations. The vertical grey line in each of the four plots, represents the number \code{nBurn} of iterations discarded before calculating maximum likelihood estimates.}
\end{figure}
Running the command \code{plot(mod)} generates the plot given in Figure \ref{fig:plot}. We can also access the cluster membership by running the command \code{clusters()}. For example, running the command \code{clusters(mod,threshold=0.7)} yields
<<>>=
clusters(mod, thresold = 0.7)
@
In the example above 2 variables, being the cancer volume (\code{lcavol}) and the prostate weight (\code{lweight}), were assigned to group 2 ($b_2=0.4737$). The other 6 variables were assigned to group 1 ($b_1=0$). Posterior probabilities of membership are available through the slot \code{P} in object of class \code{Clere}.
<<>>=
mod@P
@ 
The covariates were respectively assigned to their group with a probability larger than 0.7. Moreover, given that parameter $\gamma^2$ had very small value ($\widehat{\gamma^2} = 4.065\times10^{-8}$), we can argue that cancer volume and prostate weight are the only relevant explanatory covariates. To assess the prediction error associated with the model we can run the command \code{predict()} as follows: 
<<>>=
error <- mean((yv - predict(mod, xv))^2)
error
@
}
\subsubsection{Results of the analysis}
\par{
Table \ref{tab:realdata} summarizes the prediction errors and the
number of parameters obtained for all the methods. CLERE$_0$ had the
lowest prediction error in the analysis of the \code{Prostate} dataset
and the second best performance with the \code{eyedata} dataset. The
AVG method was also very competitive compared to variable selection
approaches stressing thus the relevance of creating groups of
variables to reduce the dimensionality (especially in the
\code{eyedata} dataset). It is worth noting that in both datasets,
imposing the constraint $b_1=0$ improved the predictive performance of CLERE.

In the \code{Prostate} dataset, CLERE robustly identified two groups
of variables representing influential ($b_2>0)$ and not relevant
variables ($b_1=0$). In the \code{eyedata} dataset in turn, AIC led to
select only one group of variables. However, this did not lessened the
predictive performance of the model since CLERE$_0$ was second best
after AVG, while needing significantly less parameters. PACS low
performed in both datasets. The Elastic net showed good predictive
performances compared to the variable selection methods like LASSO or
spike and slab model. Ridge regression and Elastic net had comparable
results in both datasets. In line with the results of the simulation
study, we observed that despite the a larger computational time (CT), CLERE
and CLERE$_0$ had a reduced mean squarred error compared to the fastest
methods. However, this improvement was less substantial than observed
in the simulation study given the differences in CT. This increased CT
may be explained by the fact that no simple stopping rule is proposed
when fitting CLERE. We may therefore consider that a smaller number of
SEM iterations could have been used to yield a similar prediction
error. Indeed, when looking at Figure \ref{fig:plot}, we see that the
convergence was achieved almost from the first 10 iterations. Still,
the observed CT for CLERE being around 22s for the \code{eyedata}
dataset and around 3s for the \code{Prostate} dataset remains
affordable in these examples.
\\

The results of this analysis on real data were made available
  as an internal dataset named \code{numExpRealData} in the \textbf{R} package
  \textbf{clere}. More details can be obtained using the command
  \code{help(numExpRealData)}.

% Directly copied from the script RealDataExample.R: %% cat(createTable2(numExpRealData))
<<echo=FALSE,results=tex>>=
library(clere)
data(numExpRealData)
pData <- numExpRealData[["ProstateData"]][, 1:27]
eData <- numExpRealData[["EyeData"]][, 1:27]
basenames <- colnames(pData)[1:9]
colnames(pData) <- paste(basenames, rep(c("cv", "df", "tt"), each = 9), sep = "_")
colnames(eData) <- paste(basenames, rep(c("cv", "df", "tt"), each = 9), sep = "_")
methname <- basenames
names(methname) <- methname
methname["ELNET"] <- "Elastic net"
methname["CLERE_s"] <- "CLERE$_0$"
methname["SS"] <- "Spike and Slab"
genLine <- function(refData, meth) {
  cv <- 100 * refData[, paste(meth, "_cv", sep = "")]
  df <- refData[, paste(meth, "_df", sep = "")]
  tt <- refData[, paste(meth, "_tt", sep = "")]
  line <- paste(
    methname[meth], "      &  ", format(mean(cv), digit = 2, nsmall = 1), " (", format(sd(cv) / sqrt(5), digit = 1), " )",
    "                     &  ", format(mean(df), digit = 2, nsmall = 1), " (", format(sd(df) / sqrt(5), digit = 1), " )",
    "                     &  ", format(mean(tt), digit = 2, nsmall = 1), " (", format(sd(tt) / sqrt(5), digit = 1), " ) \\\\\n"
  )
  return(line)
}
tab <- paste(
  "\\tiny\n",
  "\\begin{center}\n",
  "\\begin{table}[h!]\n",
  "\\begin{tabular}{lccc}\n",
  "\\toprule\n",
  "                    & 100$\\times$Averaged CV-statistic & Number of parameters & CT (seconds)\\\\\n",
  "                    & (Std. Error)                     & (Std. Error)         & (Std. Error) \\\\\n",
  "\\midrule\n",
  "\\code{Prostate dataset}                   &                       &                               &  \\\\\n",
  "\\midrule\n",
  paste(sapply(basenames, function(meth) genLine(pData, meth)), collapse = ""),
  "\\\\\n",
  "\\midrule\n",
  "\\code{eyedata}                    &                       &                               &  \\\\\n",
  "\\midrule\n",
  paste(sapply(basenames, function(meth) genLine(eData, meth)), collapse = ""),
  "\\bottomrule\n",
  "\\end{tabular}\n",
  "\\caption{\\label{tab:realdata} Real data analysis. Out-of-sample prediction error (averaged CV-statistic) was estimated using cross-validation in 100 splitted datasets. The number of parameters reported for CLERE/CLERE$_0$ was selected using AIC. CT stands for the average Computational Time.}\n",
  "\\end{table}\n",
  "\\end{center}\n",
  "\\normalsize\n"
)
cat(tab)
@
}
\section[Conclusions]{Conclusions}
\label{sec:conc}
We presented in this paper the \textbf{R} package \textbf{clere}. This
package implements two efficient algorithms for fitting the
CLusterwise Effect REgression model: the MCEM and the SEM
algorithms. If the MCEM algorithm is to be preferred when $p<n$, the
SEM algorithm is more efficient for high dimensional datasets
($n<p$). The good performances of SEM over MCEM could have been
expected regarding the computational complexities of the two
algorithms that are $\mathcal{O}\left( npg + g^3 + N_{max}ng\right)$
and $\mathcal{O}\left( M(p^2 + pg)\right)$ respectively. In fact, as
long as $p>n$, the SEM algorithm has a lower complexity. However, the
computational time to run our SEM algorithm is more variable compared
to MCEM as its \text{M step} does not have a closed form. We finally
advocate the use the MCEM algorithm only when $p\ll n$. 
Another important feature for model interpretation is proposed by
constraining the model parameter $b_1$ to equal 0, which leads to
carry out variable selection. Such constraint may also lead to a
reduced prediction error. We illustrated on a real dataset, how to run
an analysis using a detailed \textbf{R} script. Although our numerical
experiments showed that the CLERE method tended to be slower than
variable selection methods, it still had better or competitive
predictive performances. In addition, the CLERE model was often more
parsimonious than other models and was easily interpretable since
groups of regression coefficients/variables could be summarized using
a single parameter.\\
Our model can easily be extended to the analysis of binary
responses. This extension will be proposed in forthcoming version of
the package. Another direction for future research will be to develop
an efficient stoping rule for the proposed SEM algorithm, specific to
our context. Such a criterion is expected to improve the computational
performances of our estimation algorithm.

\bibliography{biblio}

\address{Loic Yengo\\
  Integrated Genomics and Metabolic Diseases Modeling\\
  CNRS UMR 8199 - Lille Institute of Biology\\
  E.G.I.D - FR3508 European Genomics Institute of Diabetes\\
  1, rue du Professeur Calmette, BP 447, 59021 Lille cedex\\
  France}
  \email{loic.yengo@cnrs.fr}

\address{Julien Jacques\\
  ERIC Laboratory\\
  University of Lyon - Lumière\\
  5 avenue Pierre Mendès France, 69676 Bron cedex\\
  France}
\email{julien.jacques@univ-lyon2.fr}

\address{Christophe Biernacki\\
  MODAL team (Inria) \& Laboratoire Paul Painlev\'e (UMR CNRS 8524)\\
  University Lille I\\
  Cité Scientifique, 59655 Villeneuve d'Ascq cedex\\
  France}
\email{christophe.biernacki@math.univ-lille1.fr}

\address{Mickael Canouil\\
  Integrated Genomics and Metabolic Diseases Modeling\\
  CNRS UMR 8199 - Lille Institute of Biology\\
  E.G.I.D - FR3508 European Genomics Institute of Diabetes\\
  1, rue du Professeur Calmette, BP 447, 59021 Lille cedex\\
  France}
  \email{mickael.canouil@cnrs.fr}

\end{article}

\end{document}
