¹

1 Introduction↩︎

In clinical care, key indicators, e.g., Hemoglobin A1c (HbA1c) for diabetes management [1], are crucial for disease monitoring and treatment adjustment [2], yet patient visits are irregular, and labs are often not ordered at every visit, resulting in pervasive missingness in longitudinal records. Recent years have seen growing interest in clinical data imputation, with methods ranging from statistical models [3]–[5] to deep generative models and representation learning [6]–[9]. Despite this progress, most approaches emphasize overall accuracy and provide limited guarantees on error control in high-stakes use, where unreliable results can directly mislead diagnosis and treatment adjustment.

Clinical deployment, therefore, raises a question beyond overall accuracy: which imputed results are reliable enough to support clinical decision-making, and how can we provide explicit control over imputation errors? This motivates reliable clinical imputation, which augments imputation with uncertainty or reliability assessment, thereby producing imputations with quantitative error control. Related quality-control frameworks have been studied in high-stakes prediction, including Bayesian methods [10], [11], imprecise probability [12], and conformal prediction [13]. Several works have begun to apply these ideas in clinical settings [14]–[17], such as conformalized models for genomic medicine [18] and SAFER [19], which uses conformal prediction to provide statistical false discovery rate (FDR) control for treatment prediction. However, risk-controlled clinical imputation remains underexplored, particularly for irregular and sparsely observed patient records [20], [21].

Building on this gap, we focus on two challenges in reliable clinical data imputation. (i) Modeling irregular, sparse, and heterogeneous longitudinal records. Real-world patient visits usually occur at uneven intervals, and lab tests are ordered opportunistically, yielding highly sparse and heterogeneous measurement patterns. As a result, intra-patient temporal evidence is often limited, while informative signals may instead come from cross-patient analogies. This calls for an imputation model that can exploit both the limited temporal trajectory and relational structure across patients under severe missingness. (ii) Reliable imputation with statistical error control. In high-stakes settings, the goal is not only to minimize average error, but to determine which imputations are reliable enough to act on, and provide a quantitative guarantee that clinically unacceptable errors are rare among the released results. This requires uncertainty quantification for each imputed entry and a statistically grounded selection procedure that can control the error rate when releasing many imputations.

To address these challenges, we propose SafeImpute, a reliable imputation framework for irregular and sparse clinical records that couples a strong event graph imputer with alse discovery rate (FDR) control. Specifically, to address the irregular and sparse nature of clinical visits and measurements, we construct an event graph with temporal edges that capture intra-patient trajectories and trend-aware value edges that connect clinically similar events to capture inter-patient analogies. On this event graph, we learn event representations with a two-relation GNN that performs relation-specific message passing and fuses temporal and value signals via adaptive gating. To ensure reliability in high-stakes settings, SafeImpute computes a proxy risk score for each imputed entry and calibrates it into conformal p-values, then applies the Benjamini–Hochberg procedure to select a subset of imputations while controlling the FDR of clinically unacceptable errors at a user-specified tolerance level. Together, SafeImpute produces accurate imputations and selectively releases them with statistical error control.

The major contributions of this paper are summarized as follows:

• We formulate the problem of reliable clinical data imputation, where the goal is not only to achieve accurate imputation but also to selectively release reliable results with statistical control over clinically unacceptable errors.

• We propose SafeImpute, a reliable imputation framework that (i) represents irregular clinical visits as an event graph, where temporal edges capture intra-patient trajectories and similarity edges connect clinically related visits across patients, and learns imputations using a two-relation GNN with adaptive fusion; and (ii) provides reliability control by calibrating a proxy risk score into conformal p-values and applying the Benjamini–Hochberg procedure to control the false discovery rate among released imputations at a user-specified tolerance.

• Extensive experiments on Mayo Clinic dataset and the MIMIC-III/IV datasets demonstrate that SafeImpute achieves strong imputation accuracy while providing reliable error control.

2 Preliminary↩︎

2.1 Clinical Imputation↩︎

In many clinical decision-making scenarios, decisions are driven by a small set of decision-critical laboratory markers that are directly relevant for diagnosis and treatment planning, rather than by completing every missing measurement. Motivated by this, we focus on imputing a designated decision-critical target, such as HbA1c for diabetes monitoring and treatment adjustment. We consider irregular, multi-lab longitudinal clinical records, where each patient has an irregular set of clinical visits over time and only a subset of laboratory tests may be measured at each visit. We treat each visit as an event node that aggregates the laboratory tests recorded during that visit. Let \(\mathcal{V}\) denote the set of all event nodes across patients. For node \(i\in\mathcal{V}\), let \(p_i\) denote its patient identity and \(t_i\) its timestamp. Let \(\mathcal{L}\) be the set of laboratory variables together with the available attributes of the patient. We denote the value of variable \(\ell\in\mathcal{L}\) at node \(i\) by \(z_{i,\ell}\). Formally, given partially observed multi-lab inputs \(\{z_{i,\ell}\}_{\ell\in\mathcal{L}}\), our goal is to impute only the designated target \(y_{i}\) at event nodes where it is missing.

2.2 Conformal Selection↩︎

Conformal selection (CS) [21] is a model-agnostic framework for selecting a subset of test instances while controlling the false discovery rate (FDR) [22] in finite samples.

CS formulates one hypothesis test for each candidate test sample. Let \(\{(x_i,y_i)\}_{i=1}^n\) denote the calibration samples and \(\{x_{n+j}\}_{j=1}^m\) denote the test candidates. For each test candidate \(j\in[m]\), CS considers the hypotheses \[H_{0j}:~ y_{n+j}\le c \label{eq:cs95hypothesis}\tag{1}\] Rejecting \(H_{0j}\) indicates that the \(j\)-th test candidate is deemed to exceed the threshold \(c\).

2.2.0.1 Nonconformity scores and conformal p-values.

CS relies on a nonconformity measure \(J(\cdot,\cdot)\) that quantifies how atypical an observation is relative to the calibration data. For calibration samples, the nonconformity scores are \[J_i = J(x_i,y_i),\quad i=1,\dots,n.\] For a test candidate \(x_{n+j}\), since \(y_{n+j}\) is unobserved, CS replaces it by the threshold \(c\) and computes \[\widehat{J}_{n+j} = J(x_{n+j},c).\] The conformal p-value [21] for testing 1 is then computed via a rank-based comparison between \(\widehat{J}_{n+j}\) and \(\{J_i\}_{i=1}^n\). Intuitively, a smaller conformal p-value provides stronger evidence against \(H_{0j}\).

