Sampling-guided Heterogeneous Graph Neural Network with Temporal Smoothing for Scalable Longitudinal Data Imputation

In practice, the process within SHT-GNN can be understood as a form of data reconstruction. From an optimization perspective, this process is conceptually similar to Variational Autoencoders (VAEs), where the goal of the reconstruction step is to minimize the reconstruction error as part of the Evidence Lower Bound (ELBO). As is widely recognized, the training objective of a standard VAE is expressed as:

where $E_{q(z|x)}\log{p(x|z)}$ represents the $**reconstruction loss**$, and $D_{\text{KL}}[q(z|x)\parallel p(z)]$ is the $**regularization term**$. In our proposed SHT-GNN, the calculation and training process can be described as follows:

where $X_{\text{obs}}$ denotes the observed values, $Z_{O}^{\text{init}}$ and $Z_{F}^{\text{init}}$ represent the initial embedding matrices of all observation and feature nodes, respectively. $Z_{O}^{\text{L}}$ and $Z_{F}^{\text{L}}$ denote the embedding matrices of all observation and feature nodes after $L$ layers of forward computation in SHT-GNN. Subsequently, the training objective in SHT-GNN is expressed as:

$\text{That is to maximize}$: $E_{q(Z^L_{O},Z^L_F|X_{obs})}\log{p(X_{obs}|Z^L_{O},Z^L_{F})} - D_{\text{KL}}[p(Z^L_O,Z^L_F|X_{obs}) \parallel p(Z^L_O,Z^L_F)]$

where $E_{q(Z_{O},Z_F|X_{obs})}\log{p(X_{obs}|Z_{O}^L,Z_{F}^L)}$ represents the $**reconstruction loss**$.

In principle, the SHT-GNN and VAE-based methods share conceptual similarities. However, the key distinctions between the SHT-GNN and VAE-based approaches in the handling of longitudinal data are as follows.

$**1. Enhanced Representation Capability through Graph Neural Networks (GNNs)**$. In SHT-GNN, the encoding and decoding processes utilize graph neural networks (GNNs) instead of the multilayer perceptrons (MLPs) employed in standard VAEs. GNNs offer significantly stronger representation capabilities by leveraging the inherent data structure and facilitating temporal smoothing across non-independent observations, which is critical in longitudinal settings.

$**2. Distinct Regularization Term for Non-Independent Observations**$. In VAE-based methods, the regularization term for the latent space distribution, $D_{\text{KL}}[q(z|x) \parallel p(z)]$, quantifies the independent distributions of individual observations. By contrast, SHT-GNN incorporates a KL divergence term, $D_{\text{KL}}[p(Z^L_O,Z^L_F|X_{obs}) \parallel p(Z^L_O,Z^L_F)]$, which measures the joint distribution of all observation representations. To account for dependencies within longitudinal data, SHT-GNN employs multi-layer temporal smoothing and introduces MADGap, enabling improved representation of the joint latent space distribution.

In summary, the foundational intuition behind SHT-GNN lies in optimizing the reconstruction error during the data reconstruction process. Moreover, in the context of longitudinal data, SHT-GNN introduces temporal smoothing and MADGap to strike a balance between temporal consistency and the diversity of observation representations within the latent space.

We have realized that the lack of straightforward presentation and theoretical intuitiveness in our paper has caused difficulties for readers. Therefore, in the revision of our paper, we will include all the aforementioned theoretical understandings as a remark within the article.

1. Problems in presentations

Thank you very much for your valuable suggestions, which are highly important for improving the quality of our paper.

1.1 Regarding the issue with the figures, we will replace them with vector graphics during the revision process, as long as they do not exceed the maximum upload capacity allowed for the PhD submission system.

1.2. For the Jaccard distance, you are correct that we made certain modifications to the original Jaccard distance. In its original form, the denominator should be the union of sets A and B, ensuring the distance is symmetric. However, since our message-passing mechanism along the timeline is directional, we aimed to measure the additional information from past observations relative to the current observation. Therefore, we adjusted the denominator to reflect the observation set corresponding to the receiving observation. We acknowledge that this deviates from the original definition, and we greatly appreciate your suggestion. We will adopt your recommendation and revise the term to "adjusted Jaccard distance."

1.3. Regarding the typos, we sincerely apologize for the errors and the inconvenience they may have caused during your reading. We will correct all the typos you identified and thoroughly review the entire manuscript for spelling and potential representation issues. Once again, thank you for your detailed review and constructive feedback on our paper.

2. Appendix Simulations and their reliance to existing studies

In the process of simulation data generation, the response values were generated using a mix of higher-order and lower-order terms, with the parameters and hyperparameters cited from the following references:

[1] Binder, H., Sauerbrei, W., Royston, P. Comparison between splines and fractional polynomials for multivariable model building with continuous covariates: a simulation study with continuous response. Statistics in Medicine, 2013, 32(13): 2262–2277.

[2] Wu, H., Zhang, J. T. Nonparametric regression methods for longitudinal data analysis: mixed-effects modeling approaches. John Wiley & Sons, 2006.

Specifically, the generation of continuous response variables followed the approach outlined in [1]. The simulation formula proposed in [1] was designed to ensure that the generated data align with a predefined Structural Causal Model. This was achieved by employing a linear or nonlinear combination of covariates’ lower-order terms (e.g., square root or original values), higher-order terms (e.g., squared or cubed values), and interaction terms.

Similar simulation methods have been widely used in studies involving clinical and behavioral data, with the primary goal of generating data that conform to a predefined Structural Causal Model. In this work, the covariate dimensionality in the simulation was relatively high. To increase the complexity of the relationship between covariates and the response variable, the final simulation employed a combination of linear and nonlinear terms.

