Unsupervised Learning of Mixture Regression Models for Longitudinal Data

TL;DR

This paper introduces a novel unsupervised penalized method for learning mixture regression models in longitudinal data, capable of selecting the number of components and estimating parameters for both continuous and discrete responses.

Contribution

It proposes a data-driven penalized approach with a modified EM algorithm for simultaneous model selection and parameter estimation in mixture regression models.

Findings

Method is consistent in component selection and parameter estimation.

Simulation studies show good finite sample performance.

Application to medical data demonstrates practical utility.

Abstract

This paper is concerned with learning of mixture regression models for individuals that are measured repeatedly. The adjective "unsupervised" implies that the number of mixing components is unknown and has to be determined, ideally by data driven tools. For this purpose, a novel penalized method is proposed to simultaneously select the number of mixing components and to estimate the mixing proportions and unknown parameters in the models. The proposed method is capable of handling both continuous and discrete responses by only requiring the first two moment conditions of the model distribution. It is shown to be consistent in both selecting the number of components and estimating the mixing proportions and unknown regression parameters. Further, a modified EM algorithm is developed to seamlessly integrate model selection and estimation. Simulation studies are conducted to evaluate the finite sample performance of the proposed procedure. And it is further illustrated via an analysis of a primary biliary cirrhosis data set.

Tables

Table 1. Table 1: The median and the 95 % percent \% confidence interval (CI) of total misclassification error rate in simulation studies. The values of median in Examples 1 and 2 are times 100. For the proposed PQL and PQL2 methods, the results below are summarized based on the models with correctly specified K 0 subscript 𝐾 0 K_{0} in 1000 replications.

(a)Example 1		Criterion	PQL	PQL2	QIFC
		median	0.000	0.000	0.058
		CI	(0.000, 0.000)	(0.000, 0.000)	(0.000, 0.000)
(b)Example 2		Criterion	PQL	PQL2	QIFC
	$\rho =0.3$	median	0.234	0.232	0.235
		CI	(0.000, 0.010)	(0.000, 0.010)	(0.000, 0.010)
	$\rho =0.6$	median	0.247	0.246	0.670
		CI	(0.000, 0.010)	(0.000, 0.010)	(0.000, 0.020)
(c)Example 3		Criterion	PQL	PQL2	QIFC
		median	0.209	0.209	0.214
		CI	(0.202, 0.218)	(0.202, 0.218)	(0.204, 0.226)

Table 2. Table 2: Estimation results in Example 1: (a) true values of mixture proportions and mixture parameters; (b) means of the parameter estimates; (c) means of the biases for the mixture proportions and mixture parameters; (d) mean squared errors (MSE) for the mixture proportions and mixture parameters. The values of bias and MSE are times 100. For the proposed PQL and PQL2 methods, the results below are summarized based on the models with correctly specified K 0 subscript 𝐾 0 K_{0} in 1000 replications.

Setting	${\beta }_{11}$	${\beta }_{12}$	${\beta }_{13}$	${\beta }_{14}$	${\beta }_{21}$	${\beta }_{22}$	${\beta }_{23}$	${\beta }_{24}$	${\sigma }_1^2$	${\sigma }_2^2$	${\pi }_1$	${\pi }_2$
$K_0=2$	True values
	0.08	-0.01	-0.4	0.06	-0.1	-0.05	3	0.3	0.5	0.8	0.5	0.5
	Mean
PQL	0.079	-0.010	-0.402	0.060	-0.099	-0.050	3.005	0.300	0.511	0.789	0.500	0.500
PQL2	0.080	-0.010	-0.402	0.060	-0.101	-0.050	3.006	0.300	0.494	0.787	0.500	0.500
QIFC	0.081	-0.010	-0.404	0.060	-0.101	-0.050	3.008	0.300	0.496	0.787	–	–
	Bias
PQL	-0.291	0.018	0.606	-0.006	-0.267	-0.008	-0.884	-0.002	-0.589	-1.090	0.047	-0.047
PQL2	-0.046	0.019	-0.209	0.018	-0.176	-0.009	0.816	-0.004	-0.557	-0.979	0.047	-0.047
QIFC	0.086	0.005	-0.109	0.008	-0.258	-0.009	0.990	-0.002	-0.548	-0.889	–	–
	MSE
PQL	0.621	0.001	4.485	0.000	1.088	0.000	3.413	0.006	2.388	0.210	0.021	0.021
PQL2	0.608	0.002	1.374	0.002	0.921	0.001	2.095	0.006	0.887	0.274	0.021	0.021
QIFC	0.558	0.001	1.281	0.000	0.949	0.001	2.131	0.002	0.098	0.253	–	–

Table 3. Table 3: Parameter estimates for primary biliary cirrhosis data.

	PQL		QIFC
Parameters	Group 0	Group 1	Group 0	Group 1
mixture proportions	0.512	0.487	–	–
	(0.129)	(0.129)	–	–
Trt	0.084	-0.097	0.055	-0.076
	(0.183)	(0.534)	(0.037)	(-0.113)
Age	-0.016	-0.051	-0.272	0.104
	(0.075)	(0.418)	(-0.261)	(0.262)
Sex	-0.366	3.220	-0.125	-0.204
	(0.097)	(1.219)	(-0.113)	(-0.064)
Time	0.068	0.313	0.093	-0.106
	(0.029)	(0.113)	(0.119)	(-0.108)
${\sigma }^2$	0.523	0.781	0.832	2.641
	(0.289)	(0.146)	(0.833)	(2.689)

Table 4. Table 4: Agreements and differences between the clinical and model classifications using the PQL and QIFC methods.

		PQL		QIFC		Total
Classify to		0	1	0	1
True	Group 0	118	54	69	103	172
	Group 1	42	98	21	119	140
	Total	160	152	90	222	312

Table 5. Table 5: Estimation results in Example 2: (a) true values of mixture proportions and mixture parameters; (b) means of the parameter estimates; (c) means of the biases for the mixture proportions and mixture parameters; (d) mean squared errors (MSE) for the mixture proportions and mixture parameters. The values of bias and MSE are times 100. For the proposed PQL and PQL2 methods, the results below are summarized based on the models with correctly specified K 0 subscript 𝐾 0 K_{0} in 1000 replications.

Setting	${\beta }_{10}$	${\beta }_{11}$	${\beta }_{12}$	${\beta }_{13}$	${\beta }_{20}$	${\beta }_{21}$	${\beta }_{22}$	${\beta }_{23}$	${\varphi }_1$	${\varphi }_2$	${\pi }_1$	${\pi }_2$
$K_0=2$	True values
	0	3	-1	1	4	-2	0	1	2	1	0.667	0.333
$\rho =0.3$	Mean
PQL	0.005	2.996	-1.005	0.995	4.001	-1.998	-0.001	0.998	1.983	0.972	0.672	0.328
PQL2	0.005	2.997	-1.006	0.995	4.001	-1.998	-0.001	0.998	2.000	0.989	0.672	0.328
QIFC	-0.013	3.014	-1.010	1.000	4.001	-2.000	-0.001	0.999	2.005	0.977	–	–
	Bias
PQL	0.584	-0.420	-0.531	-0.483	0.130	0.190	-0.329	-0.319	-1.726	-2.780	0.561	-0.561
PQL2	0.502	-0.296	-0.555	-0.456	0.149	0.180	-0.320	0.315	0.015	-1.085	0.561	-0.561
QIFC	-1.282	1.405	-0.997	0.033	0.117	0.006	-0.222	0.417	0.458	-2.254	–	–
	MSE
PQL	1.138	1.217	0.692	0.755	0.112	0.154	0.132	0.127	2.436	0.926	0.011	0.011
PQL2	1.029	1.082	0.603	0.703	0.101	0.135	0.117	0.112	2.531	0.892	0.011	0.011
QIFC	1.152	1.177	0.654	0.749	0.112	0.151	0.131	0.125	2.748	0.912	–	–
$\rho =0.6$	Mean
PQL	-0.002	2.997	-1.002	1.003	4.002	-2.000	-0.001	0.997	1.980	0.981	0.680	0.320
PQL2	-0.001	2.999	-1.001	0.998	4.002	-2.000	-0.001	0.998	2.003	1.003	0.680	0.320
QIFC	-0.023	3.020	-1.007	1.004	4.001	-2.002	-0.001	0.999	2.016	0.993	–	–
	Bias
PQL	-0.204	-0.298	-0.175	0.296	0.215	-0.033	-0.085	-0.289	-1.951	-1.878	1.334	-1.334
PQL2	-0.056	-0.131	-0.125	-0.155	0.213	-0.040	-0.088	-0.195	0.267	0.311	1.334	-1.334
QIFC	-2.346	2.003	-0.698	0.438	0.103	-0.244	-0.110	-0.102	1.557	-0.723	–	–
	MSE
PQL	1.275	1.174	0.729	0.774	0.118	0.149	0.121	0.140	3.763	1.386	0.031	0.031
PQL2	0.910	0.751	0.468	0.438	0.072	0.089	0.061	0.077	4.175	1.428	0.031	0.031
QIFC	1.101	0.907	0.518	0.491	0.083	0.100	0.071	0.089	4.489	1.410	–	–