2.2.0.2 BH procedure and FDR control.

To determine the final selected subset, CS applies the Benjamini–Hochberg (BH) [22] procedure to the set of conformal p-values \(\{p_{n+j}\}_{j=1}^m\). Let \(p_{(1)}\le \cdots \le p_{(m)}\) be the sorted p-values and define \[k^\star = \max\left\{k\in[m]:~ p_{(k)} \le \frac{k}{m}\,\alpha \right\},\] with the convention that \(k^\star=0\) if the set is empty. \(\alpha\in(0,1)\) is the target FDR level. BH rejects all hypotheses with \(p_{n+j}\le p_{(k^\star)}\), yielding the selected set. Under the standard validity conditions for conformal p-values (e.g., exchangeability), this procedure guarantees finite-sample FDR control at level \(\alpha\).

3 Method↩︎

In this section, we present SafeImpute, including event-graph construction for irregular clinical records, two-relation GNN learning over temporal and trend-aware value edges, and the selection module that provides FDR-controlled selective imputation release.

3.1 Event Graph Construction↩︎

3.1.0.1 Event Node.

We treat each patient visit as an event node. To preserve as much event-level information as possible under irregular missingness, we encode both available lab values and missingness patterns in the node feature. Specifically, for each lab \(\ell\in\mathcal{L}\) at event node \(i\), we define an observation indicator \[m_{i,\ell}=\mathbb{I}\big[z_{i,\ell}\;\text{is observed}\big].\] We then represent each event by concatenating lab values with their indicators: \[\mathbf{x}_i=\Big[\, z_{i,1}, m_{i,1},\;z_{i,2}, m_{i,2},\;\ldots,\;z_{i,|\mathcal{L}|}, m_{i,|\mathcal{L}|},\;g_i \Big],\] where unobserved \(z_{i,\ell}\) are set to \(0\) only for feature construction. To mitigate scale differences across lab values, we normalize the lab values.

With event nodes, we construct an event graph \(\mathcal{G}=(\mathcal{V},\mathcal{E}_t\cup\mathcal{E}_v)\) with two types of edges that capture two complementary signals: intra-patient trajectories and inter-patient analogies. Specifically, we include (1) intra-patient temporal edges \(\mathcal{E}_t\) that connect consecutive events of the same patient, preserving local temporal continuity; and (2) inter-patient trend-aware value edges \(\mathcal{E}_v\) that connect events among patients using lab value-level and trend-level information, enabling information sharing between clinically similar event contexts. We describe the construction of these two edge types below.

3.1.0.2 Temporal Edges.

To preserve intra-patient trajectories, for each patient \(p\) we sort its event nodes by time. Let \(\text{next}(i)\) denote the immediate subsequent event of the same patient after node \(i\). We add a temporal edge between consecutive events if the time gap is within a threshold: \[\mathcal{E}_t = \left\{(i,j)\;\middle|\;p_i=p_j,\;j=\text{next}(i),\;0 < t_j-t_i \le \Delta_t \right\}.\]

3.1.0.3 Trend-aware Value Edges.

Temporal edges capture within-patient continuity but do not facilitate information sharing across patients. To model inter-patient analogies, we additionally connect events that exhibit similar lab levels and temporal trends.

For each node \(i\) and lab \(\ell\), let \(r_{i,\ell}\) denote the most recent prior event of patient \(p_i\) before \(t_i\) where lab \(\ell\) is observed. If such an event does not exist, we set \(r_{i,\ell}=\varnothing\). We define the per-lab trend as the rate of change from the last available observation: \[\delta_{i,\ell}= \begin{cases} \dfrac{z_{i,\ell}-z_{r_{i,\ell},\ell}}{(t_i-t_{r_{i,\ell}})+\epsilon}, & r(i,\ell)\neq\varnothing,\;t_i>t_{r_{i,\ell}},\\[6pt] \text{missing}, & \text{otherwise}. \end{cases}\]

For a pair of nodes \((i,j)\), the sets of commonly observed labs for values and trends are defined as: \[\mathcal{L}_{ij} = \left\{\ell\in\mathcal{L}\mid z_{i,\ell}\;\text{and}\;z_{j,\ell}\;\text{are observed}\right\},\]

\[\mathcal{L}^{\delta}_{ij} = \left\{\ell\in\mathcal{L}\mid \delta_{i,\ell}\;\text{and}\;\delta_{j,\ell}\;\text{are observed}\right\}.\] To avoid node similarity being biased by the number of commonly observed labs, we compute RMS-normalized Euclidean distances over the common dimensions as value and trend similarities: \[\begin{align} d_v(i,j) &= \sqrt{\frac{1}{|\mathcal{L}_{ij}|}\sum_{\ell\in\mathcal{L}_{ij}}\left(\hat{z}_{i,\ell}-\hat{z}_{j,\ell}\right)^2},\\ d_\delta(i,j) &= \sqrt{\frac{1}{|\mathcal{L}^{\delta}_{ij}|}\sum_{\ell\in\mathcal{L}^{\delta}_{ij}}\left(\hat{\delta}_{i,\ell}-\hat{\delta}_{j,\ell}\right)^2}, \end{align}\]

where \(\hat{z}_{i,\ell}\) and \(\hat{\delta}_{i,\ell}\) denote the normalized value and trend features, respectively. We connect \(i\) and \(j\) if both the value and trend distance criteria are satisfied: \[\mathcal{E}_v =\{(i,j)| d_v(i,j) \le \tau_v, d_\delta(i,j)\le \tau_\delta\}.\] To avoid isolated nodes under sparse observations, for any node \(i\) with \(\mathcal{L}^{\delta}_{ij}=\emptyset\) for all \(j\), we connect \(i\) to its nearest value-based neighbor.

3.2 Two-Relation Learning↩︎

