# An Unsupervised Multivariate Time Series Kernel Approach for Identifying Patients with Surgical Site Infection from Blood Samples

## Abstract

A large fraction of the electronic health records consists of clinical measurements collected over time, such as blood tests, which provide important information about the health status of a patient. These sequences of clinical measurements are naturally represented as time series, characterized by multiple variables and the presence of missing data, which complicate analysis. In this work, we propose a surgical site infection detection framework for patients undergoing colorectal cancer surgery that is completely unsupervised, hence alleviating the problem of getting access to labelled training data. The framework is based on powerful kernels for multivariate time series that account for missing data when computing similarities. Our approach show superior performance compared to baselines that have to resort to imputation techniques and performs comparable to a supervised classification baseline.

###### keywords:

Surgical site infection, Electronic health records, Multivariate time series, Kernel methods, Missing data, Unsupervised learning## 1 Introduction

Surgical Site Infection (SSI) is one of the most common types of nosocomial infections lewis2013 (), representing up to 30% of hospital-acquired infections magill2012prevalence (); de2009surgical (). SSI can be divided into different types depending on the anatomical location of the infection ko2015american (). Superficial infections can be treated with local surgical debridement and antibiotics. On the other hand, deep infections are more complex and require lapratomies and/or percutaneous drainage and intravenous antibiotics. Recently, SSI risk factors such as advanced age, overweight, smoking, open surgery or disseminated cancer have been reported lawson2013risk (). Depending on the location of the infection, some factors contribute more to increase the risk of SSI. For example, longer lasting surgeries are associated with deep SSI, whereas a high body mass index is related with superficial SSI lawson2013risk (); blumetti2007surgical ().

Along with increased mortality rate (3%), SSI also prolongs hospitalization up to two weeks whitehouse2002 () and increases the risk of readmission shah2017evaluation (). This, in average, doubles the expenses per patient and increases the chances of readmission, with further additional costs up to 27,000 USD whitehouse2002 (); owens2014surgical (). Hence, a reduction in the number of postoperative complications like SSI will be of great benefit both for the patients and the society.

Many recent studies have focused on the analysis of SSI from blood tests, both before and after surgery silvestre2014 (); soguero2015data (); medina2016 (); angiolini2016 (). The advantage of using blood tests for this purpose is that they are recorded frequently with low burden for the patients and contain much information about their actual health status. Among blood tests results, the high predictive power of C-reactive protein (CRP) test has been evaluated and emphasized in several works. For example, authors in angiolini2016 () and gans2015diagnostic (), demonstrated the relation between different CRP cutoff values and the risk of SSI on postoperative days 3 and 4. Others have combined the results from blood tests with other structured and unstructured clinical data to predict SSI complications soguero2016support (); hu2017strategies (); Sanger2016259 ().

Since blood samples are collected over time they can naturally be represented as time series. Such clinical time series data have some special characteristics that distinguish them from time series from other domains. One key property is missing data, which can occur because of e.g. lack of documentation or lack of collection wells2013strategies (). Another characteristic is that the time series usually are multivariate. For example could the patients be described by the measurements of many different blood tests, such as e.g. hemoglobin, CRP, etc., where many of them exhibit relationships or dependencies, which can be non-trivial (non-linear).

In order to apply a machine learning algorithm on such time series, one can work directly in the input space using a similarity measure that accounts for dependencies in time and among the variables mikalsen2017time (). An advantage is that a time consuming feature learning process, or a manual feature design process that requires user intervention and domain expertise, can be then avoided. On the other hand, a key problem with classical similarity measures is that most of them are not able to cope with missing data in clinical time series data without applying a preprocessing step such as imputation in order to obtain complete data mitsa2010temporal (). However, important information about the clinical condition of the patient and the decisions of the caregiver may disappear in such a preprocessing step since the information that some data are missing is lost, and replaced by biased estimates.

In addition to the problem of analyzing multivariate time series (MTS) containing missing data, another key challenge associated with data originating from electronic health records (EHRs) is that getting access to labelled data for training the machine learning models is difficult. It is well known that manual annotation of labels, especially in the healthcare domain, is a cumbersome process that could be both time consuming and expensive MIKALSEN2017105 (), since clinical expertise is needed to create the training sets. Indeed, the workload for the clinicians is tremendous already, and with an aging and more diseased population we cannot expect these tasks to be prioritized in the future. To overcome this annotation problem, recently there have been proposed several different methods within the framework of semi-supervised learning chapelle2009semi (), of which the so-called anchor learning is an example MIKALSEN2017105 (); Halpernocw011 ().

In this work, we propose a disease classification framework for time series originating from EHR data such as e.g. blood tests. We address the key challenges described above by taking an unsupervised approach where we utilize powerful kernels for MTS containing missing data. These kernels have emerged due to recent advances in time series analysis as new family of methods that can cope with incomplete and multivariate data. Prominent examples include the learned pattern similarity (LPS) baydogan2016time () and the time series cluster kernel (TCK) mikalsen2017time (). These methods are able to handle missing data without having to resort to imputation methods. Another important property of the LPS and TCK kernels is that they can be trained in an unsupervised way and do not require tuning of critical hyperparameters. We take advantage of their robustness by using the kernels as input to an unsupervised spectral clustering framework for identifying patients with SSI. By doing so, the entire framework is grounded within the theoretically well understood kernel methods. Moreover, spectral clustering is considered a state-of-the-art clustering algorithm and has been successfully utilized in many applications ng2002spectral (); krzakala2013spectral (); lokse2017spectral ().

The proposed methodology consists of two steps, namely to first compute the kernel and then apply spectral clustering (see Fig. 1). We use the proposed framework to identify patients undergoing colorectal cancer surgery with SSI based on only blood samples, and illustrate its power by comparing to a similar framework where we use kernels that have to resort to imputation techniques. In addition, we compare to a supervised baseline for detecting SSI.

## 2 Materials and methods

### 2.1 Data description