Table 6. Table 6: Estimation results in Example 3: (a) true values of mixture proportions and mixture parameters; (b) means of the parameter estimates; (c) means of the biases for the mixture proportions and mixture parameters; (d) mean squared errors (MSE) for the mixture proportions and mixture parameters. The values of bias and MSE are times 100. For the proposed PQL and PQL2 methods, the results above are summarized based on the models with correctly specified K 0 subscript 𝐾 0 K_{0} in 1000 replications.

	True	PQL			PQL2			QIFC
	values	Mean	Bias	MSE	Mean	Bias	MSE	Mean	Bias	MSE
${\beta }_{10}$	2	1.992	-0.752	0.493	1.993	-0.698	0.368	1.996	-0.360	0.341
${\beta }_{11}$	1	0.998	-0.194	0.603	1.001	-0.014	0.327	1.000	0.107	0.357
${\beta }_{12}$	-1	-0.991	0.944	0.804	-0.989	1.113	0.439	-0.996	0.415	0.481
${\beta }_{13}$	1.5	1.496	-0.392	0.626	1.496	-0.436	0.296	1.499	-0.213	0.330
${\beta }_{14}$	1	0.998	-0.359	0.348	0.998	-0.234	0.179	0.999	-0.082	0.191
${\beta }_{20}$	-4	-3.988	1.219	0.277	-3.991	0.693	0.211	-3.999	0.071	0.215
${\beta }_{21}$	2	1.994	-0.621	0.388	1.995	-0.453	0.215	1.998	-0.350	0.227
${\beta }_{22}$	1	1.001	0.294	0.535	1.002	0.162	0.291	1.001	0.246	0.312
${\beta }_{23}$	-2	-2.001	-0.098	0.397	-2.000	0.019	0.211	-2.001	-0.112	0.226
${\beta }_{24}$	0	-0.002	-0.188	0.213	0.001	0.069	0.111	0.001	0.121	0.119
${\beta }_{30}$	-2	-1.998	0.249	0.140	-1.997	0.268	0.123	-1.999	0.235	0.135
${\beta }_{31}$	-2	-1.999	0.084	0.210	-2.000	0.046	0.169	-1.999	0.088	0.185
${\beta }_{32}$	1	0.999	-0.064	0.292	1.000	-0.009	0.245	1.001	-0.060	0.262
${\beta }_{33}$	0	-0.001	-0.065	0.199	0.000	-0.039	0.163	0.000	-0.081	0.181
${\beta }_{34}$	1	1.000	-0.090	0.106	1.000	-0.012	0.082	1.000	-0.008	0.090
${\beta }_{40}$	0	-0.006	-0.629	0.219	-0.008	-0.776	0.197	-0.001	-0.082	0.213
${\beta }_{41}$	1	0.998	-0.193	0.320	0.998	-0.194	0.260	0.999	-0.101	0.269
${\beta }_{42}$	0	-0.010	-0.988	0.417	-0.005	-0.496	0.340	0.001	0.070	0.358
${\beta }_{43}$	1	1.006	0.582	0.295	1.003	0.199	0.254	1.001	0.081	0.270
${\beta }_{44}$	1	1.001	0.140	0.163	1.002	0.174	0.135	1.001	0.109	0.144
${\beta }_{50}$	-4	-3.998	0.200	0.465	-4.001	0.001	0.463	-4.000	0.042	0.469
${\beta }_{51}$	0	-0.001	-0.050	0.883	0.000	-0.044	0.883	-0.005	-0.490	0.865
${\beta }_{52}$	-1	-0.992	0.806	1.319	-0.992	0.775	1.321	-0.995	0.882	1.265
${\beta }_{53}$	-1	-1.005	-0.416	0.958	-1.005	-0.338	0.965	-1.003	-0.338	0.936
${\beta }_{54}$	-1.5	-1.483	0.654	0.451	-1.493	0.651	0.454	-1.497	0.348	0.443
${\sigma }_1^2$	0.5	0.494	-0.592	0.189	0.493	-0.699	0.189	0.495	-0.578	0.150
${\sigma }_2^2$	0.3	0.292	-0.839	0.059	0.292	-0.824	0.059	0.297	-0.294	0.052
${\sigma }_3^2$	0.1	0.097	-0.251	0.008	0.098	-0.237	0.008	0.099	-0.137	0.008
${\sigma }_4^2$	0.15	0.143	-0.873	0.023	0.143	-0.735	0.023	0.147	-0.354	0.017
${\sigma }_5^2$	0.6	0.604	0.383	0.168	0.601	0.144	0.165	0.596	-0.396	0.146
${\pi }_1$	0.25	0.268	1.797	0.054	0.268	1.797	0.054	–	–	–
${\pi }_2$	0.25	0.264	1.418	0.037	0.264	1.418	0.037	–	–	–
${\pi }_3$	0.15	0.132	-1.791	0.050	0.132	-1.791	0.050	–	–	–
${\pi }_4$	0.15	0.132	-1.796	0.054	0.132	-1.796	0.054	–	–	–
${\pi }_5$	0.2	0.204	0.373	0.003	0.204	0.373	0.003	–	–	–

Equations

$${\mu }_{ijk}\equiv \text{E}(Y_{ij}\mid X_{ij},u_i=k)=g(X_{ij}^T{\beta }_k),$$

$${\sigma }_{ijk}^2\equiv \text{var}(Y_{ij}\mid X_{ij},u_i=k)={\varphi }_{k}V({\mu }_{ijk}),$$

$$Y_{ij}\mid X_{ij}\sim \sum _{k=1}^K{\pi }_k{f}_k(Y_{ij}\mid X_{ij}^T{\beta }_k,{\varphi }_k),$$

$$Q(\theta )=\sum _{i=1}^n\log \left[\sum _{k=1}^K{\pi }_k\exp \left\{\sum _{j=1}^{m_i}q(g(X_{ij}^T{\beta }_k);Y_{ij})\right\}\right],$$

$$Q_{\text{P}}(\theta )=Q(\theta )-n\lambda \sum _{k=1}^K\{\log (ϵ+{\pi }_k)-\log (ϵ)\},$$

$$\Psi (\theta ;Y_i|X_i)=\sum _{k=1}^K{\pi }_k\exp \left\{\sum _{j=1}^{m_i}q(g(X_{ij}^T{\beta }_k);Y_{ij})\right\}$$

$$\sqrt{n}(\hat{\theta }-{\theta }_0)\stackrel{L}{⟶}N(0,B^{-1}A{B}^{-1}).$$

$$\Psi (\eta ,\gamma ;Y_i\mid X_i)=\sum _{l=1}^{K-K_0}{\lambda }_l\eta f({\beta }_l;Y_i\mid X_i)+\sum _{k=1}^{K_0}({\pi }_{0k}+{\rho }_k\eta )f({\beta }_{0k}+\eta {\delta }_k;Y_i\mid X_i),$$

$$f({\beta }_l;Y_i\mid X_i)=\exp \left\{\sum _{j=1}^{m_i}q(g(X_{ij}^T{\beta }_l);Y_{ij})\right\}$$

$$\gamma =({\lambda }_1,\dots ,{\lambda }_{K-K_0},{\rho }_1,\dots ,{\rho }_{K_0},{\beta }_1^T,\dots ,{\beta }_{K-K_0}^T,{\delta }_1^T,\dots ,{\delta }_{K_0}^T{)}^T$$

$$Q_{\text{P}}(\eta ,\gamma )\equiv \sum _{i=1}^n\log \{\Psi (\eta ,\gamma ;Y_i\mid X_i)\}-n\lambda \sum _{k=1}^K\{\log (ϵ+{\pi }_k)-\log (ϵ)\}.$$

$$\sum _{i=1}^n\sum _{k=1}^K{u}_{ik}\left\{\log {\pi }_k+\sum _{j=1}^{m_i}q(g(X_{ij}^T{\beta }_k);Y_{ij})\right\}-n\lambda \sum _{k=1}^K\{\log (ϵ+{\pi }_k)-\log (ϵ)\}.$$

$${\hat{u}}_{ik}^{(t+1)}=\frac{{\pi }_k^{(t)}\exp \left\{{\sum }_{j=1}^{m_i}\tilde{q}(g(X_{ij}^T{\beta }_k^{(t)}),{\varphi }_k^{(t)};Y_{ij})\right\}}{{\sum }_{l=1}^K{\pi }_l^{(t)}\exp \left\{{\sum }_{j=1}^{m_l}\tilde{q}(g(X_{ij}^T{\beta }_l^{(t)}),{\varphi }_l^{(t)};Y_{ij})\right\}},$$

$$\sum _{i=1}^n\sum _{k=1}^K{\hat{u}}_{ik}^{(t+1)}\left\{\log {\pi }_k+\sum _{j=1}^{m_i}q(g(X_{ij}^T{\beta }_k);Y_{ij})\right\}-n\lambda \sum _{k=1}^K\{\log (ϵ+{\pi }_k)-\log (ϵ)\}.$$

$$\frac{\partial }{\partial \pi }\left\{\sum _{i=1}^n\sum _{k=1}^K{\hat{u}}_{ik}^{(t+1)}\log {\pi }_k-n\lambda \sum _{k=1}^K\{\log (ϵ+{\pi }_k)-\log (ϵ)\}-\xi (\sum _{k=1}^K{\pi }_k-1)\right\}=0,$$

