# Simultaneous Variable Selection, Clustering, and Smoothing in Function on
Scalar Regression

Suchit Mehrotra  Department of Statistics, North Carolina State University
Arnab Maity  Department of Statistics, North Carolina State University

###### Abstract

We address the problem of multicollinearity in a function-on-scalar regression
model by using a prior which simultaneously selects, clusters, and smooths
functional effects. Our methodology groups effects of highly correlated
predictors, performing dimension reduction without dropping relevant
predictors from the model. We validate our approach via a simulation study,
showing superior performance relative to existing dimension reduction
approaches in the function-on-scalar literature. We also demonstrate the use
of our model on a data set of age specific fertility rates from the United
Nations Gender Information database.

##  1 Introduction

Recent technological advances in computing power and data storage have
simplified the collection of vast quantities of data. Modern scientific
questions often involve statistical analyses on datasets with a large number
of predictors, many of which have no meaningful effect on the response. In
numerous situations, an analysis is further complicated by multicollinearity
in predictors. An example are microarray studies where interest lies in
understanding the relationship between a health outcome and gene expressions;
gene expression levels are often highly correlated and the number of genes is
substantially larger than the number of samples. In such situations, modern
statistical approaches have utilized the assumption that highly correlated
predictors share the same effect on the response, especially if the
correlation is due to the same underlying factor. They aim to reduce the
dimension of the problem by summing or averaging groups of highly correlated
columns and fitting a linear model with the new set of predictors, which, in
essence, is equivalent to finding clusters of regression coefficients in a
supervised manner.

Many methods exist which aim to cluster coefficients in the univariate linear
model. For example, Bondell and Reich (2008), She (2010), and Sharma et al.
(2013) all propose methods which minimize a loss function plus a penalty term
to encourage clustering in the coefficients. These methods are computationally
efficient as they are all a special case of the generalized lasso (Tibshirani
and Taylor, 2011). In the Bayesian setting, the coefficient clustering problem
is solved by utilizing mixture models as priors for the coefficients. Tadesse
et al. (2005) uses finite mixture models to cluster the coefficients, while
Nott (2008), Kim et al. (2006), Dunson et al. (2008), and Curtis and Ghosh
(2011) all use Dirichlet Process priors to estimate clusters.

In this paper we add to the existing dimension reduction literature for models
with a functional response and a set of scalar predictors, which up till now
has focused only on variable selection. At present, many group penalties exits
to select variables in a function-on-scalar regression. First, the group lasso
(Yuan and Lin, 2006) can be used to shrink basis function coefficients
belonging to the same predictor simultaneously to zero. Wang et al. (2007)
modify this approach by using a group SCAD penalty, Chen et al. (2016) use a
group MCP penalty, Fan and Reimherr (2017) propose an adaptive group lasso,
while Kowal and Bourgeois (2018) utilize a group horseshoe prior. However, an
issue with these approaches is the need to fix the number of basis functions
before fitting the model, and recent work has focused on smoothing the
coefficient functions along with selecting variables. For example, Goldsmith
and Kitago (2016) use a penalized normal prior to smooth coefficient functions
without any variable selection, while Parodi and Reimherr (2018) propose a
framework for simultaneous variable selection and smoothing.

Currently, to the best of our knowledge, no function-on-scalar approaches
exist which address multicollinearity in the predictors by exploiting a
cluster structure in the coefficient functions. Additionally, in the Bayesian
setting, methods to simultaneously select and smooth coefficient functions
have yet to be developed. Consequently, our manuscript contributes to the
current function-on-scalar literature by proposing a novel clustering based
dimension reduction technique for functional data, heretofore unexplored.

Our paper proceeds as follows. We review the relevant background literature in
Section 2, our model and relevant computational details in Sections 3 and 4, a
simulation study and its results in Section 5, along with a real data
application in Section 6.

##  2 Background

In a univariate linear model, an approach to finding clusters in a vector of
coefficients, ${{\bm{\mathbf{{\beta}}}}}$, is to assume the existence of $K$
clusters and use finite mixture models as a prior for
${{\bm{\mathbf{{\beta}}}}}$ (Tadesse et al., 2005). Using a latent class
representation for cluster membership the prior can be written as:

$\displaystyle\begin{split}(\forall p\in\\{1,\dots,P\\})\
\beta_{p}|c_{p},{{\bm{\mathbf{{\theta}}}}}&=\sum_{k=1}^{K}I(c_{p}=k)\theta_{k}\\\
(\forall p\in\\{1,\dots,P\\})\ c_{p}|{{\bm{\mathbf{{\pi}}}}}&\sim\text{
Discrete}(\pi_{1},\dots,\pi_{K})\\\ (\forall k\in\\{1,\dots,K\\})\
\theta_{k}|\lambda&\sim\mathcal{N}(0,\lambda^{-1})\end{split}$ (1)

where ${\bf c}$ is a $P\times 1$ vector which stores class labels for the
elements of ${{\bm{\mathbf{{\beta}}}}}$. The hierarchy in (1) implies that
$\beta_{p}\sim G=\sum_{k=1}^{K}\pi_{k}\delta_{\theta_{k}}(\beta_{p})$ where
$\delta_{\theta_{k}}(x)$ is 1 if $x=\theta_{k}$ and $0$ otherwise.
Consequently, each element of ${{\bm{\mathbf{{\beta}}}}}$ is a finite mixture
of K delta functions, and can only take on one of $K$ values. If we know $K$
in advance, estimating ${{\bm{\mathbf{{\theta}}}}}$ and class memberships,
${\bf c}$, provides a solution to the coefficient clustering problem.

However, this approach suffers from a major drawback in that $K$ has to be
fixed a-priori, which is information rarely available in applied settings. To
remedy this issue, one can use priors which allow for a countably infinite
number of mixture components, such as the Dirichlet Process.

###  2.1 Dirichlet Process Priors

The Dirichlet Process (DP) is indexed by a concentration parameter $\alpha$
and a base distribution $G_{0}$, denoted $DP(G_{0},\alpha)$. It has a long
history of being used as a clustering prior for a set of points,
$\\{x_{1},\dots x_{N}\\}$, and was first investigated by Ferguson (1973) and
Antoniak (1974). More recently, computationally efficient methods have been
developed for parameter estimation by Escobar and West (1995), MacEachern and
Müller (1998), and Neal (2000). Using a stick breaking construction of the
Dirichlet process, it can be shown that
$G(x)=\sum_{k=1}^{\infty}\pi_{k}\delta_{\theta_{k}}(x)$ is a sample from a
$DP(G_{0},\alpha)$, where $\pi_{k}=V_{k}\prod_{l=1}^{k-1}(1-V_{l})$,
$V_{k}\overset{iid}{\sim}\text{Beta}(1,\alpha)$, and
$\theta_{k}\overset{iid}{\sim}G_{0}$ (Sethuraman, 1994). Hence, the
distributions sampled from a Dirichlet Process are a countably infinite
mixture of point masses.

