Bayesian Clustering for Continuous-Time Hidden Markov Models

TL;DR

This paper introduces novel Bayesian clustering methods for continuous-time hidden Markov models, utilizing advanced MCMC techniques for improved inference in both finite and infinite mixture models.

Contribution

It develops new Bayesian clustering algorithms for CTHMMs using reversible-jump and split-merge MCMC methods, enhancing inference for mixture models.

Findings

Algorithms perform well on simulated data

Effective inference for finite and infinite mixtures

Demonstrated applicability on real data

Abstract

We develop clustering procedures for longitudinal trajectories based on a continuous-time hidden Markov model (CTHMM) and a generalized linear observation model. Specifically in this paper, we carry out finite and infinite mixture model-based clustering for a CTHMM and achieve inference using Markov chain Monte Carlo (MCMC). For a finite mixture model with prior on the number of components, we implement reversible-jump MCMC to facilitate the trans-dimensional move between different number of clusters. For a Dirichlet process mixture model, we utilize restricted Gibbs sampling split-merge proposals to expedite the MCMC algorithm. We employ proposed algorithms to the simulated data as well as a real data example, and the results demonstrate the desired performance of the new sampler.

Tables

Table 1. Table 1: Prior distribution for number of clusters, p 0 ​ ( M ∗ ) subscript 𝑝 0 superscript 𝑀 p_{0}(M^{*}) for Dirichlet mixture models (DMM) with α = 0.5 𝛼 0.5 \alpha=0.5 and mixture of finite mixtures (MFMs) with δ = 1 𝛿 1 \delta=1 and M ∼ P ​ o ​ i ​ s ​ s ​ o ​ n ​ ( α ​ l ​ o ​ g ​ ( N ) ) + 1 similar-to 𝑀 𝑃 𝑜 𝑖 𝑠 𝑠 𝑜 𝑛 𝛼 𝑙 𝑜 𝑔 𝑁 1 M\sim Poisson(\alpha log(N))+1 .

	$M^{\ast }$
$N$	Model	1	2	3	4	5	6	7	8	9	10	11	12	13	$\ge$14
1000	DMM	0.028	0.105	0.191	0.225	0.194	0.131	0.073	0.034	0.014	0.005	0.001	0.000	0.000	0.000
	MFM	0.037	0.119	0.188	0.193	0.183	0.126	0.090	0.038	0.021	0.004	0.001	0.000	0.000	0.000
3000	DMM	0.016	0.069	0.146	0.200	0.202	0.160	0.104	0.058	0.027	0.011	0.004	0.001	0.000	0.000
	MFM	0.024	0.075	0.160	0.174	0.198	0.153	0.104	0.058	0.030	0.015	0.008	0.001	0.000	0.000
5000	DMM	0.013	0.057	0.127	0.185	0.200	0.169	0.118	0.070	0.036	0.016	0.006	0.002	0.001	0.000
	MFM	0.012	0.068	0.126	0.181	0.177	0.169	0.122	0.076	0.045	0.012	0.008	0.002	0.002	0.000
10000	DMM	0.009	0.044	0.104	0.164	0.192	0.178	0.134	0.086	0.048	0.024	0.010	0.004	0.001	0.000
	MFM	0.005	0.044	0.098	0.168	0.195	0.168	0.129	0.085	0.065	0.023	0.013	0.006	0.001	0.000
25000	DMM	0.006	0.030	0.079	0.137	0.176	0.180	0.151	0.107	0.067	0.037	0.018	0.008	0.003	0.001
	MFM	0.009	0.034	0.082	0.139	0.169	0.166	0.154	0.106	0.074	0.031	0.022	0.012	0.002	0.000

Table 2. Table 2: Example 1: Inference for simulated data with three clusters and three latent states.

		Gaussian			Poisson
		$T=30$	$T=50$	$T=100$	$T=30$	$T=50$	$T=100$
% 3-cluster iterations	MFM-RJ	62.35%	54.40%	81.05%	66.85%	80.65%	89.95%
	MFM-SM	45.25%	46.85%	50.25%	48.65%	56.60%	58.15%
	DMM	59.65%	99.25%	99.75%	56.40%	93.05%	98.25%
Misclassification rate	MFM-RJ	1.20%	0.60%	0.00%	4.00%	2.30%	0.70%
	MFM-SM	15.60%	12.70%	9.10%	12.10%	7.30%	2.40%
	DMM	20.97%	19.16%	12.9%	14.10%	6.53%	2.40%

Table 3. Table 3: Example 2: Simulation study with three clusters. Each cluster has three latent states, the same B 𝐵 B and π 𝜋 \pi parameters, but different Q 𝑄 Q matrices.

# of cluster	$\sigma =0.5$			$\sigma =1$			$\sigma =2$
	MFM-RJ	DMM	MFM-SM	MFM-RJ	DMM	MFM-SM	MFM-RJ	DMM	MFM-SM
1	1.20%	0.05%	0.10%	0.65%	0.10%	0.10%	0.05%	0.15%	0.25%
2	38.50%	0.20%	26.05%	0.55%	0.05%	5.70%	0.75%	18.15%	13.95%
3	58.05%	80.90%	54.10 %	61.05%	95.75%	46.30%	26.15%	37.40%	37.20 %
4	1.45%	17.05%	17.25%	31.95%	4.10%	32.60%	20.95%	38.25%	28.55%
5	0.80%	1.80%	2.55%	5.80%	0.00%	11.95%	19.80%	6.05%	12.80%
6	0.00%	0.00%	2.30 %	0.00%	0.00%	7.80%	1.80%	0.00%	6.10%
$\ge 7$	0.00%	0.00%	0.00%	0.00%	0.00%	1.60%	14.20%	0.00%	1.15%
% Misclassification	4.90%	19.20%	15.40%	9.50%	22.13%	19.30%	24.20%	24.70%	25.30%

Table 4. Table 4: Example 3: Simulation study with three clusters with one time-varying factor covariate. Each cluster has three latent states and Poisson observation process.

		$T=30$	$T=50$	$T=100$
% 3-cluster iterations	MFM-RJ	97.55%	98.65%	98.55%
	MFM-SM	65.90%	65.90%	87.50%
	DMM	37.20%	85.05%	95.60%
% Misclassification	MFM-RJ	3.50%	1.90%	0.60%
	MFM-SM	12.90%	7.00%	3.60%
	DMM	11.50%	6.20%	2.90%

Table 5. Table 5: COPD data analysis: Co-clustering matrix for the two clustering methods. Cluster membership determined by the posterior modal cluster.

	Finite mixture model
		Cluster 1	Cluster 2	Cluster 3	Cluster 4
Dirichlet mixture model	Cluster 1	1002	563	584	697
	Cluster 2	322	9950	2421	1005
	Cluster 3	535	1077	2903	847
	Cluster 4	249	654	613	1290

Table 6. Table 6: Example 1: Inference for simulated data with three clusters and three latent states via a DP mixture model.

	Gaussian			Poisson
	$T=30$	$T=50$	$T=100$	$T=30$	$T=50$	$T=100$
$\Vert {\pi }_1-{\hat{\pi }}_1\Vert$	0.04	0.02	0.02	0.04	0.02	0.04
$\Vert B_1-{\hat{B}}_1\Vert$	0.47	0.43	0.41	0.22	0.08	0.04
$\Vert Q_1-{\hat{Q}}_1\Vert$	0.79	0.60	0.37	1.37	0.44	0.22
$\Vert {\pi }_2-{\hat{\pi }}_2\Vert$	0.04	0.03	0.02	0.05	0.04	0.01
$\Vert B_2-{\hat{B}}_2\Vert$	0.37	0.44	0.28	0.02	0.01	0.02
$\Vert Q_2-{\hat{Q}}_2\Vert$	0.59	0.58	0.27	0.88	0.45	0.17
$\Vert {\pi }_3-{\hat{\pi }}_3\Vert$	0.02	0.01	0.02	0.03	0.03	0.03
$\Vert B_3-{\hat{B}}_3\Vert$	0.58	0.31	0.25	0.12	0.03	0.02
$\Vert Q_3-{\hat{Q}}_3\Vert$	0.34	0.25	0.09	0.30	0.12	0.16
% of 3-cluster iterations	59.65%	99.25%	99.75%	56.40%	93.05%	98.25%
Misclassification rate	20.97%	19.16%	12.9%	14.10%	6.53%	2.40%

Table 7. Table 7: Example 1: Inference for simulated data with three clusters and three latent states via a mixture of finite mixture models using reversible-jump MCMC.

	Gaussian			Poisson
	$T=30$	$T=50$	$T=100$	$T=30$	$T=50$	$T=100$
$\Vert {\pi }_1-{\hat{\pi }}_1\Vert$	0.03	0.02	0.04	0.04	0.03	0.03
$\Vert B_1-{\hat{B}}_1\Vert$	0.10	0.03	0.04	0.13	0.20	0.13
$\Vert Q_1-{\hat{Q}}_1\Vert$	0.21	0.19	0.11	0.38	0.32	0.20
$\Vert {\pi }_2-{\hat{\pi }}_2\Vert$	0.03	0.03	0.02	0.04	0.03	0.04
$\Vert B_2-{\hat{B}}_2\Vert$	0.10	0.02	0.03	0.04	0.02	0.04
$\Vert Q_2-{\hat{Q}}_2\Vert$	0.09	0.06	0.03	0.10	0.21	0.09
$\Vert {\pi }_3-{\hat{\pi }}_3\Vert$	0.03	0.01	0.05	0.06	0.03	0.03
$\Vert B_3-{\hat{B}}_3\Vert$	0.09	0.03	0.02	0.05	0.05	0.02
$\Vert Q_3-{\hat{Q}}_3\Vert$	0.13	0.04	0.04	0.06	0.05	0.04
% of 3-cluster iterations	62.35%	54.40%	81.05%	66.85%	80.65%	89.95%
Misclassification rate	1.20%	0.60%	0.00%	4.00%	2.30%	0.70%

Table 8. Table 8: Example 1: Inference for simulated data with three clusters and three latent states via a mixture of finite mixture models using a split-merge algorithm.

	Gaussian			Poisson
	$T=30$	$T=50$	$T=100$	$T=30$	$T=50$	$T=100$
$\Vert {\pi }_1-{\hat{\pi }}_1\Vert$	0.04	0.02	0.04	0.07	0.03	0.03
$\Vert B_1-{\hat{B}}_1\Vert$	0.59	0.33	0.06	0.08	0.18	0.14
$\Vert Q_1-{\hat{Q}}_1\Vert$	0.86	0.45	0.10	0.94	0.65	0.31
$\Vert {\pi }_2-{\hat{\pi }}_2\Vert$	0.04	0.04	0.04	0.09	0.04	0.04
$\Vert B_2-{\hat{B}}_2\Vert$	0.38	0.23	0.17	0.43	0.03	0.04
$\Vert Q_2-{\hat{Q}}_2\Vert$	0.35	0.21	0.19	0.88	0.45	0.06
$\Vert {\pi }_3-{\hat{\pi }}_3\Vert$	0.08	0.02	0.03	0.12	0.05	0.03
$\Vert B_3-{\hat{B}}_3\Vert$	0.98	0.72	0.16	0.26	0.02	0.04
$\Vert Q_3-{\hat{Q}}_3\Vert$	0.64	0.45	0.05	0.54	0.12	0.08
% of 3-cluster iterations	45.25%	46.85%	50.25%	48.65%	56.60%	58.15%
Misclassification rate	15.60%	12.70%	9.10%	12.10%	7.30%	2.40%

Table 9. Table 9: Example 2: Simulation study with three clusters. Each cluster has three latent states, the same B 𝐵 B and π 𝜋 \pi parameters, but different Q 𝑄 Q matrices.

	MFM-RJ			DMM			MFM-SM
	$\sigma =0.5$	$\sigma =1$	$\sigma =2$	$\sigma =0.5$	$\sigma =1$	$\sigma =2$	$\sigma =0.5$	$\sigma =1$	$\sigma =2$
$\Vert \pi -\hat{\pi }\Vert$	0.03	0.03	0.06	0.01	0.01	0.07	0.06	0.06	0.05
$\Vert B-\hat{B}\Vert$	0.02	0.05	0.01	0.02	0.02	0.01	0.01	0.01	0.08
$\Vert Q_1-{\hat{Q}}_1\Vert$	0.30	0.26	0.78	0.32	0.51	0.77	0.79	0.51	1.08
$\Vert Q_2-{\hat{Q}}_2\Vert$	0.08	0.10	0.49	0.38	0.34	0.36	0.21	0.42	0.73
$\Vert Q_3-{\hat{Q}}_3\Vert$	0.05	0.07	0.87	0.14	0.16	0.64	0.68	0.85	0.92
% of 3-cluster iterations	58.05%	61.05%	44.10%	80.90%	95.75%	25.90%	54.10%	46.30%	37.20%
% of Misclassification	4.90%	9.50%	24.20%	19.20%	22.13%	24.70%	15.40%	19.30%	25.30%

Table 10. Table 10: Example 3: Simulation study with three clusters with one factor time-varying covariate. Each cluster has three latent states and the observation process is generated from Poisson distribution.

	MFM-RJ			DMM			MFM-SM
	$T=30$	$T=50$	$T=100$	$T=30$	$T=50$	$T=100$	$T=30$	$T=50$	$T=100$
