Covariate-Guided Bayesian Mixture of Spline Experts for the Analysis of Multivariate Time Series

We propose a tensor-product mixture of splines model for multivariate time series. For each subject $i=1,\ldots,N$, let $\boldsymbol{y}_{i}=(\boldsymbol{y}_{i1}^{\prime},\ldots,\boldsymbol{y}_{ik}^{\prime},\ldots,\boldsymbol{y}_{iK}^{\prime})^{\prime}$ be the $nK$-vector corresponding to the $K$-dimensional time series for $k=1,\ldots,K$, where $\boldsymbol{y}_{ik}=\big[y_{ik}(t_{1}),\ldots,y_{ik}(t_{j}),\ldots,y_{ik}(t_{n})\big]^{\prime}$ contains the trajectory of measurements on the $k$th entry of the time series evaluated over a grid of $n$ time points for $j=1,\ldots,n$, and $\boldsymbol{\epsilon}_{i}=(\boldsymbol{\epsilon}_{i1}^{\prime},\ldots,\boldsymbol{\epsilon}_{iK}^{\prime})^{\prime}$ is the $nK$-vector of errors. Following the model representation of Krafty and others (2017), the tensor-product model for the $K$-dimensional multivariate time series, conditional on component $g$, $g=1,\ldots,G$, can be written as:

where $\{z_{ig}\}_{g=1}^{G}$ are latent indicators as described in Section 3.3, $\boldsymbol{\alpha}_{g}=(\boldsymbol{\alpha}_{g1}^{\prime},\ldots,\boldsymbol{\alpha}_{gK}^{\prime})^{\prime}$ is a $2K$-vector of intercepts and slopes, $\boldsymbol{\beta}_{g}=(\boldsymbol{\beta}_{g1}^{\prime},\ldots,\boldsymbol{\beta}_{gK}^{\prime})^{\prime}$ is a $mK$-vector of basis function coefficients as described in Section 4.1, $\boldsymbol{I}_{K}$ is a $K\times K$ identity matrix and $\otimes$ denotes a tensor product. The matrix $\boldsymbol{X}$ is given by $\boldsymbol{X}=\begin{pmatrix}1&1&\ldots&1\\ t_{1}&t_{2}&\ldots&t_{n}\end{pmatrix}^{\prime}$ and the $m$ columns of the matrix $\boldsymbol{W}$ are smoothing splines basis functions as described in Section 4.1. We assume the error vector $\boldsymbol{\epsilon}_{i}$ follows a $\mbox{MVN}(\boldsymbol{0},\boldsymbol{\Psi}_{g}\otimes\boldsymbol{U})$ distribution, where $\boldsymbol{U}=\boldsymbol{I}_{n}$ is the $n\times n$ identity matrix, and $\boldsymbol{\Psi}_{g}=\mbox{diag}(\boldsymbol{\sigma}_{g}^{2})$ is a $K\times K$ diagonal matrix with the error variances $\boldsymbol{\sigma}_{g}^{2}=(\sigma_{g1}^{2},\ldots,\sigma_{gK}^{2})^{\prime}$. We assume each subject has a common grid of time points across all $K$ entries, such that $\boldsymbol{X}$ and $\boldsymbol{W}$ are common to all subjects, although our proposed method can be generalized to the case where subjects are observed at different grids of time points. In addition, we assume $\boldsymbol{E}(\boldsymbol{y}_{ik},\boldsymbol{y}_{ih})=\boldsymbol{0}_{n\times n}$ for $k\neq h$.

To simplify notation, we let $\boldsymbol{S}=[\boldsymbol{X}\ \boldsymbol{W}]$ and $\boldsymbol{\theta}_{g}=(\boldsymbol{\alpha}_{g1}^{\prime},\boldsymbol{\beta}_{g1}^{\prime},\ldots,\boldsymbol{\alpha}_{gK}^{\prime},\boldsymbol{\beta}_{gK}^{\prime})^{\prime}$. Equation (1) can then be rewritten as:

The mixture-of-experts model (Jacobs and others, 1991) is applied to form a covariate-guided structure for our proposed model, where the mixing weights are multinomial logits that are functions of selected covariates. As in Sun and others (2007), the mixing weights are expressed as

where $\boldsymbol{V}_{i}=(1,V_{i1},\cdots,V_{iP})^{\prime}$ is a vector of length $(P+1)$ containing values of $P$ covariates for subject $i$, and $\boldsymbol{\delta}_{g}=(\delta_{g0},\delta_{g1},\cdots,\delta_{gP})^{\prime}$ is the corresponding coefficient vector. For identifiability, we set $\boldsymbol{\delta}_{G}=\boldsymbol{0}$. Equation (3) differs slightly from the weights in the traditional mixture of experts model in that it includes a random term $\zeta_{ig}$ for each subject. This term accounts for unmeasured factors beyond the observed covariates, and enhances model performance and inference of the mixing weights.

To account for heterogeneity across subjects, we assume that the $k$th entry of the multivariate time series, $\boldsymbol{y}_{ik}$, comes from a mixture model with $G$ components, i.e.,

where $f_{gk}(\boldsymbol{y}_{ik}\mid\boldsymbol{\mu}_{gk},\sigma_{gk}^{2}\boldsymbol{I}_{n})$ is the probability density function of the multivariate normal distribution with mean vector $\boldsymbol{\mu}_{gk}=\boldsymbol{X}\boldsymbol{\alpha}_{gk}+\boldsymbol{W}\boldsymbol{\beta}_{gk}$ and covariance matrix $\sigma_{gk}^{2}\boldsymbol{I}_{n}$ for the $g$th component and the $k$th entry. The $\pi_{ig}$ are mixing weights that depend on covariates as described in Section 3.2.

As is common in mixture models, augmenting the likelihood with latent variables indicating the component from which a time series originates simplifies the computation greatly (Dempster and others, 1977). In particular, let $z_{ig}=1$ if the $i$th multivariate time series belongs to the $g$th component and $z_{ig}=0$, otherwise. Let $\boldsymbol{y}=(\boldsymbol{y}_{1},\ldots,\boldsymbol{y}_{N})^{\prime}$ be all observed multivariate time series and $\boldsymbol{\Theta}_{gk}$ be the aggregation of all parameters for component $g$ and entry $k$. The parameter vector for all components and all entries is then denoted by $\boldsymbol{\Theta}=(\boldsymbol{\Theta}_{11}^{\prime},\ldots,\boldsymbol{\Theta}_{GK}^{\prime})^{\prime}$. The augmented likelihood of all $N$ multivariate time series is given by