Ethics approval for the present study was obtained from the Data Inspectorate and the Ethics Committee at the University Hospital of North Norway (UNN) jensen2017analysis (); soguero2016predicting (). The dataset we consider consists of 7741 patients that underwent a gastrointestinal surgical procedure at UNN in the years 2004–2012. The SSI persistent in-hospital morbidity is particularly associated with colorectal cancer surgery lawson2013risk (); blumetti2007surgical (); lawson2013reliability () and therefore patients that did not undergo this type of surgery were excluded, reducing the size of the cohort to 1137 patients.

The International Classification of Diseases (ICD10) or NOMESCO Classification of Surgical Procedures (NCSP) codes related to severe postoperative complications, both superficial and deep infections, were considered to extract a cohort. For the purpose of this work and similarly to earlier studies LIMON2014127 (); Gibbons20111 (); horan_gaynes_martone_jarvis_emori_1992 (); BERGER2013974 (); soguero2015data (); Sanger2016259 (), we do not distinguish between deep and superficial SSI.

In collaboration with the clinician (author A. R.), a set of 11 different types of blood tests were defined as clinically relevant and extracted for all patients from their EHRs, namely, hemoglobine, leukocytes, CRP, potassium, sodium, creatininium, thrombocytes, albumin, carbamide, glucose and amylase. The blood tests can be considered as continuous variables over time and hence represented as MTS. For the purpose of the current analysis, we discretize time and let each time interval be one day. However, all 11 blood tests are not available every day for each patient, meaning that the dataset contains missing data. We focus on classification of SSI within 20 days after surgery and therefore define the postoperative window as the period from postoperative day 1 until day 20. Patients with less than two measurements during the postoperative window are removed from the cohort, which leads to a dataset with 232 infected patients (cases) and 651 non-infected ones (control).

### 2.2 Strategies for dealing with missing data

Since EHR data might be missing for several different reasons, in many cases the missingness mechanism cannot be described exclusively as either missing completely at random (MCAR), missing at random (MAR) or missing not at random (MNAR), but rather as a combination of all these three traditional schemes wells2013strategies (); hu2017strategies (). Indeed, one reason for missing data is lack of documentation, which occurs for example when a clinician orders a blood test but the test for some reason is not performed, or it is performed but not documented because of an error by the physician that performs the test or the data is lost during extraction. Such type of missing data is usually MCAR or MAR. However, missingness can also be caused by lack of collection. This happens, for example, when the clinician that is treating the patient thinks the health condition of the patients is so good that there is no reason to order a blood test on that particular day. In this case, data are MNAR.

Most machine learning methods, and in particular discriminative learning algorithms, work only with complete datasets little2014statistical (). However, in a clinical setting, creating a complete dataset by simply discarding patients with missing data may lead to incorrect assessments or prognostics, since the fraction of missing data is typically large soguero2016predicting (). Hence, one could end up discarding a lot of information and get a weak model. Moreover, complete case analysis as a result of discarding data only give unbiased predictions if the missingness mechanism is MCAR rubin1976inference (); schafer1997analysis ().

As an alternative, a preprocessing step involving imputation of missing values with some estimated value is common. For example one could fill the missing values for each variable of interest with the mean or median value of the observed samples. A simple, but sometimes efficient approach, is to impute all missing values with zeros. Other so-called single imputation methods include machine learning based methods such as multilayer perceptrons, self-organizing maps, k-nearest neighbors, recurrent neural networks and regression-based imputation garcia2010pattern (); RAHMAN2015198 ().

The strategies for handling missing data discussed so far are general and can be applied to both vectorial and time series data. On the other hand, methods that only apply to time series data include to impute missing values using smoothing and interpolation techniques such as the well-known last observation carried forward (LOCF) scheme that imputes the last non-missing value for the following missing values. Further approaches for MTS are linear interpolation, moving average and Kalman smoothing, to name a few ENGELS2003968 (). The list of possible imputation methods for MTS data is almost endless and a comprehensive overview of all these is beyond the scope of this paper. For a more detailed overview, we refer the interested reader to garcia2010pattern (); donders2006review (); durbin2012time ().

A drawback common to all imputation methods discussed above is that the information about which values are actually missing is lost. Moreover, imputation typically also introduces additional bias into the data due to strong assumptions made by the imputation method. For example mean imputation often leads to biased (shorter) estimates of the distances between data points than what they actually are hu2017strategies (). To resolve the bias problem it has been suggested to correct for the bias by introducing binary indicator variables to account for the missingness pattern hu2017strategies (). An additional problem with single imputation is that the uncertainty associated with the missing values is ignored since they are replaced with “certain” estimates, which in turn may lead to smaller estimated standard errors than the true standard errors. Multiple imputation white2011multiple () resolves this problem by estimating the missing values multiple times and thereby creating multiple complete datasets. Thereafter for example a classifier is applied to all datasets and the results are combined to obtain the final predictions. However, applying multiple imputation to MTS in a clustering setting is a non-trivial task that involves several challenges basagana2013framework ().

Despite that multiple imputation and other imputation methods can give satisfying results in some scenarios, these are ad-hoc solutions that lead to a multi-step procedure where the missing data are handled separately and independently from the rest of the analysis wells2013strategies (). This is not an optimal solution, and therefore several research efforts have been devoted over the last years to process incomplete temporal data without relying on imputation DBLP:journals/corr/ChePCSL16 (); bianchi2018time (); DBLP:journals/corr/abs-1711-06516 (); mikalsen2016learning (); pmlr-v56-Lipton16 (); Marlin:2012:UPD:2110363.2110408 (). In this regard, powerful kernel methods have been recently proposed, of which the TCK and LPS are prominent examples. Even though there are many similarities between these two kernels, the way missing data are dealt with is very different. In LPS the missing data handling abilities of decision trees are exploited. Along with ensemble methods, fuzzy approaches and support vector solutions, decision trees can be categorized as machine learning approaches for handling missing data garcia2010pattern (). Common to these approaches is that the missing data are handled naturally by the machine learning algorithm. One can also argue that the way missing data are dealt with in the TCK belongs to this category, since an ensemble approach is exploited. However, it can also be categorized as a likelihood-based approach since the underlying models in the ensemble are Gaussian mixture models. In the likelihood-based approaches the full, incomplete dataset is analyzed using maximum likelihood (or maximum a posteriori, equivalently), typically in combination with the expectation-maximization (EM) algorithm schafer2002missing (); little2014statistical ().