Key insights into the behaviour of the Dirichlet Process were given by
Blackwell et al. (1973) and Neal (2000). Neal (2000) derives a conditional
prior for the class labels, ${\bf c}$, by viewing the Dirichlet Process
mixture model as the limit of finite mixture models as $K\to\infty$. By
integrating out the mixing parameters ${{\bm{\mathbf{{\pi}}}}}$ he gives the
following conditional priors for the elements of ${\bf c}$

$\displaystyle\begin{split}P(c_{i}=k|{\bf
c}_{-i})=\frac{n_{-i,k}}{N-1+\alpha}\\\ P(c_{i}\neq c_{j}\text{ for all }i\neq
j|c_{-i})=\frac{\alpha}{N-1+\alpha}\end{split}$ (2)

where $n_{-i,k}$ is the number of $c_{j}$, $j\neq i$, that are equal to $k$
and ${\bf c}_{-i}$ is ${\bf c}$ without the $i^{th}$ element. This
representation allows us to bypass the sampling of mixing weights,
${{\bm{\mathbf{{\pi}}}}}$, and shows that the prior conditional probability of
an observation being assigned to an existing cluster is proportional to the
number of elements in that cluster.

For the univariate linear model, we can modify the hierarchy in (1) for the
Dirichlet Process prior as follows:

$\displaystyle\begin{split}(\forall p\in\\{1,\dots,P\\})\ \beta_{p}|G&\sim
G\\\ G&\sim DP(G_{0},\alpha)\\\
G_{0}|\lambda&=\mathcal{N}(0,\lambda^{-1})\end{split}$ (3)

Nott (2008) explored the estimation and predictive capacity of the DP prior
with a normal base measure along with its applications in penalized spline
smoothing. He also derived a Gibbs sampler for the model in (3) by leveraging
a latent class representation for group membership. His algorithm iterates
between sampling the cluster indicators for the elements of
${{\bm{\mathbf{{\beta}}}}}$, ${\bf c}$, and all other model parameters
conditioned on the number of clusters in ${\bf c}$, $K$. We leverage this
algorithm for development of our own Gibbs samplers and will discuss its
relevant details in Section 4.

###  2.2 Variable Selection and Clustering

Kim et al. (2006), Dunson et al. (2008), and Curtis and Ghosh (2011) propose
extensions to the prior in (3) for simultaneous variable selection and
clustering. They use a mixture of a point mass at zero with a Dirichlet
process prior

$\displaystyle G=\pi_{0}\delta_{0}+(1-\pi_{0})G^{*},\ G^{*}\sim
DP(G_{0},\alpha)$ (4)

where $0\leq\pi_{0}\leq 1$ is a mixture weight, $\delta_{0}$ is a point mass
probability density at $0$, and $DP(G_{0},\alpha)$ is a Dirichlet Process
prior with precision $\alpha$ and base measure $G_{0}$. Using the constructive
definition of Sethuraman (1994) the prior in (4) is again an infinite mixture
of point masses with the first point mass fixed at zero

$\displaystyle
G=\pi_{0}\delta_{0}+(1-\pi_{0})\sum_{k=1}^{\infty}\pi_{k}\delta_{\theta_{k}}=\sum_{k=0}^{\infty}\tilde{\pi}_{k}\delta_{\theta_{k}}$
(5)

where $\tilde{\pi}_{0}=\pi_{0}$, $\tilde{\pi}_{k}=(1-\pi_{0})\pi_{k}$ if
$k\geq 1$, $\theta_{0}=0$ and $\theta_{k}\overset{iid}{\sim}G_{0}$ if $k\geq
1$. Computation for this model is similar to the model discussed in (3)
because the prior in (5) can be interpreted as a clustering prior with the
restriction that the first cluster is zero. We can account for this
restriction by making a few changes to the Gibbs sampler developed by Nott
(2008).

##  3 Model

In this section we extend the methods discussed in Section 2 to models with a
functional response and a set of scalar covariates. Let $i\in\\{1,\dots,N\\}$
be the subject level index, and $p\in\\{1,\dots,P\\}$ be the predictor index.
Then the general function-on-scalar model is:

$\displaystyle y_{i}(t)=\mu(t)+\sum_{p=1}^{P}x_{ip}\beta_{p}(t)+e_{i}(t)$ (6)

where $y_{i}(t)$ is the functional response value for subject $i$ at time $t$,
$\mu(t)$ is an intercept term, $x_{ip}$ is the $p^{th}$ predictor value for
the $i^{th}$ individual, $\beta_{p}(t)$ is the functional effect for the
$p^{th}$ predictor at time $t$, and $e_{i}(t)$ is an error term. If all
observations are observed on a common grid, this model can be written in
matrix form as: ${\bf Y}={\bf X}{{\bm{\mathbf{{\beta}}}}}^{\sf T}+{\bf E}$.

Since we wish to cluster only some of the predictors (at the least the
intercept will remain unclustered), we delineate two groups of functional
effects and propose the following model for a grid of equally spaced points:

$\displaystyle{\bf Y}={\bf W}{{\bm{\mathbf{{\alpha}}}}}^{\sf T}+{\bf
X}{{\bm{\mathbf{{\beta}}}}}^{\sf T}+{\bf E}$ (7)