where $f_{gk}(\boldsymbol{y}_{ik}\mid\boldsymbol{\Theta}_{gk})$ is the probability density function as appeared in the (4). From Bayes’ rule, the distribution of the latent indicators $z_{ig}$ is given by

The conditional expectation of a mixture component in model (4) is given by $E(\boldsymbol{y}_{ik}\mid z_{ig}=1)=\boldsymbol{X}\boldsymbol{\alpha}_{gk}+\boldsymbol{W}\boldsymbol{\beta}_{gk}$. We place a smoothing spline prior on $\boldsymbol{\beta}_{gk}$ and let $\boldsymbol{\mathcal{H}}_{gk}=\boldsymbol{W}\boldsymbol{\beta}_{gk}$, where $\boldsymbol{\mathcal{H}}_{gk}=\big[\mathcal{H}_{gk}(t_{1}),\ldots,\mathcal{H}_{gk}(t_{n})\big]^{\prime}$ is a zero-mean Gaussian process with variance covariance matrix $\tau^{2}_{gk}\boldsymbol{\Phi}$ (Wahba, 1980; Wood and others, 2002), such that $\text{cov}\big[\mathcal{H}_{gk}(t_{r}),\mathcal{H}_{gk}(t_{h})\big]=\tau_{gk}^{2}\phi_{rh}$, $\tau_{gk}^{2}$ is a smoothing parameter for component $g$ and entry $k$, and the $(r,h)$th element of $\boldsymbol{\Phi}$ is given by $\phi_{rh}=\frac{1}{2}t_{r}^{2}(t_{h}-\frac{t_{r}}{3})$ for $t_{r}\leq t_{h}$. The matrix $\boldsymbol{\Phi}$ is common to all subjects since all entries of the multivariate time series are observed at common time points.

As seen above, the matrix $\boldsymbol{\Phi}$ is $n\times n$, and to avoid the computational burden for large $n$, a low-rank approximation is often adopted. To facilitate this approximation, we obtain basis functions via the spectral decomposition of $\boldsymbol{\Phi}$, as has been proposed in Wood and others (2002) and used in Rosen and others (2009); Rosen and others (2012); Krafty and others (2011). In particular, the matrix $\boldsymbol{W}$ consists of $m$ basis functions evaluated at times $t_{1},\ldots,t_{n}$, and $\boldsymbol{\beta}_{gk}$ is an $m$-dimensional vector of basis function coefficients. These basis functions are obtained by applying the spectral decomposition to $\boldsymbol{\Phi}$ such that $\boldsymbol{\Phi}=\boldsymbol{Q}\boldsymbol{\Gamma}\boldsymbol{Q}^{T}$, where $\boldsymbol{Q}$ is the matrix of eigenvectors of $\boldsymbol{\Phi}$, and $\boldsymbol{\Gamma}$ is a diagonal matrix containing the eigenvalues of $\boldsymbol{\Phi}$. We then let the design matrix $\boldsymbol{W}=\boldsymbol{Q}\boldsymbol{\Gamma}^{1/2}$ and place a normal prior $N(0,\tau^{2}_{gk}\boldsymbol{I}_{n})$ on $\boldsymbol{\beta}_{gk}$, which leads to $\boldsymbol{\mathcal{H}}_{gk}$ or $\boldsymbol{W}\boldsymbol{\beta}_{gk}\sim N(\boldsymbol{0},\tau^{2}_{gk}\boldsymbol{\Phi})$ as mentioned above.

By using the low-rank approximation, the number of columns of $\boldsymbol{W}$ is reduced from $n$ to $m$ ($m<n$), which greatly reduces the computational burden without sacrificing the model fit (Wahba, 1980; Wood, 2006). Eubank (1999) indicated that the eigenvalues in the diagonal matrix $\boldsymbol{\Gamma}$ decay rapidly as $m$ increases. Thus, we can achieve a good approximation by selecting a relatively small number $m$ of basis functions. The number of basis functions $m$ is set to $10$ in simulation studies as described in Section 6, which has been shown (Krafty and others 2011) to explain more than $98\%$ of the total variability.

The prior on $\boldsymbol{\theta}_{g}$ is thus $\boldsymbol{\theta}_{g}\sim N(\boldsymbol{0},\boldsymbol{D}_{g})$, where $\boldsymbol{D}_{g}=$ diag$(\sigma_{\alpha 1}^{2}\boldsymbol{1}_{2},\ \tau_{g1}^{2}\boldsymbol{1}_{m},\ \ldots\ ,\sigma_{\alpha K}^{2}\boldsymbol{1}_{2},\ \tau_{gK}^{2}\boldsymbol{1}_{m})$ is the covariance matrix of $\boldsymbol{\theta}_{g}$. The vector $(\sigma_{\alpha 1}^{2},\ldots,\sigma_{\alpha K}^{2})^{\prime}$ contains fixed prior variances for the regression coefficients $\boldsymbol{\alpha}_{gk}$, common to all components and entries. In particular, we fix the common prior variance $\sigma_{\alpha}^{2}=100$. The vector $\boldsymbol{\tau}_{g}^{2}=(\tau_{g1}^{2},\ldots,\tau_{gK}^{2})^{\prime}$ contains the smoothing parameters for the $g$th mixture component and $\boldsymbol{1}_{m}$ is an $m$-vector of ones. We assume independence between the regression coefficients $\boldsymbol{\alpha}_{gk}$ and the basis function coefficients $\boldsymbol{\beta}_{gk}$.