The main advantage of these methods, compared to imputation methods, is that the missing data are handled automatically and no additional tasks are left to the user. For example in multiple imputation, a careful selection of the imputation model and other variables is needed to do the imputation schafer2002missing (), which in particular in an unsupervised setting can turn out to be problematic. Moreover, similarly to multiple imputation, unbiased predictions are guaranteed if data are MAR.

A more detailed description of the TCK and LPS kernels is provided in the next subsection, along with a description of the other kernels for MTS used in this work.

### 2.3 Multivariate time series kernels

Kernel methods have been of great importance in machine learning for several decades and have applications in many different fields Jenssen2010 (); Jenssen2013 (); camps2009kernel (); soguero2016support (). Within the context of time series, a kernel is a similarity measure that also is positive semi-definite shawe2004kernel (). Once defined, such similarities between pairs of time series may be utilized in a wide range of applications, such as classification or clustering, benefiting from the vast body of work in the field of kernel methods.

#### Linear kernel

The simplest of all kernel functions is the linear kernel, which for two data points represented as vectors, and , is given by the inner product , possibly plus a constant . One can also apply a linear kernel to pairs of MTS once they are unfolded into vectors. However, by doing so the information that they are MTS and there might be inherent dependencies in time and between attributes, is then lost. Nevertheless, in some cases such a kernel can be efficient, especially if the MTS are short chen2013model (). If the MTS contain missing data, the linear kernel requires a preprocessing step involving e.g. imputation.

#### Global alignment kernel

The most widely used time series similarity measure is dynamic time warping (DTW) Berndt:1994:UDT:3000850.3000887 (), where the similarity is quantified as the alignment cost between the MTS. More specifically, in DTW the time dimension of one or both of the time series is warped to achieve a better alignment. Despite the success of DTW in many applications, similar to many other similarity measures, it is non-metric and therefore cannot non-trivially be used to design a positive semi-definite kernel marteau2015recursive (). Hence, it is not suited for kernel methods in its original formulation. Because of its popularity there have been attempts to design kernels exploiting the DTW. For example Cuturi et al. designed a DTW-based kernel using global alignments cuturi2007kernel (). An efficient version of the global alignment kernel (GAK) is provided in cuturi2011fast (). The latter has two hyperparameters, namely the kernel bandwidth and the triangular parameter. These are usually set using some heuristics. GAK does not naturally deal with missing data and incomplete datasets, and therefore also requires a preprocessing step involving imputation.

#### Time series cluster kernel

The TCK is based on an ensemble learning approach Strehl:2003 () wherein the robustness to hyperparameters is ensured by joining the clustering results of many Gaussian mixture models (GMM) to form the final kernel. Hence, no critical hyperparameters have to be tuned by the user, and the TCK can be learned in an unsupervised manner. To ensure robustness to sparsely sampled data, the GMMs that are the base models in the ensemble, are extended using informative prior distributions such that the missing data is explicitly dealt with.

More specifically, the TCK matrix is built by fitting GMMs to the set of MTS for a range of number of mixture components. The idea is that by generating partitions at different resolutions, one can capture both the local and global structure of the data. Moreover, to capture diversity in the data, randomness is injected by for each resolution (number of components) estimating the mixture parameters for a range of random initializations and randomly chosen hyperparameters. In addition, each GMM sees a random subset of attributes and segments in the MTS. The posterior distributions for each mixture component are then used to build the TCK matrix by taking the inner product between all pairs of posterior distributions. Finally, given an ensemble of GMMs, the TCK is created in an additive way by using the fact that the sum of kernels is also a kernel. In this work, we have modified the kernel slightly from the way it was originally proposed in mikalsen2017time () by normalizing the vectors of posteriors to have unit length in the -norm. This provides an additional regularization that may increase the generalization capability of the learned model. A more detailed description of the method is provided in A.

#### Learned pattern similarity

LPS is a similarity measure that satisfies the requirements of a kernel, as shown in mikalsen2017time (), which can naturally deal with MTS. Similar to the TCK, the LPS is also based on extracting random segments. Additionally, the LPS is similar to the TCK in the sense that one in an unsupervised way can learn a similarity between time series that is robust to hyperparameter choices and can deal with missing data using the missing data handling properties of tree-based learning. It generalizes the well-known autoregressive models shumway () to local autopatterns using multiple lag values for autocorrelation. These autopatterns are supposed to capture the local dependency structure in the time series and are learned using a tree-based learning strategy.

More specifically, a time series is represented as a matrix of segments. Randomness is injected to the learning process by randomly choosing time segment (column in the matrix) and lag for each tree in the random forest. A bag-of-words type compressed representation is created from the output of the leaf-nodes for each tree. The final time series representation is created by concatenating the representation obtained from the individual trees, which in turn are used to compute the similarity using a histogram intersection kernel barla2003histogram ().

Given two MTS and , a formal expression for the LPS-kernel is

(1) |

where is the th entry of the concatenated bag-of-words representation . More precisely, is a concatenation of -dimensional frequency vectors of instances in the terminal nodes from all trees.

### 2.4 Model development

The kernels that we described in the previous section are used to compute a kernel matrix on a training set, which is created by randomly selecting 80 percent of the dataset. The remaining 20 percent is used as test set. The LPS and TCK kernels are computed on the incomplete dataset containing missing data, whereas the GAK and linear kernel cannot work on incomplete datasets, and we therefore compute these on 6 different complete datasets obtained using mean imputation, LOCF imputation, 0-imputation, and replicates of these corrected for bias. In the bias corrected (BC) datasets we double the number of attributes in each MTS by stacking a binary MTS, representing imputed elements, to the imputed MTS. When using imputation of the mean, we calculate the mean for each attribute in the MTS, across all time intervals in the postoperative window and all patients in the training set. If an element is missing in the first time interval, we replace it with the mean when we do LOCF imputation.

After having computed the different kernels, we take an unsupervised approach to classifying the patients with SSI using the spectrum of the kernel matrices. We employ a variant of spectral clustering consisting of two steps, namely kPCA followed by k-means. In the first step, kPCA with the learned MTS kernel is used to compute a low dimensional representation of the MTS. Thereafter we cluster the learned representations using k-means. We assume that the number of clusters is known and set it to 2. Out-of-sample data are assigned to clusters according to the cluster labels of the k-nearest neighbors (kNN) in the training set. The processing pipeline we have described here is also illustrated in Fig. 1.