$${\pi }_k^{(t+1)}=\max \left\{0,\frac{1}{1-\lambda K}\left[\frac{1}{n}\sum _{i=1}^n{\hat{u}}_{ik}^{(t+1)}-\lambda \right]\right\},k=1,\dots ,K.$$

$$\sum _{i=1}^n\sum _{j=1}^{m_i}{\hat{u}}_{ik}^{(t+1)}{g}^{\prime }(X_{ij}^T{\beta }_k)X_{ij}\frac{Y_{ij}-g(X_{ij}^T{\beta }_k)}{V(g(X_{ij}^T{\beta }_k))}=0,$$

$${\varphi }_k^{(t+1)}=\sum _{i=1}^n\frac{{\hat{u}}_{ik}^{(t+1)}}{{\sum }_{i^{\prime }=1}^n{m}_{i^{\prime }}{\hat{u}}_{i^{\prime }k}^{(t+1)}}\sum _{j=1}^{m_i}\frac{\{Y_{ij}-g(X_{ij}^T{\beta }_k^{(t)}){\}}^2}{V(g(X_{ij}^T{\beta }_k^{(t)}))},k=1,\dots ,K.$$

$$\text{BIC}(\lambda )=-2\sum _{i=1}^n\log \left[\sum _{k=1}^{\hat{K}}{\hat{\pi }}_k\exp \left\{\sum _{j=1}^{m_i}\tilde{q}(g(X_{ij}^T{\hat{\beta }}_k),{\hat{\varphi }}_k;Y_{ij})\right\}\right]+\hat{K}(p+2)\log n,$$

$$k^{\ast }=\arg \underset{1\le k\le \hat{K}}{\max }{\hat{\pi }}_k\exp \left\{\sum _{j=1}^m\tilde{q}(g(X_{ij}^{\ast T}{\hat{\beta }}_k),{\hat{\varphi }}_k;Y_{ij}^{\ast })\right\}.$$

$$E(Y_{ij})={\beta }_{k1}{X}_{i1}+{\beta }_{k2}{X}_{i2}+{\beta }_{k3}{X}_{i3}+{\beta }_{k4}{X}_{ij4},$$

$$\log ({\mu }_{ij1})=3X_{ij1}-X_{ij2}+X_{ij3},i=1,\dots ,n_1,j=1,\dots ,m_i,$$

$$\log ({\mu }_{ij2})=4-2X_{ij1}+X_{ij3},i=1,\dots ,n_2,j=1,\dots ,m_i,$$

$$Y_i\mid u_i=k\sim MVN({\beta }_{k0}+X_i{\beta }_k,{\sigma }_k^2{R}_i^{(k)}),k=1,\cdots ,5$$

$$\text{E}(Y_{ij})={\beta }_{k1}{\text{Trt}}_{ij}+{\beta }_{k2}{\text{Age}}_{ij}+{\beta }_{k3}{\text{Sex}}_{ij}+{\beta }_{k4}{\text{Time}}_{ij},$$

$$Q_n({\hat{\theta }}^{\ast })=\sum _{i=1}^n\{\psi (n^{-1/2}{\hat{\theta }}^{\ast }+{\theta }_0;Y_i\mid X_i)-\psi ({\theta }_0;Y_i\mid X_i)\}.$$

$$Q_n({\hat{\theta }}^{\ast })$$

$${\hat{\theta }}^{\ast }\stackrel{L}{⟶}N(0,B^{-1}A{B}^{-1}).$$

$${\Psi }_0(Y\mid X)=\sum _{k=1}^{K_0}{\pi }_{0k}f({\beta }_{0k};Y\mid X).$$

$$d(\gamma ;Y\mid X)=\sum _{k=1}^{K_0}{\pi }_{0k}\sum _{j=1}^p{\delta }_{kj}\frac{D_j^{1}f({\beta }_{0k};Y\mid X)}{{\Psi }_0(Y\mid X)}+\sum _{l=1}^{K-K_0}{\lambda }_l\frac{f({\beta }_l;Y\mid X)}{{\Psi }_0(Y\mid X)}+\sum _{k=1}^{K_0}{\rho }_k\frac{f({\beta }_{0k};Y\mid X)}{{\Psi }_0(Y\mid X)},$$

Full text

Unsupervised Learning of Mixture Regression Models for Longitudinal Data

Peirong Xu1, Heng Peng2, Tao Huang3111Corresponding Author: Tao Huang, Associate Professor, School of Statistics and Management, Shanghai University of Finance and Economics, Shanghai, China, 200433. (E-mail:huang.tao@mail.shufe.edu.cn).

1 College of Mathematics and Sciences, Shanghai Normal University, Shanghai, China

2 Department of Mathematics, Hong Kong Baptist University, Hong Kong, China

3 School of Statistics and Management, Shanghai University of Finance and Economics, Shanghai, China

Abstract: This paper is concerned with learning of mixture regression models for individuals that are measured repeatedly. The adjective “unsupervised” implies that the number of mixing components is unknown and has to be determined, ideally by data driven tools. For this purpose, a novel penalized method is proposed to simultaneously select the number of mixing components and to estimate the mixture proportions and unknown parameters in the models. The proposed method is capable of handling both continuous and discrete responses by only requiring the first two moment conditions of the model distribution. It is shown to be consistent in both selecting the number of components and estimating the mixture proportions and unknown regression parameters. Further, a modified EM algorithm is developed to seamlessly integrate model selection and estimation. Simulation studies are conducted to evaluate the finite sample performance of the proposed procedure. And it is further illustrated via an analysis of a primary biliary cirrhosis data set.

Key words: Unsupervised learning, Model selection, Longitudinal data analysis, Quasi-likelihood, EM algorithm.

1 Introduction

In many medical studies, the marker of disease progression and a variety of characteristics are routinely measured during the patients’ follow-up visit to decide on future treatment actions. Consider a motivating Mayo Clinic trial with primary biliary cirrhosis (PBC), wherein a number of serological, clinical and histological parameters were recorded for each of 312 patients from 1974 to 1984. This longitudinal study had a median follow-up time of 6.3 years as some patients missed their appointments due to worsening medical condition of some labs. It is known that PBC is a fatal chronic cholesteric liver disease, which is characterized histopathologically by portal inflammation and immune-mediated destruction of the intrahepatic bile ducts (Pontecorvo, Levinson, and Roth, 1992). It can be divided into four histologic stages, but with nonuniformly affected liver. The diagnosis of PBC is important for the medical treatment with Ursodiol has been shown to halt disease progression and improve survival without need for liver transplantation (Talwalkar and Lindor, 2003). Therefore, one goal of the study was the investigation of the serum bilirubin level, an important marker of PBC progression, in relation to the time and to potential clinical and histological covariates. Another issue that should be accounted for is the unobservable heterogeneity between subjects that may not be explained by the covariates. The changes in inflammation and bile ducts occur at different rates and with varying degrees of severity in different patients, so the heterogeneous patients could potentially belong to different latent groups. To address these problems, there is a demand for mixture regression modeling for subjects on the basis of longitudinal measurements.

There are various research works on mixture regression models for longitudinal outcome data, particularly in the context of model-based probabilistic clustering (Fraely and Raftery, 2002). For example, De la Cruz-Mes$\acute{i}$a et. al. (2008) proposed a mixture of non-linear hierarchical models with Gaussian subdistributions; McNicholas and Murphy (2010) extended the Gaussian mixture models with Cholesky-decomposed group covariance structure; Kom$\acute{a}$rek and Kom$\acute{a}$rkov$\acute{a}$ (2013) introduced a generalized linear mixed model for components’ densities under the Gaussian mixture framework; Heinzl and Tutz (2013) considered linear mixed models with approximate Dirichlet process mixtures. Other relevant work includes Celeux et. al. (2005), Booth et. al. (2008), Pickles and Croudace (2010), Maroutti (2011), Erosheva et. al. (2014) and some of the references therein. Compared with heuristic methods such as the k-means method (Genolini and Falissard, 2010), issues like the selection of the number of clusters (or components) can be addressed in a principled way. However, most of them assume a parametric mixture distribution, which may be too restrictive and invalid in practice when the true data-generating mechanism indicates otherwise.

A key concern for the performance of mixture modeling is the selection of the number of components. A mixture with too many components may overfit the data and result in poor interpretations. Many statistical methods have been proposed in the past few decades by using the information criteria. For example, see Leroux (1992), Roeder and Wasserman (1997), Hennig(2004), De la Cruz-Mes$\acute{i}$a et al. (2008) and many others. However, these methods are all based on the complete model search algorithm, which result in heavy computation burden. To improve the computational efficiency, data-driven procedures are much more preferred. Recently, Chen and Khalili (2008) used the SCAD penalty (Fan and Li, 2001) to penalize the difference of location parameters for mixtures of univariate location distributions; Kom$\acute{a}$rek and Lesaffre (2008) suggested to penalize the reparameterized mixture weights in the generalized mixed model with Gaussian mixtures; Heinzl and Tutz (2014) constructed a group fused lasso penalty in linear-mixed models; Huang et. al. (2016) proposed a penalized likelihood method in finite Gaussian mixture models. Most of them are developed for independent data or based on the full likelihood. However, the full likelihood is often difficult to specify in formulating a mixture model for longitudinal data, particularly for correlated discrete data.

