Learning Multimorbidity Patterns from Electronic Health Records Using Non-negative Matrix Factorisation

TL;DR

This paper introduces a novel non-negative matrix factorisation approach for analyzing the evolution of multimorbidity patterns over time using large-scale electronic health records, providing new insights and metrics.

Contribution

The study presents a new temporal phenotyping method using NMF, along with evaluation metrics and insights into disease network evolution from extensive EHR data.

Findings

Effective identification of disease clusters and their trajectories over time

Quantitative metrics for evaluating multimorbidity patterns

Insights into the emergence and evolution of multimorbidity networks

Abstract

Multimorbidity, or the presence of several medical conditions in the same individual, has been increasing in the population, both in absolute and relative terms. However, multimorbidity remains poorly understood, and the evidence from existing research to describe its burden, determinants and consequences has been limited. Previous studies attempting to understand multimorbidity patterns are often cross-sectional and do not explicitly account for multimorbidity patterns' evolution over time; some of them are based on small datasets and/or use arbitrary and narrow age ranges; and those that employed advanced models, usually lack appropriate benchmarking and validations. In this study, we (1) introduce a novel approach for using Non-negative Matrix Factorisation (NMF) for temporal phenotyping (i.e., simultaneously mining disease clusters and their trajectories); (2) provide quantitative metrics for the evaluation of disease clusters from such studies; and (3) demonstrate how the temporal characteristics of the disease clusters that result from our model can help mine multimorbidity networks and generate new hypotheses for the emergence of various multimorbidity patterns over time. We trained and evaluated our models on one of the world's largest electronic health records (EHR), with 7 million patients, from which over 2 million where relevant to this study.

Tables

Table 1. Table 1 : Glossary of the terms frequently used in this paper, when describing the methodology.

Term	Definition
Disease cluster (DC), or multimorbidity clusters (MC)	Each DC (or MC) is a vector of which entries are weights that are assigned to each and every disease to denote the extent to which they belong to it. Note that, as implied, each disease can belong to more than one DC.
DC time course	The strength by which a DC is expressed in a given patient over time. For instance, if the time course shows a larger value at age 50 than 40, it means that its corresponding DC has a stronger expression/presence at age 50 than 40.
Multimorbidity pattern (MP) and temporal phenotype	In our study, each MP or temporal phenotype consists of a DC and its corresponding time course. For instance, our NMF approach, decomposes $𝐃$ to a matrix of DCs (i.e., $𝐁$) plus a matrix of their time courses (i.e., $𝐀$). An MP might also be referred to as component (e.g., in PCA and ICA) or factor (e.g., in NMF)

Equations

$$\mathbf{D}\approx \mathbf{A}\times \mathbf{B}\mid \mathbf{A}\ge 0,\mathbf{B}\ge 0.$$

$$\begin{array}{c}{z}_1={\sum }_{NT}I(\mathbf{v}=1,\mathbf{w}=1)\\ z_2={\sum }_{NT}I(\mathbf{v}=1,\mathbf{w}=0)\\ z_3={\sum }_{NT}I(\mathbf{v}=0,\mathbf{w}=1)\\ z_4={\sum }_{NT}I(\mathbf{v}=0,\mathbf{w}=0),\end{array}$$

$$p(z|\theta )\propto \prod _{i=1}^4{\theta }_i^{z_i},$$

$$\begin{array}{c}{\theta }_1=P(\mathbf{v}=1,\mathbf{w}=1)\\ {\theta }_2=P(\mathbf{v}=1,\mathbf{w}=0)\\ {\theta }_3=P(\mathbf{v}=0,\mathbf{w}=1)\\ {\theta }_4=P(\mathbf{v}=0,\mathbf{w}=0).\end{array}$$

$$\kappa =\frac{{\theta }_1-E}{D(max({\theta }_1)-E)+(1-D)(E-min({\theta }_1)))}$$

