In this section, we will introduce the method of constructing a graph network data structure for time-series data proposed in this paper. In the constructed graph network data, the point features and edge weights are respectively extracted and fused to construct the overall features of the patient time-series graph network. Finally, we will introduce the classifier and the overall model training framework built for the downstream multiobjective task. The overall symbols are summarized in Table 1.

To retain the fine-grained information in the patient time-series data to the greatest extent and reduce the introduction of noise information, this paper constructs the clinical time-series data as a graph network structure to input into the model calculation. In the construction of the graph network data structure, the feature space of nodes is first constructed. By disassembling the clinical time-series data, it is found that each data measurement behavior can be divided into three dimensions: measurement time (t), measurement feature (f), and measurement value (v). Based on the calculated time range (T_m), the calculated number of feature channels (F_m), the total number of measured behaviors (V_m), and the overall time-series data of the patient can be expressed as {(t1,f1,v1),(t1,f2,v2),....(ta,fb,vc)}a≤Tm,b≤Fm,c≤Vm. Among them, the measurement time and the measurement feature constitute the time node (P_t) and the feature node (P_f) in the graph network data, respectively, and the measurement value is constructed as the connection edge weight (E_w) between the time node and the feature node. The feature space of the feature node (P_f) is mapped in the feature correlation space and composed of the KL divergence values calculated between the original data within the computational time range across the F_m feature channels of each patient. This study references multiple studies on feature space construction based on KL divergence [19,20]. The use of KL divergence to measure distributional differences between distinct physiological variables (eg, heart rate and temperature) may indeed appear counterintuitive from a clinical perspective. However, it is important to clarify that KL divergence is not used here to infer causal or biological relationships between different features. Instead, it serves as a similarity metric to construct the feature node embeddings in the graph network. Specifically, it quantifies the distributional alignment between different physiological time series within the same patient, which may reflect the degree of systemic physiological integration or dysregulation—a key characteristic of critical illness progression. This approach is supported by prior work in multivariate time-series analysis, where distributional similarity is used to capture latent coordination among physiological systems. The specific construction method is as follows: first, select feature f; second, compute the KL divergence values between feature (f) and all other features to form the point feature vector for this feature; repeat this process until the point feature vectors for all features are constructed. The feature distributions used for a specific prediction time point are calculated strictly using data recorded strictly before that (T_m) time. For example, when constructing the graph network for prediction at time (T_m), only measurements with time stamps < T_m are used to compute the KL divergence and edge weights. No future information is incorporated, thereby preventing temporal data leakage. The feature space integrates the information correlation and distribution correlation between feature channels and establishes the relationship between the edge weight value and the point feature space. At the same time, the KL divergence values in the overall feature space are normalized to ensure that the feature distribution is in the feature correlation space and reduce the impact on model training. The calculation method of the feature space (P_f) is as follows:

Probability distribution construction for KL divergence under sparse data: for each feature channel, let its available measurements before prediction time be {v1, v2,...vn}. This paper constructs a discrete probability distribution over a fixed set of M bins spanning the global range of all features. When n=2, we use kernel density estimation with a Gaussian kernel; when n=1, the distribution is a Gaussian centered at that single value; and when n=0, we assign a uniform distribution. The KL divergence between 2 features is then computed on these discrete distributions. Where, S_pf represents the overall spatial feature set, and F^pf_b represents the feature vector of a feature node. Meanwhile, the feature space mapping of the time mode P_t is in the temporal difference embedding space, built by the difference between the time node and the measured time of the remaining nodes. The feature space is used to reflect the correlation of the time values in the measured time series. The calculation time point is strictly limited within the calculation time range. The information provided by the time node can assist the model in mining the time scale between different measurements and provide the change information of the time scale for the model training. The calculation method of the feature space P_t is as follows:

Where, S_pt represents the overall spatial feature set, and F^pt_a represents the feature vector of a time feature node. The connection edge E_w of the graph network is composed of v_c in the measurement behavior, and the existence of the connection edge is directly related to the existence of the measurement behavior. In this way, the fine-grained information in the data is retained to the greatest extent, and the introduction of noise information is eliminated. The construction method of the graph network is shown in Figure 1.