To learn event representations that jointly capture intra-patient progression and inter-patient similarity, we propose a two-relation graph neural network (GNN) that performs separate message passing over temporal edges \(\mathcal{E}_t\) and value edges \(\mathcal{E}_v\), followed by adaptive fusion at each layer.

Let \(\mathbf{h}_i^{(0)}=\mathbf{x}_i\). At layer \(k\), we compute relation-wise messages for each node \(i\): \[\begin{align} \mathbf{m}^{(k)}_{t,i} &= \mathrm{GNN}_t^{(k)}(\mathbf{H}^{(k-1)};\mathcal{E}_t)_i, \\ \mathbf{m}^{(k)}_{v,i} &= \mathrm{GNN}_v^{(k)}(\mathbf{H}^{(k-1)};\mathcal{E}_v)_i, \end{align}\] where \(\mathrm{GNN}_t^{(k)}\) and \(\mathrm{GNN}_v^{(k)}\) are GCN layers with independent parameters: \[\mathrm{GNN}(\mathbf{H};\mathcal{E}) = \sigma\!\left(\tilde{\mathbf{D}}^{-\frac{1}{2}}\tilde{\mathbf{A}}\tilde{\mathbf{D}}^{-\frac{1}{2}}\mathbf{H}\mathbf{W}\right),\] where \(\tilde{\mathbf{A}}\) is the adjacency with self-loops and \(\tilde{\mathbf{D}}\) is the corresponding degree matrix, \(\mathbf{W}\) is the learnable weight matrix and \(\sigma(\cdot)\) is the activation function.

We fuse the two relation-specific messages using a learnable gate computed for each node at each layer: \[\alpha_i^{(k)} = \mathrm{sigmoid}\!\left( \mathbf{w}^{(k)\top} \big(\mathbf{m}^{(k)}_{t,i}+\mathbf{m}^{(k)}_{v,i}\big) \right),\] \[\mathbf{h}_i^{(k)} = \alpha_i^{(k)}\mathbf{m}^{(k)}_{t,i} + \big(1-\alpha_i^{(k)}\big)\mathbf{m}^{(k)}_{v,i}.\]

After the final layer \(K\), we apply a multi-output linear prediction head: \[\hat{\mathbf{z}}_i = \mathbf{W}_o \mathbf{h}_i^{(K)} + \mathbf{b}_o,\] where each output dimension corresponds to one laboratory variable. The prediction for the designated target measurement is denoted by \(\tilde{y}_i\), i.e., the target dimension of \(\hat{\mathbf{z}}_i\).

Let \(\Omega_{\mathrm{train}}\) be the set of training nodes with observed target labels \(y_i\). We train the main target prediction by minimizing \[\mathcal{L}_{\mathrm{target}} = \frac{1}{|\Omega_{\mathrm{train}}|} \sum_{i\in\Omega_{\mathrm{train}}} (\tilde{y}_i-y_i)^2.\]

Although the final task is to impute the designated target measurement, we further regularize the event representations by reconstructing randomly masked non-target laboratory values. At each epoch, we randomly mask a small fraction of observed non-target labs from the input and ask the model to reconstruct them from the learned event representation. Let \(\mathcal{M}_\ell\) be the set of masked node–lab pairs for lab \(\ell\). To reduce sensitivity to outlying lab values, we use the Huber reconstruction loss: \[\mathcal{L}_{\mathrm{aux}} = \frac{1}{\sum_{\ell \ne y} |\mathcal{M}_\ell|} \sum_{\ell \ne y} \sum_{(i,\ell)\in \mathcal{M}_\ell} \rho_{\eta}(\hat{z}_{i,\ell}-z_{i,\ell}),\] where \[\rho_{\eta}(e)= \begin{cases} \frac{1}{2}e^2, & |e|\le \eta,\\ \eta\left(|e|-\frac{1}{2}\eta\right), & |e|>\eta. \end{cases}\] The final training objective is \[\mathcal{L} = \mathcal{L}_{\mathrm{target}} + \lambda \mathcal{L}_{\mathrm{aux}}.\]

3.3 Conformal Selection and FDR Control↩︎

In clinical settings, imputed laboratory values may directly influence diagnosis, prognosis assessment, and treatment planning. We therefore apply conformal selection with false discovery rate (FDR) control to retain only reliable imputations: among the released imputations, the expected fraction of clinically unacceptable errors is controlled at a user-specified level.

3.3.0.1 Risk hypothesis.

Let \(\hat{y}_i\) be the model prediction at node \(i\) and \(y_i\) be the ground-truth label (available only on calibration). We define the true error on labeled points as \[\kappa(i) = |y_i-\hat{y}_i|.\] Given a user-specified tolerance \(\delta>0\) for the target variable, we call node \(i\) risky if \(\kappa(i)\ge \delta\). For each test-time node \(j\) (with unknown \(y_j\)), we consider the null hypothesis \[H_j:\;\kappa(j)\ge \delta \qquad \text{(the imputation at j is risky)}. \label{eq:hypothesis-risky}\tag{2}\]

3.3.0.2 Proxy risk score.

Since \(\kappa(j)\) is unobservable at test time, conformal selection requires a label-free proxy score that tends to be larger for risky imputations. We therefore define a proxy risk score by combining two complementary signals: prediction instability under structure-preserving perturbations, capturing sensitivity to small changes in relational evidence; and an evidence penalty, capturing the scarcity of relational support even when predictions appear stable.

To compute prediction instability, we perform margin-weighted edge perturbations that reflect how strongly each edge is supported by the edge construction rules. For each edge \(e\), we compute a normalized margin \(m(e)\in[0,1]\) to the threshold and convert it into an edge-specific keep probability \[\pi(e) \;=\; \pi_{\min} + (\pi_{\max}-\pi_{\min})\cdot m(e)^{\gamma}, \label{eq:keep95prob}\tag{3}\] where \(\pi_{\min},\pi_{\max}\in(0,1)\) and \(\gamma>0\). Edges that barely satisfy the threshold have small \(m(e)\) and are dropped more often, while strongly supported edges have large \(m(e)\) and are kept with higher probability.

For a temporal edge \(e=(i,j)\in\mathcal{E}_t\): \[m_t(e)\;=\;1-\min\!\left\{ \frac{\Delta t_{ij}}{\Delta_t},\,1\right\}. \label{eq:temporal95margin}\tag{4}\] and for a value edge \(e=(i,j)\in\mathcal{E}_v\),