where ${\bf Y}_{N\times T}$ is a matrix of responses, ${\bf W}_{N\times
P_{f}}$ is a matrix of predictors with free effects,
${{\bm{\mathbf{{\alpha}}}}}_{T\times P_{f}}$, ${\bf X}_{N\times P_{c}}$ is a
matrix of predictors whose effects, ${{\bm{\mathbf{{\beta}}}}}_{T\times
P_{c}}$, we wish to select and cluster, and ${\bf E}_{N\times T}$ is a matrix
of errors with the ith row ${\bf
e}_{i}\sim\mathcal{N}_{T}({\bm{\mathbf{{0}}}},\tau^{-1}{\bf I})$. We expand
each coefficient function as a linear combination of $M$ b-spline basis
functions, $\\{\theta_{1},\dots,\theta_{M}\\}$, such that
${{\bm{\mathbf{{\alpha}}}}}_{T\times
P_{f}}={{\bm{\mathbf{{\Theta}}}}}_{T\times M}{\bf A}_{M\times P_{f}}$ and
${{\bm{\mathbf{{\beta}}}}}_{T\times P_{c}}={{\bm{\mathbf{{\Theta}}}}}_{T\times
M}{\bm{\tilde{\mathbf{{B}}}}}_{M\times P_{c}}$. This reduces the functional
estimation problem to estimating the elements of ${\bf A}$ and
${\bm{\tilde{\mathbf{{B}}}}}$.

In practice, care must be taken in choosing $M$, balancing the need for a
function to remain flexible with the risk of overfitting. We address this
issue by expanding the functional effects using more basis functions than we
think are necessary and use penalization to induce smoothness in the
coefficient function. We follow Goldsmith and Kitago (2016) and use a full-
rank penalty matrix ${\bf R}$ to penalize basis coefficients, setting ${\bf
R}=\eta{\bf I}+(1-\eta){\bf R}_{2}$, where ${\bf I}$ is the identity matrix,
${\bf R}_{2}$ is a second degree P-spline penalty (Eilers and Marx, 1996), and
$\eta=0.001$. Consequently, the prior for the columns of ${\bf A}$, ${\bf
a}_{p}$, is $\forall p\in\\{1,\dots,P_{f}\\}$

$\displaystyle{\bf a}_{p}$
$\displaystyle\sim\mathcal{N}_{M}({\bm{\mathbf{{0}}}},\lambda_{a_{p}}^{-1}{\bf
R}^{-1})$ (8)

where $\lambda_{a_{p}}\sim\mathcal{G}(a_{\lambda},b_{\lambda})$ and
$\mathcal{G}(a,b)$ is the Gamma distribution with shape $a$ and rate $b$.

###  3.1 Clustering Functional Effects

Many approaches exist for clustering functional data and we refer the reader
to Jacques and Preda (2014) for a review. Because we expand each coefficient
function using b-splines, our clustering approach is similar to the work
initially proposed by Abraham et al. (2003); we find clusters in
${{\bm{\mathbf{{\beta}}}}}$ by clustering the coefficients of the basis
expansion. If we assume the columns of ${{\bm{\mathbf{{\beta}}}}}$ come from
$K$ clusters, we can expand it as
${{\bm{\mathbf{{\beta}}}}}={{\bm{\mathbf{{\Theta}}}}}{\bm{\tilde{\mathbf{{B}}}}}={{\bm{\mathbf{{\Theta}}}}}{\bf
B}{\bf C}^{\sf T}$ where ${{\bm{\mathbf{{\Theta}}}}}_{T\times M}$ and ${\bf
B}_{M\times K}$ are the basis function and coefficient matrices, respectively,
and ${\bf C}_{P_{c}\times K}$ is a matrix whose rows are one-hot encoded with
class membership information. Consequently, ${\bf
X}{{\bm{\mathbf{{\beta}}}}}^{\sf T}={\bf X}{\bm{\tilde{\mathbf{{B}}}}}^{\sf
T}{{\bm{\mathbf{{\Theta}}}}}^{\sf T}={\bf X}{\bf C}{\bf B}^{\sf
T}{{\bm{\mathbf{{\Theta}}}}}^{\sf T}$, and ${\bf C}$ has the effect of summing
columns of ${\bf X}$ with the same effect on the response.

Because we are also interested in estimating smooth coefficient functions, we
give each predictor its own smoothing parameter, $\tilde{\lambda}_{b_{p}}$.
Therefore, our task is to estimate a mean and smoothing parameter for each of
our clusters. We can easily do this by modifying our base distribution to be a
multivariate normal-gamma (Ray and Mallick, 2006; Zhang and Telesca, 2014).
Consequently, for clustering functional effects, we propose the prior:

$\displaystyle\begin{split}[\forall\in\\{1,\dots,P_{c}\\}]\ \
({\bm{\tilde{\mathbf{{b}}}}}_{p},\tilde{\lambda}_{b_{p}})&\sim G\\\ G&\sim
DP(G_{0},\alpha)\\\
G_{0}&=\mathcal{N}_{M}({\bm{\mathbf{{0}}}},\tilde{\lambda}_{b_{p}}^{-1}{\bf
R}^{-1})\mathcal{G}(a_{\lambda},b_{\lambda})\end{split}$ (9)

Extending this prior for simultaneous variable selection and clustering
requires only one modification:

$\displaystyle\begin{split}[\forall\in\\{1,\dots,P_{c}\\}]\ \
({\bm{\tilde{\mathbf{{b}}}}}_{p},\tilde{\lambda}_{b_{p}})&\sim\pi_{0}\delta_{{\bm{\mathbf{{0}}}}}+(1-\pi_{0})G^{*}\\\
G^{*}&\sim DP(G_{0},\alpha)\\\
G_{0}&=\mathcal{N}_{M}({\bm{\mathbf{{0}}}},\tilde{\lambda}_{b_{p}}^{-1}{\bf
R}^{-1})\mathcal{G}(a_{\lambda},b_{\lambda})\\\
[\pi_{0},(1-\pi_{0})]&\sim\text{Dirichlet}\left(\frac{\alpha_{0}}{2},\frac{\alpha_{0}}{2}\right)\end{split}$
(10)

In all computations in this paper we give $\alpha$ a non-informative Gamma
prior and update it using the method outlined in Escobar and West (1995). We
also fix $\alpha_{0}=2$ so that $\pi_{0}\sim\text{Unif}(0,1)$.

##  4 Computation

In this section we give the details of a Gibbs sampler for estimating the
parameters of our proposed models. The crux of the sampler is the update of
the class labels for each parameter, ${\bf c}$, and we extend the approach
used by Nott (2008) to the model in (7).