We assume the smoothing parameters $\boldsymbol{\tau}_{g}^{2}=(\tau_{g1}^{2},\ldots,\tau_{gK}^{2})^{\prime}$ vary across components $g$ and entries $k$. Although the most common choice for the prior on a variance parameter is the inverse gamma distribution, Gelman (2006) and Wand and others (2011) suggested that a half-$t$ prior on the standard deviation can reflect lack of information on a scale parameter. The half-$t$ is a family of heavy-tailed distributions and has a good shrinkage performance. It can be expressed as a scale mixture of inverse gamma random variables using a latent variable which follows an inverse gamma distribution (Wand and others, 2011). Thus, we assume a half-$t$ distribution such that $\tau_{gk}\sim t_{\nu_{\tau}}^{+}(0,A_{\tau})$, where $\nu_{\tau}$ is a degrees of freedom parameter, and $A_{\tau}$ is a scale parameter. We set $\nu_{\tau}=3$ and $A_{\tau}=10$ for all components and entries.

We assume $\sigma_{gk}\stackrel{{\scriptstyle\text{i.i.d}}}{{\sim}}t_{\nu_{\sigma}}^{+}(0,A_{\sigma})$ and set $\nu_{\sigma}=3$ and $A_{\sigma}=10$ for all components and entries.

This section provides details on the prior distributions placed on the parameters of the logistic weights (3). For ease of notation, we denote $\boldsymbol{\delta}_{g}^{*}=(\boldsymbol{\delta}_{g}^{T},\boldsymbol{\zeta}_{g}^{T})^{T}$, where $\boldsymbol{\zeta}_{g}=(\zeta_{1g},\cdots,\zeta_{Ng})^{T}$, $g=1,\ldots,G$. We let $\boldsymbol{V}_{i}^{*}=(\boldsymbol{V}_{i}^{\prime},\boldsymbol{e}_{i}^{\prime})^{\prime}$ where $\boldsymbol{e}_{i}$ is a vector of all zeros except for a single $1$ in the $i$th position, and $\boldsymbol{V}^{*}$ is a matrix consisting of the rows $\boldsymbol{V}_{i}^{*T}$, $i=1,\ldots,N$. Gaussian priors are placed on the logistic parameters, i.e., $\boldsymbol{\delta}_{g}^{*}\sim N(\boldsymbol{0},\boldsymbol{B}_{g})$, where $\boldsymbol{B}_{g}=\rm{diag}(\sigma_{\delta g}^{2}\boldsymbol{1}_{P+1},\ \kappa^{2}_{\zeta g}\boldsymbol{1}_{N})$, and the priors on the random intercepts satisfy $\boldsymbol{\zeta}_{g}\sim N(\boldsymbol{0},\kappa^{2}_{\zeta g}\boldsymbol{I}_{N})$. As for the hyperparameters, we assume $\sigma_{\delta g}^{2}=10$ for all components and covariates, and $\kappa_{\zeta g}\sim t_{\nu_{\kappa}}^{+}(0,A_{\kappa})$, where $\nu_{\kappa}=3$ and $A_{\kappa}=10$ for all components.

To sample the logistic parameters, Polson and others (2013) proposed a data augmentation scheme incorporating Pólya-Gamma latent variables, which facilitates Gibbs steps. Details on sampling the logistic parameters are provided in the Supplementary Material.

Letting $\ell$ denote the current Gibbs sampling iteration, parameter values at the $(\ell+1)$th iteration are drawn according to the following steps.

1. 1.

   Draw $\boldsymbol{\theta}_{gk}^{(\ell+1)}$ from $(\boldsymbol{\theta}_{gk}^{(\ell+1)}\mid\boldsymbol{y},\boldsymbol{S},\tau_{gk}^{2(\ell)},\sigma_{gk}^{2(\ell)})\sim N(\boldsymbol{u}_{gk},\sigma_{gk}^{2}\boldsymbol{\Lambda}_{gk})$, where $\boldsymbol{u}_{gk}$ and $\boldsymbol{\Lambda}_{gk}$ are mean vectors and covariance matrices.

2. 2.

   Draw $\sigma_{gk}^{2(\ell+1)}$ from $(\sigma_{gk}^{2(\ell+1)}\mid\boldsymbol{\epsilon}_{igk}^{(\ell+1)},a_{\sigma_{gk}}^{(\ell+1)})\sim IG\Big((nN_{g}^{(\ell)}+\nu_{\sigma})/2,\sum_{i=1}^{N}z_{ig}\boldsymbol{\epsilon}_{igk}^{\prime}\boldsymbol{\epsilon}_{igk}/2+\nu_{\sigma}/a_{\sigma_{gk}}\Big)$, where $N_{g}^{(\ell)}$ is the current number of subjects in the $g$th component, $\boldsymbol{\epsilon}_{igk}$ is the error vector for the $g$th component, the $i$th subject and the $k$th entry, and $a_{\sigma_{gk}}$ is a latent variable in the $IG$ scale mixture underlying the half-$t$ distribution.

3. 3.

   Draw $\tau_{gk}^{2(\ell+1)}$ from $(\tau_{gk}^{2(\ell+1)}\mid\boldsymbol{\beta}_{gk}^{(\ell+1)},a_{\tau_{gk}}^{(\ell+1)})\sim IG\Big((\nu_{\tau}+m)/2,\boldsymbol{\beta}_{gk}^{\prime}\boldsymbol{\beta}_{gk}/2+\nu_{\tau}/a_{\tau_{gk}}\Big)$, where $a_{\tau_{gk}}$ is a latent variable as in 2.

4. 4.

   Draw $\boldsymbol{\delta}_{g}^{*(\ell+1)}$ from $(\boldsymbol{\delta}_{g}^{*(\ell+1)}\mid\boldsymbol{V^{*}},z_{ig}^{(\ell)},\omega_{ig}^{(\ell+1)},\kappa^{2(\ell)}_{\zeta g})\sim N(\boldsymbol{M}_{g},\boldsymbol{\Sigma}_{g})$, where $\omega_{ig}^{(\ell+1)}$ is a Pólya-Gamma latent variable in the augmentation described in Section 4.4.