Instead of specifying the form of distribution of the observations, a quasi-likelihood method (Wedderburn, 1974) gives consistent estimates of parameters in mixture regression models that only needs the relation between the mean and variance of each observation. Inspired by its nice property, in this paper, we propose a new penalized method based on quasi-likelihood for mixture regression models to deal with the above mentioned problems simultaneously. This would be the first attempt to handle both balanced and unbalanced longitudinal data that only requires the first two moment conditions of the model distribution. By penalizing the logarithm of mixture proportions, our approach can simultaneously select the number of mixing components and estimate the mixture proportions and unknown parameters in the semiparametric mixture regression model. The number of components can be consistently selected. And given the number of components, the estimators of mixture proportions and regression parameters can be root-$n$ consistent and asymptotically normal. By taking account of the within-component dispersion, we further develop a modified EM algorithm to improve the classification accuracy. Simulation results and the application to the motivating PBC data demonstrate the feasibility and effectiveness of the proposed method.

The rest of the paper is organized as follows. In Section 2, we introduce a new penalized method for learning semiparametric mixture regression models with longitudinal data. Section 3 presents the corresponding theoretical properties and Section 4 provides a modified EM algorithm for implementation. In Section 5, we assess the finite sample performance of the proposed method via simulation studies. We apply the proposed method to the PBC data in Section 6, and conclude the paper with Section 7. All technical proofs are provided in Appendix.

2 Learning semiparametric mixture of regressions

2.1 Model specification

In a longitudinal study, suppose $Y_{ij}$ is the response variable measured at the $j$th time point for the $i$th subject, and $X_{ij}$ is the corresponding $p\times 1$ vector of covariates, $i=1,\ldots,n$, $j=1,\ldots,m_{i}$. Let $Y_{i}=(Y_{i1},\ldots,Y_{im_{i}})^{T}$ and $X_{i}=(X_{i1},\ldots,X_{im_{i}})^{T}$. In general, the observations for different subjects are independent, but they may be correlated within the same subject. We assume that the observations of each subject belong to one of $K$ classes (components) and $u_{i}\in\{1,\ldots,K\}$ is the corresponding latent class variable. Assume that $u_{i}$ has a discrete distribution ${\rm{P}}(u_{i}=k)=\pi_{k}$, where $\pi_{k},k=1,\ldots,K$, are the positive mixture proportions satisfying $\sum^{K}_{k=1}\pi_{k}=1$. Given $u_{i}=k$ and $X_{ij}$, suppose the conditional mean of $Y_{ij}$ is

[TABLE]

where $g$ is a known link function, and $\beta_{k}$ is a $p-$dimensional unknown parameter vector. The corresponding conditional variance of $Y_{ij}$ is given by

[TABLE]

where $V$ is a known positive function and $\phi_{k}$ is a unknown dispersion parameter. In other words, conditioning on $X_{ij}$, the response variable $Y_{ij}$ follows a mixture distribution

[TABLE]

where $f_{k}(Y_{ij}\mid X^{T}_{ij}\beta_{k},\phi_{k})$’s are the components’ distributions. To avoid identifiability issues, we assume that $K$ is the smallest integer such that $\pi_{k}>0$ for $k=1,\cdots,K$, and $(\beta_{a},\phi_{a})\neq(\beta_{b},\phi_{b})$ for $1\leq a<b\leq K$. Denote $\theta=(\beta^{T}_{1},\ldots,\beta^{T}_{K},\phi^{T},\pi^{T})^{T}$ with $\beta_{k}=(\beta_{k1},\ldots,\beta_{kp})^{T}$ and $\pi=(\pi_{1},\ldots,\pi_{K-1})^{T}$, and $\phi=(\phi_{1},\ldots,\phi_{K})^{T}$.

Under the working independence correlation, the (log) quasi-likelihood of the $K$-component marginal mixture regression model is

[TABLE]

where function $q(\mu;y)$ (McCullagh and Nelder, 1989) satisfies $\frac{\partial q(\mu;y)}{\partial\mu}=\frac{y-\mu}{V(\mu)}$. It is known that, for a generalized linear model with independent data, the quasi-likelihood estimator of the regression coefficient has the same asymptotic properties as the maximum likelihood estimator. While for longitudinal data, it is equivalent to the GEE estimator (Liang and Zeger, 1986), which is consistent even when the working correlation structure is misspecified. Therefore, estimation consistency is expected to hold for the $K$-component marginal mixture regression model (2.1)-(2.2), and this will be validated in Section 3.

2.2 Penalized quasi-likelihood method

For a fixed number of $K$ components, we can maximize the quasi-likelihood function (2.3) by an expectation-maximization (EM) algorithm, which in the E-step computes the posterior probability of the class memberships and in the M-step estimates the mixture proportions and unknown parameters. However, in practice, the number of components is usually unknown and needs to be inferred from the data itself.

For the proposed marginal mixture regression model, the selection of the number of mixing components can be viewed as a model selection problem. Various conventional methods have been proposed based on the likelihood function and some information theoretic criteria. In particular, the Bayesian information criterion (BIC; Schwarz, 1978) is recommended as a useful tool for selecting the number of components (Dasgupta and Raftery, 1998; Fraley and Rafetery, 2002). Therefore, a natural idea is to propose a BIC-type criterion for selecting the number of mixing components, where the likelihood function is replaced by the quasi-likelihood function (2.3). But our simulation experience shows that it couldn’t perform as well as the traditional BIC, since (2.3) is no longer a joint density with integral equals to one.

To avoid calculating the normalizing constant, the penalization technique is preferred. By (2.3), intuitively, the $k$th component would be eliminated if $\pi_{k}=0$. But in implementation of (2.3), the quasi-likelihood function for the complete data $(u_{ik},Y_{i},X_{i})$ involves $\log{\pi_{k}}$ rather than $\pi_{k}$, where $u_{ik}$ denotes the indicator of whether $i$th subject belongs to the $k$th component (see (4.1) defined in Section 4 for details). Therefore, it is natural to penalize the logarithm of mixture proportions $\log\pi_{k}$, $k=1,\ldots,K$. Moreover, note that the gradient of $\log\pi_{k}$ increases very fast when $\pi_{k}$ is close to zero, and it would dominate the gradient of nonzero $\pi_{l}>0$. Consequently, the popular $L_{q}$ types of penalties may not able to set insignificant $\pi_{k}$ to zero. In the spirit of penalization in Huang et al. (2016), we propose the following penalized quasi-likelihood function

[TABLE]

where $\lambda$ is a tuning parameter and $\epsilon$ is a very small positive constant. Note that $\log(\epsilon+\pi_{k})-\log(\epsilon)$ is an increasing function of $\pi_{k}$ and is shrunk to zero as the mixing proportion $\pi_{k}$ goes to zero. Therefore, the proposed method (2.4) can simultaneously determine the number of mixture components and estimate mixture proportions and unknown parameters.

Remark 1**.**

The small constant $\epsilon$ is introduced to ensure the continuity of the objective function when some of mixture proportions are shrunk continuously to zero.

Remark 2**.**

The penalty $n\lambda\sum^{K}_{k=1}\{\log(\epsilon+\pi_{k})-\log(\epsilon)\}$ in (2.4) would over penalize large $\pi_{k}$ and result in a biased estimator. A more general but slightly more complicated approach is to use $n\sum^{K}_{k=1}\{\log(\epsilon+p_{\lambda}(\pi_{k}))-\log(\epsilon)\}$, where $p_{\lambda}(\cdot)$ is a penalty function that gives estimators with sparsity, unbiasedness and continuity as discussed in Fan and Li (2001).

3 Asymptotic properties

In this section, we first study the asymptotic property of the maximum quasi-likelihood estimator $\widehat{\theta}$ of (2.3) given the number of mixing components. And then, we establish the model selection consistency of the proposed method (2.4) for the general semiparametric marginal mixture regression model (2.1)-(2.2).

For a fixed number of $K$ components, denote the true value of parameter vector by $\theta_{0}$. The components of $\theta_{0}$ are denoted with a subscript, such as $\pi_{0k}$. We assume that the number of subjects $n$ increases to infinity, while the number of observations $\{m_{i}\}$ is a bounded sequence of positive integers. Let

[TABLE]

and $\psi(\theta;Y_{i}\mid X_{i})=\log(\Psi(\theta;Y_{i}\mid X_{i}))$.

We assume the following regularity conditions to derive the asymptotic properties.

1. C1

The function $g(\cdot)$ has two bounded and continuous derivatives. 2. C2

The random variables $X_{ij}$’s are bounded on the compact support $\mathcal{A}$ uniformly. For $\theta\in\Omega$, the density function of $X^{T}_{ij}\beta_{k}$ is positive and satisfies Lipschitz condition of order 1 on $\mathcal{U}_{k}=\{u=X^{T}_{ij}\beta_{k}:X_{ij}\in\mathcal{A},i=1,\ldots,n,j=1,\ldots,m_{i}\}$, $k=1,\ldots,K$. 3. C3

$\Omega$ is compact and $\theta_{0}$ is an interior point in $\Omega$. 4. C4

