Mixture-based Multiple Imputation Model for Clinical Data with a Temporal Dimension

TL;DR

This paper introduces a novel mixture-based multiple imputation model tailored for clinical multivariable time series data, effectively capturing both cross-sectional and temporal correlations to improve missing data imputation accuracy.

Contribution

The work presents a new imputation approach combining Gaussian processes with mixture models and individualized weights, specifically designed for clinical time series with missing values.

Findings

Outperforms existing imputation methods on real-world datasets.

Provides more accurate imputations across various datasets.

Effectively models both cross-sectional and temporal dependencies.

Abstract

The problem of missing values in multivariable time series is a key challenge in many applications such as clinical data mining. Although many imputation methods show their effectiveness in many applications, few of them are designed to accommodate clinical multivariable time series. In this work, we propose a multiple imputation model that capture both cross-sectional information and temporal correlations. We integrate Gaussian processes with mixture models and introduce individualized mixing weights to handle the variance of predictive confidence of Gaussian process models. The proposed model is compared with several state-of-the-art imputation algorithms on both real-world and synthetic datasets. Experiments show that our best model can provide more accurate imputation than the benchmarks on all of our datasets.

Tables

Table 1. TABLE I : Main symbols and definitions

Symbol	Definition
$𝒳$, $𝒯$	Measurement and time tensor
$x_{p,v,-b}$	Measurements in fiber $x_{p,v,:}$ excluding $x_{p,v,b}$
$t_{p,v,-b}$	Times in fiber $t_{p,v,:}$ excluding $t_{p,v,b}$
$Mi{x}_{v,b}$	Mixture model for $v$ and $b$.
$V_{v,b}$	Concatenation of $x_{:,-v,b}$ and $x_{:,v,-b}$
$𝒩$(${𝝁}_{v,b}^{\left(k\right)}$,${𝚺}_{v,b}^{\left(k\right)}$)	The $k$th prior multivariate normal distribution in $Mi{x}_{v,b}$
$m^{\left(3\right)}(\cdot )$, ${\Sigma }^{\left(3\right)}(\cdot )$	Predictive mean and variance of a GP model
${\gamma }_{v,b}$	The set of all trainable parameters of $Mi{x}_{v,b}$
$I$	The identity matrix

Table 2. TABLE II : Overall MASE by dataset and imputation model. The bold numbers are the significantly best values among all imputation models.

Dataset	GP	MTGP1	M-RNN	GMM	MICE	3D-MICE	MixMI-LL	MixMI
Real-world MIMIC ($d$=0)	0.13072	0.12599	0.12512	0.11741	0.11763	0.11186	0.09304	0.09285
Synthetic MIMIC ($d$=0.5)	0.10466	-	0.12320	0.11455	0.11573	0.09087	0.08612	0.08448
Synthetic MIMIC ($d$=1)	0.09220	-	0.12201	0.11436	0.11561	0.07715	0.08427	0.07538
NMEDW	0.18353	0.17370	0.14607	0.13593	0.13718	0.13624	0.11600	0.11589
1 MTGP cannot have multiple inputs, thus not applicable to the synthetic datasets.

Table 3. TABLE III : MASE on the real-world MIMIC dataset by variable and imputation model. The bold numbers are the best values among all imputation models.

Variable	GP	MTGP	M-RNN	GMM	MICE	3D-MICE	MixMI-LL	MixMI
Chloride	0.12993	0.12549	0.10865	0.10650	0.10575	0.10836	0.08664	0.08603
Potassium	0.11533	0.11256	0.11222	0.10963	0.10997	0.10822	0.09453	0.09442
Bicarbonate	0.13196	0.12905	0.12371	0.12231	0.12275	0.11984	0.10302	0.10254
Sodium	0.12525	0.11939	0.11557	0.10159	0.10138	0.10727	0.08787	0.08807
Hematocrit	0.11436	0.10780	0.09189	0.06518	0.06558	0.06726	0.05486	0.05482
Hemoglobin	0.14168	0.12997	0.06556	0.05774	0.05772	0.06301	0.05117	0.05103
MCV	0.14215	0.13732	0.13798	0.13391	0.13474	0.13340	0.11634	0.11657
Platelets	0.13855	0.13087	0.14369	0.14203	0.14236	0.12815	0.10090	0.10070
WBC count	0.13963	0.13614	0.14583	0.14026	0.14068	0.13060	0.10934	0.10913
RDW	0.14592	0.13668	0.15938	0.15778	0.15836	0.13897	0.11340	0.11340
Blood urea nitrogen (BUN)	0.12358	0.12720	0.16345	0.15160	0.15189	0.11814	0.09479	0.09410
Creatinine	0.13341	0.12803	0.14904	0.13211	0.13212	0.12217	0.10067	0.10014
Glucose	0.12491	0.12420	0.11677	0.11771	0.11794	0.11921	0.10493	0.10501

Equations

$$x_{:,v,b}=x_{:,-v,b}{\beta }_{v,b}^{(1)}+{ϵ}_{v,b}^{(1)},{ϵ}_{v,b}^{(1)}\sim \mathcal{N}(0,{{\sigma }_{v,b}^{(1)}}^{2}I)$$

$$x_{:,v,b}|x_{:,-v,b},{\beta }_{v,b}^{(1)},{{\sigma }_{v,b}^{(1)}}^2\sim \mathcal{N}(x_{:,-v,b}{\beta }_{v,b}^{(1)},{{\sigma }_{v,b}^{(1)}}^{2}I).$$

$$x_{:,v,b}=x_{:,v,-b}{\beta }_{v,b}^{(2)}+{ϵ}_{v,b}^{(2)},{ϵ}_{v,b}^{(2)}\sim \mathcal{N}(0,{{\sigma }_{v,b}^{(2)}}^{2}I)$$

$$x_{:,v,b}|x_{:,v,-b},{\beta }_{v,b}^{(2)},{{\sigma }_{v,b}^{(2)}}^2\sim \mathcal{N}(x_{:,v,-b}{\beta }_{v,b}^{(2)},{{\sigma }_{v,b}^{(2)}}^{2}I).$$

$$\begin{array}{rl}{x}_{p,v,b}& ={\mu }_{p,v,b}+f(t_{p,v,b}),\\ f(t_{p,v,b})& \sim \mathcal{G}\mathcal{P}(0,\mathcal{K}(t_{p,v,b},t_{p,v,b^{\prime }}))\end{array}$$

$$\begin{array}{rl}{x}_{p,v,b}|{\alpha }_{p,v,b}& \sim \mathcal{N}(m^{(3)}({\alpha }_{p,v,b}),{\Sigma }^{(3)}({\alpha }_{p,v,b}))\\ {\alpha }_{p,v,b}& =({\theta }_{v,b},x_{p,v,-b},t_{p,v,-b})\end{array}$$

$$\begin{array}{rl}p(& x_{:,v,b},V_{v,b})\\ =& {\pi }_{v,b}^{(1)}\mathcal{N}(V_{v,b}|{\delta }_{v,b}^{(1)})\mathcal{N}(x_{:,v,b}|x_{:,-v,b}{\beta }_{v,b}^{(1)},{{\sigma }_{v,b}^{(1)}}^{2}I)\\ +& {\pi }_{v,b}^{(2)}\mathcal{N}(V_{v,b}|{\delta }_{v,b}^{(2)})\mathcal{N}(x_{:,v,b}|x_{:,v,-b}{\beta }_{v,b}^{(2)},{{\sigma }_{v,b}^{(2)}}^{2}I)\\ +& {\pi }_{v,b}^{(3)}\mathcal{N}(V_{v,b}|{\delta }_{v,b}^{(3)})\mathcal{N}(x_{:,v,b}|m^{(3)}({\mathbf{\alpha }}_{v,b}),\text{diag}({\Sigma }^{(3)}({\mathbf{\alpha }}_{v,b})))\end{array}$$