$\Vert {\pi }_1-{\hat{\pi }}_1\Vert$	0.01	0.02	0.02	0.08	0.06	0.07	0.07	0.07	0.08
$\Vert B_1-{\hat{B}}_1\Vert$	0.75	0.69	0.14	0.99	0.70	0.08	1.15	1.06	0.07
$\Vert Q_1-{\hat{Q}}_1\Vert$	0.34	0.10	0.15	0.34	0.94	0.37	0.73	0.61	0.13
$\Vert {\pi }_2-{\hat{\pi }}_2\Vert$	0.06	0.03	0.02	0.10	0.05	0.04	0.09	0.08	0.03
$\Vert B_2-{\hat{B}}_2\Vert$	0.06	0.22	0.03	0.21	0.11	0.04	0.20	0.17	0.11
$\Vert Q_2-{\hat{Q}}_2\Vert$	0.11	0.14	0.16	0.47	0.47	0.39	0.48	0.31	0.26
$\Vert {\pi }_3-{\hat{\pi }}_3\Vert$	0.02	0.02	0.06	0.12	0.10	0.12	0.16	0.08	0.13
$\Vert B_3-{\hat{B}}_3\Vert$	0.41	0.30	0.12	0.42	0.36	0.19	0.51	0.68	0.09
$\Vert Q_3-{\hat{Q}}_3\Vert$	0.12	0.12	0.08	0.39	0.39	0.31	0.43	0.64	0.22
% of 3-cluster iterations	97.55%	98.65%	98.55%	37.20%	85.05%	95.60%	65.90%	65.90%	87.50%
Misclassification rate	3.50%	1.90%	0.60%	11.50%	6.20%	2.90%	12.90%	7.00%	3.60%

Table 11. Table 11: 24,712 COPD Patients: Finite mixture model analysis: exponential of coefficients for the three cluster Poisson model with 95% credible intervals in relative-risk parameterization with GP as the baseline group.

	Predictor	State 1	State 2	State 3	State 4
Cluster 1	Intercept	0.29	1.35	3.18	10.71
		(0.28,0.31)	(1.33,1.37)	(3.16,3.21)	(10.62,10.79)
	ER	1.08	1.12	1.13	1.02
		(0.96,1.21)	(1.08,1.17)	(1.11,1.16)	(1.00,1.03)
$\left(N=2108\right)$	HOSP	1.71	1.19	1.07	0.99
		(1.53,1.92)	(1.14,1.25)	(1.04,1.10)	(0.97,1.00)
	SEPC	0.72	0.93	0.99	0.99
		(0.62,0.83)	(0.90,0.97)	(0.97,1.02)	(0.97,1.01)
Cluster 2	Intercept	0.00	2.96	4.89	7.08
		(0.00,0.00)	(2.94,2.99)	(4.86,4.91)	(7.05,7.10)
	ER	1.49	0.94	1.00	1.02
		(0.76,2.76)	(0.93,0.96)	(0.99,1.00)	(1.01,1.03)
$\left(N=12244\right)$	HOSP	4.10	0.91	0.96	0.99
		(2.33,7.02)	(0.90,0.93)	(0.95,0.97)	(0.98,1.00)
	SEPC	4.04	0.96	0.95	0.98
		(2.40,6.79)	(0.95,0.97)	(0.94,0.96)	(0.97,0.98)
Cluster 3	Intercept	1.02	2.25	4.21	8.33
		(1.00,1.05)	(2.22,2.27)	(4.19,4.22)	(8.29,8.37)
	ER	0.37	1.06	1.02	1.03
		(0.35,0.40)	(1.04,1.08)	(1.01,1.03)	(1.02,1.04)
$\left(N=6521\right)$	HOSP	0.37	1.06	1.00	1.00
		(0.33,0.40)	(1.03,1.09)	(0.98,1.01)	(0.98,1.01)
	SEPC	0.82	0.97	0.99	1.00
		(0.77,0.87)	(0.95,0.99)	(0.98,1.00)	(0.99,1.01)
Cluster 4	Intercept	0.10	1.81	5.86	9.16
		(0.09,0.11)	(1.79,1.93)	(5.83,5.89)	(9.12,9.20)
	ER	0.09	0.95	0.98	1.03
		(0.06,0.12)	(0.93,0.98)	(0.97,0.99)	(1.02,1.04)
$\left(N=3839\right)$	HOSP	0.17	0.90	0.94	1.00
		(0.12,0.24)	(0.87,0.93)	(0.92,0.95)	(0.99.1.01)
	SEPC	0.44	0.99	0.98	0.99
		(0.34,0.56)	(0.96,1.02)	(0.96,0.99)	(0.98,1.01)

Table 12. Table 12: 24,712 COPD Patients: Dirichlet process mixture model analysis: exponential of coefficients for the three cluster Poisson model with 95% credible intervals in relative-risk parameterization with GP as the baseline group.

	Predictor	State 1	State 2	State 3	State 4
Cluster 1	Intercept	0.24	4.63	3.84	10.06
		(0.21,0.26)	(4.61,4.66)	(3.81,3.87)	(10.00,10.12)
	ER	0.05	1.10	1.03	1.02
		(0.02,0.08)	(1.06,1.14)	(1.02,1.05)	(1.01,1.03)
$\left(N=2846\right)$	HOSP	0.34	1.12	1.00	0.99
		(0.25,0.44)	(1.08,1.18)	(0.98,1.02)	(0.98,1.01)
	SEPC	0.41	1.01	0.99	0.98
		(0.32,0.53)	(0.97,1.04)	(0.98,1.01)	(0.97,1.00)
Cluster 2	Intercept	0.00	2.61	5.08	7.28
		(0.00,0.00)	(2.59,2.62)	(5.06,5.10)	(7.26,7.30)
	ER	0.86	0.93	1.00	1.03
		(0.36,1.80)	(0.92,0.94)	(0.99,1.01)	(1.02,1.04)
$\left(N=13698\right)$	HOSP	1.40	0.88	0.96	1.00
		(0.53,3.13)	(0.87,0.90)	(0.96,0.97)	(0.99,1.00)
	SEPC	0.38	0.94	0.97	0.98
		(0.07,1.29)	(0.93,0.95)	(0.96,0.98)	(0.98,0.99)
Cluster 3	Intercept	0.14	1.44	4.52	8.69
		(0.13,0.15)	(1.43,1.46)	(4.50,4.55)	(8.65,8.72)
	ER	0.17	1.06	1.01	1.03
		(0.13,0.22)	(1.03,1.08)	(1.00,1.02)	(1.01,1.04)
$\left(N=5362\right)$	HOSP	0.22	1.07	0.97	1.00
		(0.16,0.30)	(1.04,1.11)	(0.95,0.98)	(0.99,1.01)
	SEPC	0.33	0.96	0.96	0.99
		(0.25,0.42)	(0.93,0.98)	(0.95,0.97)	(0.98,1.00)
Cluster 4	Intercept	0.23	1.98	4.01	8.54
		(0.21,0.25)	(1.95,2.01)	(3.98,4.03)	(8.49,8.60)
	ER	0.08	0.99	1.06	1.03
		(0.05,0.13)	(0.96,1.02)	(1.04,1.08)	(1.01,1.04)
$\left(N=2806\right)$	HOSP	0.07	0.95	1.02	0.99
		(0.02,0.16)	(0.91,0.99)	(1.00,1.04)	(0.98.1.01)
	SEPC	0.33	0.96	0.98	1.01
		(0.23,0.46)	(0.92,0.99)	(0.97,1.00)	(0.99,1.02)

Table 13. Table 13: 24,712 COPD Patients: Finite mixture model, effective sample sizes (ESSs) for parameters in Q 𝑄 Q conditional on 1636 four-cluster iterations

Parameter	ESS	Parameter	ESS	Parameter	ESS	Parameter	ESS
Cluster 1
$q_{12}$	300.61	$q_{21}$	217.26	$q_{31}$	265.46	$q_{41}$	1042.76
$q_{13}$	288.48	$q_{23}$	515.78	$q_{32}$	217.66	$q_{42}$	591.93
$q_{14}$	851.18	$q_{24}$	377.74	$q_{34}$	973.16	$q_{43}$	483.49
Cluster 2
$q_{12}$	309.20	$q_{21}$	562.45	$q_{31}$	431.60	$q_{41}$	683.10
$q_{13}$	404.99	$q_{23}$	516.30	$q_{32}$	238.04	$q_{42}$	287.36
$q_{14}$	461.82	$q_{24}$	156.82	$q_{34}$	370.57	$q_{43}$	185.15
Cluster 3
$q_{12}$	384.46	$q_{21}$	217.28	$q_{31}$	345.52	$q_{41}$	1165.85
$q_{13}$	467.94	$q_{23}$	511.19	$q_{32}$	289.33	$q_{42}$	900.17
$q_{14}$	1038.50	$q_{24}$	714.76	$q_{34}$	684.81	$q_{43}$	854.58
Cluster 4
$q_{12}$	390.26	$q_{21}$	376.44	$q_{31}$	410.11	$q_{41}$	664.70
$q_{13}$	360.58	$q_{23}$	569.69	$q_{32}$	295.72	$q_{42}$	366.11
$q_{14}$	635.86	$q_{24}$	298.25	$q_{34}$	720.29	$q_{43}$	282.69

Table 14. Table 14: 24,712 COPD Patients: DP mixture model, effective sample sizes (ESSs) for parameters in Q 𝑄 Q conditional on 1894 four-cluster iterations

Parameter	ESS	Parameter	ESS	Parameter	ESS	Parameter	ESS
Cluster 1
$q_{12}$	413.62	$q_{21}$	286.26	$q_{31}$	468.32	$q_{41}$	1161.28
$q_{13}$	508.90	$q_{23}$	741.59	$q_{32}$	374.72	$q_{42}$	685.85
$q_{14}$	964.17	$q_{24}$	570.95	$q_{34}$	1271.31	$q_{43}$	1059.18
Cluster 2
$q_{12}$	876.39	$q_{21}$	715.00	$q_{31}$	398.21	$q_{41}$	819.26
$q_{13}$	452.95	$q_{23}$	750.99	$q_{32}$	411.18	$q_{42}$	366.51
$q_{14}$	586.39	$q_{24}$	298.37	$q_{34}$	1356.51	$q_{43}$	240.72
Cluster 3
$q_{12}$	306.68	$q_{21}$	292.76	$q_{31}$	407.86	$q_{41}$	1031.25
$q_{13}$	516.50	$q_{23}$	794.15	$q_{32}$	338.10	$q_{42}$	449.77
$q_{14}$	762.50	$q_{24}$	446.57	$q_{34}$	1052.53	$q_{43}$	712.01
Cluster 4
$q_{12}$	372.38	$q_{21}$	390.12	$q_{31}$	407.82	$q_{41}$	820.60
$q_{13}$	291.55	$q_{23}$	400.30	$q_{32}$	303.40	$q_{42}$	550.24
$q_{14}$	743.86	$q_{24}$	576.38	$q_{34}$	1052.53	$q_{43}$	719.27

Equations

$$f\left(O_t|X_{{\tau }_t}=k\right)=\exp \left[\left\{O_t{\theta }_{t,k}-b({\theta }_{t,k})\right\}/{\varphi }^2+c(O_t,\varphi )\right].$$

$$p\left(\mathbf{C},\Theta |\mathbf{O},\mathbf{X},M\right)\propto p_0\left(\mathbf{C}|M\right)\prod _{n=1}^N\mathcal{L}\left(O_n,X_n|{\Theta }_{C_n}\right)$$

$${\varpi }_1,\dots ,{\varpi }_M|M\sim Dirichlet\left(\delta ,\dots ,\delta \right)$$

$$p_0\left(\mathbf{C}|M,\delta \right)=\frac{\Gamma \left(M\delta \right)}{\Gamma \left(N+M\delta \right)}\prod _{m=1}^M\frac{\Gamma \left(N_m+\delta \right)}{\Gamma \left(\delta \right)}$$

$$p_0\left(\mathcal{Z}|\delta \right)=V_N\left(M^{\ast }\right)\prod _{m=1}^M{\delta }^{N_m}$$

$$V_N\left(M^{\ast }\right)=\sum _{M=1}^{\infty }\frac{M\left(M-1\right)\cdots \left(M-M^{\ast }+1\right)}{\delta M(\delta M+1)\cdots \left(\delta M+N-1\right)}{p}_0\left(M\right),\mathcal{Z}=\{E_m:\left|E_m\right|>0\}$$

$$p_0\left(\mathcal{Z}|\alpha \right)=\frac{{\alpha }^M\Gamma \left(\alpha \right)}{\Gamma \left(N+\alpha \right)}\prod _{m=1}^M\Gamma \left(N_m\right)$$

$$p_0\left(M\right)=\left|S_{N,M}\right|\frac{{\alpha }^M\Gamma \left(\alpha \right)}{\Gamma \left(N+\alpha \right)}{p}_0\left(\mathcal{Z}|M\right)=\frac{\prod _{m=1}^M\Gamma \left(N_m\right)}{\left|S_{N,M}\right|}$$

$$\sum _{i=1}^N\frac{\alpha }{\alpha +i-1}=\alpha \left(\psi (\alpha +N)-\psi (\alpha )\right)\bumpeq \alpha \log \left(1+\frac{N}{\alpha }\right)$$

$$M$$

$${\varpi }_1,\dots ,{\varpi }_M|M$$

$$\mathbb{P}\left(C_n=m|{\varpi }_1,\dots ,{\varpi }_M,M\right)$$