Subsequently, the process of averaging observations within the same time window and adding noise followed the methodology described in [2].

1. How to Apply Tabular Data Imputation to Longitudinal Data

For tabular data imputation, longitudinal data is concatenated observation-wise, regardless of the subject it originates from. All observations are separated and combined into a single tabular data structure, which is then used as the input for all missing data imputation methods.

2. How to Apply Time Series Imputation Methods to Longitudinal Data

In the tasks of missing data imputation and prediction for longitudinal data, two scenarios can be considered:

For methods that provide code specifically designed for time series across multiple subjects, these can be directly applied to longitudinal data for imputation.

For methods where the code only supports missing data imputation and prediction for a single time series, the imputation method will be applied separately to the observations corresponding to each subject.

The above summaries are drawn from both theoretical and experimental analyses. Furthermore, several time-series datasets with relevant characteristics were identified to evaluate the performance of our method in missing data imputation and prediction tasks. And we also add some other state-of-the-art baselines in longitudinal data/time series imputation for comparison. The supplementary baselines are:

[1] Transformer: Ailing Zeng, Muxi Chen, Lei Zhang, and Qiang Xu. 2022. Are transformers effective for time series forecasting? arXiv preprint arXiv:2205.13504 (2022).

[2] GP-VAE: Fortuin V, Baranchuk D, Rätsch G, et al. Gp-vae: Deep probabilistic time series imputation[C]//International conference on artificial intelligence and statistics. PMLR, 2020: 1651-1661.

[3] SAITS: Du W, Côté D, Liu Y. Saits: Self-attention-based imputation for time series[J]. Expert Systems with Applications, 2023, 219: 119619.

[4] CTA: Wi H, Shin Y, Park N. Continuous-time Autoencoders for Regular and Irregular Time Series Imputation[C]//Proceedings of the 17th ACM International Conference on Web Search and Data Mining. 2024: 826-835.

For each dataset, we introduce missing data under an MCAR mechanism with a missing ratio of 0.3. We subsequently validated the imputation performance using the Mean Absolute Error (MAE) metric.

For each dataset, we also applied MCAR with a missing ratio of 0.5, and subsequently validated the missing data imputation results using MAE.

where the AirQuality dataset consists of a single-object time series with 15 features and over 9,000 records. The Electricity dataset comprises long-term time series data for 370 subjects, containing a total of over 140,256 records, each characterized by a single-dimensional attribute. The Energy dataset captures energy usage from various electrical appliances, representing 110 subjects with more than 10,000 observations and a total of 27 features.

On the AirQuality and Electricity datasets, the SHT-GNN method demonstrates inferior performance compared to RNN and VAE-based approaches. This indicates that the SHT-GNN model is not well-suited for long-term, single-object time series, which fall outside its intended application scenario.

However, SHT-GNN achieves state-of-the-art performance on the Energy dataset. Unlike the AirQuality and Electricity datasets, the Energy dataset involves time series with multiple subjects and multidimensional features. SHT-GNN leverages its ability to effectively borrow information across features, which significantly enhances its performance, particularly in the context of irregular longitudinal data imputation.

Before presenting the supplementary datasets and baselines, I will first summarize the strengths and weaknesses of our proposed method, as demonstrated through experiments and theoretical analysis. Following this, I will present comparative results on three additional datasets and four other longitudinal and time-series datasets using cutting-edge baselines, and analyze the reasons behind these results.

What kind of missing data imputation can our method handle well?

1. MAR, MCAR, and high missing ratio: Our experiments demonstrate that the missing mechanism in the ADNI dataset aligns with MAR (Missing At Random). Through logistic modeling and testing, we identified significant associations between the missing indicator of the response variable and key covariates such as Age and APOE4. This confirms that our method is well-suited for handling both MAR and MCAR conditions. However, we currently lack a theoretical explanation for its performance under MNAR (Missing Not At Random). The ADNI dataset also presents a challenging scenario, with a missing ratio for the response variable exceeding 90%. Furthermore, in our simulations, we tested various missing ratios, and our method adapted well across these scenarios, demonstrating robustness under high missing ratios.

2. Block-wise missing: Block-wise missing data, a common issue in biomedical datasets, poses unique challenges for data imputation. Our sub-sampling method, however, remains robust and unaffected by this issue. A clear instance of block-wise missing data can be observed in the ADNI dataset. Compared to other baseline methods, addressing block-wise missing data in mini-batches typically requires additional preprocessing. For instance, in VAE-based approaches, missing value imputation often involves initializing the missing values. This process becomes particularly challenging when block-wise missing data leads to significant gaps in the data indicators within a mini-batch. In contrast, Graph Neural Networks (GNNs) bypass the need for explicit initialization of missing values. Even in the presence of block-wise missing data within mini-batches, GNNs effectively utilize information from both corresponding observations of the same subject and observations from other subjects. This ability allows GNNs to learn robust and meaningful observation embeddings, ensuring their adaptability to such challenging scenarios.

Scenarios where our method may not be advantageous?

1. Single time series and long-term memorized time series: When faced with the imputation of missing data in a single-subject time series, the inductive learning capability of our proposed method across multiple subjects cannot be fully utilized. Unlike many cutting-edge RNN-based methods that focus on imputing single time series, our method is specifically designed for longitudinal data across multiple subjects, rather than simply filling gaps within the time series of an individual subject or system. Moreover, our model currently has limited ability to capture smooth, long-term temporal features in time series. RNN-based methods achieve this through memory and update components, concentrating all computational effort on the temporal axis while neglecting the relationships between multiple subjects. As a result, our method may not perform as well as RNNs in handling missing data for datasets with very long-term memory characteristics. Furthermore, our model currently has limited capacity to capture smooth, long-term temporal dependencies in time series. RNN-based methods excel in this domain due to their memory and update mechanisms, which focus computational effort exclusively on the temporal axis while neglecting inter-subject relationships. As a result, for datasets characterized by strong long-term memory features, RNN-based methods may outperform our approach in handling missing data.