5. 5.

   Draw $\kappa^{2(\ell+1)}_{\zeta g}$ from $(\kappa^{2(\ell+1)}_{\zeta g}\mid\boldsymbol{\zeta}_{g}^{(\ell+1)},a_{\kappa_{g}}^{(\ell+1)})\sim IG\Big(\nu_{\kappa}/2,\boldsymbol{\zeta}_{g}^{\prime}\boldsymbol{\zeta}_{g}/2+(\nu_{\kappa}+N)/a_{\kappa_{g}}\Big)$, where $a_{\kappa_{g}}$ is a latent variable as in 2 and 3.

6. 6.

   The mixing weights $\pi_{ig}^{(\ell+1)}$ are obtained by computing $p(\pi_{ig}^{(\ell+1)}\mid\boldsymbol{V}^{*},\boldsymbol{\delta}_{g}^{*(\ell+1)},z_{ig}^{(\ell)})$ from Equation (3).

7. 7.

   Draw $z_{ig}^{(\ell+1)}\sim p(z_{ig}^{(\ell+1)}=1\mid\boldsymbol{y},\boldsymbol{S},\boldsymbol{\theta}_{gk}^{(\ell+1)},\sigma_{gk}^{2(\ell+1)},\pi_{ig}^{(\ell+1)})$ according to Equation (6).

Spiegelhalter and others (2002) suggested the use of the deviance information criterion (DIC) for model selection based on the effective number of parameters. Gelman and others (2003) introduced an alternative measure of effective number of parameters based on the variance of the log predictive density across MCMC iterations. This measure is robust and more accurate than the original one. Moreover, it has the advantages of always being positive and invariant to reparameterizations (Gelman and others, 2003).

In this paper, we use DIC to select the number of components for our proposed mixture model.

In this scenario, we consider the two-component trivariate model. From Equation (1), the $g$th component of the proposed mixture model is given by

where $\boldsymbol{y}(t_{j})$ is the trivariate time series evaluated at time $t_{j}$, $\boldsymbol{\alpha}_{01}=(1,-3,-2)^{\prime}$, $\boldsymbol{\alpha}_{02}=(5,4,3)^{\prime}$ and $\boldsymbol{\alpha}_{11}=(-2,2,0.5)^{\prime}$, $\boldsymbol{\alpha}_{12}=(1,-1,-0.5)^{\prime}$ are independent intercepts and slopes for each component, respectively. The vector $\boldsymbol{\beta}_{gq}$ consists of the $q$th spline coefficients of all variates for component $g$, and $w_{q}(t_{j})$ is the $q$th spline basis function evaluated at time $t_{j}$. The $\boldsymbol{\epsilon}_{gt_{j}}$ are independent zero-mean error terms, distributed as $\boldsymbol{\epsilon}_{gt_{j}}\sim\text{MVN}\Big(\boldsymbol{0},{\rm diag}(\sigma_{g1}^{2},\sigma_{g2}^{2},\sigma_{g3}^{2})\Big)$, where $\sigma_{1}^{2}=(\sigma_{11}^{2},\sigma_{12}^{2},\sigma_{13}^{2})^{\prime}=(3,5,4.5)^{\prime}$ and $\sigma_{2}^{2}=(\sigma_{21}^{2},\sigma_{22}^{2},\sigma_{23}^{2})^{\prime}=(4,3.5,4)^{\prime}$. The smoothing parameters are set to $\tau_{1}^{2}=(\tau_{11}^{2},\tau_{12}^{2},\tau_{13}^{2})^{\prime}=(3.5,5,8.5)^{\prime}$ and $\tau_{2}^{2}=(\tau_{21}^{2},\tau_{22}^{2},\tau_{23}^{2})^{\prime}=(6,2.5,1.5)^{\prime}$.

We investigate the performance of the trajectory and logistic parameter (see Equation (3)) estimates. For the former, we calculate the averaged root square error (ARSE) of each mixture component $g$

where $\mu_{gk}(t_{j})$ is the expectation of $y_{k}(t_{j})$ according to the $g$th component, and $y_{k}(t_{j})$ is the $k$th entry of the time series evaluated at time $t_{j}$. The $\hat{\mu}_{gk}(t_{j})$ are the estimated posterior means of $\mu_{gk}(t_{j})$ for $k=1,\ldots,K$ and $j=1,\ldots,n$.

To handle a potential label switching across mixture components, we compute $\text{ARSE}_{g}$ as the minimum value across all components, by using the estimate of the $g$th component and the truth of each group, $g=1,\ldots,G$. After obtaining correct component labels by evaluating ARSE, we also report the averaged bias (A-bias) and the variance of the bias (V-bias) of each mixture component $g$, where

and $\text{V-bias}_{g}$ is computed by calculating the sample variance of the bias over entries and time points.

For each replicate, time series trajectories are estimated by three methods: the proposed method, the R package gbmt (Magrini, 2022) and the TRAJ procedure in SAS (Nagin and others, 2018). Boxplots of ARSE, A-bias and V-bias of each component are given in Figure 3. Notably, TRAJ is able to fit a regression spline model by treating basis functions as time-varying covariates, while gbmt is only able to fit a cubic model. Our proposed method fits a penalized spline model under the Bayesian framework and is able to outperform both gbmt and TRAJ in terms of ARSE and V-bias for both components. A-biases are close to zero and comparable for all three methods. These findings demonstrate that all three methods are able to achieve a reasonable fit to group-based trajectories since bias over the entire time series is close to zero. Our proposed method is able to obtain more precise estimates of trajectories as is evident from the smaller V-biases.

To evaluate the performances of the logistic parameters, we compute the root mean squared error (RMSE) for each logistic parameter using the proposed method and TRAJ. Notably, gbmt is not able to incorporate covariates into the computation of mixing weights. Results of RMSEs of each logistic parameter are given in Table 1. We also compare RMSEs between the proposed method and TRAJ under four settings of different combinations of $N=150,250$ and $n=50,70$. Our proposed method yields smaller RMSEs of the logistic parameters in all cases, especially for the intercept $\delta_{0}$ and the first covariate $\delta_{1}$. This is to be expected since TRAJ uses a multinomial logistic model, which may result in inflated parameter estimates in cases of unbalanced outcomes or perfect separation, while our proposed method is able to obtain a shrinkage result using the penalization method.