$$X_n|\Theta ,C_n$$

$$O_n|X_n,\Theta ,C_n$$

$$\mathcal{L}\left(O_n,X_n|C_n,\Theta \right)=\prod _{m=1}^M{\left[{\varpi }_m\mathcal{L}\left(O_n,X_n|C_n=m,{\Theta }_m\right)\right]}^{\text{mathbbm}1\left(C_n=m\right)}.$$

$$\mathbb{P}\left(C_n=m|O_n,X_n,\Theta \right)=\frac{{\varpi }_m\mathcal{L}\left(O_n,X_n|C_n=m,{\Theta }_m\right)}{\sum _{l=1}^M{\varpi }_l\mathcal{L}\left(O_n,X_n|C_n=l,{\Theta }_l\right)}.$$

$$\mathcal{L}(\mathbf{o}|\Theta ,M)=\prod _{n=1}^N\left\{\sum _{m=1}^M{\varpi }_m\mathcal{L}(O_n|{\Theta }_m)\right\}$$

$$\alpha \left(M^{\prime },{\Theta }^{\prime };M,\Theta \right)$$

$$p\left(M,\Theta |\mathbf{o}\right)\propto \mathcal{L}(\mathbf{o}|\Theta ,M)p_0(\Theta |M)p_0(M).$$

$$\pi^{p}=\left(\pi{{}^{\prime}},\pi{{}^{\prime\prime}}\right)\qquad Q^{p}=\left(\begin{array}[]{cc}Q{{}^{\prime}}&0\\ 0&Q{{}^{\prime\prime}}\\ \end{array}\right)\qquad B^{p}=\left(B{{}^{\prime}},B{{}^{\prime\prime}}\right).$$

$$r\left(M+1,({\Theta }^p,{\varpi }^{\prime },{\varpi }^{\prime \prime });M,({\Theta }_M,{\varpi }_M)|\mathbf{o}\right)=\frac{p\left(M+1,({\Theta }^p,{\varpi }^{\prime },{\varpi }^{\prime \prime })|\mathbf{o}\right)}{p\left(M,({\Theta }_M,{\varpi }_M)|\mathbf{o}\right)}.$$

$$\min \left(1,\frac{\mathsf{q}\left(Q_M;Q^p\right)}{\mathsf{q}\left(Q^p;Q_M\right)}\frac{\mathsf{q}\left(B_M;B^p\right)}{\mathsf{q}\left(B^p;B_M\right)}\frac{\mathsf{q}\left({\pi }_M;{\pi }^p\right)}{\mathsf{q}\left({\pi }^p;{\pi }_M\right)}\frac{\mathsf{q}\left({\varpi }_M;\varpi {}^{\prime },\varpi {}^{\prime \prime }\right)}{\mathsf{q}\left(\varpi {}^{\prime },\varpi {}^{\prime \prime };{\varpi }_M\right)}r\left(M+1,({\Theta }^p,{\varpi }^{\prime },{\varpi }^{\prime \prime });M,({\Theta }_M,{\varpi }_M)|\mathbf{o}\right)\right)$$

$$=\min \left(1,\frac{d_{M+1}{\varpi }_M}{b_M{p}_{\varpi }(w)p_{\beta }({\beta }_{1,k}^{{}^{\prime \prime }})}r\left(M+1,({\Theta }^p,{\varpi }^{\prime },{\varpi }^{\prime \prime });M,({\Theta }_M,{\varpi }_M)|\mathbf{o}\right)\right)$$

$$q_{i,k}^{{}^{\prime }}=\frac{s_{M,i}}{s_{M,i}+s_{M+1,i}}\times q_{M,i,k}+\frac{s_{M+1,i}}{s_{M,i}+s_{M+1,i}}\times q_{M+1,i,k},i\ne k=1,\dots ,K$$

$${\beta }_{i,k}^{{}^{\prime }}=\frac{s_{M,i}}{s_{M,i}+s_{M+1,i}}\times {\beta }_{M,i,k}+\frac{s_{M+1,i}}{s_{M,i}+s_{M+1,i}}\times {\beta }_{M+1,i,k},i,k=1,\dots ,K$$

$${\pi }_k^{{}^{\prime }}=\frac{s_{M,i}}{s_{M,i}+s_{M+1,i}}\times {\pi }_{M,k}+\frac{s_{M+1,i}}{s_{M,i}+s_{M+1,i}}\times {\pi }_{M+1,k},k=1,\dots ,K$$

$$\min \left(1,\frac{b_M{p}_{\varpi }(w)}{d_{M+1}{\varpi }^{{}^{\prime }}}r\left(M,\Theta {}^{\prime };M+1,{\Theta }_M,{\Theta }_{M+1}|\mathbf{o}\right)\right).$$

$$G({\tilde{\Theta }}_n)=\sum _{m=1}^{\infty }{\varpi }_m\text{mathbbm}1({\tilde{\Theta }}_n={\Theta }_m)\text{with}{\varpi }_m=v_m\prod _{k<m}(1-v_k)$$

$${\pi }_C^i\sim Dirichlet\left(\sum _{n:C_n^i=C}{S}_{n,1,1}^i+{\alpha }_1,\dots ,\sum _{n:C_n^i=C}{S}_{n,1,K}^i+{\alpha }_K\right)$$

$$a\left({\mathbf{C}}^{\ast },\mathbf{C}\right)=\min \left\{1,\frac{q\left(\mathbf{C}|{\mathbf{C}}^{\ast }\right)}{q\left({\mathbf{C}}^{\ast }|\mathbf{C}\right)}\frac{{\pi }_0\left({\mathbf{C}}^{\ast }\right)}{{\pi }_0\left(\mathbf{C}\right)}\frac{\mathcal{L}\left({\mathbf{C}}^{\ast }\right)}{\mathcal{L}\left(\mathbf{C}\right)}\right\}$$

$$\frac{{\pi }_0\left({\mathbf{C}}^{\text{split}}\right)}{{\pi }_0\left(\mathbf{C}\right)}=\alpha \frac{\left(N_{C_d^{\text{split}}}-1\right)!\left(N_{C_e^{\text{split}}}-1\right)!}{\left(N_{C_d}-1\right)!}\frac{{\pi }_0\left({\mathbf{C}}^{\text{merge}}\right)}{{\pi }_0\left(\mathbf{C}\right)}=\frac{1}{\alpha }\frac{\left(N_{C_d^{\text{merge}}}-1\right)!}{\left(N_{C_d}-1\right)!\left(N_{C_e}-1\right)!}$$

$$\mathcal{L}\left(C_f\right)$$

$${\mathcal{L}}_f\left(\Theta \right)=\prod _{t=1}^{T_f}\prod _{k=1}^{K}f{\left(O_{f,t}|S_{f,t,k}\right)}^{S_{f,t,k}}\prod _{k=1}^K{\pi }_k^{S_{f,1,k}}\prod _{t=1}^{T_f}\prod _{l\ne m}\mathcal{L}\left(q_{lm}|{\Delta }_{f,t}\right).$$

Full text

Keywords: Model-based clustering; continuous-time hidden Markov models; reversible jump MCMC; split-merge proposal; nonparametric Bayesian inference;mixture models.

1 Introduction

Continuous-time Markov processes on a finite state space have been widely used to represent longitudinal data, especially if the times between successive observations are irregular. Likelihood inference for the infinitesimal generator of a continuous-time Markov jump process has been studied by Billingsley (1961). However, in reality, the process is only observed at discrete time points, and inference for the generator of the process becomes difficult. This problem arises in a range of practical settings ranging from public health surveillance (Luo et al., 2021) to molecular dynamics (Hobolth and Stone, 2009). Inference for discretely-observed continuous-time Markov processes has been explored by Bladt and Sørensen (2005) in both likelihood and Bayesian frameworks.

Continuous-time Markov processes on a finite state space have been widely used to represent longitudinal data, especially if the times between successive observations are irregular. Likelihood inference for the infinitesimal generator of a continuous-time Markov jump process has been studied by Billingsley (1961). However, in reality, the process is only observed at discrete time points, and inference for the generator of the process becomes difficult. This problem arises in a range of practical settings ranging from public health surveillance (Luo et al., 2021) to molecular dynamics (Hobolth and Stone, 2009). Inference for discretely-observed continuous-time Markov processes has been explored by Bladt and Sørensen (2005) in both likelihood and Bayesian frameworks.

In most settings, observed data are either not direct observations of the Markov process, or the process is observed with measurement errors. In such cases, we may consider a latent process as an unobserved ‘trajectory’ recording the true but unobserved state of nature. In this setting, a hidden Markov model (HMM) is suitable, and assumes that the Markov property is imposed on the unobserved process. There is a broad interest in the application of continuous-time HMMs (CTHMMs) (see for example Jackson et al., 2003; Lange et al., 2015), with the majority of research devoted to frequentist approaches; however, recently Luo et al. (2021) and Williams et al. (2020) constructed Bayesian CTHMMs using different missing data likelihood formulations for the underlying Markov chain.

If the cohort is assumed homogeneous with respect to the stochastic properties of the latent model, a common continuous-time model may be fitted to the entire data set. However, it is plausible that the study base comprises different sub-cohorts that have distinct stochastic properties. In this paper, we develop model-based clustering procedures to cluster individuals according to their sub-cohort, and to study the pattern in each cluster. In Section 2, we review the development of a continuous time HMM/generalized linear model (CTHMM-GLM) for non-equidistant longitudinal data. Section 3 discusses model-based clustering approaches; Section 3.1 compares the prior distribution of finite and infinite mixture models that are commonly used for clustering models. A reversible-jump MCMC approach to clustering under finite mixture models is presented in Section 4. Section 5 presents Dirichlet process mixture model-based clustering implemented using restricted Gibbs sampling split-merge proposals. Simulated examples and a real data example to examine the performance of model-based clustering are presented in Sections 6 and 7 respectively, and we discuss possible in Section 8.

2 The CTHMM-GLM model

We presume that, for a single individual, a sequence $\{O_{1},\ldots,O_{T}\}$ of health status-related variables is observed at time points $\{\tau_{1},\dots,\tau_{T}\}$. A latent process $\{X_{s}\}$ ($s\in\mathbb{R}^{+}$), representing the health status for a condition of interest, is assumed to be a continuous-time Markov chain (CTMC) with parameters $\left(\pi,Q\right)$, where $\pi$ is the initial distribution and $Q$ is the infinitesimal generator, taking values on the finite state space $\{1,2,\ldots,K\}$. The assumption of a finite state space for the latent process is justified in many real cases, such as the real example in Luo et al. (2021), and corresponds to conceptual stages in the progression of the underlying disease process. However such discrete models are undoubtedly an approximation in general. The observation process $O\left|{X_{s},B}\right.$ is presumed to follow an exponential family, with

[TABLE]

A GLM with link function specified as $g(u(\theta_{t,k}))={\mathbf{Z}_{{t}}}^{\top}{\beta_{k}}$ can be incorporated if there are explanatory variables $\mathbf{Z}\in\mathbb{R}^{D}$, where $u(\theta_{t,k})=\mathbb{E}\left({{O_{t}}\left|X_{\tau_{t}}=k\right.}\right)$ and $\beta_{k}$ is a coefficient vector for state $k$. Let $S_{t}=\left(S_{t,1},\ldots,S_{t,K}\right)^{\top}$ be an indicator random vector with $S_{t,k}=1$ if $X_{\tau_{t}}=k$ and 0 otherwise. We can rewrite the linear predictor into a matrix form as $g(u(\theta_{t,k}))=\mathbf{Z}_{{t}}^{\top}BS_{t}$, with $B=(\beta_{d,k})$ for $d=1,\ldots,D$ and $k=1,\ldots,K$, a coefficient matrix containing all GLM coefficients for each latent state. Luo et al. (2021) developed likelihood and Bayesian inference procedures for this model using the expectation-maximization (EM) algorithm and MCMC approaches. In the formulation, observations are indexed using an integer index (that is, $O_{t}$), and the latent process using a continuous-valued index (that is, $X_{\tau_{t}}$). In this paper, we assume that the measurement process itself (that is, the collection of times $\tau_{t},t=1,\ldots,T$) is not informative about the system either in its hidden or observed components.

This model, developed in Luo et al. (2021) (the Supplement gives details of the likelihood construction and the Bayesian hierarchical model), is a parametric model for which the likelihood is reasonably complex, albeit one that is simplified in its representation using a latent process. The model is presumed to apply to all individuals in the study, who are presumed a random sample from the target population. In this paper we consider an extension to allow for systematic heterogeneity to be exhibited by sub-populations of subjects.

3 Model-based clustering

The contribution of this paper is to develop model-based clustering for data presumed to be generated by the complex model in Section 2; specifically we develop clustering approaches for the latent trajectories based on the observed data and the presumed hidden Markov structure.

Model-based clustering is typically achieved via parametric likelihood- or density-based calculations, with the number of clusters selected using information criteria, such as AIC or BIC (Fraley and Raftery, 1998). This approach was explored extensively by Luo (2019) in the context of CTHMMs for modeling health trajectories. However, in such calculations, the number of clusters has had to be pre-specified. Bayesian model determination approaches have been a longstanding focus of interest in Bayesian inference (see, for example, Carlin and Chib, 1995; Green, 1995; Godsill, 2001). Two approaches are typically adopted to address this issue; first, reversible-jump Markov chain Monte Carlo (MCMC) (Green, 1995) exploits trans-dimensional Metropolis-Hastings (MH) moves, allowing movement across parameter spaces of different dimensions; secondly Bayesian nonparametric procedures based on the Dirichlet process are also widely used – these models are often termed as infinite mixture models, where a prior is placed on the space of discrete random measures, and where the models are limiting versions of exchangeable finite mixture models. Dirichlet process models are now widely used with implementation facilitated via MCMC; key relevant references in a well-established literature include Escobar and West (1995); MacEachern and Müller (1998); Neal (2000); Ishwaran and James (2001); Jain and Neal (2004).