2. Time series with weak feature correlations and high stability: In datasets such as those used for electricity or air quality monitoring, the data are typically stable, well-aligned, and exhibit weak correlations across features. In such scenarios, RNN-based and VAE-based methods are already highly effective. The strengths of our model lie in its flexibility to handle irregular data, its inductive ability to capture interdependencies among features, and its capacity to learn specific representations for individual observations. Although our method performs competitively on time series with weak feature correlations and high stability, it does not offer a significant advantage over existing approaches in these cases.

Thank you very much for your valuable feedback! Indeed, our method shares similarities with existing data imputation approaches based on graph neural networks (GNNs). Specifically, several prior works have utilized bipartite graphs in GNNs for data imputation, including:

[1] You, J., Ma, X., Ding, Y., et al. Handling missing data with graph representation learning. Advances in Neural Information Processing Systems, 2020, 33: 19075–19087.

[2] Zhong, J., Gui, N., Ye, W. Data imputation with iterative graph reconstruction. Proceedings of the AAAI Conference on Artificial Intelligence, 2023, 37(9): 11399–11407.

[3] Um, D., Park, J., Park, S., et al. Confidence-based feature imputation for graphs with partially known features. arXiv preprint arXiv:2305.16618, 2023.

[4] Gupta, S., Manchanda, S., Ranu, S., et al. GRAFENNE: learning on graphs with heterogeneous and dynamic feature sets. International Conference on Machine Learning, PMLR, 2023: 12165–12181.

In practice, the aforementioned methods adopt bipartite graph-based imputation frameworks, incorporating advanced modules such as confidence scoring, diffusion models, and transformer mechanisms to capture interdependencies among features. Our approach shares commonalities with these GNN-based methods, particularly in its use of edge-wise prediction for missing data imputation and its reliance on bipartite graph structures.

However, a key distinction lies in our method’s focus on longitudinal data, which is often overlooked in existing approaches. Longitudinal datasets are particularly valuable due to their temporal follow-up structure, yet they frequently encounter significant rates of missing data because of the inherent challenges of long-term studies. None of the aforementioned methods specifically address this unique context.

Our method introduces a temporal smoothing mechanism and a custom-designed loss function tailored to the characteristics of longitudinal data. These innovations enable effective modeling of temporal dependencies and inter-subject variability, setting our approach apart from existing GNN-based imputation methods.

Furthermore, our approach has demonstrated substantial advancements when applied to real-world datasets, such as those for Alzheimer's disease research, highlighting its practical utility. Below, we provide both theoretical and intuitive explanations to demonstrate the novelty of our modules and the necessity of these mechanisms for learning in longitudinal data contexts.

The content you mentioned does not match what is stated in the GRAPE paper.

In the original text of the GRAPE paper, it states: "For label prediction tasks, we use two GNN layers with 16 hidden units. $O_{edge}$ and $O_{node}$ are implemented as linear layers. The edge dropout rate is set to $r_{drop} = 0.3$. For all experiments, we run 5 trials with different random seeds and report the mean and standard deviation of the results."

Clearly, the settings we are using here for the label prediction task in GRAPE are two GNN layers with 16 hidden units, not 64.

Thank you very much for your valuable feedback, particularly for pointing out the insufficiency in our baseline and dataset comparisons. Before presenting the supplementary datasets and baselines, I will first summarize the strengths and weaknesses of our proposed method, as demonstrated through experiments and theoretical analysis. Following this, I will present comparative results on three additional datasets and four other longitudinal and time-series datasets using cutting-edge baselines, and analyze the reasons behind these results.

What kind of missing data imputation can our method handle well?

1. MAR, MCAR, and high missing ratio: Our experiments demonstrate that the missing mechanism in the ADNI dataset aligns with MAR (Missing At Random). Through logistic modeling and testing, we identified significant associations between the missing indicator of the response variable and key covariates such as Age and APOE4. This confirms that our method is well-suited for handling both MAR and MCAR conditions. However, we currently lack a theoretical explanation for its performance under MNAR (Missing Not At Random). The ADNI dataset also presents a challenging scenario, with a missing ratio for the response variable exceeding 90%. Furthermore, in our simulations, we tested various missing ratios, and our method adapted well across these scenarios, demonstrating robustness under high missing ratios.

2. Block-wise missing: Block-wise missing data, a common issue in biomedical datasets, poses unique challenges for data imputation. Our sub-sampling method, however, remains robust and unaffected by this issue. A clear instance of block-wise missing data can be observed in the ADNI dataset. Compared to other baseline methods, addressing block-wise missing data in mini-batches typically requires additional preprocessing. For instance, in VAE-based approaches, missing value imputation often involves initializing the missing values. This process becomes particularly challenging when block-wise missing data leads to significant gaps in the data indicators within a mini-batch. In contrast, Graph Neural Networks (GNNs) bypass the need for explicit initialization of missing values. Even in the presence of block-wise missing data within mini-batches, GNNs effectively utilize information from both corresponding observations of the same subject and observations from other subjects. This ability allows GNNs to learn robust and meaningful observation embeddings, ensuring their adaptability to such challenging scenarios.

Scenarios where our method may not be advantageous?

1. Single time series and long-term memorized time series: When faced with the imputation of missing data in a single-subject time series, the inductive learning capability of our proposed method across multiple subjects cannot be fully utilized. Unlike many cutting-edge RNN-based methods that focus on imputing single time series, our method is specifically designed for longitudinal data across multiple subjects, rather than simply filling gaps within the time series of an individual subject or system. Moreover, our model currently has limited ability to capture smooth, long-term temporal features in time series. RNN-based methods achieve this through memory and update components, concentrating all computational effort on the temporal axis while neglecting the relationships between multiple subjects. As a result, our method may not perform as well as RNNs in handling missing data for datasets with very long-term memory characteristics. Furthermore, our model currently has limited capacity to capture smooth, long-term temporal dependencies in time series. RNN-based methods excel in this domain due to their memory and update mechanisms, which focus computational effort exclusively on the temporal axis while neglecting inter-subject relationships. As a result, for datasets characterized by strong long-term memory features, RNN-based methods may outperform our approach in handling missing data.

2. Time series with weak feature correlations and high stability: In datasets such as those used for electricity or air quality monitoring, the data are typically stable, well-aligned, and exhibit weak correlations across features. In such scenarios, RNN-based and VAE-based methods are already highly effective. The strengths of our model lie in its flexibility to handle irregular data, its inductive ability to capture interdependencies among features, and its capacity to learn specific representations for individual observations. Although our method performs competitively on time series with weak feature correlations and high stability, it does not offer a significant advantage over existing approaches in these cases.

The above summaries are drawn from both theoretical and experimental analyses. Furthermore, several time-series datasets with relevant characteristics were identified to evaluate the performance of our method in missing data imputation and prediction tasks. And we also add some other state-of-the-art baselines in longitudinal data/time series imputation for comparison. The supplementary baselines are:

[1] Transformer: Ailing Zeng, Muxi Chen, Lei Zhang, and Qiang Xu. 2022. Are transformers effective for time series forecasting? arXiv preprint arXiv:2205.13504 (2022).

[2] GP-VAE: Fortuin V, Baranchuk D, Rätsch G, et al. Gp-vae: Deep probabilistic time series imputation[C]//International conference on artificial intelligence and statistics. PMLR, 2020: 1651-1661.

[3] SAITS: Du W, Côté D, Liu Y. Saits: Self-attention-based imputation for time series[J]. Expert Systems with Applications, 2023, 219: 119619.

[4] CTA: Wi H, Shin Y, Park N. Continuous-time Autoencoders for Regular and Irregular Time Series Imputation[C]//Proceedings of the 17th ACM International Conference on Web Search and Data Mining. 2024: 826-835.

For each dataset, we introduce missing data under an MCAR mechanism with a missing ratio of 0.3. We subsequently validated the imputation performance using the Mean Absolute Error (MAE) metric.

For each dataset, we also applied MCAR with a missing ratio of 0.5, and subsequently validated the missing data imputation results using MAE.

where the AirQuality dataset consists of a single-object time series with 15 features and over 9,000 records. The Electricity dataset comprises long-term time series data for 370 subjects, containing a total of over 140,256 records, each characterized by a single-dimensional attribute. The Energy dataset captures energy usage from various electrical appliances, representing 110 subjects with more than 10,000 observations and a total of 27 features.

On the AirQuality and Electricity datasets, the SHT-GNN method demonstrates inferior performance compared to RNN and VAE-based approaches. This indicates that the SHT-GNN model is not well-suited for long-term, single-object time series, which fall outside its intended application scenario.

However, SHT-GNN achieves state-of-the-art performance on the Energy dataset. Unlike the AirQuality and Electricity datasets, the Energy dataset involves time series with multiple subjects and multidimensional features. SHT-GNN leverages its ability to effectively borrow information across features, which significantly enhances its performance, particularly in the context of irregular longitudinal data imputation.

In the above responses, we have added the results of all supplementary experiments (including results on the PhysioNet 2012 ICU challenge dataset, results of all additional baselines on synthetic datasets and ADNI dataset, and error-bars of all these results), and we have uploaded the revision that includes all of the above supplementary results.

Once again, thank you for your valuable feedback and advice!! And we hope our responses have addressed all of your concerns.

We are grateful for the insightful feedback and constructive comments provided by the reviewers on our paper.

It is noteworthy that Yao et al. (2005) conducted a functional principal component analysis (FPCA) for sparse and irregular data. In their concluding remarks, they proposed the consideration of a regular design in which many data points may be missing for certain subjects. However, our work differs fundamentally from theirs in both scope and focus. While Yao et al. (2005) emphasized statistical theory and FPCA for longitudinal and functional data, particularly addressing the estimation of mean and covariance functions, our work prioritizes prediction accuracy in large-scale, irregular longitudinal studies characterized by missing data in both covariate and response variables. Specifically,

In functional data analysis, as described by Yao et al. (2005), observed sparse data are treated as samples from random curves that represent an underlying smooth and continuous trajectory (function) over time. The repeated observations of each subject are modeled as realizations of a continuous random function, with the primary goal being to recover and predict the entire smooth trajectory for each subject. In contrast, our approach focuses only on discrete observed time points, with particular interest in capturing trends at specific time points rather than reconstructing continuous trajectories.

Our approach combines graph-based learning with strategies to address missing data imputation challenges, whereas their methodology is rooted in functional data analysis. In the work of Yao et al. (2005), sparse and irregular observations refer to scenarios where entire observations are missing within otherwise continuous and dense data. Once the functional data are estimated, they can be used to generate dense observations to address the absence of observations. However, this approach cannot handle cases where covariates and responses are partially missing within the observed data during estimation and prediction. In this paper, we aim to address this issue, specifically focusing on the problem of partial feature missingness within sparse observations. Our work emphasizes improving the prediction accuracy and reducing the variance of missing data in each individual observation.

1. Clarify the conditions for our method