\[s_v(e)\;=\;\max\!\left\{\frac{d^{v}_{ij}}{\tau_{v}},\,\frac{d^{\delta}_{ij}}{\tau_{\delta}}\right\}, \quad m_v(e)=1-\min\{s_v(e),1\} \label{eq:value95margin}\tag{5}\] If the trend is unavailable for at least one endpoint, we use value-only strength \(s_v(e)=d^v_{ij}/\tau_v\).

Let \(\mathcal{G}^{(1)},\ldots,\mathcal{G}^{(K)}\) be \(K\) independent perturbations obtained by sampling each edge \(e\) according to \(\mathrm{Bernoulli}(\pi(e))\). Let \(\hat{y}_i^{(k)}\) be the model prediction for node \(i\) when running inference on \(\mathcal{G}^{(k)}\). We define the prediction-instability score as the dispersion of these predictions: \[S_{\mathrm{pred}}(i) =\mathrm{Std}\Big(\{\hat{y}_i^{(k)}\}_{k=1}^{K}\Big), \label{eq:spred}\tag{6}\] where \(\mathrm{Std}(\cdot)\) denotes the standard deviation across perturbations.

Beyond perturbation sensitivity, reliability also depends on the amount of relational evidence supporting a node. When the constructed relations are sparse for a node (hence small degree), there is limited information propagation from other events, and the imputation is driven by a weak neighborhood signal. To reflect this evidence scarcity, we incorporate a degree-based penalty that downweights nodes with insufficient temporal/value support, independent of the instability term.

Concretely, let \(\deg_t(i)\) and \(\deg_v(i)\) denote the degrees of node \(i\) in the temporal and value graphs, respectively. We define the evidence penalty as \[S_{\mathrm{evid}}(i) = \frac{1}{\sqrt{\deg_t(i)+1}} + \frac{1}{\sqrt{\deg_v(i)+1}}, \label{eq:sevid}\tag{7}\] so that nodes with weaker neighborhood evidence receive larger penalties.

Finally, we define the proxy risk score used for conformal testing as \[\hat{\kappa}(i) = S_{\mathrm{pred}}(i)+\beta\,S_{\mathrm{evid}}(i), \label{eq:kappa95hat}\tag{8}\] where \(\beta\ge 0\) balances instability and evidence. Intuitively, \(\hat{\kappa}(i)\) becomes large when predictions are sensitive to plausible perturbations or when the node is weakly supported by the constructed relations.

3.3.0.3 Conformal p-values.

Given the proxy risk score \(\hat{\kappa}(i)\), we test the risk hypothesis in 2 by constructing conformal p-values using the calibration set. Concretely, we partition nodes into three disjoint subsets, \(\Omega_{\mathrm{train}}\), \(\Omega_{\mathrm{cal}}\), and \(\Omega_{\mathrm{test}}\). We first define the calibration “bad” subset \[\mathcal{B}=\{i\in\Omega_{\mathrm{cal}}:\;|y_i-\hat{y}(i)|\ge \delta\}, \label{eq:bad95set}\tag{9}\] which provides an empirical reference set for the null population. A common concern for applying conformal prediction on graphs is that node dependencies may violate exchangeability. Recent studies show that conformal validity can still hold for node-level tasks on static graphs when the nonconformity score is invariant to permutations of the calibration/test nodes [23]. In our setting, calibration and test nodes are treated identically during graph construction: their features and edges are determined solely by observed covariates, with targets masked. Under this setting, we give the following exchangeability assumption between \(\Omega_{\mathrm{cal}}\) and \(\Omega_{\mathrm{test}}\).

For each node \(j\), we compute a conformal p-value [21] against \(\mathcal{B}\) as \[p_j = \frac{ 1+\sum_{i\in\mathcal{B}}\mathbb{1}\{\hat{\kappa}(i)<\hat{\kappa}(j)\} +U_j\sum_{i\in\mathcal{B}}\mathbb{1}\{\hat{\kappa}(i)=\hat{\kappa}(j)\} }{ |\Omega_{\mathrm{cal}}|+1 }, \label{eq:conformal95p}\tag{10}\] where \(U_j\sim\mathrm{Unif}[0,1]\) is used for random tie-breaking. Intuitively, a smaller \(p_j\) means that the current imputation has unusually low proxy risk compared with the calibration bad set \(\mathcal{B}\), and is therefore less likely to be risky.

3.3.0.4 FDR-controlled selection via Benjamini–Hochberg.

In clinical deployment, it is more meaningful to control the overall error burden among the imputed values that are actually presented to clinicians or downstream decision models. We therefore target control of the false discovery rate (FDR) [22], defined as the expected proportion of truly risky imputations among those we accept for deployment. To this end, we apply the Benjamini–Hochberg (BH) procedure [21], [22] to the collection of p-values \(\{p_j\}_{j\in\Omega_{\mathrm{test}}}\) at a target level \(\alpha\in(0,1)\). Let \(p_{(1)}\le \cdots \le p_{(m)}\) be the sorted p-values over the \(m=|\Omega_{\mathrm{test}}|\) test nodes, and let \[k^\star=\max\Big\{k:\;p_{(k)}\le \frac{k}{m}\,\alpha\Big\}, \label{eq:bh95cutoff}\tag{11}\] with the convention that \(k^\star=0\) if the set is empty. We then select the nodes \[\mathcal{S}(\alpha)=\{j\in\Omega_{\mathrm{test}}:\;p_j\le p_{(k^\star)}\}. \label{eq:bh95select}\tag{12}\] Since the null hypothesis \(H_j\) corresponds to “the imputation at \(j\) is risky”, rejecting \(H_j\) means selecting node \(j\) as sufficiently reliable for release. Therefore, \(\mathcal{S}(\alpha)\) is the deployed subset with a user-specified FDR control level.

The proof follows the conformal testing used in [19], [21]; for completeness we provide it in Appendix 7.

4 Results↩︎

In this section, we analyze four key aspects to demonstrate the effectiveness of SafeImpute: (i) overall imputation accuracy and reliability control of the proposed method compared with various baselines; (ii) the performance of the pure event-graph imputer without FDR control; (iii) an ablation study that quantifies the impact of major design components; and (iv) a detailed analysis of FDR-controlled selection.