In this scenario, we consider the four-component bivariate model whose $g$th component is given in Equation (7), where the values of the intercepts and slopes are $\boldsymbol{\alpha}_{01}=(1,-2)^{\prime}$, $\boldsymbol{\alpha}_{02}=(5,3)^{\prime}$, $\boldsymbol{\alpha}_{03}=(-3,5.5)^{\prime}$, $\boldsymbol{\alpha}_{04}=(4,-1)^{\prime}$, $\boldsymbol{\alpha}_{11}=(-3,0)^{\prime}$, $\boldsymbol{\alpha}_{12}=(2,-3.5)^{\prime}$, $\boldsymbol{\alpha}_{13}=(2.5,2)^{\prime}$ and $\boldsymbol{\alpha}_{14}=(-3,1.5)^{\prime}$. By analogy to the two-component trivariate model, the errors $\boldsymbol{\epsilon}_{gt_{j}}$ are independent zero-mean bivariate Gaussian random variables, distributed as $\boldsymbol{\epsilon}_{gt_{j}}\sim\text{MVN}\Big(\boldsymbol{0},{\rm diag}(\sigma_{g1}^{2},\sigma_{g2}^{2})\Big)$, where $\sigma_{1}^{2}=(\sigma_{11}^{2},\sigma_{12}^{2})^{\prime}=(6,9)^{\prime}$, $\sigma_{2}^{2}=(\sigma_{21}^{2},\sigma_{22}^{2})^{\prime}=(8,7.5)^{\prime}$, $\sigma_{3}^{2}=(\sigma_{31}^{2},\sigma_{32}^{2})^{\prime}=(10,6.5)^{\prime}$ and $\sigma_{4}^{2}=(\sigma_{41}^{2},\sigma_{42}^{2})^{\prime}=(7,8.5)^{\prime}$.

The performances of the estimated trajectories and logistic parameters for this scenario are displayed in Figure 4 and Table 2. As in the first scenario, our proposed method outperforms both gbmt and TRAJ in terms of ARSE and V-bias for all components. Notably, TRAJ fails to yield precise estimates in several replicates and thus results in larger mean ARSE and V-bias. In terms of the logistic parameters, the proposed method performs well with smaller RMSEs in almost all cases, especially for $\delta_{0}$ and $\delta_{1}$. More simulation results based on different values of $N$ and $n$ under the two scenarios considered above are presented in the Supplementary Material.

As described in Section 5 of the paper, Gibbs sampling is used to facilitate Bayesian inference. We denote by $\boldsymbol{\Theta}_{gk}=(\boldsymbol{\theta}_{gk}^{\prime},\tau_{gk}^{2},\sigma_{gk}^{2},\boldsymbol{\delta}_{g}^{*\prime},\kappa^{2}_{\zeta g})^{\prime}$ the parameters for the $g$th component and the $k$th entry, and the parameters in this vector are drawn from the corresponding conditional posterior distributions. Let $\ell$ be the current Gibbs sampling iteration; detailed Gibbs sampling steps for drawing the parameters at the $(\ell+1)$th iteration are given below.

1. 1.

   Sampling the basis function coefficients

   For each component $g$ and time series entry $k$, based on the augmented likelihood in Section 3.3 and the priors on $\boldsymbol{\theta}_{gk}=(\boldsymbol{\alpha}_{gk}^{\prime},\boldsymbol{\beta}_{gk}^{\prime})^{\prime}$ described in Section 4.1, the conditional posterior distribution of $(\boldsymbol{\theta}_{gk}^{(\ell+1)}\mid\boldsymbol{y},\boldsymbol{S},\tau_{gk}^{2(\ell)},\sigma_{gk}^{2(\ell)})$ is:

   	$\displaystyle p(\boldsymbol{\theta}_{gk}^{(\ell+1)}\mid\boldsymbol{y},\boldsymbol{S},$	$\displaystyle\tau_{gk}^{2(\ell)},\sigma_{gk}^{2(\ell)})\propto p(\boldsymbol{y}\mid\boldsymbol{S},\boldsymbol{\theta}_{gk}^{(\ell+1)},\sigma_{gk}^{2(\ell)})\cdot p(\boldsymbol{\theta}_{gk}^{(\ell+1)}\mid\tau_{gk}^{2(\ell)})$
   		$\displaystyle\propto\prod_{i=1}^{N}\Big\{(\sigma_{gk}^{2})^{-n/2}\exp\big[-\frac{1}{2\sigma_{gk}^{2}}(\boldsymbol{y}_{ik}-\boldsymbol{S}\boldsymbol{\theta}_{gk})^{\prime}(\boldsymbol{y}_{ik}-\boldsymbol{S}\boldsymbol{\theta}_{gk})\big]\Big\}^{z_{ig}}$
   		$\displaystyle\times|\boldsymbol{D}_{gk}|^{-1/2}\exp\Big(-\frac{1}{2}\boldsymbol{\theta}_{gk}\boldsymbol{D}_{gk}^{-1}\boldsymbol{\theta}_{gk}\Big)$
   		$\displaystyle\propto\exp\Big\{-\frac{1}{2\sigma_{gk}^{2}}\big[\sum_{i=1}^{N}z_{ig}(\boldsymbol{y}_{ik}-\boldsymbol{S}\boldsymbol{\theta}_{gk})^{\prime}(\boldsymbol{y}_{ik}-\boldsymbol{S}\boldsymbol{\theta}_{gk})+\boldsymbol{\theta}_{gk}^{\prime}\sigma_{gk}^{2}\boldsymbol{D}_{gk}^{-1}\boldsymbol{\theta}_{gk}\big]\Big\}$
   		$\displaystyle\propto\exp\Big[-\frac{1}{2\sigma_{gk}^{2}}(\boldsymbol{\theta}_{gk}-\boldsymbol{u}_{gk})^{\prime}(\boldsymbol{\Lambda}_{gk})^{-1}(\boldsymbol{\theta}_{gk}-\boldsymbol{u}_{gk})\Big]$
   		$\displaystyle\sim N(\boldsymbol{u}_{gk},\sigma_{gk}^{2}\boldsymbol{\Lambda}_{gk}),$

   where $\boldsymbol{\Lambda}_{gk}=(N_{g}^{(\ell)}\boldsymbol{S}^{\prime}\boldsymbol{S}+\sigma_{gk}^{2}\boldsymbol{D}_{gk}^{-1})^{-1}$, $\boldsymbol{u}_{gk}=\boldsymbol{\Lambda}_{gk}\sum_{i=1}^{N}z_{ig}\boldsymbol{S}^{\prime}\boldsymbol{y}_{ik}$, $N_{g}^{(\ell)}$ is the current number of subjects in the $g$th component, $\boldsymbol{D}_{gk}={\rm diag}(\sigma_{\alpha}^{2}\boldsymbol{1}_{2},\tau_{gk}^{2}\boldsymbol{1}_{m})$ is the prior covariance matrix for $\boldsymbol{\theta}_{gk}$. Hence, for each component $g$ and entry $k$, we draw $\boldsymbol{\theta}_{gk}^{(\ell+1)}$ from $(\boldsymbol{\theta}_{gk}^{(\ell+1)}\mid\boldsymbol{y},\boldsymbol{S},\tau_{gk}^{2(\ell)},\sigma_{gk}^{2(\ell)})\sim N(\boldsymbol{u}_{gk},\sigma_{gk}^{2}\boldsymbol{\Lambda}_{gk})$.