The principal challenge in our setting is that the core model on which the clustering will be based is driven by the unobserved continuous time trajectory represented by $\{X_{s}\}$. We meet this challenge by implementing MCMC algorithms based on the complete data likelihood described in Section 2.

3.1 Clustering via mixture models

The general principle behind model-based clustering is to cluster individuals based on the component model parameters that determine the mixture form. The basic formulation of the model envisages that the population is composed of distinct sub-populations each with distinct stochastic properties. In the case of the CTHMM formulation, this corresponds to each subpopulation having a potentially different parameter $\Theta=(\pi,Q,B)$. The estimated parameter for each cluster will provide a subpopulation-level summary, and will allow comparisons within and between clusters.

Let $C_{n}$ be the cluster membership indicator for individual $n$ with cluster-specific model parameters $\Theta_{C_{n}}=\left\{\pi_{C_{n}},Q_{C_{n}},B_{C_{n}}\right\}$. The likelihood contribution for this individual is $\mathcal{L}\left(O_{n},X_{n}\left|\Theta_{C_{n}}\right.\right)$. With $M$ denoting the number of mixture components, a prior probability, $p_{0}\left(\mathbf{C}\left|M\right.\right)$, would be assigned to each cluster membership partition $\mathbf{C}=\left\{C_{1},\ldots,C_{N}\right\}$, resulting in the posterior form

[TABLE]

where $\mathbf{C}$ is the a specific partition of the $N$ individuals into $M^{*}$ non-empty clusters. For practical inference, the key is whether the form of $\mathcal{L}\left(O_{n},X_{n}\left|\Theta_{C_{n}}\right.\right)$ can be easily computed, and the cluster-specific parameter $\Theta_{C_{n}}$ can be integrated out. This would greatly simplify the MCMC algorithm. In our model specification, conditional on a proposed clustering model, we will adopt the formulation of Luo et al. (2021) for the CTHMM, so therefore differences in inference will be driven by the specific clustering model adopted. In next sections, we compare aspects of the finite mixture and Dirichlet process mixture models, noting their similarities and where they differ.

3.2 Mixture of finite mixtures

Conventionally, in a finite mixture model, the number of clusters is determined by fitting an $M$-component mixture model with various $M$ and then using information criteria to choose among them. A natural extension is to regard the number of components, $M$, as a random variable, specifying the prior as $M\sim p_{0}\left(M\right)$, a mass function on $\left\{1,2,3,\ldots\right\}$, resulting in the joint prior distribution of the form $p_{0}\left(\mathbf{C},M\right)=p_{0}\left(\mathbf{C}|M\right)p_{0}\left(M\right)$. However, there is a crucial distinction between the number of components $M$ in the mixture model and the number of clusters $M^{*}$ in the data which is defined as the number of components used to generate the observed data, or the number of “filled” mixture components. By specifying a prior on the number of components $M$, we implicitly place a prior on $M^{*}$. We typically specify mixture weights

[TABLE]

with $\delta=1$, making the weight distribution uniform, and then draw cluster labels $\mathbf{C}$ independently from the multinomial distribution with $\mathbb{P}\left({C_{n}=m}\left|M,\varpi_{1},\ldots,\varpi_{M}\right.\right)=\varpi_{m}$. The resulting conditional distribution of $\mathbf{C}$ given that there are $M$ components has the form

[TABLE]

where $N_{m}$ is the number of subjects placed in component $m$ by this procedure.

Miller and Harrison (2018) discussed similarities and differences between the Dirichlet mixture model and a mixture of finite mixtures (MFMs). They showed that when $N\to\infty$, the distribution of the number of clusters $M^{*}$ behaves similarly to the number of components $M$, and under the posterior they also behave similarly. In addition, they also illustrated that this prior is a simple exchangeable partition distribution closely paralleling that of the Dirichlet process, with the distribution of the random partition

[TABLE]

where

[TABLE]

and $E_{m}=\{n:C_{n}=m\}$ for $m\in\{1,2,\ldots\}$. Therefore, posterior samples can be generated via a direct application of MCMC for the Dirichlet process mixture, and the only difference is that the new element is placed in an existing cluster $m$ with probability $\propto\left(N_{m}+\delta\right)$ or a new cluster with probability $\propto\delta\frac{V_{N}\left(M^{*}+1\right)}{V_{N}\left(M^{*}\right)}$. In Section 4, we will explicitly introduce this model in the context of CTHMM and inference based on reversible-jump MCMC (Green, 1995) with an efficient proposal. In addition, we can also directly apply the algorithm developed in Section 5 for Dirichlet mixture models to this model by replacing $\alpha$ with $\delta\frac{V_{N}\left(M^{*}+1\right)}{V_{N}\left(M^{*}\right)}$ and $N_{m}$ with $N_{m}+\delta$.

3.3 Dirichlet mixture models

In the finite mixture formulation, allowing $M\longrightarrow\infty$, given a finite sample size $N$, will still yield a valid prior distribution albeit with some $N_{m}$’s equal zero. This leads to another choice of a prior distribution, the Dirichlet process (Ferguson, 1973; Antoniak, 1974; Lo, 1984). Specifically, for a Dirichlet process prior, the prior distribution of the partition has the form

[TABLE]

and the probability mass function of $M$ and the conditional distribution of $\mathbf{C}$ given $M$ are (Green and Richardson, 2001)

[TABLE]

where $\left|S_{N,M}\right|=\sum_{\text{all }\mathcal{Z}}\prod\limits_{m=1}^{M}\Gamma\left(N_{m}\right)$, is the absolute value of a Stirling number of the first kind, the number of ways of partitioning $N$ items into $M$ non-empty subsets.

In Section 5, we will use the Dirichlet process prior in the context of CTHMMs. This representation marginalizes out parameters of the model, which is particularly useful in terms of the computation. Specifically, a split-merge MCMC for this model will be introduced to facilitate the computation of updating cluster memberships.

3.4 Calibrating the prior on the number of clusters

In order to make our analysis under the finite and Dirichlet mixture models as comparable as possible, we aim to match their cluster specifications. For the parameters in the CTHMM itself, it is possible to match the specification exactly, so we focus on implied the cluster structures. The Dirichlet process prior on $M$, given by the first expression in (4), can be computed directly for any $N$; this is numerically challenging when $N$ is large, so rather than compute it exactly we adopt a Monte Carlo strategy, as sampling the number of clusters is straightforward: under this prior $M^{*}\stackrel{{\scriptstyle d}}{{=}}\sum_{i=1}^{N}\mathbbm{1}\left\{U_{i}>\frac{\alpha}{\alpha+i-1}\right\}$, where $U_{1},\ldots,U_{N}$ are independent $Uniform(0,1)$ random variables, and $\mathbbm{1}\{.\}$ is the indicator function. Table 1 shows $p_{0}(M^{*})$ computed for $\alpha=0.5$ and $N=1000,3000,5000,10000,25000$, with Monte Carlo simulation being used to compute for the two largest sample sizes.

The expected number of clusters under this prior is

[TABLE]

where $\psi(.)$ is the Digamma function. Thus the bulk of the mass in this distribution is concentrated on low values of $M^{*}$ even when $N$ is very large, and the expectation increases very slowly with $N$. To make approaches comparable, in the simulation and real data analysis, we will use this $p_{0}(M^{*})$ as the prior in the finite mixture model discussed in Section 3.2.

While under MFMs, we can use the similar procedure to sample the number of clusters by replacing $\alpha$ with $\delta\frac{V_{N}\left(M^{*}+1\right)}{V_{N}\left(M^{*}\right)}$ and $N_{m}$ with $N_{m}+\delta$. We use the Monte Carlo simulation to generate the empirical distribution of $p_{0}\left(M^{*}\right)$ (shown in Table 1). Under the prior of $M\sim Poisson(\alpha log(N))+1$, the prior of $M^{*}$ is similar to $Poisson(\alpha log(N))+1$ as well as with Dirichlet mixture models. In terms of cluster sizes, as is evident from the form in (4), the Dirichlet process prior favours few, large clusters (see, for example, Green and Richardson, 2001; Miller and Harrison, 2018), whereas under the symmetric prior with $\delta=1$, the finite mixture model specification puts higher mass on clusters of the same order of magnitude.

4 Computation for the Mixture of Finite Mixtures CTHMM

We deploy reversible-jump MCMC procedures to allow the number of clusters to be inferred via the posterior distribution under a finite mixture model formulation. Specifically, we use a split-combine move to update the number of clusters and then implement fixed dimension MCMC.

We start with finite mixture model-based clustering for CTHMM in this section. We have the following hierarchy:

[TABLE]

The complete-data likelihood for subject $n$ is

[TABLE]

A subject is assigned to cluster $m$ according to the probability

[TABLE]

In reality, the model parameters, $\Theta$, and the values of the latent states, $X_{n}$, are not known, and must be inferred from the observed data.

4.1 Reversible-jump MCMC

The reversible-jump algorithm (Green, 1995) can be deployed to implement trans-dimensional MCMC for the finite mixture model with ${M^{*}}$ unknown. We study split/combine moves for pairs of states similar in spirit to the split/merge moves of Richardson and Green (1997) and Dellaportas and Papageorgiou (2006). If the model consists of $M$ components with $K$ latent states ($K$ is assumed fixed) in each component, then the model can be viewed as a CTHMM-GLM with $K\times{M}$ states. Infinitesimal generator $Q$ can be expressed as a block diagonal matrix with diagonal blocks $Q_{1},\dots,Q_{m},\dots,Q_{M}$, where $Q_{m}$ is the $K\times K$ infinitesimal generator matrix for cluster $m$ $\left(m=1,\dots,{M}\right)$. This constraint prevents the transition of subjects between clusters across time. One iteration of the algorithm includes

a move that considers splitting a cluster into two, or combining two clusters into one; 2. 2.

an update of the cluster label for each individual according to the posterior probability (5), given parameters in each cluster, and the latent states; 3. 3.

an update of the model parameters using standard MCMC moves for each cluster with the number of clusters $M$ fixed, specifically proposals that

• •

update latent state indicators $S_{n,t}$;

• •

update the parameters associated with the observation process $B$;

• •

update the initial distribution $\pi$;

• •

update the infinitesimal generator $Q$.

For any empty cluster from Step 2, we generate model parameters from prior distributions.

For the split and combine moves, we will implement a reversible jump algorithm (Green, 1995). We carry out this move on the marginalized model, where the cluster labels and latent processes are marginalized out from the calculation, and use the likelihood

[TABLE]

where $\Theta$ denotes the collection of $\varpi$ and CTMC parameters.

Consider a proposal from the current state $(M,\Theta)$ to a new state $(M^{\prime},\Theta^{\prime})$ using the proposal density $\mathsf{q}\left(M^{\prime},\Theta^{\prime};M,\Theta\right)=\mathsf{q}_{1}\left(M^{\prime};M\right)\mathsf{q}_{2}\left(\Theta^{\prime};\Theta\right)$, that is, using independent proposals for the two components. The acceptance probability for this proposal is given by

[TABLE]

where $p\left(M,\Theta\left|\mathbf{o}\right.\right)$ is the posterior distribution of $(M,\Theta)$ given the observed data $\mathbf{o}$, which can, up to proportionality, be decomposed into the marginal (or ‘incomplete data’) likelihood of the data $\mathcal{L}(\mathbf{o}|\Theta,M)$ times the prior distribution for $(M,\Theta)$, with prior distribution as $p_{0}$;

[TABLE]

Our algorithm relies upon the ability to compute the marginal likelihood efficiently for any $\Theta$; however, this is a standard ‘forward’ calculation for continuous-time hidden Markov models. The reversible-jump algorithm is for the most part standard, and the only noteworthy elements are the trans-dimensional moves. We discuss these in more detail below.

4.2 Split and combine move

To construct efficient split and combine moves, we adopt the idea of centered proposals (Brooks et al., 2003) designed to produce similar likelihood contributions for the current and proposed parameters. The combine move is designed to choose a cluster, $m$ say, at random and select another cluster $i$ such that $\left\|B_{i}-B_{m}\right\|_{2}$ is smallest for $i\neq m$. The reverse split move is to randomly select a cluster, $m$ to split into two clusters, say $m$ and $m^{*}$, and check if the condition, $\left\|B_{m*}-B_{m}\right\|_{2}<\left\|B_{j}-B_{m}\right\|_{2}$ for $j\neq m$. If this condition is not met, then the split move is rejected directly.

4.2.1 Split move