After constructing a graph network with time-series data, this paper proposes the construction of a graph convolutional network (GCN) by using the point feature space and the edge weight connection relationship to extract the overall features in the point feature space. At the same time, a multihead attention mechanism network with a gated mechanism connecting multiple layers is proposed for feature extraction in the edge weight feature space. Finally, a learnable weight head is added to the features of the 2 dimensions for feature fusion to meet the downstream multiobjective tasks.

First, in the GCN, the feature spaces S_pf and S_pt are jointly constructed as the feature matrix X_f of the nodes. The collected values of each feature for an individual patient at different time points are normalized, and the connection values corresponding to time nodes and feature nodes form the edge weight matrix A_w, and the learnable parameter Wis added. The GCN [21] network with L) layer input edge weights is used to construct and calculate the overall feature F_p of the node feature space. The calculation is as follows:

Second, during the computation process of the GCN network, the model primarily focuses on the intrapatient variations, and the edge features tend to gradually weaken as the training depth increases. To fully express the features in the edge weight space, this paper proposes a targeted approach. First, each feature of the edge weight (E_w) is independently normalized across the entire set of patient samples. This normalized edge weight data is then used to construct the feature value matrix (W_E). Finally, (W_E) is input into the convolutional gated attention mechanism network.

During the computation process of the convolutional gated attention mechanism network, first, W_E is fed into the forgetting gate to construct the partially forgotten feature matrix (W_EF). Then, W_EF and W_E are input into the 2D convolution layer at the same time, and the multidimensional data (C_EF) and (C_E) are generated after convolution and input into the multihead attention mechanism layer [22] to obtain A_E(1). It is constructed as shown in the following formula:

Here, (h₁) represents the weight matrix of the first layer and (b₁) represents the bias term of the first layer. The feature A_E(1) obtained from the first layer is reinput into the convolutional gated network to generate the C_AE after passing through the average pooling layer. C_AE, A_E(1), and C_E are again input into the multihead attention mechanism layer. Repeat the operation of feature extraction through (L) layers, and finally obtain the extracted feature weight (F_e) of the edge feature space. The overall calculation formula is as follows:

Finally, the feature (F_p) obtained from the point feature space and the feature (F_e) obtained from the edge weight space were added to the learnable weight head (H_fp) and (H_fe), respectively, to form the overall feature extracted from the time-series data graph network space. The features are computed as follows:

After obtaining the overall features extracted from the temporal graph network, we build different decoders to satisfy the downstream multiobjective task. For the classification prediction task, it is necessary to trace the overall feature changes of patients and conduct a comprehensive analysis. For the 2 different classification tasks of patient mortality and patient clinical representation, the overall feature matrix (F_total) is decomposed and arranged by time and mapped to the corresponding dimension through the feedforward neural network for training. The binary cross-entropy loss function is used to calculate the task loss value, and the calculation formula is as follows:

For the regression task, the key consideration is the feature vector of patients at the final time point. For regression prediction tasks such as predicting the future length of hospital stay or the future change of different physical sign data, this paper uses the features of the last dimension of the point feature space and edge weight space, which are weighted by H_fp and H_fe and then mapped to the corresponding dimension by the feedforward neural network for training. The calculation formula is as follows:

The length of hospital stay was normalized before model training using minimum-maximum scaling. Specifically, let (y_st) denote the original stay time in hours. The normalized value (y_norm) is computed as:

The overall model construction and training are shown in Figures 2A and 2B.

This research was conducted in accordance with the Declaration of Helsinki and was approved by the Ethics Committee of the First Affiliated Hospital of Wenzhou Medical University (approval KY2023-R061). The requirement for written informed consent was waived by the Institutional Review Board due to the retrospective nature of the study. The study data are anonymous or deidentified. Meanwhile, no identification of individual images in the paper or additional material is possible.