2. 2.

   Sampling the error variances

   Gelman (2006) proposed using the half-$t$ distribution as the prior on scale parameters. We follow Wand and others (2011) and express the half-$t$ prior of Section 4.3 as a scale mixture of inverse Gamma distributions as follows

   $$(\sigma_{gk}^{2}\mid a_{\sigma_{gk}})\sim IG\Big(\frac{\nu_{\sigma}}{2},\frac{\nu_{\sigma}}{a_{\sigma_{gk}}}\Big),a_{\sigma_{gk}}\sim IG\Big(\frac{1}{2},\frac{1}{A_{\sigma}^{2}}\Big).$$

   Therefore, the conditional posterior distribution of the latent variable $a_{\sigma_{gk}}$ is

   $$p(a_{\sigma_{gk}}^{(\ell+1)}\mid\sigma_{gk}^{2(\ell)})\propto\exp\Big[-\frac{1}{a_{\sigma_{gk}}}\Big(\frac{\nu_{\sigma}}{\sigma_{gk}^{2}}+\frac{1}{A_{\sigma}^{2}}\Big)\Big]\times(a_{\sigma_{gk}})^{-(\frac{1}{2}+1+\frac{\nu_{\sigma}}{2})},$$

   which is $IG\Big(\frac{\nu_{\sigma}+1}{2},\frac{\nu_{\sigma}}{\sigma_{gk}^{2}}+\frac{1}{A_{\sigma}^{2}}\Big)$. Denoting by $\boldsymbol{\epsilon}_{igk}$ the error vector of time series $\boldsymbol{y}_{ik}$ for component $g$, we have $\boldsymbol{\epsilon}_{igk}=\boldsymbol{y}_{ik}-\boldsymbol{S}\boldsymbol{\theta}_{gk}$, where $\boldsymbol{\epsilon}_{igk}\sim N(\boldsymbol{0},\sigma_{gk}^{2}\boldsymbol{I}_{n})$. The conditional distribution of the error variance is

   	$\displaystyle p(\sigma_{gk}^{2(\ell+1)}$	$\displaystyle\mid\boldsymbol{\epsilon}_{igk}^{(\ell+1)},a_{\sigma_{gk}}^{(\ell+1)})\propto p(\boldsymbol{\epsilon}_{igk}^{(\ell+1)}\mid\sigma_{gk}^{2(\ell+1)})\cdot p(a_{\sigma_{gk}}^{(\ell+1)}\mid\sigma_{gk}^{2(\ell+1)})\cdot p(\sigma_{gk}^{2(\ell+1)})$
   		$\displaystyle\propto\prod_{i=1}^{N}\Big[(\sigma_{gk}^{2})^{-\frac{n}{2}}\exp\Big(-\frac{1}{2\sigma_{gk}^{2}}\boldsymbol{\epsilon}_{igk}^{\prime}\boldsymbol{\epsilon}_{igk}\Big)\Big]^{z_{ig}}\times(\sigma_{gk}^{2})^{-(\frac{\nu_{\sigma}}{2}+1)}\exp\Big(-\frac{\nu_{\sigma}}{\sigma_{gk}^{2}a_{\sigma_{gk}}}\Big)$
   		$\displaystyle\propto(\sigma_{gk}^{2})^{-(\frac{n}{2}N_{g}^{(\ell)}+\frac{\nu_{\sigma}}{2}+1)}\exp\Big[-\frac{1}{\sigma_{gk}^{2}}\Big(\frac{\sum_{i=1}^{N}z_{ig}\boldsymbol{\epsilon}_{igk}^{\prime}\boldsymbol{\epsilon}_{igk}}{2}+\frac{\nu_{\sigma}}{a_{\sigma_{gk}}}\Big)\Big],$

   which is $IG\Big(\frac{nN_{g}^{(\ell)}+\nu_{\sigma}}{2},\frac{\sum_{i=1}^{N}z_{ig}\boldsymbol{\epsilon}_{igk}^{\prime}\boldsymbol{\epsilon}_{igk}}{2}+\frac{\nu_{\sigma}}{a_{\sigma_{gk}}}\Big)$. The sampling scheme proceeds by first sampling $(a_{\sigma_{gk}}^{(\ell+1)}\mid\sigma_{gk}^{2(\ell)})$ and then $(\sigma_{gk}^{2(\ell+1)}\mid\boldsymbol{\epsilon}_{igk}^{(\ell+1)},a_{\sigma_{gk}}^{(\ell+1)})$.