For each $\theta\in\Omega$, $\psi(\theta;Y_{i}\mid X_{i})$ admits third order partial derivatives with respect to $\theta$. And there exist functions $M_{l}(X_{i},Y_{i})$, $l=0,1,2,3$ such that for $\theta$ in a neighborhood of $\theta_{0}$, $|\psi(\theta;Y_{i}\mid X_{i})-\psi(\theta_{0};Y_{i}\mid X_{i})|\leq M_{0}(X_{i},Y_{i})$, $|\partial\psi(\theta;Y_{i}\mid X_{i})/\partial\theta_{j}|\leq M_{1}(X_{i},Y_{i})$, $|\partial^{2}\psi(\theta;Y_{i}\mid X_{i})/\partial\theta_{j}\partial\theta_{k}|\leq M_{2}(X_{i},Y_{i})$, and $|\partial^{3}\psi(\theta_{0};Y_{i}\mid X_{i})/\partial\theta_{j}\partial\theta_{k}\partial\theta_{l}|\leq M_{3}(X_{i},Y_{i})$ with $E\{M_{l}(X_{i},$ $Y_{i})\}<\infty$, for all $i=1,\ldots,n$. 5. C5

$\theta_{0}$ is the identifiably unique maximizer of $E\{Q(\theta)\}$. 6. C6

Let $A=\textup{var}\{\partial\psi(\theta_{0};Y_{i}\mid X_{i})/\partial\theta\}$. The second derivative matrix $B=\textup{E}\{-\partial^{2}\psi(\theta_{0};Y_{i}\mid X_{i})/\partial\theta\partial\theta^{T}\}$ is positive definite.

Conditions C1-C2 are typical assumptions in the estimation literature, which are also found in Xu and Zhu (2012) and Xu et. al. (2016). Conditions C3-C6 are mild conditions in the literature of mixture models, which are used for the proof of weak consistency and asymptotic normality.

Theorem 1**.**

Under conditions C1-C6, the maximum quasi-likelihood estimator $\widehat{\theta}$ of (2.3) given the number of components is consistent and has the asymptotic normality

[TABLE]

Next, we study the model selection consistency of the proposed method (2.4) for the marginal mixture regression model (2.1)-(2.2). We assume that there are $K_{0}$ mixture components, $K_{0}\leq K$ with $\pi_{l}=0$, for $l=1,\ldots,K-K_{0}$ and $\pi_{l}=\pi_{0k}$ for $l=K-K_{0}+1,\ldots,K$, $k=1,\ldots,K_{0}$. In the spirit of locally conic parametrization (Dacunha-Castelle and Gassiat, 1997), define $\pi_{l}=\lambda_{l}\eta$, $l=1,\ldots,K-K_{0}$ and $\pi_{l}=\pi_{0k}+\rho_{k}\eta$, $l=K-K_{0}+1,\ldots,K$, $k=1,\ldots,K_{0}$. Then, the function (3.1) can be rewritten as

[TABLE]

where

[TABLE]

and

[TABLE]

with restrictions $\lambda_{l}\geq 0,\beta_{l}\in R^{p},l=1,\ldots,K-K_{0},$ $\delta_{k}\in R^{p},\rho_{k}\in R,k=1,\ldots,K_{0},$ $\sum^{K-K_{0}}_{l=1}\lambda_{l}+\sum^{K_{0}}_{k=1}\rho_{k}=0$ and $\sum^{K-K_{0}}_{l=1}\lambda^{2}_{l}+\sum^{K_{0}}_{k=1}\rho^{2}_{k}+\sum^{K_{0}}_{k=1}\|\delta_{k}\|^{2}=1$. By the permutation, such a parametrization is locally conic and identifiable. And then, the penalized quasi-likelihood function (2.4) can be rewritten as

[TABLE]

To establish the model selection consistency of the proposed method, we need the following additional conditions:

1. C7

There exists a positive constant $\varepsilon$ such that $g(X^{T}_{ij}\beta_{k})$ and $V(g(X^{T}_{ij}\beta_{k}))$ are bounded on $\mathcal{B}=\{\beta:\|\beta-\beta_{0}\|\leq\varepsilon\}$ uniformly in $i=1,\ldots,n$, $j=1,\ldots,m_{i}$, $k=1,\ldots,K$. 2. C8

Let $\Sigma_{ik}=\textup{cov}(Y_{i}\mid X_{i},u_{i}=k)$, and $V_{ik}$ be a $m_{i}\times m_{i}$ diagonal matrix with $j$th element $\sigma^{2}_{ijk}$. Both the eigenvalues of $\Sigma_{ik}$ and $V_{ik}$ are uniformly bounded away from 0 and infinity.

Condition C7 is analogous to conditions (A2) and (A6) in Wang (2011), which is generally satisfied for marginal models. For example, when the marginal model follows a Poisson regression, $V(g(X^{T}_{ij}\beta_{k}))=g(X^{T}_{ij}\beta_{k})=\exp(X^{T}_{ij}\beta_{k})$’s are uniformly bounded around on $\mathcal{B}$. Condition C8 is similar to conditions (C3) and (C4) in Huang et. al. (2007), which ensures the non-singularity of the covariance matrices and the working covariance matrices.

Theorem 2**.**

Under conditions C1-C8, if $\lim_{n\rightarrow\infty}\sqrt{n}\lambda=a$ and $\epsilon=o(n^{-1/2}/\log{n})$, where $a$ is a constant, there exists a local maximizer $(\eta,\gamma)$ of (3.2) such that $\eta=O_{p}(n^{-1/2})$, and for such local maximizer, we have $\widehat{K}\stackrel{{\scriptstyle P}}{{\longrightarrow}}K_{0}$.

Theorem 2 indicates that by choosing an appropriate tuning parameter $\lambda$ and a small constant $\epsilon$, the proposed method (2.4) can select the number of mixing components consistently.

4 Implementation and tuning parameter selection

In this section, we propose an algorithm to implement the proposed method (2.4) and a procedure to select the tuning parameter $\lambda$.

4.1 Modified EM Algorithm

Since the membership of each subject is unknown, it is natural to use EM algorithm to implement (2.4). But note that the criterion (2.4) is a function unrelated to different dispersion parameters $\phi_{k},k=1,\ldots,K$, therefore, the naive EM algorithm may decrease the classification accuracy for the observations by ignoring the within-component dispersion. Therefore, we here propose a modified EM algorithm in consideration of different component dispersion.

Let $u_{ik}$ denote the indicator of whether the $i$th subject is in the $k$th class. That is, $u_{ik}=1$ if the $i$th subject belongs to the $k$th component, and $u_{ik}=0$ otherwise. If the missing data $\{u_{ik},i=1,\ldots,n,k=1,\ldots,K\}$ were observed, the penalized quasi-likelihood function for the complete data is given by

[TABLE]

Denote $\Theta=(\pi^{T},\beta^{T},\phi^{T})^{T}$ as the vector of all parameters in the $K$-component marginal mixture regression model (2.1)-(2.2) with $\beta=(\beta^{T}_{1},\ldots,\beta^{T}_{K})^{T}$. In the E-step, given the current estimate $\Theta^{(t)}=(\pi^{(t)T},\beta^{(t)T},\phi^{(t)T})^{T}$, we impute values for the unobserved $u_{ik}$ by

[TABLE]

where $\tilde{q}(\mu,\phi;y)=\int^{\mu}_{y}\frac{y-t}{\phi V(t)}dt.$ Plugging them into (4.1), we obtain the function

[TABLE]

In the M-step, the goal is to update $\pi^{(t)}$ and $\beta^{(t)}$ by maximizing (4.2) with the constraint $\sum^{K}_{k=1}\pi_{k}=1$ and update $\phi^{(t)}$ by the residual moment method. Specifically, to update $\pi^{(t)}$, we solve the following equations

[TABLE]

where $\xi$ is the Lagrange multiplier. Then, when $\epsilon$ is very close to zero, it gives

[TABLE]

$\beta^{(t)}_{k}$ can be updated by solving the following equations

[TABLE]

where $g^{\prime}(\cdot)$ is the first derivative of $g$, $k=1,\ldots,K$. And using the residual moment method, we update $\phi^{(t)}$ as follows

[TABLE]

Remark 3**.**

In the initial step, we pre-specify a large number of components, and once a mixing proportion is shrunk to zero by (4.3), the corresponding parameters in this component are set to zero and fewer components are kept for the remaining EM iterations. Here we use the same notation $K$ for the whole process. In practice, during the iterations, $K$ becomes smaller and smaller until the algorithm converges.

Remark 4**.**

Although in theory we require $\epsilon=o(n^{-1/2}/\log{n})$, we can update $\pi$ using (4.3) without choosing $\epsilon$ in practice.

4.2 Turning Parameter Selection and Classification Rule

In terms of selecting the tuning parameter $\lambda$, we follow the suggestion in Wang, Li, and Tsai (2007) and use a BIC-type criterion:

[TABLE]

where $\widehat{K}$ and $\widehat{\beta}$ are estimators of $K_{0}$ and $\beta_{0}$ by maximizing (2.4) for a given $\lambda$.

Let $\widehat{K}$, $\widehat{\pi}$, $\widehat{\beta}$, and $\widehat{\phi}$ be the final estimators of the number of components, the mixture proportions and unknown parameters, respectively. Then, in the sense of clustering, a subject can be assigned to the class whose empirical posterior is the largest. For example, a subject $(Y^{*},X^{*})$ with $m$ times observations is assigned to the class

[TABLE]

Consequently, a nature predictor of $Y^{*}$ is given by $g(X^{T}\widehat{\beta}_{k^{*}})$.

Remark 5**.**

One may claim that $\widehat{\beta}$ would loss some efficiency if the within-subject correlation is strong. It would be better to incorporate correlation information to gain estimation efficiency. However, a correlation analysis would lead to additional computational cost and increase the chance of the convergence problem for the proposed modified EM algorithm. In practice, we suggest to estimate $\beta$ once again given the component information derived from (2.4). Specifically, we first fit the mixture regression model (2.1)-(2.2) and cluster samples into $\widehat{K}$ classes by (4.5); then, in each class, the marginal generalized linear model is estimated by applying GEE with a working correlation structure. It is expected that this two-step technique may improve the estimation efficiency if the correlation of the longitudinal data is strong and the working structure is correctly specified.

5 Simulation studies

In this section, we conduct a set of Monte Carlo simulation studies to assess the finite sample performance of the proposed method. The maximum initial number of clusters is set to be ten, the initial value for the modified EM algorithm is estimated by K-means clustering and the tuning parameter $\lambda$ is selected by the proposed BIC criterion (4.4). To test the classification accuracy and estimation accuracy, we conduct 1000 replications and compare the method (2.4) with the two-step method mentioned in Remark 5 and the QIFC method proposed by Wang and Qu (2014). QIFC is a supervised classification technique for longitudinal data. To permit comparison, we assume that the true number of components, the true class label and the true within-subject correlation are known for the QIFC method. We denote the proposed method and the two-step method as PQL and PQL2 in the following, respectively.

Example 1. Motivated by the real data application, we simulate PBC data from a two-component normal mixture as follows. We set $n=300$, $K=2$, $m_{i}=6$, and $\pi_{1}=\pi_{2}=0.5$. For $k$th component, the mean structure of each response is set as

[TABLE]

and the marginal variance is assumed as $\sigma^{2}_{k}$. The true values of the regression parameters $\beta_{kj}$’s and the marginal variances $\sigma^{2}_{k}$’s are given in Table 2. Covariates $X_{i1}$ are generated independently from Bernoulli distribution $B(1,0.5)$ with 0 for placebo and 1 for D-penicillamine. Covariates $X_{i2}$, representing the age of $i$th patient at entry in years, are generated independently from uniform distribution $U(30,80)$. Covariates $X_{i3}$ are randomly sampled from Bernoulli distribution $B(1,0.5)$ with 0 for male and 1 for female. For each subject, $m_{i}=6$ visit times $Z_{ij}$’s are generated, with the first time being equal to 0 and the remaining five visit times being generated from uniform distributions on intervals (350, 390), (710, 770), (1080, 1160), (1450, 1550), and (1820, 1930) days, respectively. Then, let $X_{ij4}=Z_{ij}/30.5$ be the $j$th visit time of $i$th subject in months. Further, for each subject, we assume the within correlation structure is AR(1) with correlation coefficient $0.6$.

To measure the performance of the proposed tuning parameter selector (4.4), we show the histograms of the estimated component numbers and report the percentage of selecting the correct number of components. To check the convergence of the proposed modified EM algorithm, we draw the evolution of the penalized quasi-likelihood (2.4) in one run. With respect to classification, we generate 100 new subjects from each component with the same setting as in each configuration and measure the performance in terms of the misclassification error rate. We summarize the median and the 95$\%$ confidence interval of misclassification error rate from a model with correctly identified $K_{0}$ for PQL and PQL2 and report these quantities for QIFC as well. To measure the performance of the proposed estimators, the mean values of the estimators, the means of the biases, and the mean squared errors (MSE) for the mixture proportions and regression parameters are reported when the number of components $K_{0}$ is correctly identified. Correspondingly, the mean values, the means of biases, and the MSE of the estimated QIFC estimators are also summarized as a benchmark for comparison. Note that the label switching might arise in practice. Yao (2015) and Zhu and Fan (2016) proposed many feasible labeling methods and algorithms. In our simulation studies, we solve the label switching by putting an order constraint on components’ mean parameters.

Figure 1(a) draws the histogram of the estimated component numbers. It shows that the proposed PQL method with the BIC tuning parameter selector can identify the correct number of components at least with probability 0.962, which is in accordance with the model selection consistency in Theorem 2. Figure 1(b) depicts the evolution of the penalized quasi-likelihood function (2.4) for the simulated data set in one run, showing that how our proposed modified EM algorithm converges numerically.

When the number of components is correctly identified, Table 1(a) reports the median and the 95$\%$ confidence interval of the misclassification error rate from the model-based clustering. We can see that the proposed methods perform better than QIFC with relatively smaller misclassification error rate. Since QIFC is proved asymptotically optimal in terms of misclassification error rate (see Theorem 1 in Wang and Qu, 2014), the observations in Table 1(a) imply the optimality of the proposed methods in terms of misclassification error rate numerically. Further, in terms of parameter estimation, we summarize the estimation of mixture proportions and regression parameters in Table 2. The means of the PQL estimators seem to provide consistent estimates of the regression parameters. It is not surprising that, for regression parameters, the PQL approach performs not as well as the QIFC method with larger bias (in absolute value) and MSE, since QIFC estimators are oracle by assuming the known class memberships and the true within-subject correlation structure. It implies that ignoring the working correlation would affect the efficiency of parameter estimation. However, we can improve the estimating efficiency if the correct correlation information is incorporated. This is reflected in the PQL2 estimators that have much smaller biases (in absolute value) and MSEs compared with the PQL estimators. Indeed, the PQL2 method performs similarly to the QIFC approach.

In addition, combining Table 1(a) and Table 2, we can observe that the two-step technique is able to improve the estimation efficiency for the mean regression parameters without reducing the classification accuracy, which validate our guess in Remark 5 numerically. In general, when the within-subject correlation is strong, it is recommended to use PQL2 to provide more predictive power by utilizing the within-subject correlation information.

Example 2. By design, the application of the proposed method is not restricted to continuous responses, and we next evaluate the performance of PQL and PQL2 on count responses. We generate correlated count outcomes from a two-component overdispersed Poisson mixture with mixture proportions $\pi_{1}=1/3$ and $\pi_{2}=2/3$. For component 1, the mean function of repeated measurements $Y_{ij}$ is

[TABLE]

and the marginal variance is $\phi_{1}\mu_{ij1}=\textup{var}(Y_{ij})=2\mu_{ij1}$. The correlation structure within a subject is AR(1) with correlation coefficient $\rho$. For component 2, $Y_{i}$ has the same correlation matrix as in component 1, except that

[TABLE]

with dispersion parameter $\phi_{2}=1$. The number of repeated measurements $m_{i}$ is randomly drawn from a Poisson distribution with mean 3 and increased by 2, and the sample size is $n=150$. Covariates $X_{ijp}$ are generated independently from uniform distribution $U(0,1)$. Two values of $\rho$ are considered, $\rho=0.3$ and 0.6, to represent different correlation magnitude.

Figure 1(c) and (d) depict the histograms of the estimated component numbers with different correlation magnitude. It shows that our proposed PQL method can identify the correct model in more than 95$\%$ cases. Even with large within-subject correlation, Figure 3 (b) in Appendix B shows that the modified EM algorithm converges numerically with the maximum number of components as 10. Once the model is correctly selected, the classification accuracy is quite satisfactory. Table 1(b) implies that PQL and PQL2 provide more predictive power, especially for large within-subject correlation. With respect to bias and MSE in the estimation of parameters, Table 5 in Appendix B indicates that our modified EM algorithm gives consistent estimates for parameters and mixture proportions by considering the within-class dispersions. Similar to that in Example 1, when the within-subject correlation is large, the PQL2 approach enhances the estimation efficiency by incorporating the correlations within each subject while retaining the class membership prediction accuracy.

Example 3. In the third example, we consider a five-component Gaussian mixture of AR(1), exchangeable (CS), and independence (IND). This is a more challenging example with more components but having different correlation structures. Specifically, we generate 500 samples with mixture proportions $\pi_{1}=\pi_{2}=0.25$, $\pi_{3}=\pi_{4}=0.15$ and $\pi_{5}=0.2$. Conditional on the class label $u_{i}$, the response vector $Y_{i}$ is generated from five multivariate normal distributions:

[TABLE]

where the within-subject correlation structures are set as $R_{i}^{(1)}=R_{iAR(1)}(0.6)$, $R_{i}^{(2)}=R_{iAR(1)}(0.6)$, $R_{i}^{(3)}=R_{iCS}(0.3)$, $R_{i}^{(4)}=R_{iCS}(0.3)$, $R_{i}^{(5)}=R_{iIND}$, and the true values of the regression parameters $(\beta_{k0},\beta_{k})$’s and the variance parameters $\sigma^{2}_{k}$’s are given in Table 6. The number of repeated measurements $m_{i}$ and the covariates are generated as in Example 2.