We first list the model hierarchy conditioned on the cluster indicator matrix,
${\bf C}$, and assume that it has $K$ columns. We also modify (7) and work
with the transposed equation, ${\bf Y}^{\sf T}={{\bm{\mathbf{{\Theta}}}}}{\bf
A}{\bf W}^{\sf T}+{{\bm{\mathbf{{\Theta}}}}}{\bf B}({\bf X}{\bf C})^{\sf
T}+{\bf E}^{\sf T}$. This, combined with the identity $\text{vec}({\bf A}{\bf
B}{\bf C}^{\sf T})=({\bf C}\otimes{\bf A})\text{vec}({\bf B})$, implies that
our model can be written as ${\bf y}={\bm{\tilde{\mathbf{{W}}}}}{\bf
a}+{\bm{\tilde{\mathbf{{X}}}}}{\bf b}+{\bf e}$ where ${\bf y}=\text{vec}({\bf
Y}^{\sf T})$, ${\bf e}=\text{vec}({\bf E}^{\sf T})$, ${\bf a}=\text{vec}({\bf
A})$, ${\bf b}=\text{vec}({\bf B})$, ${\bm{\tilde{\mathbf{{X}}}}}=({\bf X}{\bf
C})\otimes{{\bm{\mathbf{{\Theta}}}}}$ and ${\bm{\tilde{\mathbf{{W}}}}}={\bf
W}\otimes{{\bm{\mathbf{{\Theta}}}}}$. Therefore, conditioned on ${\bf C}$, our
model hierarchy is given by:

$\displaystyle\begin{split}{\bf
y}&\sim\mathcal{N}_{NT}({\bm{\tilde{\mathbf{{W}}}}}{\bf
a}+{\bm{\tilde{\mathbf{{X}}}}}{\bf b},\tau^{-1}{\bf I})\\\ {\bf
a}&\sim\mathcal{N}_{MP_{f}}({\bm{\mathbf{{0}}}},{{\bm{\mathbf{{\Lambda}}}}}_{a}^{-1}\otimes{\bf
R}^{-1})\\\ {\bf
b}&\sim\mathcal{N}_{MK}({\bm{\mathbf{{0}}}},{{\bm{\mathbf{{\Lambda}}}}}_{b}^{-1}\otimes{\bf
R}^{-1})\\\ [\forall p\in\\{1,\dots,P_{c}\\}]\ \
\lambda_{a_{p}}&\sim\mathcal{G}(a_{\lambda},b_{\lambda})\\\ [\forall
p\in\\{1,\dots,K\\}]\ \
\lambda_{b_{p}}&\sim\mathcal{G}(a_{\lambda},b_{\lambda})\\\
\tau&\sim\mathcal{G}(a_{\tau},b_{\tau})\end{split}$ (11)

where
${{\bm{\mathbf{{\Lambda}}}}}_{a}=\text{diag}(\lambda_{a_{1}},\dots,\lambda_{a_{P_{f}}})$
and
${{\bm{\mathbf{{\Lambda}}}}}_{b}=\text{diag}(\lambda_{b_{1}},\dots,\lambda_{b_{K}})$.
From (11) it is easy to see that the updates of the parameters are similar to
those for a normal linear model, with the posteriors of ${\bf a}$ and ${\bf
b}$ being multivariate normal, and the posteriors of the $\lambda$’s and
$\tau$ being Gamma. Therefore, the rest of this section will focus on
estimation of the cluster indicator matrix, first for the clustering only
Dirichlet Process prior (9) followed the variable selection and clustering
prior (10).

###  4.1 Dirichlet Process

Updating ${\bf c}$ is the most computationally intensive step of the Gibbs
sampler because its elements have to be updated sequentially. To update the
class label for a predictor, $c_{p}$, we need to calculate posterior
probabilities of the predictor belonging to each cluster. In the linear model
setting, Nott (2008) updates $c_{p}$ by integrating out the parameter vector,
which is important for the mixing of the Markov chains, but requires $K+1$
matrix inversions for each predictor. Here, $K$ is the number of unique
elements in ${\bf c}_{-p}$ and the additional inversion is for the proposal of
a new cluster. Hence, his algorithm is a modification of the collapsed Gibbs
sampler; Algorithm 3 in Neal (2000).

For our purposes, the inclusion of a prior on the smoothing parameters in (9)
complicates the integral required to collapse the over all cluster parameters.
Hence, we combine Algorithms 3 and 8 from Neal (2000) and collapse only over
the cluster means. Our update of the elements of ${\bf c}$ proceeds as
follows. For each $c_{i}$, $i\in\\{1,\dots,P_{c}\\}$:

  * •

Let $k^{-}$ be the number of distinct $c_{j}$ for $j\neq i$ and $h=k^{-}+1$.

  * •

If $c_{i}=c_{j}$ for some $j\neq i$ draw from
$\mathcal{G}(a_{\lambda},b_{\lambda})$ for $\lambda_{b_{h}}$. If $\forall
j\neq i$, $c_{i}\neq c_{j}$, set $\lambda_{b_{h}}=\tilde{\lambda}_{b_{i}}$,
where $\tilde{\lambda}_{b_{i}}$ is the smoothing parameter for the $i^{th}$
predictor.

  * •

For each $k\in\\{1,\dots,h\\}$

$\displaystyle P(c_{p}=k|{\bf c}_{-p},{\bm{\mathbf{{\lambda}}}}_{b},\cdot)$
$\displaystyle\propto f({\bf y}|{\bf
c}^{\prime},{\bm{\mathbf{{\lambda}}}}_{b},\cdot)P(c_{p}=k|{\bf
c}_{-p},\alpha)$ $\displaystyle=\int f({\bf y}|{\bf b},{\bf
c}^{\prime},{\bm{\mathbf{{\lambda}}}}_{b},{\bf X},{\bf W},{\bf a},\tau)f({\bf
b}|{\bm{\mathbf{{\lambda}}}}_{b},{\bf R})\ d{\bf b}\times P(c_{p}=k|{\bf
c}_{-p},\alpha)$

where

$\displaystyle\begin{split}&\int f({\bf y}|{\bf b},{\bf
c}^{\prime},{\bm{\mathbf{{\lambda}}}}_{b},{\bf X},{\bf W},{\bf a},\tau)f({\bf
b}|{\bm{\mathbf{{\lambda}}}}_{b},{\bf R})\ d{\bf b}=\\\ &\ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ \
(2\pi)^{-\frac{NT}{2}}\tau^{\frac{NT-M\tilde{K}}{2}}|{{\bm{\mathbf{{\Lambda}}}}}_{b}^{\prime}|^{\frac{M}{2}}|{\bf
R}|^{\frac{\tilde{K}}{2}}|{\bf
G}|^{-\frac{1}{2}}\exp\left\\{-\frac{\tau}{2}({{\bm{\hat{\mathbf{{e}}}}}}_{W}^{\sf
T}{{\bm{\hat{\mathbf{{e}}}}}}_{W}-{\bf g}^{\sf T}{\bf G}^{-1}{\bf
g})\right\\}\end{split}$ (12)

