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

## Abstract

Abstract In this paper, we propose a novel framework, the Sampling-guided Heterogeneous Graph Neural Network (SHT-GNN), to effectively tackle the challenge of missing data imputation in longitudinal studies. Unlike traditional methods, which often require extensive preprocessing to handle irregular or inconsistent missing data, our approach accommodates arbitrary missing data patterns while maintaining computational efficiency.
SHT-GNN models both observations and covariates as distinct node types, connecting observation nodes at successive time points through subject-specific longitudinal subnetworks, while covariate-observation interactions are represented by attributed edges within bipartite graphs. By leveraging subject-wise mini-batch sampling and a multi-layer temporal smoothing mechanism, SHT-GNN efficiently scales to large datasets, while effectively learning node representations and imputing missing data.
Extensive experiments on both synthetic and real-world datasets, including the Alzheimer’s Disease Neuroimaging Initiative (ADNI) dataset, demonstrate that SHT-GNN significantly outperforms existing imputation methods, even with high missing data rates (e.g., 80%). The empirical results highlight SHT-GNN’s robust imputation capabilities and superior performance, particularly in the context of complex, large-scale longitudinal data.

## 1 Introduction

With the advancement of data collection and analysis techniques, longitudinal data has gained increasing importance across various scientific fields, such as biomedical science, economics, and e-commerce (Pekarčík et al., ; Sadowski et al., ; Mundo et al., ). Observation schedules in longitudinal studies vary widely as some follow regular intervals for each subject while others follow irregular patterns. For example, clinical visits scheduled every two months represent a regular observation schedule, whereas follow-ups at 6 months, 1 year, and 2 years post-treatment is an irregular schedule. Additionally, observation schedules may be either consistent or inconsistent across different subjects—some studies may impose uniform intervals for all participants, while others allow for variability due to individual health conditions or other circumstances.

Imputation of missing values poses a major challenge in longitudinal data analysis, both from theoretical and practical standpoints (Daniels & Hogan, ; Little & Rubin, ). In fact, in complex longitudinal data analysis, such as those in neuroimaging or electronic health record studies, patient follow-up data is collected over time, often resulting in missing values across variables (cross-sectional missingness) and across time (longitudinal missingness). Ideally, as illustrated in Figure 1, each subject would have observations taken at consistent intervals, with an equal number of observations per subject, and each observation would contain fully observed covariates and response variable. However, when missing data arises at specific time points, this misalignment of observations across time can occur. Furthermore, the absence of entire observations at certain time steps can transform regular observation schedules into irregular ones and consistent schedules into inconsistent ones, significantly complicating longitudinal data analysis.

Given the frequent occurrence of missing measurements in studies with extensive data collection, addressing the missingness of both covariates and response variables is crucial for accurate predictions of the target response variable. For instance, in Alzheimer’s disease research, the prediction of key biomarkers often relies on incomplete datasets to assess disease progression. Despite methodological advancements, managing missing data in longitudinal studies presents several ongoing challenges. Firstly, how can we effectively model longitudinal data that is irregularly and inconsistently observed? Secondly, how can we design a cohesive forward pass process for data imputation that accommodates varied observation schedules across subjects? Thirdly, how can we accurately predict target response values in the presence of missing covariate data and seamlessly integrate this process into model training? Fourthly, in the context of large-scale longitudinal studies, how can we ensure the scalability of the missing data imputation method?

In this paper, we address the key challenges by introducing the Sampling-guided Heterogeneous Graph Neural Network (SHT-GNN). Unlike traditional imputation methods, our GNN leverages scalable modeling of complex data structures to effectively learn from irregular and inconsistent longitudinal observations. Our SHT-GNN incorporates three key innovations:

- 1.
We handle longitudinal data by sharing trainable parameters across sampled graphs constructed from subject-wise mini-batches. This sampling-guided process ensures scalability for massive longitudinal datasets.
- 2.
Subject-wise longitudinal subnetworks are constructed by connecting adjacent observations with directed edges, while covariate values are transformed into attributed edges, linking observation nodes with covariate nodes. This structure allows SHT-GNN to flexibly model longitudinal data under arbitrary missing data conditions.
- 3.
We introduce a temporal smoothing mechanism across observations for each subject using a multi-layer longitudinal subnetwork. The novel MADGap statistic controls the smoothness within the subnetwork, balancing temporal smoothing with the specificity of observation node representations.

We conducted extensive experiments on both real data and synthetic data. The results demonstrate that SHT-GNN consistently achieves state-of-the-art performance across different temporal characteristics and performs exceptionally well even under high missing rates for both covariates and response variable. Furthermore, we applied the proposed SHT-GNN model to the Alzheimer’s Disease Neuroimaging Initiative (ADNI) dataset to predict the critical biomarker $A\beta 42/40$. SHT-GNN’s strong performance in predicting $A\beta 42/40$ is confirmed through validation with ground truth data and in downstream analyses. Additionally, we conduct extensive ablation studies, demonstrating the longitudinal network’s ability to perform temporal smoothing through multi-layer longitudinal subnetworks and the critical role of MADGap in guaranteeing variance between observations.

Figure: Figure 1: The ideal format of longitudinal data without missing (left); The regularly and consistently observed longitudinal data with missing in covariate variables and response variable (middle). And the irregular and inconsistent observation schedule in longitudinal data due to missing data (right).
Refer to caption: extracted/5985358/Figure1.png

## 2 Related Work

### 2.1 Statistical imputation methods

Statistical methods for imputing missing values in longitudinal data are often derived from multiple imputation techniques originally developed for cross-sectional data. For example, the 3D-MICE method extends the MICE (Multivariate Imputation by Chained Equations) framework to account for both cross-sectional and longitudinal dependencies in data (Luo et al., ). Some methods like trajectory mean and last observation carried forward (LOCF) leverage the time-series nature for missing data imputation (Lane, ). However, these basic models often fail to capture the complex, nonlinear spatio-temporal dependencies inherent in longitudinal data and lack the scalability to manage both continuous and discrete covariates simultaneously. The multi-directional multivariate time series (MD-MTS) method seeks to address these limitations by integrating temporal and cross-sectional covariates into a unified imputation framework through extensive feature engineering (Xu et al., ). A notable limitation of MD-MTS is its reliance on manually crafted features can be labor-intensive and prone to errors. To overcome such challenges, the time-aware dual-cross-visit (TA-DualCV) method leverages both longitudinal dependencies and covariate correlations using Gibbs sampling for imputation (Gao et al., ). Nonetheless, when applied to large-scale longitudinal datasets with numerous repeated observations, those methods involving chained equations and Gibbs sampling can become computationally expensive.

### 2.2 Machine learning imputation methods

Machine learning methods for imputing missing values in longitudinal data have primarily utilized recurrent neural networks (RNNs) and their advanced variants, such as long short-term memory (LSTM) networks and gated recurrent units (GRUs). BRITS, for example, leverages bidirectional LSTMs to capture both longitudinal and cross-sectional dependencies by utilizing past and future trends (Cao et al., ). Similarly, NAOMI employs a recursive divide-and-conquer approach with bidirectional RNNs to handle missing data imputation (Woillez et al., ). CATSI enhances this by combining bidirectional LSTMs with a context vector to account for patient-specific temporal dynamics (Kazijevs & Samad, ). A significant limitation of all RNN-based methods, however, is their reliance on consistent time step lengths across different subjects, often requiring padding or truncation to achieve this uniformity. When applied to large datasets, where the number of observations can vary greatly between subjects, padding and truncation may lead to excessive and inefficient computations due to the presence of numerous invalid time steps. Beyond RNN-based approaches, the GP-VAE model integrates deep variational autoencoders with Gaussian processes to address missing data in time series. However, the computational complexity of Gaussian processes, which scales at $O(N^{3})$, introduces significant computational bottlenecks when processing large-scale longitudinal data (Fortuin et al., ).

### 2.3 Graph-based imputation methods

Numerous advanced graph-based methods have been developed for data imputation tasks (Zhang & Chen, ; Li et al., ; Zhang et al., ). However, these methods face limitations when handling mixed discrete and continuous features and often fail to capture the temporal dependencies inherent in longitudinal data. For example, RAINDROP employs a graph-guided network to handle irregularly observed time series with varying intervals, effectively addressing missing data imputation and outperforming other methods in downstream prediction tasks (Zhang et al., ). Similarly, TGNN4I, which integrates GNNs with gated recurrent units, has shown strong performance in imputing missing data for longitudinal graph data (Oskarsson et al., ).
Despite their strengths, both RAINDROP and TGNN4I rely on predefined node connections within graph datasets, limiting their applicability to longitudinal tabular data where such connections between subjects do not inherently exist. Additionally, GRAPE demonstrated success using a bipartite graph for feature imputation, but its edge size increases linearly with the number of observations, posing scalability challenges (You et al., ). IGRM (Zhong et al., ), an extension of GRAPE, constructs a friend network among sample nodes to capture similarities and improve imputation. However, IGRM requires continuous updating of a fully connected graph among all samples, making it computationally expensive and lacking scalability for large-scale longitudinal data.

## 3 Problem Formulation