$$D=\{\begin{array}{c}\frac{{\theta }_1-E}{2(max({\theta }_1)-E)}+0.5\text{, if }{\theta }_1\ge E\\ 0.5-\frac{{\theta }_1-E}{2(E-min({\theta }_1))}\text{, otherwise}.\end{array}$$

$${\tau }_{\mathbf{v},\mathbf{w}}=\{\begin{array}{c}1-\frac{{\theta }_1+{\theta }_3}{{\theta }_1+{\theta }_2)}\text{, if }{\theta }_2\ge {\theta }_3\\ \frac{{\theta }_1+{\theta }_2}{{\theta }_1+{\theta }_3}-1\text{, otherwise}\end{array}$$

Full text

††journal: arXiv

1 Introduction

Multimorbidity is generally defined as the presence of two or more chronic conditions in an individual [1]. There is growing evidence that the number of people with multimorbidity has been increasing in many populations, both in relative and absolute terms. This increasing burden has been attributed to a number of factors, including the trend towards an ageing population demographics, as well as factors relating to changes in lifestyle, health-seeking behaviour and the environment [2]. Research in this area has been growing, but most investigations have focused on predicting, preventing and managing disorders in isolation from one another. Therefore, more research is needed for a better understanding of this emerging burden, its underlying patterns and mechanisms, in order to anticipate its consequences for the health services and the provision of appropriate care [3].

Different types of studies in the past have tried several methods to investigate multimorbidity. We refer to the first group of such studies as “pairwise methods”; in these methods, disease pairs that show co-occurrence frequencies that are different from what their individual frequencies in the population would predict, are considered to be “connected” [4, 5, 6, 7]. Treating diseases as nodes and connectedness as edges, these studies formed networks of which the properties were then used to characterise various multimorbidity pathways. Liu et al. [8] took the idea further and formed a network where nodes were medications, tests and diagnoses. While pairwise methods are valuable in generating comorbidity hypotheses for disease pairs, their inability to address conditional independence (i.e., where the correspondence between disease $c_{i}$ and $c_{j}$ is due to a disease $c_{k}$ they both are linked to, i.e., $P(c_{i}|c_{j},c_{k})=P(c_{i}|c_{k})$, or $c_{i}\perp\!\!\!\perp c_{j}~{}|~{}c_{k}$ [9]) can make the multi-disease networks resulting from them misleading.

This has led to the rise of alternative methods for disease phenotyping and the study of multimorbidity, that can deal with multiple diseases simultaneously. Except a limited number of studies that used “probabilistic methods” such as latent class growth analysis [7] and Hidden Markov Model [10], majority of recent studies in the field have relied on “factorisation methods”. The earlier factorisation methods in this field started by forming a matrix, where the entry $i,j$ denotes a metric related to disease $i$ (or other concepts in EHR, such as medication and/or clinical measurements) in patient $j$; factor methods decompose such a matrix into $P$ multimorbidity patterns (MPs), each consisting of a disease cluster (DC) and the expression of that DC across patients. While relatively effective, such matrix factorisation approaches did not take into account the temporal aspect of MPs (i.e., disease-based temporal phenotyping) and have been mostly limited to mining static disease (and sometimes medication) clusters [11, 12, 13, 14, 15]. As an exception, Zhou et al. [16] introduced a matrix factorisation method that considers the temporal patterns in EHR data; while the method can be employed for multimorbidity analyses, the study’s primary focus was on the prediction of future diseases (as opposed to multimorbidity).

On the other hand, solutions for joint phenotyping (based on diseases plus other concepts, for instance) have been achieved through the use of tensor factorisation. For instance, both [17] and [18], used non-negative tensor factorisation for joint phenotyping of diagnoses and medications, through adding medications as a new dimension to the input matrix (and hence the tensor). Following a similar approach Perros et al. [19], introduced a time dimension and used tensor factorisation for joint phenotyping of diseases and time (i.e., disease-based temporal phenotyping); instead of explicitly dealing with time, however, their model deals with the chronological order/index of encounters. In another recent study of temporal phenotyping, Zhao et al. [20] used tensor factorisation and showed the effectiveness of their resulting phenotypes for stratifying the risk of cardiovascular diseases. In summary, previous temporal phenotyping studies relied on relatively complex tensor-based factorisation approaches, which makes the use of simpler matrix factorisation techniques for this goal under-explored.

In this study, we introduce a novel design that can enable matrix factorisation techniques, such as NMF [21, 22], for disease-based temporal phenotyping (i.e., to mine the DCs and their expression over time in patients). This was facilitated by an assumption that is key to our study: the underlying DCs are the same for all patients, but their expression pattern (i.e., the strength of a DC’s expression in a given year of one’s life) varies from one person to the other. We trained our approach on one of the world’s largest EHR datasets – consisting of more than 7 million individuals (2 million were appropriate for this study) – and evaluated its results for the study of multimorbidity. Given that the past studies of multimorbidity are often cross-sectional, use arbitrary narrow age ranges, and lack appropriate benchmarking and validations (of their methodology and results), this paper can introduce a new methodology pipeline to the field of multimorbidity research that can employ the temporal information – hidden in longitudinal datasets such as EHR – to generate and evaluate new hypotheses on how diseases occur over time.

2 Materials and methods

In this section, we explain the source of our EHR data, our approach in using NMF, and the network analyses that we carried out to show the strength of temporal disease-based phenotyping for the study of MPs. Note that, for the rest of the paper, matrices will be denoted by upper case bold fonts (e.g., $\mathbf{A}$), vectors will be denoted by lower case bold fonts (e.g., $\mathbf{a}$), and everything else (scalar and indices) will be denoted with no bolding of the fonts.

2.1 EHR Data

An important factor that can help address the aforementioned issues in the current state of multimorbidity studies, is the rapid growth in the development of healthcare information systems, and the growing interest in utilising EHR. Particularly, the longitudinal nature of EHR and its richness (e.g., containing diagnoses, medications, and tests/measurement) can provide a unique opportunity to study temporal multi-modal phenotyping.

In this study, we used the Clinical Practice Research Datalink (CPRD) [23]; its longitudinal primary care data from a network of 674 general practices (GP) in the UK, and is linked to secondary care (i.e., Hospital Episode Statistics, HES) and other health and administrative databases (e.g., Office for National Statistics’ death registration). One in ten GPs in the UK contributes data to CPRD; it covers 35 million patients, among whom nearly 7 million currently registered patients. CPRD is broadly representative of the population by age, sex, and ethnicity [24]. It has been extensively validated and is considered as the most comprehensive longitudinal primary care database [25], with several large-scale epidemiological reports [26, 27, 28] adding to its credibility.

HES, on the other hand, contains data on hospitalisations, outpatient visits, accident and emergency for all admissions to National Health Service (NHS) hospitals in England [29]. Approximately 75% of the CPRD GPs in England (58% of all UK CPRD GPs) participate in patient-level record linkage with HES, which is performed by the Health and Social Care Information Centre [30]. In this study, we only considered the data from GPs that consented to (and hence have) record linkage with HES. The importance of primary care at the centre of the National Health System in the UK, the additional linkages, and all the aforementioned properties, make CPRD one of the most suitable EHR datasets in the world for data-driven clinical/medical discovery and machine learning.

2.2 Study Population

In order to have a comprehensive coverage of patients’ health journey, we only considered patients with at least 5 years of follow-up; this resulted in a total number of 2,204,178 patients with a total number of 25,791,493 clinical encounters. Note that, in this study, we are interested in “incident cases” and not in the “prevalent cases” of diseases, i.e. new occurrences of diseases rather than diseases carried over. Therefore, we only considered the first occurrence of each disease happening after 1 year of any patient’s registration date with the general practice clinic (as the first year after registration is likely to contain diseases carried over rather than new occurrences of diseases).

Another important step in processing CPRD was to create consistent disease classifications between GP and HES data and choosing the appropriate level of granularity in diseases’ hierarchy . In HES, diseases are coded using ICD-10 (International Classification of Diseases [31]), whereas in the GP records diseases are coded using Read Code [32]. ICD-10 offers a hierarchical form that makes its use for data mining and machine learning much more convenient. Therefore, we decided to map the diseases to the ICD-10 domain, i.e., we mapped the Read Codes to ICD-10 codes using the mapping provided by NHS Digital [33]. When no direct mapping was available, Read Codes were first mapped to SNOMED-CT codes [34] (also provided by NHS Digital [35]), and the latter are then mapped to ICD-10 codes using the mapping provided by the US National Library of Medicine [36].

As mentioned earlier, ICD-10 codes are organised in a tree-like hierarchy. Working at the highest level of the ICD-10 hierarchy will result in only a few diseases and hence is not likely to help unravel the complex underlying dependencies among diseases. Conversely, working at the lowest level of the hierarchy will generate thousands of diseases, each with small number of occurrences (and even a smaller number of co-occurrences) and hence an overall difficulty of mining useful inter-disease patterns. In this study, similarly to the work in [15], we chose to work at the ICD-10 block level, which provides a good trade off between granularity and co-occurrence, and can lead to medically interpretable results. Furthermore, for the purpose of this study, we eliminate all the diagnoses relating to pregnancy, general symptoms, external causes and administration (i.e., ICD-10 chapters XV, XVI, XVIII, XIX, XX and XXI); similar to [6]. The resulting 142 ICD-10 blocks are what we will refer to as diseases from here on – more details on our inclusion/exclusion process is shown in Figure 1.

At the end of this process, our cut of CPRD dataset includes nearly 2.2M adult patients (aged 16 years and over) and 65M events; from here, this will be referred to as the data in this study. Figure 2 shows some of the key characteristics of the data. When compared to the past studies, in terms of number of patients, length of follow-up and being representative of the population, our data is comparable to, if not better than, the best-case studied so far.

2.3 Non-negative Matrix Factorisation

Non-negative matrix factorisation (NMF, or NNMF) refers to a group of algorithms that decompose a matrix $\mathbf{D}$ into (usually) two matrices $\mathbf{A}$ and $\mathbf{B}$, with the property that all three matrices have no negative elements, i.e.,

[TABLE]

This non-negativity makes the resulting matrices easier to inspect and interpret. Also, in applications such as processing of count data (i.e., the starting point of our multimorbidity analysis), non-negativity is inherent to the data being considered. Since the NMF problem does not have an exact analytical solution in general, there has been a range of numerical approximations for it [37, 38, 39, 40, 41, 42, 43, 44, 45, 46]. In this paper, we use the Kullback-Leibler divergence and simple multiplicative updates [37, 39], enhanced to avoid numerical underflow [41] as implemented by Nimfa python package [44].

2.4 Modelling Pipeline

We start our modelling pipeline by forming a “disease matrix” $\mathbf{D}_{p}$ for each patient $p$, where $\mathbf{D}_{p}(i,j)=1$ if patient $p$ had the first incidence of disease $i$ at age (in years) $j$; $\mathbf{D}_{p}(i,j)=0$ otherwise. This makes $\mathbf{D}_{p}$ a $T\times C$ matrix, where $C$ is the number of conditions/diseases (i.e., 142 in our case) and $T$ is the maximum age we track a patient for (i.e., 114 years in our case). Denoting the total number patients by $N$, this process will result in $N$ such $\mathbf{D}_{p}$ matrices. We carry out two processing steps to create a final disease matrix $\mathbf{D}$ for the NMF analysis.

Given the variability in disease prevalences (i.e., some diseases are more common than others), the counts corresponding to rare diseases (such as tuberculosis) are expected to be much lower than the counts corresponding to more common diseases (such as respiratory infections). Therefore, when solving the NMF model, the results can be biased towards explaining the more frequent diseases (i.e., the higher counts). In order to correct for this, in our first processing step, we used an adjustment inspired by TF-IDF (term frequency-inverse document frequency), which is commonly used in natural language processing and information retrieval [47]. More specifically, we introduce DF-IPF (disease frequency - inverse patient frequency), which equals $DF\times IPF$, where we set $DF=1$ as we only considered the first occurrence of each disease. Then for each disease $i$, we defined $IPF(i)=\log(N/N_{i})$, where $N$ is the total number of patients and $N_{i}$ is the number of patients who had disease $i$. The DF-IPF adjustment for $\mathbf{D}_{p}$ will simply result from the multiplication of its entries with the appropriate inverse patient frequency, i.e., $\mathbf{D}_{p}(i,j)\leftarrow\mathbf{D}_{p}(i,j)*IPF(i)$.

As each patient will only have a relatively small number of diseases, $\mathbf{D}_{p}$ is expected to be sparse. On the other hand, NMF does not explicitly model age as a temporal concept (i.e., no explicit model for the relationship among rows in $\mathbf{D}$). In medicine, however, one will not see a meaningful difference between a disease happening at age $a$ or $a\pm$1 or 2 years; in factorisation of a matrix like $\mathbf{D}$ using standard NMF, this property will not exist and the two scenarios will not be seen as similar. Furthermore, we know from the practice of medicine that the date at which certain chronic disease gets recorded is a noisy concept: For instance, one’s hypertension diagnosis can be delayed by months or years due to not noticing or ignoring the symptoms, and/or delaying a doctor visit; or, diseases such dementia are known to have long preclinical periods, where patients who visit their doctors less regularly can, on average, be more likely to have their diagnosis delayed. Therefore, in a second adjustment of $\mathbf{D}_{p}$, and in order to address NMF’s lack of temporal regularisation, and take into account the noise and/or uncertainty we have around the age that a particular disease is recorded at, we smooth each column of $\mathbf{D}_{p}$, using a Gaussian filter $g_{\sigma}$; the optimal value of $\sigma$ will be determined empirically.

As mentioned earlier, to the best of authors’ knowledge, the only MP analyses that took the time into account used tensor factorisation methods. Therefore, one of the key contributions of this paper is to enable matrix factorisation techniques to result in temporal phenotyping. Assuming that DCs are the same across all patients – and what varies from patient to patient is the strength of these DCs expressing themselves at different ages – enables us to use NMF in a unique way that enables temporal phenotyping. That is, before we carry out the NMF analysis, we perform one last step: Concatenating the $\mathbf{D}_{p}$ matrices across the age dimension to form $\mathbf{D}$ (i.e., $\mathbf{D}=[\mathbf{D}_{1},\mathbf{D}_{2},...,\mathbf{D}_{N}]$); a $(T*N)\times C$ matrix. Decomposing $\mathbf{D}$ to $Q$ DCs (or factors), makes $\mathbf{A}$ a $(T*N)\times Q$ matrix and $\mathbf{B}$ a $Q\times C$ matrix. Each row in $\mathbf{B}$ is a DC (i.e., what we assumed are the same for all patients ); a column in $\mathbf{A}$, on the other hand, is the time-course associated with a DC across a number of patients (i.e., how strongly a DC is expressed in a patient at a given age). See Figure 3 for an illustration of our NMF analysis, and Table 1 for a glossary of some of the terms introduced in this paper and frequently used in for describing the methodology.

Note that, in this study we focus on disease-based temporal phenotyping, but the approach we introduce can easily incorporate additional phenotypes (e.g., medications) to the input matrix and result in temporal phenotypes that are not disease-only.

2.5 Ascendency Analysis

While mining disease clusters is a key objective of multimorbidity research, providing a view into how such clusters evolve over time is another important goal of multimorbidity research that is less studied. Given that our approach results in a time course for each cluster, we aim to treat these clusters as nodes and use their time courses to find their connectivity. More specifically, we employed a simple network modelling technique, which was originally introduced by Patel et al [48] for the study of the networks in the brain; it was shown to outperform many network modelling techniques in both finding edges and their directions in a range of settings where ground truth for network edges are available (i.e., a simulation study) [49].

Imagine we have two DCs, with binary time courses $\mathbf{v}$ and $\mathbf{w}$ (each T*N dimensional); in our method, if $\mathbf{v}$ and $\mathbf{w}$ are active together and inactive together, we consider them associated/connected. This concept is measured by $\kappa_{\mathbf{v},\mathbf{w}}\in[-1,1]$, and can be seen as an undirected edge in a network; it monotonically increases with joint activation of $\mathbf{v}$ and $\mathbf{w}$, as their individual activations remain fixed; it monotonically decreases with $\mathbf{v}$’s activation as $\mathbf{w}$’s activation and the joint activation remain fixed, or conversely, it monotonically decreases with $\mathbf{w}$’s activation as $\mathbf{v}$’s activation and the joint activation remain fixed; it will be 0 when $\mathbf{v}$ and $\mathbf{w}$ are statistically independent. More formally, in this approach a bivariate Bernoulli Bayesian model is constructed for the joint activation of each pair of time courses using a multinomial likelihood with a Dirichlet prior. The data to model the joint activation/inactivation probabilities for the two time courses be:

[TABLE]

The multinomial likelihood of our data takes the form:

[TABLE]

where each parameter $\theta$ is defined as

[TABLE]

$\kappa$ is defined as a ratio, with its numerator measuring the difference between the joint activation probability and the expected joint activation probability under independence, and its denominator forcing it to range from -1 to 1, i.e.,

[TABLE]

where $E=(\theta_{1}+\theta_{2})(\theta_{1}+\theta_{3})$, $max(\theta_{1})=min(\theta_{1}+\theta_{2},\theta_{1}+\theta_{3})$, $min(\theta_{1})=max(0,2\theta_{1}+\theta_{2}+\theta_{3}-1)$ and

[TABLE]

Given the connectivity between two nodes with $\mathbf{v}$ and $\mathbf{w}$ time courses, if $\mathbf{v}$ exhibits elevated activity for a subset of the period in which $\mathbf{w}$ exhibits elevated activity, we consider $\mathbf{w}$ to be ascendant to $\mathbf{v}$ and vice versa. This is measured by $\tau_{\mathbf{v},\mathbf{w}}\in[-1,1]$, and $\tau_{\mathbf{v},\mathbf{w}}>0$ means that $\mathbf{v}$ is ascendant to $\mathbf{w}$, whereas $\tau_{\mathbf{v},\mathbf{w}}<0$ means that $\mathbf{w}$ is ascendant to $\mathbf{v}$. More formally:

[TABLE]

2.6 Benchmarks and Evaluations

It is always challenging to assess the quality of unsupervised/exploratory models, particularly when there is limited access to relevant objective data on their goodness. Therefore, in this study, we aim to provide the field with an objective approach to evaluate the goodness of the DCs that result from multimorbidity research. We used two different sources of objective medical knowledge, each providing a list of “comorbid” disease pairs: (1) Jensen et al. [6] concluded a list of 1,556 comorbid disease pairs based on the analysis of a large national health dataset, followed by thorough medical due diligence; and (2) Beam et al. [50] studied a large body of medical literature to conclude their list of 113 comorbid condition pairs. Furthermore, through the study of large body of medical literature Beam et al. [50] also concluded a list of 156 causal pairs, where one disease causes the other. Note that, these are the numbers after converting their disease notations to ICD-10 blocks (i.e., the level at which our analysis has been carried out).

Assume that we denote the $q$’th DC as $\mathbf{q}_{k}=\{w_{qi}\}_{i=1}^{C}$ (where $w_{qi}$ is the weight for disease $c_{i}$ in this cluster), and each disease pair in the look up as $p_{ij}=(c_{i},c_{j})$. We can define a score $C_{L}=a/b$, where $b$ is the number of $p_{ij}$’s such that either $c_{i}$ or $c_{j}$ belong to the top $L$ diseases in any cluster (i.e., diseases corresponding to the top $L$ weights in a cluster), and $a$ is the number of $p_{ij}$’s such that both $c_{i}$ and $c_{j}$ belong to the top $L$ diseases in the same cluster. Using $C_{L}$ one can evaluate the goodness of DCs. Furthermore, given that the optimal values of our modelling pipeline are expected to lead to better $C_{L}$, it can be used to guide the grid search for $\sigma$ and $Q$. We defined a similar score based on causal pairs, to assess the goodness of our ascendancy analysis. Assume that we denote a causal pair from Beam et al. as $p_{i\rightarrow j}$ (meaning $c_{i}$ causes $c_{j}$), and an ascendancy edge from our analyses as $o_{m\rightarrow n}$ (meaning an arrow from cluster $m$ to $n$). We defined $A_{L}=a/b$, where $b$ is the number of $p_{i\rightarrow j}$’s that either $i$ or $j$ belong to the top $L$ diseases in either $\mathbf{q}_{m}$ or $\mathbf{q}_{n}$, and $a$ is the number of $p_{i\rightarrow j}$’s such that diseases $i$ and $j$ belong to the top $L$ diseases in clusters $\mathbf{q}_{m}$ and $\mathbf{q}_{n}$, respectively.

In order to check if the $C_{L}$ and $A_{L}$ values we observe are higher than what would be expected by chance alone, we sampled from their distribution under the null hypothesis ($\mathcal{H}_{0}$) that diseases co-occur at random. Under this $\mathcal{H}_{0}$, we can permute the disease labels in our NMF results (i.e., matrix $\mathbf{B}$), and calculate the corresponding $C_{L}$ and $A_{L}$. We carried out this label permutation for 1,000 times to create a nonparametric null distribution for $C_{L}$ and $A_{L}$, and used this to assess the statistical significance of our observed $A_{L}$ and $C_{L}$ values.

3 Results

Our NMF analysis relies on some key free parameters: $Q$ (i.e., the number of DCs) and $\sigma$. We computed $C_{L}$ for a range of of $Q$ and $\sigma$ values on Jensen disease pairs, which consists of a more comprehensive list of disease pairs than Beam, and hence more suitable for the overall parameter tuning. The resulting $C_{L}$ scores are shown in Figure 4.a; using these results, we see a better performance at $\sigma=3$ and $Q\approx 36$. In order to find the best value for $Q$, we calculated the cophenetic correlation coefficient, which measures the stability of clustering derived from NMF results at a given rank [41] for a range of $Q$ values. According to the results shown in Figure 4.b, we see that the factorisation at $Q=34$ results in the highest score. Thus, the rest of our analyses and results will be based on 34 DCs.

Another free parameter of our analysis is the binarisation threshold that we used to prepare the time courses for ascendancy analysis. In order to set this parameter, we calculated $A_{L}$ for a range of $L$ values at different binarisation thresholds and saw the best results to be at 75%. Therefore, for the rest of our analysis we thresholded/binarised each time course $\mathbf{w}$ at $0.75*max(\mathbf{w}$).

Given that our approach to NMF does not include additional variables such as gender, we decided to carry out the analyses separately for men and women. In our data, we have 994,563 male and 1,209,615 female patients; Figure LABEL:fig:nmf_components shows the heatmap in which each column corresponds to DCs that result from NMF (i.e., matrix $\mathbf{B}$). Note that, corresponding to each one of these clusters we have a time-course that has one value for every patient-age, which shows how that DC is expressed in that patient at that age.

Lastly, to demonstrate how the time-courses that results from our model can be used for the analysis of disease trajectories, we carried out an ascendancy analysis, where each DC is treated as a node and the binarised time courses for each disease pair are used to calculate their corresponding $\kappa$ and $\tau$. The resulting $\kappa$ and $\tau$ can be seen in Figure LABEL:fig:mean_tau_kappa. Connecting the pairs with high $\kappa$ values, and defining the direction of the resulting edges based on their corresponding $\tau$, a network map can be derived that is shown in Figures 5 and 6, for male and female patients, respectively. In both these figures, we chose a threshold for $\kappa$, so that the graph has 60 edges to simplify the visual investigation. Note that for every cluster we show the top 3 diseases (using a bar plot showing their weights). We can see for instance in the network corresponding to male patients (Figure 5) that $DC20_{m}$ in which the top disease is “Diabetes mellitus” leads to $DC24_{m}$ in which the top disease is “melabolic disorders”. We can also see in the network corresponding to female patients (Figure 6) that $DC27_{f}$ in which the top disease is “Disorders of gallbladder, biliary tract and pancreas” leads to $DC18_{f}$ in which the top disease is “Disorders of thyroid gland”.

4 Conclusions and Discussion

In this study, we employed a well-known matrix factorisation technique called NMF to mine the MPs using one of the largest EHR datasets in the world (i.e., CPRD). The key reason behind this research was to provide a simple and effective solution for multimorbidity research. To be more specific, our study attempted to build on the past studies’ learnings, while addressing some of their limitations (e.g., using relatively small data, relying on a narrow observation windows, focusing on a small number of diseases, solely extracting the DCs (instead of temporal phenotyping), and lack of appropriate quantitative benchmarking and evaluation of the results). To the best of the authors’ knowledge, our approach is the first one in the literature that has used a matrix (as opposed to tensor) factorisation technique for temporal phenotyping and the study of MPs’ temporal patterns.

Enabling NMF to result in temporal phenotyping was achieved through a simple assumption, which is fairly common across similar studies: Every disease belongs to a number of DCs (with a degree of membership); a combination of these clusters for each patient at each year of life explains the observed patterns of diagnoses. This is very similar to a common approach in the study of functional MRI data in neuoimaging known as multi-session ICA through temporal concatenation [51]. In cases where one is looking for DCs that are common across population, without assuming them having a consistent expression across different people’s lifetime, this approach can be even preferred over the tensor-based factorisation techniques, which assumes the same cluster and time course for everyone.

Given a number of DCs, the next important objective of multimorbidity analyses is to conclude a disease network, which summarises how diseases interact with one another and influence each other’s occurrence. While network models have the potential to solve such a problem (when given the time courses for diseases), the definition of nodes can be a challenge for researchers. That is, if we operate in the ICD-10 universe, we can see scenarios where the network can have anything ranging from 22 nodes (at chapter level) to 10,138 nodes (at level 4). We know from the network modelling literature that the search space for finding the best network is of a super-exponential size on the number of nodes (i.e., $O(n!2^{\binom{n}{2}})$) [52]. This makes the optimisation for learning a network of 10K nodes a huge challenge; both in terms of data (relatively small number of patients, and low prevalence and hence co-occurrence of most diseases at this level) and computing. On the other hand, using ICD-10 chapters, which will result in network with 22 nodes, is likely to lead to results that are hard to be clinically meaningful and interpretable; due to the heterogeneity of the diseases they each contain. Operating at levels such as ICD-10 blocks and 2-digit ICD-10 codes, while not suffering from too many nodes, is still likely to have many highly correlated/co-occurring nodes that might make sense to be combined (particularly, given the data and computing challenges that we face) when learning large networks.

In this study, we introduced a new concept for the node (i.e., the DCs resulting from the NMF) and a new framework to mine these nodes’ relationships with each other. This definition of nodes has a few advantages. Firstly, as our analysis suggests, it leads to a relatively small number of nodes, for which the corresponding network will be easier to learn. Secondly, from an empirical point of view, given that such nodes are driven by diseases that usually co-occur, splitting them into sub-nodes is not likely to be the source of any advantage (specially that such a split will make the network more complex and hence more difficult to learn). And lastly, from a clinical perspective, we are implying that diseases tend to happen in clusters and what the network will tell us is the influence of one cluster on another, given the rest of the clusters. This is in correspondence with what many in the clinical world have been arguing for that the definition of diseases today might not be the most accurate one (and hence various research on phenomapping of the diseases) [53]. Based on all these, we carried out our analysis and derived a network, using a simple and yet powerful technique known as ascendancy analysis. Of course, the data from our NMF approach is equally useful for any other network model (e.g., Bayesian networks), and hence there is need for follow-up research on the use of alternative network modelling techniques (an exhaustive list of such techniques have been used and compared in Smith et al [49]).

As in most such analyses, our modelling pipeline relied on a mix of choices and assumptions we made; from preprocessing to factorisation technique and beyond. For instance, NMF has a range of different implementations; a follow-up research on comparing different implementations of NMF can surely improve the approach we introduced. Furthermore, the data preprocessing we introduced, from smoothing to adjustment for count can all be done with some differences that can be subject of future research. Lastly, our approach was focused on $\mathbf{D}$ as defined in Section 2. The earlier works, however, have introduced more concepts (e.g., measurements and medications) to the starting matrix. Given that our approach can simply accommodate additional phenotypes to the input matrix, without necessarily needing to add a new dimension to it (not needing to go from matrix to tensor space) there is an opportunity that can lead to richer analyses that not only can take time/age into account, but also can benefit from additional concepts and interventions that do influence the DCs and their trajectories. However, for the multimorbidity analysis, disease-only analysis was a necessary starting point.

While we showed some of the results on $C_{L}$ for parameter tuning earlier, here we also show both $A_{L}$ and $C_{L}$ for all the disease pairs, using the optimal values of the free parameters. Note that, while ascendancy analysis does not necessarily mean causal relationship, it has empirically been shown to have high correspondence with it. For instance, Smith et al. [49] have shown how ascendancy metrics are better than almost all well-known causal modelling tools, when assessed on simulated data with known ground truth (i.e., simulated data). Furthermore, disease trajectories and various diseases’ relationship with each other are complex; there will be cases where different experts might not agree on whether a relationship is causal or comorbid. Therefore, showing $A_{L}$ and $C_{L}$ for all pairs will provide additional insight about the strengths and weaknesses of our approach. Figure 7 shows these results for a number of $L$ values. For instance, at $L=15$, which represents 17% of our diseases, the results show that our analysis has successfully identified a considerable number of disease pairs in our lookups either as part of the same cluster ($>33\%$), or as linked via ascendancy ($>61\%$).

Lastly, there has been various developments in methods related to our study that can provide multiple new directions for future works. For instance, deep learning’s success in the past few years has led to “deep phenotyping” research on EHR; while such models can help the study of MPs, their use has been limited to learning disease representations (or embeddings) for disease/event predictions. The earlier works in this space, despite not taking time into account [54, 55], have shown that meaningful DCs can be learned from EHR. Similar results have been shown using CNN [56] and RNN [57, 58, 59]. Furthermore, there has been multiple neural solutions for matrix factorisation [60, 61]. Another novel methodology that has the potential to improve such research is Temporal Regularised Matrix Factorisation (TRMF) [62], which has the ability to regularise the temporal aspect of the factors so it is influenced by prior knowledge/assumptions such as smoothness.

Acknowledgements

This research was funded by the Oxford Martin School (OMS) and supported by the National Institute for Health Research (NIHR) Oxford Biomedical Research Centre (BRC). The views expressed are those of the authors and not necessarily those of the OMS, the UK National Health Service (NHS), the NIHR or the Department of Health and Social Care. This work uses data provided by patients and collected by the NHS as part of their care and support and would not have been possible without access to this data. The NIHR recognises and values the role of patient data, securely accessed and stored, both in underpinning and leading to improvements in research and care.