Our method is capable of handling missing data under MCAR (Missing Completely at Random) and MAR (Missing at Random) assumptions, even in scenarios with a high missing ratio. Our experiments demonstrate that the missing mechanism in the ADNI dataset aligns with MAR (Missing At Random). Through logistic modeling and testing, we identified significant associations between the missing indicator of the response variable and key covariates such as Age and APOE4. This confirms that our method is well-suited for handling both MAR and MCAR conditions. However, we currently lack a theoretical explanation for its performance under MNAR (Missing Not At Random). The ADNI dataset also presents a challenging scenario, with a missing ratio for the response variable exceeding 90%. Furthermore, in our simulations, we tested various missing ratios, and our method adapted well across these scenarios, demonstrating robustness under high missing ratios.

2. Longitudinal observations often become inconsistent due to the following circumstances:

2.1. Missing Data at Specific Time Points: Observations for certain variables may be unavailable at specific time points, resulting in incomplete datasets. This leads to misalignment across time and variables, as not all measurements are available for every subject at each observation time.

2.2. Absence of Entire Observations at Certain Time Steps: Complete data points, including all variables for a specific subject at a given time, may be missing. This disrupts the regular observation schedule, converting a consistent timeline into an irregular one.

2.3. Variability in Observation Frequencies Between Subjects: Subjects may have differing numbers of observations or data collected at varying intervals. For example, one subject might be observed weekly, while another is observed monthly, resulting in inconsistent observation schedules across the dataset.

2.4. Irregular Data Collection Schedules Due to External Factors: Practical constraints, such as missed appointments or technical issues during data collection, can lead to observations being recorded at irregular intervals, further complicating the consistency of longitudinal data.

In the last revision, we have emphasized these inconsistencies in the introduction to underscore their significance in longitudinal data analysis.

1. Does the temporal smoothing method work for longitudinal data with only a few measurements per subject?

Our method is effective in handling scenarios with limited measurements per subject.

Methodology: The foundation of our approach lies in the representation power of graph neural networks (GNNs), which serve as the foundation for both missing data imputation and response prediction. Each observation derives its embedding from the observed values associated with its attributed edges. Leveraging the strong expressive capabilities of GNNs, our temporal smoothing mechanism achieves robust performance without relying on a large number of measurements.

Experimental Evidence: In the ADNI dataset, the number of measurements varies considerably across subjects, ranging from as few as 10–15 follow-ups to as many as 30–40 follow-ups. Our approach is specifically designed to accommodate such variability, demonstrating strong adaptability to the structure of the ADNI data. Meanwhile, in our simulation studies, the GLOBEM human behavior dataset features subjects with only 10–20 measurements. Despite this limitation, our method has consistently exhibited strong performance in processing such longitudinal data, further underscoring its robustness and capacity to handle sparse observations.

2. Can your method handle noise in longitudinal measurements?

Indeed, our method is capable of handling noise in longitudinal measurements and has demonstrated excellent performance on simulation datasets with noisy data. In particular, we introduced noise during the generation of response variables, as detailed in Section 5.2 of our paper, which provides a comprehensive description of the process used to incorporate noise into the simulations.

We are able to provide a complete theoretical understanding of our method and explain its applicability to MAR from a theoretical perspective. For all of the following content, we will include it in our revision:

From an optimization perspective, the missing data imputation process in our method is conceptually similar to Variational Autoencoders (VAEs), where the goal of the reconstruction step is to minimize the reconstruction error as part of the Evidence Lower Bound (ELBO). As is widely recognized, the training objective of a standard VAE for missing data imputation is expressed as:

where $E_{q(z|x)}\log{p(x|z)}$ represents the reconstruction loss, and $D_{KL}[q(z|x)\parallel p(z)]$ is the regularization term. When maximizing $\log{p(x_{obs})}$, it is guaranteed that the estimated results for missing data will be consistent in both MCAR and MAR scenarios, a point that has been emphasized in many studies [1][2].

[1] Mattei P A, Frellsen J. MIWAE: Deep generative modelling and imputation of incomplete data sets[C]//International conference on machine learning. PMLR, 2019: 4413-4423.

[2] Collier M, Nazabal A, Williams C K I. VAEs in the presence of missing data[J]. arXiv preprint arXiv:2006.05301, 2020.

In our proposed SHT-GNN, we are also theoretically optimizing $\log{p(x_{obs})}$. Specifically, the calculation and training process can be described as follows:

$\text{That is to maximize}$: $E_{q(Z^L_{O},Z^L_F|X_{obs})}\log{p(X_{obs}|Z^L_{O},Z^L_{F})} - D_{\text{KL}}[p(Z^L_O,Z^L_F|X_{obs}) \parallel p(Z^L_O,Z^L_F)]$

where $X_{obs}$ denotes the observed values, $Z_{O}^{\text{init}}$ and $Z_{F}^{\text{init}}$ represent the initial embedding matrices of all observation and feature nodes, respectively. $Z_{O}^{\text{L}}$ and $Z_{F}^{\text{L}}$ denote the embedding matrices of all observation and feature nodes after $L$ layers of forward computation in SHT-GNN. Subsequently, the training objective in SHT-GNN is expressed as:

where $E_{q(Z_{O},Z_F|X_{obs})}\log{p(X_{obs}|Z^L_{O},Z^L_{F})}$ represents the reconstruction loss. Here, $E_{q(Z_{O},Z_F|X_{obs})}\log{p(X_{obs}|Z^L_{O},Z^L_{F})}$ is the joint distribution over all observations, which differs from the $E_{q(z|x_{obs})}\log{p(x_{obs}|z)}$ in VAE.

Previously, in the case of VAE, their expectation is calculated on the sample level, typically using the MSE over the observed values of all samples to approximate the reconstruction loss. The specific form is:

$\text{Loss} = \frac{1}{N} {\sum_{i=1}^N \sum_{j=1}^p {m_{ij}}\cdot(\hat{X_{ij}}-X_{ij})^2}$

where $m_{ij}$ is the missing indicator for the $j$-th feature in the $i$-th observation.

In contrast, in SHT-GNN, our reconstruction loss is in the form of a joint distribution, and it is not possible to estimate it by averaging over the samples. This is why we use edge dropout in missing data imputation, because directly using the loss over all observed edges will not provide an effective estimate of $E_{q(Z_{O},Z_F|X_{obs})}\log{p(X_{obs}|Z^L_{O},Z^L_{F})}$. Specifically, we estimate the overall reconstruction loss by randomly calculating the loss on some edges in each batch of different heterogeneous graphs.

In principle, the SHT-GNN and VAE-based methods share conceptual similarities. Furthermore, the reconstruction loss used in both methods, along with the maximized target $p(X_{obs})$, theoretically indicates that our method is also capable of handling missing data in the MAR scenario, in the same way as VAE-based methods.

Thank you for your valuable reply and suggestions!

We have addressed your questions in the above answers and revision, and made changes according to your requirements. We look forward to receiving your feedback!

We thank the Reviewer for the valuable comment. We acknowledge that our paper is structured in multiple steps, which may have caused some difficulty in understanding. To address this, we offer the following intuitive explanation to provide a clearer understanding of our method.

In practice, the process within SHT-GNN can be understood as a form of data reconstruction. From an optimization perspective, this process is conceptually similar to Variational Autoencoders (VAEs), where the goal of the reconstruction step is to minimize the reconstruction error as part of the Evidence Lower Bound (ELBO). As is widely recognized, the training objective of a standard VAE is expressed as:

where $E_{q(z|x)}\log{p(x|z)}$ represents the $**reconstruction loss**$, and $D_{\text{KL}}[q(z|x)\parallel p(z)]$ is the $**regularization term**$. In our proposed SHT-GNN, the calculation and training process can be described as follows:

where $X_{\text{obs}}$ denotes the observed values, $Z_{O}^{\text{init}}$ and $Z_{F}^{\text{init}}$ represent the initial embedding matrices of all observation and feature nodes, respectively. $Z_{O}^{\text{L}}$ and $Z_{F}^{\text{L}}$ denote the embedding matrices of all observation and feature nodes after $L$ layers of forward computation in SHT-GNN. Subsequently, the training objective in SHT-GNN is expressed as:

$\text{That is to maximize}$: $E_{q(Z^L_{O},Z^L_F|X_{obs})}\log{p(X_{obs}|Z^L_{O},Z^L_{F})} - D_{\text{KL}}[p(Z^L_O,Z^L_F|X_{obs}) \parallel p(Z^L_O,Z^L_F)]$

where $E_{q(Z_{O},Z_F|X_{obs})}\log{p(X_{obs}|Z_{O}^L,Z_{F}^L)}$ represents the $**reconstruction loss**$.

In principle, the SHT-GNN and VAE-based methods share conceptual similarities. However, the key distinctions between the SHT-GNN and VAE-based approaches in the handling of longitudinal data are as follows.

$**1. Enhanced Representation Capability through Graph Neural Networks (GNNs)**$. In SHT-GNN, the encoding and decoding processes utilize graph neural networks (GNNs) instead of the multilayer perceptrons (MLPs) employed in standard VAEs. GNNs offer significantly stronger representation capabilities by leveraging the inherent data structure and facilitating temporal smoothing across non-independent observations, which is critical in longitudinal settings.

$**2. Distinct Regularization Term for Non-Independent Observations**$. In VAE-based methods, the regularization term for the latent space distribution, $D_{\text{KL}}[q(z|x) \parallel p(z)]$, quantifies the independent distributions of individual observations. By contrast, SHT-GNN incorporates a KL divergence term, $D_{\text{KL}}[p(Z^L_O,Z^L_F|X_{obs}) \parallel p(Z^L_O,Z^L_F)]$, which measures the joint distribution of all observation representations. To account for dependencies within longitudinal data, SHT-GNN employs multi-layer temporal smoothing and introduces MADGap, enabling improved representation of the joint latent space distribution.

In summary, the foundational intuition behind SHT-GNN lies in optimizing the reconstruction error during the data reconstruction process. Moreover, in the context of longitudinal data, SHT-GNN introduces temporal smoothing and MADGap to strike a balance between temporal consistency and the diversity of observation representations within the latent space.

We have realized that the lack of straightforward presentation and theoretical intuitiveness in our paper has caused difficulties for readers. Therefore, in the revision of our paper, we will include all the aforementioned theoretical understandings as a remark within the article.

Thank you for your valuable feedback! I will clarify which modules in our proposed method contribute the most to the performance of our model:

$**1. The inductive learning and representation capability of GNN as the framework’s foundation**$: Graph Neural Networks (GNNs) provide the foundational strength of this framework, particularly through their ability to effectively handle longitudinal data with irregular and inconsistent observation schedules. Unlike RNN-based methods, GNNs demonstrate superior representation capabilities, which underpin the accuracy of SHT-GNN in both missing data imputation and longitudinal data prediction. At each computational step, operations such as concatenation, aggregation, and updates adhere to the principles of GNNs, enabling the effective learning of observation embeddings. This robust representation capability is crucial for capturing complex relationships within irregular longitudinal datasets, ensuring the framework’s strong performance.

$**2. Temporal smoothing mechanism and MADGap as key components for data reconstruction in longitudinal data**$: The superior performance of this method arises from the multi-layer temporal smoothing mechanism, specifically designed for longitudinal data, and the incorporation of MADGap in the loss function. This combination facilitates a trainable trade-off between temporal smoothing and the preservation of the unique characteristics of observation representations, which is critical for effective data reconstruction in longitudinal datasets.

$**3. Mini-batch subject-wise sampling is not critical to the performance of SHT-GNN**$. The primary role of mini-batch subject-wise sampling is to address the memory storage challenges posed by massive datasets through the use of sharable trainable parameters. While it does not directly enhance the model’s accuracy or generalization capabilities, it demonstrates the feasibility of distributed training for the proposed model on large-scale datasets. For instance, this approach enables parallel training on multiple subgraphs using a gradient-sharing strategy, ensuring scalability for large-scale implementations. }

The above summaries are drawn from both theoretical and experimental analyses. Furthermore, several time-series datasets with relevant characteristics were identified to evaluate the performance of our method in missing data imputation and prediction tasks. And we also add some other state-of-the-art baselines in longitudinal data/time series imputation for comparison. The supplementary baselines are:

[1] Transformer: Ailing Zeng, Muxi Chen, Lei Zhang, and Qiang Xu. 2022. Are transformers effective for time series forecasting? arXiv preprint arXiv:2205.13504 (2022).

[2] GP-VAE: Fortuin V, Baranchuk D, Rätsch G, et al. Gp-vae: Deep probabilistic time series imputation[C]//International conference on artificial intelligence and statistics. PMLR, 2020: 1651-1661.

[3] SAITS: Du W, Côté D, Liu Y. Saits: Self-attention-based imputation for time series[J]. Expert Systems with Applications, 2023, 219: 119619.

[4] CTA: Wi H, Shin Y, Park N. Continuous-time Autoencoders for Regular and Irregular Time Series Imputation[C]//Proceedings of the 17th ACM International Conference on Web Search and Data Mining. 2024: 826-835.

For each dataset, we introduce missing data under an MCAR mechanism with a missing ratio of 0.3. We subsequently validated the imputation performance using the Mean Absolute Error (MAE) metric.

For each dataset, we also applied MCAR with a missing ratio of 0.5, and subsequently validated the missing data imputation results using MAE.

where the AirQuality dataset consists of a single-object time series with 15 features and over 9,000 records. The Electricity dataset comprises long-term time series data for 370 subjects, containing a total of over 140,256 records, each characterized by a single-dimensional attribute. The Energy dataset captures energy usage from various electrical appliances, representing 110 subjects with more than 10,000 observations and a total of 27 features.

On the AirQuality and Electricity datasets, the SHT-GNN method demonstrates inferior performance compared to RNN and VAE-based approaches. This indicates that the SHT-GNN model is not well-suited for long-term, single-object time series, which fall outside its intended application scenario.

However, SHT-GNN achieves state-of-the-art performance on the Energy dataset. Unlike the AirQuality and Electricity datasets, the Energy dataset involves time series with multiple subjects and multidimensional features. SHT-GNN leverages its ability to effectively borrow information across features, which significantly enhances its performance, particularly in the context of irregular longitudinal data imputation.

In addition, regarding the PhysioNet 2012 ICU dataset, we also completed experiments on this dataset.

It can be observed that our method demonstrates a significant advantage on clinical longitudinal follow-up data, including the ADNI and PhysioNet datasets.

Here, we present the performance of all the added baselines on the synthetic data:

On the ADNI dataset, their performance is as follows:

Concerning the missing mechanism:

In both the ADNI dataset and the PhysioNet 2012 ICU dataset, there is a significant presence of a MAR (Missing At Random) mechanism:

1. Existing studies indicate that the missing data in these datasets follows a MAR mechanism. The strong performance of our methods in these datasets reflects our approach's capability to effectively handle MAR mechanisms.

   [1] Campos S, Pizarro L, Valle C, et al. Evaluating imputation techniques for missing data in ADNI: a patient classification study[C]//Progress in Pattern Recognition, Image Analysis, Computer Vision, and Applications: 20th Iberoamerican Congress, CIARP 2015, Montevideo, Uruguay, November 9-12, 2015, Proceedings 20. Springer International Publishing, 2015: 3-10.

   [2] Lo R Y, Jagust W J, Aisen P, et al. Predicting missing biomarker data in a longitudinal study of Alzheimer disease[J]. Neurology, 2012, 78(18): 1376-1382.

2. Regarding experiments under the MNAR (Missing Not At Random) perspective, we will add these if time permits.

Single time series and weak feature correlation:

1. Single time series: The AirQuality dataset consists of a single object time series with over 9,000 records and 15 features.

2. Weak feature correlation: The Electricity dataset contains long-term time series data for 370 locations, totaling over 140,256 records, with each record including data collected at 370 locations at the current time point. There is weak correlation between different locations (i.e., different features).

From the results presented above, it is evident that our method does not perform notably on the AirQuality and Electricity datasets. Regarding the consideration of new baselines, we have added many cutting-edge baselines for longitudinal data imputation from the years 2022 to 2024, specifically including:

[1] Transformer: Ailing Zeng, Muxi Chen, Lei Zhang, and Qiang Xu. 2022. Are transformers effective for time series forecasting? arXiv preprint arXiv:2205.13504 (2022).

[2] GP-VAE: Fortuin V, Baranchuk D, Rätsch G, et al. Gp-vae: Deep probabilistic time series imputation[C]//International conference on artificial intelligence and statistics. PMLR, 2020: 1651-1661.

[3] SAITS: Du W, Côté D, Liu Y. Saits: Self-attention-based imputation for time series[J]. Expert Systems with Applications, 2023, 219: 119619.

[4] CTA: Wi H, Shin Y, Park N. Continuous-time Autoencoders for Regular and Irregular Time Series Imputation[C]//Proceedings of the 17th ACM International Conference on Web Search and Data Mining. 2024: 826-835.

Problems about large tables

Finally, regarding the final presentation style, thank you for your suggestion! We will change the large tables to graphical form in the final revision.

Seems like you have added some additional experiments to address my other comments. However, it is unclear to me based on the additional experiments if you have addressed these:

1. Varying degrees of those types of missingness, to empirically evaluate against other methods?
2. Single-time series and weak feature correlations, to show where the methods break down?

Also, in general, I find large tables incomprehensible, and prefer graphics. See various blog posts from Andrew Gelman on the topic.

In the above responses, we clarify the novelty of our work (What did we create? What did we apply? What did we modify?). We have uploaded the revision that includes all of these clarifications. Specifically, we have added Section 5.4, which explains in which situations our method excels, what types of data it is applicable to, and how it adapts to different missing data types.

To support these clarifications, we have also included the results of all additional baselines on the newly added datasets in the added appendix. In the appendix, we present the experiments using graphics as you suggested, rather than large tables. Finally, we have included a theoretical understanding of the method in the appendix as well.

Once again, we thank you for your patience and feedback!! We hope that we have addressed all of your concerns.

Thank you for your valuable reply and suggestions!

We have addressed your questions in the above answers and revision, and made changes according to your requirements. We look forward to receiving your feedback!

We thank the Reviewer for the insightful comment. In practice, the process within SHT-GNN can be understood as a form of data reconstruction. From an optimization perspective, this process is conceptually similar to Variational Autoencoders (VAEs), where the goal of the reconstruction step is to minimize the reconstruction error as part of the Evidence Lower Bound (ELBO). As is widely recognized, the training objective of a standard VAE is expressed as:

where $E_{q(z|x)}\log{p(x|z)}$ represents the $**reconstruction loss**$, and $D_{\text{KL}}[q(z|x)\parallel p(z)]$ is the $**regularization term**$. In our proposed SHT-GNN, the calculation and training process can be described as follows:

where $X_{\text{obs}}$ denotes the observed values, $Z_{O}^{\text{init}}$ and $Z_{F}^{\text{init}}$ represent the initial embedding matrices of all observation and feature nodes, respectively. $Z_{O}^{\text{L}}$ and $Z_{F}^{\text{L}}$ denote the embedding matrices of all observation and feature nodes after $L$ layers of forward computation in SHT-GNN. Subsequently, the training objective in SHT-GNN is expressed as:

$\text{That is to maximize}$: $E_{q(Z^L_{O},Z^L_F|X_{obs})}\log{p(X_{obs}|Z^L_{O},Z^L_{F})} - D_{\text{KL}}[p(Z^L_O,Z^L_F|X_{obs}) \parallel p(Z^L_O,Z^L_F)]$

where $E_{q(Z_{O},Z_F|X_{obs})}\log{p(X_{obs}|Z_{O}^L,Z_{F}^L)}$ represents the $**reconstruction loss**$.

In principle, the SHT-GNN and VAE-based methods share conceptual similarities. However, the key distinctions between the SHT-GNN and VAE-based approaches in the handling of longitudinal data are as follows.

$**1. Enhanced Representation Capability through Graph Neural Networks (GNNs)**$. In SHT-GNN, the encoding and decoding processes utilize graph neural networks (GNNs) instead of the multilayer perceptrons (MLPs) employed in standard VAEs. GNNs offer significantly stronger representation capabilities by leveraging the inherent data structure and facilitating temporal smoothing across non-independent observations, which is critical in longitudinal settings.

$**2. Distinct Regularization Term for Non-Independent Observations**$. In VAE-based methods, the regularization term for the latent space distribution, $D_{\text{KL}}[q(z|x) \parallel p(z)]$, quantifies the independent distributions of individual observations. By contrast, SHT-GNN incorporates a KL divergence term, $D_{\text{KL}}[p(Z^L_O,Z^L_F|X_{obs}) \parallel p(Z^L_O,Z^L_F)]$, which measures the joint distribution of all observation representations. To account for dependencies within longitudinal data, SHT-GNN employs multi-layer temporal smoothing and introduces MADGap, enabling improved representation of the joint latent space distribution.

In summary, the foundational intuition behind SHT-GNN lies in optimizing the reconstruction error during the data reconstruction process. Moreover, in the context of longitudinal data, SHT-GNN introduces temporal smoothing and MADGap to strike a balance between temporal consistency and the diversity of observation representations within the latent space.

We have realized that the lack of straightforward presentation and theoretical intuitiveness in our paper has caused difficulties for readers. Therefore, in the revision of our paper, we will include all the aforementioned theoretical understandings as a remark within the article.

1. $**Ignored baseline: “On the unreasonable effectiveness of feature propagation in learning on graphs with missing node features”**$

Although this paper leverages Graph Neural Networks, it does not address longitudinal or dynamic time-series data. Instead, its primary focus is on static graph data, such as social interactions and traffic flows. While the study incorporates missing data imputation within a GNN framework, it is not designed for temporal data. Similar works include:

• [1] Um D, Park J, Park S, et al. Confidence-based feature imputation for graphs with partially known features [J]. arXiv preprint arXiv:2305.16618, 2023.

• [2] Telyatnikov L, Scardapane S. EGG-GAE: scalable graph neural networks for tabular data imputation [C]//International Conference on Artificial Intelligence and Statistics. PMLR, 2023: 2661-2676.

These studies also focus on missing data imputation in static graph or tabular datasets. However, they do not address the challenges of imputing longitudinal or time-series data, which is the core problem our work seeks to tackle.

Thanks the authors for the detailed response. I have a few additional questions, but overall I will maintain my current rating and recommendation.