In the experimental data part of this paper, the time-series data of patients extracted from the public datasets Medical Information Mart for Intensive Care (MIMIC) III [23,24] and MIMIC-IV [25] are used. Under the guidance of the clinical practice experience of clinicians, this paper selects 42 important features of concern in clinical practice for model construction. Meanwhile, the patient’s survival, length of hospital stay, and clinical phenotypes are constructed as the outcome variables of the model. Since the recording time of the MIMIC-IV dataset is in the transition period of characterizing the coding rules, it is split into two parts: ICD-9 (International Classification of Diseases, 9th Revision) and ICD-10 (International Classification of Diseases, 10th Revision). The data of patients collected by the Critical Care Medicine Department of the First Affiliated Hospital of Wenzhou Medical University from 2012 to 2021 are also used to verify the effectiveness of the model. The data of 2000 patients in each dataset are selected in a 1:1 ratio of death to survival to explore the effectiveness of the proposed model in a balanced and fair way. Overall statistics for the data are provided in Table 2. Among them, the format of categorical data is expressed as “number of patients to number of nonpatients,” and the format of numerical data is expressed as “mean (SD).” To verify the stability and effectiveness of the model, the proposed model was evaluated on an expanded dataset of 4000 samples with a 1:1 ratio of death to survival, analyzing temporal feature data more than 6-hour and 12-hour time windows. The results were compared with those of state-of-the-art methods, and the detailed comparison is provided in Multimedia Appendix 1.

In this paper, the proposed model is compared with advanced time-series prediction models in two single- and multi-objective tasks: predicting the final death of a patient and predicting whether multiple phenotypes of a patient occur. The study explores the advantages and scope of application of the proposed model. The experimental equipment adopted in this paper includes an Intel(R) Xeon(R) Silver 4210R CPU running at 2.40 GHz, an NVIDIA GeForce RTX 3090 graphics card with 24 GB of video memory, and 16 GB of system memory. The proposed model in this paper is compared with state-of-the-art time-series prediction models designed for irregularly sampled data, including time-aware long short-term memory [26], gated recurrent unit with decay [27], multitime attention networks for irregularly sampled time series [28], and ContiFormer [29]. Additionally, we adopt time-series analysis methods such as fully convolutional networks [30], temporal convolutional networks [31], neural basis expansion analysis for interpretable time series forecasting [32], Crossformer [33], and DSformer [34], which are trained on regularly sampled data with missing values imputed by the mean. The final comparison results provide a more multidimensional perspective for evaluating time-series prediction performance. For the time range of data length, this paper refers to the experience in the clinical practice process to construct 5 important time period lengths for training and comparison. The hyperparameter tuning of the model refers to multiple state-of-the-art hyperparameter optimization methods [35-39], and the optimal training parameters for each model are obtained. Detailed data are provided in Multimedia Appendix 2. The accuracy and area under the curve for predicting death for each model are provided in Tables 3 and 4.

From the results, it can be seen that the proposed model is in a leading position in most tasks of the 4 datasets, and the results of the proposed method are equal to or slightly inferior to the results of the comparison methods in some time periods of some datasets. Meanwhile, the proposed model demonstrates accuracy in predicting mortality on a larger-scale dataset. The detailed comparison results are provided in Multimedia Appendix 1. From the analysis, it can be found that death is a simple binary classification task, and most of the time-series analysis methods can effectively capture the change law of time-series data to make the correct judgment. The wrong part of the judgment may be due to changes in the patient’s condition, resulting in the patient’s condition getting better or suddenly worse. Moreover, in the public datasets, the clinical time-series data of patients are collected with high frequency and high data density, and the results of various time-series prediction methods are basically similar. However, in the private dataset, due to different clinical policies, the data collection frequency of patients is relatively sparse, and there are many missing values in the normal time-series arrangement, which leads to poor and unstable prediction results. The proposed algorithm can predict stably and efficiently even when there are many vacant values.

For the multitarget classification task of predicting whether a patient‘s phenotype occurs, the average accuracy, area under the curve, and precision-recall area under the curve (PRAUC) of each target are used as evaluation metrics. PRAUC is particularly suitable for imbalanced classification tasks, as it focuses on the model’s performance on the minority class—a common scenario in clinical phenotype prediction where positive cases are often rare. The prediction performance of each model is summarized in Tables 5-7.