4.1 Experimental Settings↩︎

4.1.1 Datasets↩︎

We evaluate SafeImpute  on HbA1c (A1c) imputation for diabetes patients, where HbA1c is a key indicator for long-term glycemic control and diabetes management, using three datasets: Mayo Clinic, MIMIC-III [25], and MIMIC-IV [26]. Mayo Clinic dataset contains real-world records collected in 2022–2023 from patients diagnosed with type 1 diabetes, where each patient has multiple visits; for each visit, we extract {A1c, glucose, cholesterol, gender} as features, and the lab variables exhibit naturally occurring missingness because tests are not ordered at every visit. To create a test set for evaluation, we additionally mask 10% of observed A1c values as held-out targets. For MIMIC-III and MIMIC-IV, we construct diabetes cohorts by filtering patients with diabetes diagnoses and treat each A1c measurement as an anchor visit; for each anchor, we retrieve the most recent glucose and cholesterol measurements recorded at or before the A1c timestamp (within 7 days for glucose and 90 days for cholesterol), and use gender as the demographic feature, yielding irregular longitudinal event sequences with naturally occurring missingness across patients. We randomly mask 40% of observed A1c values as held-out targets. Apart from the held-out test nodes, we split the remaining nodes into training/validation/calibration sets with a 70%/15%/15% ratio. Statistics of processed datasets and detailed analysis of Mayo Clinic data are provided in the appendix 8. Definitions of the evaluation metrics are provided in Appendix 9.

4.1.2 Baselines↩︎

Existing imputation baselines generally do not provide native uncertainty estimation or selective-release mechanisms. Since prior work on uncertainty estimation for imputation [27] suggests that methods such as MICE [28] and MIWAE [29] can reflect both calibration and accuracy, we include a diverse set of baselines spanning statistical, deep learning, and sequential methods. Statistical methods: Mean imputes each feature by the training-set mean. KNN imputes from the k nearest samples based on similarity over observed dimensions. MICE [28] iteratively fits feature-wise conditional regressors to update missing entries. MissForest [30] performs iterative imputation using random-forest predictors. HyperImpute [31] selects a strong imputer via automated model and hyperparameter search. Deep learning methods: GAIN [32] uses adversarial training to generate realistic imputations. GRAPE [6] leverages graph neural message passing to propagate information for missing values. MIWAE [29] performs imputation under a variational generative framework. TDM [33] models complex tabular dependencies with a deep imputation architecture. ReMasker [7] learns masked reconstruction with a denoising objective. DiffImputer [8] uses diffusion-based denoising to generate completed samples. Sequential methods: LSTM [34] imputes using recurrent hidden states over irregular event sequences. TRANS [35] is a transformer-based sequential model for electronic health records data prediction; we adapt it for imputation by casting missing-value recovery as a prediction problem over masked targets. We evaluate all methods on the same set of predictions selected by the proxy risk score and report each method’s empirical Precision based on its own imputation outputs.

4.1.3 Implementation↩︎

We choose the target FDR level \(\alpha\) according to the empirical FDR–acceptance trade-off on each dataset. Specifically, we set \(\alpha{=}0.15\) for Mayo Clinic and \(\alpha{=}0.35\) for MIMIC-III/IV. Very strict target levels can yield an extremely small or empty accepted set, making the selective results unstable and less practically useful. We set the clinically interpretable threshold \(\delta\) to \(0.6\) percentage points for Mayo Clinic and \(1.0\) percentage point for MIMIC-III/IV. For all datasets, we use an auxiliary loss weight \(\lambda{=}0.1\) and set \(\beta{=}0.1\) when computing the proxy risk score. All experiments are run on a single NVIDIA H200 GPU. For baselines, we use the default PyPI implementation for MICE; for HyperImpute and MissForest, we use their official packages; and for KNN, we set \(k{=}10\). For other baselines, we adopt the original implementations and follow the hyperparameter recommendations in the corresponding papers, with minor tuning for each dataset, and report the best-performing configuration. Each experiment is repeated with five random seeds, and we report the mean performance.

4.2 Overall Performance↩︎

Table 1 shows the overall performance of SafeImpute against a broad set of baselines on Mayo Clinic, MIMIC-III, and MIMIC-IV. Overall, SafeImpute is consistently top-performing across datasets, demonstrating strong performance with FDR control and robustness under irregular and sparse missing data. Beyond the overall gains, comparing baseline behaviors to our method yields several findings. (i) Among non-neural baselines that are still commonly used in practice, MissForest and HyperImpute are relatively strong competitors, but they become less competitive in more challenging settings (e.g., when missingness is severe); notably, after applying our FDR-controlled selection, these simple methods also improve compared to their uncontrolled results (Table 2), suggesting that risk-aware selection can be a broadly useful add-on for simple imputers. (ii) LSTM becomes less competitive, indicating that when sufficient temporal context is not available, classic time-series models remain sensitive to irregular long gaps and sparse supervision. (iii) GRAPE, as a purely graph-based method, performs strongly, reinforcing the effectiveness of capturing intra-patient and inter-patient relations via the graph. (iv) TRANS achieves comparable performance on Mayo Clinic, implying that combining a temporal Transformer with graph modeling can be effective when sequences are short and irregular; at the same time, our comparison suggests that high-quality error control often requires method-aware design to support more reliable selection.

3.0pt

4.3 Event Graph Evaluation↩︎

We then evaluate the pure imputation performance on the event graph without any error control. In this setting, we construct edges only using the proposed temporal and trend-aware value criteria, and we do not add extra value-only edges when no eligible trend-based neighbor exists, to avoid noisy connections. Table 2 reports results on Mayo Clinic, MIMIC-III, and MIMIC-IV. The proposed event-graph imputer achieves the best or second best performance across all three datasets, highlighting the benefit of event-graph construction and a two-relation GNN that combines intra-patient temporal continuity with inter-patient value-based similarity under irregular and sparse visit patterns. Sequential baselines are less competitive, suggesting that irregular visits and sparse observations limit sequence-only modeling, while classical imputers remain behind because feature-conditional modeling cannot fully capture latent relational signals under severe missingness. Additionally, the pure event-graph imputer sometimes slightly outperforms the FDR-controlled version because the latter adds value-only fallback edges to support reliability estimation, which may introduce noisy cross-patient connections. This reflects a small accuracy and reliability trade-off for selective release.