### 2.5 Model evaluation

The different models are evaluated both on the training and test set. Because of the imbalanced classes we decide to use F1-score hripcsak2005agreement () instead of accuracy as performance measure. F1-score is a function of two metrics, namely precision and recall. These two metrics are also commonly referred to as positive predictive value and sensitivity, respectively. Precision is the fraction of true positives (have infection) among all those that are classified (clustered) as positive cases, whereas recall is the fraction of positive cases in the gold standard classified as positive. F1-score can be expressed in terms of these two metrics as follows:

(2) |

In order to adapt this to an unsupervised regime, we define clustering F1-score similarly to how clustering accuracy is defined. i.e. we use the permutation of the labels provided by the clustering algorithm that gives the highest score. In the following, we refer to both clustering F1-score and classification F1-score simply as F1-score.

The procedure described in the previous section is repeated 10 times such that we can compute both mean and standard errors for the F1-score, i.e. we randomly select 10 different training and test sets (80 and 20 percent, respectively), and repeat the same process on all of them. In addition, in order to study how stable and robust the different methods are to varying length of the MTS, we vary the size of the postoperative window from 7 to 20 days.

## 3 Results

The TCK and LPS are run using default hyperparameters MikalsenTCK (); baydogan2016time (), with the exception for the LPS that we increase the minimal segment length from to percent of the length of the MTS to account for the short time series. In accordance with Cuturi (), for GAK we set the bandwidth to two times the median distance of all MTS in the training set scaled by the square root of the median length of all MTS, and the triangular parameter Cuturi () to 0.2 times the median length of all MTS. Distances are measured using the canonical metric induced by the Frobenius norm. In the linear kernel we set the constant to 0.

Fig. 2 shows mean F1-score over 10 runs on the training (left) and test set (right), and standard errors, obtained using LPS and TCK kernels, followed by kPCA to 10 dimensions and k-means, where test data are clustered using a kNN classifier with and the cluster assignments as labels. Initial experiments showed that the clustering results are stable to varying values of these hyperparameters. For easier comparison, in the same figure we have also added results obtained with the GAK and linear kernel on the imputed dataset that gives the highest F1-score, namely 0-imputation, whereas the results for all 6 complete datasets are shown in Fig. 3.

We observe that the two kernels, TCK and LPS, that can explicitly deal with the missing data, give very similar results both on the training and test set. When the postoperative window is 7 days the two methods yield an F1-score of approximately 0.63, then the performance increases almost linearly from day 7 to 15 where it stabilizes around a F1-score of 0.80. Further, it can be seen that TCK, LPS and GAK perform worse than the linear kernel when the postoperative window is short ( days). However, as the size of the postoperative window increases, the relative performance of the TCK and LPS with respect to both GAK and the linear kernel improve. Even if the differences between LPS and TCK are small, it can be seen that the latter yields slightly better performance (in average) for all sizes of the postoperative window both in training and test. Moreover, the standard errors with the TCK are smaller, in particular in training when the postoperative window is short. This is particularly important in unsupervised frameworks like the one we are proposing.

For the two kernels that work on the imputed data, we observe that GAK (0 and 0+BC) and Linear (0, 0+BC, LOCF+BC) perform quite similarly, and these are the five imputed data methods that give the best performance. A common pattern for these five methods, and especially for GAK (0 and 0+BC), is that the F1-score increases when the postoperative window is increased from 7 to around 11 days (at least on the training set), but then the F1-score slowly starts decrease after that. The performance of the seven other imputation methods is considerably worse.

To further investigate the differences between LPS and TCK, beside F1-score, in Fig. 4(a) and Fig. 4(b) we show the kPCA representation corresponding to the two largest eigenvalues obtained using these two kernels on the training set with postoperative window size equal to 20. Interestingly the representations created by the LPS and TCK are very different. The LPS has a clearer manifold structure, whereas the TCK is more spread out in the plane. Even though it is difficult to argue that one of these two representations is superior to the other, the TCK at least more clearly reflects that the cohort of patients is very diverse and complex because of large individual differences. To better understand the performance of the GAK and linear kernel we also plot the 2D kPCA representation obtained on 0-imputed data in Fig. 4(c) and Fig. 4(d). We note that these are heavily influenced by outliers. Apart from the outliers, the other datapoints become very compact and close to each other, and it is therefore not strange that the clustering algorithm does not identify the groups correctly.

As a baseline to compare the performance of our method, we follow the idea proposed by hu2017strategies () and compute higher level non-temporal features from the longitudinal data. Specifically, in addition to the average value of each of the blood tests, we compute two extreme values (the maximum and minimum) over the postoperative window. In this setting we consider a missing value as a specific blood test that is missing entirely during the postoperative window. This baseline cannot naturally deal with missing data and therefore we use the same six missing data imputation strategies as described for the temporal features. We apply a linear kernel to the manual features and then we follow the same scheme as in the temporal case. Fig. 5 shows mean F1-score over 10 runs on training (left) and test set (right) and standard errors obtained using this baseline on the six imputed datasets. Similarly to the results obtained with the linear kernel and GAK on the temporal data, the results obtained with this baseline heavily depend on the type of imputation method, and are also very fluctuating as function of length of the postoperative window.

In addition to comparing different kernels and methods for dealing with missing data, to see how well the unsupervised classification part works, we also benchmark the proposed unsupervised framework against a supervised baseline where the two first steps are the same, namely to compute the kernel and thereafter do kPCA. However, on the 10 dimensional kPCA representation we employ a kNN classifier with using the true labels. Figure 6 shows mean classification F1-score and standard errors over the same 10 randomly drawn training (left) and test sets (right). Also with this baseline the performance of the TCK and LPS kernel is very similar. Moreover, in general, the supervised baseline performs better on the training set than the unsupervised method. This is expected. However, on the test set, the F1-scores obtained using the supervised baseline is almost identical to the proposed method for all sizes of the postoperative window. This implies that labelling is not a required task for identifying patients with surgical site infection from blood samples.