Figure 1(e) draws the histogram of estimated numbers of components and Figure 4 depicts the evolution of the penalized quasi-likelihood function (2.4) in one run. Though PQL uses a single correlation structure (IND), it is able to identify the correct number of components with high probability, and the corresponding modified EM algorithm converges numerically. Further, the classification results summarized in Table 1(c) shows that PQL gives more accurate prediction of the class’s membership compared with QIFC, which is oracle by assuming the known class memberships and the true different within-subject correlation structures. Table 6 in Appendix B indicates, across different finite mixture correlation models, PQL estimators are still consistent. It may loss some efficiency, but can be improved by PQL2.

6 Application to primary biliary cirrhosis data

In this section, we apply the proposed method to study a doubled-blinded randomized trail in primary biliary cirrhosis (PBC) conducted by the Mayo Clinic between 1974 and 1984 (Dickson, Grambsch, Fleming, Fisher, and Langworthy, 1989).

This data set consists of 312 patients who consented to participate in the randomized placebo-controlled trial with D-penicillamine for treating primary biliary cirrhosis until April 1988. Each patient was supposed to have measurements taken at 6 months, 1 year, and annually thereafter. However, 125 of the original 312 patients had died at updating of follow-up in July 1986. Of the remainder, a sizable portion of patients missed their measurements because of worsening medical condition of some labs, which resulted in an unbalanced data structure. A number of variables were recorded for each patient including ID number, time variables such as age and number of months between enrollment and this visit date, categorical variables such as drug, gender and status, continuous measurement variables such as the serum bilirubin level. PBC is a rare but fatal chronic cholestatic liver disease, with a prevalence of about 50-cases-per-millon-population. Affected patients are typically middle-aged women. As in this data set, the sex ratio is 7.2 : 1 (women to men), where the median age of women patients is 49 years old. Identification of PBC is crucial to balancing the need for medical treatment to halt disease progression and extend survival without need for liver transplantation, while minimizing drug-induced toxicities. Biomedical research indicates that serum bilirubin concentration is a primary indicator to help evaluate and track the absence of liver diseases. It is generally normal at diagnosis (0.1$\sim$1 mg/dl) but rise with histological disease progression (Talwalkar and Lindor, 2003). Therefore, we concentrate on modeling the relationship between marker serum bilirubin and other covariates of interest.

We set the log-transformed serum bilirubin level (lbili) as the response variable, since the original level has positive observed values (Murtaugh, Dickson, van Dam, Malinchoc, Grambsch, Langworthy, and Gips, 1994). Figure 2(a) depicts the plot of a set of observed transformed longitudinal profiles of serum bilirubin marker. It shows that the trend of profiles vary over time and the variability may be large for different patients. The median age of 312 patients is 50 years, but varies between 26 and 79 years. The two sample t-test indicates that there exists significant difference in means of age between male and female groups (p-value = 0.001). Therefore, we consider the marginal semiparametric mixture regression model (2.1)-(2.2) with the identity link. The mean structure in the $k$th component takes the form

[TABLE]

and the marginal variance is assumed as $\textup{var}(Y_{ij})=\sigma^{2}_{k}$, $i=1,\ldots,312$, $j=1,\ldots,m_{i}$, $k=1,\ldots,K$, where variable Trt is a binary variable with 0 for placebo and 1 for D-penicillamine, variable Sex is binary with 0 for male and 1 for female, and variable Time is the number of months between enrollment and this visit date.

We first standardize data so that there is no intercept term in model (6.1). Then, we apply the proposed method to simultaneously select the number of components and to estimate the mixture proportions and unknown parameters. As in the simulation studies, the maximum initial number of clusters is set to be ten and the initial value for the modified EM algorithm is estimated by K-means clustering. For comparison purposes, the standard linear mixed-effects model (LMM) with heterogeneity (Verbeke and Lesaffre, 1996; De la Cruz-Mes$\acute{i}$a, Quintana, and Marshall, 2008) is also considered for continuous response variable lbili. The R package mixAK (Komárek and Komárková, 2013) is used to estimate the model and select the number of groups. The proposed method detects 2 groups, which is same as the clinical classification, while LMM favors 3 groups. Figure 5 in Appendix B depicts the boxplots of residuals in these three groups. The boxplots exhibit the heavy-tailed phenomenon for residuals, especially for those patients in Group 1. It implies that the normality assumptions for the random effects and errors appear inappropriate for modeling this data set. A misspecified distribution of random quantities in the model can seriously influence parameter estimates as well as their standard errors, subsequently leading to invalid statistical inferences. Therefore, it is better to use the proposed semiparametric mixture regression model that only requiring the first two moment conditions of the model distribution. To check the stability of the proposed method, we run our method 100 replications. To be specific, the variable “status” is a triple variable with 0 for censored, 1 for liver transplanted and 2 for dead. It describes the status of a patient at the endpoint of the cohort study. For each run, we randomly draw 80$\%$ of patients for each of these three status without replacement. Figure 1(f) shows that our proposed method selects two groups with high probability.

The resulting estimators of parameters and mixture proportions along with the standard deviations are shown in Table 3. One scientific question of this cohort study is that whether the drug D-pencillamine has effective impact on slowing the rate of increase in serum bilirubin level. According to the estimates and standard deviations with respect to covariate “Trt” in Table 3, it implies that there is little benefit of D-pencillamine to lowering the rate of increase in serum bilirubin level even harmful effect, which is in accordance with findings in other literatures (eg Hoofnagle, David, Schafer, Peters, Avigan, Pappas, Hanson, Minuk, Dusheiko, and Campbell, 1986; Pontecorvo, Levinson, and Roth, 1992).

Another goal of this study is to identify groups of patients with similar characteristics by using the values of the marker serum bilirubin and to see how the bilirubin levels evolve over time. Figure 6 in Appendix B depicts the fitted mean profiles in identified two groups, showing the increasing trend of bilirubin levels in both groups. According to the estimates and standard deviations of parameters in Table 3, it implies that the covariate “Time” is significant and bilirubin level increases over time in both treatment and control arms. Moreover, note that $\widehat{\beta}_{14}=0.068<\widehat{\beta}_{24}=0.313$, which implies that the bilirubin level increases more slowly over time in Group 0. Therefore, from the clinical point of view, Group 0 should correspond to patients with a better prognosis compared to Group 1. To confirm this conclusion, Kaplan-Meier estimates of the survival probabilities are calculated based on data from patients classified in each group. We can see that, from Figure 2(b), the survival prognosis of Group 0 is indeed much better than that of Group 1 with the estimated 5-year survival probability in Group 0 of 0.926 compared to 0.729 in Group 1, and the 10-year survival probabilities 0.771 and 0.310 in Groups 0 and 1, respectively. The p-value of the log rank test is near 0, which implies that the survival distributions corresponding to identified groups are quite different. Further, according to the variable “status”, the group levels for 312 patients are predefined. At the endpoint of the cohort study, 140 of the patients had died, Group 1, while 172 were known to be alive, Group 0. Therefore, it is of interest to compare the classification results using the fitted semiparametric two-component mixture models shown in Table 4. For comparison purposes, the fitting results and classification results of the QIFC method are presented in Tables 3 and  4, respectively. It can be observed that the proposed method provides more accurate classification performance than the QIFC.

7 Conclusion

In this paper, we have proposed a penalized method for learning mixture regression models from longitudinal data which is able to select the number of components in an unsupervised way. The proposed method only requires the first two moment conditions of the model distribution, and thus is suitable for both the continuous and discrete responses. It penalizes the logarithm of mixing proportions, which allows one to simultaneously select the number of components and to estimate the mixture proportions and unknown parameters. Theoretically, we have shown that our proposed approach can select the number of components consistently for general marginal semiparametric mixture regression models. And given the number of components, the estimators of mixture proportions and regression parameters are root-$n$ consistent and asymptotically normal.

To improve the classification accuracy, a modified EM algorithm has been proposed by considering the within-component dispersion. Simulation results and the real data analysis have shown its convergence, but further theoretical investigation is needed. And we have introduced a BIC-type method to select the tuning parameter automatically. Numerical studies show it works well, while the theoretical consistency deserves a further study.

Another issue is the consideration of the within-subject correlation. The proposed penalized approach is introduced under the working independence correlation. Simulation results have implied that it may lose some estimation efficiency, especially when the within-subject correlation is large. Therefore, we suggest a two-step technique to refine the estimates. Simulations show that the efficiency improvement is significant if the correlation information is incorporated and the working structure is correctly specified. It would be worthwhile to systematically study the unsupervised learning of mixtures by incorporating correlations.

Finally, in the presence of missing data at some time points, our implicit assumption is missing completely at random, under which the quasi-likelihood method yield consistent estimates (Liang and Zeger, 1986). Such an assumption is applicable to our motivating example, as patients missed their measurements due to administrative reasons. However, when the missing values are informative, the proposed method has to be modified so as to incorporate missing mechanisms. This is beyond the current scope of the work and would warrant further investigations.

Appendix

A. Proofs of Theorems

Proof of Theorem 1