From the Tables 5 and 6, it can be seen that the proposed model has achieved superior results in different datasets and different time periods, which fully reflects the advantages of the proposed algorithm in the face of complex downstream tasks. As shown in Table 7, the proposed method consistently achieves the highest PRAUC across most datasets and time windows, indicating its robustness in handling severe class imbalance. Meanwhile, the proposed model demonstrates accuracy in predicting multiple clinical phenotypes on a larger-scale dataset. The detailed comparison results are provided in Multimedia Appendix 1. Comparing different time-series analysis algorithms, the sampling-and-filling method does cause the loss of fine-grained information. However, the data input structure of the graph network proposed in this paper can meet the required input of the model while retaining the fine-grained information in the time-series data to the greatest extent. At the same time, the experimental results show that the simpler the model structure may produce better results when computing time-series data with sampling and filling.

In the experimental architecture of this paper, predicting the length of hospital stay of a patient is set as a regression prediction analysis task to explore the prediction accuracy of the model in the face of different types of downstream tasks. At the same time, the model can provide multidimensional information required by clinical practice. In this paper, the normalized length of hospital stay of patients is set as the prediction target, and the mean absolute error (MAE) value between the true value and the predicted value is used as the standard to evaluate the quality of the model. The comparison of regression predictions is shown in Table 8. For each dataset, the corresponding minimum and maximum values were used for scaling. Detailed dataset-specific normalization parameters are provided in Multimedia Appendix 3.

From the analysis of the table, it can be found that the prediction accuracy of the model proposed in this paper is in the leading position in most datasets and time periods and is slightly inferior to the best prediction results in a small number of time periods. This finding fully reflects the stability of the model proposed in this paper in the face of different downstream tasks. Meanwhile, the proposed model demonstrates accuracy in predicting the length of hospital stay on a larger-scale dataset. The detailed comparison results are provided in Multimedia Appendix 1. Compared with other models, it can be found that some time-series analysis models have good accuracy in classification experiments, but have large deviation values in regression prediction tasks. In this paper, the prediction indicators of different time periods of each dataset are averaged, and the regression prediction and classification prediction are projected into the same set of score space to verify the stability of the model in the face of multitask learning.

In the ablation experiment of the model, the point feature space and edge weight space proposed in this paper are experimented on using the above 3 regression and prediction analysis tasks (death, phenotype, and stay time) on 4 datasets with different time lengths, and the contribution of feature vectors extracted from different feature spaces to the overall model is explored. Meanwhile, this paper also analyzed the connection results of 2 feature spaces under the learnable head and the direct connection results of 2 feature spaces. Results from top to bottom, the accuracy for predicting death, the average accuracy for predicting the occurrence of different representations of patients, and the MAE for predicting the length of hospital stay were arranged. The values are given in Table 9.

From the analysis of the table, it can be found that both the point feature space and the edge weight space can achieve certain prediction accuracy in the classification task. However, due to the characteristics of the graph network space and the fact that the GCN network is easy to fall into local graph space, the features extracted from the point feature space lead to unstable prediction results. The lack of measurement time information in the edge weight space also leads to low prediction accuracy of the extracted features. The combination of the two can effectively construct a stable and high-accuracy prediction model. In regression prediction analysis, the information provided by the point feature space is limited, and the results of regression prediction have large deviations. Under the action of the learnable head, it can effectively fuse with the feature vector extracted from the edge weight space. It improves the prediction accuracy and prediction stability of the model. Meanwhile, a single-head attention mechanism network and a single-layer GCN are supplemented to verify the effectiveness of hyperparameter configuration and the construction of the multilevel network framework. The results are shown in Table 10.

Aiming at the problem of uneven sampling time intervals and a large number of missing values in medical time-series data, this paper proposes to transform medical time-series data into a graph network structure for input. The feature vectors are extracted by constructing algorithms respectively for the point feature space and the edge weight space in the graph network data structure and are fused to form the feature vector of the graph network space under the action of the learnable head. The fused features are then decoded by the decoder to cope with different downstream tasks. In the exploration of model interpretability, this study extracts each feature at different time points from the trained model. The extracted features are fitted to different target tasks using logistic regression, and Shapley Additive Explanations analysis is used to obtain the changes in the contribution of features to the final target task across various time points. Figure 3 illustrates the variation of feature contributions to mortality prediction at different time scales based on 48-hour data from the MIMIC-IV ICD-10 dataset. Results for the entire dataset are provided in Multimedia Appendix 4.