## 4 Discussion

In this work, we have proposed a framework for SSI identification based on secondary use of EHR data. Our first objective was to alleviate the problem of getting access to labelled EHR data. The results presented in the previous section clearly show that the proposed framework can detect accurately SSI, without relying on supervised training. In fact, the supervised baseline did not improve the performance compared to the proposed method. The second objective has been to deal with missing data effectively. We tackled this problem by introducing robust time series kernels as a main component in our framework and in the following we discuss our findings.

For the two kernels that work on the imputed data, GAK and linear, the choice of imputation method heavily affects the performance more than the choice of kernel itself. Both the linear kernel and GAK give very bad results especially with mean imputation, but also with LOCF imputation. In those cases, the mean F1-score is far from smooth over time. On the other hand, 0-imputation gives good results compared to the other imputation methods for both kernels. This can be surprising, since 0-imputation introduces a strong bias because blood test values are positive and therefore the value 0 is very rarely a good estimate. However, a possible explanation might be that 0-imputation provides some auto-correction for the bias, since the 0s now describe the missingness pattern. This result is in accordance with what Hu et al. found in hu2017strategies () when taking a manual feature design approach to the problem.

Regarding bias correction, it is maybe not so surprising that it does not lead to improved performance on the 0-imputed dataset for any of the two kernels, since 0-imputation in itself seems to have the same effect (ref. previous paragraph). On the other hand, since the performance obtained using the linear kernel combined with mean or LOCF imputation is substantially improved, whereas nothing happens for the GAK kernel, it is difficult to draw conclusions about the effect of bias correction.

We also note that the best imputed data methods all share the same pattern; the performance peaks around postoperative window size 11. The natural behaviour should be that the F1-score should increase as more information is added when the length of the postoperative window is increased. However, when the opposite happens, this indicates that all these imputation methods introduce a bias into the data that affects the performance. Hence, due to the high variance in the performance across different imputation methods and length of postoperative window, and since is difficult to do model selection in the unsupervised setting of our framework, we conclude that the linear kernel and GAK are not suitable for the task under analysis.

Not surprisingly, the time series kernels, TCK, LPS and GAK perform worse than the linear kernel when the postoperative window is short ( days). When the MTS are shorter than 10 days the time dependency and time structure is probably not clear enough to be fully utilized by the specialized time series kernels. However, as the size of the postoperative window increases, the F1-score obtained using LPS and TCK increase smoothly, before it stabilizes around a F1-score of 0.80 when the length is around 14-15 days, considerably higher than the F1-score obtained using the GAK and linear kernel. This behaviour verifies the robustness of these two kernels with respect to varying length of postoperative window and missingness patterns.

We note that for the TCK and LPS kernels an underlying assumption is that values are MAR, whereas the missingness for EHR-data could be partly due to MNAR. However, it has been demonstrated by several authors that a slightly wrong assumption of MAR in many realistic scenarios do not have a big impact in terms of biased predictions collins2001comparison (); schafer2002missing (); hu2017strategies (). In our case we do not know by how much the assumption of MAR is broken. However, experiments in mikalsen2017time () demonstrated these kernels’ (and in particular TCK’s) robustness to large fractions of missing data – also in the case of MNAR data, whereas the imputation methods suffered in cases with much missing data. Even though, the data considered in this paper is completely different, it is interesting to observe a similar behaviour here.

### 4.1 Limitations and further work

Although the results obtained in this paper are promising just based on blood results, previous studies soguero2016predicting () have shown that the combination of heterogeneous data sources (e.g. free text, drugs, ICD-9 or vital signs) from the EHR might provide better performance. Including more data, however, comes at the cost of a more complicated and computationally demanding analysis. In addition, we did not differentiate between deep and superficial SSI in this work. These issues are subject to further work. Moreover, in this work, we focused on identifying patients with SSI based on postoperative data. In further work we would like to do prediction of SSI based on preoperative data and investigate if we can predict that the patient gets SSI before he or she is diagnosed with the complication, which is a framework that would be very valuable to operationalize in the clinic. The framework used in this work is general, and not restricted to SSI, it will therefore be interesting to also see how well it works on detecting other complications or diseases.

We could have benchmarked the proposed method against other kernels as well, but chose to not do that for the clarity of the presentation, and because the results obtained using GAK and the linear kernel indicate that results are more dependent on the imputation method than the type of kernel. We also tried with different clustering algorithms in the last step such as hierarchical clustering and kNN mode seeking ensemble clustering NordhaugMyhre2018491 (), but initial experiments did not improve performance, and we have therefore, for the conciseness of the presentation, not included them.

## 5 Conclusions

Hospital acquired infections in general, and surgical site infection in particular, are major problems at modern hospitals nowadays. To be able to reduce this problem, accurate prediction of SSI is of utmost importance. In this study, we showed that analyzing EHR data as MTS within a kernel framework can be very powerful in that respect. In particular, the LPS and TCK kernels that explicitly can deal with the missing data, turn out to be more robust and work better than those kernels that require the incomplete data to be pre-processed using some imputation method.

Moreover, because of the two kernels’ robustness to hyperparameters, we showed that we can completely unsupervised identify patients with SSI and perform similarly to a supervised baseline, hence alleviating the problem of a time consuming and expensive manual label annotation process, often unfeasible for large datasets. Worth mentioning in that respect, is that in this paper we have also illustrated the power of using only blood tests as the data source, hence also reducing the burden for the patients and data engineers.

## Conflict of interest

The authors have no conflict of interest related to this work.

## Acknowledgement

This work was partially funded by the Norwegian Research Council FRIPRO grant no. 239844 on developing the Next Generation Learning Machines. Cristina Soguero-Ruiz is partially supported by project TEC2016-75361-R from Spanish Government and by project DTS17/00158 from Institute of Health Carlos III (Spain).

The authors would like to thank Kristian Hindberg from UiT The Arctic University of Norway for his assistance on preprocessing and extracting the data from the EHR system. We would also like to thank Rolv-Ole Lindsetmo and Knut Magne Augestad from the University Hospital of North Norway, Fred Godtliebsen from UiT, together with Stein Olav Skrøvseth from the Norwegian Centre for E-health Research for helpful discussions throughout the study and manuscript preparation.