In the longitudinal data structure depicted in Figure 1, the observed variables are divided into two categories: covariates and the response variable. Assume there are $n$ subjects, with the $k$-th subject having $n_{k}$ observations measured. The total number of observations across all subjects is $N=\sum_{k=1}^{n}n_{k}$. Typically, the set of observed covariates varies across observations in longitudinal studies. For simplicity, we denote the complete set of covariates that can potentially be observed as $X=\{x_{1},\ldots,x_{p}\}$, where $p$ represents the total number of covariates. Any missing values within this set are treated as missing data. The covariate data for all observations can thus be represented as a matrix $\mathbf{D}=({D}_{il})\in\mathbb{R}^{N\times p}$, where $l\in\{1,\ldots,p\}$ and $i=\sum_{k=1}^{j-1}n_{k}+m$ indexes the $m$-th observation of the $j$-th subject. Similarly, the response variable for all observations is denoted as $\mathbf{Y}=(Y_{i})\in\mathbb{R}^{N\times 1}$. Given the presence of missing data, we introduce a mask matrix for the covariate data $\mathbf{M}^{O}=({M}^{O}_{{i}l})\in\{0,1\}^{N\times p}$, where ${D}_{{i}l}$ is observed if ${M}^{O}_{{i}l}=1$. Likewise, the mask matrix for the response variable is denoted by $\mathbf{M}^{Y}=({M}^{Y}_{{i}})\in\{0,1\}^{N\times 1}$, where the value of ${Y}_{{i}}$ is observed if ${M}^{Y}_{{i}}=1$. Within this irregular or inconsistent longitudinal data structure as depicted in Figure 1, the goal is to predict the response variable ${Y}_{{i}}$ for all ${i}$ where ${M}^{Y}_{{i}}=0$, leveraging the available observation information despite missingness.

## 4 Method

We propose SHT-GNN to address the challenges of missing data in irregular or inconsistent longitudinal data mentioned in Section 1. SHT-GNN models the observations and covariates as nodes, connecting them with attributed edges that represent the observed covariate values. Additionally, it employs specially designed longitudinal subnetworks to perform temporal smoothing among observations within the same subject. Initially, observations in longitudinal data can be represented as a graph $\mathcal{G}_{L}(\mathcal{V}_{O},\mathcal{E}_{OO})$ with subject-wise longitudinal subnetworks, where the observation node set $\mathcal{V}_{O}=\{u_{1},u_{2},\ldots,u_{N}\}$ represents all observations. The observation nodes that belong to the same subject form their own independent longitudinal subnetwork. In each longitudinal subnetwork, observation nodes at adjacent time points are connected by directed edges. The directed edge set in $\mathcal{G}_{L}$ is denoted as $\mathcal{E}_{OO}=\left\{(u_{i},u_{i^{\prime}},e_{i\rightarrow{i^{\prime}}})|u_
{i},u_{i^{\prime}}\in\mathcal{V}_{O},S(u_{i})=S(u_{i^{\prime}}),u_{i}\prec u_{
i^{\prime}}\right\}$, where $S(u)$ denotes the subject that observation $u$ belong to, and $u_{i}\prec u_{i^{\prime}}$ represents $u_{i}$ is the direct predecessor of $u_{i^{\prime}}$ in terms of time.

In irregularly and inconsistently observed longitudinal observation data, each observation may have measured different covariates. SHT-GNN establishes a covariate node set $\mathcal{V}_{C}=\{v_{1},v_{2},\ldots,v_{p}\}$, and construct edges between observation nodes and covariate nodes to represent the observed covariate values. The longitudinal data matrix $\mathbf{D}\in\mathbb{R}^{N\times p}$ and the missing indicator matrix $\mathbf{M}^{O}\in\{0,1\}^{N\times p}$ can be represented as a bipartite graph $\mathcal{G}_{B}=(\mathcal{V},\mathcal{E}_{OC})$. The undirected edge set $\mathcal{E}_{OC}=\left\{(u_{i},v_{l},e_{u_{i}v_{l}})|u_{i}\in\mathcal{V}_{O},v
_{l}\in\mathcal{V}_{C},{M}^{O}_{il}=1\right\}$, where the edge feature $e_{u_{i}v_{l}}$ takes the value of the corresponding feature $e_{u_{i}v_{l}}={D}_{il}$. Then a longitudinal dataset can be represented as the union of one undirected bipartite graph and multiple subject-wise longitudinal subnetworks: $\mathcal{G}(\mathcal{V},\mathcal{E})=\mathcal{G}_{B}\cup\mathcal{G}_{L}=
\mathcal{G}(\mathcal{V}_{O}\cup\mathcal{V}_{C},\mathcal{E}_{OO}\cup\mathcal{E}
_{OC})$. After constructing SHT-GNN, the task of response variable prediction under missing covariate imputation can be represented as learning the mapping: ${Y}_{i}=[f(\mathcal{G})]_{i}$ by minimizing the difference between ${Y}_{i}$ and $\widehat{{Y}}_{i}$, for all $i$ where ${M}^{Y}_{i}=1$.

Figure: Figure 2: The flow chart for Sampling-guided Heterogeneous Graph Neural Network.
Refer to caption: extracted/5985358/Algorithm.jpg

### 4.1 Learning in Sampling-guided Heterogeneous Graph Neural Network

#### 4.1.1 Subject-wise Mini-Batch Sampling

It is evident that the spatial complexity of graph $\mathcal{G}(\mathcal{V},\mathcal{E})$ increases with the number of subjects. When working with longitudinal data, it is common to encounter a large number of repeated observations or subjects. When dealing with millions or even tens of millions of observations, directly training a huge graph composed of all subjects is impractical and inefficient (You et al., ; Zhong et al., ). SHT-GNN performs subject-wise mini-batch sampling at the beginning of each training phase. The corresponding graph is then constructed based on the sampled subjects for the current training phase. In SHT-GNN, the graphs associated with each training phase share the same trainable parameters.

Specifically, assume there are $n$ subjects, represented by the subject set $S=\{1,2,\ldots,n\}$, where the $k$-th subject corresponds to an observation set $O_{k}$. Across all observation sets, there are a total of $N$ observations. At the beginning of each training phase, $s$ subjects are randomly sampled from $n$ subjects to form a subject-wise mini-batch $S_{0}=\text{Sample}(S,s)$, yielding the corresponding observation-wish mini-batch ${O}_{s}=\bigcup_{k\in S_{0}}{O}_{k}$. The observations in ${O}_{s}$ are then employed to construct the graph $\mathcal{G}_{s}=(\mathcal{V}_{s},\mathcal{E}_{s})$, where $\mathcal{V}_{s}=\mathcal{V}_{O_{s}}\cup\mathcal{V}_{C}$, with $\mathcal{V}_{O_{s}}=\{u_{i}|i\in{O}_{s}\}$. The edge set $\mathcal{E}_{s}$ is constructed within $\mathcal{V}_{s}$ according to the definitions provided above. Assuming the trainable parameters in SHT-GNN are denoted by $\theta$, and the loss function under input graph $\mathcal{G}_{\text{input}}$ is $L(\mathcal{G}_{\text{input}},\theta)$. The parameter learning process on the sampled graph in each sampling phase can be expressed as:

$$ $\theta^{(t+1)}=\theta^{(t)}-\eta_{t}\nabla_{\theta}L_{\mathcal{S}_{0}}( \mathcal{G}_{s},\theta^{(t)}),$ (1) $$

where $L_{\mathcal{S}_{0}}(\mathcal{G}_{s},\theta^{(t)})$ represents the loss function obtained from the forward pass in graph $\mathcal{G}_{s}$, and $\eta_{t}$ represents the learning rate. Here the forward pass in $\mathcal{G}_{s}$ begins with inductive learning within a multi-layer bipartite graph to derive representations of observation and covariate nodes. Afterward, temporal smoothing is conducted within the longitudinal subnetwork of each subject.

#### 4.1.2 Inductive learning in multi-layer bipartite graph

In SHT-GNN, all the information from the observed data is derived from the attributed edges connecting observation nodes and covariate nodes. The representations of all nodes need to be inductively learned from these edges. Inspired by GraphSAGE, which is a variant of GNNs known for its strong inductive learning capabilities. SHT-GNN modifies GraphSAGE architecture by introducing edge embeddings in the bipartite graph.
At the $l$-th layer in $\mathcal{G}_{B}$, the message generation function takes the concatenation of the embedding of the source node $\mathbf{h}_{v}^{(l-1)}$ and the edge embedding $\mathbf{e}^{(l-1)}_{uv}$ as input:

$$ $\mathbf{m}_{v}^{(l)}=\text{A}\text{\small GG}_{l}\left(\sigma\left(\mathbf{P}^ {(l)}\cdot\text{C}\text{\small ONCAT}(\mathbf{h}_{u}^{(l-1)},\mathbf{e}_{uv}^{ (l-1)})\right)\mid\forall u\in\mathcal{N}(v,\mathcal{E}_{s})\right)~{}~{}~{}~{ }\forall v\in\mathcal{V}_{s},$ (2) $$

where $\text{A}\text{\small GG}_{l}$ is message aggregation function for all node, $\sigma$ is the non-linearity activation function, $\mathbf{P}^{(l)}$ is the trainable weight matrix and $\mathcal{N}$ is the node neighborhood function. Subsequently, all nodes update their embeddings based on the concatenation of aggregated messages and their local representations:

$$ $\mathbf{h}_{v}^{(l)}=\sigma\left(\mathbf{Q}^{(l)}\cdot\text{C}\text{\small ONCAT }(\mathbf{h}_{v}^{(l-1)},\mathbf{m}_{v}^{(l)})\right)~{}~{}~{}~{}~{}~{}~{} \forall v\in\mathcal{V}_{s},$ (3) $$

where $\mathbf{Q}^{(l)}$ is the trainable weight matrix. Subsequently, the edge embeddings are then updated based on the updated node embeddings at both ends of each edge:

