Mixture of hidden Markov models for accelerometer data

TL;DR

This paper introduces a finite mixture of hidden Markov models tailored for accelerometer data, enabling characterization of physical activity patterns and subpopulation identification without predefined activity levels.

Contribution

It presents a novel mixture HMM approach that estimates activity levels from data and accounts for population heterogeneity, with theoretical guarantees on classification accuracy.

Findings

Misclassification probability decreases exponentially with data length.

Model effectively captures activity transition dynamics.

Numerical simulations support theoretical results.

Abstract

Motivated by the analysis of accelerometer data, we introduce a specific finite mixture of hidden Markov models with particular characteristics that adapt well to the specific nature of this type of data. Our model allows for the computation of statistics that characterize the physical activity of a subject (\emph{e.g.}, the mean time spent at different activity levels and the probability of the transition between two activity levels) without specifying the activity levels in advance but by estimating them from the data. In addition, this approach allows the heterogeneity of the population to be taken into account and subpopulations with homogeneous physical activity behavior to be defined. We prove that, under mild assumptions, this model implies that the probability of misclassifying a subject decreases at an exponential decay with the length of its measurement sequence. Model identifiability is also investigated. We also report a comprehensive suite of numerical simulations to support our theoretical findings. Method is motivated by and applied to the PAT study.

Tables

Table 1. Table 1: Convergence of estimators when 1000 replicates are drawn from case medium: ARI between estimated and true partition, ARI between estimated and true latent states and MSE between the MLE and the true parameters

		ARI (latent variables)		MSE (model parameters)
$n$	$T$	partition	states	${𝑨}_k$	${\epsilon }_h$	$a_h$	$b_h$	${\delta }_k$
10	100	0.995	0.621	0.021	0.001	0.088	0.024	0.047
10	500	1.000	0.632	0.007	0.000	0.020	0.005	0.048
100	100	0.996	0.630	0.004	0.000	0.011	0.003	0.005
100	500	1.000	0.634	0.003	0.000	0.005	0.002	0.005

Table 2. Table 2: Convergence of estimators obtained over 1000 replicates with and without missing data when data are sampled from case medium: ARI between estimated and true partition, ARI between estimated and true latent states and MSE between the MLE and the true parameters

			Adjusted Rand index		Mean square error
$n$	$T$	missingness	partition	states	${𝑨}_k$	${\epsilon }_h$	$a_h$	$b_h$	${\delta }_k$
10	100	no missingness	0.995	0.621	0.021	0.001	0.088	0.024	0.047
		MCAR-1	0.991	0.613	0.024	0.001	0.102	0.028	0.047
		MCAR-2	0.987	0.605	0.028	0.001	0.113	0.032	0.047
		MNAR	0.934	0.497	0.051	0.003	0.398	0.050	0.050
10	500	no missingness	1.000	0.632	0.007	0.000	0.020	0.005	0.048
		MCAR-1	1.000	0.631	0.007	0.000	0.020	0.005	0.048
		MCAR-2	1.000	0.631	0.007	0.000	0.019	0.005	0.048
		MNAR	0.999	0.516	0.021	0.003	0.233	0.028	0.048
100	100	no missingness	0.996	0.630	0.004	0.000	0.011	0.003	0.005
		MCAR-1	0.994	0.624	0.004	0.000	0.013	0.003	0.005
		MCAR-2	0.989	0.618	0.005	0.000	0.014	0.004	0.005
		MNAR	0.951	0.512	0.014	0.002	0.200	0.026	0.005
100	500	no missingness	1.000	0.634	0.003	0.000	0.005	0.002	0.005
		MCAR-1	1.000	0.633	0.002	0.000	0.006	0.002	0.005
		MCAR-2	1.000	0.632	0.002	0.000	0.005	0.002	0.005
		MNAR	1.000	0.520	0.011	0.002	0.198	0.026	0.005

Table 3. Table 3: Parameters and mean time per states for the three subjects

State name	${\epsilon }_h$	$a_h$	$b_h$	mean	sd
intensive-level	0.00	98.94	0.65	152.76	15.36
moderate-level	0.00	11.09	0.11	99.34	29.84
low-level	0.00	2.32	0.11	20.98	13.79
sleeping	0.22	1.48	0.72	2.06	1.70

Table 4. Table 4: Information criteria obtained on PAT data with four levels of activity (minima are in bold)

$K$	1	2	3	4	5	6	7
BIC	-2953933	-2952313	-2951809	-2951705	-2951308	-2951364	-2951696
ICL	-2953933	-2952313	-2951810	-2951707	-2951309	-2951364	-2951697

Table 5. Table 5: Parameters describing the four activity levels for PAT data and statistics of the distributions

Name of the activity level	Parameters			Statistics
	${\epsilon }_h$	$a_h$	$b_h$	mean	marginal cutoffs