## Appendix A Tck

### Notation

The following notation is used. A multivariate time series (MTS) is defined as a (finite) sequence of univariate time series (UTS), where each attribute, , is a UTS of length . The number of UTS, , is the dimension of . The length of the UTS is also the length of the MTS . Hence, a –dimensional MTS, , of length can be represented as a matrix in . Given a dataset of MTS, we denote the -th MTS. An incompletely observed MTS is described by the pair , where is a binary MTS with entry if the realization is missing and if it is observed.

### DiagGMM

To build the TCK kernel matrix, we first fit different diagonal covariance GMM (DiagGMM) to the MTS dataset. In the DiagGMM one assumes time-dependent means, expressed by , where is a UTS, and a time-constant covariance matrix is , being the variance of attribute . Moreover, the data is assumed to be missing at random (MAR), i.e. the missing elements are only dependent on the observed values. Under these assumptions, missing data can be analytically integrated away, such that imputation is not needed rubin1976inference (), and the pdf for the incompletely observed MTS is given by

(3) |

The conditional probability of given , can be found using Bayes’ theorem,

(4) |

The parameters of the DiagGMM are learned using a maximum a posteriori expectation maximization algorithm, as described in mikalsen2017time ().

### Ensemble strategy

To ensure diversity, each GMM model uses a number of components from the interval , where is the maximal number of mixture components. For each number of components, we apply different random initial conditions and hyperparameters. We let be the index set keeping track of initial conditions and hyperparameters (), and the number of components (). Moreover, each model is trained on a random subset of MTS, accounting only a random subset of variables , with cardinality , over a randomly chosen time segment . The inner products of the posterior distributions from each mixture component are then added up to build the TCK kernel matrix, according to the ensemble strategy ensemble (). Algorithm 1 describes the details of the method.

### Method details

Algorithm 1 describes the details of the method. is the index set keeping track of initial conditions and hyperparameters (), and the number of components ().

In order to be able to compute similarities with MTS not available at the training phase, one needs to store the time segments , subsets of attributes , DiagGMM parameters and posteriors . Then, the TCK for such out-of-sample MTS is evaluated according to Algorithm 2.

## References

### References