3. 3.

   Sampling the smoothing parameters

   The smoothing parameters $\tau_{gk}^{2}$ are drawn by analogy to the error variances. We first draw $(a_{\tau_{gk}}^{(\ell+1)}\mid\tau_{gk}^{2(\ell)})\sim IG\Big(\frac{\nu_{\tau}+1}{2},\frac{\nu_{\tau}}{\tau_{gk}^{2}}+\frac{1}{A_{\tau}^{2}}\Big)$. The conditional posterior distribution of the smoothing parameters is

   	$\displaystyle p(\tau_{gk}^{2(\ell+1)}\mid\boldsymbol{\beta}_{gk}^{(\ell+1)},a_{\tau_{gk}}^{(\ell+1)})$	$\displaystyle\propto p(\boldsymbol{\beta}_{gk}^{(\ell+1)}\mid\tau_{gk}^{2(\ell+1)})\cdot p(a_{\tau_{gk}}^{(\ell+1)}\mid\tau_{gk}^{2(\ell+1)})\cdot p(\tau_{gk}^{2(\ell+1)})$
   		$\displaystyle\propto(\tau_{gk}^{2})^{-\frac{m+\nu_{\tau}}{2}}\exp{\Big[-\frac{1}{\tau_{gk}^{2}}\Big(\frac{\nu_{\tau}}{a_{\tau_{gk}}}+\frac{\boldsymbol{\beta}_{gk}^{\prime}\boldsymbol{\beta}_{gk}}{2}\Big)\Big]},$

   which is $IG\Big(\frac{\nu_{\tau}+m}{2},\frac{\boldsymbol{\beta}_{gk}^{\prime}\boldsymbol{\beta}_{gk}}{2}+\frac{\nu_{\tau}}{a_{\tau_{gk}}}\Big)$. The sampling scheme proceeds by first sampling $(a_{\tau_{gk}}^{(\ell+1)}\mid\tau_{gk}^{2(\ell)})$ and then $(\tau_{gk}^{2(\ell+1)}\mid\boldsymbol{\beta}_{gk}^{(\ell+1)},a_{\tau_{gk}}^{(\ell+1)})$.

4. 4.

   Sampling the logistic parameters

   Let $\boldsymbol{\delta}_{g}^{*}=(\boldsymbol{\delta}_{g}^{T},\boldsymbol{\zeta}_{g}^{T})^{T}$ be the aggregation of the logistic parameters and all random intercepts for the $g$th component. Based on the logits of Section 3.2 and the corresponding priors described in Section 4.4, the conditional posterior distribution of $(\boldsymbol{\delta}_{g}^{*(\ell+1)}\mid\boldsymbol{V}^{*},z_{ig}^{(\ell)},\kappa^{2(\ell)}_{\zeta g})$ is

   	$\displaystyle p(\boldsymbol{\delta}_{g}^{*(\ell+1)}\mid\boldsymbol{V}^{*},z_{ig}^{(\ell)},\kappa^{2(\ell)}_{\zeta g})$	$\displaystyle\propto p(z_{ig}^{(\ell)}=1\mid\boldsymbol{V}^{*},\boldsymbol{\delta}_{g}^{*(\ell+1)})\cdot p(\boldsymbol{\delta}_{g}^{*(\ell+1)}\mid\kappa^{2(\ell)}_{\zeta g})$
   		$\displaystyle=\prod_{i=1}^{N}\Big[\frac{\exp(\boldsymbol{V}_{i}^{*\prime}\boldsymbol{\delta}_{g}^{*})}{\sum_{h=1}^{G}\exp(\boldsymbol{V}_{i}^{*\prime}\boldsymbol{\delta}_{h}^{*})}\Big]^{z_{ig}}p(\boldsymbol{\delta}_{g}^{*(\ell+1)}\mid\kappa^{2(\ell)}_{\zeta g}),$

   where $\boldsymbol{V}^{*}=(\boldsymbol{V}_{1}^{*},\ldots,\boldsymbol{V}_{N}^{*})^{\prime}$ is a $N\times(P+1)$ matrix with $\boldsymbol{V}_{i}^{*}$ representing all covariates (including intercepts) for subject $i$. To sample from the posterior distribution of $p(\boldsymbol{\delta}_{g}^{*(\ell+1)}\mid\boldsymbol{V}^{*},z_{ig}^{(\ell)},\kappa^{2(\ell)}_{\zeta g})$, we adopt the Póyla-Gamma data augmentation strategy of Polson and others (2013) by introducing a latent variable $\omega_{ig}$ coming from the Pólya-Gamma distribution. Thus, the conditional posterior distributions of the logistic parameters are

   	$\displaystyle p(\boldsymbol{\delta}_{g}^{*(\ell+1)}\mid\boldsymbol{V}^{*},z_{ig}^{(\ell)},\omega_{ig}^{(\ell+1)},\kappa^{2(\ell)}_{\zeta g})$	$\displaystyle\propto p(z_{ig}^{(\ell)}=1\mid\boldsymbol{V}^{*},\omega_{ig}^{(\ell+1)},\boldsymbol{\delta}_{g}^{*(\ell+1)})\cdot p(\omega_{ig}^{(\ell+1)}\mid\boldsymbol{V}^{*},\boldsymbol{\delta}_{g}^{*(\ell)})$
   		$\displaystyle\cdot p(\boldsymbol{\delta}_{g}^{*(\ell+1)}\mid\kappa^{2(\ell)}_{\zeta g})$
   		$\displaystyle\propto\exp\Big(-\frac{\omega_{ig}\eta_{ig}^{2}}{2}\Big)\cdot p(\omega_{ig}\mid 1,0)|\boldsymbol{B}_{g}|^{-P/2}\exp\Big(-\frac{1}{2}\boldsymbol{\delta}_{g}^{*\prime}\boldsymbol{B}_{g}^{-1}\boldsymbol{\delta}_{g}^{*}\Big),$

   where $\eta_{ig}=\boldsymbol{V}_{i}^{*\prime}\boldsymbol{\delta}_{g}^{*}-C_{ig}$ and $C_{ig}=\log\sum_{h\neq j}\exp(\boldsymbol{V}_{i}^{*\prime}\boldsymbol{\delta}_{h}^{*})$, $p(\omega_{ig}\mid 1,0)$ is the Pólya-gamma distribution $PG(b,c)$ with $b=1$ and $c=0$, $\boldsymbol{B}_{g}$ is the prior covariance matrix of Section 4.4 and $\boldsymbol{B}_{g}={\rm diag}(\sigma_{\delta g}^{2}\boldsymbol{1}_{P+1},\kappa^{2}_{\zeta g}\boldsymbol{1}_{N})$. By assuming the conjugate prior $N(\boldsymbol{0},\boldsymbol{B}_{g})$ on $\boldsymbol{\delta}_{g}^{*}$, the posterior distribution of the Pólya-gamma latent variable is

   $$(\omega_{ig}^{(\ell+1)}\mid\boldsymbol{V}^{*},\boldsymbol{\delta}_{g}^{*(\ell)})\sim PG(1,\eta_{ig}).$$

   Thus, the conditional distributions of the logistic parameters (including the random intercepts) are

   $$(\boldsymbol{\delta}_{g}^{*(\ell+1)}\mid\boldsymbol{V}^{*},z_{ig}^{(\ell)},\omega_{ig}^{(\ell+1)},\kappa^{2(\ell)}_{\zeta g})\sim N(\boldsymbol{M}_{g},\boldsymbol{\Sigma}_{g}),$$

   where $\boldsymbol{\Sigma}_{g}=(\boldsymbol{V}^{*\prime}\boldsymbol{\Omega}_{g}\boldsymbol{V}^{*}+\boldsymbol{B}_{g}^{-1})^{-1}$, $\boldsymbol{M}_{g}=\boldsymbol{\Sigma}_{g}\big[\boldsymbol{V}^{*\prime}(\boldsymbol{\Omega}_{g}\boldsymbol{C}_{g}+\boldsymbol{\xi}_{g})\big]$, $\boldsymbol{\Omega}_{g}={\rm diag}(\omega_{1g},\cdots,\omega_{Ng})$, $\boldsymbol{C}_{g}=(C_{1g},\cdots,C_{Ng})^{\prime}$, and $\boldsymbol{\xi}_{g}=(\xi_{1g},\cdots,\xi_{Ng})^{\prime}$, with $\xi_{ig}=z_{ig}-\frac{1}{2}$. Thus, $\boldsymbol{\delta}_{g}^{*(\ell+1)}$ is drawn by first sampling $(\omega_{ig}^{(\ell+1)}\mid\boldsymbol{V}^{*},\boldsymbol{\delta}_{g}^{*(\ell)})$ and then $(\boldsymbol{\delta}_{g}^{*(\ell+1)}\mid\boldsymbol{V}^{*},z_{ig}^{(\ell)},\omega_{ig}^{(\ell+1)},\kappa^{2(\ell)}_{\zeta g})$.