$$\begin{array}{rl}{p}_1& (x_{:,v,b},V_{v,b})\\ & =p(x_{:,v,b}|V_{v,b})p(V_{v,b})\\ & =\mathcal{N}(x_{:,v,b}|{\mu }_1(V_{v,b}),{\sigma }_1(V_{v,b}))\mathcal{N}(V_{v,b}|{\delta }_{v,b}^{(1)})\\ & =\mathcal{N}(x_{:,v,b}|x_{:,-v,b}{\beta }_{v,b}^{(1)},{{\sigma }_{v,b}^{(1)}}^{2}I)\mathcal{N}(V_{v,b}|{\delta }_{v,b}^{(1)})\end{array}$$

$$p(V_{v,b})=\sum _{k=1}^3{\pi }_{v,b}^{(k)}\mathcal{N}(V_{v,b}|{\delta }_{v,b}^{(k)})$$

$$\begin{array}{rl}p& (x_{:,v,b}|V_{v,b})\\ & =\frac{p(x_{:,v,b},V_{v,b})}{p(V_{v,b})}\\ & =C_{v,b}^{(1)}\mathcal{N}(x_{:,v,b}|x_{:,-v,b}{\beta }_{v,b}^{(1)},{{\sigma }_{v,b}^{(1)}}^{2}I)\\ & +C_{v,b}^{(2)}\mathcal{N}(x_{:,v,b}|x_{:,v,-b}{\beta }_{v,b}^{(2)},{{\sigma }_{v,b}^{(2)}}^{2}I)\\ & +C_{v,b}^{(3)}\mathcal{N}(x_{:,v,b}|m^{(3)}({\mathbf{\alpha }}_{v,b}),\text{diag}({\Sigma }^{(3)}({\mathbf{\alpha }}_{v,b})))\end{array}$$

$${\hat{\gamma }}_{v,b}=\arg \underset{{\gamma }_{v,b}}{\max }\ln p(x_{:,v,b}^{obs},V_{v,b}^{tr};{\gamma }_{v,b})$$

$${\beta }_{v,b}^{(1)},{{\sigma }_{v,b}^{(1)}}^2,{\beta }_{v,b}^{(2)},{{\sigma }_{v,b}^{(2)}}^2,{\theta }_{v,b}\text{and}{\pi }_{v,b}^{(k)},{\mathbf{\mu }}_{v,b}^{(k)},{\mathbf{\Sigma }}_{v,b}^{(k)}\text{for}k=1,2,3.$$

$${\Pi }_{p,v,b}^{(1)}{x}_{p,-v,b}{\hat{\beta }}_{v,b}^{(1)}+{\Pi }_{p,v,b}^{(2)}{x}_{p,v,-b}{\hat{\beta }}_{v,b}^{(2)}+{\Pi }_{p,v,b}^{(3)}{m}^{(3)}({\hat{\alpha }}_{p,v,b})$$

$${\Pi }_{p,v,b}^{(k)}=\frac{{\hat{\pi }}_{v,b}^{(k)}\mathcal{N}(V_{p,v,b}|{\hat{\mathbf{\mu }}}_{v,b}^{(k)},{\hat{\mathbf{\Sigma }}}_{v,b}^{(k)})}{{\sum }_{j=1}^3{\hat{\pi }}_{v,b}^{(j)}\mathcal{N}(V_{p,v,b}|{\hat{\mathbf{\mu }}}_{v,b}^{(j)},{\hat{\mathbf{\Sigma }}}_{v,b}^{(j)})},k=1,2,3$$

$$\begin{array}{rl}\mathcal{L}& ({\gamma }_{v,b})\\ & =\sum _{p\in P_{v,b}^{tr}}\sum _{q_{v,b}}{Q}_p(q_{v,b})\ln \frac{p(x_{p,v,b},V_{p,v,b},q_{v,b};{\gamma }_{v,b})}{Q_p(q_{v,b})}\\ & \le \sum _{p\in P_{v,b}^{tr}}\ln \sum _{k=1}^{3}p(q_{v,b}=k)p(x_{p,v,b},V_{p,v,b};{\gamma }_{v,b}|q_{v,b}=k)\\ & =\ell ({\gamma }_{v,b})\end{array}$$

$$Q_p(q_{v,b}=k)=\frac{p(q_{v,b}=k)p(x_{p,v,b},V_{p,v,b}|q_{v,b}=k)}{{\sum }_{j=1}^{3}p(q_{v,b}=j)p(x_{p,v,b},V_{p,v,b}|q_{v,b}=j)}.$$

$$\begin{array}{rl}\Delta {\tilde{x}}_i& =\frac{|x_i-x_{i-1}|}{{\sum }_{i=2}^B|x_i-x_{i-1}|}\\ \Delta {\tilde{t}}_i& =\frac{|t_i-t_{i-1}|}{{\sum }_{i=2}^B|t_i-t_{i-1}|}.\end{array}$$