and ${\bf G}=[({\bf X}{\bf C}^{\prime})\otimes{{\bm{\mathbf{{\Theta}}}}}]^{\sf
T}[({\bf X}{\bf
C}^{\prime})\otimes{{\bm{\mathbf{{\Theta}}}}}]+\frac{1}{\tau}({{\bm{\mathbf{{\Lambda}}}}}_{b}^{\prime}\otimes{\bf
R})$, ${\bf g}=[({\bf X}{\bf
C}^{\prime})\otimes{{\bm{\mathbf{{\Theta}}}}}]^{\sf
T}{{\bm{\hat{\mathbf{{e}}}}}}_{W}$,
${{\bm{\hat{\mathbf{{e}}}}}}_{W}=\text{vec}({\bf Y}^{\sf
T}-{{\bm{\mathbf{{\Theta}}}}}{\bf A}{\bf W}^{\sf T})$, ${\bf c}^{\prime}$ is
the vector of class labels with $c_{p}=k$, ${\bf C}^{\prime}$ is the
corresponding one-hot encoded matrix, $\tilde{K}$ is the number of clusters in
${\bf c}^{\prime}$, ${\bf b}$ is a $M\tilde{K}\times 1$ vector,
${{\bm{\mathbf{{\Lambda}}}}}_{b}^{\prime}$ is a $\tilde{K}\times\tilde{K}$
dimensional diagonal matrix of smoothing parameters for clusters present in
${\bf c}^{\prime}$, and $P(c_{p}=k|{\bf c}_{-i},\alpha)$ can be calculated by
using the formulas in (2). It should be noted that while $h$ inverses have to
be computed for each predictor, they can be done in parallel to speed up the
algorithm.

###  4.2 Dirichlet Process + Point Mass

Computation for this model is similar to the model discussed in Section 4.1
because the prior in (9) is a clustering prior with the restriction that the
first cluster is zero. We can account for this by making a few changes to the
Gibbs sampler outlined above. To update ${\bf c}$, conditional prior
probabilities of the form in (2) need to be determined and the integral in
(12) has to be modified by dropping, from ${\bf X}$, the columns in the null
(zero) cluster. The update of ${\bm{\tilde{\mathbf{{B}}}}}$ also relies on
dropping the columns in the null cluster because their effects are restricted
to be zero; the rest of the update proceeds as before with the smaller
predictor matrix.

If we let $P_{0}$ denote the number of variables which are in the null cluster
and $P_{nz}$ denote the number of variables which are non-zero
$(P_{c}=P_{0}+P_{nz})$, we can calculate conditional prior probabilities of
the form in (2) by thinking of the prior in two levels: the first being a two
component finite mixture model and the second a Dirichlet Process. To begin,
we need to determine the prior probabilities of the predictor $p$ being equal
to zero, $P(c_{p}=0|{\bf c}_{-p},\alpha_{0})$, where $[\pi_{0},(1-\pi_{0})]$
is given a
$\text{Dirichlet}\left(\frac{\alpha_{0}}{2},\frac{\alpha_{0}}{2}\right)$
prior. This is given directly by Neal (2000) to be

$\displaystyle\pi^{*}_{0}=\pi(c_{p}=0|{\bf c}_{-p},\alpha_{0})$
$\displaystyle=\frac{n_{-p,0}+\alpha_{0}/2}{P_{c}-1+\alpha_{0}}$ (13)

Because only a subset of the predictors are non zero, the Dirichlet Process
prior can be thought of as active for only $P_{nz}$ predictors. Therefore, the
probabilities from (13) can be combined with the probabilities from (2) to get
the conditional prior for cluster membership:

$\displaystyle\begin{split}\pi^{*}_{0}=P(c_{p}=0|{\bf
c}_{-p},\alpha_{0})&=\frac{n_{-p,0}+\alpha_{0}/2}{P_{c}-1+\alpha_{0}}\\\
\text{For $k\geq 1$, }P(c_{p}=k|{\bf
c}_{-p},\alpha)&=(1-\pi_{0}^{*})\left(\frac{n_{-p,k}}{P_{nz}-1+\alpha}\right)\\\
P(c_{p}\neq 0\text{ and }c_{p}\neq c_{j}\text{ for all }p\neq j|{\bf
c}_{-p},\alpha)&=(1-\pi_{0}^{*})\left(\frac{\alpha}{P_{nz}-1+\alpha}\right)\end{split}$
(14)

##  5 Simulations

| N = 30 | N = 60 | N = 120 | N = 240
---|---|---|---|---
Design 1 |  |  |  |
FOSR | 3.02 (0.16) | 1.19 (0.08) | 0.64 (0.04) | 0.3 (0.02)
FOSR-PM | 2.12 (0.16) | 1.2 (0.08) | 0.71 (0.05) | 0.29 (0.02)
FOSR-DP | 0.94 (0.05) | 0.4 (0.03) | 0.16 (0.01) | 0.07 (0.01)
FOSR-DPPM | 1.12 (0.06) | 0.43 (0.03) | 0.17 (0.01) | 0.06 (0.01)
Design 2 |  |  |  |
FOSR | 4.61 (0.22) | 2.8 (0.15) | 1.38 (0.08) | 0.73 (0.04)
FOSR-PM | 4.1 (0.22) | 3.49 (0.2) | 1.94 (0.11) | 0.99 (0.07)
FOSR-DP | 1.49 (0.09) | 0.89 (0.06) | 0.37 (0.03) | 0.16 (0.01)
FOSR-DPPM | 2.1 (0.12) | 1.13 (0.07) | 0.48 (0.03) | 0.19 (0.02)
Design 3 |  |  |  |
FOSR | 2.75 (0.1) | 1.18 (0.03) | 0.64 (0.02) | 0.31 (0.01)
FOSR-PM | 2.49 (0.08) | 1.53 (0.04) | 0.81 (0.04) | 0.36 (0.01)
FOSR-DP | 1.85 (0.05) | 1.18 (0.03) | 0.75 (0.03) | 0.41 (0.01)
FOSR-DPPM | 2.16 (0.06) | 1.27 (0.04) | 0.78 (0.03) | 0.42 (0.02)
Design 4 |  |  |  |
FOSR | 1.73 (0.07) | 0.7 (0.03) | 0.31 (0.01) | 0.16 (0.01)
FOSR-PM | 1.09 (0.03) | 0.67 (0.03) | 0.29 (0.01) | 0.14 (0.01)
FOSR-DP | 0.95 (0.03) | 0.62 (0.02) | 0.3 (0.01) | 0.17 (0.01)
FOSR-DPPM | 1.16 (0.04) | 0.62 (0.02) | 0.3 (0.01) | 0.17 (0.01)
Table 1: Mean pointwise MSE for 100 datasets for each model. Standard errors
(in parentheses) are estimated using the bootstrap with 100 repetitions.