5. 5.

   Sampling the variances of the random intercepts

   By analogy with sampling the error variances and sthe moothing parameters, we first draw $(a_{\kappa_{g}}^{(\ell+1)}\mid\kappa^{2(\ell)}_{\zeta g})\sim IG\Big(\frac{\nu_{\kappa}+1}{2},\frac{\nu_{\kappa}}{\kappa^{2}_{\zeta g}}+\frac{1}{A_{\kappa}^{2}}\Big)$. The conditional posterior distributions of the variances of the random intercepts are

   	$\displaystyle p(\kappa^{2(\ell+1)}_{\zeta g}\mid\boldsymbol{\zeta}_{g}^{(\ell+1)},a_{\kappa_{g}}^{(\ell+1)})$	$\displaystyle\propto p(\boldsymbol{\zeta}_{g}^{(\ell+1)}\mid\kappa^{2(\ell+1)}_{\zeta g})\cdot p(a_{\kappa_{g}}^{(\ell+1)}\mid\kappa^{2(\ell+1)}_{\zeta g})\cdot p(\kappa^{2(\ell+1)}_{\zeta g})$
   		$\displaystyle\propto(\kappa^{2}_{\zeta g})^{-\frac{N+\nu_{\kappa}}{2}+1}\exp{\Big[-\frac{1}{\kappa^{2}_{\zeta g}}\Big(\frac{\nu_{\kappa}}{a_{\kappa_{g}}}+\frac{\boldsymbol{\zeta}_{g}^{T}\boldsymbol{\zeta}_{g}}{2}\Big)\Big]},$

   which is $IG\Big(\frac{\nu_{\kappa}+N}{2},\frac{\boldsymbol{\zeta}_{g}^{T}\boldsymbol{\zeta}_{g}}{2}+\frac{\nu_{\kappa}}{a_{\kappa_{g}}}\Big)$. The sampling scheme proceeds by first sampling $(a_{\kappa_{g}}^{(\ell+1)}\mid\kappa^{2(\ell)}_{\zeta g})$ and then $(\kappa^{2(\ell+1)}_{\zeta g}\mid\boldsymbol{\zeta}_{g}^{(\ell+1)},a_{\kappa_{g}}^{(\ell+1)})$.

6. 6.

   Computing the mixing weights

   After drawing the $\boldsymbol{\delta}_{g}^{*}$, the mixing weights $\pi_{ig}^{(\ell+1)}$ for each component, given the design matrix $V_{i}^{*}$, are computed by

   $$p(\pi_{ig}^{(\ell+1)}\mid V_{i}^{*},\boldsymbol{\delta}_{g}^{*(\ell+1)})=\frac{\exp(\boldsymbol{V}_{i}^{*T}\boldsymbol{\delta}_{g}^{*})}{\sum_{h=1}^{G}\exp(\boldsymbol{V}_{i}^{*T}\boldsymbol{\delta}_{h}^{*})}.$$

7. 7.

   Sampling the latent indicators

   After sampling all parameters and computing the mixing weights, the final Gibbs step is to allocate subjects to different components by drawing the latent indicators $z_{ig}$. As in Section 3.3, the conditional posterior of these indicators is

   $$p(z_{ig}^{(\ell+1)}=1\mid\boldsymbol{y},\boldsymbol{S},\boldsymbol{\Theta}^{(\ell+1)},\pi_{ig}^{(\ell+1)})=\frac{\pi_{ig}\prod_{k=1}^{K}f_{gk}(\boldsymbol{y}_{ik}\mid\boldsymbol{\Theta}_{gk})}{\sum_{h=1}^{G}\pi_{ih}\prod_{k=1}^{K}f_{hk}(\boldsymbol{y}_{ik}\mid\boldsymbol{\Theta}_{hk})},$$

   and the indicators are drawn from the multinomial distribution.