4.4 Ablation Study↩︎

We conduct ablation studies on Mayo Clinic and MIMIC-III to assess both the event-graph imputer and the reliability-control module. Removing temporal edges or trend-aware value edges degrades performance, showing that intra-patient temporal continuity and inter-patient clinical similarity provide complementary signals. Removing the adaptive fusion gate causes a substantial drop, highlighting the importance of dynamically balancing the two relations. The auxiliary reconstruction loss improves performance on Mayo Clinic and serves as an additional regularizer for learning event representations from sparse non-target labs.

We also ablate the proxy risk score for FDR-controlled selection. Using prediction instability alone leads to much higher FDR than the full proxy score, increasing FDR from 0.0000 to 0.3575 on Mayo Clinic and from 0.3667 to 0.6304 on MIMIC-III. This confirms that the evidence penalty is important for reliable release under sparse relational support.

4.5 A Close Look at FDR Control↩︎

4.5.0.1 FDR and Power.

Recall that our test-time null hypothesis is \(H_j:\;|y_j-\hat{y}(j)|\ge \delta\), i.e., node \(j\) has a risky imputation. Our goal is to reject these risky-imputation nulls and retain nodes whose imputations are sufficiently reliable for release. Figure 2 sweeps the BH target level \(\alpha\) and reports the FDR among the selected set, defined as the fraction of selected nodes that are in fact risk points (with \(|y_j-\hat{y}(j)|\ge\delta\)). We also report power, measured as the fraction of non-risk nodes (with \(|y_j-\hat{y}(j)|<\delta\)) that are successfully retained. Across \(\delta\in\{0.5,0.6,0.7,0.8\}\), the achieved FDR closely follows the target, indicating conservative control. As \(\alpha\) increases, more hypotheses are rejected and more nodes are selected, leading to higher power and eventual saturation when the BH cutoff becomes non-restrictive. Larger tolerances \(\delta\) make the risky-imputation null less likely to hold, resulting in fewer true risk points; consequently, the retained set attains higher power and lower FDR, consistent with the intended screening behavior. FDR curves on MIMIC-III and MIMIC-IV are provided in Appendix 10.

4.5.0.2 Selective Release Summary

Table 5 summarizes the selective release results under the chosen operating point for each dataset. The retained subset is selected by the pre-specified proxy risk score, conformal calibration, and the Benjamini–Hochberg procedure, without access to test labels during selection. The results show that SafeImpute retains a nontrivial fraction of imputations under the target FDR constraint, rather than achieving error control by selecting only a very small subset. The acceptance fraction reflects the coverage of released imputations under the target FDR constraint, while power measures how many truly reliable imputations are successfully retained. Across datasets, the selected operating points illustrate the trade-off between empirical FDR and retained coverage: stricter selection yields lower FDR but smaller acceptance fractions, while more permissive target levels retain more imputations at the cost of higher empirical FDR.

4.5.0.3 Quality of the Proxy Risk Score.

Our FDR-controlled selection relies on conformal p-values constructed from the proxy risk score \(\hat{\kappa}\), which in turn determines the BH rejection set. A well-behaved proxy should ideally exhibit a monotonic relationship with the realized imputation error: nodes with larger \(\hat{\kappa}\) should, on average, have larger realized errors \(|y-\hat{y}|\). Such alignment is crucial because it makes the ranking induced by \(\hat{\kappa}\) informative for separating high-error nodes from low-error ones. Here, we evaluate this property on both the Mayo Clinic and MIMIC-III datasets. For each dataset, we partition evaluation nodes into quantile groups (5 for Mayo Clinic and 10 for MIMIC-III) according to \(\hat{\kappa}\) and plot the mean absolute error within each group against the mean proxy score. Figure 3 shows a clear increasing trend on both datasets: groups with higher average \(\hat{\kappa}\) consistently have higher mean absolute residuals. This indicates that \(\hat{\kappa}\) is a meaningful surrogate of the unobserved error and is suitable for conformal selection when test labels are unavailable.

5 Related Work↩︎

5.1 Data Imputation↩︎

Missing data is a pervasive challenge in many applications [36]. Classical methods, e.g., simple heuristics (mean/median substitution) and machine learning methods such as KNN imputation [4] offer straightforward solutions but rely on restrictive assumptions and ignore feature relationships. Deep learning approaches leverage various techniques and training paradigms [37]–[40] (e.g., GANs, GNNs, and diffusion models). Representative examples include GAIN [32], which trains a generator to impute missing entries conditioned on observed entries while a discriminator distinguishes observed entries from imputed ones; GRAPE [6], which represents tabular data [41] as a bipartite graph with observed values treated as edges and formulates imputation as an edge-level prediction task; and DiffPuter [8], which combines diffusion models with an EM-style training procedure to learn the joint distribution of complete data and perform conditional sampling for missing entries. Autoencoder-based methods [42], [43], such as ReMasker [7], adapt masked autoencoding to tabular imputation by masking entries and training a deep learning-based model to reconstruct entries from context. Data imputation is particularly important in healthcare settings, where missing data is common due to irregular measurements and clinical workflows. Many of these ideas have also been adapted to clinical records. For example, GRU-D [44] and BRITS [45] leverage RNNs to model temporal dynamics. Attention-based models such as SAITS [46] use self-attention with self-supervised masking to better capture long-range dependencies of clinical variables. GRIN [47] exploits GNNs [48], [49] and performs message passing across variables and time. More recently, diffusion-based approaches such as CSDI [50] perform conditional denoising to iteratively refine missing values and have shown strong performance on ICU benchmarks. Although these methods are effective, they primarily optimize the reconstruction error and overlook the reliability of imputations.

5.2 Uncertainty Quantification in Healthcare Models↩︎