(a) RAND Index | N = 30 | N = 60 | N = 120 | N = 240
---|---|---|---|---
Design 1 |  |  |  |
FOSR-DP | 0.69 (0.01) | 0.79 (0.01) | 0.92 (0.01) | 0.96 (0.01)
FOSR-DPPM | 0.65 (0.01) | 0.8 (0.01) | 0.92 (0.01) | 0.97 (0.01)
Design 2 |  |  |  |
FOSR-DP | 0.69 (0.01) | 0.78 (0.01) | 0.87 (0.01) | 0.93 (0.01)
FOSR-DPPM | 0.63 (0.01) | 0.76 (0.01) | 0.84 (0.01) | 0.92 (0.01)

(b) Adjusted RAND Index | N = 30 | N = 60 | N = 120 | N = 240
---|---|---|---|---
Design 1 |  |  |  |
FOSR-DP | 0.41 (0.02) | 0.59 (0.02) | 0.83 (0.02) | 0.93 (0.01)
FOSR-DPPM | 0.47 (0.03) | 0.64 (0.02) | 0.84 (0.02) | 0.94 (0.01)
Design 2 |  |  |  |
FOSR-DP | 0.43 (0.02) | 0.57 (0.03) | 0.78 (0.02) | 0.87 (0.02)
FOSR-DPPM | 0.22 (0.02) | 0.45 (0.03) | 0.64 (0.02) | 0.87 (0.02)

Table 2: Mean RAND and Adjusted RAND index for 100 datasets is reported for
each model. Standard errors (in parentheses) are estimated using the bootstrap
with 100 repetitions.

In this section we evaluate the efficacy of our models using a simulation
study. Our simulations consist of four designs and compare the performance of
the proposed priors with two priors from the literature; the smoothing prior
proposed by Goldsmith and Kitago (2016) and a modification of the variable
selection prior proposed by Goldsmith and Schwartz (2017). For clarity, the
priors and the corresponding names used in our tables are listed below:

  * •

FOSR: Goldsmith and Kitago (2016)

$\displaystyle({\bm{\tilde{\mathbf{{b}}}}}_{p},\tilde{\lambda}_{b_{p}})$
$\displaystyle\sim\mathcal{N}_{M}({\bm{\mathbf{{0}}}},\tilde{\lambda}_{b_{p}}^{-1}{\bf
R}^{-1})\mathcal{G}(a_{\lambda},b_{\lambda})$

  * •

FOSR-PM: Modification of Goldsmith and Schwartz (2017) to include smoothing

$\displaystyle({\bm{\tilde{\mathbf{{b}}}}}_{p},\tilde{\lambda}_{b_{p}})$
$\displaystyle\sim\pi_{0}\delta_{{\bm{\mathbf{{0}}}}}+(1-\pi_{0})G\text{ with
}G=\mathcal{N}_{M}({\bm{\mathbf{{0}}}},\tilde{\lambda}_{b_{p}}^{-1}{\bf
R}^{-1})\ \mathcal{G}(a_{\lambda},b_{\lambda})$

  * •

FOSR-DP: Clustering prior given in (9)

  * •

FOSR-DPPM: Variable selection and clustering prior given in (10)

###  5.1 Design

All of the methods considered are fit under the assumption of iid errors, even
though our simulations use an exponential covariance matrix with a nugget
term. Using the notation outlined in Equation (7), we set $T=15$, $P_{f}=5$,
and $P_{c}=15$ for all simulations, and use three Fourier basis functions to
generate the coefficient functions. Additionally, we generate correlated
values for ${\bf X}$ with
$Cov(x_{i,p},x_{i,p^{\prime}})=0.75^{|p-p^{\prime}|}$, while
$w_{ij}\overset{iid}{\sim}\mathcal{N}(0,1)$. Finally, we simulate each dataset
with
${{\bm{\mathbf{{\Sigma}}}}}={{\bm{\mathbf{{\Sigma}}}}}^{\prime}+\sigma^{2}{\bf
I}$ where ${{\bm{\mathbf{{\Sigma}}}}}^{\prime}_{i,j}=\exp(-10*(i-j)^{2})$ and
$\sigma^{2}$ was chosen to set the signal-to-noise ratio in the data to 1\.
The major variation in our simulation designs is the cluster membership
pattern of ${{\bm{\mathbf{{\beta}}}}}$; it is varied as follows:

  * •

Design 1: ${{\bm{\mathbf{{\beta}}}}}$ has three clusters with ${\bf
c}=({\bm{\mathbf{{1}}}}_{7}^{\sf T},{\bm{\mathbf{{2}}}}_{4}^{\sf
T},{\bm{\mathbf{{3}}}}_{4}^{\sf T})^{\sf T}$. The coefficients of cluster 1
are set to zero.

  * •

Design 2: ${{\bm{\mathbf{{\beta}}}}}$ has three clusters with the same pattern
as Design 1. However, all of the effects are non-zero.

  * •

Design 3: ${{\bm{\mathbf{{\beta}}}}}$ has no clusters and none of the
coefficients are zero.

  * •

Design 4: Similar to Design 3 but the first seven elements of
${{\bm{\mathbf{{\beta}}}}}$ are set to zero.

###  5.2 Results

Each model was estimated using the Gibbs sampler described in Section 4, run
for 5000 iterations with 2500 used as burn-in. Additionally, we used eight
b-spline basis functions to estimate each coefficient function. We evaluate
our methods using two criteria: parameter estimation using the pointwise mean
squared error (MSE) of the functional effects, and clustering ability using
the RAND and adjusted RAND indices.