$$ $\mathbf{e}_{uv}^{(l)}=\sigma\left(\mathbf{W}^{(l)}\cdot\text{C}\text{\small ONCAT }(\mathbf{e}_{uv}^{(l-1)},\mathbf{h}_{u}^{(l-1)},\mathbf{h}_{v}^{(l-1)})\right )~{}~{}~{}~{}~{}~{}~{}\forall\mathbf{e}_{uv}\in\mathcal{E}_{s},$ (4) $$

where $\mathbf{W}^{(l)}$ is the trainable weight matrix.

#### 4.1.3 Temporal smoothing in multi-layer longitudinal subnetworks

Temporal smoothing is a crucial technique in the imputation of longitudinal data because it leverages the temporal correlation within the observations. After $L$ layers of inductive learning in the bipartite graph, SHT-GNN innovatively performs temporal smoothing through $K$ layers of message passing and representation updates within the subject-wise longitudinal subnetworks. At the $k$-th layer of each subject-wise longitudinal subnetwork, the message passing function for observations $u_{i}$ to $u_{i^{\prime}}$ takes the source node embedding $\mathbf{h}_{u_{i}}^{(L+k)}$ and the edge weight $w_{u_{i}\rightarrow{u_{i^{\prime}}}}^{(L+k)}$ as input:

$$ $\mathbf{m}_{u_{i}\rightarrow{u_{i^{\prime}}}}^{(L+k)}=\mathbf{h}_{u_{i}}^{(L+k )}\cdot w_{u_{i}\rightarrow{u_{i^{\prime}}}}^{(L+k)},\ \text{where}~{}~{}S(u_{ i})=S(u_{i^{\prime}}),u_{i}\prec u_{i^{\prime}},\forall u_{i},u_{i^{\prime}} \in\mathcal{V}_{O_{s}},$ (5) $$

Figure: Figure 3: The variation of information in the representation of observation nodes during multi-layer message passing and representation updates for temporal smoothing within longitudinal subnetworks.
Refer to caption: extracted/5985358/G.png

in which $S(u)$ denotes the subject that observation $u$ belong to and $u_{i}\prec u_{i^{\prime}}$ represents that $u_{i}$ is the direct predecessor of $u_{i^{\prime}}$ in terms of time. For each pair of message passing described in (5), the target observation node $u_{i^{\prime}}$ updates its representation as follows:

$$ $\mathbf{h}_{u_{i}}^{(L+k+1)}=\sigma\left(\mathbf{U}^{(L+k)}\cdot\text{C}\text{ \small ONCAT}(\mathbf{h}_{u_{i^{\prime}}}^{(L+k)},\ \mathbf{m}_{u_{i} \rightarrow{u_{i^{\prime}}}}^{(L+k)})\right),$ (6) $$

where $\mathbf{U}^{L+k}$ is the trainable parameter matrix for embedding updates. In (5), the weight for edges connecting observation nodes needs to be calculated and updated during training. The edge weight is computed as follows:

$$ $w^{(L+k)}_{u_{i}\rightarrow{u_{i^{\prime}}}}={D_{u_{i}\rightarrow{u_{i^{\prime }}}}}\cdot{J_{{u_{i}\rightarrow{u_{i^{\prime}}}}}}\cdot{\text{C}\text{\small os }\left(\mathbf{h}^{(L+k)}_{u_{i}},\mathbf{h}^{(L+k)}_{u_{i^{\prime}}}\right)}.$ (7) $$

Here, $D_{u_{i}\rightarrow{u_{i^{\prime}}}}$ represents the time decay weight in SHT-GNN. When facing irregular observation schedules in longitudinal data, the time interval between observations $u_{i}$ and $u_{i^{\prime}}$ is not fixed. We employ exponential decay functions to compute the time decay weights in (7). We define
$D_{{u_{i}}\rightarrow u_{i^{\prime}}}=\exp\left({{-\frac{|T(u_{i})-T(u_{{i^{
\prime}}})|}{\Delta_{max}}}}\right)$
(8)
$\mbox{and}~{}~{}~{}~{}~{}J_{{u_{i}}\rightarrow u_{i^{\prime}}}=\gamma-\frac{|
\mathbb{A}\cap\mathbb{B}|}{|\mathbb{B}|}.$
(9)

In (8), $T(u)$ represents the time step associated with observation $u$, and $\Delta_{max}$ represents the longest time interval between adjacent observations for the current subject. We further introduce the Jaccard distance between the sets of observed covariates in the weight calculation for message passing. Suppose the covariate sets observed by $\mathbf{h}_{u_{i}}$ and $\mathbf{h}_{u_{i^{\prime}}}$ are denoted as $\mathbb{A}$ and $\mathbb{B}$, respectively. As shown in (9), a larger $J_{{u_{i}}\rightarrow u_{i^{\prime}}}$ implies that observation $u_{i}$ contains more covariates that are not observed in observation $u_{i^{\prime}}$, and $\gamma$ is a hyperparameter constant for preventing the weights from becoming zero. Therefore, $\mathbf{h}_{u_{i^{\prime}}}$ need to borrow more information from $\mathbf{h}_{u_{i}}$. Beyond the time decay weight and Jaccard distance, the strength of message passing for a specific observation $\mathbf{h}_{u_{i^{\prime}}}$ from its immediate predecessor $\mathbf{h}_{u_{i}}$ is intuitively determined by the similarity between the currently observed covariates of $\mathbf{h}_{u_{i}}$ and $\mathbf{h}_{u_{i^{\prime}}}$. SHT-GNN then introduces the cosine similarity between $\mathbf{h}_{u_{i}}$ and $\mathbf{h}_{u_{i^{\prime}}}$ in the edge weight as discribed in (7).

Intuitions on why longitudinal subnetworks work. In the SHT-GNN, multi-layer message passing and representation updates within longitudinal subnetworks enable observations from the same subject to leverage information from previous time steps. As shown in Figure 3, in a two-layer longitudinal subnetwork, assume there are time-ordered observations $u_{t_{1}},u_{t_{2}},u_{t_{3}}andu_{t_{4}}$. The message passing and representation update in the first layer allow ${u_{t_{m}}}$ to directly draw information from ${u_{t_{m-1}}}$. Subsequently, while the temporal smoothing in the second layer enable ${u_{t_{m}}}$ to indirectly capture information from ${u_{t_{m-2}}}$ via the representation of ${u_{t_{m-1}}}$. This iterative process can be extended with additional layers, allowing observation nodes to integrate information from progressively earlier time steps, thereby complete temporal smoothing in subject-wise longitudinal subnetworks.

#### 4.1.4 Covariate imputation in S HT-GNN

After conducting $L$ layers of inductive learning on the bipartite graph and $K$ layers of temporal smoothing within longitudinal subnetworks, edge-level predictions are made at the $(L+K)$-th layer:

$$ $\widehat{{D}}_{il}=\mathbf{O}_{\text{impute}}\left({\text{C}\text{\small ONCAT }(\mathbf{h}^{(L+K)}_{u_{i}},\mathbf{h}^{(L+K)}_{v_{l}})}\right)~{}~{}~{}~{} \forall u_{i}\in{V}_{O_{s}},\ v_{l}\in{V_{F}},$ $$

where $\mathbf{O}_{\text{impute}}$ is a multilayer perceptron (MLP). Here, $\widehat{{D}}_{il}$ represents the imputation output for the $l$-th covariate for the $i$-th observation.

#### 4.1.5 Response variable prediction

Finally, we complete the prediction of the response variable based on the imputed covariates:

$$ $\widehat{{Y}}_{i}=\mathbf{O}_{\text{predict}}\left({\text{C}\text{\small ONCAT }({\widehat{D}}_{i1},{\widehat{D}}_{i2},...,{\widehat{D}}_{ip})}\right)~{}~{}~ {}~{}~{}\forall u_{i}\in\mathcal{V}_{O_{s}},$ $$

where $\mathbf{O}_{\text{predict}}$ is a multilayer perceptron and $\widehat{{Y}}_{i}$ represents the predicted response.

#### 4.1.6 Loss function for S HT-GNN

In both the response variable prediction and covariate imputation tasks, the loss function takes the following form:

$$ $\text{L\text{\small{oss}}}=\text{M\text{\small{SE}}}-\lambda\cdot\text{M\text{ \small{ADGap}}}.$ $$

Here, ${\text{M}\text{\small SE}}=(\sum_{i=1}^{N}{M}^{Y}_{i})^{-1}\sum_{i=1}^{N}{M}^{
Y}_{i}\cdot(\widehat{{Y}}_{i}-{Y}_{i})^{2}$, where ${M}^{Y}_{i}$ is the missing indicator. MADGap (Mean Average Distance Gap) is a statistical measure used to quantify the degree of over-smoothing in GNNs (Chen et al., ). A large MADGap value indicates that the node receives more useful information than noise. In SHT-GNN, the multi-layer longitudinal subnetwork between observation nodes represents the process of temporal smoothing. However, multi-layer message passing can lead to over-smoothing, causing the embedding representations of different observations within the same subject to become overly similar. SHT-GNN addresses this by maximizing MADGap to mitigate over-smoothing in the GNN.
MADGap is defined for individual nodes as follows:

$$ $\text{M\text{\small{ADGap}}}=\text{M\text{\small{AD}}}_{\text{remote}}-\text{M \text{\small{AD}}}_{\text{neighbour}}.$ $$

In SHT-GNN, for the $k$-th subject’s subnetwork with $n_{k}$ time ordered observations $\{u_{1},u_{2},...,u_{n_{k}}\}$, MADGap is calculated as:

$$ $\text{M\text{\small{ADGap}}}_{k}=\frac{1}{n_{k}}\sum_{m=1}^{n_{k}}\textbf{1}_{ m>2}\cdot\left(\frac{1}{m-2}\sum_{{m^{\prime}}=1}^{m-2}{\text{C\text{\small{os }}}(\mathbf{h}_{u_{m^{\prime}}},\mathbf{h}_{u_{m}})-\text{C\text{\small{os}}}( \mathbf{h}_{u_{m}},\mathbf{h}_{u_{m-1}})}\right).$ $$

Here $\frac{1}{m-2}\sum_{m^{\prime}=1}^{m-2}\text{Cos}(\mathbf{h}_{u_{m}},\mathbf{h}
_{u_{m^{\prime}}})$ denotes the similarity between the representations of the $m$-th observation and its past ancestor observation nodes in the longitudinal subnetwork, where $\text{Cos}(\mathbf{h}_{u_{m}},\mathbf{h}_{u_{m^{\prime}}})$ represents the cosine similarity between $\mathbf{h}_{u_{m}}$ and $\mathbf{h}_{u_{m^{\prime}}}$. Then for all $n$ subjects:

$$ $\small\text{M\text{\small{ADGap}}}=\frac{1}{n}\sum^{n}_{k=1}\left[\frac{1}{n_{ k}}\sum_{m=1}^{n_{k}}\textbf{1}_{m>2}\cdot\left(\frac{1}{m-2}\sum_{{m^{\prime} }=1}^{m-2}{\text{C\text{\small{os}}}(\mathbf{h}_{u_{m^{\prime}}},\mathbf{h}_{u _{m}})-\text{C\text{\small{os}}}(\mathbf{h}_{u_{m}},\mathbf{h}_{u_{m-1}})} \right)\right].$ $$

## 5 Experiment

### 5.1 Baselines

We consider eight baseline methods as follows.
1. Mean (Huque et al., ): imputes missing values using the covariate-wise mean.
2. Copy-mean Last Observation Carried Forward (LOCF) (Jahangiri et al., ): first impute missing values using the LOCF method, then refine the results based on the population’s mean trajectory.
3. Multivariate Imputation by Chained Equations (MICE)(Van Buuren & Groothuis-Oudshoorn, ): employs multiple regressions to model each missing value conditioned on other observed covariate values.
4. 3D-MICE (Kazijevs & Samad, ): combines MICE and Gaussian processes for longitudinal data imputation.
5. Decision Tree (DT) (Xia et al., ): a commonly used statistical method that can handle reponse prediction under missing data.
6. GRAPE: handles feature imputation through graph representation learning. (You et al., )
7. CASTI (Yin et al., ): handles missing data in longitudinal data by employing bidirectional LSTM and MLPs.
8. IGRM (Zhong et al., ): enhances feature imputation by leveraging the similarity between observations in the GNN.

### 5.2 Experiment on Synthetic Data Simulated from Real Data

Real data introduction. To comprehensively evaluate the performance of various methods on longitudinal data with diverse temporal characteristics, we first conduct experiments using synthetic data simulated from real-world datasets. We selected the longitudinal behavior modeling (GLOBEM) dataset as the basis for our synthetic data. The GLOBEM dataset is an extensive, multi-year collection derived from mobile and wearable sensing technologies. It contains data from 497 unique subjects, with over 50,000 observations and more than 2,000 covariates, including phone usage, Bluetooth scans, physical activity, and sleep patterns. For simulating synthetic datasets, all covariate values are sampled directly from the original data. Since the GLOBEM dataset does not include a predefined response variable, we simulate response values based on the observed covariates to construct complete synthetic datasets.

Response Variable simulation based on real data. In line with common assumptions in longitudinal data studies, when simulating the response values for the $m$-th observation of the $k$-th subject, assume the response value $Y_{i}$, where $i=\sum_{k=1}^{j-1}n_{k}+m$, lies within the $t$-th temporal smoothing window. For simplicity of expression, denote it as $Y^{t}_{km}$, which is assumed to follow a normal distribution $N(\mu^{t}_{k},{\sigma^{t}_{k}}^{2})$. Then $\mu_{k}^{t}$ and ${\sigma_{k}^{t}}^{2}$ respectively represent the mean and variance of the response variable for the $t$-th temporal smoothing window for the $k$-th subject. In practice, suppose the $k$-th subject has $n_{k}$ observations within the time step set $T=\{1,\ldots,n_{k}\}$. We first divide $T$ into $W$ temporal smoothing windows. For the $m$-th observation, which occurs within the $t$-th temporal smoothing window, the response value ${Y}^{t}_{km}$ is modeled as follows:
${Y}^{t}_{km}=\mu^{t}_{k}+\epsilon_{km},$
where $\mu^{t}_{k}$ represents the mean response value for the $k$-th subject within the $t$-th temporal smoothing window, and $\epsilon_{km}$ denotes the random fluctuations for the $m$-th observation.

Response variable simulation execution.
In our experiment, we run multiple random trials under a set of fixed parameters, reporting the average performance of all methods across these trials. In each random trial, we first randomly select $p$ dimensions from the 2000 covariates in the GLOBEM dataset, denoted as $X=\{x_{1},x_{2},...,x_{p}\}$. We then simulate the response values for each observation based on a specified model $f(X)=f(x_{1},x_{2},...,x_{p})$, which incorporates both linear and nonlinear elements. The details about two different specified models $f(X)$ are provided in Appendix A.1. According to the longitudinal effect model ${Y}^{t}_{km}=\mu^{t}_{k}+\epsilon_{km}$ described above, we first impute the mean value of the response variable for observations belonging to the same subject and within the same temporal smoothing window. Next, random fluctuations $\epsilon_{km}$ are sampled from the distribution $N(0,\epsilon)$ and incorporated into the formula ${Y}^{t}_{km}=\mu^{t}_{k}+\epsilon_{km}$ to simulate the final response values.

Experimental Procedure.
After obtaining complete synthetic datasets, we first split all subjects into a training set and a test set according to ratio $r$. By extracting the observations of all subjects, we obtain the covariate matrix $\mathbf{D}^{\text{train}}\in\mathbb{R}^{N_{\text{train}}\times p}$ and response vector $\mathbf{Y}^{\text{train}}\in\mathbb{R}^{N_{\text{train}}}$ for training, as well as $\mathbf{D}^{\text{test}}\in\mathbb{R}^{N_{\text{test}}\times p}$ and $\mathbf{Y}^{\text{test}}\in\mathbb{R}^{N_{\text{test}}}$ for testing. We respectively generate missing indicator matrices $\mathbf{M}^{\text{train}}$ and $\mathbf{M}^{\text{test}}$ for $\mathbf{D}^{\text{train}}$ and $\mathbf{D}^{\text{test}}$ according to the missing ratio $r_{X}$, along with the missing indicator vector $\mathbf{V}^{\text{train}}$ for $\mathbf{Y}^{\text{train}}$ based on $r_{Y}$. For all methods, we employ $\{{D}_{il}^{\text{train}}|{M}^{\text{train}}_{il}=1\}$ and $\{{Y}_{i}^{\text{train}}|{V}^{\text{train}}_{i}=1\}$ as input for training. In the testing phase, we use $\{{D}_{il}^{\text{test}}|{M}^{\text{test}}_{il}=1\}$ as input to predict all missing response values $\{\widehat{{Y}}_{i}^{\text{test}}|{V}^{\text{test}}_{i}=0\}$. Finally, we evaluate the response variable using root mean square error (RMSE) between ${Y}_{i}^{\text{test}}$ and $\widehat{{Y}}_{i}^{\text{test}}$ for all $i$ under ${V}^{\text{test}}_{i}=0$.

SHT-GNN configurations. We train SHT-GNN for 20 sampling phases with a sampling size of 200. For each sampled graph, we run 1500 training epochs using the Adam optimizer with a learning rate of 0.001(Kingma & Ba, ). We employ a three-layer bipartite graph and two-layer longitudinal subnetworks for all subjects. We use the ReLU activation function as the non-linear activation function. The dimensions of both node embeddings and edge embeddings are set to 32. The message aggregation function $\text{A\small{GG}}_{l}$ is implemented as a mean pooling function MEAN(·). Both $\textbf{O}_{impute}$ and $\textbf{O}_{predict}$ are implemented as multi-layer perceptrons (MLP) with 32 hidden units. The $\lambda$ in loss function is set to 0.001.

Baseline implementation. For Dicision Tree, GRAPE and IGRM, we directly conduct response prediction under missing covariate matrix. For other baselines, as no end-to-end response prediction approach is available, we first perform covariate imputation using the baselines, followed by utilizing a Multilayer Perceptron (MLP) as the prediction model. To ensure a fair comparison, we apply the same dimensional settings for these MLPs as those used in SHT-GNN.