$$\begin{array}{rl}{t}_i^{\prime }& =\{\begin{array}{ll}{t}_1,& \text{if }i=1\\ t_i+{\sum }_{j=2}^i[d(\Delta {\tilde{x}}_j-\Delta {\tilde{t}}_j)S],& \text{otherwise}\end{array}\\ S& =\sum _{j=2}^B(|t_j-t_{j-1}|).\end{array}$$

$$\begin{array}{r}MASE(v)=\frac{1}{{\sum }_{\bar{p}}{I}_{\bar{p},v}}\sum _p\frac{{\sum }_{i\in mas{k}_{p,v}}|{\tilde{x}}_{p,v,i}-x_{p,v,i}^{obs}|}{\frac{J_{p,v}}{J_{p,v}-1}{\sum }_{j=2}^{J_{p,v}}|Y_{p,v,j}-Y_{p,v,j-1}|}.\end{array}$$

$$\begin{array}{rl}[& C_{p,v,b}^{(1)}{]}^{(j)}\\ & =[{\pi }_{v,b}^{(1)}{]}^{(j)}[D_{p,v,b}^{(1)}{]}^{(j)}\mathcal{N}(x_{p,v,b}|x_{p,-v,b}[{\beta }_{v,b}^{(1)}{]}^{(j)},[{{\sigma }_{v,b}^{(1)}}^2{]}^{(j)})\\ [& C_{p,v,b}^{(2)}{]}^{(j)}\\ & =[{\pi }_{v,b}^{(2)}{]}^{(j)}[D_{p,v,b}^{(2)}{]}^{(j)}\mathcal{N}(x_{p,v,b}|x_{p,v,-b}[{\beta }_{v,b}^{(2)}{]}^{(j)},[{{\sigma }_{v,b}^{(2)}}^2{]}^{(j)})\\ [& C_{p,v,b}^{(3)}{]}^{(j)}\\ & =[{\pi }_{v,b}^{(3)}{]}^{(j)}[D_{p,v,b}^{(3)}{]}^{(j)}\mathcal{N}(x_{p,v,b}|m^G([{\alpha }_{p,v,b}{]}^{(j)}),{\Sigma }^G([{\alpha }_{p,v,b}{]}^{(j)}))\\ [& D_{p,v,b}^{(k)}{]}^{(j)}=\mathcal{N}(V_{p,v,b}|[{\mathbf{\mu }}_{v,b}^{(k)}{]}^{(j)},[{\mathbf{\Sigma }}_{v,b}^{(k)}{]}^{(j)}),k=1,2,3\\ [& w_{p,v,b}^{(k)}{]}^{(j)}=\frac{[C_{p,v,b}^{(k)}{]}^{(j)}}{{\sum }_{i=1}^3[C_{p,v,b}^{(i)}{]}^{(j)}},k=1,2,3\end{array}$$

$$\begin{array}{rl}[{\pi }_{v,b}^{(k)}{]}^{(j+1)}& =\frac{1}{|P_{v,b}^{tr}|}\sum _{p\in P_{v,b}^{tr}}[w_{p,v,b}^{(k)}{]}^{(j)},k=1,2,3\\ [{\mathbf{\mu }}_{v,b}^{(k)}{]}^{(j+1)}& =\frac{{\sum }_{p\in P_{v,b}^{tr}}[w_{p,v,b}^{(k)}{]}^{(j)}{V}_{p,v,b}}{{\sum }_{p\in P_{v,b}^{tr}}[w_{p,v,b}^{(k)}{]}^{(j)}}\\ [{\mathbf{\Sigma }}_{v,b}^{(k)}{]}^{(j+1)}& =\frac{1}{{\sum }_{p\in P_{v,b}^{tr}}[w_{p,v,b}^{(k)}{]}^{(j)}}\sum _{p\in P_{v,b}^{tr}}[w_{p,v,b}^{(k)}{]}^{(j)}[U_{p,v,b}^{(k)}{]}^{(j+1)}\\ [U_{p,v,b}^{(k)}{]}^{(j+1)}& =\{V_{p,v,b}-[{\mathbf{\mu }}_{v,b}^{(k)}{]}^{(j+1)}\}\{V_{p,v,b}-[{\mathbf{\mu }}_{v,b}^{(k)}{]}^{(j+1)}{\}}^{\prime }\\ [{\beta }_{v,b}^{(1)}{]}^{(j+1)}& =\{\{Z_{v,b}^{\prime }[{\mathbf{w}}_{v,b}^{(1)}{]}^{(j)}{Z}_{v,b}{\}}^{-1}{Z}_{v,b}^{\prime }[{\mathbf{w}}_{v,b}^{(1)}{]}^{(j)}{x}_{:,v,b}^{obs}{\}}^{\prime }\\ [{{\sigma }_{v,b}^{(1)}}^2{]}^{(j+1)}& =\frac{1}{{\sum }_{p\in P_{v,b}^{tr}}[w_{p,v,b}^{(1)}{]}^{(j)}}[S_{v,b}^{(1)}{]}^{(j+1)}\\ [S_{v,b}^{(1)}{]}^{(j+1)}& =\sum _{p\in P_{v,b}^{tr}}[w_{p,v,b}^{(1)}{]}^{(j)}\{x_{p,v,b}-x_{p,-v,b}[{\beta }_{v,b}^{(1)}{]}^{(j+1)}{\}}^2\\ [{\beta }_{v,b}^{(2)}{]}^{(j+1)}& =\{\{Y_{v,b}^{\prime }[{\mathbf{w}}_{v,b}^{(2)}{]}^{(j)}{Y}_{v,b}{\}}^{-1}{Y}_{v,b}^{\prime }[{\mathbf{w}}_{v,b}^{(2)}{]}^{(j)}{x}_{:,v,b}^{obs}{\}}^{\prime }\\ [{{\sigma }_{v,b}^{(2)}}^2{]}^{(j+1)}& =\frac{1}{{\sum }_{p\in P_{v,b}^{tr}}[w_{p,v,b}^{(2)}{]}^{(j)}}[S_{v,b}^{(2)}{]}^{(j+1)}\\ [S_{v,b}^{(2)}{]}^{(j+1)}& =\sum _{p\in P_{v,b}^{tr}}[w_{p,v,b}^{(2)}{]}^{(j)}\{x_{p,v,b}-x_{p,v,-b}[{\beta }_{v,b}^{(2)}{]}^{(j+1)}{\}}^2\\ [{\theta }_{v,b}{]}^{(j+1)}& =G([{\mathbf{w}}_{v,b}^{(3)}{]}^{(j)},[{\theta }_{v,b}{]}^{(j)},x_{:,v,b}^{obs},Y_{v,b},t_{:,v,:})\end{array}$$

$$\begin{array}{rl}{x}_{p,t}& ={\mu }_{p,t}+f(t),\\ f(t)& \sim \mathcal{G}\mathcal{P}(0,\mathcal{K}(t,t^{\prime }))\end{array}$$

$$\begin{array}{rl}{m}^{(3)}(\theta ,x_{p,-t},\bar{t})& =(\frac{1-r^T{R}^{-1}{1}_n}{1_n^T{R}^{-1}{1}_n}{1}_n^T+r^T)R^{-1}{x}_{p,-t}\\ {\Sigma }^{(3)}(\theta ,x_{p,-t},\bar{t})& ={\sigma }_f^2[1-r^T{R}^{-1}r+\frac{(1-1_n^T{R}^{-1}r{)}^2}{1_n{R}^{-1}{1}_n}]\end{array}$$

$$R_{t,t^{\prime }}=\exp (-\theta |t-t^{\prime }{|}^2).$$

$${\sigma }_f^2=\frac{C^T{R}^{-1}C}{n},C=x_{p,-t}-1_n(1_n^T{R}^{-1}{1}_n{)}^{-1}(1_n^T{R}^{-1}{x}_{p,-t})$$

$$\frac{\partial L}{\partial \theta }=\sum _{p=1}^{|p^{tr}|}{w}_p\frac{\partial }{\partial \theta }\ln \mathcal{N}(x_{p,t};m^{(3)}(\theta ,x_{p,-t},\bar{t}),{\Sigma }^{(3)}(\theta ,x_{p,-t},\bar{t})).$$

$$\begin{array}{rl}\frac{\partial L}{\partial \theta }& =\sum _{p=1}^{|p^{tr}|}{w}_p\frac{\partial }{\partial \theta }\ln \mathcal{N}(x_{p,t};g_p(\theta ),h_p(\theta ))\\ & =\sum _{p=1}^{|p^{tr}|}{w}_p\frac{\partial }{\partial \theta }\{\ln \frac{1}{\sqrt[]{2\pi h_p(\theta )}}-\frac{[x_{p,t}-g_p(\theta ){]}^2}{2h_p(\theta )}\}\\ & =\sum _{p=1}^{|p^{tr}|}{w}_p\{-\frac{1}{2h_p(\theta )}\frac{\partial h_p(\theta )}{\partial \theta }-\frac{\partial }{\partial \theta }\frac{[x_{p,t}-g_p(\theta ){]}^2}{2h_p(\theta )}\}\\ & =\sum _{p=1}^{|p^{tr}|}{w}_p\{-\frac{1}{2h_p(\theta )}\frac{\partial h_p(\theta )}{\partial \theta }\\ & -\frac{1}{2h_p^2(\theta )}\{2[x_{p,t}-g_p(\theta )][-\frac{\partial g_p(\theta )}{\partial \theta }]h_p(\theta )\\ & -\frac{\partial h_p(\theta )}{\partial \theta }[x_{p,t}-g_p(\theta ){]}^2\}\}\\ & =\sum _{p=1}^{|p^{tr}|}{w}_p\{-\frac{1}{2h_p(\theta )}\frac{\partial h_p(\theta )}{\partial \theta }\\ & +\frac{[x_{p,t}-g_p(\theta )]\frac{\partial g_p(\theta )}{\partial \theta }}{h_p(\theta )}+\frac{\frac{\partial h_p(\theta )}{\partial \theta }[x_{p,t}-g_p(\theta ){]}^2}{2h_p^2(\theta )}\}.\end{array}$$

$$\begin{array}{rl}\frac{\partial g_p(\theta )}{\partial \theta }& =(\frac{\partial H_1}{\partial \theta }{R}^{-1}+H_1\frac{\partial R^{-1}}{\partial \theta })x_{p,-t}\\ \frac{\partial h_p(\theta )}{\partial \theta }& ={\sigma }_f^2\frac{\partial H_3}{\partial \theta }+\frac{\partial {\sigma }_f^2}{\partial \theta }{H}_3\end{array}$$

$$\begin{array}{rl}{H}_1& =\frac{[1-(r{R}^{-1}{1}_n)]}{1_n^T{R}^{-1}{1}_n}{1}_n^T+r\\ \frac{\partial H_1}{\partial \theta }& =\frac{-(\frac{\partial r}{\partial \theta }{R}^{-1}+r\frac{\partial R^{-1}}{\partial \theta })1_n(1_n^T{R}^{-1}{1}_n)}{1_n^T{R}^{-1}{1}_n^2}\\ & -\frac{(1_n^T\frac{\partial R^{-1}}{\partial \theta }{1}_n)[1-(r{R}^{-1}{1}_n)]}{1_n^T{R}^{-1}{1}_n^2}{1}_n^T+\frac{\partial r}{\partial \theta }\\ fc& =(1-1_n^T{R}^{-1}{r}^T{)}^2\\ gc& =1_n^T{R}^{-1}{1}_n\\ \frac{\partial fc}{\partial \theta }& =2(1-1_n^T{R}^{-1}{r}^T)[-1_n^T(\frac{\partial R^{-1}}{\partial \theta }{r}^T+R^{-1}\frac{\partial r^T}{\partial \theta })]\\ \frac{\partial gc}{\partial \theta }& =1_n^T\frac{\partial R^{-1}}{\partial \theta }{1}_n\\ H_2& =\frac{(1-1_n^T{R}^{-1}{r}^T{)}^2}{1_n^T{R}^{-1}{1}_n}\\ \frac{\partial H_2}{\partial \theta }& =\frac{\frac{\partial fc}{\partial \theta }gc-\frac{\partial gc}{\partial \theta }fc}{g{c}^2}\end{array}$$

$$\begin{array}{rl}{H}_3& =1-(r{R}^{-1}{r}^T)+H_2\\ \frac{\partial H_3}{\partial \theta }& =-(\frac{\partial r}{\partial \theta }{R}^{-1}{r}^T+r\frac{\partial R^{-1}}{\partial \theta }{r}^T+r{R}^{-1}\frac{\partial r^T}{\partial \theta })+\frac{\partial H_2}{\partial \theta }\\ H_4& =x_{p,-t}-1_n\frac{(1_n^T{R}^{-1}{x}_{p,-t})}{1_n^T{R}^{-1}{1}_n}\\ \frac{\partial H_4}{\partial \theta }& =-1_n\frac{1}{(1_n^T{R}^{-1}{1}_n{)}^2}[(1_n^T\frac{\partial R^{-1}}{\partial \theta }{x}_{p,-t})(1_n^T{R}^{-1}{1}_n)\\ & -(1_n^T\frac{\partial R^{-1}}{\partial \theta }{1}_n)(1_n^T{R}^{-1}{x}_{p,-t})]\\ \frac{\partial {\sigma }_f^2}{\partial \theta }& =\frac{1}{n}[(\frac{\partial H_4}{\partial \theta }{)}^T{R}^{-1}{H}_4+H_4^T\frac{\partial R^{-1}}{\partial \theta }{H}_4+H_4^T{R}^{-1}\frac{\partial H_4}{\partial \theta }].\end{array}$$

Full text

©2019 IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works. DOI: 10.1109/BigData47090.2019.9005672

I Introduction

The computational modeling in clinical applications attracts growing interest with the realization that the quantitative understanding of patient pathophysiological progression is crucial to clinical studies [1]. With a comprehensive and precise modeling, we can have a better understanding of a patient’s state, offer more precise diagnosis and provide better individualized therapies [2]. Researchers are increasingly motivated to build more accurate computational models from multiple types of clinical data. However, missing values in clinical data challenge researchers using analytic techniques for modeling, as many of the techniques are designed for complete data.

Traditional strategies used in clinical studies to handle missing values include deleting records with missing values and imputing missing entries by mean values. However, deleting records with missing values and some other filtering strategies can introduce biases [3] that can impact modeling in many ways, thus limiting its generalizability. Mean imputation is widely used by researchers to handle missing values. However, it is shown to yield less effective estimates than many other modern imputation techniques [4, 5, 6, 7], such as maximum likelihood approaches and multiple imputation methods (e.g. multiple imputation by chained equations (MICE) [8]). The maximum likelihood and multiple imputation methods are based on solid statistical foundations and become standard in the last few decades [9, 10].

Multiple imputation is originally proposed by Rubin [11] and with the idea of replacing each missing value with multiple estimates that reflect variation within each imputation model. The classic multiple imputation method treats the complete variables as predictors and incomplete variables as outcomes, and uses multivariate regression models to estimate missing values [12, 13]. Another imputation approach relies on chained equations where a sequence of regression models is fitted by treating one incomplete variable as the outcome and other variables with possible best estimates for missing values as predictors [14, 15, 16]. Multiple imputation has its success in healthcare applications, where many imputation methods designed for various tasks are based on multiple imputation [17, 18, 19, 20, 21, 22, 23].

In recent years, additional imputation methods are proposed. Although many imputation methods [7, 8, 24, 14, 25, 26, 27, 28] show their effectiveness in many applications, few of them are designed for time series-based clinical data. These clinical data are usually multivariable time series, where patients have measurements of multiple laboratory tests at different times. Many methods are designed for cross-sectional imputation (measurements taken at the same time point) and do not consider temporal information that is useful in making predictions or imputing missing values. Ignoring informative temporal correlations and only capturing cross-sectional information may yield less effective imputation.

In order to address the limitations mentioned above, we present a mixture-based multiple imputation model (MixMI) for clinical time series 111 The implementations of our model is available at: https://github.com/y-xue/MixMI.. Our model captures both cross-sectional information and temporal correlations to estimate missing values using mixture models. We model the distribution of measurements using a mixture model. The mixture is composed of linear regression to model cross-sectional correlations and Gaussian processes (GPs) to capture temporal correlations. The problem of integrating GP within a standard mixture model is that GP models in all patient cases get the same mixing weights, while the confidence of predictions by GP models can vary largely across different patient cases. We overcome this problem by introducing individualized mixing weights for each patient cases, instead of assigning a fixed weight. We train our model using the Expectation-Maximization (EM) algorithm. We demonstrate the effectiveness of our model by comparing it with several state-of-the-art imputation algorithms on multiple clinical datasets.

Our main contributions are summarized as follows.

1. To the best of our knowledge, we are the first to build imputation model for time series by integrating GP within mixture models. We address the problem that all GP models in all patient cases get a fixed mixing weight by introducing individualized mixing weights.

2. We test the performance of our model on two real-world clinical datasets and several synthetic datasets. Our best model outperforms all comparison models including several state-of-the-art imputation models. Using synthetic datasets, we also explore and discover the properties of the data that benefit our model and/or comparison models. Experiments show that our best model is robust to the variation of these properties and outperforms comparison models on all synthetic datasets.

The remainder of this paper is structured as follows. Section II discusses related work while in Section III, the proposed method is described. The experimental setup, including dataset collection and evaluation procedure, is described in Section IV. Section V discusses the computational results and underlying analyses. The conclusions are drawn in Section VI.

II Related work

Research in designing imputation methods for multivariable time series attracts growing interest in recent decades. Previous studies generally fall into two categories. One comes from methods using linear or other simple parametric functions to estimate missing values. The other is the methods treating time series as smooth curves and estimating missing values using GP or other nonparametric methods.

In the first category, multivariable time series are modeled based on either linear models [29], linear mixed models [30, 31] or autoregressive models [32, 33]. However, in these methods, the potential trajectories of variables are only limited to linear or other simple parametric functions. Alternatively, many authors choose GPs or other nonparametric functions to model time series. Compared to linear models, GPs only have locality constraints in which close time points in a time series usually have close values. Therefore, GPs bring in more flexibility in capturing temporal trajectories of variables.

However, directly applying GPs on imputation for multivariable time series usually yields less effective imputation due to the lack of considerations for the similarities and correlations across time series. Futoma et al. [34] extends the GP-based approach to imputation in the multivariable setting by applying the multi-task Gaussian Processes (MTGP) prediction [35]. However, the quality of estimating missing values relies on the estimation of covariance structure among variables when using MTGP or other multi-task functional approaches [36, 37, 38]. To make a confident estimation of the covariance, a large amount of time points with shared observations of these variables are often required by these multi-task approaches [35, 39]. Due to the fact that many patients only have records with a limited number of time points, time series of inpatient clinical laboratory tests fall short of such a requirement. Hori et al. [40] extend MTGP to imputation on multi-environmental trial data in a third-order tensor. However, this approach is not applicable to clinical data with a large number of patients, since the covariance matrix of observed values needs to be computed together with its inverse, which is intractable for our datasets.

Recently, Luo et al. [41] explore the application of GPs in clinical data and propose an algorithm, 3-dimensional multiple imputation with chained equations (3D-MICE), that combines the imputation from GP and traditional MICE based on weighting equations. However, the weighting equations are calculated only based on the standard deviations of GP and MICE predictions for missing values. The weighting strategy is static and not optimized. We postulate that calculating weights through an optimization problem can help to improve the imputation quality. In our work, instead of the predictive mean matching used in [41], we choose linear regression as one component of our model. Our model is also based on a statistical model and thus statistically justified which is not the case for [41].

The mixture model proposed in this work differs from the mixtures of Gaussian processes (MGP) [42] in several ways. The mixing components in our model are not limited to Gaussian process, thus bringing in possibilities of capturing correlations with various choices of parametric functions. More importantly, our model is designed for multivariable time series and captures both temporal and cross-sectional correlations. The latter is ignored in [42] when MGP is applied to multivariable time series, which yields less effective predictions since some clinical tests have strong cross-sectional relations. Without a Gaussian process component, the standard Gaussian mixture model (GMM) is well studied and has been extended to imputation tasks [43, 44, 45, 46]. These approaches still suffer from losing informative temporal correlations when imputing missing values in longitudinal data.

In order to effectively model the interaction between different aspects, we represent the data as a tensor with each aspect being one mode. Therefore, our method is also considered as a tensor completion approach. Tensor completion problems are extensively studied. However, the classic tensor completion methods [47, 48, 49] focus on general tensors and usually do not consider temporal aspects. In recent years, many studies explore the application of temporal augmented tensor completion on imputing missing values in time series [50, 51, 52, 53, 54]. These methods discretize time into evenly sampled intervals. However, due to the fact that inpatient clinical laboratory tests are usually measured at varying intervals, assembling clinical data over regularly sampled time periods might have several drawbacks, such as leading to sparse tensors if discretizing time at fine granularity (e.g. every minute) while some laboratory tests are measured less frequently (e.g. daily). Furthermore, extending these methods to the case, where time is not regularly sampled, is not easy and straightforward, requiring changing design details and the objective functions to be optimized. Recently, Yang et al. [55] propose a tensor completion method that can deal with irregularly sampled time. However, most components of this approach are tailored to the characteristics of septic patients. We implement this approach with only the component of time-aware mechanism that is general and applicable to our experimental settings. However, this approach is not so effective in our experiments, thus not included as a benchmark in this paper. Lately, Zhe et al. [56] propose a Bayesian nonparametric tensor decomposition model that captures temporal correlations between interacting events. However, this approach is not directly applicable to continuous multivariable time series because it focuses on discrete events and captures temporal correlations between the occurrences of events.

Recurrent Neural Networks (RNNs) are another choice to capture temporal relationships and have been studied in imputation tasks. Early studies [57, 58, 59] in RNN-based imputation models capture temporal correlations either within each time series or across all time series [60]. Yoon et al. [60] recently propose a novel multi-directional RNN that captures temporal correlations both within and across time series. Compared with traditional statistical models, neural network-based models are usually harder to train and provide limited explainability, while the performance is not always necessarily better. We use [60] as a benchmark since it is the most recent model. Besides imputation, researchers also attempt to build classification models with RNNs that utilize the missing data patterns [61, 62, 63]. Instead of relying on the estimation of missing values, these models perform predictions by exploring the missingness itself and no separate imputation is required. These models can not be a substitute for imputation methods since they do not provide imputed values.

III Methodology

III-A Imputation Framework

In many predictive tasks on temporal clinical data, time series are often aligned into the same-length sequences to derive more robust patient phenotypes through matrix decomposition or discover feature groups by applying sequence/graph mining techniques [55, 64]. We model under this assumption. In this work, we use tensor representation, in which patients have the same number of time points. We represent the data collected from $P$ patients with $V$ laboratory tests and $B$ time points as two 3D tensors $\mathcal{X}\in\mathbb{R}^{P\times V\times B}$ and $\mathcal{T}\in\mathbb{R}^{P\times V\times B}$, shown in Fig. 1. Each laboratory test measurement $x_{p,v,b}$ is stored in the measurement tensor $\mathcal{X}$. Each time $t_{p,v,b}$, when $x_{p,v,b}$ is measured, is stored in the time tensor $\mathcal{T}$.

Table I lists the main symbols we use throughout the paper. Missing values in the measurement tensor are denoted as $x^{mis}_{p,v,b}$, showing that the value of test $v$ at time index $b$ for patient $p$ is missing. Correspondingly, $x^{obs}_{p,v,b}$ denotes an observed value. The time tensor $\mathcal{T}$ is complete, since we only collect patient records at the time when at least one laboratory test measurement is available. We assume we know the “prescribed” time a missing measurement should have been taken. In the matrix $x_{:,:,b}$ at time index $b$, the measurement time $t_{p,v,b}$ and $t_{q,v,b}$ can be different when $p\neq q$, whereas for a given patient $p$, we have $t_{p,v,b}=t_{p,u,b}$ for $v,u\in[1:V]$. That is, all tests for a particular patient are taken at the same time.

Disregarding the temporal dimension, the imputation problem is well studied. If one dimension is time, in order to apply imputation methods that are not designed for time series, we need to disregard the temporal aspect or ignore temporal correlations of the data. However, temporal trajectories can reveal patient’s underlying pathophysiological evolution, the modeling of which can help to better estimate missing values. For the reason that both cross-sectional information and temporal aspects can impact the estimation of missing measurements, we explore mixture models, which are composed of several base models through either a cross-sectional or temporal view. We introduce these base models in Section III-B.

In our imputation framework, a mixture model is trained for each variable and time index. We use $Mix_{v,b}$ to denote the mixture model to impute missing values of variable $v$ at time index $b$. The missing values $x^{mis}_{:,v,b}$ in the fiber $x_{:,v,b}$ are imputed by the optimized $Mix_{v,b}$. In each iteration of the algorithm, it is assumed that all other values are known and only $x^{mis}_{:,v,b}$ is imputed. This is wrapped in an outer loop. We call this procedure simple-pass tensor imputation, i.e., one pass through all $v,b$. Since several simple-pass tensor imputations are conducted, our approach is also considered as an iterative imputation [65], which can also be regarded as a sampling-based approach where a Gibbs sampler is used to approximate convergence. The convergence of iterative imputation methods can be quite fast with a few iterations [8].

In detail, the iterative imputation approaches start by replacing all missing data with values from simple guesses; we fill in all missing values with initial estimates by taking random draws from observed values. This procedure is called an initial imputation. Then we perform iterative tensor imputation on each copy separately. The training procedure and imputation for fibers are introduced in Section III-C and III-D. We also rely on the concept of multiple imputations, where several iterative imputations are performed and the imputed values are averaged at the end. Each iterative imputation starts with a different iterative imputation tensor and/or uses a different order of $v$,$b$.

In summary, the algorithm creates $M$ different copies $\mathcal{X}^{1},\ldots,\mathcal{X}^{M}$ of $\mathcal{X}$, each one filled with different random $\mathcal{X}^{mis}$. For each $i=1,\ldots,M$, we then perform $K$ simple-pass tensor imputations. Each simple-pass has a loop over all $v$,$b$, which uses $Mix_{v,b}$ to adjust $\mathcal{X}^{i,mis}$. At the end of the $M$ imputations, $\mathcal{X}^{i,mis}$ are averaged across all $i$ to yield the final imputed tensor.

The whole imputation process involves $M\times K\times V\times B$ imputation models. We next first focus on the base models behind $Mix_{v,b}$ and then on the actual mixture model.

III-B Base Models

Our mixture models are composed of three components that are derived from two base models, linear regression and Gaussian processes. One component consists of GP models and the other two components are linear models through two different views of the measurement tensor. Through a cross-sectional view, the tensor can be considered as a vector of patient-by-variable matrices at different time indices. Through a temporal view, we can view the tensor as a vector of patient-by-time matrices for different variables.

III-B1 Linear model through cross-sectional view

We can view the measurement tensor $\mathcal{X}$ as a vector of patient-by-variable matrices. On the slice $x_{:,:,b}$, we use a linear regression model to fit the target variable $v$ as a function of the other variables except $v$. The target values $x_{:,v,b}$ are modeled as

[TABLE]

where $\beta^{(1)}_{v,b}$ is the column vector of coefficients and $\sigma^{(1)}_{v,b}$ is the standard deviation of the error $\epsilon^{(1)}_{v,b}$, regarding the regression model through cross-sectional view for variable $v$ and time index $b$.

The likelihood distribution of $x_{:,v,b}$ is then given by

[TABLE]

The training data consists of observed target values $(x_{p,v,b})_{p\in P^{tr}_{v,b}}$ and the input data $(x_{p,-v,b})_{p\in P^{tr}_{v,b}}$, where $P^{tr}_{v,b}$ is the training patient set and includes $p$ only if $x_{p,v,b}$ is observed.

III-B2 Linear model through temporal view

In addition to the cross-sectional view, we can also view the measurement tensor $\mathcal{X}$ as a vector of patient-by-time matrices. On matrix $x_{:,v,:}$, we use linear regression to model the measurements at time index $b$ against those at other indices.

The target values $x_{:,v,b}$ are modeled as

[TABLE]

where $\beta^{(2)}_{v,b}$ is the column vector of coefficients and $\sigma^{(2)}_{v,b}$ is the standard deviation of the error $\epsilon^{(2)}_{v,b}$, regarding the linear regression model though temporal view for variable $v$ and time index $b$.

The likelihood distribution of $x_{:,v,b}$ is given by

[TABLE]

The training data consists of observed target values $(x_{p,v,b})_{p\in P^{tr}_{v,b}}$ and input data $(x_{p,v,-b})_{p\in P^{tr}_{v,b}}$.

III-B3 Gaussian processes through temporal view

Gaussian processes are commonly used to capture trajectories of variables, thus used in our mixture model to capture temporal correlations. Through the same temporal view as introduced above, on matrix $x_{:,v,:}$, we fit GPs on time series for each patient.

The target value $x_{p,v,b}$ is modeled as

[TABLE]

where $\mu_{p,v,b}$ is the overall mean of the model, $f(\cdot)$ is a Gaussian process with mean of [math] and a covariance matrix $\mathcal{K}(t,t^{\prime})$ of time pairs ($t$,$t^{\prime}$). Then the likelihood distribution of $x_{p,v,b}$ is written as

[TABLE]

where $\theta_{v,b}$ are the kernel parameters of the GP models, and the predictive mean and variance are given by $m^{(3)}(\cdot)$ and $\Sigma^{(3)}(\cdot)$; see more details in Appendix B. For a certain $v$ and $b$, all GP models share the same kernel parameters $\theta_{v,b}$.

III-C The Mixture Model

Given the likelihood distribution of all three components, we model the joint mixture distribution, regarding the variable $v$ and time index $b$, in the following way

[TABLE]

where we define $V_{v,b}=[x_{:,-v,b}\ x_{:,v,-b}]$, $\boldsymbol{\alpha}_{v,b}=(\alpha_{p,v,b})_{p\in P}$, $\delta^{(k)}_{v,b}=(\boldsymbol{\mu}^{(k)}_{v,b},\boldsymbol{\Sigma}^{(k)}_{v,b})$ and $\pi^{(k)}_{v,b}$ is the mixing weight for the $k$th component. This model can be interpreted as the joint distribution between observed data $V_{v,b}$ and missing values $x_{:,v,b}$, which consists of a mixture of three distributions. The first one $p_{1}(x_{:,v,b},V_{v,b})$ is modeled as

[TABLE]

the remaining two follow the same logic.

By marginalizing over $x_{:,v,b}$, the prior probability distribution $p(V_{v,b})$ is written as

[TABLE]

which is a mixture of Gaussians. It also follows that

[TABLE]

where we define $C^{(k)}_{v,b}=\frac{\pi^{(k)}_{v,b}\mathcal{N}(V_{v,b}|\delta^{(k)}_{v,b})}{\sum_{j=1}^{3}\pi^{(j)}_{v,b}\mathcal{N}(V_{v,b}|\delta^{(j)}_{v,b})}$.

We train our mixture model on observed target values by maximizing the log likelihood of the joint mixture distribution

[TABLE]

where $V^{tr}_{v,b}$ is the training input data and defined as the concatenation of $(x_{p,-v,b})_{p\in P^{tr}_{v,b}}$ and $(x_{p,v,-b})_{p\in P^{tr}_{v,b}}$, and $\gamma_{v,b}$ is the set of all trainable parameters

[TABLE]

After training, the missing values are imputed using individualized (per-patient) mixing weights that are derived from the conditional distribution. The missing value $x^{mis}_{p,v,b}$ is imputed by

[TABLE]

where the individualized mixing weight of the $k$th component for patient $p$, variable $v$ and time index $b$ is defined as

[TABLE]

where $V_{p,v,b}$ is the observed data and defined as the concatenation of $x_{p,-v,b}$ and $x_{p,v,-b}$.

III-D Mixture Parameter Estimation

Let $\ell(\gamma_{v,b})$ be the log likelihood $\ln{p(x^{obs}_{:,v,b},V^{tr}_{v,b};\gamma_{v,b})}$. Explicitly maximizing $\ell(\gamma_{v,b})$ is hard. Instead, we use the EM algorithm to repeatedly construct a lower-bound on $\ell(\gamma_{v,b})$ and then optimize that lower-bound $\mathcal{L}(\gamma_{v,b})$. We first define a latent indicator variable $q_{v,b}\in\{1,2,3\}$ that specifies which mixing component that data points come from. Then we use Jensen’s inequality to get the lower-bound $\mathcal{L}(\gamma_{v,b})$, which is given by

[TABLE]

where

[TABLE]

In (12) and (13), the marginal distribution $p(q_{v,b}=k)$ over $q_{v,b}$ is specified by the mixing coefficients $\pi^{(k)}_{v,b}=p(q_{v,b}=k)$. We can view $Q_{p}(q_{v,b}=k)$ as the responsibility that component $k$ of the mixture model $Mix_{v,b}$ takes to “explain” $x_{p,v,b}$. We use the standard EM algorithm to maximize the lower-bound $\mathcal{L}(\gamma_{v,b})$; see more details about the estimation of the parameters in Appendix A.

III-E The Imputation Model

The proposed mixture model provides flexibility in changing base models and can be solved under various assumptions. For example, a full mixture model consists of all three base models introduced earlier. However, under the assumption of only linear correlations within time series, mixture models composed of only the two linear components (denoted as MixMI-LL) can be used. Although the choice of mixtures can be determined empirically or based on expert knowledge of the data, we propose an imputation model (denoted as MixMI) that chooses the type of mixture automatically.

In MixMI, for each variable and time index, we train both the MixMI-LL and the full mixture model, then select the one with lower training error as the final mixture model that is used to perform imputation in inference. We use the absolute training error, instead of likelihood, as the selection criterion because the likelihood of a mixture model with two linear models is not of the same scale as the full mixture model.

IV Experimental Setup

IV-A Real-world Datasets

We collect two real-world datasets from the Medical Information Mart for Intensive Care (MIMIC-III) database [66] and the Northwestern Medicine Enterprise Data Warehouse (NMEDW). Each dataset contains inpatient test results from 13 laboratory tests. These tests are quantitative, frequently measured on hospital inpatients and commonly used in downstream classification tasks [67, 63]. They are the same as those used in [41] in their imputation study. We organize the data by unique admissions. We distinguish multiple admissions of the same patient. Each admission consists of time series of the 13 laboratory tests. Although our model is evaluated only on continuous variables, they are also applicable on categorical variables with a simple extension (e.g. through predictive mean matching [68]).

In both MIMIC-III and NMEDW datasets, the length of time series varies across admissions. To apply our imputation model on these datasets, we truncate time series so that they have the same length. The length is the average number of time points across all admissions. Before truncating, the average number of time points in MIMIC-III dataset is 11. We first exclude admissions that have less than 11 time points, and then we truncate time series by removing measurements taken after the 11-th time point. We also exclude admissions that contain time series that have no observed values. Our MIMIC-III dataset includes 26,154 unique admissions and the missing rate is about 28.71%. The same data collection procedure is applied on the NMEDW dataset (approved by Northwestern IRB) where we end up with 13,892 unique admissions that have 7 time points. The missing rate of the NMEDW dataset is 24.22%.

IV-B Synthetic Datasets

We create synthetic datasets to explore and discover the properties of the data that might benefit our model and/or comparison models. In synthetic datasets, we augment the correlation between measurements and times. We do not augment correlations by imposing strong constraints on time series where closer measurements have closer values. Instead, we generate synthetic times by altering real times so that the constraints in synthetic data are “slightly” stronger than real-world data. We also introduce a scaling factor $d$ to control the strength of the constraints in synthetic data.

The synthetic datasets are generated based on the real-world MIMIC-III dataset. We move two consecutive times of a time series closer, if the relative difference $\Delta\tilde{x}$ in two consecutive measurements is smaller than the relative difference $\Delta\tilde{t}$ in two consecutive times. The relative differences $\Delta\tilde{x}$ and $\Delta\tilde{t}$ of a time series are given by

[TABLE]

The scaling factor $d\in(0,1)$ controls how much we move times. If $d=0$, we do not move times. In other words, the synthetic dataset at $d=0$ is the same as the real-world MIMIC-III dataset. As $d$ increases, stronger constraints are introduced to synthetic data. The synthetic time $t^{\prime}$ for a time series is generated as follows:

[TABLE]

If a time series has missing values, we first calculate the synthetic times for the observed measurements. Then we perform a linear interpolation between real times and synthetic times for observed measurements to generate synthetic times for missing measurements.

IV-C Evaluation of Imputation Quality

We randomly mask 20% observed measurements in a data set as missing and treat the masked values as the test set. The remaining observed values are used in training. We impute originally missing and masked values together, and compare the imputed values with the ground truth for masked data to evaluate imputation performance. The Mean Absolute Scaled Error (MASE) [69] is used to measure the quality of imputation on the test set. MASE is a scale-free measure of the accuracy of predictions and recommended by Hyndman et al [69, 70] to measure the accuracy of predictions for series. We calculate MASE for all variables and take a weighted average according to the number of masked values of a variable to get an overall MASE per dataset.

Let $mask_{p,v}$ be the set of cardinality $I_{p,v}$ of all time indices that have been masked for patient $p$ and variable $v$. Also let $Y_{p,v}=(x^{obs}_{p,v,j})_{j}$ be the sequence of length $J_{p,v}$ of all observed values for patient $p$ and variable $v$, and let $\tilde{x}_{p,v,i}$ represent the imputed value. The MASE for variable $v$ is defined as

[TABLE]

To show the effectiveness of our imputation model MixMI and its variant MixMI-LL, we compare MASE scores of MixMI and MixMI-LL with other six imputation methods: (a) MICE [8] with 100 imputations, where the average of all imputations are used for evaluation; (b) the pure Gaussian processes, where a GP model is fitted to the observed data of each time series using GPfit [71] in R and missing values are replaced with the predictions from the fitted GP models; (c) 3D-MICE, a state-of-the-art imputation model [41] for which we obtain their code and adapt it to account for our use of the tensor representation; (d) multi-task Gaussian process prediction (MTGP) [35]; (e) a Gaussian mixture model for imputation (GMM) [43, 72] and (f) M-RNN, a state-of-the-art RNN-based imputation model [60].

To tune hyperparameters if any in these models, we mask out 20% observed measurements in the training set as a validation set and tune hyperparameters on the validation set. We introduce $\pi^{(1)},\pi^{(2)}$ and $\pi^{(3)}$ as hyperparameters of our imputation model. In our experiments, all mixture models start with the same mixing weights, i.e., $\pi^{(k)}_{v,b}=\pi^{(k)}$ for $k=1,2,3$. Other trainable parameters are initialized directly from the data after initial imputation and do not require a manual specification.

We also evaluate our imputation model under a classification task on our MIMIC dataset and compare it with GRU-D [63], a state-of-the-art RNN model for classifications with multivariable time series, which does not require a separate data imputation stage. We build a RNN model with the standard Gated Recurrent Unit as our classification model (GRU-MixMI) trained on data imputed by MixMI. GRU-D is trained on data without imputation. We predict whether a patient dies within 30 days after admission, instead of 48 hours as in [63], because of severe imbalance in our cohort, where only 0.57% patients have positive labels within 48 hours.

V Results

V-A Performance Comparison

V-A1 Imputation task

Table II and Fig. 2 compare all imputation models on 4 datasets. Table II shows the overall MASE score of each imputation model on all datasets. Fig. 2 provides a comparison for all imputation models in MASE score over 3D-MICE by showing the percentage deviation against 3D-MICE. We select 3D-MICE since it is a state-of-the-art benchmark model. We observe that MixMI outperforms all comparison models on all 4 datasets and is significantly better than the second best model ($p$=.001, permutation test with 1000 replicates) on all 4 datasets.

Table III shows a variable-wise comparison of the imputation models on the real-world MIMIC ($d$=0) dataset. Our two imputation models outperform all comparison models on all variables. MixMI is better than MixMI-LL on most variables. All models except GP and MTGP achieve a much lower error on Hematocrit and Hemoglobin than on other variables. The reason is that these two variables are highly correlated. Those methods that capture the correlation between variables can reasonably infer missing values for Hematocrit from observed measurements of Hemoglobin, and vice versa. Compared to MICE, MixMI achieves even lower errors on these two variables, which indicates that temporal correlations captured by our model help to make better estimation of missing values, even when there is a more dominant cross-sectional correlation.

As shown in Table II, all models except MTGP benefit from the increment of $d$, the scaling factor used to generate synthetic data. The reason is that all models take into account temporal aspects and the measurements in the synthetic time series have stronger temporal correlations as $d$ increases. MICE and GMM also benefit from the temporal correlations because we include time as a feature in addition to variable features. We also try to exclude times, however, experiments show that these two models perform better when times are included. MixMI-LL is outperformed by 3D-MICE when $d$ increases to $1$, while MixMI shows its robustness to the variation of $d$ in our current experimental settings.

We compare the running times of all models on the real-world MIMIC dataset. MixMI-LL (taking 4.2 hours), GP (1.1 hours), MTGP (1.2 hours) GMM (2.2 hours) and M-RNN (0.9 hours) are the fastest models; 3D-MICE (156.1 hours) is the slowest; MICE (77.5 hours) and MixMI (109.5 hours) are in the middle.

V-A2 Classification task

Our classification model GRU-MixMI is compared with GRU-D on data with different folds and different seeds. First, we make a copy of the original data and split each copy into the same training and test set. One copy is then imputed by MixMI and used to train GRU-MixMI. The other is used by GRU-D directly without imputation. Both models are trained and tested 5 times with different random seeds pertaining to algorithms. The area under the ROC curve (AUC) scores on test sets are reported. This procedure repeats 3 times with different splits. On average, GRU-MixMI ties with GRU-D in terms of AUC scores ($0.7857$ vs. $0.7853$, $p>0.5$). However, compared with models that perform imputation and prediction together, our imputation model, as a pure imputation method, offers more opportunities for other supervised and unsupervised tasks, and modeling choices and effective feature engineering.

V-B Component Weights

The key information our imputation model uses to estimate missing values is the component weight $\pi$, which quantifies the interaction of cross-sectional and temporal correlations. We study MixMI on the real-world MIMIC dataset and visualize how the interaction is captured by our model. Fig. 3 shows a comparison of all component weights for variables across all patients and times. When imputing Hematocrit and Hemoglobin, our model relies mostly (79.6% and 84.8%) on the linear model through the cross-sectional view because of the strong correlation between these two variables. Interestingly, most values of $\pi^{(3)}$ are lower than 10%, which indicates that temporal correlations are not strong in these variables in our dataset. However, our model is still able to detect the week temporal correlations and utilizes them to improve imputation. Additionally, we observe that when predicting the missing values at the beginning and the end of a time series, our model reasonably uses a lower value of $\pi^{(3)}$ than in the middle. On average, $\pi^{(3)}$ of the first and the last time indices is 13.9% lower than those in the middle. This is due to the usage of GPs, which usually produce less confident estimates at the end points of series. Furthermore, we observe an increment of $\pi^{(3)}$ as we impose stronger temporal correlations in synthetic datasets, which further validates the ability of our model in capturing such interaction.

V-C Individualized Weights

By introducing individualized (per patient) mixing weights $\Pi$ defined in (11), we improve the performance in MASE score from 0.08351 to 0.07538, an improvement of 9.73% compared against the model where each mixture component has a fixed weight for all patient cases. The reason that individualized weights are better than fixed weights in our model might be that they better approximate the responsibilities $Q$ defined in (13).

In training, we can optimize the responsibility a component should take to “explain” an observed target value $x_{p,v,b}$ for $p\in P^{tr}_{v,b}$. However, when making inference, the responsibility each component should take to “explain” a missing value is unknown, because responsibilities depend on observed target values, according to (13). We have to use $\Pi$, the individualized mixing weights, as an approximation of the responsibilities in inference. As defined in (11), $\Pi$ only depends on the inputs, therefore, we can calculate them when making inferences on the test set.

In a standard mixture model, we could use $\pi^{(k)}_{v,b}$, which is the average of responsibilities of the $k$th component across all training patients, as a fixed weight that the $k$th component should contribute to impute missing values $x^{mis}_{:,v,b}$ for all test patients. However, patient time series can be very different and the confidence of predictions by the GP component can vary largely across different patient cases. A fixed weight can not reflect such variation in prediction confidence.

We shall view an individualized mixing weight as an approximation of how much responsibility a component should take to impute the missing value for a particular patient case. It is tailored for each patient. To visualize how individualized mixing weights help to produce better estimates, in Fig. 4, we plot the distribution of the individualized weights $\Pi$ of the GP component in the training set and compare it with the distribution of the optimized responsibility values $Q$. The responsibilities that the GP component should take can vary a lot in different patient cases, especially on the synthetic dataset, which implies that it is more reasonable for patients to get individualized mixing weights than a fixed weight. We also observe that the individualized mixing weights reasonably mimic the distribution of the optimized responsibilities on the training set. The improvement of our model on the test set attests that the individualized weights approximate the responsibilities better than fixed weights.

VI Conclusions

We present and demonstrate a mixture-based imputation model MixMI for multivariable clinical time series. Our model can capture both cross-sectional and temporal correlations in time series. We integrate Gaussian processes with mixture models and introduce individualized mixing weights to further improve imputation accuracy. We show that MixMI can provide more accurate imputation than several state-of-the-art imputation models. Although in this work our model is tested on inpatient clinical data, the proposed idea is general and should work for all multivariable time series where interactions of cross-sectional and temporal correlations exist.

Appendix A Parameter Estimation in EM

In the E (Expectation) step, we calculate the responsibilities $w^{(k)}_{p,v,b}=Q_{p}(q_{v,b}=k)$ for $p\in P^{tr}_{v,b}$ using the current values of the parameters in iteration $j$:

[TABLE]

Let $Z_{v,b}=(x_{p,-v,b})_{p\in P^{tr}_{v,b}}$ and $Y_{v,b}=(x_{p,v,-b})_{p\in P^{tr}_{v,b}}$. In the M (Maximization) step, we re-estimate the parameters in iteration $(j+1)$ using the $j$th responsibilities:

[TABLE]

where $[\mathbf{w}^{(k)}_{v,b}]^{(j)}$ is the vector of $[w^{(k)}_{p,v,b}]^{(j)}$ for $p\in P^{tr}_{v,b}$ in iteration $j$. The kernel parameters $\theta_{v,b}$ of GP models are evaluated by function $G$, a gradient descent method that calculates the estimates of $[\theta_{v,b}]^{(j+1)}$ to maximize $\mathcal{L}_{v,b}(\gamma)$, using $[\theta_{v,b}]^{(j)}$ as the starting point. The first order derivatives of $\mathcal{L}_{v,b}(\gamma)$ with respect to $\theta_{v,b}$ that are used in $G$ are given in Appendix C.

Appendix B GP Model

We assume the GP model discussed here in a mixture model for a certain variable and time, and thus we exclude the subscripts $v$ and $b$. We use $x_{p,t}$ to denote a measurement of the time series $x_{p}$ at time $t$ for patient $p$ of a certain variable. We use $x_{p,-t}$ to denote a time series without the measurement at time $t$. The GP model is given by

[TABLE]

where $\mu_{p,t}$ is the overall mean of the model and $f(t)$ is a Gaussian process with mean of [math] and covariance of $\mathcal{K}(t,t^{\prime})$. Following the maximum likelihood approach, the best linear unbiased predictor (BLUP) 222Sacks, Jerome, et al. ”Design and analysis of computer experiments.” Statistical science (1989): 409-423. at $t$ and the mean squared error are

[TABLE]

where $r_{t}(t^{\prime})=corr(f(t),f(t^{\prime}))$, $r$ is the vector of $r_{t}(t^{\prime})$ for all possible $t$, $\bar{t}$ is a vector of time except for time $t$, $R$ is the $(B-1)\times(B-1)$ correlation matrix and the correlation function is given by

[TABLE]

The estimator $\sigma^{2}$ is given by

[TABLE]

where $1_{n}$ is a vector with length $(B-1)$ of all ones.

Appendix C Partial Derivatives in GP

To simplify the notations, we assume that the likelihood function $L$ under consideration is for a mixture model for a certain variable and time. The partial derivative with respect to Gaussian process parameters $\theta$ is

[TABLE]

Letting $g_{p}(\theta)=m^{(3)}(\theta,x_{p,-t},\bar{t})$ and $h_{p}(\theta)=\Sigma^{(3)}(\theta,x_{p,-t},\bar{t})$, we have

[TABLE]

Then $\frac{\partial g_{p}(\theta)}{\partial\theta}$ and $\frac{\partial h_{p}(\theta)}{\partial\theta}$ are given by

[TABLE]

where $H_{1}$, $\frac{\partial H_{1}}{\partial\theta}$, $H_{3}$ and $\frac{\partial H_{3}}{\partial\theta}$ are given as follows:

[TABLE]

[TABLE]