Given the high stakes of the domain, modern healthcare machine learning models increasingly incorporate uncertainty quantification (UQ) techniques to improve trustworthiness and safety [51], [52]. Bayesian approaches [53], [54] place probabilistic priors on model parameters to capture epistemic uncertainty and have been used in medical applications [55], [56] to indicate low-confidence predictions. For example, Qiu et al. [57] applied a Bayesian neural network to clinical records and found that instances with high predictive uncertainty were harmful to overall model performance. Related approximations include Monte Carlo dropout [58], which performs stochastic inference by enabling dropout at test time; TrUE-Net [59] applies this strategy to AD genomic variant classification. Deep ensembles are another empirical way to improve predictive performance and quantify the uncertainty via inter-model variability. DEGU [60] integrates ensemble learning to improve the robustness and explainability of genomic DNNs. Another complementary direction is conformal prediction (CP), which offers distribution-free uncertainty quantification with formal coverage guarantees through calibrated prediction sets. Building on CP, conformal selection [21] accepts predictions only when the conformal prediction set is sufficiently informative, yielding an explicit and risk-controllable mechanism for accepting a prediction [61]. This paradigm has been increasingly adopted in recent clinical studies [18], [21], [62]. SAFER [19] uses conformal inference to provide statistical guarantees while filtering out uncertain treatment recommendations. Together, these advances motivate the integration of uncertainty quantification to support more reliable and risk-aware clinical AI.

6 Conclusion↩︎

We propose SafeImpute, a framework for reliable clinical data imputation that models irregular and sparse longitudinal records by constructing an event graph capturing both intra-patient temporal structure and inter-patient clinical similarity. For reliability control, SafeImpute calibrates a proxy risk score into conformal p-values and applies the Benjamini–Hochberg procedure to control the false discovery rate of clinically unacceptable imputation errors. Extensive experiments demonstrate that SafeImpute achieves strong imputation performance while providing meaningful reliability guarantees, outperforming selected baselines in both standard and FDR-controlled imputation evaluation. Future work includes: (i) incorporating multi-modal information and multi-source clinical data; and (ii) moving beyond imputation metrics to evaluate the effect of reliable imputations on downstream tasks, such as risk prediction and treatment recommendation.

Ethical Use of Data↩︎

To protect patient privacy, all data were anonymized during collection and processing. The use of this dataset has been reviewed and approved by the Institutional Review Board (IRB) to ensure compliance with ethical standards and participant protections. We strictly adhere to relevant laws and regulations regarding data use and protection to ensure the legal and compliant use of the data.

Acknowledgments↩︎

This work is supported by National Science Foundation under Award No. IIS-2117902. The views and conclusions are those of the authors and should not be interpreted as representing the official policies of the funding agencies or the government.

7 Proof of FDR Control↩︎

7.0.0.1 Goal.

Let \(\kappa(i)=|y_i-\hat{y}(i)|\) be the true imputation error on labeled nodes, and fix a clinical tolerance \(\delta>0\). For each test node \(j\in\Omega_{\mathrm{test}}\) we test the null \[H_j:\;\kappa(j)\ge \delta,\] i.e., the imputation at \(j\) is risky.

Let \(\hat{\kappa}(x)\) be our label-free proxy risk score in 8 . To apply the selective conformal testing argument in [21], we define a nonconformity score that is monotone in the binary label \(\mathbb{1}\{\kappa(i)\ge \delta\}\). Let \(M\ge \sup_i \hat{\kappa}(i)\) (finite since \(\hat{\kappa}\) is bounded in practice). Define, for any node \(i\) and candidate risk label \(b\in\{0,1\}\), \[J(x,b)\;=\;\hat{\kappa}(x)\;+\;2M\cdot \mathbb{1}\{b=1\}. \label{eq:J95def}\tag{13}\] Thus, if the node is not risky (\(b=0\)), then \(J(x,b)=\hat{\kappa}(x)\); if it is risky (\(b=1\)), then \(J(x,1)=\hat{\kappa}(x)+2M\). The score \(J(i,b)\) is monotone in terms of the risk label \(b\).

For each calibration node \(i\in\Omega_{\mathrm{cal}}\), define its observed risk label \[B_i=\mathbb{1}\{\kappa(i)\ge \delta\},\] and for each test node \(j\in\Omega_{\mathrm{test}}\), under the null \(H_j\) we have \(B_j=1\). Let \[J_i = J(i,B_i)\quad (i\in\Omega_{\mathrm{cal}}), \qquad \widehat{J}_j = J(j,1)\quad (j\in\Omega_{\mathrm{test}}),\] i.e., we evaluate the test score at the null label \(b=1\).

7.0.0.2 Conformal p-values.

Let \(n=|\Omega_{\mathrm{cal}}|\) and index calibration nodes as \(\{1,\ldots,n\}\), and test nodes as \(\{n+1,\ldots,n+m\}\) with \(m=|\Omega_{\mathrm{test}}|\). Using the standard conformal rank with random tie-breaking, define for each test node \(j\in\{1,\ldots,m\}\), \[p_j = \frac{ \sum_{i=1}^{n}\mathbb{1}\{J_i<\widehat{J}_{n+j}\} \;+\; U_j\Big(1+\sum_{i=1}^{n}\mathbb{1}\{J_i=\widehat{J}_{n+j}\}\Big) }{ n+1 } \label{eq:p95def}\tag{14}\] where \(U_j\sim\mathrm{Unif}[0,1]\). This is exactly the conformal p-value construction, as defined in [21], for testing the null label \(b=1\), with a nonconformity score monotone in \(b\).

7.0.0.3 Super-uniformity.

By exchangeability of calibration and test nodes (Assumption 1) and monotonicity of \(J\) in 13 , the p-values in 14 are selectively super-uniform with respect to the nulls \(H_j\): for every \(\alpha\in[0,1]\), \[\mathbb{P}\big(\,(j\in S)\;\wedge\;(p_j\le \alpha)\,\big)\;\le\;\alpha, \label{eq:selective95superuniform}\tag{15}\] where \(S\) denotes the set of tested hypotheses produced by the procedure [21].

7.0.0.4 BH step and FDR control.