Table 1 shows the pointwise MSE of our proposed methods in comparison to the
other priors considered. It can be seen that when cluster structure is present
in the coefficients, as in Designs 1 and 2, the FOSR-DP and FOSR-DPPM methods
outperform FOSR and FOSR-PM handily. This is to be expected since our priors
are designed to exploit this structure. What is interesting however, is the
out-performance of our methods in the small sample setting for Designs 3 and
4. The data in Design 3 is simulated to give FOSR an advantage, however both
FOSR-DP and FOSR-DPPM outperform FOSR when $N=30$. Since the total number of
predictors in the simulations is fixed at 20, and none are set to zero in
Design 3, the high dimensionality of problem seems to lend itself well to the
use of a dimension reduction method which borrows information from nearby
covariates. We see something similar in Design 4, where FOSR-DP outperforms
the FOSR-PM prior.

We also evaluated the clustering capacity of our methods using the RAND and
Adjusted RAND indices, with results in Table 2. As can be seen, our methods
are able to learn cluster structure in the coefficients despite miss-
specification of the covariance matrix and a low signal-to-noise ratio.

##  6 Application

We demonstrate the use of our methods by analyzing age-specific fertility
rates from the United Nations Gender Information (UNGEN) database. We use data
from surveys conducted between 2000-2005, estimates of which are plotted in
Figure 1. Due to the nature of the survey, fertility curves are only observed
at 7 points, where each time point corresponds to an age group for the
participants in question. Each curve seems to be a relatively smooth function
of time, making a function-on-scalar analysis of this data appropriate. A
similar data set is explored by Kowal (2018) but our analysis considers a
larger population of countries and a greater number of demographic and
socioeconomic predictors.

Figure 1: Age-specific fertility rates for 92 countries from the United
Nations Gender Information Database

The response information from the UNGEN database is combined with fifteen
country level covariates available on Gapminder. Each covariate included in
the final model is a complete case average of Gapminder data from 2000-2005.
This allows us to see how average values of demographic and socioeconomic
variables during the years of the UN survey affect the fertility curve. The
variables we consider include some potentially causal factors such as age at
first marriage, something which should increase fertility in the early parts
of the curve, along with seemingly inconsequential variables such as the
proportion of dollar billionaires in a country and the amount of alcohol its
inhabitants consume. Unfortunately, many countries had no covariate
information for the years of the survey and were dropped from the analysis for
simplicity. The final data set consists of 15 covariates for 92 countries
around the world.

Figure 2 shows many highly correlated predictors among consideration. For
example, Under-5 Mortality is highly positively correlated with Maternal
Deaths (correlation coefficient of 0.85), but is highly negatively correlated
with Contraception Prevalence (correlation of -0.85). Contraception Prevalence
on the other hand is positively correlated with many socioeconomic factors
such as life expectancy and the number of births attended by trained birth
staff.

Figure 2: Plot of the correlation matrix of predictors used to fit the FOSR-DP
and FOSR-DPPM models

We fit the FOSR-DP and FOSR-DPPM models on the data set. Each method was fit
using 10,000 iterations with 5,000 used for burn-in. Figure 3 show the cluster
dendogram for the analyses. Viewing both dendograms we see that covariates
with high positive correlation tend to be clustered together. Additionally,
the structure of the tree seems similar across the two models, with small
differences in the order that they are joined in a particular node. Finally,
the variable that we thought a-priori should most clearly impact the age-
specific fertility curve, age at first marriage, is a standalone covariate in
both dendograms. Table 3 shows the percent of iterations after burn-in that
each covariate was set to zero by the FOSR-DPPM model; of the fifteen
covariates considered, only six are estimated to have non-zero effects if a 5%
cutoff is used.

Figure 3: Dendograms from the DP and DPPM clust models | Percent Zero
---|---
Contraception Prevalence (15-49) | $0.000$
Maternity Related Deaths per 1K Women | $0.000$
Births Atten. by Exp. Staff (% of Total) | $0.000$
Age at First Marriage | $0.006$
Male to Female Sex Ratio (15-49) | $0.007$
U5 Mortality | $0.012$
Cervical Cancer Deaths (Per 100K) | $0.148$
Health Expenditures (% of GDP) | $0.291$
Female BMI | $0.556$
Female Labor Force Participation | $0.579$
GDP Per Capita | $0.814$
Mean Yrs. of School (Women % Men) (25-34) | $0.861$
Life Expectancy | $0.901$
Dollar Billionares (per 1M) | $0.948$
Alcohol Consumption per Adult | $0.959$
Table 3: Shows the proportion of iterations that the variables listed were set
to zero by the FOSR-DPPM model

Finally, the estimated effects, on a standardized scale, of the non-zero
coefficients are shown in Figure 4. From this we can see that the age at first
marriage has a time varying effect; a decrease of one standard deviation of
age increases the number of births in the 15-19 and 20-24 age groups, but the
effect disappears as age increases. Additionally, U5 mortality and Maternal
Deaths per 1-K are highly correlated variables with clustered coefficient
curves. Their effect is positive throughout, but decreases as age increases.
Finally, an interesting finding in our estimates is that contraception
prevalence primarily impacts fertility in the middle of the curve; it is not
impactful in the 15-19 age group.

Figure 4: Plots showing the estimated coefficient curves from the FOSR-DP and
FOSR-DPPM models. The solid lines are the posterior means of 5000 samples from
the Gibbs sampler, while the dotted lines are the 2.5 and 97.5 percentile
respectively

##  7 Conclusion

In this paper we extend methodology for addressing multicollinearity in the
predictors to function-on-scalar regression. We contribute to the literature
by proposing a prior which simultaneously selects, clusters, and smooths the
coefficient functions using a Bayesian approach, while the current Bayesian
literature only focuses independently on selection or smoothing. Our model
allows coefficient estimates to remain flexible while controlling overfitting,
and performs dimension reduction by summing columns in the covariate matrix
which have the same effect on the response. We also develop a Gibbs sampler
which converges to the posterior distribution quickly in practice by combining
multiple sampling algorithms developed in the clustering literature.

Our work can be extended in many ways. First, the Gibbs sampler proposed in
Section 4 is not scalable to problems with large $N$ or $P$, because
collapsing the sampler over the basis function coefficients forces us to solve
many linear systems in each iteration. Faster methods utilizing other
algorithms for Dirichlet Process priors, with a focus on Variational Bayes
approaches, have the potential to make our methodology applicable to larger
data sets. Second, our methods assume that the coefficient curves are globally
clustered, which can be a strong assumption in practice. Finally, if point
estimation is the only goal, methods such as OSCAR, the Clustered Lasso, and
PACS can be extended to the function-on-scalar regression framework, possibly
providing another computationally efficient alternative to our approach.