- S. S. Lewis, R. W. Moehring, L. F. Chen, D. J. Sexton, D. J. Anderson, Assessing the relative burden of hospital-acquired infections in a network of community hospitals, Infection Control & Hospital Epidemiology 34 (11) (2013) 1229–1230.
- S. S. Magill, W. Hellinger, J. Cohen, R. Kay, C. Bailey, B. Boland, D. Carey, J. d. Guzman, K. Dominguez, J. Edwards, et al., Prevalence of healthcare-associated infections in acute care hospitals in Jacksonville, Florida, Infection Control 33 (03) (2012) 283–291.
- G. de Lissovoy, K. Fraeman, V. Hutchins, D. Murphy, D. Song, B. B. Vaughn, Surgical site infection: incidence and impact on hospital utilization and treatment costs, American Journal of Infection Control 37 (5) (2009) 387–397.
- C. Y. Ko, B. L. Hall, A. J. Hart, M. E. Cohen, D. B. Hoyt, et al., The American college of surgeons national surgical quality improvement program: achieving better and safer surgery, The Joint Commission Journal on Quality and Patient Safety 41 (5) (2015) 199–AP1.
- E. H. Lawson, B. L. Hall, C. Y. Ko, Risk factors for superficial vs deep/organ-space surgical site infections: implications for quality improvement initiatives, JAMA surgery 148 (9) (2013) 849–858.
- J. Blumetti, M. Luu, G. Sarosi, K. Hartless, J. McFarlin, B. Parker, S. Dineen, S. Huerta, M. Asolati, E. Varela, et al., Surgical site infections after colorectal surgery: do risk factors vary depending on the type of infection considered?, Surgery 142 (5) (2007) 704–711.
- J. D. Whitehouse, N. D. Friedman, K. B. Kirkland, W. J. Richardson, D. J. Sexton, The impact of surgical-site infections following orthopedic surgery at a community hospital and a university hospital adverse quality of life, excess length of stay, and extra cost, Infection Control & Hospital Epidemiology 23 (04) (2002) 183–189.
- R. Shah, E. Pavey, M. Ju, R. Merkow, R. Rajaram, M. W. Wandling, M. E. Cohen, A. Dahlke, A. Yang, K. Bilimoria, Evaluation of readmissions due to surgical site infections: A potential target for quality improvement, The American Journal of Surgery.
- P. L. Owens, M. L. Barrett, S. Raetzman, M. Maggard-Gibbons, C. A. Steiner, Surgical site infections following ambulatory surgery procedures, JAMA 311 (7) (2014) 709–716.
- J. Silvestre, J. Rebanda, C. Lourenço, P. Póvoa, Diagnostic accuracy of C-reactive protein and procalcitonin in the early detection of infection after elective colorectal surgery–a pilot study, BMC infectious diseases 14 (1) (2014) 444.
- C. Soguero-Ruiz, W. M. Fei, R. Jenssen, K. M. Augestad, J.-L. R. Álvarez, I. M. Jiménez, R.-O. Lindsetmo, S. O. Skrøvseth, Data-driven temporal prediction of surgical site infection, in: AMIA Annual Symposium Proceedings, Vol. 2015, American Medical Informatics Association, 2015, p. 1164.
- F. J. Medina-Fernández, D. J. Garcilazo-Arismendi, R. García-Martín, L. Rodríguez-Ortiz, J. Gómez-Barbadillo, J. M. Gallardo-Valverde, J. L. Martínez-Dueñas, E. Navarro-Rodríguez, E. Torres-Tordera, C. A. Díaz-López, et al., Validation in colorectal procedures of a useful novel approach for the use of C-reactive protein in postoperative infectious complications, Colorectal Disease 18 (3) (2016) O111–O118.
- M. R. Angiolini, F. Gavazzi, C. Ridolfi, M. Moro, P. Morelli, M. Montorsi, A. Zerbi, Role of C-reactive protein assessment as early predictor of surgical site infections development after pancreaticoduodenectomy, Digestive surgery 33 (4) (2016) 267–275.
- S. L. Gans, J. J. Atema, S. Van Dieren, B. G. Koerkamp, M. A. Boermeester, Diagnostic value of C-reactive protein to rule out infectious complications after major abdominal surgery: a systematic review and meta-analysis, International journal of colorectal disease 30 (7) (2015) 861–873.
- C. Soguero-Ruiz, K. Hindberg, J. L. Rojo-Álvarez, S. O. Skrøvseth, F. Godtliebsen, K. Mortensen, A. Revhaug, R.-O. Lindsetmo, K. M. Augestad, R. Jenssen, Support vector feature selection for early detection of anastomosis leakage from bag-of-words in electronic health records, IEEE journal of biomedical and health informatics 20 (5) (2016) 1404–1415.
- Z. Hu, G. B. Melton, E. G. Arsoniadis, Y. Wang, M. R. Kwaan, G. J. Simon, Strategies for handling missing clinical data for automated surgical site infection detection from the electronic health record, Journal of Biomedical Informatics 68 (2017) 112–120.
- P. C. Sanger, G. H. van Ramshorst, E. Mercan, S. Huang, A. L. Hartzler, C. A. Armstrong, R. J. Lordon, W. B. Lober, H. L. Evans, A prognostic model of surgical site infection using daily clinical wound assessment, Journal of the American College of Surgeons 223 (2) (2016) 259 – 270.e2.
- B. J. Wells, K. M. Chagin, A. S. Nowacki, M. W. Kattan, Strategies for handling missing data in electronic health record derived data, eGEMs 1 (3).
- K. Ø. Mikalsen, F. M. Bianchi, C. Soguero-Ruiz, R. Jenssen, Time series cluster kernel for learning similarities between multivariate time series with missing data, Pattern Recognition 76 (2018) 569–581.
- T. Mitsa, Temporal data mining, Chapman & Hall/CRC, 2010.
- K. Ø. Mikalsen, C. Soguero-Ruiz, K. Jensen, K. Hindberg, M. Gran, A. Revhaug, R.-O. Lindsetmo, S. O. Skrøvseth, F. Godtliebsen, R. Jenssen, Using anchors from free text in electronic health records to diagnose postoperative delirium, Computer Methods and Programs in Biomedicine 152 (Supplement C) (2017) 105 – 114.
- O. Chapelle, B. Schlkopf, A. Zien, Semi-Supervised Learning, 1st Edition, The MIT Press, 2010.
- Y. Halpern, S. Horng, Y. Choi, D. Sontag, Electronic medical record phenotyping using the anchor and learn framework, Journal of the American Medical Informatics Association.
- M. G. Baydogan, G. Runger, Time series representation and similarity based on local autopatterns, Data Mining and Knowledge Discovery 30 (2) (2016) 476–509.
- A. Y. Ng, M. I. Jordan, Y. Weiss, On spectral clustering: Analysis and an algorithm, in: Advances in neural information processing systems, 2002, pp. 849–856.
- F. Krzakala, C. Moore, E. Mossel, J. Neeman, A. Sly, L. Zdeborová, P. Zhang, Spectral redemption in clustering sparse networks, Proceedings of the National Academy of Sciences 110 (52) (2013) 20935–20940.
- S. Løkse, F. M. Bianchi, A.-B. Salberg, R. Jenssen, Spectral clustering using PCKID–a probabilistic cluster kernel for incomplete data, in: Scandinavian Conference on Image Analysis, Springer, 2017, pp. 431–442.
- K. Jensen, C. Soguero-Ruiz, K. Ø. Mikalsen, R.-O. Lindsetmo, I. Kouskoumvekaki, M. Girolami, S. O. Skrovseth, K. M. Augestad, Analysis of free text in electronic health records for identification of cancer patient trajectories, Scientific Reports 7 (2017) 46226.
- C. Soguero-Ruiz, K. Hindberg, I. Mora-Jiménez, J. L. Rojo-Álvarez, S. O. Skrøvseth, F. Godtliebsen, K. Mortensen, A. Revhaug, R.-O. Lindsetmo, K. M. Augestad, et al., Predicting colorectal surgical complications using heterogeneous clinical data and kernel methods, Journal of Biomedical Informatics 61 (2016) 87–96.
- E. H. Lawson, C. Y. Ko, J. L. Adams, W. B. Chow, B. L. Hall, Reliability of evaluating hospital quality by colorectal surgical site infection type, Annals of surgery 258 (6) (2013) 994–1000.
- E. Limón, E. Shaw, J. Badia, M. Piriz, R. Escofet, F. Gudiol, M. Pujol, Post-discharge surgical site infections after uncomplicated elective colorectal surgery: impact and risk factors. the experience of the VINCat program, Journal of Hospital Infection 86 (2) (2014) 127 – 132.
- C. Gibbons, J. Bruce, J. Carpenter, A. Wilson, J. Wilson, A. Pearson, D. Lamping, Z. Krukowski, B. Reeves, Identification of risk factors by systematic review and development of risk-adjusted models for surgical site infection, Health Technol Assess 15 (2011) 1–156.
- T. C. Horan, R. P. Gaynes, W. J. Martone, W. R. Jarvis, T. G. Emori, CDC definitions of nosocomial surgical site infections, 1992: A modification of CDC definitions of surgical wound infections, Infection Control Hospital Epidemiology 13 (10) (1992) 606–608.
- R. L. Berger, L. T. Li, S. C. Hicks, J. A. Davila, L. S. Kao, M. K. Liang, Development and validation of a risk-stratification score for surgical site occurrence and surgical site infection after open ventral hernia repair, Journal of the American College of Surgeons 217 (6) (2013) 974 – 982.
- R. J. Little, D. B. Rubin, Statistical analysis with missing data, John Wiley & Sons, 2014.
- D. B. Rubin, Inference and missing data, Biometrika 63 (3) (1976) 581–592.
- J. L. Schafer, Analysis of incomplete multivariate data, CRC press, 1997.
- P. J. García-Laencina, J.-L. Sancho-Gómez, A. R. Figueiras-Vidal, Pattern classification with missing data: a review, Neural Computing and Applications 19 (2) (2010) 263–282.
- S. A. Rahman, Y. Huang, J. Claassen, N. Heintzman, S. Kleinberg, Combining Fourier and lagged k-nearest neighbor imputation for biomedical time series data, Journal of Biomedical Informatics 58 (2015) 198 – 207.
- J. M. Engels, P. Diehr, Imputation of missing longitudinal data: a comparison of methods, Journal of Clinical Epidemiology 56 (10) (2003) 968 – 976.
- A. R. T. Donders, G. J. van der Heijden, T. Stijnen, K. G. Moons, Review: a gentle introduction to imputation of missing values, Journal of clinical epidemiology 59 (10) (2006) 1087–1091.
- J. Durbin, S. J. Koopman, Time series analysis by state space methods, Vol. 38, OUP Oxford, 2012.
- I. R. White, P. Royston, A. M. Wood, Multiple imputation using chained equations: issues and guidance for practice, Statistics in medicine 30 (4) (2011) 377–399.
- X. Basagaña, J. Barrera-Gómez, M. Benet, J. M. Antó, J. Garcia-Aymerich, A framework for multiple imputation in cluster analysis, American journal of epidemiology 177 (7) (2013) 718–725.
- Z. Che, S. Purushotham, K. Cho, D. Sontag, Y. Liu, Recurrent neural networks for multivariate time series with missing values, CoRR abs/1606.01865. arXiv:1606.01865.
- F. M. Bianchi, L. Livi, A. Ferrante, J. Milosevic, M. Malek, Time series kernel similarities for predicting paroxysmal atrial fibrillation from ECGs, arXiv preprint arXiv:1801.06845.
- A. S. Strauman, F. M. Bianchi, K. Ø. Mikalsen, M. Kampffmeyer, C. Soguero-Ruiz, R. Jenssen, Classification of postoperative surgical site infections from blood measurements with missing data using recurrent neural networks, CoRR abs/1711.06516. arXiv:1711.06516.
- K. Ø. Mikalsen, F. M. Bianchi, C. Soguero-Ruiz, S. O. Skrøvseth, R.-O. Lindsetmo, A. Revhaug, R. Jenssen, Learning similarities between irregularly sampled short multivariate time series from EHRs, 3rd ICPR International Workshop on Pattern Recognition for Healthcare Analytics, Cancun, Mexico, 2016.
- Z. C. Lipton, D. Kale, R. Wetzel, Directly modeling missing data in sequences with RNNs: Improved classification of clinical time series, in: Proceedings of the 1st Machine Learning for Healthcare Conference, Vol. 56 of Proceedings of Machine Learning Research, PMLR, Children’s Hospital LA, Los Angeles, CA, USA, 2016, pp. 253–270.
- B. M. Marlin, D. C. Kale, R. G. Khemani, R. C. Wetzel, Unsupervised pattern discovery in electronic health care data using probabilistic clustering models, in: Proc. of 2nd ACM SIGHIT Int. Health Informatics Symposium, 2012, pp. 389–398.
- J. L. Schafer, J. W. Graham, Missing data: our view of the state of the art, Psychological methods 7 (2) (2002) 147.
- R. Jenssen, Kernel entropy component analysis, IEEE Trans Pattern Anal Mach Intell 33 (5) (2010) 847–860.
- R. Jenssen, Entropy-relevant dimensions in the kernel feature space: Cluster-capturing dimensionality reduction, IEEE Signal Processing Magazine 30 (4) (2013) 30–39.
- G. Camps-Valls, L. Bruzzone, Kernel methods for remote sensing data analysis, John Wiley & Sons, 2009.
- J. Shawe-Taylor, N. Cristianini, Kernel methods for pattern analysis, Cambridge university press, 2004.
- H. Chen, F. Tang, P. Tino, X. Yao, Model-based kernel for efficient time series analysis, in: Proceedings of the 19th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, ACM, 2013, pp. 392–400.
- D. J. Berndt, J. Clifford, Using dynamic time warping to find patterns in time series, in: Proceedings of the 3rd International Conference on Knowledge Discovery and Data Mining, AAAI Press, 1994, pp. 359–370.
- P.-F. Marteau, S. Gibet, On recursive edit distance kernels with application to time series classification, IEEE Transactions on Neural Networks and Learning Systems 26 (6) (2015) 1121–1133.
- M. Cuturi, J.-P. Vert, O. Birkenes, T. Matsui, A kernel for time series based on global alignments, in: Acoustics, Speech and Signal Processing, 2007. ICASSP 2007. IEEE International Conference on, Vol. 2, IEEE, 2007, pp. II–413.
- M. Cuturi, Fast global alignment kernels, in: Proceedings of the 28th International Conference on Machine Learning, 2011, pp. 929–936.
- A. Strehl, J. Ghosh, Cluster ensembles – a knowledge reuse framework for combining multiple partitions, Journal of Machine Learning Research 3 (2003) 583–617.
- R. H. Shumway, D. S. Stoffer, Time series analysis and its applications: with R examples, Springer Science & Business Media, 2010.
- A. Barla, F. Odone, A. Verri, Histogram intersection kernel for image classification, in: Proceedings of International Conference on Image Processing, Vol. 3, IEEE, 2003, pp. III–513.
- G. Hripcsak, A. S. Rothschild, Agreement, the f-measure, and reliability in information retrieval, Journal of the American Medical Informatics Association 12 (3) (2005) 296–298.
- K. Ø. Mikalsen, Time series cluster kernel (TCK) Matlab implementation, http://site.uit.no/ml (2017).
- Fast global alignment kernel Matlab implementation, http://www.marcocuturi.net/GA.html, accessed: 2017-12-02.
- L. M. Collins, J. L. Schafer, C.-M. Kam, A comparison of inclusive and restrictive strategies in modern missing data procedures., Psychological methods 6 (4) (2001) 330.
- J. N. Myhre, K. Ø. Mikalsen, S. Løkse, R. Jenssen, Robust clustering using a kNN mode seeking ensemble, Pattern Recognition 76 (2018) 491 – 505.
- A. L. Fred, A. K. Jain, Combining multiple clusterings using evidence accumulation, IEEE Transactions on Pattern Analysis and Machine Intelligence 27 (6) (2005) 835–850.