Apply Benjamini–Hochberg (BH) [22] at level \(\alpha\) to \(\{p_j\}_{j=1}^m\). Let \(p_{(1)}\le \cdots \le p_{(m)}\) and \[k^\star=\max\Big\{k:\;p_{(k)}\le \frac{\alpha k}{m}\Big\},\qquad \mathcal{R}=\{j:\;p_j\le p_{(k^\star)}\},\] where \(\mathcal{R}\) is the BH rejection set. In our setting, rejecting \(H_j\) means declaring node \(j\) reliable (i.e., deploying its imputation). Let \(R=|\mathcal{R}|\) and define the number of false discoveries (deployed but truly risky) as \[V=\sum_{j=1}^m \mathbb{1}\{j\in\mathcal{R},\;H_j\;\text{true}\}.\] The FDR is \(\mathrm{FDR}=\mathbb{E}\big[\frac{V}{\max\{R,1\}}\big]\).

To obtain FDR control after selecting a deployed subset via Benjamini–Hochberg, we further assume the null p-values are independent or satisfy PRDS [24] [20], [24]. Define indicator variables \(R_j=\mathbb{1}\{j\in\mathcal{R}\}\) and \(T_j=\mathbb{1}\{H_j\;\text{true}\}\). Then \[\mathrm{FDR} = \mathbb{E}\!\left[\frac{\sum_{j=1}^m T_j R_j}{\max\{1,\sum_{j=1}^m R_j\}}\right] = \mathbb{E}\!\left[\frac{1}{\max\{1,\sum_{j=1}^m R_j\}}\sum_{j=1}^m T_j R_j\right].\] From the classical BH guarantee under PRDS super-uniform p-values [19], [21], [22], [24], we have for each \(j\), \[\mathbb{E}[T_j R_j]\;\le\;\alpha\,\mathbb{E}[R_j].\] Summing over \(j\) and applying the standard BH bounding argument yields \[\mathrm{FDR} = \mathbb{E}\!\left[\frac{\sum_{j=1}^m T_j R_j}{\max\{1,\sum_{j=1}^m R_j\}}\right] \le \alpha.\] Therefore, among the deployed imputations (BH rejections), the expected fraction of truly risky ones is controlled at level \(\alpha\). This completes the proof.

8 Datasets Analysis↩︎

8.0.0.1 Dataset statistics.

Table 6 reports the basic statistics of the three preprocessed datasets used in our study. The Mayo Clinic dataset was collected from patients diagnosed with type 1 diabetes who had multiple clinical visits at Mayo Clinic during 2022–2023. We focus on routinely used metabolic and anthropometric measurements, including HbA1c (A1C), glucose, cholesterol and lipid panel components, and basic body measurements (weight, height, BMI), together with basic demographics. MIMIC-III and MIMIC-IV are publicly available critical-care EHR datasets from the Beth Israel Deaconess Medical Center. MIMIC-III covers ICU stays between 2001 and 2012, while MIMIC-IV contains more recent ICU and emergency department records from 2008 to 2019.

8.0.0.2 Missingness patterns across measurements.

Figure 4 visualizes missingness across visits and lab measurements in the Mayo Clinic data. Rows correspond to visits sorted by patient and time, and columns correspond to measurements. Missingness is highly heterogeneous: some measurements are recorded consistently, whereas lipid panel components such as LDL, HDL, and triglycerides show substantially higher missing rates. This pattern reflects opportunistic test ordering in routine care.

8.0.0.3 Heterogeneous longitudinal follow-up.

Figure 5 (left) shows the distribution of the number of visits per patient in Mayo Clinic data. Most patients have a moderate number of visits, while a small fraction have substantially fewer or more visits. This heterogeneity means that some patients provide limited intra-patient history, motivating the use of cross-patient information when individual trajectories are short.

8.0.0.4 Irregular inter-visit intervals.

Figure 5 (right) reports inter-visit time gaps, in days, between consecutive visits within each patient in Mayo Clinic data. The distribution is highly non-uniform and is consistent with periodic diabetes monitoring, while also showing substantial irregularity across patients.

8.0.0.5 Visit-wise missingness dynamics.

Figure 6 summarizes visit-wise missing ratios for 10 randomly sampled patients in Mayo Clinic data. Missingness varies both across visits within the same patient and across patients at the same visit index, further highlighting the irregular and opportunistic nature of measurement collection in real-world clinical practice.

9 Evaluation Metrics↩︎

Let \(\hat{y}_i\) denote the imputed value and \(y_i\) denote the ground-truth value for an evaluated node \(i\). The mean absolute error (MAE) is defined as \[\mathrm{MAE} = \frac{1}{|\Omega_{\mathrm{test}}|} \sum_{i \in \Omega_{\mathrm{test}}} |\hat{y}_i - y_i| .\] The root mean squared error (RMSE) is \[\mathrm{RMSE} = \sqrt{ \frac{1}{|\Omega_{\mathrm{test}}|} \sum_{i \in \Omega_{\mathrm{test}}} (\hat{y}_i - y_i)^2 }.\]

For selective-release evaluation, let \(\mathcal{S}\) denote the set of released imputations selected by the FDR-control procedure. Given the clinically meaningful error threshold \(\delta\), an imputation is considered reliable if \(|\hat{y}_i-y_i| < \delta\) and risky otherwise. We define precision as the fraction of released imputations that are truly reliable: \[\mathrm{Precision} = \frac{ |\{i \in \mathcal{S}: |\hat{y}_i-y_i| < \delta\}| }{ \max\{|\mathcal{S}|,1\} }.\] Equivalently, since empirical FDR measures the fraction of risky imputations among the released set, precision corresponds to \(1-\mathrm{FDR}\). Higher precision indicates that a larger fraction of released imputations are within the clinically acceptable error range.

10 FDR Curves on MIMIC-III and MIMIC-IV↩︎

Figure 7 provides additional FDR–power trade-off curves on MIMIC-III and MIMIC-IV. We fix the clinically interpretable tolerance to \(\delta=1.0\) and sweep the target FDR level \(\alpha\). These results show how the BH target level controls the selective release behavior on public clinical datasets: increasing \(\alpha\) makes the selection rule less restrictive and generally increases power, while the achieved FDR remains interpretable relative to the target level.

References↩︎

---

1. Accepted at the 32nd ACM SIGKDD Conference on Knowledge Discovery and Data Mining (KDD 2026). This is the authors’ accepted manuscript. The final version of record will appear in the ACM Digital Library with DOI: 10.1145/3770855.3817967.↩︎