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:
 
Y= WT+ XT+ E
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, \({}_{T \times P_f}\) , \({\bf X}_{N \times P_c}\)  is a matrix of predictors whose effects, \({}_{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(, \tau ^{-1}{\bf I})\) . We expand each
coefficient function as a linear combination of \(M\)  b-spline basis functions,
\(\lbrace \theta _1, \dots , \theta _M \rbrace \) , such that \({}_{T \times P_f} = {{{\mathbf {\Theta }}}}_{T\times M} {\bf A}_{M \times P_f}\)  and \({}_{T \times P_c} = {{{\mathbf {\Theta }}}}_{T \times M}{{\tilde{\mathbf {}{B}}}_{M \times P_c}. This reduces the functional estimation problem toestimating the elements of {\bf A} and {{\tilde{\mathbf {}{B}}}.}In practice, care must be taken in choosing M, balancing the need for afunction to remain flexible with the risk of overfitting. We address thisissue by expanding the functional effects using more basis functions than wethink are necessary and use penalization to induce smoothness in the coefficientfunction. We follow \cite {goldsmith2016assessing} and use a full-rank penaltymatrix {\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-splinepenalty \cite {eilers1996flexible}, and\eta = 0.001. Consequently, the prior for the columns of {\bf A}, {\bf a}_p, is\forall p \in \lbrace 1, \dots , P_f \rbrace \begin{align} {\bf a}_p & \sim \mathcal {N}_{M}(, \lambda _{a_p}^{-1}{\bf R}^{-1})\end{align}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.}\subsection {Clustering Functional Effects}\) Many approaches exist for clustering functional data and we refer the reader to
[1]} for a review. Because we expand each coefficient
function using b-splines, our clustering approach is similar to the work
initially proposed by [2]}; we find clusters in
\({}\)  by clustering the coefficients of the basis expansion. If we assume the
columns of \({}\)  come from \(K\)  clusters, we can expand it as \({}={{{\mathbf {\Theta }}}}{{\tilde{\mathbf {}{B}}}= {{{\mathbf {\Theta }}}}{\bf B}{\bf C}^{\sf T} where {{{\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 classmembership information. Consequently, {\bf X}{}^{\sf T}= {\bf X}{{\tilde{\mathbf {}{B}}}^{\sf T}{{{\mathbf {\Theta }}}}^{\sf T}= {\bf X}{\bf C}{\bf B}^{\sf T}{{{\mathbf {\Theta }}}}^{\sf T}, and {\bf C} has the effect of summingcolumns of {\bf X} with the same effect on the response.}Because we are also interested in estimating smooth coefficient functions, wegive each predictor its own smoothing parameter, \tilde{\lambda }_{b_p}.Therefore, our task is to estimate a mean and smoothing parameter for each ofour clusters. We can easily do this by modifying our base distribution to be amultivariate normal-gamma \cite {ray2006functional, zhang2014joint}. Consequently, for clustering functional effects, wepropose the prior:\begin{align} \begin{split}[\forall \in \lbrace 1, \dots , P_c \rbrace ] \ \ ({{\tilde{\mathbf {}{b}}}_p, \tilde{\lambda }_{b_p}) & \sim G \\G & \sim DP(G_0, \alpha ) \\G_0 & = \mathcal {N}_{M}(, \tilde{\lambda }_{b_p}^{-1}{\bf R}^{-1}) \mathcal {G}(a_{\lambda }, b_{\lambda })}\end{split}Extending this prior for simultaneous variable selection and clustering requiresonly one modification:\begin{align} \begin{split}[\forall \in \lbrace 1, \dots , P_c \rbrace ] \ \ ({{\tilde{\mathbf {}{b}}}_p, \tilde{\lambda }_{b_p}) & \sim \pi _0 \delta _{} + (1 - \pi _0) G^* \\G^* & \sim DP(G_0, \alpha ) \\G_0 & = \mathcal {N}_{M}(, \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}In all computations in this paper we give \alpha  a non-informative Gamma priorand update it using the method outlined in \cite {escobar1995bayesian}. We alsofix \alpha _0 = 2 so that \pi _0 \sim \text{Unif}(0, 1).\end{align}\section {Computation} \end{align}In this section we give the details of a Gibbs sampler for estimating theparameters of our proposed models. The crux of the sampler is the update of theclass labels for each parameter, {\bf c}, and we extend the approach used by\cite {nott2008predictive} to the model in (\ref {eqn:dense_model}).}We first list the model hierarchy conditioned on the cluster indicatormatrix, \)C\(, and assume that it has \) K\( columns. We also modify(\ref {eqn:dense_model}) and work with the transposed equation, \)YT= AWT+ B(XC)T+ ET\(. This, combined with the identity\) vec(ABCT) = (CA) vec(B)\(, implies that ourmodel can be written as \)y= Wa+ Xb+ e\( where \)y=
vec(YT)\(, \)e= vec(ET)\(, \)a= vec(A)\(, \)b=
vec(B)\(, \) X= (XC) \( and \) W= W\(. Therefore, conditioned on \)C\(, our model hierarchy is given by:\begin{align} \begin{split}{\bf y}& \sim \mathcal {N}_{NT}({{\tilde{\mathbf {}{W}}}{\bf a}+ {{\tilde{\mathbf {}{X}}}{\bf b}, \tau ^{-1}{\bf I}) \\{\bf a}& \sim \mathcal {N}_{MP_f}(, {{{\mathbf {\Lambda }}}}_a^{-1}\otimes {\bf R}^{-1}) \\{\bf b}& \sim \mathcal {N}_{MK}(, {{{\mathbf {\Lambda }}}}_b^{-1}\otimes {\bf R}^{-1}) \\[\forall p \in \lbrace 1, \dots , P_c \rbrace ] \ \ \lambda _{a_p} & \sim \mathcal {G}(a_{\lambda }, b_{\lambda }) \\[\forall p \in \lbrace 1, \dots , K \rbrace ] \ \ \lambda _{b_p} & \sim \mathcal {G}(a_{\lambda }, b_{\lambda }) \\\tau & \sim \mathcal {G}(a_{\tau }, b_{\tau })}}where {{{\mathbf {\Lambda }}}}_a = \text{diag}(\lambda _{a_1}, \dots , \lambda _{a_{P_f}}) and{{{\mathbf {\Lambda }}}}_b = \text{diag}(\lambda _{b_1}, \dots , \lambda _{b_K}). From (\ref {eqn:cond_hier}) it is easy to see that the updates of theparameters are similar to those for a normal linear model, with the posteriorsof {\bf a} and {\bf b} being multivariate normal, and the posteriors of the\lambda ^{\prime }s and \tau  being Gamma. Therefore, the rest of this section willfocus on estimation of the cluster indicator matrix, first for the clusteringonly Dirichlet Process prior (\ref {eqn:prior_dp}) followed the variableselection and clustering prior (\ref {eqn:prior_dppm}).\end{split}\subsection {Dirichlet Process} \end{align}Updating \)c\( is the most computationally intensive step of the Gibbssampler because its elements have to be updated sequentially. To update theclass label for a predictor, \) cp\(, we need to calculate posterior probabilitiesof the predictor belonging to each cluster. In the linear model setting,\cite {nott2008predictive} updates \) cp\( by integrating out the parametervector, 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 uniqueelements in \)c-p\( and the additional inversion is for the proposal of a newcluster. Hence, his algorithm is a modification of the collapsed Gibbs sampler;Algorithm 3 in \cite {neal2000markov}.\) For our purposes, the inclusion of a prior on the smoothing parameters in
() complicates the integral required to collapse the over all
cluster parameters. Hence, we combine Algorithms 3 and 8 from
[3]} and collapse only over the cluster means. Our update of
the elements of \({\bf c}\)  proceeds as follows. For each \(c_i\) , \(i \in \lbrace 1, \dots , P_c \rbrace \) :