We consider an update that changes the number of clusters from $M\to M+1$. Without loss of generality, we aim to split the $M^{\text{th}}$ cluster with parameters $\Theta_{M}=\{\pi_{M},Q_{M},B_{M}\}$ into two clusters, with corresponding parameters $\Theta{{}^{\prime}}=\{\pi{{}^{\prime}},Q{{}^{\prime}},B{{}^{\prime}}\}$ and $\Theta{{}^{\prime\prime}}=\{\pi{{}^{\prime\prime}},Q{{}^{\prime\prime}},B{{}^{\prime\prime}}\}$. To implement the idea of centering proposals, we use a deterministic proposal for $Q$ and $\pi$, and let $Q{{}^{\prime}}=Q{{}^{\prime\prime}}=Q_{M}$ and $\pi^{{}^{\prime}}=\pi^{{}^{\prime\prime}}=\pi_{M}$. For observation model parameter $B$, we can use a similar proposal: for $k=1,\ldots,K$, let $\beta_{1,k}^{{}^{\prime}}=\beta_{M,1,k}$ and $\beta_{1,k}^{{}^{\prime\prime}}\sim\mathcal{N}\left(\beta_{M,1,k},c^{2}\right)$, with $\beta_{j,k}^{{}^{\prime}}=\beta_{j,k}^{{}^{\prime\prime}}=\beta_{M,j,k}$ for $j=2,\ldots,D$. For $\varpi$, let $w\sim Beta(2,2)$ and set $\varpi{{}^{\prime}}=w\varpi_{M}$ and $\varpi{{}^{\prime\prime}}=(1-w)\varpi_{M}$. The proposed move is from $M$ to $M+1$ clusters, with the new parameters $\Theta^{p}=\left\{\pi^{p},Q^{p},B^{p}\right\}$ with

[TABLE]

If we denote the posterior ratio as

[TABLE]

Then, the acceptance probability for this proposal is

[TABLE]

where $b_{M}$ is the probability of choosing the split move and $d_{M+1}=1-b_{M}$ is the probability of choosing the combine move, and $p_{\beta}(\cdot)$ is the Normal density with mean $\beta_{1,K}$ and variance $c^{2}$ and $p_{\varpi}(\cdot)$ is the $Beta(2,2)$ density.

4.2.2 Combine move

For the combine move, we consider an update from $M+1\to M$ clusters. Again, without loss of generality, we consider combine the $(M+1)^{\text{th}}$ and $M^{\text{th}}$ clusters into one cluster. We first find the stationary probabilities, $s_{M}$ and $s_{M+1}$ associated with $Q_{M}$ and $Q_{M+1}$ respectively. To combine $Q_{M}$ and $Q_{M+1}$ into $Q^{{}^{\prime}}$, the operation is as follows:

[TABLE]

and $q_{M,k,k}=-\sum_{i\neq k}q_{M,i,k}$ for $k=1,\ldots,K$. For the observation process parameter $B$,

[TABLE]

For the initial distribution $\pi$,

[TABLE]

and rescale the sum to 1. For mixture weight update, we set $\varpi^{{}^{\prime}}=\varpi_{M}+\varpi_{M+1}$. The acceptance probability for the move from $M+1$ to $M$ clusters is

[TABLE]

5 Dirichlet Process Mixture CTHMM

In the Dirichlet process mixture model, for some $\alpha>0$ and distribution $G_{0}$, we assume that $\Theta_{m}=\left\{\pi_{m},Q_{m},B_{m}\right\}\sim G_{0}$ for $m=1,2,\ldots$ again represent the component-specific model parameters. We then assume $\widetilde{\Theta}_{n}\sim G\left(.\right)$, and if $v_{m}\sim Beta(1,\alpha)$, $m=1,2,\ldots$

[TABLE]

We then suppose that cluster label $C_{n}$ is defined so that $C_{n}=m$ implies that $\widetilde{\Theta}_{n}=\Theta_{m}$. Within cluster $m$, for the latent process, $X_{n}|\widetilde{\Theta}_{n},C_{n}=m\sim\text{CTMC}(\pi_{m},Q_{m})$, and for the observations $O_{n}|{X_{n},\Theta_{n},C_{n}=m}\sim{\text{Exponential Family}}({{B_{m}}})$. We assume a priori for the equilibrium probabilities, $\pi\sim Dirichlet\left(\alpha_{1},\ldots,\alpha_{K}\right)$, and for the off-diagonal elements in $Q$, $q_{lm}\sim Gamma\left(a_{lm},b_{l}\right)$ for $1\leq l\neq m\leq K$; this prior is conjugate with the complete data likelihood – the representation using the latent trajectory as auxiliary data. Conjugacy can be relaxed using the approaches outlined in Neal (2000), but here we restrict attention to the conjugate case.

5.1 Split-merge MCMC algorithm for DMM

Inference for the Dirichlet process mixture model is often carried out using MCMC (Escobar and West, 1995; MacEachern and Müller, 1998; Neal, 2000; Jain and Neal, 2004, 2007). For this model, Gibbs sampling via a Pólya urn scheme is often used, but this can be slow mixing when clustering a large volume of data, and it is non-trivial to parallelize. This limitation has motivated the development of MCMC algorithms which partially address the inherently sequential nature of the updates. Celeux et al. (2000) noted that the standard MCMC sampler tends to stay within the neighborhood of the local mode. As a consequence, it is less likely to move to a new mixture component even with well-separated components because of the low probability of moving to an intermediate state: updating a group of subjects simultaneously can help resolve this problem. Green and Richardson (2001) introduced a split-merge update under the reversible-jump MCMC framework, and also showed how to construct the split proposal. Subsequently, Jain and Neal (2004, 2007) extended this to a Metropolis-Hastings (MH) sampling scheme with split-merge updates: the extended approach suggested splitting the component in a deterministic manner by employing restricted Gibbs sampling, which would increase the probability of forming a new component.

We will use the split-merge approach to obtain the posterior samples. We first focus on the case where the observation model is covariate-free and then the time-varying covariates are categorical only, and use the original parameterization $\theta$ in the exponential family as in (1) instead of $B$. The algorithm proceeds as follows, with superscript $i=0,1,\ldots$ denoting iteration number.

• •

Initialization: Randomly sample cluster labels $\mathbf{C}^{0}$ from $\{1,\dots,M^{0}\}$ for each subject, where $M^{0}$ is an arbitrary positive integer with $\left|\mathbf{C}^{0}\right|=M^{0}$. Starting with initial values (for example, those values can be obtained from the EM algorithm) $\pi^{0}_{C^{0}}$, $Q^{0}_{C^{0}}$, $\theta^{0}_{C^{0}}$ and $\phi^{0}_{C^{0}}$, compute the ‘forward’ and ‘backward’ quantities $a_{n,t,k}$ and $b_{n,t,k,j}$ using the forward-backward algorithm (see the Supplement for definitions of these quantities).

• •

Update latent state indicators: For each $n$ and $t$, generate the random vector $S_{n,t}^{i}$ from the multinomial distribution with the parameter set $a_{n,t}^{i}=(a_{n,t,1}^{i},\ldots,a_{n,t,K}^{i})$ where $S_{n,t}=\left(S_{n,t,1},\ldots,S_{n,t,K}\right)$ is an indicator random vector with $S_{n,t,k}=1$ if $X_{\tau_{t}}=k$ and 0 otherwise.

• •

Simulate the path of the latent process: For each individual, the path simulation follows in two steps. First, draw the current and next state ($X_{n,\tau_{n,t}}$ and $X_{n,\tau_{n,t+1}}$) from a multinomial distribution with the parameter matrix $b_{n,t,k,j}$, then simulate ${N_{n,l,m}}\left({{\Delta_{n,t}}}\right)$ and ${R_{n,l}}\left({{\Delta_{n,t}}}\right)$ from the Markov processes with infinitesimal generator $Q_{C_{n}}^{i-1}$ through the intervals $\left[\tau_{n,t},\tau_{n,t+1}\right)$ initiated at $X_{n,\tau_{n,t}}$ and end point $X_{n,\tau_{n,t+1}}$ sampled previously.

• •

Update $\mathbf{C}$ by split-merge (Jain and Neal, 2004): Details of this procedure can be found in the Supplement.

• •

**Update the component parameters $\theta$, $Q$ and $\pi$ **: Assume that $\mathbf{C}^{i}=\left\{C_{1}^{i},\ldots,C_{N}^{i}\right\}$.

Update $\mathbf{Q}$: Update the ${N_{l,m}}\left({{\Delta_{n,t}}}\right)$ and ${R_{l}}\left({{\Delta_{n,t}}}\right)$ from the updated label component generator $Q_{C}$, and then $q_{l,m\left|C\right.}^{i}$ associated with component $C$ from a Gamma distribution with shape parameter $\Lambda_{l,m}^{C}$ and rate parameter $\Upsilon_{l}^{C}$ where $\Lambda_{l,m}^{C}=\sum\limits_{n:C_{n}^{i}=C}{\sum\limits_{t=1}^{T_{n}}{{N_{l,m}}\left({{\Delta_{n,t}}}\right)}}+a_{lm}$ and $\Upsilon_{l}^{C}=\sum\limits_{n:C_{n}^{i}=C}{\sum\limits_{t=1}^{T_{n}}{{R_{l}}\left({{\Delta_{n,t}}}\right)}}+b_{l}$. 2. 2.

Update $\mathbf{\theta}$: For prior $\pi_{0}\left(\theta\right)$, generate the $\theta_{C}^{i}$ from the conditional posterior distribution using only the individuals with cluster label $C$ via a Gibbs/MH step. 3. 3.

Update $\mathbf{\pi}$: For conjugate $Dirichlet\left(\alpha_{1},\ldots,\alpha_{K}\right)$ prior,

[TABLE]

5.2 Acceptance probability for the split-merge proposal

The acceptance probability for the proposal to update $\mathbf{C}$ for split and merge moves takes the form

[TABLE]

The prior distribution of $\mathbf{C}$ is the product over the partition of $\left\{C_{1},\ldots,C_{N}\right\}$. Therefore

[TABLE]

where $N_{j}$ denotes the count of the number of subjects with label $j$ in the configuration. For MFMs, $\alpha$ and $N_{j}$ will be replaced by $\delta\frac{V_{N}\left(M^{*}+1\right)}{V_{N}\left(M^{*}\right)}$ and $N_{j}+\delta$, respectively.

As there is only one way to assign all $k\in S$ into one component, then the proposal density $q(\mathbf{C}^{\text{merge}}\left|\mathbf{C}\right.)=q(\mathbf{C}\left|\mathbf{C}^{\text{split}}\right.)=1$. For $q\left(\mathbf{C}^{\text{split}}\left|\mathbf{C}\right.\right)$, the probability is the product of transition probabilities from the last launch state $\mathbf{C}^{l}$ to the final proposed state $\mathbf{C}^{\text{split}}$, that is, $q(\mathbf{C}^{\text{split}}\left|\mathbf{C}\right.)=\prod_{k\in S}\mathbb{P}(C_{k}^{\text{split}}\left|C_{-k}\right.)$ which is calculated by the final Pólya urn scan in the split update . Each time, the $C_{k}$ is incrementally modified during the Pólya urn scan and the updated $C_{k}$ is used in the subsequent restricted Gibbs sampling computation.

Given the current labels for $f\in\mathcal{M}\cup\left\{d,e\right\}$, the likelihood contribution is given by

[TABLE]

where $H_{f,j}$ is the posterior distribution of $\Theta$ based on $G_{0}$ and all subjects such that their label $C_{g}=j$ but $g<f$. As parameters $\theta$, $\pi$ and $Q$ are separable in $\mathcal{L}_{f}\left(\Theta\right)$, $\mathcal{L}_{f}\left(\Theta\right)$ can be written as

[TABLE]

In the split step, the likelihood of the cluster label configuration is

[TABLE]

For the merge step, the likelihood of the cluster label configuration is

[TABLE]

The term $q\left(\mathbf{C}\left|\mathbf{C}^{\text{merge}}\right.\right)$ is the product of the transition probabilities from the last launch state to the original ‘split’ state, but there is no actual sampling step since the ‘split’ states are already known. Full details of the calculation of a Poisson example are given in the Supplement.

6 Examples

We examine the performance of the proposed algorithm in simulation. In all examples, data are generated from the finite mixture model with a fixed number of clusters. We place $Gamma(2.5,5)$ priors on the off-diagonal elements of $Q$, $\mathcal{N}\left(-2,1\right)$, $\mathcal{N}\left(0,1\right)$ and $\mathcal{N}\left(2,1\right)$ priors for states 1,2, and 3 respectively in the Gaussian case and $Gamma\left(5,10\right)$, $Gamma\left(10,10\right)$ and $Gamma\left(20,10\right)$ respectively in the Poisson case for the components of $B$, and a $Dirichlet\left(10,\ldots,10\right)$ prior for $\pi$. For MFMs, we take $\delta=1$ and $M\sim Poisson\left(0.5log(N)\right)+1$.

6.1 Example 1

In the first example, a three-cluster CTHMM-GLM is considered, with each cluster having three latent states and transition matrices

[TABLE]

with associated coefficient matrices

• •

Gaussian case: $B_{1}=\left(-4,0,5\right)$, $B_{2}=\left(-5.5,0.5,5.5\right)$, $B_{3}=\left(-5,1,4.8\right)$.

• •

Poisson case: $B_{1}=\left(-2,1.2,3\right)$, $B_{2}=\left(-1,1,2.5\right)$, $B_{3}=\left(-1.5,1.1,2.8\right)$.