As analyzed from Figure 3, heart rate and glucose indicators are important factors affecting patient survival in the early stage of hospitalization. With the passage of time, carbohydrate metabolites and patient age gradually become the dominant factors.

In this paper, the proposed model is tested in 5 different data time ranges in 4 datasets, and it is found that the proposed model has good prediction accuracy and prediction stability in classification and regression prediction tasks. Meanwhile, paired DeLong tests were performed on the mortality probability prediction results of the proposed model and comparative models in the MIMIC-IV ICD-9 and private datasets. The results showed that the P value of the test in MIMIC-IV ICD-9 was <0.001 and that in private data was <0.05. Additionally, paired U-tests were conducted on the MAE values of the regression predictions for length of hospital stay, revealing that P values in both MIMIC-IV ICD-9 and private data were <0.05. The specific value of P is provided in Multimedia Appendix 5. These findings fully demonstrate the effectiveness of the proposed model in improving prediction performance. Moreover, this paper conducts ablation experiments on the point feature space model and the edge weight space model and finds that the accuracy and stability of the classification prediction of the model come from the simultaneous contributions of the 2 spaces. The accuracy of regression prediction comes from the adjustment and fusion of the 2 spatial feature vectors by the learnable head.

In comparing our proposed model with existing methods, it is important to acknowledge a confounding factor in the baseline experimental design. Our graph-based model operates directly on raw, irregularly sampled time-series data without any imputation, whereas standard deep learning models, such as a fully convolutional network and temporal convolutional network, are trained on data that have been preprocessed with mean imputation to enforce regular sampling intervals. This design choice conflates 2 distinct sources of potential performance gain: the architectural advantage of the graph network and the representational benefit of avoiding information loss due to imputation. The primary rationale for this comparison was to reflect a realistic clinical deployment scenario, where missing values are prevalent, and imputation remains the default preprocessing step for conventional models. As such, the current benchmark captures the end-to-end practical advantage of a model that inherently accommodates irregular data without the need for imputation. However, we acknowledge that this setup does not isolate the contribution of the graph architecture alone. Despite this limitation, we believe the current comparison remains clinically meaningful, as it reflects the performance gap between models that require imputation and those that do not in real-world settings.

The model proposed in this paper also has the following limitations: (1) since the method proposed in this study fully uses the information in the data to construct nodes and edges and uses both nodes and edges for computation simultaneously, it imposes certain requirements on the training space; and (2) as the number of computational features and the scope of computation time increase, the demand for computational space grows exponentially.

In the future, this study intends to develop a graph network construction paradigm with a more compact structure, which can be compatible with longer time-series data and more features, thereby enabling the tracking and monitoring of patients’ overall disease progression. While the proposed model demonstrates strong predictive performance on irregularly sampled medical time-series data, its scalability to large-scale clinical databases remains a challenge. As noted in the Limitations subheading, the current graph construction and joint computation over node and edge feature spaces impose significant memory overhead that grows exponentially with the number of features and time points. To facilitate real-world deployment in hospital informatics environments that routinely handle tens of thousands of patients, several optimization strategies can be explored. First, sparse graph representations could be adopted to retain only clinically meaningful edges, reducing the memory footprint. Second, feature or temporal down-sampling techniques may help control the node count while preserving essential information. Third, distributed computing frameworks or graph mini-batching strategies could partition the computation across multiple devices, enabling parallel processing of large patient cohorts. Finally, memory-efficient graph convolution operators, such as those based on neighborhood sampling or low-rank approximations, offer promising avenues for scaling graph networks to larger datasets. In future work, we aim to investigate these directions to enhance the scalability and clinical utility of the model without compromising predictive accuracy.

Based on the model’s prediction results, analyzing patients’ clinical time-series data from the early time period can reveal the disease progression for many patients. In addition to carefully examining the characteristic representations of patients identified by the model as having higher mortality risk, it is also possible to focus on predicting treatment regimens that are likely to be misapplied to a given patient. Through the similarity of the patient’s graph network to form the edge and through the hierarchical clustering method, a good treatment plan is mined, and the wrong treatment plan is excluded. Overall, this paper proposes a deep learning–based temporal data analysis model for accurately predicting the survival outcomes of critically ill patients, while also providing critical features that require close attention at various time points.