Recall $\Psi(\theta;Y_{i}|X_{i})=\sum^{K}_{k=1}\pi_{k}\exp\left\{\sum^{m_{i}}_{j=1}q(g(X^{T}_{ij}\beta_{k});Y_{ij})\right\}$, and $\psi(\theta;Y_{i}\mid X_{i})=\log(\Psi(\theta;Y_{i}\mid X_{i}))$. Under Condition C5, $\theta_{0}$ is a maximizer of $n^{-1}\sum^{n}_{i=1}E\{\psi(\theta;Y_{i}\mid X_{i})-\psi(\theta_{0};Y_{i}\mid X_{i})\}$. Then, $\theta_{0}$ is identifiability unique. Therefore, in the spirit of Theorem 17 in Ferguson (1996) and Theorem 2.1 in Bollerslev and Wooldridge (1992), $\widehat{\theta}$ is weak consistent under Conditions C1-C5. Let $\widehat{\theta}^{*}=\sqrt{n}(\widehat{\theta}-\theta_{0})$. Then, $\widehat{\theta}^{*}$ maximizes

[TABLE]

An application of Taylor expansion yields that

[TABLE]

where $\mathbf{1}$ is a $pK+(K-1)$ dimensional all-ones vector. It can be shown that $B_{n}\stackrel{{\scriptstyle P}}{{\longrightarrow}}-B$. Then, by (7.1) and quadratic approximation lemma, we have $\widehat{\theta}^{*}=B^{-1}D_{n}+o_{p}(\mathbf{1})$. Note that $\textup{var}(D_{n})=A$. And under the regularity conditions, we have $D_{n}\stackrel{{\scriptstyle L}}{{\longrightarrow}}N(0,A)$. Hence,

[TABLE]

In order to establish Theorem 2, we need the following lemma first, which can be derived using similar arguments as the proof of Proposition A.1 of Huang et al. (2016).

For a data pair $(Y,X)$ with $m$ times observations, define

[TABLE]

Let $\mathcal{D}$ be the subset of functions of form

[TABLE]

where $D^{1}_{j}f(\beta_{0k};Y\mid X)$ is the first derivative of $f(\beta_{0k};Y\mid X)$ for the $j$th component of $\beta_{0k}$.

Lemma 7.1**.**

Under conditions C1-C6, $\mathcal{D}$ is a Donsker class.

Proof: Under conditions C1-C6, it is straightforward to show that $\mathcal{D}$ satisfies conditions P0 and P1 in Dacunha-Castelle and Gassiat (1999) as in Keribin (2000). Then, there exists a $\Psi_{0}\nu$-square integrable envelope function $\bar{d}(\cdot)$ such that $|d(\gamma;Y|X)|\leq\bar{d}(Y|X)$. On the other hand, the sequences of coefficients of $d(\gamma;Y|X)$ are bounded under the restrictions imposed on $\gamma$. Hence, similar to the proof of Proposition 3.1 in Dacunha-Castelle and Gassiat (1999), we can show that $\mathcal{D}$ has the Donsker property with the bracketing number $N(\varepsilon)=\varepsilon^{-pK}$.

Proof of Theorem 2

In the spirit of proof of Theorem 3.2 in Huang et al. (2016), we divide our proof into two parts. First, we show that there exists a maximizer $(\eta,\gamma)$ such that $\eta=O_{p}(n^{-1/2})$ when $\lambda=a/\sqrt{n}$. It is sufficient to show that, for a large constant $C$, $Q_{\textup{P}}(\eta,\gamma)<Q_{\textup{P}}(0,\gamma)$ when $\eta=C/\sqrt{n}$. Let $\eta=C/\sqrt{n}$, and note that

[TABLE]

Then,

[TABLE]

For $S_{1}$, an application of Taylor expansion yields

[TABLE]

for $\eta=C/\sqrt{n}$, where $|t_{i}|\leq 1$. By Taylor expansion again for $\Psi(\eta,\gamma;Y\mid X)$ at $\eta=0$, we have $\Psi(\eta,\gamma;Y\mid X)=\Psi_{0}(Y\mid X)+\eta\Psi^{\prime}(0,\gamma;Y\mid X)+\frac{\eta^{2}}{2}\Psi^{\prime\prime}(\tilde{\theta},\gamma;Y\mid X),$ for a $\tilde{\theta}\leq\theta$. Then, by conditions C1-C5, we have

[TABLE]

By Lemma 7.1 for the class $\mathcal{D}$, we know that $\frac{1}{\sqrt{n}}\sum^{n}_{i=1}\frac{\Psi^{\prime}(0,\gamma;Y_{i}\mid X_{i})}{\Psi_{0}(Y_{i}\mid X_{i})}$ converges uniformly in distribution to a Gaussian process and $\sum^{n}_{i=1}\left\{\frac{\Psi^{\prime}(0,\gamma;Y_{i}\mid X_{i})}{\Psi_{0}(Y_{i}\mid X_{i})}\right\}^{2}=O_{p}(n)$ by the law of large numbers. Therefore,

[TABLE]

For $S_{2}$, we know that $|\pi_{l}-\pi_{0(l-K+K_{0})}|=|\eta\rho_{l-K+K_{0}}|\leq\frac{C}{\sqrt{n}},\quad l=K-K_{0}+1,\ldots,K,$ by the restriction condition on $\rho_{k}$, $k=1,\ldots,K_{0}$. Thus, by Taylor expansion, we have

[TABLE]

if $\sqrt{n}\lambda\rightarrow a$. Therefore, when $C$ is large enough, the second term in $S_{1}$ dominates $S_{2}$ and other terms in $S_{1}$. Consequently, we have $Q_{\textup{P}}(\eta,\gamma)-Q_{\textup{P}}(0,\gamma)<0$ with probability tending to one. Hence, there exists a maximizer $(\eta,\gamma)$ such that $\eta=O_{p}(n^{-1/2})$ with probability tending to one.

Then, we show that $\widehat{K}=K_{0}$ or $\widehat{\pi}_{k}=0,k=1,\ldots,K-K_{0}$ when the maximizer $(\eta,\gamma)$ satisfies $\eta=O_{p}(n^{-1/2})$. We first show that, for any maximizer $Q_{\textup{P}}(\eta^{*},\gamma^{*})$ with $|\eta^{*}|\leq Cn^{-1/2}$, if there is a $k\leq K-K_{0}$ such that $Cn^{-1/2}\geq\pi^{*}_{k}>n^{-1/2}/\log n$, there exists another maximizer of $Q_{\textup{P}}(\eta,\gamma)$ in the area of $|\eta|\leq Cn^{-1/2}$. It is equivalent to show that $Q_{\textup{P}}(\eta^{*},\gamma^{*})<Q_{\textup{P}}(0,\gamma^{*})$ holds with probability tending to one for any such kind maximizer $Q_{\textup{P}}(\eta^{*},\gamma^{*})$ with $|\eta^{*}|\leq Cn^{-1/2}$. For any $k<K-K_{0}+1$, we have

[TABLE]

As shown before, we have $S_{1}+S_{2}=O_{p}(C^{2})$. For $S_{3}$, because $\epsilon=o(n^{-1/2}/\log n)$ and $\pi_{k}<n^{-1/2}/\log n$, we have

[TABLE]

which implies that $S_{3}$ dominates $S_{1}$ and $S_{2}$, and hence $Q_{\textup{P}}(\eta^{*},\gamma^{*})<Q_{\textup{P}}(0,\gamma^{*})$. So, in the following step, we only need to consider the maximizer $Q_{\textup{P}}(\widehat{\eta},\widehat{\gamma})$ with $|\widehat{\eta}|\leq Cn^{-1/2}$ and $\widehat{\pi}_{k}<n^{-1/2}/\log n$ for $k\leq K-K_{0}$.

Let $Q^{*}(\theta)=Q_{\textup{P}}(\theta)-\xi(\sum^{K}_{k=1}\pi_{k}-1)$, where $\xi$ is a Lagrange multiplier. Then, it is sufficient to show that, for the maximizer $(\eta,\gamma)$,

[TABLE]

with probability tending to one. For $k=1,\ldots,K$, note that $\widehat{\pi}_{k}$ satisfies

[TABLE]

where $f_{l}(\beta_{l};Y_{i}\mid X_{i})=\exp\left\{\sum^{m_{i}}_{j=1}q(g(X^{T}_{ij}\beta_{l});Y_{ij})\right\}$. By the law of large numbers, the first term of (7.3) is of order $O_{p}(n)$. If $k>K-K_{0}$ and $\eta=O_{p}(n^{-1/2})$, we have that

[TABLE]

Hence, the second term of (7.3) is of order $O_{p}(n\lambda)=o_{p}(n)$. Thus, $\xi=O_{p}(n)$. If $k\leq K-K_{0}$, since $\widehat{\pi}_{k}=O_{p}(n^{-1/2}/\log n)$, $\lambda=a/\sqrt{n}$ and $\epsilon=o(n^{-1/2}/\log n)$, we have

[TABLE]

with probability tending to one. Hence, the second term in (7.3) dominates the first and third terms when $k\leq K-K_{0}$ and $\widehat{\pi}_{k}<n^{-1/2}/\log n$, which implies that (7.2) holds or equivalently, $\widehat{\pi}_{k}=0$, $k=1,\ldots,K-K_{0}$ with probability tending to one. This completes the proof of Theorem 2.

B. Tables and Graphs