The initial distributions for three clusters are $\pi_{1}=\left(0.5,0.4,0.1\right)^{\top}$, $\pi_{2}=\left(0.3,0.5,0.2\right)^{\top}$ and $\pi_{3}=\left(0.45,0.45,0.1\right)^{\top}$. In the Gaussian outcome model, the residual error standard deviation $\sigma$ is set to equal 1. In the split-merge Gibbs sampler, the intermediate restricted sampling scanned three times before performing the actual split or merge update. We initiate the model with one cluster. Data are generated by constructing the continuous Markov chain from the generator $Q_{i}$ for cluster $i=1,2,3$, a continuous-time realization $\left\{X_{s},0\leq s\leq 15\right\}$, and uniformly sampled $T-1$ time points between 0 and 15. Data are generated with 1000 subjects for each cluster.

Results are shown in Table 2 (an extended version is given in the Supplement), and are based on the three-cluster iterations taken from a total of 2000 iterations. MFM-RJ represents MFMs using reversible jump MCMC, and MFM-SM represents MFMs using split-merge algorithm. For Dirichlet mixture models (DMMs), in all cases, the posterior modal number of clusters was three, and this is especially clear when $T=50,100$ over 90% of iterations resided in the three-cluster model. For a summary output, cluster membership is assigned to the subject according to its posterior mode conditional on the three-cluster iterations. When $T=30$ for Normal case, there were 23.70% and 14.85% of two and four clusters respectively, while for the Poisson case, 43.40% of iterations contained four clusters. In general, restricted Gibbs sampling using split-merge proposals performed well; The norm differences in parameters between true values and the posterior means were small, and misclassification rates were low. It was more difficult to cluster trajectories, and the misclassification rate increased, and as sample size decreased, the Gaussian cases had higher misclassification than the Poisson cases. For MFMs, the posterior mode of the number of clusters is three for all cases. MFM-SM has a lower percentage of 3-cluster iterations as the split-merge algorithm encourages more trans-dimensional moves to explore the parameter space, while our reversible jump MCMC yields more accurate in terms of the cluster membership because it is updated in each sweep. In Miller and Harrison (2018), they included both the split-merge algorithm and the Pólya urn scheme in each MCMC step to encourage trans-dimensional moves as well as increase the accuracy of clustering. For a small sample size, this approach is feasible, but in our case when the sample size is large and the model structure is complex, it will be extremely computationally intensive to include both steps.

6.2 Example 2

In the second example, data are generated from a finite mixture model with different infinitesimal generators $Q$ within each cluster, but with other parameters $B$ and $\pi$ identical to Cluster 1 from Example 1. The observation process is taken to be the Gaussian distribution only. Data are generated with 1000 subjects in each cluster with 3000 subjects in total. The restricted Gibbs sampling procedure with split-merge proposals (with five intermediate Gibbs steps) for component parameters $Q$ only was implemented (see Supplement for the algorithm).

Results for this simulation are shown in Table 3 with an extended version in the Supplement: again, results presented are conditional on three-cluster iterations with total 2000 iterations. As this is a more difficult clustering problem, the percentage of three-cluster iterations decreased compared to the previous example. However, for DMMs, posterior modes of the number of clusters were still 3 apart from when $\sigma=2$, in which case the posterior was more diffuse; for the DMM also, the posterior distribution of the number of clusters was also more dispersed in this case. For $\sigma=0.5$, 17.05% of iterations had four clusters, while for $\sigma=1$, there were 4.10%. When $\sigma=2$, 18.15%, 38.25% and 6.05% of iterations had two, four and five clusters respectively. The misclassification rate for $\sigma=2$ is around 40%, but the posterior mean estimates were still close to the true values. While for MFMs, the posterior modes of the number of clusters were 3 for all cases. Specifically for MFM-RJ with $\sigma=0.5$, the model spent more iterations in two clusters before stabilizing while the model stabilized in three clusters faster, as for $\sigma=0.5$ cluster separation is greater. When $\sigma=2$, the algorithm tends to move to higher dimensions, with 20.95%, 19.80% 18.10% and 10.15% of iteration in four to seven clusters respectively. Similar to the previous example, MFM-RJ had the smallest misclassification rates conditional on three-cluster iterations, but MFM-SM yielded more trans-dimensional moves. This example illustrates that separating clusters via $Q$ is more challenging than for the other parameters.

6.3 Example 3

In the third example, again using a three cluster generating model, we add a three level, time-varying factor covariate $Z_{1}\sim{Multinomial}\left(1;1/4,1/4,1/2\right)$ at each observation time point to modify a Poisson outcome model, with coefficient matrices

[TABLE]

Parameters $Q$ and $\pi$ are as in Example 1. Data are generated with 300, 500, 200 subjects in Cluster 1, 2, and 3 respectively with 1000 subjects in total. The intermediate restricted Gibbs sampling scanned twice before performing the actual split or merge update. In this case, we place $Gamma\left(11,10\right)$, $Gamma\left(55,10\right)$ and $Gamma\left(165,10\right)$ for state 1,2,3 as prior distributions for $B$.

Results summarizing the three-cluster iterations for this example are shown in Table 4 (full version in the Supplement). We observed similar results in terms of misclassification rates for the Poisson cases in Example 1. For DMMs, as we have more parameters in this case, the posterior distribution of the number of clusters was again more dispersed (see plots in the Supplement). For $T=30$, the number of clusters fluctuated between three and four clusters, with 37.2% and 56.8% of total iterations respectively. When $T$ increases to 50 and 100, the mode becomes three with less posterior variation observed – for $T=50,100$, only 13.55% and 2.95% of 2000 iterations had four clusters respectively. For $T=100$, there were a few iterations with more than five clusters initially, but the number of clusters soon stabilized at three. For MFMs, all cases have posterior modal three. MFM-RJ cases tended to stay in one model for a longer time, which leads to more 3-cluster iterations. While for MFM-SM cases, they had more dimension changes but suffered from higher misclassification rates, which is consistent with previous examples.

7 Real example: Health surveillance of COPD patients

Our real health trajectories example relates to healthcare surveillance for the chronic condition COPD in greater Montreal, Québec, Canada. In 1998, a 25% random sample was drawn from the registry of the Régie de l’assurance maladie du Québec (RAMQ, the Québec provincial health authority) with a residential postal code in the census metropolitan area of Montreal. Subsequently at each year, a 25% random sample of residents new to Montreal within the previous year were sampled to maintain the representative cohort. If people died or changed their residential address outside of Montreal, the follow-up was ended. The data include outpatient diagnoses and procedures submitted through RAMQ billing claims, and procedures and diagnoses from inpatient claims. Using established case-definitions based on medical diagnostic codes (Blais et al., 2014; Lix et al., 2018), 76,888 COPD patients were enrolled with an incident event occurring after a minimum of two years at risk with no events. Patients were followed from January 1998, starting from the time of their first diagnosis, until December 2014. Physicians only observed these patients during medical visits, which occurred when patients chose to interact with the healthcare system, and at which information, including the number of prescribed medications, is collected. However, as this information was only available for patients with drug insurance, we restrict the cohort to patients over 65 years old with COPD, as prescription data are available for all of these patients.

In our analysis, we fit and compare MFMs and DMMs. Specifically, we use reversible-jump MCMC for MFMs, as from simulation studies it yielded more accurate results with a reasonable amount of trans-dimensional moves. The number of states in each cluster was fixed to four to match the convention in the COPD literature (GOLD Executive Committee, 2017), with the states termed mild, moderate, severe and very severe. We fit the model with parameters $Q$, $\theta_{\text{HOSP}}$, $\theta_{\text{SPEC}}$, $\theta_{\text{GP}}$, $\theta_{\text{ER}}$ and $\pi$ in each cluster, where $\theta_{U}$ ($U=$ HOSP, SPEC, GP, ER) represents the log mean parameter for the number of drugs prescribed in a Poisson model for each healthcare utilization. Each algorithm was initiated in the one-cluster model. Using elicitation via simulating trajectories, subjective prior distributions for off-diagonal elements of $Q$ is $Gamma\left(20,500\right)$, chosen with consideration for plausible holding and transition times within each state, and $Gamma\left(1.5,10\right)$, $Gamma\left(30,10\right)$, $Gamma\left(60,10\right)$ and $Gamma\left(100,10\right)$ for state 1, 2, 3, and 4 respectively for $B$ and $Dirichlet\left(10,\ldots,10\right)$ for $\pi$. Finally, we place $p_{0}(M)$ in Table 1 with $N=25000$ as the prior distribution on the number of clusters for the finite mixture model.

In the analysis to compare the performance of the two models and demonstrate the feasibility of clustering trajectories, we implemented this method on a data set comprising of 24,712 COPD patients. We present results for both models based on 2000 MCMC samples with 100 burn-in iterations after initialization using an EM algorithm fit of the one-cluster model as described in Luo et al. (2021). For both models, the mode of the posterior distribution of the number of clusters was four, and hence the results presented are from the four-cluster iterations. The Supplement gives details of the effective sample sizes of the posterior samples for $Q$ and $B$ in each cluster: these were satisfactory. Parameters $\theta_{\text{HOSP}}$, $\theta_{\text{SPEC}}$, $\theta_{\text{GP}}$ and $\theta_{\text{ER}}$ are converted to a contrast (relative-risk) parameterization in each cluster for both models. The posterior mean of the exponential of the coefficients for the four-cluster model with the number of patients in each cluster are shown in Tables in the Supplement. The numbers of drugs across different healthcare utilizations are different within each state. All four clusters have distinct numbers of prescribed drugs in each state. In general, specialist visits on average has smaller number of drugs prescribed compared to general practitioner visits; Cluster 1 has on average a greater number of drugs prescribed. Cluster 2 is the largest cluster, and has on average a greater number of drugs prescribed than the other two clusters. The posterior distribution for the Dirichlet mixture model varied slightly more than for finite mixture models as a result of the nature of the assumptions on cluster structure.

Recall that one of the principal computational challenges in this problem is that at each iteration we essentially need to sample and store complete trajectories for each of the almost 25,000 patients in the study; however, these posterior sampled trajectories can themselves be useful in gathering inference concerning time spent in each state, time to transition into more severe states and so on. In the Supplement, we show figures of posterior mean transition patterns over 5 years for the four clusters. In both plots, we observed similar transition patten; we summarize as follows: for State 1, patients in Cluster 2 are the most likely to hold in State 1, while Cluster 4 has the greatest probability for patients to exit this state; for State 2, all the clusters have a similar transition pattern; in State 3, patients in Cluster 2 have a relatively higher probability transition to State 4 compared to other clusters; Cluster 2 has notably different patterns from other clusters as patients tend to stay in State 4 but are more likely to transition to State 4 from State 3 – it is believed that those patients are likely in a progression from middle to end stage COPD.

Due to the data structure, with different start times and observations times for each subject, and the latent nature of the modeled process, it is challenging to demonstrate the clustering of individual trajectories by using conventional plotting of longitudinal curves so other summaries must suffice. The transition probabilities Figures discussed previously reveal some but not all of the differences between clusters. To examine further the distinction between clusters, Figure 1 display the real component of the three (non-trivial) eigenvalues of the posterior sampled $Q$ matrices in the four cluster realizations. In both plots, it clear that although there is some overlap for the second and third eigenvalues between Clusters 1 and 3, they are distinct. This evidence, coupled with the differences in outcome response levels indicated by the results in Tables of the posterior mean of the exponential of the coefficients – Cluster 1 has on average a greater number of drugs prescribed, indicating that patients in Cluster 1 have more severe health condition than those in Cluster 3 – confirms the presence of population substructure.

Table 5 displays the co-clustering matrix for two methods. The two methods are consistent in grouping the majority of patients into the same cluster. Clusters 1 and 3 showed some overlap in Figure 1, while the Dirichlet mixture model has a more distinct separation. This is because a high percentage of subjects in Cluster 1 from the Dirichlet mixture model is grouped in Cluster 3 in the finite mixture model, which brings the separation less distinct in the left panel of Figure 1. As discussed in the simulation example, the discrepancy in cluster memberships of two methods mainly comes from the updating mechanism of the number of clusters and the cluster membership of the two algorithms. Reversible-jump MCMC updates the cluster membership each iteration after the split/combine move of the number of clusters. However, split-merge MCMC updates both the cluster membership and the number of clusters simultaneously, which results in the cluster membership changing less often as reversible-jump MCMC.

8 Discussion

We implemented model-based clustering procedures for health trajectories based on continuous-time hidden Markov models based on finite and infinite mixture models. The methodology was applied to simulated examples, where the Markov transition rate matrices dictate that each cluster has its own transition characteristics and observation process. The posterior distribution for the cluster labels was computed by reversible-jump MCMC (Green, 1995) for finite mixture models and via Pólya urn schemes (Neal, 2000) and more efficient split-merge updates (Jain and Neal, 2004) for Dirichlet mixture models. Simulation studies demonstrated that both mixture models can identify the correct number of clusters and sample the target posterior distribution. Only the conjugate DP mixture case was explored in this paper; for handling the non-conjugate case, an auxiliary variable method (Neal (2000), Algorithm 8) or non-conjugate split-merge proposals (Jain and Neal, 2007) could be used. This would allow us to incorporate continuous time-varying covariates in the observation process and baseline covariates of general form in $Q$. Alternatively, if log-linear specifications are used in latent and observation models, an approximate conjugate analysis can be carried out using Gaussian approximations. In our illustrative analyses, we have examined a cohort of over 24,000 patients which comprised around one million records. This demonstrates the computing capability of our proposed algorithm. The main computational obstacle is the imputation of the latent process which is needed to facilitate all aspects of the likelihood calculation, and this may be in part overcome by parallelization of the individual likelihood calculations.