sleeping	0.988	7.470	7.470	0.012	$[0,\mathrm{\hspace{0.33em}0}]$
low-level	0.260	0.974	0.020	36.926	$]0,\mathrm{\hspace{0.33em}97.7}]$
moderate-level	0.025	1.408	0.004	329.249	$]97.7,\mathrm{\hspace{0.33em}614.4}]$
intensive-level	0.007	2.672	0.002	1696.935	$]614.4,+\infty [$

Table 6. Table 6: Mean time spent at the different activity levels for the five classes

	sleeping	low-level	moderate-level	intensive-level
Class active	0.306	0.284	0.338	0.072
Class sedentary	0.467	0.209	0.263	0.061
Class moderate	0.304	0.411	0.225	0.060
Class very sedentary	0.504	0.366	0.124	0.006
Class very active	0.189	0.351	0.316	0.143

Table 7. Table 7: Transition matrix for the five classes

	Class moderate
	sleeping	low-level	moderate-level	intensive-level
sleeping	0.76	0.21	0.03	0.00
low-level	0.16	0.73	0.11	0.00
moderate-level	0.03	0.20	0.73	0.04
intensive-level	0.01	0.04	0.16	0.80
	Class very sedentary
	sleeping	low-level	moderate-level	intensive-level
sleeping	0.85	0.08	0.06	0.00
low-level	0.20	0.67	0.13	0.01
moderate-level	0.10	0.11	0.76	0.03
intensive-level	0.01	0.04	0.14	0.82
	Class very active
	sleeping	low-level	moderate-level	intensive-level
sleeping	0.80	0.14	0.05	0.01
low-level	0.08	0.74	0.17	0.01
moderate-level	0.03	0.18	0.69	0.10
intensive-level	0.01	0.05	0.21	0.74

Table 8. Table 8: Convergence of the estimators with 1000 replicates: ARI between estimated and true partition, ARI between estimated and true latent states and MSE between the MLE and the true parameters

		Adjusted Rand index		Mean square error
$n$	$T$	partition	states	${𝑨}_k$	${\epsilon }_h$	$a_h$	$b_h$	${\delta }_k$
Case hard ($e=0.75$ and $a_2=3$)
10	100	0.812	0.350	0.097	0.002	0.285	0.085	0.057
10	500	1.000	0.385	0.031	0.000	0.077	0.024	0.052
100	100	0.870	0.377	0.032	0.000	0.087	0.028	0.007
100	500	1.000	0.391	0.018	0.000	0.045	0.017	0.005
Case medium-easy ($e=0.75$ and $a_2=5$)
10	100	0.996	0.661	0.018	0.001	0.253	0.031	0.052
10	500	0.999	0.669	0.004	0.000	0.050	0.006	0.052
100	100	0.998	0.669	0.002	0.000	0.027	0.003	0.005
100	500	1.000	0.671	0.000	0.000	0.005	0.001	0.005
Case easy ($e=0.90$ and $a_2=5$)
10	100	1.000	0.832	0.006	0.000	0.159	0.017	0.048
10	500	1.000	0.839	0.001	0.000	0.030	0.003	0.046
100	100	1.000	0.836	0.001	0.000	0.015	0.002	0.005
100	500	1.000	0.839	0.000	0.000	0.003	0.000	0.005

Table 9. Table 9: Convergence of the estimators obtained on 1000 replicates with and without missing data when data are sampled from case hard: ARI between estimated and true partition, ARI between estimated and true latent states and MSE between the MLE and the true parameters

			Adjusted Rand index		Mean square error
$n$	$T$	missingness	partition	states	${𝑨}_k$	${\epsilon }_h$	$a_h$	$b_h$	${\delta }_k$
10	100	no missingness	0.812	0.350	0.097	0.002	0.285	0.085	0.057
		MCAR-1	0.753	0.342	0.103	0.002	0.360	0.092	0.061
		MCAR-2	0.690	0.333	0.120	0.003	0.320	0.091	0.068
		MNAR	0.542	0.303	0.152	0.004	0.456	0.109	0.080
10	500	no missingness	1.000	0.385	0.031	0.000	0.077	0.024	0.052
		MCAR-1	1.000	0.386	0.029	0.000	0.072	0.024	0.052
		MCAR-2	0.999	0.385	0.032	0.000	0.074	0.024	0.052
		MNAR	0.991	0.347	0.042	0.003	0.177	0.076	0.052
100	100	no missingness	0.870	0.377	0.032	0.000	0.087	0.028	0.007
		MCAR-1	0.837	0.373	0.033	0.000	0.091	0.030	0.008
		MCAR-2	0.801	0.369	0.035	0.000	0.099	0.032	0.008
		MNAR	0.664	0.337	0.038	0.003	0.099	0.083	0.011
100	500	no missingness	1.000	0.391	0.018	0.000	0.045	0.017	0.005
		MCAR-1	1.000	0.391	0.018	0.000	0.050	0.017	0.005
		MCAR-2	1.000	0.391	0.018	0.000	0.044	0.017	0.005
		MNAR	0.996	0.355	0.023	0.003	0.109	0.073	0.005

Table 10. Table 10: Convergence of the estimators obtained on 1000 replicates with and without missing data when data are sampled from case medium-easy: ARI between estimated and true partition, ARI between estimated and true latent states and MSE between the MLE and the true parameters

			Adjusted Rand index		Mean square error
$n$	$T$	missingness	partition	states	${𝑨}_k$	${\epsilon }_h$	$a_h$	$b_h$	${\delta }_k$
10	100	no missingness	0.996	0.661	0.018	0.001	0.253	0.031	0.052
		MCAR-1	0.993	0.659	0.020	0.001	0.276	0.034	0.052
		MCAR-2	0.986	0.657	0.023	0.001	0.308	0.038	0.051
		MNAR	0.966	0.621	0.052	0.004	0.561	0.100	0.052
10	500	no missingness	0.999	0.669	0.004	0.000	0.050	0.006	0.052
		MCAR-1	0.999	0.669	0.004	0.000	0.051	0.006	0.052
		MCAR-2	0.999	0.668	0.004	0.000	0.052	0.007	0.052
		MNAR	1.000	0.636	0.021	0.003	0.217	0.085	0.052
100	100	no missingness	0.998	0.669	0.002	0.000	0.027	0.003	0.005
		MCAR-1	0.996	0.668	0.002	0.000	0.030	0.003	0.005
		MCAR-2	0.993	0.667	0.002	0.000	0.033	0.004	0.005
		MNAR	0.978	0.641	0.013	0.003	0.160	0.074	0.005
100	500	no missingness	1.000	0.671	0.000	0.000	0.005	0.001	0.005
		MCAR-1	1.000	0.671	0.000	0.000	0.005	0.001	0.005
		MCAR-2	1.000	0.671	0.000	0.000	0.005	0.001	0.005
		MNAR	1.000	0.644	0.011	0.003	0.140	0.076	0.005

Table 11. Table 11: Convergence of the estimators obtained on 1000 replicates with and without missing data when data are sampled from case easy: ARI between estimated and true partition, ARI between estimated and true latent states and MSE between the MLE and the true parameters

			Adjusted Rand index		Mean square error
$n$	$T$	missingness	partition	states	${𝑨}_k$	${\epsilon }_h$	$a_h$	$b_h$	${\delta }_k$
10	100	no missingness	1.000	0.832	0.006	0.000	0.159	0.017	0.048
		MCAR-1	1.000	0.829	0.007	0.000	0.185	0.019	0.048
		MCAR-2	1.000	0.826	0.008	0.001	0.207	0.022	0.048
		MNAR	1.000	0.771	0.017	0.003	0.315	0.042	0.048
10	500	no missingness	1.000	0.839	0.001	0.000	0.030	0.003	0.046
		MCAR-1	1.000	0.838	0.001	0.000	0.031	0.003	0.046
		MCAR-2	1.000	0.837	0.001	0.000	0.031	0.003	0.046
		MNAR	1.000	0.781	0.007	0.003	0.119	0.030	0.046
100	100	no missingness	1.000	0.836	0.001	0.000	0.015	0.002	0.005
		MCAR-1	1.000	0.833	0.001	0.000	0.018	0.002	0.005
		MCAR-2	1.000	0.831	0.001	0.000	0.019	0.002	0.005
		MNAR	1.000	0.779	0.006	0.003	0.084	0.026	0.005
100	500	no missingness	1.000	0.839	0.000	0.000	0.003	0.000	0.005
		MCAR-1	1.000	0.838	0.000	0.000	0.003	0.000	0.005
		MCAR-2	1.000	0.838	0.000	0.000	0.003	0.000	0.005
		MNAR	1.000	0.783	0.005	0.003	0.072	0.025	0.005

Table 12. Table 12: Transition matrix for the two remaining classes

	Class active
	sleeping	low-level	moderate-level	intensive-level
sleeping	0.87	0.12	0.01	0.00
low-level	0.17	0.73	0.10	0.00
moderate-level	0.04	0.30	0.66	0.01
intensive-level	0.08	0.08	0.18	0.66
	Class sedentary
	sleeping	low-level	moderate-level	intensive-level
sleeping	0.79	0.16	0.05	0.00
low-level	0.17	0.66	0.16	0.01
moderate-level	0.05	0.14	0.79	0.03
intensive-level	0.01	0.02	0.15	0.82

Equations

$$p({\mathbf{y}}_i;\mathbf{\theta })=\sum _{k=1}^K{\delta }_{k}p({\mathbf{y}}_i;{\mathbf{\pi }}_k,{\mathbf{A}}_k,\mathbf{\lambda },\mathbf{\epsilon }),$$

$$p({\mathbf{y}}_i;{\mathbf{\pi }}_k,{\mathbf{A}}_k,\mathbf{\lambda },\mathbf{\epsilon })=\sum _{{\mathbf{x}}_i\in \mathcal{X}}p({\mathbf{x}}_i;{\mathbf{\pi }}_k,{\mathbf{A}}_k)p({\mathbf{y}}_i\mid {\mathbf{x}}_i;\mathbf{\lambda },\mathbf{\epsilon }).$$

$$p({\mathbf{x}}_i;{\mathbf{\pi }}_k,{\mathbf{A}}_k)=\prod _{h=1}^{\ell }{\pi }_{kh}^{x_{i(0)h}}\prod _{t=1}^T\prod _{h=1}^M\prod _{\ell =1}^M({\mathbf{A}}_k[h,\ell ]{)}^{x_{i(t-1)h}{x}_{i(t)\ell }},$$

$$p({\mathbf{y}}_i\mid {\mathbf{x}}_i;\mathbf{\lambda },\mathbf{\epsilon })=\prod _{t=0}^T\prod _{h=1}^{M}g(y_{i(t)};{\mathbf{\lambda }}_h,{\epsilon }_h{)}^{x_{i(t)h}},$$

$$g(y_{i(t)};{\mathbf{\lambda }}_h,{\epsilon }_h)=(1-{\epsilon }_h)g_c(y_{i(t)};{\mathbf{\lambda }}_h)+{\epsilon }_h{1}_{\{y_{i(t)}=0\}},$$

$$\forall k\in \{1,\dots ,K\},{\nu }_2({\mathbf{A}}_k)\le \nu ,$$

$$\forall h\in \{1,\dots ,M-1\},\underset{y_{i(1)}\to \rho }{\lim }\frac{g_c(y_{i(1)};{\mathbf{\lambda }}_{h+1})}{g_c(y_{i(1)};{\mathbf{\lambda }}_h)}=0.$$

$$\forall k\in \{1,\dots ,K\},\forall k^{\prime }\in \{1,\dots ,K\}\setminus \{k\},\sum _{h=1}^M\sum _{\ell =1}^M{\pi }_{kh}\log \frac{{\mathbf{A}}_k[h,\ell ]}{{\mathbf{A}}_{k^{\prime }}[h,\ell ]}>\zeta .$$

$$\forall {\mathbf{y}}_i,p({\mathbf{y}}_i;\mathbf{\theta })=p({\mathbf{y}}_i;\tilde{\mathbf{\theta }})\Rightarrow \mathbf{\theta }=\tilde{\mathbf{\theta }}.$$

$${\mathbb{P}}_0\left[\frac{\mathbb{P}(Z_{ik}=1\mid {\mathbf{y}}_i)}{\mathbb{P}(Z_{i{k}_0}=1\mid {\mathbf{y}}_i)}>a\right]\le \mathcal{O}(e^{-cT}),$$

$$\frac{\mathbb{P}(Z_{ik}=1\mid {\mathbf{y}}_i)}{\mathbb{P}(Z_{i{k}_0}=1\mid {\mathbf{y}}_i)}\underset{T\to +\infty }{\stackrel{a.s.}{⟶}}0,\mathbb{P}(Z_{i{k}_0}=1\mid {\mathbf{y}}_i)\underset{T\to +\infty }{\stackrel{a.s.}{⟶}}1\text{ and }\mathbb{P}(Z_{ik}=1\mid {\mathbf{y}}_i)\underset{T\to +\infty }{\stackrel{a.s.}{⟶}}0.$$

$$p({\mathbf{y}}_i;{\mathbf{\pi }}_k,{\mathbf{A}}_k,\mathbf{\lambda },\mathbf{\epsilon })\simeq \prod _{s=1}^{S_i}p({\mathbf{y}}_{is};{\mathbf{\pi }}_k,{\mathbf{A}}_k,\mathbf{\lambda },\mathbf{\epsilon }),$$

$${\mathbf{y}}_i=(y_{i1(0)},\dots ,y_{i1(T_{i1})},y_{i2(0)},\dots ,y_{i2(T_{i2})},\dots ,y_{i{S}_i(0)},\dots ,y_{i{S}_i(T_{i{S}_i})}),$$

$$d_{min}\ge \frac{1}{1-{\nu }_k^{\ast }}\log \frac{1}{\eta {\min }_h{\pi }_{kh}},$$

$$\underset{h}{\max }\Vert A_k^D[h,\cdot ]-{\pi }_k{\Vert }_{TV}\le \eta .$$

$${\ell }_K(\mathbf{\theta };\mathbf{y})=\sum _{i=1}^n\log \left(\sum _{k=1}^K{\delta }_k\prod _{s=1}^{S_i}p({\mathbf{y}}_{is};{\mathbf{\pi }}_k,{\mathbf{A}}_k,\mathbf{\lambda },\mathbf{\epsilon })\right).$$

$${\ell }_K(\mathbf{\theta };\mathbf{y},\mathbf{z})=\sum _{i=1}^n\sum _{k=1}^K{z}_{ik}\log {\delta }_k+\sum _{i=1}^n\sum _{k=1}^K{z}_{ik}\left(\sum _{s=1}^{S_i}\log p({\mathbf{y}}_{is};{\mathbf{\pi }}_k,{\mathbf{A}}_k,\mathbf{\lambda },\mathbf{\epsilon })\right).$$

$${\gamma }_{ikhs(t)}(\mathbf{\theta })=\mathbb{P}(X_{is(t)}=h\mid {\mathbf{y}}_{is},Z_{ik}=1;\mathbf{\theta })=\frac{{\alpha }_{ikhs(t)}(\mathbf{\theta }){\beta }_{ikhs(t)}(\mathbf{\theta })}{{\sum }_{\ell =1}^M{\alpha }_{ik\ell s(t)}(\mathbf{\theta }){\beta }_{ik\ell s(t)}(\mathbf{\theta })}.$$

$$\begin{array}[]{ll}\xi_{ikh\ell s(t)}(\boldsymbol{\theta})&\vspace{1em}=\mathbb{P}(X_{is(t)}=\ell,X_{is(t-1)}=h\mid\boldsymbol{y}_{is},Z_{ik}=1;\boldsymbol{\theta})\\ &=\displaystyle\frac{\alpha_{ikhs(t)}(\boldsymbol{\theta})\boldsymbol{A}_{k}[h,\ell]g(y_{is(t)};\boldsymbol{\lambda}_{\ell},\varepsilon_{\ell})\beta_{ik\ell s(t)}(\boldsymbol{\theta})}{\sum_{h^{\prime}=1}^{M}\sum_{\ell^{\prime}=1}^{M}\alpha_{ikh^{\prime}s(t)}(\boldsymbol{\theta})\boldsymbol{A}_{k}[h^{\prime},\ell^{\prime}]g(y_{is(t)};\boldsymbol{\lambda}_{\ell^{\prime}},\varepsilon_{\ell^{\prime}})\beta_{ik\ell^{\prime}s(t)}(\boldsymbol{\theta})}.\end{array}$$

$${\tau }_{ik}(\mathbf{\theta })=\mathbb{P}(Z_{ik}=1\mid {\mathbf{y}}_i,\mathbf{\theta })=\frac{{\prod }_{s=1}^{S_i}{\sum }_{h=1}^M{\alpha }_{ikhs(T_{is})}(\mathbf{\theta })}{{\sum }_{k^{\prime }=1}^K{\prod }_{s=1}^{S_i}{\sum }_{h=1}^M{\alpha }_{i{k}^{\prime }hs(T_{is})}(\mathbf{\theta })}.$$

$${\eta }_{ihs(t)}(\mathbf{\theta })=\mathbb{P}(X_{is(t)}=h\mid {\mathbf{y}}_i,\mathbf{\theta })=\sum _{k=1}^K{\tau }_{ik}(\mathbf{\theta }){\gamma }_{ikhs(t)}(\mathbf{\theta }).$$

$${\tau }_{ik}({\mathbf{\theta }}^{[r-1]}),{\gamma }_{ikhs(t)}({\mathbf{\theta }}^{[r-1]}),{\eta }_{ihs(t)}({\mathbf{\theta }}^{[r-1]},\text{ and }{\xi }_{ikh\ell s(t)}({\mathbf{\theta }}^{[r-1]}).$$

$${\delta }_k^{[r]}=\frac{n_k({\mathbf{\theta }}^{[r-1]})}{n},{\pi }_{kh}^{[r]}=\frac{n_{kh(0)}({\mathbf{\theta }}^{[r-1]})}{n_k({\mathbf{\theta }}^{[r-1]})},{\mathbf{A}}_k[h,\ell {]}^{[r]}=\frac{n_{kh\ell }({\mathbf{\theta }}^{[r-1]})}{n_{kh}({\mathbf{\theta }}^{[r-1]})},{\epsilon }_h^{[r]}=\frac{w_h({\mathbf{\theta }}^{[r-1]})}{n_{kh}({\mathbf{\theta }}^{[r-1]})},$$

$$\text{and }{\mathbf{\lambda }}_h^{[r]}=\underset{{\lambda }_h}{\mathrm{argmax}}\sum _{i=1}^n\sum _{s=1}^{S_i}\sum _{t=0}^{T_{is}}{\eta }_{ihs(t)}({\mathbf{\theta }}^{[r-1]})g_c(y_{is(t)};{\mathbf{\lambda }}_h),$$

$$n_k(\mathbf{\theta })=$$

$$n_{kh\ell }(\mathbf{\theta })$$

$${\epsilon }_1={\epsilon }_2=0.1,a_1=1,b_1=b_2=1,{\mathbf{A}}_1=\left[\begin{array}{cc}e& 1-e\\ 1-e& e\end{array}\right]\text{ and }{\mathbf{A}}_2=\left[\begin{array}{cc}1-e& e\\ e& 1-e\end{array}\right].$$

$${\hat{\mathbf{A}}}_1=\left[\begin{array}{cccc}0.86& 0.14& 0.00& 0.00\\ 0.12& 0.81& 0.06& 0.01\\ 0.00& 0.07& 0.79& 0.14\\ 0.00& 0.00& 0.13& 0.87\end{array}\right].$$

$$\mathrm{BIC}(K)={\ell }_K(\mathbf{\theta };\mathbf{y})-\frac{{\nu }_K}{2}\log (\sum _{i=1}^n\sum _{s=1}^{S_i}{T}_{is}+1),$$

$$\mathrm{ICL}(K)=\mathrm{BIC}(K)+\sum _{i=1}^n\sum _{k=1}^K{\hat{z}}_{ik}(\hat{\mathbf{\theta }})\log {\tau }_{ik}(\hat{\mathbf{\theta }}),$$

$$\hat{z}_{ik}(\hat{\boldsymbol{\theta}})=\left\{\begin{array}[]{rl}1&\text{if }\tau_{ik}(\hat{\boldsymbol{\theta}})=\mathrm{argmax}_{\ell=1,\ldots,K}\;\tau_{i\ell}(\hat{\boldsymbol{\theta}})\\ 0&\text{otherwise}\end{array}\right..$$

Full text

Keywords: Accelerometer data; Hidden Markov model; Longitudinal data; Missing data; Mixture models

1 Introduction

Inadequate sleep and physical inactivity affect physical and mental well-being while often exacerbating health problems. They are currently considered major risk factors for several health conditions (see, for instance Kimm et al., (2005); Taheri et al., (2004); Lee et al., (2012); Grandner et al., (2013); McTiernan, (2008)). Therefore, appropriate assessment of activity and sleep periods is essential in disciplines such as medicine and epidemiology. The use of accelerometers to evaluate physical activity—by measuring the acceleration of the part of the body to which they are attached—is a classic method that has become widespread in public health research. Indeed, since the introduction in 2003 of the first objective assessment of physical activity using accelerometers, as part of the National Health and Nutrition Examination Survey (NHANES), the analysis of actigraphy data has been the subject of extensive studies over the past two decades. Recently, the New York City (NYC) Department of Health and Mental Hygiene conducted the 2010-2011 Physical Activity and Transit (PAT) Survey111NYC Department of Health and Mental Hygiene. Physical Activity and Transit Survey 2010-2011; public use datasets accessed on May 10, 2019. The data are freely accessible on this page: https://www1.nyc.gov/site/doh/data/data-sets/physical-activity-and-transit-survey-public-use-data.page, a random survey of adult New Yorkers that tracked levels of sedentary behavior and physical activity at work, at home, and for leisure. A subset of interviewees was also invited to participate in a follow-up study to measure objectively their activity level using an accelerometer. One of the objectives of this study is to describe measured physical activity levels and to compare estimates of adherence to recommended physical activity levels, as assessed by accelerometer, with those from self-reports. In contrast to NHANES accelerometer data, PAT data still seem relatively unexplored in the statistical literature.

This paper is motivated by the analysis of the accelerometer data worn by 133 individuals aged at least of 65 years who responded PAT survey. Our objective is to propose a model adapted to the specificities of these data and study its properties. Indeed, this data set raises various challenges, such as managing the heterogeneity of the population or missing data of different natures. In order to motivate the development of a new model, we present an overview of the literature on accelerometer data analysis.

Pioneer approaches used for analyzing accelerometer data have focused on automatic detection of the sleep and wake-up periods (Cole et al.,, 1992; Sadeh et al.,, 1994; Pollak et al.,, 2001; Van Hees et al.,, 2015). More recent developments are interested in the classification of different levels of activity (see Yang and Hsu, (2010) for a review). These methods provide summary statistics like the mean time spent at different activity levels. In epidemiological studies, time spent by activity level is often used as a covariate in predictive models (see, for instance, the works of Noel et al., (2010); Palta et al., (2015); Innerd et al., (2018) where the links between physical activity and obesity are investigated). These statistics can be computed using deterministic cutoff levels (Freedson et al.,, 1998). However, with such an approach, the dependency in time is neglected and the cutoff levels are pre-specified and not estimated from the data.

Accelerometer data are characterized by a time dependency between the different measures. They can be analyzed by methods developed for functional data or by Hidden Markov Models (HMM). Methods for functional data need the observed data to be converted into a function of time (Morris et al.,, 2006; Xiao et al.,, 2014; Gruen et al.,, 2017). For instance, Morris et al., (2006) use wavelet basis for analyzing accelerometer profiles. The use of a function basis reduces the dimension of the data, and therefore the computing time. However, these methods do not define levels of activity and thus cannot directly provide the time spent at different activity levels.

When considering a discrete latent variable to model time dependence, HMM are appropriate for adjusting sequence data (Scott et al.,, 2005; Altman,, 2007; Gassiat et al.,, 2016). Titsias et al., (2016) expand the amount of information which can be obtained from HMM including a procedure to find maximum a posteriori (MAP) of the latent sequences and to compute posterior probabilities of the latent states. HMM are used on activity data for monitoring circadian rythmicity (Huang et al., 2018b, ) or directly for estimating the sequence of activity levels from accelerometer data (Witowski et al.,, 2014). For simulated data, Witowski et al., (2014) established the superiority of different HMM models, in terms of classification error, over traditional methods based on a priori fixed thresholds. While the simplicity of implementing threshold-based methods is an obvious advantage, they have some significant disadvantages compared to the HMM methods, particularly for real data. Indeed, the variation in counts and the resulting dispersion is large, leading to considerable misclassification of counts recorded in erroneous activity ranges. The approach of Witowski et al., (2014) assumes homogeneity of the population and does not consider missingness within the observations. However, heterogeneity in physical activity behaviors is often present (see, for instance, Geraci, (2018)) and the use of more than one HMM allows it to be taken into account (see, e.g., Van de Pol and Langeheine, (1990)). Clustering enables the heterogeneity of the population to be addressed by grouping observations into a few homogeneous classes. Finite mixture models (McLachlan and Peel,, 2004; McNicholas,, 2016) permit to cluster different types of data like: continuous (Banfield and Raftery,, 1993), integer (Karlis and Meligkotsidou,, 2007), categorical (Goodman,, 1974), mixed (Hunt and Jorgensen,, 2011; Kosmidis and Karlis,, 2015), network (Hoff et al.,, 2002; Hunter et al.,, 2008; Matias et al.,, 2018) and sequence data (Wong and Li,, 2000). Recent methods use clustering for accelerometer data analysis. For instance, Wallace et al., (2018) use a specific finite mixture to identify novel sleep phenotypes, Huang et al., 2018a perform a matrix-variate-based clustering on accelerometer data while Lim et al., (2019) use a clustering technique designed for functional data. Mixed Hidden Markov Models (MHMM) are a combination of HMM and Generalized Linear Mixed Models (Van de Pol and Langeheine,, 1990; Bartolucci et al.,, 2012). These models consider one (or more) random effect(s) coming from either a continuous distribution (Altman,, 2007) or a discrete distribution (Bartolucci et al.,, 2011; Maruotti,, 2011). Note that a MHMM with a single discrete random effect distribution, having a finite number of states, is a finite mixture of HMM. Such a model allows to estimate a partition among the population and to consider the population heterogeneity. The impact of the random effect can be on the measurement model or on the latent model.

This paper focuses on the analysis of PAT data with a two-fold objective: obtain summary statistics about physical activity of the subjects without pre-specifying cutoff levels and obtain a partition which groups subjects into homogeneous classes. We define a class as homogeneous if its subjects have similar average times spent into the different activity levels and similar transition probabilities between activity levels. To achieve this goal, we introduce a specific finite mixture of HMM for analyzing accelerometer data. This model considers two latent variables: a categorical variable indicating each subject’s class membership and a sequence of categorical variables indicating the subject’s level of activity each time its acceleration is measured. At time $t$, the measure is independent of the class membership conditionally on the activity level (i.e., the latent state) and follows a zero-inflated distribution—a distribution that allows for frequent zero-valued observations. The activity level defines the parameter of this distribution. The use of zero-inflated distribution is quite common for modeling accelerometer data (Ae Lee and Gill,, 2018; Bai et al.,, 2018), as the acceleration is measured every second, many observations are zero. Note that the definitions of the activity levels are equal among the mixture components. This is an important point for the use of the summary statistics (e.g., time spent at different activity levels, probabilities of transition between levels) in a future statistical study. The model we consider is thus a specific MHMM with a finite-states random effect that only impacts the distribution of latent physical activity levels. MHMM with a finite-states random effect have few developments in the literature (Bartolucci et al.,, 2011; Maruotti,, 2011), especially when the random effects only impact the latent model (and not the measurement model). We propose to theoretically study the model properties by showing that the probability of misclassifying an observation decreases at an exponential rate. In addition, since the distribution given the latent state is itself a bi-component mixture (due to the use of zero-inflated distributions), we investigate sufficient conditions for model identifiability.

In practice, the data collected often include missing intervals due to non-compliance by participants (e.g., if the accelerometer is removed). Geraci and Farcomeni, (2016) propose to identify different profiles of physical activity behaviors using a principal component analysis that allows for missing values. The PAT data contain three types of missing values corresponding to periods when the accelerometer is removed, making statistical analysis more challenging. First, missingness occurs at the beginning and at the end of the measure sequences due to the installation and the removing of the accelerometer. Second, subjects are asked to remove the accelerometer when they sleep at night. Third, missing values appear during the day (e.g., due to a shower period, nap, …). We remove missing values which occur at the begin and at the end of the sequence. For missingness caused by night sleep, we consider that the different sequences describing different days of observations of a subject are independent and that the starting point (e.g., first observed measure of the accelerometer of the day) is drawn from the stationary distribution. For missing values measured during the day, the model and the estimation algorithm can handle these data. Moreover, we propose an approximation of the distribution that avoids the computation of large powers of the transition matrices in the algorithm used for parameter inference and thus reduces computation time. Theoretical guarantees and numerical experiments show the relevance of our proposition.

The R package MHMM which implements the method introduced in this paper is available on CRAN (Du Roy de Chaumaray et al.,, 2020). It permits to analyze other accelerometer data and thus it is complementary to existing packages for MHMM. Indeed, it takes into account the specificities of accelerometer data (the class membership only impacts the transition matrices, the emission distributions are zero-inflated gamma (ZIG) distributions). Among the R packages implementing MHMM methods, one can cite the R packages LMest (Bartolucci et al.,, 2017) and seqHMM (Helske and Helske,, 2019) which focus on univariate longitudinal categorical data and the R package mHMMbayes (Aarts,, 2019) which focuses on multivariate longitudinal categorical data.

This paper is organized as follows. Section 2 presents the PAT data and the context of the study. Section 3 introduces our specific mixture of HMM and its justification in the context of accelerometer data analysis. Section 4 presents the model properties (model identifiability, exponential decay of the probabilities of misclassification and a result for dealing with the non-wear periods). Section 5 discusses the maximum likelihood inference and Section 6 illustrates the model properties on both simulated and real data. Section 7 illustrates the approach by analyzing a subset of the PAT accelerometer data. Section 8 discusses some future developments. Proofs and technical lemmas are postponed in Appendix.

2 PAT data description

We consider a subset of the data from the PAT survey, the subjects who participated in the follow-up study to measure objectively their activity level using an accelerometer. A detailed methodological description of the study and an analysis of the data is provided in Immerwahr et al., (2012). Note that the protocols for accelerometer data for the PAT survey and NHANES were identical. One of the objectives of the PAT study is to investigate the relationships between self-reported physical activity and physical activity measured by the accelerometer in order to provide best practice recommendations for the use of self-reported data Wyker et al., (2013). Indeed, self-reported data may be subject to overreporting. This is particularly the case among less active people, due in particular to a social desirability bias or the cognitive challenge associated with estimating the frequency and duration of daily physical activity (see, e.g., Slootmaker et al., (2009); Dyrstad et al., (2014); Lim et al., (2015)). The results of Wyker et al., (2013) show that males tend to underreport their physical activity, while females and older adults (65 years and older) overreported (see also Troiano et al., (2008) for a detailed study of the differences between self-reported physical activity and accelerometer measurements in NHANES 2003-2004). Consequently, the study of data measured by accelerometer for these specific populations makes it possible to determine methods for correcting estimates from self-reported data, such as stratification by gender and/or age when comparing groups.

In this work, we are particularly interested in the age category above 65 years old ($n=133$). We present some characteristics related to PAT data and refer to Immerwahr et al., (2012) for a full description222Raw accelerometer data, covariates allowing the selection of the subset of the population, as well as a detailed dictionary are freely accessible here: https://www1.nyc.gov/site/doh/data/data-sets/physical-activity-and-transit-survey-public-use-data.page. Accelerometers were worn for one week (beginning on Thursday and ending on Wednesday) and measured the activity minute-by-minute. The trajectory associated with each subject is therefore of length 10080. In addition, a participant’s day spans from 3am-3am (and not a calendar day) in order to record late night activities and transit and contains missing data sequences of variable length at the beginning and end of the measurement period (these missing data sequences were excluded from the analysis). This length varies from one subject to another, and the mean and minimum trajectory length for the population under consideration (after excluding those missing at the edges) are 9474 and 5199 respectively (with a total number of observations equal to 1259981). The model of the accelerometer used is Actigraph GT3X, it is worn on the hips (which results in the fact that certain activities, such as lifting weights or biking, cannot be measured). In addition, participants were also asked to remove it when sleeping, swimming or bathing, hence the data contains approximately 44$\%$ of missing values that appear mainly in sequence, appearing at night but also during the day. Figure 1 gives an example of accelerometer data measured on one subject (i.e., patcid:1200255) for one week where the three types of missing data can be seen. The four levels of physical activity based on the classification established by the US Department of Health and Human Services, (2008) in the Physical Activity Guidelines for Americans (PAGA) report is also shown in the Figure 1. Specifically, the PAT protocol for accelerometer data has established a classification according to PAGA, characterizing each minute of activity. Activity minutes with less than 100 activity counts were classified as Sedentary, minutes with 100-2019 counts were classified as Light, the class Moderate corresponds to 2020-5998 counts/minute and Vigorous 5999 and above counts/minute. A comparison between our method and this traditional threshold-based approach is provided in Section 7.3.

3 Mixture of hidden Markov models for accelerometer data

In this section we present the proposed model and the application context for which it has been defined.

3.1 The data

Observed data $\boldsymbol{y}=(\boldsymbol{y}_{1}^{\top},\ldots,\boldsymbol{y}_{n}^{\top})$ are composed of $n$ independent and identically distributed sequences $\boldsymbol{y}_{i}$. Each sequence $\boldsymbol{y}_{i}=(y_{i(0)},\ldots,y_{i(T)})$ which contains the values measured by the accelerometer at times $t\in\{0,1,\ldots,T\}$ for subject $i$, with $y_{i(t)}\in\mathbb{R}^{+}$. Throughout the paper, index $i$ is related to the label of the subject and index $(t)$ is related to the time of measurement.

The model considers $M$ different activity levels (which are unobserved). These levels impact the distribution of the observed sequences of accelerometer data. The sequence of the hidden states $\boldsymbol{x}_{i}$ indicates the activity level of subject $i$ at the different times. Thus, $\boldsymbol{x}_{i}=(\boldsymbol{x}_{i(0)},\ldots,\boldsymbol{x}_{i(T)})\in\mathcal{X}$ and the activity level (among the $M$ possible levels) of subject $i$ at time $t$ is defined by the binary vector $\boldsymbol{x}_{i(t)}=(x_{i(t)1},\ldots,x_{i(t)M})$ where $x_{i(t)h}=1$ if subject $i$ is at state $h$ at time $t$ and $x_{i(t)h}=0$ otherwise.

The heterogeneity (in the sense of different physical activity behaviors) between the $n$ subjects can be addressed by grouping subjects into $K$ homogeneous classes. This is achieved by clustering that assesses a partition $\boldsymbol{z}=(\boldsymbol{z}_{1},\ldots,\boldsymbol{z}_{n})$ among the $n$ subjects based on their accelerometer measurements. Thus the vector $\boldsymbol{z}_{i}=(z_{i1},\ldots,z_{iK})$ indicates the class membership of subject $i$, as $z_{ik}=1$ if observation $i$ belongs to class $k$ and $z_{ik}=0$ otherwise. Throughout the paper, index $k$ refers to the label of a class grouping homogeneous subjects.

Each subject $i$ is described by three random variables: one unobserved categorical variable $\boldsymbol{z}_{i}$ (defining the membership of the class of homogeneous physical activity behaviors for subject $i$), one unobserved categorical longitudinal data $\boldsymbol{x}_{i}$ (a univariate categorical discrete-time time series which defines the activity level of subject $i$ at each time) and one observed positive longitudinal data $\boldsymbol{y}_{i}$ (a univariate positive discrete-time time series which contains the values of the accelerometer measured on subject $i$ at each time).

3.2 Generative model

The model described below considers that the observations are independent between the subjects and identically distributed. It is defined by the following generative model and summarized by Figure 2 (note that this figure is similar to Figure 6.2 of Bartolucci et al., (2012)):

sample class membership $\boldsymbol{z}_{i}$ from a multinomial distribution; 2. 2.

sample the sequence of activity levels $\boldsymbol{x}_{i}$ from a Markov model whose transition matrix depends on class membership; 3. 3.

sample the accelerometer measurement sequence given the activity levels (each $\boldsymbol{Y}_{i(t)}$ follows a ZIG distribution whose parameters are defined only by $\boldsymbol{x}_{i(t)}$).

3.3 Finite mixture model for heterogeneity

The sequence of accelerometer measures obtained on each subject is assumed to independently arise from a mixture of $K$ parametric distributions, so that the probability distribution function (pdf) of the sequence $\boldsymbol{y}_{i}$ is

[TABLE]

where $\boldsymbol{\theta}=\{\boldsymbol{\lambda},\boldsymbol{\varepsilon}\}\cup\{\delta_{k},\boldsymbol{\pi}_{k},\boldsymbol{A}_{k};k=1,\ldots,K\}$ groups the model parameters, $\delta_{k}$ is the proportion of components $k$ with $\delta_{k}>0$, $\sum_{k=1}^{K}\delta_{k}=1$, and $p(\cdot;\boldsymbol{\pi}_{k},\boldsymbol{A}_{k},\boldsymbol{\lambda},\boldsymbol{\varepsilon})$ is the pdf of component $k$ parametrized by $(\boldsymbol{\pi}_{k},\boldsymbol{A}_{k},\boldsymbol{\lambda},\boldsymbol{\varepsilon})$ defined below. Thus, $\delta_{k}$ is the marginal probability that a subject belongs to class $k$ (i.e., $\delta_{k}=\mathbb{P}(Z_{ik}=1)$). Moreover, $p(\cdot;\boldsymbol{\pi}_{k},\boldsymbol{A}_{k},\boldsymbol{\lambda},\boldsymbol{\varepsilon})$ defines the distribution of a sequence of values measured by the accelerometer on a subject belonging to class $k$ (i.e., $p(\cdot;\boldsymbol{\pi}_{k},\boldsymbol{A}_{k},\boldsymbol{\lambda},\boldsymbol{\varepsilon})$ is the pdf of $\boldsymbol{y}_{i}$ given $Z_{ik}=1$).

3.4 Hidden Markov model for activity levels

The model assumes that the distribution of the hidden state sequence depends on the class membership, and that the distribution of activity measurements depends on the state at time $t$ but not on the component membership given the state (i.e., $\boldsymbol{X}_{i}\not\perp\boldsymbol{Z}_{i}$, $Y_{i(t)}\not\perp\boldsymbol{X}_{i(t)}$ and $Y_{i(t)}\perp\boldsymbol{Z}_{i}\mid\boldsymbol{X}_{i(t)}$). It is crucial that the distribution of $Y_{i(t)}$ given $\boldsymbol{X}_{i(t)}$ is independent to $\boldsymbol{Z}_{i}$. Indeed, each activity level is defined by the distribution of $Y_{i(t)}$ given the state. Therefore, to extract summary statistics on the whole population (as the average time spent by level of activity) the definition of the activity levels (and the distribution of $y_{i(t)}$ given the state) must to be the same among the mixture components.

The pdf of $\boldsymbol{y}_{i}$ for components $k$ (i.e., given $Z_{ik}=1$) is

[TABLE]

The Markov assumption implies that

[TABLE]

where $\boldsymbol{\pi}_{k}=(\pi_{k1},\ldots,\pi_{kM})$ defines the initial probabilities so that $\pi_{kh}=\mathbb{P}(X_{i(1)h}=1\mid Z_{ik}=1)$, $\boldsymbol{A}_{k}$ is the transition matrix so that $\boldsymbol{A}_{k}[h,\ell]=\mathbb{P}(X_{i(t)\ell}=1\mid X_{i(t-1)h}=1,Z_{ik}=1)$. Finally, we have

[TABLE]

where $g(\cdot;\boldsymbol{\lambda}_{h},\varepsilon_{h})$ is the pdf of a zero-inflated distribution defined by

[TABLE]

where $g_{c}(\cdot;\boldsymbol{\lambda}_{h})$ is the density of a distribution defined on a positive space and parametrized by $\boldsymbol{\lambda}_{h}$. The choice of considering zero-inflated distributions is motivated by the large number of zero in the accelerometer data (see Figure 1). For the application of Section 7, we use a gamma distribution of $g_{c}(\cdot;\boldsymbol{\lambda}_{h})$. However, model properties and inference are discussed for a large family of densities $g_{c}(\cdot;\boldsymbol{\lambda}_{h})$.

4 Model properties

In this section, we present the properties of the mixture of parametric HMM. We start with a discussion of three assumptions. Then, model identifiability is proved. It is shown that the probability of making an error in the partition estimation exponentially decreases with $T$, when the model parameters are known. Finally, the analysis of missing data is discussed.

4.1 Assumptions

Assumption 1**.**

For each component $k$, the Markov chain is irreducible. Moreover, we assume that the sequence is observed at its stationary distribution (i.e., $\boldsymbol{\pi}_{k}$ is the stationary distribution so $\boldsymbol{\pi}_{k}^{\top}\boldsymbol{A}_{k}=\boldsymbol{\pi}_{k}^{\top}$). Therefore, there exists $0\leq\nu<1$ such that

[TABLE]

where $\nu_{2}(\boldsymbol{A}_{k})$ is the second-largest eigenvalue of $\boldsymbol{A}_{k}$. Finally, we denote by $\bar{\nu}_{2}(\boldsymbol{A}_{k})=\max(0,\nu_{2}(\boldsymbol{A}_{k}))$.

Assumption 2**.**

The hidden states define different distributions for the observed sequence. Therefore, for $h\in\{1,\ldots,M\}$, $h^{\prime}\in\{1,\ldots,M\}\setminus\{h\}$, we have $\boldsymbol{\lambda}_{h}\neq\boldsymbol{\lambda}_{h^{\prime}}$. Moreover, the parametric family of distributions defining $g_{c}(\cdot;\boldsymbol{\lambda}_{1}),\ldots,g_{c}(\cdot;\boldsymbol{\lambda}_{M})$ permits to consider an ordering such that for a fix value $\rho\in\mathbb{R}^{+}\setminus\{0\}$, we have

[TABLE]

Assumption 3**.**

The transition probabilities are different over the mixture components and are not zero. Therefore, for $k\in\{1,\ldots,K\}$, $k^{\prime}\in\{1,\ldots,K\}\setminus\{k\}$, we have $\forall(h,\ell)$, $\boldsymbol{A}_{k}[h,\ell]\neq\boldsymbol{A}_{k^{\prime}}[h,\ell]$. Moreover, there exists $\zeta>0$ such that

[TABLE]

Finally, without loss of generality, we assume that $A_{k}[1,1]>A_{k+1}[1,1]$.

Assumption 1 states that the state at time 1 is drawn from the stationary distribution of the component that the observation belongs to. To obtain the model identifiability we do not need the assumption that the stationary distribution is different over the mixture components. As a result, two components having the same stationary distribution but different transition matrices can be considered. Assumption 2 and Assumption 3 are required to obtain the model identifiability. Assumption 3 can be interpreted as the Kullback-Leibler divergence between the distribution of the states under component $k$ and their distribution under component $k^{\prime}$. This constraint is required for model identifiability because it is related to the definition of the classes. Consequently, the matrices of the transition probability must be different among components.

4.2 Identifiability

Model identifiability is crucial for interpreting the estimators of the latent variables and of the parameters. It has been studied for some mixture models (Teicher,, 1963, 1967; Allman et al.,, 2009; Celisse et al.,, 2012) and HMM (Gassiat et al.,, 2016), but not for the mixture of HMM. Generic identifiability (up to switching of the components and of the states) of the model defined in (1) implies that

[TABLE]

The following theorem states this property.

Theorem 1**.**

If Assumptions 1, 2 and 3 hold, then the model defined in (1) is generically identifiable (up to switching of the components and of the states) if $T>2K$.

Proof of Theorem 1 is given in Appendix, Section A. The model defined by the marginal distribution of a single $y_{i(t)}$ is not identifiable. Indeed, the marginal distribution of $y_{i(t)}$ is a mixture of zero-inflated distributions and such mixture is not identifiable (i.e., different class proportions and inflation proportions can define the same distribution). It is therefore this dependency over time that makes the proposed mixture generically identifiable. Note that such statement has been made by Gassiat et al., (2016) when they discuss the case where the emission distribution for an HMM follows a mixture model.

4.3 Probabilities of misclassification

In this section, we examine the probability that an observation will be misclassified when the model parameters are known. We consider the ratio between the probability that subject $i$ belongs to class $k$ given $\boldsymbol{y}_{i}$ and the probability that this subject belongs to its true class, and we quantify the probability of it being greater than some positive constant $a$. Let $\theta_{0}$ be the true model parameter and $\mathbb{P}_{0}=\mathbb{P}(\cdot\mid Z_{ik_{0}}=1,\theta_{0})$ denote the true conditional distribution (true label of subject $i$ and parameters are known).

Theorem 2**.**

Assume that Assumptions 1 and 3 hold. If $a>0$ is such that Assumption 4 (defined in Appendix Section B) holds, then for every $k\neq k_{0}$

[TABLE]

where $c>0$ is a positive constant

Moreover, the exponential bounds of Theorem 2 allows to use the Borel-Cantelli’s lemma to obtain the almost sure convergence.

Corollary 1**.**

Assume that Assumptions 1 and 3 hold. If $\boldsymbol{y}_{i}$ is generated form component $k_{0}$ (i.e., $Z_{ik_{0}}=1$), then for every $k\neq k_{0}$

[TABLE]

Therefore, by considering $a=1$, Theorem 2 and Corollary 1 show that the probability of misclassifying the subject $i$ based on the observation $\boldsymbol{y}_{i}$, using the maximum a posteriori rule, tends to zero when $T$ increases, if the model parameters are known. Proof of Theorem 2 and a sufficient condition that allows to consider $a=1$ (value of interest when the partition is given by the MAP rule) are given in Appendix, Section B. It should be noted that it is not so common to have an exponential rate of convergence for the ratio of the posterior probability of classification. Similar results are obtained for network clustering using the stochastic block model (Celisse et al.,, 2012) or for co-clustering (Brault and Mariadassou,, 2015). For these two models, the marginal distribution of a single variable provides information about the class membership. For the proposed model, this is the dependency between the different observed variables which is the crucial point for recovering the true class membership.

4.4 Dealing with missing values

Due to the markovian character of the states, missing values can be handled by iterating the transition matrices. In our particular context, missing values appear when the accelerometer is not worn (see Section 2 for explanations of the reasons of missingness). We will not observe isolated missing values but rather wide ranges of missing values. Let $d$ be the number of successive missing values, we thus have to compute the matrix $A_{k}^{d+1}$ to obtain the distribution of the state at time $t+d$ knowing the state at time $t-1$. These powers of transition matrices should be computed many times during the algorithm used for inference (see Section 5). Moreover, after $d+1$ iterations with $d$ large enough, the transition matrix can be considered sufficiently close to stationarity (e.g., for any $(h,\ell)$, $A_{k}^{d+1}[h,\ell]\simeq\pi_{k\ell}$), which has actually been chosen as the initial distribution. Therefore, for numerical reasons, we will avoid computing the powers of the transition matrices and we will make the following approximation. An observation $\boldsymbol{y}_{i}$ with $S_{i}$ observed sequences split with missing value sequences of size at least $d$ are modeled as $S_{i}$ independent observed sequences with no missing values, all belonging to the same component $k$. Namely, for each individual $i$, the pdf $p(\boldsymbol{y}_{i};\boldsymbol{\pi}_{k},\boldsymbol{A}_{k},\boldsymbol{\lambda},\boldsymbol{\varepsilon})$ of component $k$ is approximated by the product of the pdf of the $S_{i}$ observed sequences $\boldsymbol{y}_{i1},\boldsymbol{y}_{i2},\ldots,\boldsymbol{y}_{iS_{i}}$:

[TABLE]

where, for each $s$, $\boldsymbol{y}_{is}$ is an observed sequence of length $T_{is}+1$: $\boldsymbol{y}_{is}=(y_{is(0)},\ldots,y_{is(T_{is})})$ and $p(\boldsymbol{y}_{is};\boldsymbol{\pi}_{k},\boldsymbol{A}_{k},\boldsymbol{\lambda},\boldsymbol{\varepsilon})$ is defined as in (2). We note that the observation $\boldsymbol{y}_{i}$ can thus be rewritten as follows

[TABLE]

with $y_{i2(0)}=y_{i(T_{i1}+d_{i1}+1)}$ where the $d_{i1}$ values $y_{i(T_{i1}+1)},\ldots y_{i(T_{i1}+d_{i1})}$ correspond to the first sequence of missing values, and more generally, for each $s=2,\ldots,S_{i}$, $y_{is(0)}=y_{i(\sum_{j=1}^{s-1}(T_{ij}+d_{ij}+1))}$, with $d_{ij}$ being the number of missing values between the observed sequences $\boldsymbol{y}_{is_{j}}$ and $\boldsymbol{y}_{is_{j+1}}$.

Once the estimation of the parameters has been done, we make sure that this assumption was justified by verifying that the width of the smallest range $d_{min}=\min\left\{d_{i1},\ldots,d_{i\,S_{i}-1}\right\}$ of missing values is sufficiently large to be greater than the mixing time of the obtained transition matrix. To do so, we use an upper bound for the mixing time given by Levin and Peres, (2017, Theorem 12.4, p. 155). For each component $k$, we denote by $\nu_{k}^{*}$ the second maximal absolute eigenvalue of $\boldsymbol{A}_{k}$. For any positive $\eta$, if for each $k$

[TABLE]

then for any integer $D\geq d_{min}$, the maximum distance in total variation satisfies

[TABLE]

5 Maximum likelihood inference

This section presents the methodology used to estimate the model parameters.

5.1 Inference

We proposed to estimate the model parameters by maximizing the log-likelihood function where missing values are managed as in Section 4.4 and we recall that the log-likelihood is also approximated for numerical reasons, to avoid computing large powers of the transition matrices. We want to find $\hat{\boldsymbol{\theta}}$ which maximizes the following approximated log-likelihood function

[TABLE]

This maximization is achieved via an EM algorithm (Dempster, A.P. and Laird, N.M. and Rubin, D.B.,, 1977) which considers the complete-data log-likelihood defined by

[TABLE]

5.2 Conditional probabilities

Let $\alpha_{ikhs(t)}(\boldsymbol{\theta})$ be the probability of the partial sequence $y_{is(0)},\ldots,y_{is(t)}$ and ending up in state $h$ at time $t$ under component $k$. Moreover, let $\beta_{ikhs(t)}(\boldsymbol{\theta})$ be the probability of the ending partial sequence $y_{is(t+1)},\ldots,y_{is(T_{is})}$ given a start in state $h$ at time $t$ under component $k$. These probabilities can be easily obtained by the forward/backward algorithm (see Appendix, Section C). We deduce that the probability $\gamma_{ikhs(t)}(\boldsymbol{\theta})$ of being in state $h$ at time $t\in\{0,\ldots,T_{is}\}$ for $\boldsymbol{y}_{i}$ under component $k$ is

[TABLE]

The probability $\xi_{ikh\ell s(t)}(\boldsymbol{\theta})$ of being in state $\ell$ at time $t\in\Omega_{i}$ and in state $h$ at time $t-1$ for observation $\boldsymbol{y}_{i}$ under component $k$ is

[TABLE]

The probability $\tau_{ik}$ that one observation arises from component $k$ is

[TABLE]

The probability $\eta_{ihs(t)}$ that observation $i$ is at state $h$ at time $t$ of sequence $s$ is

[TABLE]

5.3 EM algorithm

The EM algorithm is an iterative algorithm randomly initialized at the model parameter $\boldsymbol{\theta}^{[0]}$. It alternates between two steps: the Expectation step (E-step) consisting in computing the expectation of the complete-data likelihood under the current parameters, and the maximization step (M-step) consisting in maximizing this expectation over the model parameters. Iteration $[r]$ of the algorithm is defined by

E-step Conditional probability computation, updating of

[TABLE]

M-step Parameter updating

[TABLE]

[TABLE]

where

[TABLE]

6 Numerical illustrations

This section aims to highlight the main properties of the model on numerical experiments. First, simulated data are used to illustrate the exponential decay of the probabilities of misclassification (given by Theorem 2), the convergence of estimators and the robustness of the approach to missingness. Second, our approach is applied to the data from the PAT study. All the experiments are conducted with the R package MHMM available on CRAN.

6.1 Simulated data

Simulation design

All the simulations are performed according to the same model. This model is a bi-components mixture of HMM with two states (i.e., $K=M=2$) and equal proportions (i.e., $\delta_{1}=\delta_{2}=1/2$). The distribution of $Y_{i(t)}$ conditionally on the state $h$ is a ZIG distribution. We have

[TABLE]

The parameter $a_{2}>1$ controls the separation of the distribution of $Y_{i(t)}$ given the state. The parameter $e$ controls the separation of the distribution of $X$ given the class (when $e$ increases, the constant $c$ in Theorem 2 increases). We consider four cases: hard ($e=0.75$ and $a_{2}=3$), medium-hard ($e=0.90$ and $a_{2}=3$), medium-easy ($e=0.75$ and $a_{2}=5$) and easy ($e=0.90$ and $a_{2}=5$).

Illustrating the exponential rate of the probabilities of misclassification

Theorem 2 states that the probabilities of misclassification decrease at an exponential rate with $T$. To illustrate this property, 1000 sequences are generated for $T=1,\ldots,100$ and the four cases. For each sequence $\boldsymbol{y}_{i}$, we compute $\log({\mathbb{P}(Z_{ik}=1\mid\boldsymbol{y}_{i})}/{\mathbb{P}(Z_{ik_{0}}=1\mid\boldsymbol{y}_{i})})$ when $k_{0}$ is the true class, $k$ the alternative and the true model parameters are used. Figure 3(a) shows the behavior of $\log({\mathbb{P}(Z_{ik}=1\mid\boldsymbol{y}_{i})}/{\mathbb{P}(Z_{ik_{0}}=1\mid\boldsymbol{y}_{i})})$ (the median of this log ratio is plotted in plain and a 90$\%$ confidence interval is plotted in gray). Note that this log ratio of probabilities linearly decreases with $T$ which illustrates the exponential decay of the probabilities of misclassification. Moreover, Figure 3(b) presents the empirical probabilities of misclassification and thus also illustrates Theorem 2. As expected, this shows that the decay of the probabilities of misclassification is faster as the overlaps between class decreases.

Illustrating the convergence of the estimators

We illustrate the convergence of the estimators (partition, latent states and parameters) when the model parameters are estimated by maximum likelihood (see Section 5). We compute the mean square error (MSE) between the model parameters and their estimators. Moreover, we compute the adjusted Rand index (ARI; Hubert and Arabie, (1985)) between the true partition and the partition given by the MAP rule, and between the true state sequences and the estimated state sequences given by the MAP rule (obtained with the Viterbi algorithm (Viterbi,, 1967)). Table 1 shows the results obtained with two different sample sizes $n$ and two different lengths of sequences $T$, considering the case medium-hard. It can be seen that the partition and the model parameters are well estimated. Indeed, the MLE converge to the true parameters as $T$ or $n$ increases, except for the proportion of each component $\delta_{k}$. The convergence of the estimator of the proportions depends mainly on the sample size $n$. We notice that the partition obtained by our estimation procedure corresponds to the true partition (for $n$ and $T$ large enough) even if we are not under the true parameters but under the MLE, which is not an immediate consequence of Theorem 2. On the contrary, we do not find the true state sequences a.s., as the number of states to be estimated is also growing with $n$ and $T$. This result was expected because the number of latent states increases with $T$ and $n$ while the number of parameters and the dimension of the partition does not increase with $T$. Results obtained for the three other cases are similar and are presented in Appendix, Section D.1.

Illustrating the robustness to missingness

We now investigate the robustness of the proposed method with missingness. We compare the accuracy of the estimators (ARI for the latent variables and MSE for the parameters) obtained on samples without missingness to the accuracy of the estimators obtained when missingness is added to the samples. Three situations of missingness are considered: missing completely at random-1 (MCAR-1) (i.e., one sequence of 10 missing values is added to each sequence $\boldsymbol{y}_{i}$, the location of the sequence follows a uniform distribution), MCAR-2 (i.e., two sequences of 20 missing values are added for each sequence $\boldsymbol{y}_{i}$, the location of the sequences follows a uniform distribution) and missing not at random (MNAR) (i.e., the probability to observe the value $y_{i(t)}$ is equal to ${e^{y_{i(t)}}}/{(1+e^{y_{i(t)}})}$). Note that the last situation adds many missing values when the true value of $y_{i(t)}$ is close to zero, so the occurrence of missing values depends on the latent states. Table 2 compares the results obtained with and without missingness, considering case medium-hard. It shows that estimators are robust to missingness. Results obtained for the other three cases are similar and are reported in Appendix, Section D.1.

6.2 Using the approach on classical accelerometer data

We consider the accelerometer data measured on three subjects available from Huang et al., 2018b . The accelerometer measures the activity every five minutes for one week. Note that the first subject has $2\%$ of missing values. The purpose of this section is to illustrate the differences between the method of Huang et al., 2018b and the method proposed in this paper.

Huang et al., 2018b consider one HMM per subject with three latent states. This model is used for monitoring the circadian rhythmicity, subject by subject. Because they fit one HMM per sequence measured by the accelerometer of a subject, the definition of the activity level is different for each subject (see, Huang et al., 2018b (, Figure 4)). This is not an issue for their study because the analysis is done subject by subject. However, the mean time spent by activity levels cannot be compared among the subjects.

The method proposed here makes this comparison possible. Figure 4 depicts the activity data of the three subjects, the expected value of $Y_{i(t)}$ conditionally to the most likely state and on the most likely component and the probability of each state. Based on the QQ-plot (see, Appendix, Section D.2), we consider $M=4$ activity levels. These levels can be easily characterized with the model parameters presented in Table 3. Moreover, the transition matrices also make sense. For instance, class 1 (subjects 9 and 20) has an almost tri-diagonal transition matrix (by considering an order between the states given through the activity levels per state) and class-2 (subject 2) is composed of a subject with low-overall activity

[TABLE]

7 Analysis of PAT data

In this section, we analyze the data presented in Section 2.

7.1 Experimental conditions

In order to compare our approach to the cuts defined a priori in the PAT study (see Section 2), the model was fitted with four activity levels. Note that selecting the number of states in HMM stays a challenging problem (see the discussion in the conclusion). However, approaches considering four activity levels are standard for accelerometer data. The number of components (i.e., the number of classes) is estimated, using an information criterion unlike the PAT study where it is arbitrarily set at 3 or 4. For each number of components, 5000 random initializations of the EM algorithm are performed. The analysis needs about one day of computation on a 32-Intel(R) Xeon(R) CPU E5-4627 v4 @ 2.60GHz.

7.2 Model selection

To select the number of components, we use two information criteria which are generally used in clustering: the BIC (Schwarz,, 1978) and the ICL (Biernacki et al.,, 2000) defined by

[TABLE]

and

[TABLE]

where $\nu_{K}=(K-1)+K(M+M^{2})+3M$ is the number of parameters for a model with $K$ components and $M$ states and $\hat{z}_{ik}(\hat{\boldsymbol{\theta}})$ defines the partition by the MAP rule associated to the MLE such that

[TABLE]

The ICL is defined according as the integrated complete-data likelihood computed with the partition given by the MAP rule with the MLE. The values of the information criteria are given in Table 4, for different number of classes. Both criteria select five components. The values of $\mathrm{ICL}(K)$ are close to those of the $\mathrm{BIC}(K)$, implying that the entropy $\sum_{i=1}^{n}\sum_{k=1}^{K}\hat{z}_{ik}\log\tau_{ik}(\hat{\boldsymbol{\theta}})\approx 0$. This is a consequence of Theorem 2 (see also numerical experiments in Section 6). In the following, we interpret the results obtained with $M=4$ activity levels and $K=5$ classes.

7.3 Description of the activity levels

The parameters of the ZIG distributions are presented in Table 5. The four distributions are ordered by the value of their means. The sleeping state is characterized by a large probability of observing zero (i.e., $\varepsilon_{h}$ is close to one). However, $\varepsilon_{h}$ is not equal to zero for the other states, but the more active the state is, the smaller $\varepsilon_{h}$ is. We also compute the marginal cutoffs (i.e., the cutoffs by considering the MAP of $\mathbb{P}(X_{i(t)}\mid Y_{i(t)})$). These cutoffs neglect the time dependency due to the Markov structure, but can be compared to the cutoffs proposed by the PAT study. Indeed, according to the PAT study, minutes with $<100$ counts are assigned to Sedentary activity, minutes with 100-2019 counts were classified as Light, the class Moderate corresponds to 2020-5998 counts/minute and Vigorous 5999 and above counts/minute. The marginal cutoff associated with the low-level state is very close to that of the Sedentary class of the PAT. We find, however, that our marginal cutoffs are more accurate for higher levels of activity. PAT cutoffs do not adequately characterize the activity level of the study population. Finally, contrary to classical thresholds, our modeling approach allows to capture and characterize the variability associated with the different levels of activity, variability which seems important (see Figure 5 and Table 5).

7.4 Description of the classes

Classes can be described using their proportions and their associated parameters presented in Table 5. The data are composed of a majority class ($\delta_{1}=0.518$). Three other classes are composed of more sedentary individuals (e.g., their marginal probabilities of being in states 1 and 2 are higher). Finally, there is a small class ($\delta_{5}=0.045$) which contains the most active subjects (i.e., $\pi_{k4}=0.143$). For three of the five classes, Figure 5 presents a characteristic subject of each class and the probabilities of the activity levels (the associated graphs for the two remaining classes are given in Appendix, Section D.3). Classes can be interpreted from the mean time spent at different activity levels presented in Table 6 and from transition matrices presented in Table 7 which are almost tri-diagonal. This could be expected because it seems relevant to obtain a low probability of jumping between the sleeping state and the intensive state. Additionally, the approximation made for efficiently handling the missingness (see Section 4.4) turns out to be relevant. The minimal range of missing values is indeed equal to $d_{\min}=60$ which leads to a distance in total variation between the $d_{\min}$-power of the transition matrices and the stationary distribution being less than $5.10^{-4}$ for any component.

8 Conclusion

A specific mixture of HMM has been introduced to analyze accelerometer data. It avoids the traditional cutoff point method and provides a better characterization of activity levels for the analysis of these data, while adapting to the population. The proposed model could be applied to a population with different characteristics (e.g., younger) which would lead to different definitions of activity levels. In addition, the use of several HMMs make to take into account dependency over time and thus improve the traditional method based on cutoff points (Witowski et al.,, 2014). This approach also allows us to take into account the heterogeneity of the population (in the sense of physical activity).

An interesting perspective is to consider adjusting for confusing factors (e.g., gender or age). These confusing factors could impact the probabilities of transition between the latent spaces (e.g., using a generalized linear model approach) and/or the definition of the accelerometer measurement given a state (e.g., linear regression on some parameters of the ZIG distribution).

In the application, the number of activity levels was not estimated but fixed at a common value for accelerometer data. Estimating the number of states for a mixture of HMM is an interesting but complex topic: for instance, the use of BIC is criticized (see, Cappé et al., (2005, Chapter 15)). This makes the study of relevant information criteria for selecting the number of states an interesting topic for future work. Pseudo-likelihood based criteria could be used (Gassiat,, 2002; Csiszár and Talata,, 2006) but the fact that the marginal distribution of one $Y_{i(t)}$ is not identifiable limits this approach. A more promising approach could be to use cross-validated likelihood (Celeux and Durand,, 2008) but it would be computationally intensive because accelerometer data provides a large number of observations.

Appendix A Model identifiability

The proof of Theorem 1 is split in two parts:

Identifiability of the parameters of the specific distribution per state is obtained using the approach of Teicher, (1963). Hence $\forall h=1,\ldots,M$

[TABLE] 2. 2.

Identifiability of the transition matrices and of the $\varepsilon$ is shown using properties of Vandermonde matrices. Hence,

[TABLE]

A.1 Identifiability of the parameters of the specific distribution per state

Considering the marginal distribution at time $t=0$, we have

[TABLE]

Note that $g(y_{i(0)};\boldsymbol{\lambda}_{h},\varepsilon_{h})=(1-\varepsilon_{h})g_{c}(y_{i(0)};\boldsymbol{\lambda}_{h})+\varepsilon_{h}\mathbf{1}_{\{y_{i(0)}=0\}}$ is a pdf of a zero-inflated distribution, so it is a pdf of a bi-component mixture. We now use the same reasoning as Teicher, (1963). We have

[TABLE]

Considering $y_{i(0)}\to\rho$, by Assumption 2, we have

[TABLE]

Repeating the previous argument with $h=2,\ldots,M$, we conclude that, for $h\in\{1,\ldots,M\}$,

[TABLE]

A.2 Identifiability of the transition matrices

First, we introduce two technical lemmas of which proofs are discussed in the next subsection. Second, we show that $\boldsymbol{A}_{k}[1,1]=\tilde{\boldsymbol{A}}_{k}[1,1]$ then we extend the results to the whole transition matrices.

Lemma 1**.**

Let $N_{0}$, $N_{1}$, $\tilde{N}_{0}$ and $\tilde{N}_{1}$ be four definite positive matrices of size $K\times K$ such that for $u\in\{1,\ldots,K\}$ and $k\in\{1,\ldots,K\}$,

[TABLE]

with $a_{k}>a_{k+1}>0$, $\tilde{a}_{k}>\tilde{a}_{k+1}>0$ and $a_{1}\geq\tilde{a}_{1}$. If for any $\tilde{w}\in\mathbb{R}^{K}_{+}$ there exists $w\in\mathbb{R}^{K}_{+}$ $N_{0}w=\tilde{N}_{0}\tilde{w}$ and $N_{1}w=\tilde{N}_{1}\tilde{w}$ then for $k\in\{1,\ldots,K\}$ $a_{k}=\tilde{a}_{k}$ and $w=\tilde{w}$.

Lemma 2**.**

Let $N_{0}$, $\tilde{N}_{0}$ be two definite positive matrices of size $K\times K$ such that for $u\in\{1,\ldots,K\}$ and $k\in\{1,\ldots,K\}$,

[TABLE]

with $a_{k}>a_{k+1}>0$, $\tilde{a}_{k}>\tilde{a}_{k+1}>0$ and $a_{1}\geq\tilde{a}_{1}$. Let $D_{u}=\mathrm{diag}(a_{1}^{Ku},\ldots,a_{K}^{Ku})$ and $\tilde{D}_{u}=\mathrm{diag}(\tilde{a}_{1}^{Ku},\ldots,\tilde{a}_{K}^{Ku})$. If there exist $\alpha\in]0,1[$, $\tilde{\alpha}\in]0,1[$, $w\in\mathbb{R}^{K}_{+}$ and $\tilde{w}\in\mathbb{R}^{K}_{+}$ such that for $u\in\{0,\ldots,K-1\}$, we have

[TABLE]

then for $k\in\{1,\ldots,K\}$ $a_{k}=\tilde{a}_{k}$ and $w=\tilde{w}$.

We consider the marginal distribution of $(y_{i(0)},\ldots,y_{i(t-1)})$ with $t=1,\ldots,2K$, where $y_{i(0)}=y_{i(t^{\prime})}$ for each $t^{\prime}=1,\ldots,t-3$, $y_{i(t-2)}=y_{i(0)}^{\tau_{1}}$ , $y_{i(t-1)}=y_{i(0)}^{\tau_{2}}$ and $y_{i(t)}=y_{i(0)}^{\tau_{3}}$. Therefore, taking $\tau_{1}=\tau_{2}=\tau_{3}=1$ and letting $y_{i(0)}$ tend to $\rho$ (see Assumption A), we obtain, for $t=1,\ldots,2K$, that

[TABLE]

Because, we consider $2K$ marginal distributions, we can use Lemma 2 by setting $\alpha=(1-\varepsilon)$, $\tilde{\alpha}=(1-\tilde{\varepsilon})$, $a_{k}=\boldsymbol{A}_{k}[1,1](1-\varepsilon_{1})$, $\tilde{a}_{k}=\tilde{\boldsymbol{A}}_{k}[1,1](1-\tilde{\varepsilon}_{1})$, $w_{k}=\delta_{k}\pi_{k1}$ and $w_{k}=\tilde{\delta}_{k}\tilde{\pi}_{k1}$. Therefore, we have $\varepsilon=\tilde{\varepsilon}$, $\boldsymbol{A}_{k}[1,1]=\tilde{\boldsymbol{A}}_{k}[1,1]$ and $\delta_{k}\pi_{k1}=\tilde{\delta}_{k}\tilde{\pi}_{k1}$. Using the previous approach, with $\tau_{1}=\tau_{2}=1$ and $\tau_{3}<1$, with $h=2,\ldots,M$, we have for $t=1,\ldots,K$

[TABLE]

and thus $\boldsymbol{A}_{k}[1,h]=\tilde{\boldsymbol{A}}_{k}[1,h]$ and $\varepsilon_{h}=\tilde{\varepsilon}_{h}$. Similarly, taking $\tau_{2}<1$ and $\tau_{1}=\tau_{3}=1$, we have $\boldsymbol{A}_{k}[h,1]=\tilde{\boldsymbol{A}}_{k}[h,1]$. Finally, we have $\boldsymbol{A}_{k}[h,h^{\prime}]=\tilde{\boldsymbol{A}}_{k}[h,h^{\prime}]$ by increasing $h$ and $h^{\prime}$, by noting that with suitable choices of $\tau_{1}$, $\tau_{2}$ and $\tau_{3}$, we have for $t=1,\ldots,K$

[TABLE]

A.3 Proofs of the two technical lemmas

Proof of Lemma 1.

Since $a_{k}\neq a_{k^{\prime}}$ and $\tilde{a}_{k}\neq\tilde{a}_{k^{\prime}}$, then $N_{0}$, $N_{1}$, $\tilde{N}_{0}$ and $\tilde{N}_{1}$ are Vandermonde matrices and thus are invertible. Therefore, we have $w=N_{0}^{-1}\tilde{N}_{0}\tilde{w}=N_{1}^{-1}\tilde{N}_{1}\tilde{w},$ and thus

[TABLE]

or similarly for $u\in\{1,\ldots,K\}$,

[TABLE]

Since the previous equation holds for any $\tilde{w}$, we have $N_{0}^{-1}\tilde{N}_{0}=N_{1}^{-1}\tilde{N}_{1}$. Moreover, we have $N_{1}=N_{0}D$ and $\tilde{N}_{1}=\tilde{N}_{0}\tilde{D}$ where $D=\text{diag}(a_{1}^{K},\ldots,a_{K}^{K})$ and $\tilde{D}=\mathrm{diag}(\tilde{a}_{1}^{K},\ldots,\tilde{a}_{K}^{K})$. Denoting $R=N_{0}^{-1}\tilde{N}_{0}$, $DR=R\tilde{D}$ and then for $u\in\{1,\ldots,K\}$ and $k\in\{1,\ldots,K\}$

[TABLE]

We now show that $D=\tilde{D}$ and, $w=\tilde{w}$, and hence $R=I_{K}$ and $\tilde{N}_{0}=N_{0}$, where $I_{K}$ is the identity matrix of size $K$. First we show that $a_{1}=\tilde{a}_{1}$ and $w_{1}=\tilde{w}_{1}$.

• •

If $R[1,j]\neq 0$, (3) implies that $a_{1}^{K}R[1,j]=\tilde{a}_{j}^{K}R[1,j]$ and thus $a_{1}=\tilde{a}_{j}$. However, this is impossible because $a_{1}\geq\tilde{a}_{1}>a_{j}$ for $j\in\{2,\ldots,K\}$. Hence, we have $R[1,j]=0$ for $j=2,\ldots K$.

• •

Noting that $R$ is a product of two invertible matrices, $R$ is invertible. Therefore, $R[1,1]\neq 0$ because $R[1,j]=0$ for $j=2,\ldots K$. Hence, we have $a_{1}=\tilde{a}_{1}$.

• •

Note that $R[1,1]=\sum_{k=1}^{K}(N_{0}^{-1})[1,k]\tilde{N}_{0}[k,1]$ and that $\tilde{N}_{0}[k,1]=\tilde{a}_{1}^{k}=a_{1}^{k}=N_{0}[k,1]$. Therefore, we have $R[1,1]=\sum_{k=1}^{K}(N_{0}^{-1})[1,k]N_{0}[k,1]=(N_{0}^{-1}N_{0})[1,1]=1$.

• •

For $j=2,\ldots K$, $a_{1}>a_{j}$ so we have $R[j,1]=0$, because $a_{1}=\tilde{a}_{1}$.

• •

Because $w=R\tilde{w}$, we have $w_{1}=\tilde{w}_{1}$.

Equality of $a_{k}=\tilde{a}_{k}$ and $w_{k}=\tilde{w}_{k}$ can be shown recursively for $k=2,\ldots,K$ using the same reasoning. ∎

Appendix B Probabilities of misclassification

B.1 Technical lemmas

This section presents some notations and three lemmas which are used for the proof of Theorem 2. The technical lemmas discuss the concentration of the frequency of the observation $y_{i(t)}$ in a region of interest $W$, give an upper bound of $p(\boldsymbol{y}_{i}\mid Z_{ik}=1)$ and a concentration result of the ratio of $\frac{p(\tilde{\boldsymbol{x}}_{ik},\boldsymbol{y}_{i}\mid Z_{ik}=1)}{p(\tilde{\boldsymbol{x}}_{ik_{0}},\boldsymbol{y}_{i}\mid Z_{ik_{0}}=1)}$, where $\tilde{x}_{ik}=\operatorname*{argmax}_{\boldsymbol{x}_{i}\in\mathcal{X}}p(\boldsymbol{x}_{i},\boldsymbol{y}_{i}\mid Z_{ik}=1)$ is the estimator of the latent states conditionally on the observation $\boldsymbol{y}_{i}$ and on component $k$ obtained by applying the maximum a posteriori rule with the true parameter $\boldsymbol{\theta}$. The proof of the lemmas uses two concentration results for hidden Markov chains given by Kontorovich and Weiss, (2014) and by León and Perron, (2004).

Preliminaries Let $\boldsymbol{v}_{ik(t)}=(v_{ik(t)h\ell};h=1,\ldots,M;\ell=1,\ldots,M)$ with $\sloppy v_{ik(t)h\ell}=x_{ik(t-1)h}x_{ik(t)\ell}$ and $\tilde{\boldsymbol{v}}_{ik(t)}=(\tilde{v}_{ik(t)h\ell};h=1,\ldots,M;\ell=1,\ldots,M)$ with $\tilde{v}_{ik(t)h\ell}=\tilde{x}_{ik(t-1)h}\tilde{x}_{ik(t)\ell}$. In the following, $\mathbb{P}_{0}\left(\cdot\right)=\mathbb{P}\left(\cdot\mid Z_{ik_{0}}=1\right)$ by considering the true parameters.

Remark 1**.**

For any $k=1,\ldots,g$, $\boldsymbol{V}_{ik(t)}$ is a finite, ergodic and reversible Markov chain with $M^{2}$ states and transition matrix $\boldsymbol{P}_{k}$ with general term defined for any $(h_{1},h_{2},h_{3},h_{4})\in M^{4}$ by

[TABLE]

Moreover, the non-zero eigenvalues of $\boldsymbol{P}_{k}$ are the non-zero eigenvalues of $\boldsymbol{A}_{k}$ and the eigenvectors of $\boldsymbol{P}_{k}$ are obtained from the eigenvectors of $\boldsymbol{A}_{k}$.

Theorem 3** (Kontorovich and Weiss, (2014)).**

Let $U_{(1)},U_{(2)},\ldots$ be a stationary $\mathbb{N}$-valued $(G,\eta)$-geometrically ergodic Markov or hidden Markov chain, and consider the occupation frequency

[TABLE]

When $\sum_{u\in\mathbb{N}}\sqrt{\rho_{u}}<\infty$ with $\rho_{u}=\mathbb{P}(U_{(1)}=u)$, then for any $\varepsilon>0$

[TABLE]

where

[TABLE]

Theorem 4** (León and Perron, (2004)).**

For all pairs $(V,f)$, such that $V=(V_{(1)},\ldots,V_{(T)})$ is a finite, ergodic and reversible Markov chain in stationary state with the second-largest eigenvalue $\lambda$ and $f$ is a function taking values in $[0,1]$ such that $\mathbb{E}[f(V_{(t)})]<\infty$, the following bounds, with $\lambda_{0}=\max(0,\lambda)$, hold for all $s>0$ such that $\mathbb{E}[f(V_{(t)})]+s<1$ and all time $T$

[TABLE]

Concentration of the frequency of the observations in $W$ Let $W\subset\mathbb{R}^{+}$ be the subset of $\mathbb{R}^{+}$ where the estimator of $x_{i(t)}$ obtained by the maximum a posteriori rule is sensitive to $x_{i(t-1)}$ and $x_{i(t)}$ conditionally on $y_{i(t)}$ and component $k$. We define

[TABLE]

where

[TABLE]

and

[TABLE]

Lemma 3**.**

Let $\rho_{k_{0}}=\mathbb{P}_{0}(Y_{i(2)}\in W)$ and $\hat{\rho}_{k_{0}}=\sum_{t=1}^{T}\mathbf{1}_{\{y_{i(t)}\in W\}}$. For any $\delta_{1}>\frac{1}{\sqrt{2T}}$,

[TABLE]

$\delta_{1}=\varepsilon+\frac{1}{\sqrt{2T}}$* and $c_{1}=\frac{1}{2}(\delta_{1}-\frac{1}{\sqrt{2T}})^{2}>0$. Moreover, $\hat{\rho}_{k_{0}}$ is a consistent estimate of $\rho_{k_{0}}$ because the marginal distribution of $Y_{i(t)}$ is the same for any $t$, and thus $\rho_{k_{0}}=\mathbb{P}_{0}(Y_{i(t)}\in W)$ for any $t$.*

Proof of Lemma 3.

We have,

[TABLE]

Let $U_{(t)}=\mathbf{1}_{\{y_{i(t)}\in W\}}$. Then, for any $k=1,\ldots,g$, $U_{(1)},\ldots,U_{(T)}$ is a stationary $\{0,1\}$-valued $(1,0)$-geometrically ergodic hidden Markov chain conditionally on component $k$. Hence, by Theorem 3,

[TABLE]

We can conclude that

[TABLE]

$\delta_{1}=\varepsilon+\frac{1}{\sqrt{2T}}$ and $c_{1}=\frac{1}{2}(\delta_{1}-\frac{1}{\sqrt{2T}})^{2}$. ∎

Upper-bound of the conditional probability of $\boldsymbol{y}_{i}$ given $Z_{ik}=1$ Let $\gamma$ and $\bar{\gamma}$ be upper-bounds of the ratio $\frac{p(\tilde{x}_{i(t-1)},x_{i(t)},\tilde{x}_{i(t+1)},y_{i(t)}\mid Z_{ik}=1)}{p(\tilde{x}_{i(t-1)},\tilde{x}_{i(t)},\tilde{x}_{i(t+1)},y_{i(t)}\mid Z_{ik}=1)}$ when $y_{i(t)}\in W$ and $y_{i(t)}\notin W$ respectively. Thus, $\gamma=\max_{k}\max_{h_{1},h_{2},h_{3},h_{4}}\frac{A_{k}[h_{1},h_{2}]}{A_{k}[h_{3},h_{4}]}$ and $\bar{\gamma}$ permit to upper bound the ratio between the likelihood computed for any $(\boldsymbol{x}_{i},\boldsymbol{y}_{i})$ given $Z_{ik}=1$ and the likelihood computed with $(\tilde{\boldsymbol{x}}_{ik},\boldsymbol{y}_{i})$ given $Z_{ik}=1$. We have, if $y_{i(t)}\in W$,

[TABLE]

where $\mathbf{e}_{k}(u;h_{2})=\{(h_{1},h_{3}):h_{2}=\operatorname*{argmax}e_{k}(u;h_{1},h_{2},h_{3})\}$. Moreover, we have, if $y_{i(t)}\not\in W$,

[TABLE]

Note that $\gamma\geq 1$ and $\bar{\gamma}<1$.

Lemma 4**.**

We have, for any $k=1,\ldots,g$,

[TABLE]

where $c_{2}=1+\frac{\gamma}{1+\bar{\gamma}}$ and $\tilde{\gamma}=\max(2,\gamma)$.

Proof of Lemma 4.

By definition, we have

[TABLE]

Let $B_{p}(\tilde{\boldsymbol{x}}_{ik})=\{\boldsymbol{x}:\;\mid\mid\boldsymbol{x}-\boldsymbol{x}_{ik}\mid\mid_{0}=p\}$, then

[TABLE]

Remark that

[TABLE]

Moreover, we observe $T_{W}=T\hat{\rho}_{k_{0}}$ elements of the sequence $y_{i(1)},\ldots,y_{i(T)}$ which belongs to $W$. We have

[TABLE]

We have

[TABLE]

and

[TABLE]

where $\tilde{\gamma}=\max(2,\gamma)$. Noting that $1+\bar{\gamma}<\tilde{\gamma}+\bar{\gamma}$ and $1+\frac{\gamma}{1+\bar{\gamma}}>\frac{\tilde{\gamma}}{\bar{\gamma}+\tilde{\gamma}}$, we have

[TABLE]

where $c_{2}=1+\frac{\gamma}{1+\bar{\gamma}}$. Note that $\gamma+1>c_{2}>1$.

∎

Concentration of the ratio of complete-data likelihood

Lemma 5**.**

For any $k\neq k_{0}$ and for any $\delta_{3}$ such that $-\zeta<\delta_{3}<u_{kk_{0}}$, we have

[TABLE]

where $c_{3}=2\frac{1-\bar{\nu}_{2}(\boldsymbol{A}_{k_{0}})}{1+\bar{\nu}_{2}(\boldsymbol{A}_{k_{0}})}s^{2}>0$ and $s=\frac{\delta_{3}}{\omega_{kk_{0}}}+\frac{1}{\omega_{kk_{0}}}\sum_{h=1}^{M}\sum_{\ell=1}^{M}\pi_{k_{0}h}A_{k_{0}}[h,\ell]\log\left(\frac{A_{k_{0}}[h,\ell]}{A_{k}[h,\ell]}\right)$.

Proof of Lemma 5.

Let $f(\cdot)\in[0,1]$ defined by

[TABLE]

where $\omega_{kk_{0}}=u_{kk_{0}}+u_{k_{0}k}$, $u_{kk_{0}}=\max_{(h,\ell)}\log\left(\frac{A_{k}[h,\ell]}{A_{k_{0}}[h,\ell]}\right)$. Denoting $\mathbb{E}_{0}[\cdot]=\mathbb{E}[\cdot\mid Z_{ik_{0}}=1]$ the conditional expectation computed with the true parameters, we have, for $t=1,\ldots,T$,

[TABLE]

Therefore, we have

[TABLE]

where $s=\frac{\delta_{2}}{T\omega_{kk_{0}}}+\frac{1}{\omega_{kk_{0}}}\sum_{h=1}^{M}\sum_{\ell=1}^{M}\pi_{k_{0}h}A_{k_{0}}[h,\ell]\log\left(\frac{A_{k_{0}}[h,\ell]}{A_{k}[h,\ell]}\right)$.

Note that $\omega_{kk_{0}}>0$ and that, by Assumption 3, $\sum_{h=1}^{M}\sum_{\ell=1}^{M}\pi_{k_{0}h}A_{k_{0}}[h,\ell]\log\left(\frac{A_{k_{0}}[h,\ell]}{A_{k}[h,\ell]}\right)>\zeta>0$ because it is a weighted sum of $M$ Kullback-Leibler divergences. Thus, if $-T\zeta<\delta_{2}$ then $s>0$. Moreover, if $\delta_{2}<Tu_{kk_{0}}$, then $\mathbb{E}[f(\boldsymbol{V}_{i(1)})]+s<1$. Assumption 1 and Remark 1 imply that $\bar{\nu}_{2}(\boldsymbol{A}_{k_{0}})$ is the maximum between zero and the second-largest eigenvalue of reversible Markov chain of $\boldsymbol{V}_{i(t)}$. Therefore, using Theorem 4, we have for any $\delta_{3}$ such that $-\zeta<\delta_{3}<u_{kk_{0}}$,

[TABLE]

where $c_{3}=2\frac{1-\bar{\nu}_{2}(\boldsymbol{A}_{k_{0}})}{1+\bar{\nu}_{2}(\boldsymbol{A}_{k_{0}})}s^{2}$ and $s=\frac{\delta_{3}}{\omega_{kk_{0}}}+\frac{1}{\omega_{kk_{0}}}\sum_{h=1}^{M}\sum_{\ell=1}^{M}\pi_{k_{0}h}A_{k_{0}}[h,\ell]\log\left(\frac{A_{k_{0}}[h,\ell]}{A_{k}[h,\ell]}\right)$. ∎

B.2 Proof of Theorem 2

Noting that $\mathbb{P}(Z_{ik}=1\mid\boldsymbol{y}_{i})\propto\delta_{k}p(\boldsymbol{y}_{i}\mid Z_{ik}=1)$ and using Lemma 4, we have

[TABLE]

Moreover,

[TABLE]

where

[TABLE]

and

[TABLE]

Therefore, we have

[TABLE]

with $c_{4}=\log\frac{\delta_{k}}{a\delta_{k_{0}}}+\log\left(2M\tilde{\gamma}\gamma\max_{h,\ell}\frac{\pi_{kh}}{\pi_{k\ell}}\right)+\max_{k,k_{0},h,\ell}\log\frac{\pi_{kh}}{\pi_{k_{0}\ell}}$. By definition of $W$, we have $\tilde{v}_{ik(t)h\ell}=\tilde{v}_{ik_{0}(t)h\ell}$ if $y_{i(t)}\notin W$. Moreover, because $\tilde{\boldsymbol{v}}_{ik_{0}}$ is the maximum a posteriori rule, if $y_{i(t)}\in W$, then $d_{k1(t)}<\gamma$. Therefore, we have

[TABLE]

Hence, we have,

[TABLE]

Using Lemma 3, if $\delta_{1}>\frac{1}{\sqrt{2T}}$, the first term of the right side of the previous equation can be upper bounded by $e^{-Tc_{1}}$ with $c_{1}=\frac{1}{2}(\delta_{1}-\frac{1}{\sqrt{2T}})^{2}$.

Using Lemma 5, the second term of the right-hand side of the previous equation can be upper bounded by $e^{-Tc_{3}}$ with $c_{3}=2\frac{1-\bar{\nu}_{2}(\boldsymbol{A}_{k_{0}})}{1+\bar{\nu}_{2}(\boldsymbol{A}_{k_{0}})}s^{2}$, where $s=\frac{\delta_{3}}{\omega_{kk_{0}}}+\frac{1}{\omega_{kk_{0}}}\sum_{h=1}^{M}\sum_{\ell=1}^{M}\pi_{k_{0}h}A_{k_{0}}[h,\ell]\log\left(\frac{A_{k_{0}}[h,\ell]}{A_{k}[h,\ell]}\right)$ and $\delta_{3}=-\frac{c_{4}}{T}-\log(\tilde{\gamma}+\bar{\gamma})-(\gamma+c_{2})(\rho_{k_{0}}+\delta_{1})$, if $\delta_{3}$ is such that $-\zeta<\delta_{3}<u_{kk_{0}}$. Thus, we have the following condition on $\delta_{1}$

[TABLE]

Noting that $\gamma+c_{2}<1+2\gamma$, the previous upper bound can be satisfied under the following assumption

Assumption 4**.**

It holds that

[TABLE]

with $c_{4}=\log\frac{\delta_{k}}{a\delta_{k_{0}}}+\log\left(2M\tilde{\gamma}\gamma\max_{h,\ell}\frac{\pi_{kh}}{\pi_{k\ell}}\right)+\max_{k,k_{0},h,\ell}\log\frac{\pi_{kh}}{\pi_{k_{0}\ell}}$.

For any $a$ such that Assumption 4 holds and for any $\delta_{1}$ with $\frac{1}{\sqrt{2T}}<\delta_{1}<\frac{\zeta-c_{4}-\log(\tilde{\gamma}+\bar{\gamma})}{\gamma+c_{2}}-\rho_{k_{0}}$, we have

[TABLE]

and

[TABLE]

with $\delta_{3}=-\frac{c_{4}}{T}-\log(\tilde{\gamma}+\bar{\gamma})-(\gamma+c_{2})(\rho_{k_{0}}+\delta_{1})$. Therefore, there exists $c>0$ such that

[TABLE]

If the misclassification error is studied, we should consider $a=1$. Then, a sufficient condition to have the exponential decreasing of the probability of misclassifying an observation is obtained on the basis of Assumption 4 with $a=1$.

Appendix C Details about the conditional distribution

Forward formula We define

[TABLE]

which measures the probability of the partial sequence $y_{is(0)},\ldots,y_{is(t)}$ and ending up in state $h$ at time $t$ under component $k$. For any $(i,k,h,s)$, we can define $\alpha_{ikhs(t)}$ recursively, as follows,

[TABLE]

Considering independence between the $S_{i}$ sequences $\boldsymbol{y}_{is}$, the pdf of $\boldsymbol{y}_{i}$ under component $k$ is

[TABLE]

Therefore,

[TABLE]

Backward formula We define

[TABLE]

which measures the probability of the ending partial sequence $y_{is(t+1)},\ldots,y_{is(T_{is})}$ given a start in state $h$ at time $t$ under component $k$. We can define $\beta_{ikhs(t)}(\boldsymbol{\theta})$ recursively, for any $(i,k,h,s)$, as

[TABLE]

Considering independence between the $S_{i}$ sequences $\boldsymbol{y}_{is}$, the pdf of $\boldsymbol{y}_{i}$ under component $k$ is

[TABLE]

[TABLE]

Appendix D Additional simulation results

D.1 Additional tables for the analysis of simulated data

D.2 Additional figures for the analysis of the data from Huang et al., 2018b

D.3 Additional figures for the analysis of the PAT data