**Table 1: Performance comparison with different methods under varying temporal smoothing window sizes and covariate missing ratios. All RMSE values are 0.1 of the actual values.**
| Missing ratio | $r_{X}=0.3,r_{Y}=0.3$ | $r_{X}=0.3,r_{Y}=0.5$ |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- |
| Window size<br>Variance | $w=3$<br>$\sigma=0.1$ | $w=5$<br>$\sigma=0.15$ | $w=7$<br>$\sigma=0.2$ | $w=3$<br>$\sigma=0.1$ | $w=5$<br>$\sigma=0.15$ | $w=7$<br>$\sigma=0.2$ |
| Mean | 0.693$\pm$0.009 | 0.826$\pm$0.011 | 0.936$\pm$0.011 | 0.820$\pm$0.012 | 0.833$\pm$0.012 | 1.051$\pm$0.016 |
| LOCF | 0.786$\pm$0.021 | 0.813$\pm$0.017 | 0.903$\pm$0.028 | 0.772$\pm$0.024 | 0.787$\pm$0.018 | 0.920$\pm$0.021 |
| MICE | 0.724$\pm$0.051 | 0.851$\pm$0.042 | 0.978$\pm$0.038 | 0.825$\pm$0.051 | 0.863$\pm$0.052 | 0.920$\pm$0.041 |
| 3D MICE | 0.689$\pm$0.031 | 0.753$\pm$0.037 | 0.847$\pm$0.029 | 0.741$\pm$0.035 | 0.785$\pm$0.037 | 0.883$\pm$0.021 |
| DT | 0.653$\pm$0.012 | 0.768$\pm$0.020 | 0.871$\pm$0.025 | 0.771$\pm$0.017 | 0.791$\pm$0.027 | 0.881$\pm$0.019 |
| GRAPE | 0.671$\pm$0.013 | 0.786$\pm$0.020 | 0.865$\pm$0.034 | 0.765$\pm$0.013 | 0.799$\pm$0.027 | 0.935$\pm$0.037 |
| CATSI | 0.701$\pm$0.034 | 0.732$\pm$0.026 | 0.832$\pm$0.023 | 0.749$\pm$0.047 | 0.748$\pm$0.038 | 0.885$\pm$0.027 |
| IGRM | 0.682$\pm$0.015 | 0.768$\pm$0.011 | 0.874$\pm$0.013 | 0.786$\pm$0.034 | 0.831$\pm$0.020 | 0.928$\pm$0.012 |
| Our Method | 0.552$\pm$0.011 | 0.650$\pm$0.014 | 0.759$\pm$0.018 | 0.623$\pm$0.013 | 0.653$\pm$0.018 | 0.818$\pm$0.020 |

Results in medium-dimensional covariates and moderate missing ratios. Following the experimental procedure described earlier, we set the covariate dimension $p$ to 50. For the window size $w$ and the variance of random fluctuation $\sigma(\epsilon)$, we use the parameter combinations $\{w=3,\sigma(\epsilon)=0.1\}$, $\{w=5,\sigma(\epsilon)=0.1\}$, and $\{w=7,\sigma(\epsilon)=0.2\}$. For each configuration, we first split the subjects with a test ratio of $r=0.2$, then apply two missing ratio settings to the covariate matrix and response vector: $\{r_{X}=0.3,r_{Y}=0.3\}$ and $\{r_{X}=0.3,r_{Y}=0.5\}$. We run All methods for 5 random trials per setting, and report the average RMSE along with its standard deviation for the response prediction on the test set. As shown in Table 1, SHT-GNN outperforms all baselines across all settings, achieving an average reduction of 17.5% in prediction RMSE compared to the best baseline. Across all methods, performance noticeably decline as the temporal smoothing window size and variance increase, particularly for LOCF, Decision Tree, and IGRM. This indicates that the various temporal smoothing characteristics in longitudinal data pose significant challenges for these methods. By integrating inductive learning with temporal smoothing, SHT-GNN demonstrate stable and superior performance in these scenarios.

Results in high-dimensional covariates and high missing ratios. We also compare the performance of different methods under higher covariate dimensions and missing ratios. As shown in Table 4 (Appendix A.3), SHT-GNN consistently outperforms all baselines across different settings, achieving an average reduction of 18% in prediction RMSE compared to the best baseline.

### 5.3 Alzheimer’s Disease Neuroimaging Initiative study dataset

ADNI dataset introduction. We apply SHT-GNN to the real data from Alzheimer’s Disease Neuroimaging Initiative (ADNI) study. We propose our model to predict the CSF biomarker Amyloid beta 42/40 ($A\beta 42/40$), which has been demonstrated as a crucial biomaker in ADNI study. ADNI dataset containing 1,153 subjects and 10,033 observations. The covariate matrix has 83 dimensions, with a missing ratios of 0.32 for covariate matrix and a missing ratio of 0.83 for $A\beta 42/40$. More details about ADNI dataset can be found in Appendix A.5,

SHT-GNN configurations for ADNI dataset. For the ADNI dataset, the configurations of SHT-GNN are detailed in Appendix A.6.

**Table 2: Performance comparison in ADNI dataset. For all methods, we conduct 5-fold cross-validation and report the mean values and standard deviations of the results.**
| Method | RMSE | AUC | Accuracy |
| --- | --- | --- | --- |
| Mean | 0.112$\pm$0.002 | 0.671$\pm$0.009 | 0.682$\pm$0.008 |
| LOCF | 0.108$\pm$0.003 | 0.706$\pm$0.008 | 0.717$\pm$0.015 |
| MICE | 0.109$\pm$0.004 | 0.717$\pm$0.021 | 0.701$\pm$0.021 |
| 3D MICE | 0.103$\pm$0.004 | 0.731$\pm$0.017 | 0.721$\pm$0.018 |
| DT | 0.107$\pm$0.001 | 0.703$\pm$0.011 | 0.709$\pm$0.009 |
| GRAPE | 0.106$\pm$0.002 | 0.724$\pm$0.010 | 0.714$\pm$0.011 |
| CATSI | 0.104$\pm$0.003 | 0.713$\pm$0.017 | 0.708$\pm$0.013 |
| IGRM | 0.105$\pm$0.002 | 0.721$\pm$0.013 | 0.713$\pm$0.009 |
| Our Method | 0.090$\pm$0.001 | 0.829$\pm$0.012 | 0.782$\pm$0.011 |

Validation on observed response value. We use RMSE to validate the prediction accuracy. As shown in Table 2, our proposed SHT-GNN significantly outperforms other methods in predicting $A\beta 42/40$, with an mean 15.5% improvement in accuracy compared to the best baseline.

Validation on diagnostic labels. In the previous section, we evaluated the $A\beta 42/40$ prediction results based on the observed ground truth. However, most observations in ADNI do not have recorded $A\beta 42/40$ values, making it impossible to validate the predictions for these observations using RMSE. In the ADNI study, the ultimate goal is to predict disease progression, so we validate the predictions by employing them to classify the diagnostic labels. Each observation is associated with a binary diagnostic label: either AD(Alzheimer’s Disease) or no-AD. We combine the predicted $A\beta 42/40$ values with fully observed demographic features (40 genetic data dimensions, sex, and age) and build a logistic model for classification. The classification performance can reflect the performance of $A\beta 42/40$ predictions. As shown in Table 2, using SHT-GNN-predicted $A\beta 42/40$ combined with basic individual features yields an AUC of 82.92% for the no-AD vs. AD classification. Notably, 35% of subjects in the ADNI dataset exhibit label transitions across different observations, indicating that SHT-GNN effectively capturing observation variance while conducting temporal smoothing.

### 5.4 Ablation Study on the Layer Count of Longitudinal Subnetworks

**Table 3: Performance comparison of longitudinal subnetworks with different layer count across different window sizes. All RMSE values are 0.1 of the actual values.**
| Window size | $w=3$ | $w=5$ | $w=7$ |
| --- | --- | --- | --- |
| None | 0.673$\pm$0.019 | 0.757$\pm$0.015 | 0.933$\pm$0.018 |
| One-Layer | 0.621$\pm$0.012 | 0.739$\pm$0.015 | 0.889$\pm$0.017 |
| Two-Layer | 0.593$\pm$0.015 | 0.687$\pm$0.019 | 0.827$\pm$0.020 |
| Three-Layer | 0.552$\pm$0.011 | 0.659$\pm$0.014 | 0.763$\pm$0.018 |

As shown in Figure 3, multi-layer longitudinal subnetworks allow observation nodes to borrow information from previous time steps. To validate that longitudinal subnetworks can borrow longer-term historical information by stacking layers, we applied SHT-GNN with one, two and three layers’ longitudinal subnetworks under data simulated under temporal smoothing window sizes $w=3$, $w=5$, and $w=7$. In all cases, we set the subject test ratio $r$ to 0.2, along with the missing rate of both the covariate matrix and the response vector to 0.3. As shown in Table 3, as $w$ increases, SHT-GNN with three-layer subnetworks outperform the one-layer and two-layer versions by a larger margin. This suggests that additional layers in longitudinal subnetworks enhance performance in data with longer-span temporal smoothing, which demonstrating the multi-layer longitudinal networks designed for temporal smoothing are effectve. Additionally, we conduct an ablation study for the MADGap component in SHT-GNN. The results and more details can be found in Appendix A.4.

### 5.5 Scalability of S HT-GNN

Figure: Figure 4: The comparison of scalability under different observation sizes across GNN-based methods.
Refer to caption: extracted/5985358/scalability_plot.png

In longitudinal data, a large number of repeated observations often lead to massive data size, making scalability a critical concern for missing data imputation methods. Specifically, we compare the scalability of SHT-GNN with two cutting-edge GNN-based imputation methods: GRAPE and IGRM. For experimental simplicity, we conduct imputation on datasets with varying observation sizes of $N\times 10$ (where $N=50,500,1000,$ and $5000)$. In SHT-GNN, we fix the subsample size at 500 subjects, resulting in a total of 5000 observations. For all three methods, we report the memory consumption and forward computation time per epoch during training. The results show that both GRAPE and IGRM exhibit significant increases in memory usage and computation time as observation size grows. However, by fixing the sampled subject batch size, SHT-GNN achieves consistent memory usage and computation time per epoch, regardless of the observation size.