Focusing on the number of clusters, the Dirichlet precision parameter $\alpha$ and the prior distribution on the component specific parameters both influence the number of inferred clusters – this is inevitable in the context of the unsupervised learning problem to which model-based clustering corresponds. For example, a more precise prior on the components of $Q$ will encourage more clusters for the same observed data. Similarly, larger values of $\alpha$ will encourage more clusters; it is also straightforward to treat $\alpha$ as an unknown parameter in the MCMC algorithm, although we did not pursue that approach here.

In this paper, we have not explored an approach that allows the number of states to vary, but that is also possible. An efficient construction of proposal distributions is required in order to allow the sampler to move around and explore the parameter space for both the number of clusters and states (Brooks et al., 2003). We have extended our algorithm to allow the number of states to be inferred using the reversible jump approach, and will report on this elsewhere.

In Luo et al. (2021), other potential extensions or the model are discussed, and we recap them here. We have not addressed issues of informative dropout; in principle, dropout can be handled using standard Bayesian missing data procedures once a suitable missingness mechanism has been proposed, although the specification of a realistic model may be challenging to construct, and for the COPD data, the influence of informative dropout is likely to be minimal, as the principal cause of dropout is subjects leaving the province, which mainly affects younger subjects. Secondly, censoring due to death of the subject can be handled within the framework of the latent HMM by the inclusion of an absorbing state: however, information on time of death is not available in our real data set as the cohort contains only data on interactions with the healthcare system. This could be addressed by linking the RAMQ data to provincial death records. We assume here that the sequence of observation times $\tau_{t},t=1,\ldots,T$ is not informative about the underlying Markov chain or the outcome process, but the model could be extended to account for informative observation times by introducing a further stochastic process, for example an inhomogeneous Poisson process with rate dependent on the latent HMM – this would complicate the computation and be inhibitive for large data sets, but would still be feasible in principle.

Appendix A The CTHMM-GLM Model

This section outlines the formulation of Luo et al. (2021). To recap the model for a single individual, a sequence $\{O_{1},\ldots,O_{T}\}$ of variables is observed at time points $\{\tau_{1},\dots,\tau_{T}\}$. A latent process $\{X_{s}\}$ ($s\in\mathbb{R}^{+}$), representing the health status for a condition of interest, is assumed to be a continuous-time Markov chain (CTMC) with parameters $\left(\pi,Q\right)$, where $\pi$ is the initial distribution and $Q$ is the infinitesimal generator, taking values on the finite state space $\{1,2,\ldots,K\}$. The transition probability from state $i$ to $j$ in the time interval of length $\Delta_{t}=\tau_{t+1}-\tau_{t}$ between observation times $t$ and $t+1$ takes the form

[TABLE]

where $t=1,\ldots,T-1$, $Q=\left(q_{ij}\right)$ for $1\leq i,j\leq K$ is the infinitesimal generator for the continuous-time Markov process $\{X_{s}\}$, and $\text{expm}\left(A\right)$ is the matrix exponential of matrix $A$. The initial state distribution $\pi$ for $\{X_{s}\}$ is $\pi_{i}=\mathbb{P}\left(X_{0}=i\right)$ for $i=1,\ldots,K$. The model is parameterized by $\Theta=\left\{q_{ij},1\leq i,j\leq K,\pi,B,\phi\right\}$. The following diagram provides a schematic of the presumed data generating structure for one subject.

12321\tau_{1}=0$$\tau_{2}$$\tau_{3}$$\tau_{4}$$\tau_{5}$$\tau_{T}$$O_{1}$$O_{2}$$O_{3}$$O_{4}$$O_{5}$$O_{T}

$s$

$O_{s}\left|C\right.$

$X_{s}\left|C\right.$

In this diagram, outcome $O_{s}$ is observed at $s=\tau_{t},t=1,\ldots,T$. The underlying trajectory determines that the subject begins in state 1, then progresses through states 2,3,2 etc. until finally reverting to state 1 after the final observation. Observations are made at times that do not coincide with the transition times between states, and $\Delta_{t}=\tau_{t+1}-\tau_{t}$ records the interval between observations $O_{{t+1}}$ and $O_{t}$. $C$ represents the cluster label of this individual. In one cluster case, $C$ will be the same across all the subjects. In this paper, we assume that the measurement process itself (that is, the collection of times $\tau_{t},t=1,\ldots,T$) is not informative about the system either in its hidden or observed components. To illustrate the challenging nature of clustering trajectories, Figure 2 shows simulations from three groups with different parameters (taken from the Gaussian example in Example 1); Cluster 1 has in general smaller variation, but it is not easy to distinguish between Clusters 2 and 3. As differences among clusters are also coming from the underlying continuous-time Markov process $X_{s}$, it is necessary to extend the basic model to allow for heterogeneity in observation and latent processes.

Likelihood: Using standard arguments, the transition probability from state $i$ to $j$ in the time interval of length $\Delta_{t}=\tau_{t+1}-\tau_{t}$ between observation times $t$ and $t+1$ takes the form

[TABLE]

where $t=1,\ldots,T-1$, $Q=\left(q_{ij}\right)$ for $1\leq i,j\leq K$ is the infinitesimal generator for the continuous-time Markov process $\{X_{s}\}$ with $\sum_{j\neq i}q_{ij}=-q_{ii}>0$, and $\exp\left(A\right)$ is the matrix exponential of matrix $A$. The initial state distribution $\pi$ for $\{X_{s}\}$ is $\pi_{i}=\mathbb{P}\left(X_{0}=i\right)$ for $i=1,\ldots,K$. The model is parameterized by $\Theta=\left\{q_{ij},1\leq i,j\leq K,\pi,B,\phi\right\}$. For any individual, as $\{X_{s}\}$ is a Markov process with infinitesimal generator $Q$, if $\{X_{s}\}$ has been observed continuously in the time interval $\left[0,\tau\right]$, the likelihood function of $Q$ is

[TABLE]

where $N_{l,m}(\tau)$ is the number of transitions from state $l$ to state $m$ in the time interval $\left[0,\tau\right]$ and $R_{l}\left({{\tau}}\right)$ is the total time that the process has spent in state $l$ in $\left[0,\tau\right]$,

[TABLE]

Note that the quantities, $N_{l,m}\left(\tau\right)$ and $R_{l}\left({{\tau}}\right)$, are unobserved, but can be computed given a realization of the latent process on $[0,\tau]$.

For a random sample of $N$ subjects, let $O_{n,t}$ $\left(t=1,\ldots,T_{n}\right)$ be the $t^{\text{th}}$ observation for subject $n$ with the associated observation time $\tau_{n,t}$. The complete data likelihood derived from $\left\{O_{n}\right\}$ and $\left\{X_{n,\tau_{n}}\right\}$ can be factorized $\mathcal{L}(\Theta)\equiv\mathcal{L}(\mathbf{O},\mathbf{X}|\Theta)=\mathcal{L}(\mathbf{X}|\Theta)\mathcal{L}(\mathbf{O}|\mathbf{X},\Theta)$ where $\mathbf{O}=\{O_{n,t}\}$ and $\mathbf{X}=\{X_{n,\tau_{n,t}}\}$ for $n=1,\ldots,N,t=1,\ldots,T_{n}$ and

[TABLE]

The log-likelihood written in terms of the latent state indicator random vectors $\left\{S_{k}\right\}_{k=1}^{K}$ is

[TABLE]

where

[TABLE]

records the probability of transition from state $k$ to state $j$ in the interval $\Delta_{n,t}=\tau_{n,t+1}-\tau_{n,t}$. $N_{n,l,m}^{j,k}\left(\Delta_{n,t}\right)$ and $R_{n,l}^{j,k}\left(\Delta_{n,t}\right)$ are the amended versions of $N_{l,m}$ and $R_{l}$ computed conditional on starting in state $j$ and ending in state $k$ over the interval $\Delta_{n,t}$.

A.1 Definition of $a_{n,t,k}$ and $b_{n,t,k,j}$

For $1\leq k\leq K$, let

[TABLE]

say, where $b_{n,t,k,j}=\mathbb{P}(S_{n,t,k}=S_{n,{t+1},j}=1|O;\Theta)$. The $\left(a_{n,t,k},b_{n,t,k,j}\right)$ can be obtained by the forward-backward algorithm (Baum and Eagon, 1967; Baum and Sell, 1968). The forward variable $\alpha_{n,t,k}$ is defined as $\alpha_{n,t,k}=\mathbb{P}\left(O_{n,1},\ldots,O_{n,t},X_{n,\tau_{n,t}}=k\right)$. Let $O_{n}^{(t)}=(O_{n,1},\ldots,O_{n,t})$. Then

[TABLE]

The initial value $\alpha_{n,1,k}=\pi_{k}\times f\left(O_{n,1}\left|X_{n,\tau_{n,1}}=k\right.\right)$. The backward variable $\gamma_{n,t,k}$ is defined as $\gamma_{n,t,k}=\mathbb{P}\left(O_{n,T},\ldots,O_{n,t+1}\left|X_{\tau_{n,t}}=k\right.\right)$. Let $\overleftarrow{O}_{n}^{(t)}=(O_{n,T},\ldots,O_{n,t+1})$. Then

[TABLE]

The first backward value $\gamma_{n,T,k}$ is initialized to 1 for all $n$ and $k$. Then, for $t=1,\ldots,T_{n}-1$, $b_{n,t,k,j}=\mathbb{P}\left(S_{n,t,k}=S_{n,{t+1},j}=1\left|O\right.\right)$. Define $\overrightarrow{O}_{n}^{(t)}=(O_{n,t},\ldots,O_{n,T})$. Then

[TABLE]

with $a_{n,t,k}=\sum\limits_{j=1}^{K}{b_{n,t,k,j}}$.

For $t=T_{n}$,

[TABLE]

Appendix B Update $\mathbf{C}$ by split-merge by Jain and Neal (2004)

We layout the split-merge update for the cluster membership by Jain and Neal (2004) in this section. Denote by $M^{i-1}$ the number of components in the current label configuration. Select two distinct subjects $d$ and $e$ and denote their cluster labels $j_{d}=C_{d}^{i-1}$ and $j_{e}=C_{e}^{i-1}$.

Let $\mathcal{M}=\{f:f\neq d,e,\textrm{but }C_{f}^{i-1}=j_{d}\textrm{ or }C_{f}^{i-1}=j_{e}\}\subseteq\{1,2,\ldots,N\}$. 2. 2.

Define the launch state, $\mathbf{C}^{l}$ by running a Gibbs sampler scan restricted to the labels of subjects $f\in\mathcal{M}$.

1. (i)

If $j_{d}\neq j_{e}$: for $f\in\mathcal{M}$, $\mathbb{P}\left(C_{f}=j\left|C_{-f}\right.\right)$ for $j\in\{j_{d},j_{e}\}$ is given by

[TABLE]

where

• •

$N_{-f,j}=\sum\limits_{m\neq f}{\mathbbm{1}\left({{C_{m}}=j}\right)}$. For MFM, $N_{-f,j}$ will be replaced by $N_{-f,j}+\delta$;

• •

$\mathcal{L}_{f}\left(\Theta\right)$ is defined in the main paper Section 5.2;

• •

$H_{-f,j}\left(\Theta\right)$ is the posterior distribution of $\Theta$ based on the prior $G_{0}$ and subjects $g\in\mathcal{M}\cup\left\{d,e\right\}$ but $g\neq f$, such that their labels $C_{g}=j$. 2. (ii)

If $j_{d}=j_{e}$: for $f\in\mathcal{M}$, the conditional probability $\mathbb{P}\left(C_{f}=j\left|C_{-f}\right.\right)$ for $j\in\{j_{e},M^{i-1}+1\}$ is calculated in the same fashion in (13). The value $M^{i-1}+1$ is a new component label, not represented in the current label set. 3. 3.

Split-merge step:

1. (i)

If $j_{d}=j_{e}$, propose a split configuration, $\mathbf{C}^{\text{split}}$.

• •

Set $C_{d}^{\text{split}}=M^{i-1}+1$ and $C_{e}^{\text{split}}=C_{e}^{i-1}=j_{e}$.

• •

For $f\notin\mathcal{M}$, let $C_{f}^{\text{split}}=C_{f}^{i-1}$; for $f\in\mathcal{M}$, set $C_{f}^{\text{split}}$ by performing one more Pólya urn scan from the launch state.

• •

Perform a MH update with acceptance probability $a\left(\mathbf{C}^{\text{split}},\mathbf{C}^{i-1}\right)$. If the proposal is accepted, then set $\mathbf{C}^{i}=\mathbf{C}^{\text{split}}$; if rejected, set $\mathbf{C}^{i}=\mathbf{C}^{i-1}$. 2. (ii)

If $j_{d}\neq j_{e}$, propose the merge configuration, $\mathbf{C}^{\text{merge}}$.

• •

Set $C_{d}^{\text{merge}}=C_{e}^{\text{merge}}=C_{e}^{i-1}=j_{e}$.

• •

For $f\in\mathcal{M}$, set $C_{f}^{\text{merge}}=j_{e}$; for $f\notin\mathcal{M}$, let $C_{f}^{\text{merge}}=C_{f}^{i-1}$.

• •

Perform a MH update with acceptance probability $a\left(\mathbf{C}^{\text{merge}},\mathbf{C}^{i-1}\right)$. If the proposal is accepted, then set $\mathbf{C}^{i}=\mathbf{C}^{\text{merge}}$; if rejected, set $\mathbf{C}^{i}=\mathbf{C}^{i-1}$.