## References

  * Abraham et al. (2003) Christophe Abraham, Pierre-André Cornillon, ERIC Matzner-Løber, and Nicolas Molinari.  Unsupervised curve clustering using b-splines.  _Scandinavian journal of statistics_ , 30(3):581–595, 2003.
  * Antoniak (1974) Charles E. Antoniak.  Mixtures of dirichlet processes with applications to bayesian nonparametric problems.  _Ann. Statist._ , 2(6):1152–1174, 11 1974.
  * Blackwell et al. (1973) David Blackwell, James B MacQueen, et al.  Ferguson distributions via pólya urn schemes.  _The annals of statistics_ , 1(2):353–355, 1973\.
  * Bondell and Reich (2008) Howard D Bondell and Brian J Reich.  Simultaneous regression shrinkage, variable selection, and supervised clustering of predictors with oscar.  _Biometrics_ , 64(1):115–123, 2008.
  * Chen et al. (2016) Yakuan Chen, Jeff Goldsmith, and R Todd Ogden.  Variable selection in function-on-scalar regression.  _Stat_ , 5(1):88–101, 2016.
  * Curtis and Ghosh (2011) S McKay Curtis and Sujit K Ghosh.  A bayesian approach to multicollinearity and the simultaneous selection and clustering of predictors in linear regression.  _Journal of Statistical Theory and Practice_ , 5(4):715–735, 2011.
  * Dunson et al. (2008) David B Dunson, Amy H Herring, and Stephanie M Engel.  Bayesian selection and clustering of polymorphisms in functionally related genes.  _Journal of the American Statistical Association_ , 103(482):534–546, 2008.
  * Eilers and Marx (1996) Paul HC Eilers and Brian D Marx.  Flexible smoothing with b-splines and penalties.  _Statistical science_ , pages 89–102, 1996.
  * Escobar and West (1995) Michael D Escobar and Mike West.  Bayesian density estimation and inference using mixtures.  _Journal of the american statistical association_ , 90(430):577–588, 1995.
  * Fan and Reimherr (2017) Zhaohu Fan and Matthew Reimherr.  High-dimensional adaptive function-on-scalar regression.  _Econometrics and statistics_ , 1:167–183, 2017.
  * Ferguson (1973) Thomas S. Ferguson.  A bayesian analysis of some nonparametric problems.  _Ann. Statist._ , 1(2):209–230, 03 1973.
  * Goldsmith and Kitago (2016) Jeff Goldsmith and Tomoko Kitago.  Assessing systematic effects of stroke on motor control by using hierarchical function-on-scalar regression.  _Journal of the Royal Statistical Society: Series C (Applied Statistics)_ , 65(2):215–236, 2016.
  * Goldsmith and Schwartz (2017) Jeff Goldsmith and Joseph E Schwartz.  Variable selection in the functional linear concurrent model.  _Statistics in medicine_ , 36(14):2237–2250, 2017\.
  * Jacques and Preda (2014) Julien Jacques and Cristian Preda.  Functional data clustering: a survey.  _Advances in Data Analysis and Classification_ , 8(3):231–255, 2014.
  * Kim et al. (2006) Sinae Kim, Mahlet G Tadesse, and Marina Vannucci.  Variable selection in clustering via dirichlet process mixture models.  _Biometrika_ , 93(4):877–893, 2006.
  * Kowal (2018) Daniel R Kowal.  Dynamic function-on-scalars regression.  _arXiv preprint arXiv:1806.01460_ , 2018.
  * Kowal and Bourgeois (2018) Daniel R Kowal and Daniel C Bourgeois.  Bayesian function-on-scalars regression for high dimensional data.  _arXiv preprint arXiv:1808.06689_ , 2018.
  * MacEachern and Müller (1998) Steven N MacEachern and Peter Müller.  Estimating mixture of dirichlet process models.  _Journal of Computational and Graphical Statistics_ , 7(2):223–238, 1998.
  * Neal (2000) Radford M Neal.  Markov chain sampling methods for dirichlet process mixture models.  _Journal of computational and graphical statistics_ , 9(2):249–265, 2000.
  * Nott (2008) David J Nott.  Predictive performance of dirichlet process shrinkage methods in linear regression.  _Computational Statistics & Data Analysis_, 52(7):3658–3669, 2008.
  * Parodi and Reimherr (2018) Alice Parodi and Matthew Reimherr.  Simultaneous variable selection and smoothing for high-dimensional function-on-scalar regression.  _Electron. J. Statist._ , 12(2):4602–4639, 2018\.  doi: 10.1214/18-EJS1509.  URL https://doi.org/10.1214/18-EJS1509.
  * Ray and Mallick (2006) Shubhankar Ray and Bani Mallick.  Functional clustering by bayesian wavelet methods.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 68(2):305–332, 2006.
  * Sethuraman (1994) Jayaram Sethuraman.  A constructive definition of dirichlet priors.  _Statistica sinica_ , pages 639–650, 1994.
  * Sharma et al. (2013) Dhruv B. Sharma, Howard D. Bondell, and Hao Helen Zhang.  Consistent group identification and variable selection in regression with correlated predictors.  _Journal of Computational and Graphical Statistics_ , 22(2):319–340, 2013.
  * She (2010) Yiyuan She.  Sparse regression with exact clustering.  _Electronic Journal of Statistics_ , 4:1055–1096, 2010\.
  * Tadesse et al. (2005) Mahlet G Tadesse, Naijun Sha, and Marina Vannucci.  Bayesian variable selection in clustering high-dimensional data.  _Journal of the American Statistical Association_ , 100(470):602–617, 2005.
  * Tibshirani and Taylor (2011) Ryan J. Tibshirani and Jonathan Taylor.  The solution path of the generalized lasso.  _Ann. Statist._ , 39(3):1335–1371, 06 2011.  doi: 10.1214/11-AOS878.  URL https://doi.org/10.1214/11-AOS878.
  * Wang et al. (2007) Lifeng Wang, Guang Chen, and Hongzhe Li.  Group scad regression analysis for microarray time course gene expression data.  _Bioinformatics_ , 23(12):1486–1494, 2007.
  * Yuan and Lin (2006) Ming Yuan and Yi Lin.  Model selection and estimation in regression with grouped variables.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 68(1):49–67, 2006.
  * Zhang and Telesca (2014) Yafeng Zhang and Donatello Telesca.  Joint clustering and registration of functional data.  _arXiv preprint arXiv:1403.7134_ , 2014.