## 6 Conclusion

In this paper, we present SHT-GNN, a scalable and accurate framework for longitudinal data imputation. Our method combines a sampling-guided training policy with inductive learning, temporal smoothing, and a custom-designed loss function to address the challenges posed by irregular and inconsistent longitudinal missing data. Compared to state-of-the-art imputation techniques, SHT-GNN consistently improves response prediction across both synthetic and real-world datasets. Extensive experiments confirm the efficacy of our multi-layer longitudinal subnetwork for temporal smoothing, and the loss function with MADGap effectively guarantees observation variance, enhancing model performance.

## 7 Reproducibility Statement

In this paper, we used the GLOBEM and ADNI datasets for our experiments. For the GLOBEM dataset, access can be requested at https://physionet.org/content/globem/, which is divided into four groups of subjects. In this study, we use the original data without any preprocessing. For the ADNI dataset, access can be requested at https://ida.loni.usc.edu/. The ADNI study is divided into different phases, and we select those subjects with complete basic genetic observations in the ADNI 1, ADNI 2, and ADNI GO phases. In the supplementary materials, we provide the complete code for the proposed SHT-GNN. The original codes for other baselines can be found through the links to the referenced papers in the Section. In this study, all models are trained on a Windows 10 64-bit OS (version 19045) with 32GB of RAM (AMD Ryzen 7 4800H CPU @ 2.9GHz) and 4 NVIDIA GeForce RTX 2060 GPUs with Max-Q Design.

## Appendix A Appendix

### A.1 Response vector simulation

In the experiment involving medium-dimensional covariates and moderate missing ratios, we simulate response variables from a 50-dimensional covariate matrix for each observation. Specifically, we use a response simulation model that includes both linear and nonlinear components. The default model for simulating the response variables is as follows:

$$ $\displaystyle y=0.25\cdot x_{1}+2\cdot\left(\frac{\log(x_{2}+10)}{25}\right)^{ 2}-0.4\cdot x_{3}-0.15\cdot\left(x_{4}+5\cdot e^{-5\left(1.5-\log(x_{5})\right )^{2}/2}\right)$ $\displaystyle-0.25\cdot\log(x_{6}+1)+0.4\cdot x_{7}+0.021\cdot\sin(x_{8})+0.04 \cdot\sqrt{x_{9}}+0.1\cdot e^{x_{10}}$ $\displaystyle+0.05\cdot\log(x_{11}+1)+0.02\cdot\tan(x_{12})+0.015\cdot\cos(x_{ 13})$ $\displaystyle+0.07\cdot\log(x_{14}+1)+3.5\cdot\sqrt{x_{15}}+\epsilon$ $$

Each $x_{i}$ is randomly selected without repetition from all the covariates. Before simulating the response values based on the observed covariates, all covariate values are normalized using the MinMax Scaler (Rajaraman, ). The term $\epsilon$ represents random noise following a normal distribution $N(0,0.175)$. In the experiment involving high-dimensional covariates and moderate missing ratios, we simulate response variables from a 100-dimensional covariate matrix for each observation. We use a response simulation model with higher-order inputs and more complex expressions, as shown below:

$$ $\displaystyle y$ $\displaystyle=0.3\cdot\sqrt{x_{1}}-0.4\cdot x_{2}^{2}+0.15\cdot\log(x_{3}+10^{ -6})+0.2\cdot\exp(0.5\cdot x_{4})-0.1\cdot x_{5}$ $\displaystyle\quad+0.05\cdot\sin(2\pi\cdot x_{6})+0.25\cdot\log 1p(x_{7})-0.1 \cdot\cos(2\pi\cdot x_{8})+0.35\cdot\tan(\text{clip}(x_{9},-0.5,0.5))$ $\displaystyle\quad+0.05\cdot\arcsin(\text{clip}(x_{10},-1,1))+0.2\cdot x_{11}^ {3}-0.3\cdot\sqrt{x_{12}}+0.4\cdot\frac{\log(x_{13}+1)}{10}$ $\displaystyle\quad+0.15\cdot\sin(2\pi\cdot x_{14})-0.1\cdot\log 1p(x_{15})+0.1 \cdot\exp(x_{16})-0.05\cdot\log(x_{17}+1)$ $\displaystyle\quad+0.2\cdot x_{18}^{2}+0.3\cdot\cos(x_{19})-0.07\cdot\tan(x_{2 0})+0.05\cdot\frac{\sin(x_{21})}{x_{22}+1}$ $\displaystyle\quad+0.25\cdot\log 1p(x_{23})+0.15\cdot\arcsin(\text{clip}(x_{24 },-1,1))+0.1\cdot x_{25}^{3}-0.05\cdot\sqrt{x_{26}}$ $\displaystyle\quad+0.07\cdot\log(x_{27}+1)+0.2\cdot\frac{\tan(x_{28})}{1+x_{29 }^{2}}-0.1\cdot\exp(x_{30})$ $\displaystyle\quad+0.3\cdot\log(x_{31}+10)+0.25\cdot x_{32}+\epsilon.$ $$

Similarly, each $x_{i}$ is randomly selected without repetition from all the covariates. Before simulating the response values based on the observed covariates, all covariate values are normalized using the MinMax Scaler. The term $\epsilon$ represents random noise following a normal distribution $N(0,0.2)$.

### A.2 Details for baselines

Mean: For each covariate with missing, we fill in the missing values using the mean of the observed values for that covariate across all observations.

Copy-mean LOCF: Following the process described in (Jahangiri et al., ), missing values are initially imputed using the LOCF (Last Observation Carried Forward) method within each subject to provide an approximation. Next, the population’s mean trajectory is used to further refine these imputed values.

MICE: As outlined in (Van Buuren & Groothuis-Oudshoorn, ), MICE (Multiple Imputation by Chained Equations) performs multiple imputations by modeling each missing value conditioned on the non-missing values in the data. A maximum of 20 iterations is used during the imputation process.

3D-MICE: Following the procedure described in (Luo et al., ), MICE is configured to perform cross-sectional imputation with a maximum of 20 iterations. And Gaussian Process Regression is applied longitudinally to the time-indexed data for each feature. The Gaussian Process uses an RBF kernel combined with a constant kernel ($\text{Kernel}=\text{C}(1.0)\times\text{RBF}(1.0)$), and the predictions from the Gaussian Process are averaged with the MICE-imputed values to capture both temporal and cross-sectional patterns.

Decision Tree: Following the procedure described in (Twala et al., ), Decision tree can handle the missing values by using surrogate splits, which are alternative features that closely correlate with the primary splitting feature. And the constructed tree is used for continuous response prediction.

GRAPE: Following the setup in (You et al., ), GRAPE is trained for 20,000 epochs using the Adam optimizer with a learning rate of 0.001. We employ two GNN layers with 16 hidden units and ReLU activation. The $\text{A\small{GG}}_{l}$ function is implemented as a mean pooling function MEAN(·). Both the edge imputation and response prediction neural networks are implemented as linear layers.

CATSI: Following the setup in (Yin et al., ), CATSI is trained for 3,000 epochs using the Adam optimizer with a learning rate of 0.001. We employ the default settings for MLPs and LSTM models in CASTI to impute covariate values.

IGRM: Following the default settings described in (Zhong et al., ), IGRM employs three GraphSAGE layers with 64 hidden units for bipartite Graph Representation Learning (GRL), and one GraphSAGE layer for friend network GRL. The Adam optimizer with a learning rate of 0.001 and ReLU activation function is used. For initializing the friend network, observation nodes are randomly connected with $|U|$ edges to form the initial network structure, which is updated every 100 epochs during bipartite graph training.

### A.3 Results in high-dimensional covariates and high missing ratios

To further compare the performance of different methods under higher covariate dimensions and higher missing ratios, we set the covariate dimension $p$ to 100. The same settings for temporal smoothing window size $w$ and the variance of random fluctuation $\sigma(\epsilon)$ were applied as before. Additionally, two higher missing ratio settings were employed: $\{r_{X}=0.5,r_{Y}=0.5\}$ and $\{r_{X}=0.5,r_{Y}=0.7\}$. All methods were run for 5 random trials per setting, and the average RMSE of response prediction on the test set are recorded. As shown in Table 4, SHT-GNN consistently outperforms all baselines across all settings, achieving an average reduction of 18% in prediction RMSE compared to the best baseline.

**Table 4: Performance comparison with different methods under varying temporal smoothing windows and covariate missing ratios. All RMSE values are 0.1 of the actual values.**
| Missing ratio | $r_{X}=0.5,r_{Y}=0.5$ | $r_{X}=0.5,r_{Y}=0.7$ |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- |
| Window size<br>Variance | $w=3$<br>$\sigma=0.1$ | $w=5$<br>$\sigma=0.15$ | $w=7$<br>$\sigma=0.2$ | $w=3$<br>$\sigma=0.1$ | $w=5$<br>$\sigma=0.15$ | $w=7$<br>$\sigma=0.2$ |
| Mean | 0.803$\pm$0.015 | 0.907$\pm$0.022 | 1.054$\pm$0.023 | 0.908$\pm$0.019 | 0.893$\pm$0.020 | 0.975$\pm$0.027 |
| LOCF | 0.751$\pm$0.029 | 0.813$\pm$0.035 | 0.938$\pm$0.035 | 0.825$\pm$0.023 | 0.823$\pm$0.019 | 0.937$\pm$0.019 |
| MICE | 0.797$\pm$0.043 | 0.848$\pm$0.041 | 1.014$\pm$0.059 | 0.837$\pm$0.036 | 0.893$\pm$0.041 | 0.974$\pm$0.065 |
| 3D MICE | 0.745$\pm$0.021 | 0.773$\pm$0.047 | 0.948$\pm$0.045 | 0.795$\pm$0.036 | 0.793$\pm$0.046 | 0.902$\pm$0.046 |
| DT | 0.763$\pm$0.017 | 0.798$\pm$0.018 | 0.981$\pm$0.020 | 0.781$\pm$0.019 | 0.794$\pm$0.027 | 0.891$\pm$0.031 |
| GRAPE | 0.724$\pm$0.026 | 0.793$\pm$0.015 | 0.952$\pm$0.021 | 0.809$\pm$0.021 | 0.745$\pm$0.021 | 0.944$\pm$0.021 |
| CATSI | 0.831$\pm$0.029 | 0.725$\pm$0.021 | 0.945$\pm$0.041 | 0.791$\pm$0.031 | 0.755$\pm$0.041 | 0.903$\pm$0.031 |
| IGRM | 0.795$\pm$0.010 | 0.831$\pm$0.013 | 0.904$\pm$0.014 | 0.785$\pm$0.012 | 0.815$\pm$0.019 | 0.895$\pm$0.021 |
| Our Method | 0.632$\pm$0.010 | 0.672$\pm$0.017 | 0.821$\pm$0.014 | 0.658$\pm$0.021 | 0.641$\pm$0.013 | 0.819$\pm$0.020 |

### A.4 Ablation Study for M ADGap in S HT-GNN

In SHT-GNN, the multi-layer longitudinal subnetworks are designed for temporal smoothing. However, the degree of smoothing among observations for the same subject may vary across different longitudinal studies. The loss function incorporates MADGap to promote greater variance amonng observation representations for a given subject, enabling the model to effectively trade off between temporal smoothing and representation diversity. To evaluate the impact of MADGap, we test SHT-GNN with and without MADGap under different temporal smoothing window sizes. As shown in Table 5, the results shows that incorporating MADGap enhances SHT-GNN’s performance by an average of 6%. Moreover, as the degree of temporal smoothing increases, the design incorporating MADGap exhibited a more substantial advantage over the one without it. This highlights MADGap’s capacity to help SHT-GNN capture the unique characteristics of each observation during imputation.

**Table 5: Performance comparison of SHT-GNN with and without MADGap across different temporal smoothing windows and missing ratios. All RMSE values are 0.1 of the actual values.**
|  | SHT-GNN without MADGap | SHT-GNN with MADGap |  |  |  |
| --- | --- | --- | --- | --- | --- |
|  | $r_{X}=0.3$ | $r_{X}=0.5$ | $r_{X}=0.3$ | $r_{X}=0.5$ | Enhancement |
|  | $r_{Y}=0.3$ | $r_{Y}=0.5$ | $r_{Y}=0.3$ | $r_{Y}=0.5$ |  |
| Window size 3 | 0.591$\pm$0.011 | 0.691$\pm$0.013 | 0.552$\pm$0.011 | 0.632$\pm$0.010 | 9.68% |
| Window size 5 | 0.687$\pm$0.012 | 0.718$\pm$0.011 | 0.651$\pm$0.014 | 0.672$\pm$0.017 | 5.57% |
| Window size 7 | 0.801$\pm$0.011 | 0.865$\pm$0.011 | 0.769$\pm$0.018 | 0.821$\pm$0.014 | 4.91% |

### A.5 A DNI dataset introduction

We apply SHT-GNN to the real data from Alzheimer’s Disease Neuroimaging Initiative (ADNI) study. ADNI is a multi-centre longitudinal neuroimaging study with the aim of developing effective treatments that can slow or halt the progression of Alzheimer’s Disease (AD). The ADNI participants were followed prospectively, with follow-up time points at 3 months, 6 months, then every 6 months until up to 156 months. The ADNI study includes a wide range of clinical data such as cognitive assessments, magnetic resonance imaging (MRI) and cerebrospinal fluid (CSF) biomarkers. Numerous studies show that CSF biomarkers are strong indicators of AD progression, but collecting CSF requires invasive procedures like lumbar puncture, leading to high missing data rates. We propose the SHT-GNN model to predict the CSF biomarker Amyloid beta 42/40 ($A\beta 42/40$), which has been a key biomaker. ADNI dataset containing 1,153 subjects and 10,033 observations. The covariate matrix has 83 dimensions, with a missing ratios of 0.32 for covariate matrix and a missing ratio of 0.83 for $A\beta 42/40$.

### A.6 S HT-GNN configurations for ADNI dataset

We train SHT-GNN for 10 sampling phases with a sampling size of 200. For each sampled graph, we run 1500 training epochs using the Adam optimizer with a learning rate of 0.001. We employ a three-layer bipartite graph and two-layer longitudinal subnetworks for all subjects. We use the ReLU activation function as the non-linear activation function. The dimensions of both node embeddings and edge embeddings are set to 32. The message aggregation function $\text{A\small{GG}}_{l}$ is implemented as a mean pooling function MEAN(·). Both $\textbf{O}_{impute}$ and $\textbf{O}_{predict}$ are implemented as multi-layer perceptrons (MLP) with 32 hidden units. The $\lambda$ in loss function is set to 0.001.

## Peer Reviews

confidence: 3
contribution: 3
initial_rating: 6
presentation: 3
questions: - Is there any transformer based model for longitudinal data imputation? If there is work with sequence models like RNN and LSTM, I would think there should be transformer based models too. If so, they should also be baselines, and they could potentially mitigate the computation burden by parallel computing
- The response simulation model in Appendix 1 seems strange. Are you just combining different linear/non-linear functions together to make it complex or is there any intuition or real world model to guide the design?
- One suggestion is to change the abbreviation to SH-GNN
soundness: 3
strengths: - The paper is well written and easy to follow. Illustrations make it clear to understand the framework.
- The proposed framework, especially the way to construct the graph is novel.
- Experiments show promising results on the provided framework, especially on the real ADNI dataset, which could have significant impact on biomedical science.
summary: The paper proposes Sampling-guided Heterogeneous Graph Neural Network (SHT-GNN) for Longitudinal data imputation. SHT-GNN treats both observations and covariates as nodes, and connects observations at successive time points, while keep covariate-observation interactions in a bipartite graph. Three key component in this model: subject-wise mini-batch sampling to reduce computation burden, inductive learning from observed edges and temporal smoothing in longitudinal subnetworks. With the embeddings generated with GNN, the authors employs MLPs to get the imputed covariates and predict the response. The loss function also takes consideration of a term Mean Average Distance Gap to prevent over-smoothing of GNN. Extensive experiment results and ablation studies are provided with semi-simulated and real data to demonstrate the effectiveness of the proposed framework.
weaknesses: - The paper is mainly heuristic on proposing a framework and show experiment results, without theoretical understanding of why it works. Although I do not think this is a deal breaker, it would be nice if the authors could provide more intuitions on why they design certain module the way it is. For example, it is unclear whether the way to get the weight with equation (7) is the best option, are there any alternative ways? How stable it is in practice and how to set $\gamma$ in practice?
- For covariate imputation, it seems to me one important baseline is missing, which is to do feature propagation with the proposed graph, similar to the way label propagation works. A recent work is [1]

[1] Rossi, Emanuele, et al. "On the unreasonable effectiveness of feature propagation in learning on graphs with missing node features." Learning on graphs conference. PMLR, 2022.

---

confidence: 3
contribution: 2
initial_rating: 5
presentation: 2
questions: 1. The method contains many steps, I did not follow them all, to be completely honest, there is a lot of new notation and terminology.  It would be great to simply the exposition of the method, if possible.

2. Of all the fancy things included in the method, including temporal smoothing, subject-wise mini-batch sampling, and concatening lots of stuff, experiments that illustrate which of those steps matters most, and how much, would be much more informative.

3. A simulation setting in which something trivial and straightforward outperforms your method, would be instructive.
soundness: 2
strengths: The paper focuses on a timely topic, dealing with longitudinal health data.  Moreover, they apply it to a real-world dataset of interest, ADNI. The developed a method that, at least in theory, overcomes limitations of existing tools.
summary: This paper introduces a method to impute data in multivariate longitudinal data.
weaknesses: There are 4-5 pages of methods.  Many steps were taken, each with its own few decisions.  It was quite difficult for me to ascertain which steps were important, and how much each step mattered. The ablation study did not help me much, I was not sure precisely whether this was on the real data or the synthetic data. 

The paper only shows experiments in which the SH-GNN method outperforms everything else, for all possible parameter values of the simulation, and also the real data.  That is a bit fishy.  No method outperforms all other methods in all other scenarios.  Especially one that runs faster and takes less memory.  I understand that there are likely no 'benchmark data' for this kind of problem.  And yet, my interests when I read a paper are to get very clear when I *should* use a method, and when I should not.  I get no real insight into that from this manuscript.  The claim seems to be that I should use it when I have multi-subject longitudinal data, with cross-sectional and cross-temporal missingness.  But, there are certainly conditions on the distribution and missingness under which other things are better, both in theory and practice.  What are the conditions?

I find it a bit weird that decision tree was included, rather than random forest, or some other decision forest based approach, since forests are known to be far superior to trees.

---

confidence: 2
contribution: 3
initial_rating: 6
presentation: 3
questions: 1. In what contexts are longitudinal observations inconsistent? Can your method handle noise in longitudinal measurements?

2. Does the temporal smoothing method work for longitudinal data with only a few measurements per subject?
soundness: 3
strengths: The topic of imputation for longitudinal data is significant, and the proposed solution is appealing.
summary: The goal of this paper is to present a scalable imputation method for handling missing data in longitudinal studies. A sampling-guided heterogeneous graph neural network (SHT-GNN) is proposed and it is both scalable and capable of learning effectively from irregular and inconsistent longitudinal observations. The proposed method is evaluated against eight baseline methods on synthetic data generated from real data, and it outperforms all baseline methods.

Missing data is a common issue in longitudinal studies, making this an important problem to address. However, imputation methods for longitudinal covariates have been extensively studied, and a key method by Yao et al. (2005) was overlooked and not included in the comparison study. Although this method was designed for Gaussian data, it performs well for non-Gaussian data in practice. It would be valuable to see how this method compares to SHT-GNN.

The authors claim that their method can accommodate arbitrary missing data. However, this may be an overstatement, as it seems likely that their approach is effective only for data missing completely at random. If the method is indeed applicable to missing-at-random or informatively missing schemes, this should be clarified.

Reference: 

Yao, Müller and Wang (2005). Functional data analysis for sparse longitudinal data.
Journal of the American Statistical Association.
weaknesses: An important method was omitted from the comparison study.

---

confidence: 4
contribution: 3
initial_rating: 5
presentation: 2
questions: Could the authors consider including more datasets, for instance as used previously for imputation of irregular data as referenced above? See in particular datasets and methodologies used here:
https://www.sciencedirect.com/science/article/pii/S0167947317300403
https://www.sciencedirect.com/science/article/pii/S0925231221003003
https://dl.acm.org/doi/10.1145/3616855.3635831
In particular, the second reference list many benchmark data sets (Table 3) as well as existing methodologies here used that are not currently compared against. Given the vast literature addressing imputation of irregular time-series data I think the paper needs to establish results much more extensively in terms of data and compared methods and the study is currently in its experimentation rather limited and in my eyes too limited. This makes it hard to judge the impact and utility of the approach compared to this rather large literature. I therefore strongly encourage the authors to include additional experimentation on well-established irregular temporal imputation datasets also considering additional well established methodologies for imputation in such irregular data (as given in the recent references above).

In summary, I think the authors’ methodology is sound and potentially indeed meritable and worth publication, but at this point I find this too unclear as the experimentation in terms of datasets and compared methodologies is too limited.
soundness: 3
strengths: The paper is well written , easy to follow, and with nice illustrations explaining the approach.

The methodology is sound and the included components including temporal smoothing and regularization for oversmoothing in GNN well motivated.

The approach have merits both in compute and ability to impute exploring temporal and multivariate/multimodal structure of data.
summary: The authors propose the Sampling-guided Heterogeneous Graph Neural Network (SHT-GNN) approach that explores GNN structures to learn to impute longitudinal irregular sampled multivariate/modal data. The approach relies on GNNs for linking observations through time exploring a subjects minibatched sampling stategy for scalable inference with edge weights defined in terms of temporal smoothing. The latter defined using the product of an exponential decay in time, overlap in missing pattern and cosine similarity of the associated nodes embeddings. To prevent overfitting the Mean Average Distance Gap is used as regularization during training. The approach is contrasted simple to more advanced approaches for irregular data imputation including some existing GNN based imputation methodologies on a synthetic (based on simulated response variables from a real GLOBEM dataset) and real longitudinal dataset (ADNI) finding superior imputation performance than the compared baselines as well as improved computational efficiency when compared to GRAPE and iGRM.
weaknesses: Whereas the approach is sound and the experimental comparison reasonable, it is somewhat limited. I.e. two datasets of which one is with simulated responses.

Furthermore, the literature on imputation of time-series data is vast and in this space, the paper covers only some of these works, see for instance:
https://paperswithcode.com/task/imputation

Arguably many of these listed procedures do not handle irregular sampled time-series, however, there are still many relevant methods and datasets for time-series imputation of irregular data that it would be highly interesting and relevant to see the performance of the current methodology against beyond the considered baselines. In particular the authors should consult Table 3 in the following recent survey of imputation of irregular time-series data for methods and datasets:
https://www.sciencedirect.com/science/article/pii/S0925231221003003
as well as the recently published article on imputation of irregular time-series data using autoencoders:
https://dl.acm.org/doi/10.1145/3616855.3635831

In particular, health care data such as MIMIC-III and the Physionet 2012 ICU challenge could here be relevant to consider as previously used and potentially also the other ICU datasets here listed. 

Methodology-wise it would be interesting to compare against the recent autoencoder framework of:
https://dl.acm.org/doi/pdf/10.1145/3616855.3635831
as well as, BRITS, GP-VAE, SAITS also here compared against.

Whereas the approach is sound and I believe with merits the experimentation is currently too limited to fully see how meritable the approach is and its impact upon this very large body of existing literature.


Minor:
Dicision Tree -> Decision Tree

---

confidence: 4
contribution: 2
initial_rating: 5
presentation: 3
questions: Q1: How to apply these tabular data imputation like GRAPE and IGRM for longitudinal data.

Q2: How to apply the time-series data imputation method for multiple observations in longitudinal data.

Q3: How do the state-of-the-art tabular imputation methods and time-series imputation methods perform on the longitudinal data?
soundness: 2
strengths: S1: The introduction of the method is thoroughly detailed, with complete and comprehensive mathematical formulas provided for each step.

S2: The motivation behind each step is clearly explained and sufficient.
summary: This manuscript proposes a new framework, SHT-GNN, for longitudinal data imputation based on GNN networks. 
It uses the bipartite graph to model the relationship between the covariates and observations and uses a direct graph to model the sequential relations of observations. 
Modeling the imputation as link prediction in the edge utilizes the graph to aggregate the information for an accurate imputation.
Experiments on both synthetic data and one realistic dataset validate its high performance.
weaknesses: W1: The novelty of the proposed method is limited.

W1.1. The method utilizes a bipartite graph to model the relationship between observations and covariates, addressing the task of imputation through link prediction. While similar approaches have been widely used for tabular imputation, as seen in models like GRAPE and IGRM, it appears that the author has merely adapted these methods for longitudinal data, which is essentially not very different from standard tabular data. GRAPE and IGRM also utilize the bipartite graph with U to denote the instance and V to denote the feature. The method proposed in this paper also follows this pipeline besides a directed graph used for instances to model the temporal relationship.

W1.2. It utilizes sampling methods for training the graph neural network, a widely used technique for GNNs. However, the sampling strategy employed here does not appear to be specifically designed for longitudinal data.

W1.3. A temporal smoothing module is proposed to connect different observations; however, its effectiveness has not been studied through ablation analysis. Although some experimental results are provided in section 5.4, these experimental results are conducted on synthetic datasets, which already imply strong temporal relationships. In real-world datasets, this relationship may not exist. 

=======

W2: The presentations of the work can be enhanced.

W2.1. The figures (Figures 1, 2, 3) are not vector images, causing them to become blurry when zoomed in.

W2.2. In Equation 9, the Jaccard distance (or similarity) is not correctly defined according to the standard formula, so it is inappropriate to call it the “Jaccard distance.” In the Jaccard definition, the denominator should be the union of sets A and B. Referring to it as a “Modified Jaccard Distance” may be more appropriate. You can either correct the formula or explicitly state and justify their modified version of the Jaccard distance.

W2.3: In Equation 6, the left side seems to be $h_{u_{i\prime}}$ rather than $h_{u_{i}}$.

W2.3: There are some typos, Row 223 typo, “observation-wish” -> “observation-wise”; Row 303 typo, “and” inside the formula; Row 416 typo, “dicision tree” -> “decision tree”.

======

W3: Some experiments can be enhanced to better highlight its performance.

W3.1: The simulation process in the generation of the datasets is complex and lacks motivation. As provided in the appendix, a 50-dimensional covariate matrix is constructed for each observation. They use a fixed equation with a lot of fixed numbers (hyperparameters). How these numbers or hyper-parameters are selected and while the usage of this fixed equation needs further explanation or citation.
  
W3.2: Some default settings of baselines are changed. For example, in the hidden dimension of the GRAPE, it is 64 by default in its original setting but the authors set it to 16 in the experiments (Line 701)

W3.3: The authors adopt two sets of parameters for different datasets, which may affect the model's adaptability and performance across varying data distributions.

W3.4: The authors include some imputation methods for tabular data (e.g., GRAPE and IGRM) and time series (e.g., CASTI) but they are not SOTA. Some more recent imputation methods such as ReMasker for tabular data [1] and [2] for time series are proposed. Coud you please why these specific baselines were chosen, or to include comparisons with the more recent methods like in [1] and [2].

[1]: ReMasker: Imputing Tabular Data with Masked Autoencoding

[2]: Mining of Switching Sparse Networks for Missing Value Imputation in Multivariate Time Series

Meta-review: This work proposed a GNN model for longitudinal data imputation, where the observations, covariates, temporal and subject information are jointly considered in a heterogeneous graph. The reviewers agree that the paper addressed an important problem in dealing with longitudinal data, and that the revision has improved the empirical evaluation. The reviewers share the sentiment that the overall novelty can be better clarified to position this work in the literature. Despite the authors made efforts in explaining the mechanism based on graph VAE, the current writing lacks motivations and justifications from a formal perspective, and on how the authors' approach addresses limitations of existing methods. The generalizability to other types of missing data and more general data types should be thoroughly explained.

Decision: Reject