The acceptance probability for this proposal is computed in the following subsection.

Appendix C Split-Merge Proposals for the Conjugate Dirichlet Process Mixture Model with Component Parameter $Q$

We now propose an algorithm to cluster individual trajectories via the DP mixture for the CTHMM-GLM with component parameters $Q$.

• •

Initialization: Randomly sample cluster label $\mathbf{C}^{0}$ from $\{1,\dots,M^{0}\}$ for each subject, where $M^{0}$ is an arbitrary positive integer with $\left|\mathbf{C}^{0}\right|=M^{0}$. Starting with initial values $\pi^{0}$, $Q^{0}_{C^{0}}$, $\theta^{0}$ and $\phi^{0}$, compute $a_{n,t,k}$ and $b_{n,t,k,j}$ using the forward-backward algorithm.

• •

Update latent state indicators: For each $n$ and $t$, generate the random vector $S_{n,t}^{i}$ from the multinomial distribution with the parameter set $a_{n,t}^{i}=(a_{n,t,1}^{i},\ldots,a_{n,t,K}^{i})$ where $S_{n,t}=\left(S_{n,t,1},\ldots,S_{n,t,K}\right)^{\top}$ is an indicator random vector with $S_{n,t,k}=1$ if $X_{\tau_{t}}=k$ and 0 otherwise.

• •

Update $B$ and $\phi$: Sample coefficient matrix $B^{i}$ and scale parameter $\phi^{i}$ given $S_{n,t}^{i}$ via the MH algorithm as their conditional posterior distributions are not of standard form.

• •

Update $\pi$: Under prior $Dirichlet\left(\alpha_{1},\ldots,\alpha_{K}\right)$, sample $\pi^{i}$ from a Dirichlet distribution with parameters

[TABLE]

• •

Simulate the path for the latent process:

• –

Sample the current state and next state ($X_{n,\tau_{n,t}}$ and $X_{n,\tau_{n,t+1}}$) from a multinomial distribution with the parameter matrix $b_{n,t,k,j}$.

• –

Simulate ${N_{n,l,m}}\left({{\Delta_{n,t}}}\right)$ and ${R_{n,l}}\left({{\Delta_{n,t}}}\right)$ from the Markov jump processes step-by-step with infinitesimal generator $Q_{C_{n}}^{i-1}$ through the intervals $\left[\tau_{n,t},\tau_{n,t+1}\right)$ initiated at $X_{n,\tau_{n,t}}$ and end point $X_{n,\tau_{n,t+1}}$ sampled previously with current label $C^{i-1}$.

• •

Update label $\mathbf{C}$ by split-merge: Denote by $M^{i-1}$ the number of clusters in the current label configuration. Select two distinct subjects $d$ and $e$ and denote their cluster labels $j_{d}=C_{d}^{i-1}$ and $j_{e}=C_{e}^{i-1}$.

Let $\mathcal{M}=\{f:f\neq d,e,\textrm{but }C_{f}^{i-1}=j_{d}\textrm{ or }C_{f}^{i-1}=j_{e}\}\subseteq\{1,2,\ldots,N\}$. 2. 2.

Define the launch state, $\mathbf{C}^{l}$ by running a Gibbs sampler scan restricted to the labels of subjects $f\in\mathcal{M}$.

1. (i)

If $j_{d}\neq j_{e}$: for $f\in\mathcal{M}$, the conditional probability $\mathbb{P}\left(C_{f}=j\left|C_{-f}\right.\right)$ for $j\in\{j_{d},j_{e}\}$ is given by

[TABLE]

where

[TABLE]

and where $H_{-f,j}\left(q_{lm}\right)$ denotes the posterior distribution of $q_{lm}$ based on the prior $G_{0}$ and subjects $g\in\mathcal{M}\cup\left\{d,e\right\}$ but $g\neq f$ such that $C_{g}=j$. 2. (ii)

If $j_{d}=j_{e}$: for $f\in\mathcal{M}$, the conditional probability $\mathbb{P}\left(C_{f}=j\left|C_{-f}\right.\right)$ for $j\in\{j_{e},M^{i-1}+1\}$ is calculated in the same fashion. The value $M^{i-1}+1$ is a new component label, not represented in the current label set. 3. 3.

Split-merge step:

1. (i)

If subjects $d$ and $e$ are in the same mixture component, i.e., $C_{d}^{i-1}=C_{e}^{i-1}$, propose the split procedure, $\mathbf{C}^{\text{split}}$.

• –

Set $C_{d}^{\text{split}}=M^{i-1}+1$ and $C_{e}^{\text{split}}=C_{e}^{i-1}=j_{e}$.

• –

For $f\notin\mathcal{M}$, let $C_{f}^{\text{split}}=C_{f}^{i-1}$; for $f\in\mathcal{M}$, modify $C_{f}^{\text{split}}$ by performing one more Pólya urn scan from the launch state label $C^{l}_{f}$.

• –

Perform the MH update with acceptance probability $a\left(\mathbf{C}^{\text{split}},\mathbf{C}^{i-1}\right)$. If the proposal is accepted, then set $\mathbf{C}^{i}=\mathbf{C}^{\text{split}}$; if rejected, set $\mathbf{C}^{i}=\mathbf{C}^{i-1}$. 2. (ii)

If subjects $d$ and $e$ are in different mixture components, i.e., $C_{d}^{i-1}\neq C_{e}^{i-1}$, propose the merge procedure, $\mathbf{C}^{\text{merge}}$.

• –

Set $C_{d}^{\text{merge}}=C_{e}^{\text{merge}}=C_{e}^{i-1}=j_{e}$.

• –

For $f\in\mathcal{M}$, let $C_{f}^{\text{merge}}=C_{e}^{i-1}$; for $f\notin\mathcal{M}$, let $C_{f}^{\text{merge}}=C_{f}^{i-1}$.

• –

Perform the MH update with acceptance probability $a\left(\mathbf{C}^{\text{merge}},\mathbf{C}^{i-1}\right)$. If the proposal is accepted, then set $\mathbf{C}^{i}=\mathbf{C}^{\text{merge}}$; if rejected, set $\mathbf{C}^{i}=\mathbf{C}^{i-1}$.

• •

Update the component parameter $Q$: For all $C$ in $\mathbf{C}^{i}=\left\{C_{1}^{i},\ldots,C_{N}^{i}\right\}$, update ${N_{l,m}}\left({{\Delta_{n,t}}}\right)$ and ${R_{l}}\left({{\Delta_{n,t}}}\right)$ from the updated label component generator $Q_{C}$ (obtained when calculating $a\left(\mathbf{C}^{*},\mathbf{C}^{i-1}\right)$). Then $q_{l,m\left|C\right.}^{i}$ associated with component $C$ is generated from a Gamma distribution with shape parameter $\Lambda_{l,m}^{C}$ and rate parameter $\Upsilon_{l}^{C}$.

In terms of calculating the acceptance probability, the procedure is the same as the general algorithm except for the likelihood ratio. In the general algorithm, the likelihood of the component parameter $Q$ is based on the unobserved Markov process, namely, ${N_{l,m}}\left({{\Delta_{n,t}}}\right)$ and ${R_{l}}\left({{\Delta_{n,t}}}\right)$. Therefore, when calculating the likelihood ratio for split-merge procedures, the unobserved Markov process should be modified as well. The likelihood of the current label for $f\in\mathcal{M}\cup\left\{d,e\right\}$ is

[TABLE]

where ${N_{l,m}}\left({{\Delta_{f,t}}}\right)$ and ${R_{l}}\left({{\Delta_{f,t}}}\right)$ are generated from the previous step. In the split step, the likelihood of the split label is calculated by

[TABLE]

where

[TABLE]

and ${N^{e,\text{split}}_{l,m}}\left({{\Delta_{f,t}}}\right)$, ${R^{e,\text{split}}_{l}}\left({{\Delta_{f,t}}}\right)$ are simulated from the Markov jump processes step-by-step with infinitesimal generator $Q^{e,\text{split}}$ through the intervals $\left[\tau_{n,t},\tau_{n,t+1}\right)$ initiated at $X_{\tau_{n,t}}$ and end point $X_{\tau_{n,t+1}}$ respectively and each entry in $Q^{e,\text{split}}$ is generated from a Gamma distribution with shape parameter $\Lambda_{l,m}^{j_{e}}$ and rate parameter $\Upsilon_{l}^{j_{e}}$. Similarly, ${N^{d,\text{split}}_{l,m}}\left({{\Delta_{f,t}}}\right),{R^{d,\text{split}}_{l}}\left({{\Delta_{f,t}}}\right)$ are sampled associated with $Q^{d,\text{split}}$ where each entry in $Q^{d,\text{split}}$ is sampled from a Gamma distribution with shape parameter $\Lambda_{l,m}^{C_{d}^{\text{split}}}$ and rate parameter $\Upsilon_{l}^{C_{d}^{\text{split}}}$.

If the split step is accepted, $Q^{d,\text{split}}$ and $Q^{e,\text{split}}$ will be carried into the next step to update $Q$. Similarly in the merge step, the likelihood of the merge label is calculated by

[TABLE]

where $\mathcal{L}^{e,\text{merge}}\left(q_{lm}\left|\Delta_{f,t}\right.\right)={{q_{l,m}}^{{N^{e,\text{merge}}_{l,m}}\left({{\Delta_{f,t}}}\right)}e^{\left({-{q_{l,m}}{R^{e,\text{merge}}_{l}}\left({{\Delta_{n,t}}}\right)}\right)}}$ and ${N^{e,\text{merge}}_{l,m}}\left({{\Delta_{f,t}}}\right)$, ${R^{e,\text{merge}}_{l}}\left({{\Delta_{f,t}}}\right)$ are simulated associated with $Q^{e,\text{merge}}$ where each entry in $Q^{e,\text{merge}}$ is sampled from a Gamma distribution with shape parameter $\Lambda_{l,m}^{j_{e}}$ and rate parameter $\Upsilon_{l}^{j_{e}}$.

If the merge step is accepted, $Q^{e,\text{merge}}$ will be carried into the next step when updating $Q$.

Appendix D Split-merge illustration: Poisson model

Suppose $O_{t}\left|X_{\tau_{t}}=k\right.\sim\text{Poisson}\left(\theta_{k}\right)$ and that there are no covariates in outcome model or latent model. Data $\mathbf{o}=\left\{o_{n,t}\right\}$ for $n=1,\ldots,N,t=1,\ldots,n$, have the complete data log-likelihood

[TABLE]

where $p_{n,t}^{k,j}$ is the probability of transition from state $k$ to state $j$, defined in (12). If a priori $\pi\sim Dirichlet\left(\alpha_{1},\ldots,\alpha_{K}\right)$, $q_{lm}\sim Gamma\left(a_{lm},b_{l}\right)$ for $1\leq l\neq m\leq K$ and $\theta_{k}\sim Gamma\left(a_{k},b_{k}\right)$ independently, then we have a conjugate model, and the integrals defining the split and merge probabilities are analytically tractable; we can integrate out the model parameters separately. To update $C_{f}$ via intermediate Gibbs sampling, a new value of $C_{f}$ is drawn from (13).

**Outcome model parameters: ** For $\theta_{k}$ ($k=1,\ldots,K$), given latent state indicators, since this is a conjugate Poisson-Gamma model, $H_{-f,C}\left(\theta_{k}\right)\sim Gamma(a_{k}^{\prime},b_{k}^{\prime})$ where

[TABLE]

Therefore, as

[TABLE]

we have that

[TABLE]

Initial state parameters: Within each cluster, $H_{-f,j}\left(\pi\right)\sim Dirichlet\left(\alpha_{1}^{\prime},\ldots,\alpha_{K}^{\prime}\right)$, where $\alpha_{k}^{\prime}=\alpha_{k}+\sum\limits_{n\neq f:C_{n}=j}S_{n,1,k}$. Then,

[TABLE]

Transition model parameters: For the components $q_{lm}$ of $Q$, given $N_{lm}\left(\tau\right)$ and $R_{l}\left(\tau\right)$ as defined previously, $H_{-f,C}\left(q_{lm}\right)\sim Gamma(a_{lm}^{\prime},b_{lm}^{\prime})$, where

[TABLE]

Then,

[TABLE]

Therefore,

[TABLE]

and the probability in (13) readily computed.

The MH acceptance probability also involves calculating the marginal likelihood of the cluster membership, $\mathcal{L}\left(\mathbf{C}\right)$ defined in the main paper Section 5.2. The integrals can be calculated in the same fashion above and the only difference is to change $H_{-f,j}$ to $H_{f,j}$. Specifically,

[TABLE]

A similar conjugate approach can be used if factor predictors are incorporated into the models for the outcome or the transition infinitesimal generator.

Appendix E Example Results

The following results demonstrate that the algorithm successfully recovers the simulating parameters in the posterior distributions. In each case, the posterior mean was used to calculate the norm difference from the true values.

E.1 Example 1

E.2 Example 2

E.3 Example 3

Appendix F Real example: Health surveillance of COPD patients

The following figures and tables demonstrate adequacy of the MCMC approaches for the Bayesian posterior calculations for the real data.

24,712 COPD Patients: Finite mixture model analysis, effective sample sizes (out of 1636) from the MCMC analysis.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

24,712 COPD Patients: DP mixture model analysis, effective sample sizes (out of 1894) from the MCMC analysis.

[TABLE]

[TABLE]

[TABLE]

[TABLE]