TASTE: Temporal and Static Tensor Factorization for Phenotyping Electronic Health Records

3.3.1. Solution for factor matrix Q_k.

We can rewrite objective function (3) with respect to Q_k based on trace properties [22] as:

Removing the constant terms and applying the trace property Trace(ABC) = Trace(CAB) yields the following new objective:

The optimal value of Q_k can then be computed via the Orthogonal Procrustes problem [21] which has the closed form solution $Q_k=B_k{C}_k^T\text{where}{B}_k\in R^{I_k\times R}\text{and}{C}_k\in R^{R\times R}$ are the right and left singular vectors of ${\mu }_k{U}_k{H}^T$. Note that each Q_k can contain negative values.

3.3.2. Solution for factor matrix H.

The objective function with respect to H can be rewritten as an unconstrained problem:

Note that Equation 6 is different than the original formulation introduced in Equation 3, where the Frobenius norm contains the term U_k − Q_kH. Through this reformulation, TASTE can utilize the least square solver to efficiently update H.

To obtain the new objective function, we observe that $Q_k\in {ℝ}^{I_k\times R}$ is a rectangular orthogonal matrix $(Q_k^T{Q}_k=I\in {ℝ}^{R\times R})$. We introduce a new orthogonal matrix, $\widetilde{Q_k}\in {ℝ}^{I_k\times \left(I_k-R\right)},\text{where}{\widetilde{Q_k}}^T\widetilde{Q_k}=I\in {ℝ}^{I_k-R\times I_k-R}\text{and}{\widetilde{Q_k}}^T{Q}_k=0$. This can be used to produce a new square orthogonal matrix $[\begin{array}{ll}{Q}_k\hfill & \widetilde{Q_k}\hfill \end{array}]$.

Since $[\begin{array}{ll}{Q}_k\hfill & \widetilde{Q_k}\hfill \end{array}]$ is a square orthogonal matrix (shown in Equation (7)), we can now demonstrate that Equation (6) and Equation (3) are equivalent objectives for H.

where ${\sum }_{k=1}^K{\Vert {\widetilde{Q_k}}^T{U}_k\Vert }_F^2$ is a constant and independent of the parameter under minimization. Therefore, the value of H that minimizes ${\sum }_{k=1}^K\frac{{\mu }_k}{2}{\Vert Q_{k}H-U_k\Vert }_F^2$ also minimizes ${\sum }_{k=1}^K\frac{{\mu }_k}{2}{\Vert H-Q_k^T{U}_k\Vert }_F^2$ and the update rule for factor matrix H is based on the least square solution and has the following form:

3.3.3. Solution for phenotype evolution matrix U_k.

After up-dating the factor matrices Q_k, H, we focus on solving for U_k. In classic PARAFAC2 [12, 13], this factor is retrieved through the simple multiplication U_k = Q_kH. However, for improved interpretability, we prefer temporal factor matrix U_k to be non-negative because the temporal phenotype evolution for patient k (U_k) should not be negative. As shown in the empirical results, a naive enforcement of non-negativity (max (0, Q_k H)) violates the important uniqueness property of PARAFAC2. Therefore, we consider U_k as an additional factor matrix, constrain it to be non-negative, and minimize its difference to Q_k H.

The objective function with respect to U_k can be combined into the following NNLS form:

As mentioned earlier, factor matrix U_k is updated based on block principal pivoting method.

3.3.4. Solution for temporal phenotype definition V.

Factor matrix V defines the participation of temporal features in different phenotypes. In Equation (3), the factor matrix V participates in the PARAFAC2 part with the non-negativity constraint. Therefore, the objective function for factor matrix V has the following form:

To update V based on block principal pivoting, the algorithm calculates ${\left(U_k{S}_k\right)}^T\left(U_k{S}_k\right)\text{and}{U}_k{S}_k{X}_k$ for all K samples which can be done in an embarrassingly parallel fashion.

3.3.5. Solution for factor matrix W or {S_k}.

The objective function with respect to W yields the following format:

Factor matrices {S_k} are shared between the PARAFAC2 input and matrix A where W(k,:)=diag(S_k). Since vec(U_kS_kV^T) = (V ⊙ U_k)W (k, :)^T, Equation (11) can be rewritten in the following NNLS form:

where $\odot$ denotes Khatri-Rao product. Each row of factor matrix W (W(k, :) or diag(S_k)) can be solved separately and in parallel. Unfortunately, the update for each factor matrix S_k involves computing two time-consuming operations: $1){\left(V\odot U_k\right)}^T\left(V\odot U_k\right)$ and $2){\left(V\odot U_k\right)}^T\mathrm{vec}\left(X_k\right)$. Instead of explicitly forming the Khatri-Rao product, both operations can be replaced with more efficient counterparts. The first operation can be replaced with $V^{T}V*U_k^T{U}_k$ where * denotes the element-wise (Hadamard) product [23]. The second operation also can be replaced with diag $(U_k{X}_k{V}^T)$ [23]. Thus, each row of W can be efficiently updated via block principal pivoting.

3.3.6. Solution for static phenotype definition F.

Finally, factor matrix F represents the participation of static features for the phenotypes. The objective function for factor matrix F has the following NNLS form:

which can be easily updated via block principal pivoting.

Sutter:

This dataset is from Sutter Palo Alto Medical Foundation, a large primary care and multispecialty group practice. The data set contains the EHRs for patients with new onset of heart failure and matched controls (matched by encounter time, and age). It includes 5912 cases and 59300 controls. For all patients, encounter features (e.g., medication orders, diagnoses) were extracted from the electronic health records. We use standard medical concept groupers to convert the available ICD-9 or ICD-10 diagnosis codes to Clinical Classification Software (CCS level 3) [24]. We also group the normalized drug names based on unique therapeutic sub-classes using the Anatomical Therapeutic Chemical (ATC) Classification System. Static patient information includes their gender, age, race, smoking status, alcohol status and BMI.

Centers for Medicare and Medicaid (CMS):²

The second data set is CMS 2008–2010 Data Entrepreneurs’ Synthetic Public Use File (DE-SynPUF). The goal of CMS data set is to provide a set of realistic data by protecting the privacy of Medicare beneficiaries by using 5% of real data to synthetically construct the whole dataset. We extract the ICD-9 diagnosis codes and convert them to CCS diagnostic categories as in the case of Sutter dataset.

RMSE:

Accuracy is evaluated as the Root Mean Square Error (RMSE) which is a standard measure used in coupled matrix-tensor factorization literature [25, 26]. Given an input collection of matrices $X_k\in {ℝ}^{I_k\times J},\forall k=1,\dots ,K$ and a static input matrix $A\in {ℝ}^{K\times P}$, we define

X_k (i, j) denotes the (i, j) element of input matrix X_k and ${\widehat{X}}_k(i,j)$ its approximation through a model’s factors (the (i, j) element of the product U_kS_kV^T in the case of TASTE). Similarly, A(i, j) is the (i, j) element of input matrix A and $\widehat{A}(i,j)$ is its approximation (in TASTE, this is the (i, j) element of WF^T).

Cross-Product Invariance (CPI):

We use CPI to assess the solution’s uniqueness, since this is the core constraint promoting it [13]. In particular we check whether $U_k^T{U}_k$ is close to constant $(H^{T}H)\forall k\in \{1,\dots ,K\})$. The cross-product invariance measure is defined as:

The range of CPI is between [−∞, 1], with values close to 1 indicating unique solutions(i.e., ${\text{U}}_k^T{\text{U}}_k$ is close to constant).

Area Under the ROC Curve (AUC):

We examine the classification model’s performance when the data is imbalanced by comparing the actual and estimated labels. We use AUC on the test set to evaluate predictive model performance.

4.4.1. Baseline Approaches:

In this experiment, we compare TASTE with methods that incorporate non-negativity constraint on all factor matrices. Note that SPARTan [10] and COPA [11] are not able to incorporate non-negativity constraint on factor matrices {U_k}.

4.4.2. Setting hyper-parameters:

We perform a grid search for $\lambda \in \{0.01,0.1,1\}\text{and}{\mu }_1=\cdots ={\mu }_K=\mu \in \{0.01,0.1,1\}$ for TASTE and Cohen+ for different target ranks (R = {5, 10, 20, 40}). Each method is run with the specific parameter for 5 random initializations and the best values of λ and μ are selected based on the lowest average RMSE value. For COPA+, we search for the best value of λ ∈ {0.1, 1, 10} since it does not have a μ parameter.

4.4.3. Results:

Apart from purely evaluating the RMSE and the computational time achieved, we assess to what extent the cross-product invariance constraint is satisfied [13]. Therefore, in Figure 2 we present the average and standard deviation of RMSE, CPI, and the computational time for both the Sutter and CMS data sets for four different target ranks (R ∈ {5, 10, 20, 40}). In Figures 2a, 2d, we compare the RMSE for all three methods. We observe that all methods achieve comparable RMSE values on the two different data sets. On the other hand, Figures 2b, 2e show the cross-product invariance (CPI) for Sutter and CMS respectively. COPA+ achieves poor values of CPI for both data sets. This indicates that the output factors violate model constraints and do not satisfy the uniqueness property [13]. Also TASTE significantly outperforms Cohen on CPI in Figures 2b and 2e. Finally, Figures 2c, 2f show the running time comparison for all three methods where TASTE is up to 4.5× and 2× faster than Cohen on Sutter and CMS data sets. Therefore, our approach is the only one that achieves a fast and accurate solution (in terms of RMSE) and preserves model uniqueness (in terms of CPI).

4.6.1. Evaluation Metric: Similarity between two factor matrices:

We define the cosine similarity between two vectors $x_i,y_j\text{as}{C}_{ij}=\frac{x_i^T{y}_j}{\Vert x_i\Vert \Vert y_j\Vert }$. Then the similarity between two factor matrices $X\in {ℝ}^{I\times R},Y\in {ℝ}^{I\times R}$ can be computed as (similar to [13]):

The range of Sim is between [0,1] and values near 1 indicate higher similarity.

4.6.2. Synthetic Data Construction:

We construct the ground-truth factor matrices $\tilde{H}\in {ℝ}^{R\times R},\tilde{V}\in {ℝ}^{J\times R},\tilde{W}\in {ℝ}^{K\times R},\tilde{F}\in {ℝ}^{P\times R}$ by drawing a number uniformly at random between (0,1) to each element of each matrix. For each factor matrix ${\tilde{Q}}_k$, we create a binary non-negative matrix such that ${\tilde{Q}}_k^T{\tilde{Q}}_k=I$ and then compute ${\tilde{U}}_k={\tilde{Q}}_k\tilde{H}$. After constructing all factor matrices, we compute the input based on $X_k={\tilde{U}}_k\mathrm{diag}(\tilde{W}(k,:)){\tilde{V}}^T\text{and}A=\tilde{W}{\tilde{F}}^T$. We set K = 100, J = 30, P = 20, I_k = 100, and R = 4. We then add Gaussian normal noise to varying percentages of randomly-drawn elements ({5%, 10%, 15%, 20%, 25%, 30%, 35%, 40%}) of X_k, ∀k = 1, …, K and A input matrices.

4.6.3. Results:

All three methods achieve the same value for RMSE, therefore, we omit the RMSE versus different noise levels plot. We assess the similarity measure between each ground truth latent factor and its corresponding estimated one (e.g., Sim $(\tilde{V},V)$), and consider the average (·, ·) measure across all output factors as shown in Figure 4a. We also measure CPI and provide the results in 4b for different levels of noise. We observe that despite achieving comparable RMSE, COPA+ scores the lowest on the similarity between the true and the estimated factors. On the other hand, our model achieves the highest amount of recovery, in accordance to the fact that it achieves the highest CPI among all approaches. Overall, we demonstrate how promoting uniqueness (by enforcing the CPI measure to be preserved [13]) leads to more accurate parameter recovery, as suggested by prior work [13, 28].

4.7.1. Cohort Construction:

After applying the preprocessing steps (i.e. removing sparse features and eliminating patients with less than 5 clinical visits), we create a data set from Sutter with 35,113 patients where 3,244 of them are cases and 31,869 are controls (prevalence of 9.2 %). For case patients, we know the date that they are diagnosed with heart failure (HF dx). Control patients also have the same index dates as their corresponding cases. We extract 145 medications, 178 diagnosis codes, and 22 static features from a 2-year observation window and set the prediction window length to 6 months. Figure 5 depicts the observation and prediction windows in more detail.

4.7.2. Baselines:

We assess the performance of TASTE with 6 different baselines.

4.7.3. TrainingDetails:

To calculate the AUC score for our model, we extended 5-fold cross-validation processes (described below) to access how to use phenotyping models to perform HF prediction, by calculating AUC score in the cross validation. At each fold, we take 80% percent of patients as the training set and the remaining 20% as the test set. Figure 6 depicts our heart failure prediction framework which contains 5 steps including:

1. First, we apply TASTE on case patients in the training set and extract the HF phenotypes (V_cases, F_cases) and calculate phenotype score values for them $\left(\{S_{(K_{\mathit{\text{cases}}})}\}\right)$.

2. Second, we assign the existing HF phenotypes (V_cases, F_cases) to the control patient information $(\{X_{(K_{(controls)}})\})$ from training set (based on section 3.4) and extract personalized phenotype scores for control patients $(\{S_{(K_{controls)})}\})$.

3. We train a regularized logistic regression classifier on personalized phenotype scores for all patients in the training set $(\{S_{(K_{cases})}\},\{S_{(K_{controls})}\})$.

4. We assign the existing HF phenotypes (V_cases, F_cases) to the patient information from test set (including cases and controls) and extract their personalized phenotype scores $(\{S_{(K_{test})}\})$.

5. Finally, we predict HF (AUC score) for patients in the test set based on the classifier model trained in step 3. we pick the best parameters (C, λ, μ) based on the highest average AUC score on the test set.

All the other tensor baselines have the same training strategy as TASTE. For all the baselines under comparison, we apply 5-fold cross-validation processes and train a Lasso Logistic Regression ⁴. Lasso Logistic Regression has regularization parameter (C = [1e − 2, 1e − 1, 1, 10, 100, 1000, 10000]). For all 6 baselines, we just need to tune parameter C. However, for TASTE we need to perform a 3-D grid search over λ ∈ {0.01, 0.1, 1} and μ₁ = μ₂ = … = μ_k = μ = ∈ {0.01, 0.1, 1} and C.

Results:

Figure 7 shows the average of AUC for all baselines and TASTE. For COPA, COPA(+static), CNTF and TASTE we report the AUC score for different values of R ({5, 10, 20, 40, 60}). TASTE improves the AUC score over a simple non-negative PARAFAC2 model (COPA and COPA(+static)) and CNTF which suggests: 1) incorporating static features with dynamic ones will increase the predictive power (comparison of TASTE with COPA and CNTF); and 2) incorporating static features using a coupled matrix improves predictive power (comparison of TASTE and COPA(+static)). We also observe that TASTE with R=60 (AUC=0.7687) performs slightly better than the RNN baseline model. Moreover, TASTE offers interpretability as the phenotype definitions can be readily extracted. RNNs require additional mechanisms to explain the model [30].

4.8.1. Cohort Construction:

We select the patients diagnosed with HF from the EHRs in Sutter dataset. We extract 145 medications and 178 diagnosis codes from a 2-year observation window which ends 6 months before the heart failure diagnosis date (HFdx). ⁶ The total number of patients (K) is 3,244 (the HF case patients of Sutter dataset) same as section 4.7.

4.8.2. Pure PARAFAC2 cannot handle static feature integration.

In this experiment, we further analyze the results of the naive way of incorporating static feature information into a simpler PARAFAC2-based framework [11]. We posit that this results in less interpretable phenotypes. We incorporate the static features into PARAFAC2 input by repeating the value of static features on all clinical visits of the patients in the same fashion as COPA(+static). For instance, if the male feature of patient k has value 1, we repeat the value 1 for all the clinical visits of that patient. Then we compare the phenotype definitions discovered by TASTE (matrices V, F) and by COPA (matrix V). Table 3 contains two sample phenotypes discovered by this baseline, using the same truncation threshold that we use throughout this work (we only consider features with values greater than 0.1). We observe that the static features introduce a significant amount of bias into the resulting phenotypes: the phenotype definitions are essentially dominated by static features, while the values of weights corresponding to dynamic features are close to 0. This suggests that pure PARAFAC2-based models such as the work in [11] are unable to produce meaningful phenotypes that handle both static and dynamic features. Such a conclusion extends to other PARAFAC2-based work which does not explicitly model side information [10, 13, 27].

4.8.3. TASTE Findings of HF Phenotypes.

Based on Figure 7, we present the top 5 phenotypes extracted from TASTE using R = 40 due to space limitations⁷. This rank is selected as outperforms all but the RNN baseline and is comparable to R = 60 in terms of performance. The 5 phenotypes are all confirmed and annotated by an expert cardiologist. Table 4 provides the details of these phenotypes. The clinical description of the 5 phenotypes as provided by the cardiologist are: