Paper deep dive
Patients-like-me: A Variational LM--GNN Framework for Explainable Clinical Prediction
Xinyu Wang, Yixuan Li, Hanwei Wu, Qincheng Lu, Chi-Kuang Yeh, Xiao-Wen Chang, Ziyang Song
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 93%
Last extracted: 8/9/2026, 3:08:51 AM
Summary
The paper introduces Patients-like-me (PLM), a unified framework combining Language Models (LMs) and Graph Neural Networks (GNNs) for clinical prediction. PLM uses a Variational Expectation-Maximization (VEM) algorithm to iteratively refine patient representations by integrating local textual semantics from LMs with global cohort structure from GNNs. It demonstrates superior performance on MIMIC-III and MIMIC-IV datasets for tasks like readmission prediction and drug recommendation, while providing explainable reference-patient attributions.
Entities (9)
Relation Signals (8)
Patients-like-me → evaluatedon → MIMIC-IV
confidence 99% · Extensive experiments on MIMIC-III and MIMIC-IV show that PLM consistently outperforms...
Patients-like-me → evaluatedon → MIMIC-III
confidence 99% · Extensive experiments on MIMIC-III and MIMIC-IV show that PLM consistently outperforms...
Patients-like-me → uses → Language Models
confidence 95% · PLM couples a patient-level LM with a cohort-level GNN within a VEM framework.
Patients-like-me → uses → Graph Neural Networks
confidence 95% · PLM couples a patient-level LM with a cohort-level GNN within a VEM framework.
Patients-like-me → trainedwith → Variational Expectation-Maximization
confidence 92% · we introduce a Variational Expectation-Maximization algorithm that alternates LM and GNN updates...
Patients-like-me → performstask → Readmission Prediction
confidence 90% · Experiments on MIMIC-III and MIMIC-IV across readmission prediction... show that PLM consistently outperforms...
Patients-like-me → performstask → Drug Recommendation
confidence 90% · Experiments on MIMIC-III and MIMIC-IV across readmission prediction... and drug recommendation show that PLM consistently outperforms...
Patients-like-me → provides →
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Language models (LMs) offer strong textual representations for electronic health records (EHRs), but they encode patient sequences in isolation and provide limited explainability. Graph neural networks (GNNs) complement LMs by incorporating inter-patient relationships and enabling reference-patient attribution, yet they rely on high-quality patient representations. We propose Patients-like-me (PLM), a unified LM--GNN framework that integrates local patient semantics with global cohort structure. To train PLM efficiently, we introduce a Variational Expectation-Maximization algorithm that alternates LM and GNN updates under a supervised variational objective. Extensive experiments on MIMIC-III and MIMIC-IV show that PLM consistently outperforms state-of-the-art methods, with improvements generalizing across encoder-only and decoder-only LM backbones. These gains are achieved with only modest additional computational overhead. PLM also provides reference-patient explanations by retrieving influential similar patients, while edge-masking experiments confirm that the highest-ranked references have the greatest impact on model predictions.
Tags
Links
- Source: https://arxiv.org/abs/2608.04193v1
- Canonical: https://arxiv.org/abs/2608.04193v1
PDF not stored locally. Use the link above to view on the source site.
Full Text
103,074 characters extracted from source content.
Expand or collapse full text
Patients-like-me: A Variational LM–GNN Framework for Explainable Clinical Prediction Xinyu Wang1 Yixuan Li111footnotemark: 1 Hanwei Wu2 Qincheng Lu1 Chi-Kuang Yeh3 Xiao-Wen Chang1 Ziyang Song4 1McGill University 2Université de Montréal 3Georgia State University 4Ohio University Equal contribution. Emails: xinyu.wang5@mail.mcgill.ca and yixuan.li2@mail.mcgill.ca.Corresponding author: ziyangs@ohio.edu. Abstract Language models (LMs) offer strong textual representations for electronic health records (EHRs), but they encode patient sequences in isolation and provide limited explainability. Graph neural networks (GNNs) complement LMs by incorporating inter-patient relationships and enabling reference-patient attribution, yet they rely on high-quality patient representations. We propose Patients-like-me (PLM), a unified LM–GNN framework that integrates local patient semantics with global cohort structure. To train PLM efficiently, we introduce a Variational Expectation-Maximization algorithm that alternates LM and GNN updates under a supervised variational objective. Extensive experiments on MIMIC-I and MIMIC-IV show that PLM consistently outperforms state-of-the-art methods, with improvements generalizing across encoder-only and decoder-only LM backbones. These gains are achieved with only modest additional computational overhead. PLM also provides reference-patient explanations by retrieving influential similar patients, while edge-masking experiments confirm that the highest-ranked references have the greatest impact on model predictions. Patients-like-me: A Variational LM–GNN Framework for Explainable Clinical Prediction Xinyu Wang1†thanks: Equal contribution. Emails: xinyu.wang5@mail.mcgill.ca and yixuan.li2@mail.mcgill.ca. Yixuan Li111footnotemark: 1 Hanwei Wu2 Qincheng Lu1 Chi-Kuang Yeh3 Xiao-Wen Chang1 Ziyang Song4†thanks: Corresponding author: ziyangs@ohio.edu. 1McGill University 2Université de Montréal 3Georgia State University 4Ohio University 1 Introduction Electronic health records (EHRs) contain heterogeneous patient records spanning one or more visits, including structured codes (e.g., diagnoses, procedures, medications) and unstructured free-text notes (Xu et al., 2022). The accumulation of large-scale EHR data enables representation learning that supports accurate supervised clinical prediction for clinical decision support, including readmission and length-of-stay (LOS) prediction as well as personalized medication recommendation (Wang et al., 2024; Shang et al., 2019b; Yang et al., 2021b). Language models (LMs) are a natural backbone for EHR modeling because they can encode structured codes and free-text notes into textual representations (Wornow et al., 2025). Recent work therefore leverages LMs as patient-level encoders, mapping each patient’s EHR sequence to a representation for downstream predictions (Makarov et al., 2025; Hegselmann et al., 2025). However, LM-based encoders are inherently local, representing each patient sequence in isolation, thus overlooking global cohort structure such as shared phenotypes and comorbidities. This hurt tasks that benefit from cohort-level signals, such as personalized drug recommendation, where clinically similar reference patients (i.e., “patients like me”) provide evidence for prediction (Alsentzer et al., 2025; Kauffman et al., 2025). For example, among patients with congestive heart failure, medication selection is often influenced by shared comorbidities such as kidney dysfunction or diabetes. In addition, LMs offer limited explainability, making it difficult to provide explanations linked to reference cases. Graph neural networks (GNNs) complement LMs by modeling inter-patient relations on a text-attributed graph (TAG), where nodes represent patient sequences and edges reflect cohort similarity (Yan et al., 2023; Boll et al., 2024). These structural relationships preserve patient proximity, such that connected patients are more likely to exhibit similar clinical profiles. GNN message passing propagates cohort-level signals to refine patient representations, addressing the locality limitation of LMs. In addition, graph-based attribution methods can link predictions to influential neighboring cases, yielding reference-patient explanations for prediction (Ying et al., 2019). This motivates a unified LM–GNN framework in which LMs encode patient representations and GNNs refine them using global cohort structure. Existing LM-GNN methods mainly adopt two-stage, end-to-end, or alternating optimization. Two-stage methods train the LM and GNN separately, offering scalability but preventing graph information from iteratively refining the LM (Li et al., 2024). End-to-end methods jointly optimize both components, enabling tighter interaction but incurring higher memory and computational costs (Huang et al., 2024; Hu et al., 2025). Alternating methods, such as GLEM (Zhao et al., 2023), update the LM and GNN in turn to balance interaction and efficiency, but are primarily developed for semi-supervised node classification with missing labels. Their formulation therefore does not address fully supervised patient-level clinical prediction. In this work, we propose Patients-like-me (PLM), a unified LM-GNN framework for clinical prediction. PLM formulates a Variational Expectation-Maximization (VEM) algorithm that alternates LM-GNN optimization under a supervised variational objective, allowing local patient semantics and global cohort structure to iteratively refine one another with modest computational overhead. Experiments on MIMIC-I and MIMIC-IV across readmission prediction, length-of-stay prediction, and drug recommendation show that PLM consistently outperforms existing clinical prediction models, LM-based baselines, and alternative LM-GNN training strategies. PLM is backbone-agnostic, improving both encoder-only and decoder-only models with only modest computational overhead. Beyond predictive performance, PLM provides reference-patient explanations that identify influential similar patients, with edge-masking experiments showing that the highest-ranked references have the strongest effect on predictions. 2 Related Work 2.1 Language Modeling for EHRs Modeling EHR sequences is central to clinical prediction. EHRs consist of heterogeneous patient records, where each visit may contain structured clinical codes and free-text notes (Xu et al., 2022). Early probabilistic models are employed to uncover latent structure from EHRs (Song et al., 2022). Subsequent work uses deep neural networks to learn representations from patient records, including multi-layer perceptrons (MLPs) (Miotto et al., 2016), recurrent networks (Choi et al., 2016a; Pham et al., 2016), convolutional networks (Ma et al., 2020), and attention-based models (Choi et al., 2016b, 2017). Recently, Transformer-based LMs have become the dominant paradigm for encoding EHR data into contextual patient representations (Lee et al., 2020; Huang et al., 2019; Rasmy et al., 2021; Labrak et al., 2024). These models span both encoder-only architectures (e.g., BioBERT (Lee et al., 2020), BioClinical ModernBERT (Sounack et al., 2025)) and decoder-only LMs (e.g., BioMistral, Meerkat) (Labrak et al., 2024; Kim et al., 2025). However, these biomedical LMs primarily encode local textual representations (Defilippo et al., 2024), motivating the incorporation of inter-patient structure to complement these representations (Alsentzer et al., 2025; Kauffman et al., 2025). Figure 1: Overview of Patients-like-me. (a) The LM encodes each patient sequence into an embedding, which the GNN refines using patient graph A for prediction and reference patient explanations. (b) The patient graph A retains edges with Aij≥τA_ij≥τ. (c) The variational EM algorithm alternates between an E-step that updates the LM while freezing the GNN and an M-step that updates the GNN while freezing the LM. 2.2 Graph-Based Patient Modeling and LMs GNNs on patient graphs are widely used to exploit cohort structure for clinical prediction (Boll et al., 2024; Vaida and Huang, 2025). Applying GNNs to EHRs requires mapping heterogeneous patient records to continuous node features (Gilmer et al., 2017), motivating TAGs and textual representation learning for node initialization (Yan et al., 2023). Prior work computes node embeddings from hand-crafted features or shallow text features (e.g., bag-of-words) (Kipf and Welling, 2017a). Sequential models have been adopted to map patient sequences into node embeddings for GNNs (Liu et al., 2021; Rocheteau et al., 2021; Dong et al., 2023). Recent work uses biomedical LMs to produce semantically rich embeddings for GNNs (Chen et al., 2024; Jiang et al., 2024). Most LM–GNN pipelines rely on two-stage or end-to-end optimization methods, but the former is often ineffective and the latter is not scalable to large data (Golmaei and Luo, 2021; Hu et al., 2025). Alternating methods iteratively update LMs and GNNs, but are tailored to semi-supervised node classification with partial labels (Zhao et al., 2023). Beyond training, prior GNN explainability methods such as GNNExplainer (Ying et al., 2019), PGExplainer (Luo et al., 2020), and SubgraphX (Yuan et al., 2021) provide edge- or subgraph-level attributions. In PLM, we leverage graph-based neighbor attribution to provide reference-patient explanations for clinical prediction. 3 Methodology 3.1 Patients-like-me Architecture For clinical prediction, each patient n is associated with an EHR sequence xnx_n and a target label yny_n. PLM couples a patient-level LM with a cohort-level GNN within a VEM framework (Fig. 1a). We first serialize each EHR into a text prompt by concatenating demographics with chronologically ordered visits, where each visit includes co-occurring diagnosis, procedure, and medication codes. The LM encodes the prompt into a patient embedding ono_n, which initializes the GNN node feature as hn(0)=onh_n^(0)=o_n. The GNN then propagates information over the patient graph to refine node representations hn(L)h_n^(L) for downstream prediction. Task-specific prompt templates are in Appendix B.1. As illustrated in Fig. 1b, we construct a weighted patient graph from the cohort. We represent it as a TAG =(,A,x)G=(N,A,x), where N is the set of patient nodes, A∈ℝN×NA ^N× N is the adjacency matrix, and each node is associated with a patient EHR sequence. We construct the patient-patient graph using diagnosis-code overlap, following prior work (Lu and Uddin, 2021). To prevent temporal leakage, edges for each prediction instance are computed only from diagnosis codes available within the corresponding input prefix. Specifically, the edge weight AijA_ij is defined as the number of diagnosis codes shared by patients i and j. This graph captures coarse inter-patient comorbidity structure, while patient-specific details are encoded by the LM. To maintain sparsity and suppress noisy links, we retain an edge only if Aij≥τA_ij≥τ. We select τ on the validation set and use it across all tasks, as it provides the best trade-off between predictive performance and graph sparsity (Section 5.4). We model the resulting weighted patient graph with a Graph Convolutional Network (GCN) (Kipf and Welling, 2017b), which supports weighted message passing under standard GCN normalization. After the patient-level split, we construct separate training, validation, and test graphs within each split. Details on split-specific graph construction, test-time message passing, and edge-weight analysis are provided in Appendix C.6. 3.2 VEM Algorithm As illustrated in Fig. 1c, we present a pseudo-likelihood VEM algorithm for PLM. We consider a probabilistic model in which patient sequences x are observations and target labels y are latent variables. Our goal is to infer the posterior distribution p(y∣x)p(y x), which is generally intractable (Fox and Roberts, 2012). Variational inference approximates this posterior with a tractable family qθ(y∣x)q_θ(y x) and optimizes the evidence lower bound (ELBO) (Blei et al., 2017): ELBO=qθ(y∣x)[logpϕ(x,y)−logqθ(y∣x)],ELBO=E_q_θ(y x) [ p_φ(x,y)- q_θ(y x) ], where θ and ϕφ denote the variational and model parameters, respectively. Following amortized variational inference (Kingma and Welling, 2022), we approximate the intractable posterior with a variational distribution qθ(y∣x)q_θ(y x) parameterized by an LM, which maps each patient sequence xnx_n to a distribution over yny_n by a single forward pass. We use a mean-field variational family that factorizes over the latent variables, qθ(y∣x)=∏n=1Nqθ(yn∣xn),q_θ(y x)= _n=1^Nq_θ(y_n x_n), where the factorization is imposed only on the variational family qθq_θ for scalable amortized inference. We instantiate the distribution pϕ(y∣x)p_φ(y x) with a GNN and omit the explicit dependence on the patient graph A. We optimize the ELBO using a VEM algorithm that alternates between E-step and M-step. In the E-step, we fix the GNN and update the LM by minimizing the KL divergence KL(qθ(y∣x)∥pϕ(y∣x))=qθ(y∣x)[logqθ(y∣x)−logpϕ(y∣x)].KL (q_θ(y x)\,\|\,p_φ(y x) )=E_q_θ(y x)\! [ q_θ(y x)- p_φ(y x) ]. This injects graph structural signals into the LM representations and yields a tighter ELBO. In the M-step, we fix the LM and update the GNN by maximizing a pseudo-likelihood objective (Besag, 1975), qθ(y∣x)[logpϕ(y∣x)]≈qθ(y∣x)[∑n=1Nlogpϕ(yn∣x,y∖n)],E_q_θ(y x) [ p_φ(y x) ] _q_θ(y x) [ _n=1^N p_φ (y_n x,y_N n ) ], where y∖ny_N n denotes the labels of all nodes except node n. It improves the GNN’s modeling of global structural signals using the LM representations. 3.3 E-Step: LM Optimization In the E-step, we fix the GNN and update the LM by maximizing the ELBO, thereby injecting global relational structure into the LM. Under the mean-field factorization, the optimal update for the factor of qθ(yn∣xn)q_θ(y_n x_n) for patient n is given by the following theorem. Theorem 3.1 (Optimal mean-field update of qθq_θ). Given the fixed pϕ(yn∣xn)p_φ(y_n x_n), the local optimum of the mean-field variational distribution qθ(yn∣xn)q_θ(y_n x_n), denoted by qθ∗(yn∣xn)q_θ^*(y_n x_n), satisfies logqθ∗(yn∣xn) q_θ^*(y_n x_n) =qθ(y∣x)[logpϕ(yn∣xn,yNB(n))] =E_q_θ(y x) [ p_φ(y_n x_n,y_NB(n)) ] +const, +const, (1) where the constant term is independent of yny_n and can be omitted during optimization. The proof is provided in Appendix A.2. Based on Theorem 3.1, computing the local optimum qθ∗(yn∣xn)q_θ^*(y_n x_n) requires evaluating an expectation under the current variational distribution qθ(y∣x)q_θ(y x). We estimate this expectation with a single-sample Monte Carlo estimator (Kingma and Welling, 2022). As shown in Appendix C.8, a single sample yields performance comparable to that of multiple samples, so we use a single sample in all experiments for efficiency. This Monte Carlo approximation is implemented using a single forward pass of the LM to obtain the patient embedding ono_n for each node, which is then used as the input node feature for the GNN, i.e., hn(0)=onh_n^(0)=o_n, yielding qθ(y∣x)[logpϕ(yn∣xn,yNB(n))]≃logpϕ(yn∣on,yNB(n)). split&E_q_θ(y x) [ p_φ(y_n x_n,y_NB(n)) ]\\ & p_φ(y_n o_n,y_NB(n)). split (2) Given equation 3.1 and equation 2, the local optimum qθ∗(yn∣xn)q_θ^*(y_n x_n) can thus be approximated by pϕ(yn∣on,yNB(n))p_φ(y_n o_n,y_NB(n)). This optimality condition indicates that, during the E-step, the LM is optimized through the fixed GNN, allowing its patient representations to incorporate both local semantics and neighborhood interactions. The E-step seeks to maximize the ELBO, equivalently minimizing KL[qθ(y∣x)||pϕ(y∣x)]KL [q_θ(y x)||p_φ(y x) ]. However, directly optimizing this KL is typically intractable as it requires computing over all possible qθ(y∣x)q_θ(y x). To address this issue, we adopt the wake–sleep algorithm, which updates the inference model qθq_θ by minimizing the reverse KL divergence (Hinton et al., 1995), θ∗=argminθKL(pϕ(y∣x)∥qθ(y∣x))=argmaxθpϕ(x,y)[logqθ(y∣x)]. splitθ^*&= _θKL (p_φ(y x)\,\|\,q_θ(y x) )\\ &= _θE_p_φ(x,y) [ q_θ(y x) ]. split (3) We omit the constant term independent of θ. Combining the above results, we obtain the supervised E-step objective for updating qθq_θ by fine-tuning the LM for label prediction, ℒq=pϕ(x,y)[logqθ(y∣x)]=pϕ(x,y)[∑n=1Nlogqθ(yn∣on)]≈∑n=1Nlogqθ(yn∣hn(L)(o;ϕ¯)), splitL_q&=E_p_φ(x,y) [ q_θ(y x) ]\\ &=E_p_φ(x,y) [ _n=1^N q_θ\! (y_n o_n ) ]\\ &≈ _n=1^N q_θ\! (y_n h_n^(L)(o; φ) ), split (4) where hn(L)(o;ϕ¯)h_n^(L)(o; φ) denotes the deterministic GNN representation computed from the LM embeddings o under the frozen GNN parameters ϕ¯ φ. Since expectation over pϕ(x,y)p_φ(x,y) is intractable, we estimate it by a single-sample Monte Carlo sampling. In the E-step, the LM first generates embeddings ono_n, which are passed through the frozen GNN to obtain hn(L)(o,ϕ¯)h_n^(L)(o, φ). Because the GNN parameters are fixed, the loss is computed on hn(L)(o,ϕ¯)h_n^(L)(o, φ) in place of ono_n. Gradients backpropagate through the fixed GNN to ono_n, updating only the LM parameters θ. 3.4 M-Step: GNN Optimization Note that logpϕ(x,y)=logpϕ(y∣x)+logpϕ(x) p_φ(x,y)= p_φ(y x)+ p_φ(x). In the M-step, we fix the LM and update the GNN to maximize the expected log-likelihood qθ(y∣x)[logpϕ(x,y)]=qθ(y∣x)[logpϕ(y∣x)]+logpϕ(x),E_q_θ(y x) [ p_φ(x,y) ]=E_q_θ(y x) [ p_φ(y x) ]+ p_φ(x), where the second term on the right hand side is independent of ϕφ since the GNN only parameterizes pϕ(y∣x)p_φ(y x), and is therefore omitted when optimizing ϕφ. Directly maximizing logpϕ(y∣x) p_φ(y x) is often intractable because the partition function requires summing over the exponentially large label space to normalize the distribution (Wainwright and Jordan, 2008). Following the pseudo-likelihood approximation (Besag, 1975), we instead maximize a product of tractable local conditional likelihoods, where each term yny_n is conditioned on the remaining labels y∖ny_N n and the full observed evidence x: qθ(y∣x)[∑nNlogpϕ(yn∣x,y∖n)]=qθ(y∣x)[∑nNlogpϕ(yn∣xn,yNB(n))],E_q_θ(y x) [ _n^N p_φ(y_n x,y_N n) ]=E_q_θ(y x) [ _n^N p_φ(y_n x_n,y_NB(n)) ], where the graphical Markov property implies that yny_n is conditionally independent of all non-neighbor labels given the evidence and its neighbors, i.e., pϕ(yn∣x,y∖n)=pϕ(yn∣on,yNB(n))p_φ(y_n x,y_N n)=p_φ(y_n o_n,y_NB(n)). We estimate the expectation under qθ(y∣x)q_θ(y x) using a single Monte Carlo sample, and represent the node xnx_n by its LM embedding ono_n. The M-step updates the GNN by minimizing the supervised cross-entropy loss using the LM-encoded representations, ℒp=∑nNlogpϕ(yn∣on,yNB(n)).L_p= _n^N p_φ(y_n o_n,y_NB(n)). (5) Input: Patient sequences x, target labels y, LM qθq_θ, GNN pϕp_φ, and patient graph A Output: LM embeddings onn=1N\o_n\_n=1^N, GNN node embeddings hn(L)n=1N\h_n^(L)\_n=1^N, and predictions y^nn=1N\ y_n\_n=1^N while not converged do // E-step: update the LM while fixing the GNN update LM qθq_θ by minimizing ℒqL_q in equation 4; // M-step: update the GNN while fixing the LM update GNN pϕp_φ by minimizing ℒpL_p in equation 5; Algorithm 1 VEM Optimization 3.5 Optimization Algorithm 1 summarizes our VEM algorithm, which alternates between updating the LM qθq_θ and the GNN pϕp_φ while fixing the other component. In the E-step, we freeze the GNN parameters ϕφ and update the LM by minimizing the loss in equation 4 through the GNN output hn(L)h_n^(L). Gradients update only the LM parameters θ. In the M-step, we freeze the LM parameters θ, use the LM representations ono_n as fixed node features, and update the GNN by minimizing the loss in equation 5. Gradients update only the GNN parameters ϕφ. During inference, the LM first encodes xnx_n into ono_n and initializes the node feature as hn(0)=onh_n^(0)=o_n, after which the GNN computes hn(L)h_n^(L) for prediction. At validation and test time, both modules are fixed and applied to the corresponding cohort graph for prediction. 3.6 Reference-Patient Attribution PLM provides reference-patient explanations by retrieving reference patients for a target patient i from the patient graph A. For each target patient, we retrieve up to 10 reference patients. When more than 10 neighbors are available, we select the top 10 based on cosine similarity of their embeddings. For each reference patient j, we quantify its contribution to the prediction for patient i via a gradient-based importance score, Imp(i,j)=∂logiti/∂AijImp(i,j)=∂\,logit_i/∂\,A_ij (Simonyan et al., 2014). This is defined as the gradient of the target logit with respect to the edge weight AijA_ij, which reflects the contribution of j. We rank reference patients by Imp(i,j)Imp(i,j) and report the top peers as explanations for prediction. To quantitatively assess the faithfulness of our reference-patient attribution, we perform a masking-based evaluation by removing selected reference edges and measuring the resulting prediction drop. 4 Experiments 4.1 Experimental Setting Datasets and Preprocessing. We use two large-scale, de-identified EHR datasets, MIMIC-I and MIMIC-IV (Johnson et al., 2016, 2023). Each patient is represented as a temporally ordered visit sequence, with diagnoses, procedures, and medications treated as unordered sets within each visit. Data statistics and preprocessing are deferred to Appendix C.1 and C.2. Tasks. We evaluate three tasks: readmission within 15 days (AUPRC/AUROC), length-of-stay prediction over 10 classes (AUPRC/F1), and drug recommendation (AUPRC/F1/Jaccard). We use BCE for readmission and drug recommendation, and CE for LOS. Detailed task settings and metric definitions are provided in Appendix C.3. Baselines. We compare against representative clinical prediction models, namely Deepr (Nguyen et al., 2017), RETAIN (Choi et al., 2016b), GRAM (Choi et al., 2017), StageNet (Gao et al., 2020), AdaCare (Ma et al., 2020), and GRASP (Zhang et al., 2021). We further compare with LM-graph models that incorporate external medical knowledge or graph, including G-BERT (Shang et al., 2019a), LEADER (Liu et al., 2024), GLEM (Zhao et al., 2023), GraphCare (Jiang et al., 2024), KARE (Jiang et al., 2025), ColaCare (Wang et al., 2025). For drug recommendation, we also include SafeDrug (Yang et al., 2021b), MICRON (Yang et al., 2021a), GAMENet (Shang et al., 2019b), MoleRec (Yang et al., 2023), UDC (Zhao et al., 2025), and KEHGCN (Zhang et al., 2026). Details about baseline implementations are provided in Appendix C.4. Table 1: Clinical prediction performance on MIMIC-I and IV. We report the mean (standard deviation) performance (%) over 10 runs. The best results are highlighted for both datasets, and the second-best results are underlined. N/A indicates not applicable. Model (%) Task 1: Readmission Task 2: LOS MIMIC-I MIMIC-IV MIMIC-I MIMIC-IV AUPRC AUROC AUPRC AUROC AUPRC F1 AUPRC F1 PLM (BCMBERT-396M) 51.7 ± 0.5 81.2 ± 0.4 48.6 ± 0.4 80.6 ± 0.4 83.6 ± 0.2 37.8 ± 0.4 82.8 ± 0.2 35.0 ± 0.3 PLM (Meerkat-8B) 52.2 ± 0.4 81.4 ± 0.4 49.7 ± 0.2 80.3 ± 0.3 84.5 ± 0.3 38.6 ± 0.2 85.3 ± 0.1 35.6 ± 0.2 Deepr 44.8 ± 0.9 76.1 ± 0.4 42.6 ± 0.2 74.7 ± 0.3 77.9 ± 0.1 35.0 ± 0.4 79.5 ± 0.3 32.3 ± 0.1 RETAIN 40.6 ± 1.0 71.0 ± 0.7 38.5 ± 0.4 69.5 ± 0.5 78.2 ± 0.1 34.9 ± 0.4 78.9 ± 0.3 32.0 ± 0.2 GRAM 43.8 ± 0.7 74.1 ± 0.4 41.4 ± 0.2 73.0 ± 0.3 78.2 ± 0.1 34.5 ± 0.2 78.8 ± 0.2 31.9 ± 0.3 StageNet 42.5 ± 0.6 73.8 ± 0.4 40.2 ± 0.1 72.6 ± 0.1 78.3 ± 0.2 34.4 ± 0.4 79.2 ± 0.3 31.3 ± 0.3 AdaCare 44.0 ± 0.6 74.1 ± 0.3 41.6 ± 0.1 72.9 ± 0.1 N/A N/A N/A N/A GRASP 43.1 ± 0.4 73.4 ± 0.6 40.9 ± 0.3 72.3 ± 0.2 N/A N/A N/A N/A G-BERT 49.6 ± 0.6 79.2 ± 0.8 46.3 ± 0.4 77.1 ± 0.7 79.7 ± 0.4 35.2 ± 0.3 78.6 ± 0.2 32.0 ± 0.1 LEADER 45.8 ± 0.5 76.3 ± 0.4 44.7 ± 0.2 75.3 ± 0.2 79.3 ± 0.3 35.4 ± 0.2 78.2 ± 0.2 31.5 ± 0.2 GLEM 45.7 ± 0.8 75.1 ± 0.9 42.6 ± 0.6 73.9 ± 0.7 78.5 ± 0.3 35.2 ± 0.4 79.0 ± 0.2 32.5 ± 0.3 GraphCare 48.8 ± 0.5 77.6 ± 0.6 45.0 ± 0.4 76.3 ± 0.5 81.4 ± 0.3 37.5 ± 0.5 81.7 ± 0.2 34.2 ± 0.2 KARE 49.5 ± 0.3 78.8 ± 0.3 48.1 ± 0.2 78.2 ± 0.3 80.7 ± 0.4 32.9 ± 0.2 81.3 ± 0.3 33.5 ± 0.2 ColaCare 50.1 ± 0.4 79.0 ± 0.3 47.5 ± 0.4 78.2 ± 0.3 82.3 ± 0.4 35.5 ± 0.4 82.5 ± 0.2 33.7 ± 0.2 Model (%) Task 3: Drug Recommendation MIMIC-I MIMIC-IV AUPRC F1 Jaccard AUPRC F1 Jaccard PLM (BCMBERT-396M) 79.6 ± 0.2 68.4 ± 0.3 51.5 ± 0.3 75.4 ± 0.1 64.8 ± 0.3 48.9 ± 0.3 PLM (Meerkat-8B) 80.3 ± 0.2 69.4 ± 0.3 53.0 ± 0.3 77.1 ± 0.2 65.7 ± 0.4 50.4 ± 0.4 Deepr 72.3 ± 0.1 60.3 ± 0.4 44.7 ± 0.3 63.7 ± 0.1 53.1 ± 0.4 39.8 ± 0.4 RETAIN 75.1 ± 0.3 65.2 ± 0.2 49.4 ± 0.2 65.7 ± 0.6 56.9 ± 0.4 41.5 ± 0.4 GRAM 74.7 ± 0.1 62.9 ± 0.3 47.9 ± 0.3 65.3 ± 0.2 53.1 ± 0.2 40.3 ± 0.3 StageNet 73.4 ± 0.1 61.4 ± 0.3 45.8 ± 0.4 63.1 ± 0.1 52.2 ± 0.3 37.5 ± 0.4 SafeDrug 75.8 ± 0.2 66.2 ± 0.2 50.5 ± 0.2 67.1 ± 0.3 58.2 ± 0.3 43.0 ± 0.2 MICRON 76.5 ± 0.3 67.4 ± 0.3 51.1 ± 0.2 66.6 ± 0.4 59.5 ± 0.3 44.1 ± 0.3 GAMENet 76.1 ± 0.1 66.0 ± 0.1 50.2 ± 0.2 67.2 ± 0.3 58.7 ± 0.3 43.4 ± 0.3 MoleRec 72.8 ± 0.1 61.1 ± 0.3 46.1 ± 0.3 62.3 ± 0.1 56.3 ± 0.4 41.2 ± 0.3 UDC 78.5 ± 0.3 67.5 ± 0.3 51.7 ± 0.2 71.7 ± 0.2 62.7 ± 0.3 47.8 ± 0.3 KEHGCN 78.8 ± 0.2 68.1 ± 0.3 52.2 ± 0.2 72.3 ± 0.4 64.5 ± 0.3 49.1 ± 0.3 G-BERT 69.0 ± 0.2 61.9 ± 0.2 45.8 ± 0.3 62.4 ± 0.3 51.7 ± 0.2 37.3 ± 0.2 LEADER 78.0 ± 0.2 67.4 ± 0.3 51.8 ± 0.1 71.2 ± 0.2 63.0 ± 0.2 47.8 ± 0.2 GLEM 73.9 ± 0.1 62.1 ± 0.4 46.1 ± 0.4 66.0 ± 0.5 54.5 ± 0.4 41.4 ± 0.3 GraphCare 78.5 ± 0.2 66.2 ± 0.3 49.8 ± 0.4 70.7 ± 0.5 60.4 ± 0.3 45.7 ± 0.4 4.2 Implementation Details To demonstrate that our method is backbone-agnostic and consistently improves performance across backbones, we instantiated the PLM framework with both encoder-only and decoder-only biomedical LMs. For encoder-only backbone, we used BioClinical ModernBERT-396M (BCMBERT) (Sounack et al., 2025) and fine-tuned the last six Transformer layers with a task-specific prediction head. For decoder-only backbone, we used a Llama-based LLM Meerkat-8B, obtained sequence embeddings via LLM2Vec (BehnamGhader et al., 2024), and fine-tuned the model using LoRA adapters on attention projection layers (Hu et al., 2022), while freezing the remaining parameters. We used a standard three-layer GCN on the patient graph. Details about model architecture and hyperparameters are provided in Appendix C.5. 5 Results Figure 2: Interpretability of PLM. (a) Ego-graph for reference-patient attribution given a diabetes patient. Nodes are positioned closer and colored darker for higher importance scores. Direct diabetes evidence is highlighted in red, indirectly related endocrine conditions and diabetes complications are also annotated. (b) Quantitative evaluation of importance scores. We assess explanation faithfulness using a masking-based evaluation, where more negative Δ indicates that the masked reference patients were more influential. Figure 3: Sensitivity analyses. (a) Sensitivity to the graph sparsification threshold τ. Increasing τ produces a sparser graph with lower average degree. The readmission prediction performance peaks around τ=8τ=8. (b) Drug recommendation performance stratified by visit number on MIMIC-IV. Performance remains stable for one to three visits and declines noticeably for more than three visits, indicating PLM does not rely on long visit histories. Table 2: Ablation studies on MIMIC-IV comparing PLM with four training strategies under the BCMBERT and Meerkat backbones. VEM consistently yields the best performance across tasks. Backbone Training Readmission LOS Drug Recommendation AUPRC AUROC AUPRC F1 AUPRC F1 Jaccard BCMBERT-396M VEM 48.6 80.6 82.8 35.0 75.4 64.8 48.9 Alternating 45.6 76.9 79.8 33.4 69.2 58.2 44.8 E2E 44.2 76.0 78.4 31.6 71.5 62.5 46.4 2-stage 45.5 77.4 78.7 32.0 72.4 62.8 46.3 LM-only 43.6 74.5 77.0 30.7 70.2 61.4 45.5 Meerkat-8B VEM 49.7 80.3 85.3 35.6 77.1 65.7 50.4 Alternating 46.4 77.6 81.2 33.1 72.4 60.9 47.1 E2E 44.4 77.2 78.0 30.8 70.4 60.2 45.4 2-stage 46.2 77.3 81.8 32.9 73.2 63.6 47.9 LM-only 43.9 75.4 79.1 30.9 71.3 60.7 45.8 5.1 Clinical Predictions on MIMIC Datasets As shown in Table 1, PLM achieves the strongest overall performance across datasets and tasks. Both PLM variants consistently outperform prior clinical prediction models across readmission, LOS, and drug recommendation, with the Meerkat variant performing best in most settings. Meanwhile, PLM with BCMBERT remains highly competitive and still surpasses the strongest prior baselines on MIMIC-IV. The gains are particularly consistent for drug recommendation, suggesting that combining patient-level semantics with cohort structure benefits multi-label clinical prediction. Table 7 reports runtime and peak GPU memory on MIMIC-IV with the Meerkat-8B backbone, showing that VEM achieves substantially better performance than all baselines. VEM outperforms alternating optimization while requiring close training time and GPU memory. Compared with LM-only, VEM adds only 4.1% runtime and 4.7% memory overhead, while remaining more efficient than E2E and two-stage training. Further details are provided in Section C.10. 5.2 Generalization Across LM Backbones To assess the generalization of VEM across backbones, we further evaluated PLM with four backbones on MIMIC-IV spanning both encoder-only models (BCMBERT and BioBERT) and decoder-only models (Meerkat and BioMistral). VEM consistently improves over each corresponding LM-only baseline across all three tasks, demonstrating that the proposed framework generalizes across different LM backbones rather than depending on a specific backbone. Detailed experimental results are provided in Appendix C.9. 5.3 Interpretability Analysis We assessed interpretability of PLM on insulin prediction in drug recommendation. Fig. 2.a shows an ego graph for a target diabetes patient 57139 with ten retrieved reference cases. PLM assigns each reference patient an importance score (Imp); higher scores denote closer distance to the target and darker color. We highlighted in red the CCS codes directly associated with insulin (e.g., 50: Diabetes mellitus with complications) and annotated indirectly related CCS codes such as endocrine conditions and diabetes complications. Patient 305 with a direct diabetes CCS code has the highest Imp 0.41. Patients 107 and 384 also have high scores (0.33 and 0.29) given their multiple diabetes-related comorbidities. Overall, the retrieved references support insulin prediction by highlighting both direct diabetes cases and cases with diabetes-related complications. We further compared against a baseline that ranks reference patients using only weighted edges AijA_ij (Fig. 5). This edge-weight baseline can miss clinically salient references. Despite being the most influential reference with a direct diabetes code, patient 305 has an edge weight of 9, only slightly above the neighborhood threshold of 8. Therefore, edge weights mainly reflect shared-code overlap, whereas PLM leverages both textual semantics and graph structure to yield relevant cases as explanations. Additional details and interpretability analyses are provided in Appendix C.7. To quantitatively evaluate the faithfulness of our reference-patient attributions, we conducted a masking-based evaluation on the readmission, LOS, and drug recommendation tasks. For each test patient, we removed the three retrieved reference edges selected by a given masking strategy, ran inference on the updated test graph, and measured the resulting change in AUPRC. We considered three selection strategies: top-3 by our importance scores, top-3 by edge weights, and random selection. Fig. 2.b shows that masking the top-3 references ranked by our importance scores consistently yields the largest performance drop across all three tasks. This suggests that the computed importance scores more faithfully identify the reference patients that drive the model’s predictions. 5.4 Sensitivity Analysis We study readmission sensitivity to the graph sparsification threshold τ, where an edge is retained only if Aij≥τA_ij≥τ. As shown in Figure 3.a, increasing τ prunes weaker links, lowers average degree, and improves performance up to τ=8τ=8. Beyond that, the graph becomes overly sparse and performance degrades. We therefore use τ=8τ=8, which balances efficiency and predictive performance. We also stratify MIMIC-IV drug recommendation performance by visit count. Figure 3.b shows stable performance for one to three visits and a clearer decline beyond three, suggesting that PLM does not rely on long histories for strong performance; this is consistent with prior work on increasing patient complexity (Meng et al., 2024). 5.5 Ablation Studies We ablated PLM under both BCMBERT and Meerkat backbones by comparing it with LM-only and LM-GNN variants trained using alternating, E2E or in two-stage methods. As shown in Table 2, PLM with VEM achieves the best performance across all tasks under both backbones. These results indicate that the improvements stem from VEM’s probabilistic training scheme, rather than simply adding a GNN, switching LM backbones, or using various optimization methods. 6 Conclusion We propose PLM, a unified LM-GNN framework for clinical prediction that integrates patient textual semantics with graph relational structure through a VEM algorithm. Across MIMIC-I and MIMIC-IV, PLM consistently improves readmission prediction, LOS forecasting, and drug recommendation, and these gains generalize across both encoder-only and decoder-only LM backbones. Compared with other LM-GNN training strategies, VEM achieves these improvements with only modest computational overhead. Qualitative case studies show that PLM identifies influential similar patients as reference-patient explanations, while quantitative edge masking confirms that the highest-ranked references have the strongest effect on predictions. Limitations PLM currently constructs patient graphs using diagnosis-code overlap, which may not capture richer clinical relationships available from multimodal EHR data or external knowledge. In addition, our reference-patient attributions are evaluated quantitatively for faithfulness, while further clinician assessment is needed before real-world deployment. Future work will explore richer graph construction and evaluate these attributions with clinical stakeholders. References E. Alsentzer, M. M. Li, S. N. Kobren, A. Noori, Undiagnosed Diseases Network, I. S. Kohane, and M. Zitnik (2025) Few shot learning for phenotype-driven diagnosis of patients with rare genetic diseases. NPJ Digit. Med. 8 (1), p. 380 (en). Cited by: §1, §2.1. P. BehnamGhader, V. Adlakha, M. Mosbach, D. Bahdanau, N. Chapados, and S. Reddy (2024) LLM2Vec: large language models are secretly powerful text encoders. In First Conference on Language Modeling, External Links: Link Cited by: §C.5, §4.2. J. Besag (1975) Statistical analysis of non-lattice data. Statistician 24 (3), p. 179–195. Cited by: §3.2, §3.4. D. M. Blei, A. Kucukelbir, and J. D. McAuliffe (2017) Variational inference: a review for statisticians. J. Am. Stat. Assoc. 112 (518), p. 859–877 (en). Cited by: §3.2. H. O. Boll, A. Amirahmadi, M. M. Ghazani, W. O. de Morais, E. P. de Freitas, A. Soliman, F. Etminani, S. Byttner, and M. R. Mendoza (2024) Graph neural networks for clinical risk prediction based on electronic health records: A survey. J. Biomed. Informatics 151, p. 104616. Cited by: §1, §2.2. J. Chen, C. Yin, Y. Wang, and P. Zhang (2024) Predictive modeling with temporal graphical representation on electronic health records. In IJCAI, p. 5763–5771. Cited by: §2.2. E. Choi, M. T. Bahadori, A. Schuetz, W. F. Stewart, and J. Sun (2016a) Doctor ai: predicting clinical events via recurrent neural networks. In Proceedings of the 1st Machine Learning for Healthcare Conference, Proceedings of Machine Learning Research, Vol. 56, p. 301–318. External Links: Link Cited by: §2.1. E. Choi, M. T. Bahadori, L. Song, W. F. Stewart, and J. Sun (2017) GRAM: graph-based attention model for healthcare representation learning. In KDD, p. 787–795. Cited by: §C.4, §2.1, §4.1. E. Choi, M. T. Bahadori, J. Sun, J. Kulas, A. Schuetz, and W. F. Stewart (2016b) RETAIN: an interpretable predictive model for healthcare using reverse time attention mechanism. In NIPS, p. 3504–3512. Cited by: §C.4, §2.1, §4.1. A. Defilippo, P. Veltri, P. Lió, and P. H. Guzzi (2024) Leveraging graph neural networks for supporting automatic triage of patients. CoRR abs/2403.07038. Cited by: §2.1. A. P. Dempster, N. M. Laird, and D. B. Rubin (1977) Maximum likelihood from incomplete data via the em algorithm. Journal of the Royal Statistical Society. Series B (Methodological) 39 (1), p. 1–38. External Links: ISSN 00359246, Link Cited by: §A.3. X. Dong, R. Wong, W. Lyu, K. Abell-Hart, J. Deng, Y. Liu, J. G. Hajagos, R. N. Rosenthal, C. Chen, and F. Wang (2023) An integrated lstm-heterorgnn model for interpretable opioid overdose risk prediction. Artificial intelligence in medicine 135, p. 102439. Cited by: §2.2. C. W. Fox and S. J. Roberts (2012) A tutorial on variational bayesian inference. Artif. Intell. Rev. 38 (2), p. 85–95 (en). Cited by: §3.2. J. Gao, C. Xiao, Y. Wang, W. Tang, L. M. Glass, and J. Sun (2020) StageNet: stage-aware neural networks for health risk prediction. In Proceedings of The Web Conference 2020, W ’20, New York, NY, USA, p. 530–540. External Links: ISBN 9781450370233, Link, Document Cited by: §C.4, §4.1. J. Gilmer, S. S. Schoenholz, P. F. Riley, O. Vinyals, and G. E. Dahl (2017) Neural message passing for quantum chemistry. In ICML, Proceedings of Machine Learning Research, Vol. 70, p. 1263–1272. Cited by: §2.2. S. N. Golmaei and X. Luo (2021) DeepNote-gnn: predicting hospital readmission using clinical notes and patient network. In BCB, p. 19:1–19:9. Cited by: §2.2. S. Hegselmann, G. von Arnim, T. Rheude, N. Kronenberg, D. Sontag, G. Hindricks, R. Eils, and B. Wild (2025) Large language models are powerful electronic health record encoders. External Links: 2502.17403, Link Cited by: §1. G. E. Hinton, P. Dayan, B. J. Frey, and R. M. Neal (1995) The “wake-sleep” algorithm for unsupervised neural networks. Science 268 (5214), p. 1158–1161 (en). Cited by: §A.1, §3.3. E. J. Hu, Y. Shen, P. Wallis, Z. Allen-Zhu, Y. Li, S. Wang, L. Wang, and W. Chen (2022) LoRA: low-rank adaptation of large language models. In International Conference on Learning Representations, External Links: Link Cited by: §4.2. S. Hu, G. Zou, S. Yang, S. Lin, Y. Gan, B. Zhang, and Y. Chen (2025) Large language model meets graph neural network in knowledge distillation. In Proceedings of the AAAI Conference on Artificial Intelligence, Vol. 39, p. 17295–17304. Cited by: §1, §2.2. K. Huang, J. Altosaar, and R. Ranganath (2019) ClinicalBERT: modeling clinical notes and predicting hospital readmission. CoRR abs/1904.05342. Cited by: §2.1. X. Huang, K. Han, Y. Yang, D. Bao, Q. Tao, Z. Chai, and Q. Zhu (2024) Can gnn be good adapter for llms?. In Proceedings of the ACM Web Conference 2024, p. 893–904. Cited by: §1. P. Jiang, C. Xiao, A. Cross, and J. Sun (2024) GraphCare: enhancing healthcare predictions with personalized knowledge graphs. In International Conference on Learning Representations, Cited by: §C.2, §C.2, §C.3, §C.4, §2.2, §4.1. P. Jiang, C. Xiao, M. Jiang, P. Bhatia, T. Kass-Hout, J. Sun, and J. Han (2025) Reasoning-enhanced healthcare predictions with knowledge graph community retrieval. In The Thirteenth International Conference on Learning Representations, External Links: Link Cited by: §C.4, §4.1. A. E. Johnson, L. Bulgarelli, L. Shen, A. Gayles, A. Shammout, S. Horng, T. J. Pollard, S. Hao, B. Moody, B. Gow, et al. (2023) MIMIC-iv, a freely accessible electronic health record dataset. Scientific data 10, p. 1. Cited by: §C.1, §4.1. A. E. Johnson, T. J. Pollard, L. Shen, L. H. Lehman, M. Feng, M. Ghassemi, B. Moody, P. Szolovits, L. Anthony Celi, and R. G. Mark (2016) MIMIC-i, a freely accessible critical care database. Scientific data 3, p. 1–9. Cited by: §C.1, §4.1. J. Kauffman, E. Holmes, A. Vaid, A. W. Charney, P. Kovatch, J. Lampert, A. Sakhuja, M. Zitnik, B. S. Glicksberg, I. Hofer, and G. N. Nadkarni (2025) InfEHR: clinical phenotype resolution through deep geometric learning on electronic health records. Nat. Commun. 16 (1), p. 8475 (en). Cited by: §1, §2.1. H. Kim, H. Hwang, J. Lee, S. Park, D. Kim, T. Lee, C. Yoon, J. Sohn, J. Park, O. Reykhart, T. Fetherston, D. Choi, S. H. Kwak, Q. Chen, and J. Kang (2025) Small language models learn enhanced reasoning skills from medical textbooks. NPJ Digit. Med. 8 (1), p. 240 (en). Cited by: §2.1. D. P. Kingma and M. Welling (2022) Auto-encoding variational bayes. External Links: 1312.6114, Link Cited by: §3.2, §3.3. T. N. Kipf and M. Welling (2017a) Semi-supervised classification with graph convolutional networks. In ICLR (Poster), Cited by: §2.2. T. N. Kipf and M. Welling (2017b) Semi-supervised classification with graph convolutional networks. In International Conference on Learning Representations (ICLR), Cited by: §C.5, §3.1. Y. Labrak, A. Bazoge, E. Morin, P. Gourraud, M. Rouvier, and R. Dufour (2024) BioMistral: A collection of open-source pretrained large language models for medical domains. In ACL (Findings), Findings of ACL, Vol. ACL 2024, p. 5848–5864. Cited by: §2.1. P. Langley (2000) Crafting papers on machine learning. In Proceedings of the 17th International Conference on Machine Learning (ICML 2000), P. Langley (Ed.), Stanford, CA, p. 1207–1216. Cited by: §C.10. J. Lee, W. Yoon, S. Kim, D. Kim, S. Kim, C. H. So, and J. Kang (2020) BioBERT: a pre-trained biomedical language representation model for biomedical text mining. Bioinform. 36 (4), p. 1234–1240. Cited by: §2.1. Y. Li, Z. Li, P. Wang, J. Li, X. Sun, H. Cheng, and J. X. Yu (2024) A survey of graph meets large language model: progress and future directions. In Proceedings of the Thirty-Third International Joint Conference on Artificial Intelligence, IJCAI ’24. External Links: ISBN 978-1-956792-04-1, Link, Document Cited by: §1. Q. Liu, X. Wu, X. Zhao, Y. Zhu, Z. Zhang, F. Tian, and Y. Zheng (2024) Large language model distilling medication recommendation model. arXiv preprint arXiv:2402.02803. Cited by: §C.4, §4.1. Y. Liu, S. Yang, Y. Zhang, C. Miao, Z. Nie, and J. Zhang (2021) Learning hierarchical review graph representations for recommendation. IEEE Transactions on Knowledge and Data Engineering 35, p. 658–671. Cited by: §2.2. H. Lu and S. Uddin (2021) A weighted patient network-based framework for predicting chronic diseases using graph neural networks. Sci. Rep. 11 (1), p. 22607 (en). Cited by: §C.6, §3.1. D. Luo, W. Cheng, D. Xu, W. Yu, B. Zong, H. Chen, and X. Zhang (2020) Parameterized explainer for graph neural network. In Advances in Neural Information Processing Systems, Vol. 33. Cited by: §2.2. L. Ma, J. Gao, Y. Wang, C. Zhang, J. Wang, W. Ruan, W. Tang, X. Gao, and X. Ma (2020) AdaCare: explainable clinical health status representation learning via scale-adaptive feature extraction and recalibration. In AAAI, p. 825–832. Cited by: §C.4, §2.1, §4.1. N. Makarov, M. Bordukova, P. Quengdaeng, D. Garger, R. Rodriguez-Esteban, F. Schmich, and M. P. Menden (2025) Large language models forecast patient health trajectories enabling digital twins. npj Digit. Medicine 8 (1). Cited by: §1. F. Meng, X. Li, X. Hou, M. Lu, and Y. Zhang (2024) MGRN: toward robust drug recommendation via multi-view gating retrieval network. Bioinformatics 40 (10) (en). Cited by: §5.4. R. Miotto, L. Li, B. A. Kidd, and J. T. Dudley (2016) Deep patient: an unsupervised representation to predict the future of patients from the electronic health records. Sci. Rep. 6 (1) (en). Cited by: §2.1. P. Nguyen, T. Tran, N. Wickramasinghe, and S. Venkatesh (2017) Deepr: A convolutional net for medical records. IEEE J. Biomed. Health Informatics 21 (1), p. 22–30. Cited by: §C.4, §4.1. T. Pham, T. Tran, D. Q. Phung, and S. Venkatesh (2016) DeepCare: A deep dynamic memory model for predictive medicine. In PAKDD (2), Lecture Notes in Computer Science, Vol. 9652, p. 30–41. Cited by: §2.1. L. Rasmy, Y. Xiang, Z. Xie, C. Tao, and D. Zhi (2021) Med-bert: pretrained contextualized embeddings on large-scale structured electronic health records for disease prediction. NPJ digital medicine 4, p. 86. Cited by: §2.1. E. Rocheteau, C. Tong, P. Veličković, N. Lane, and P. Liò (2021) Predicting patient outcomes with graph representation learning. arXiv preprint arXiv:2101.03940. Cited by: §2.2. J. Shang, T. Ma, C. Xiao, and J. Sun (2019a) Pre-training of graph augmented transformers for medication recommendation. In Proceedings of the Twenty-Eighth International Joint Conference on Artificial Intelligence, IJCAI-19, p. 5953–5959. External Links: Document, Link Cited by: §C.4, §4.1. J. Shang, C. Xiao, T. Ma, H. Li, and J. Sun (2019b) Gamenet: graph augmented memory networks for recommending medication combination. In proceedings of the AAAI Conference on Artificial Intelligence, Vol. 33, p. 1126–1133. Cited by: §C.4, §1, §4.1. K. Simonyan, A. Vedaldi, and A. Zisserman (2014) Deep inside convolutional networks: visualising image classification models and saliency maps. External Links: 1312.6034, Link Cited by: §3.6. Z. Song, Y. Hu, A. Verma, D. L. Buckeridge, and Y. Li (2022) Automatic phenotyping by a seed-guided topic model. In Proceedings of the 28th ACM SIGKDD Conference on Knowledge Discovery and Data Mining, KDD ’22, New York, NY, USA, p. 4713–4723. External Links: ISBN 9781450393850, Link, Document Cited by: §2.1. T. Sounack, J. Davis, B. Durieux, A. Chaffin, T. J. Pollard, E. Lehman, A. E. W. Johnson, M. McDermott, T. Naumann, and C. Lindvall (2025) BioClinical modernbert: a state-of-the-art long-context encoder for biomedical and clinical nlp. External Links: 2506.10896, Link Cited by: §2.1, §4.2. M. Vaida and Z. Huang (2025) Multimodal graph neural networks in healthcare: a review of fusion strategies across biomedical domains. Frontiers in Artificial Intelligence 8, p. 1716706. Cited by: §2.2. M. J. Wainwright and M. I. Jordan (2008) Graphical models, exponential families, and variational inference. Found. Trends® Mach. Learn. 1 (1-2), p. 1–305 (en). Cited by: §3.4. J. Wang, J. Luo, M. Ye, X. Wang, Y. Zhong, A. Chang, G. Huang, Z. Yin, C. Xiao, J. Sun, and F. Ma (2024) Recent advances in predictive modeling with electronic health records. In IJCAI, p. 8272–8280. Cited by: §1. Z. Wang, Y. Zhu, H. Zhao, X. Zheng, D. Sui, T. Wang, W. Tang, Y. Wang, E. Harrison, C. Pan, J. Gao, and L. Ma (2025) ColaCare: enhancing electronic health record modeling through large language model-driven multi-agent collaboration. In Proceedings of the ACM on Web Conference 2025, W ’25, New York, NY, USA, p. 2250–2261. External Links: ISBN 9798400712746, Link, Document Cited by: §C.4, §4.1. M. Wornow, S. Bedi, M. A. F. Hernandez, E. Steinberg, J. A. Fries, C. Ré, S. Koyejo, and N. Shah (2025) Context clues: evaluating long context models for clinical prediction tasks on EHR data. In International Conference on Learning Representations, Cited by: §1. J. Xu, X. Xi, J. Chen, V. S. Sheng, J. Ma, and Z. Cui (2022) A survey of deep learning for electronic health records. Applied Sciences 12 (22), p. 11709. Cited by: §1, §2.1. H. Yan, C. Li, R. Long, C. Yan, J. Zhao, W. Zhuang, J. Yin, P. Zhang, W. Han, H. Sun, W. Deng, Q. Zhang, L. Sun, X. Xie, and S. Wang (2023) A comprehensive study on text-attributed graphs: benchmarking and rethinking. In NeurIPS, Cited by: §1, §2.2. C. Yang, C. Xiao, L. Glass, and J. Sun (2021a) Change matters: medication change prediction with recurrent residual networks. In Proceedings of the Thirtieth International Joint Conference on Artificial Intelligence, IJCAI-21, Z. Zhou (Ed.), p. 3728–3734. Note: Main Track External Links: Document, Link Cited by: §C.4, §4.1. C. Yang, C. Xiao, F. Ma, L. Glass, and J. Sun (2021b) SafeDrug: dual molecular graph encoders for recommending effective and safe drug combinations. In Proceedings of the Thirtieth International Joint Conference on Artificial Intelligence, IJCAI-21, Z. Zhou (Ed.), p. 3735–3741. Note: Main Track External Links: Document, Link Cited by: §C.4, §1, §4.1. N. Yang, K. Zeng, Q. Wu, and J. Yan (2023) Molerec: combinatorial drug recommendation with substructure-aware molecular representation learning. In Proceedings of the ACM web conference 2023, p. 4075–4085. Cited by: §C.4, §4.1. R. Ying, D. Bourgeois, J. You, M. Zitnik, and J. Leskovec (2019) GNNExplainer: generating explanations for graph neural networks. In Proceedings of the 33rd International Conference on Neural Information Processing Systems, Red Hook, NY, USA, p. 3594–3603. Cited by: §1, §2.2. H. Yuan, H. Yu, J. Wang, K. Li, and S. Ji (2021) On explainability of graph neural networks via subgraph explorations. In Proceedings of the 38th International Conference on Machine Learning, M. Meila and T. Zhang (Eds.), Proceedings of Machine Learning Research, Vol. 139, p. 12241–12252. Cited by: §2.2. C. Zhang, X. Gao, L. Ma, Y. Wang, J. Wang, and W. Tang (2021) GRASP: generic framework for health status representation learning based on incorporating knowledge from similar patients. In Proceedings of the AAAI conference on artificial intelligence, Vol. 35, p. 715–723. Cited by: §C.4, §4.1. Z. Zhang, H. Liu, X. Guo, T. Sun, and Z. Wu (2026) Knowledge-enhanced explainable HyperGraph convolution network for medication recommendation. Proc. Conf. AAAI Artif. Intell. 40 (19), p. 16424–16432. Cited by: §4.1. C. Zhao, H. Tang, J. Zhang, and X. Li (2025) Unveiling discrete clues: superior healthcare predictions for rare diseases. In Proceedings of the ACM on Web Conference 2025, W ’25, New York, NY, USA, p. 1747–1758. External Links: ISBN 9798400712746, Link, Document Cited by: §C.4, §4.1. J. Zhao, M. Qu, C. Li, H. Yan, Q. Liu, R. Li, X. Xie, and J. Tang (2023) Learning on large-scale text-attributed graphs via variational inference. In International Conference on Learning Representations, Cited by: §C.4, §1, §2.2, §4.1. Appendix A Detailed Mathematical Derivations A.1 Derivation of Variational EM Algorithm We define the ELBO objective of the VI as ℒ(θ):=qθ(y∣x)[logpϕ(x,y)−logqθ(y∣x)].L(θ):=E_q_θ(y x) [ p_φ(x,y)- q_θ(y x) ]. For the variational distribution qθ(y∣x)q_θ(y x), we adopt the mean-field assumption by factorizing latent variable z independently: qθ(y∣x)=∏n=1Nqθ(yn∣xn).q_θ(y x)= _n=1^Nq_θ(y_n x_n). Formally, maximizing the ELBO with respect to the variational distribution qθ(y∣x)q_θ(y x) is equivalent to minimizing the KL divergence between the posterior distribution and the variational distribution. The KL divergence is defined as KL[qθ(y∣x)||pϕ(y∣x)] [q_θ(y x)||p_φ(y x) ] =qθ(y∣x)[logqθ(y∣x)−logpϕ(y∣x)] =E_q_θ(y x) [ q_θ(y x)- p_φ(y x) ] =∑yqθ(y∣x)logqθ(y∣x)pϕ(y∣x). = _yq_θ(y x) q_θ(y x)p_φ(y x). Therefore, in the E-step, we update the variational distribution qθ(y∣x)q_θ(y x) (i.e., fine-tuning the LM) by minimizing KL[qθ(y∣x)||pϕ(y∣x)]KL [q_θ(y x)||p_φ(y x) ]. However, this objective is typically intractable because it requires to compute the expectations under qθ(⋅)q_θ(·) over all possible values. To address this, we use the wake-sleep algorithm (Hinton et al., 1995) to minimize the reverse KL divergence, where the sleep phase minimizes the KL divergence the wrong way round: θ∗ θ^* :=argminθKL[pϕ(y∣x)∥qθ(y∣x)] = _θKL [p_φ(y x)\|q_θ(y x) ] =argminθpϕ(y∣x)[logpϕ(y∣x) = _θE_p_φ(y x) [ p_φ(y x) −logqθ(y∣x)] - q_θ(y x) ] =argmaxθpϕ(y∣x)[logqθ(y∣x)]+const = _θE_p_φ(y x) [ q_θ(y x) ]+const =argmaxθpϕ(x,y)[logqθ(y∣x)], = _θE_p_φ(x,y) [ q_θ(y x) ], where const=pϕ(y∣x)[logpϕ(y∣x)]const=E_p_φ(y x) [ p_φ(y x) ] does not depend on θ and is therefore omitted during optimization. ℒq _q =pϕ(x,y)[logqθ(y∣x)] =E_p_φ(x,y) [ q_θ(y x) ] =pϕ(x,y)[∑nNlogqθ(yn∣xn)] =E_p_φ(x,y) [ _n^N q_θ(y_n x_n) ] =pϕ(x,y)[∑nNlogqθ(yn∣on)] =E_p_φ(x,y) [ _n^N q_θ(y_n o_n) ] ≈∑nNlogqθ(yn∣hn(L)(o,ϕ¯)), ≈ _n^N q_θ (y_n h_n^(L)(o, φ) ), (6) where hn(L)(o;ϕ¯)h_n^(L)(o; φ) denotes the deterministic GNN representation obtained from the LM embeddings o under frozen GNN parameters ϕ¯ φ. Since the expectation over pϕ(x,y)p_φ(x,y) is intractable, we estimate it with a single-sample Monte-Carlo sampling. In the E-step, LM first generates embeddings ono_n, which are passed through the frozen GNN to obtain hn(L)h_n^(L). Because the GNN parameters are fixed, the loss is computed on hn(L)(o,ϕ¯)h_n^(L)(o, φ) in place of ono_n. Gradients then backpropagate through the fixed GNN to ono_n, updating only the LM parameters θ. Thus, optimizing ℒqL_q amounts to maximizing the label log-likelihood, which we approximate by replacing the LM embedding ono_n with the fixed GNN representation hn(L)(o,ϕ¯)h_n^(L)(o, φ). In the M-step, we then update the GNN model pϕ(y∣x)p_φ(y x) towards maximizing the following pseudo-likelihood: qθ(y∣x)[logp(x,y)] _q_θ(y x) [ p(x,y) ] =Eqθ(y∣x)[logpϕ(y∣x)+logp(x)] =E_q_θ(y x) [ p_φ(y x)+ p(x) ] =qθ(y∣x)[logpϕ(y∣x)]+const =E_q_θ(y x) [ p_φ(y x) ]+const ≈qθ(y∣x)[∑nNlogpϕ(yn∣x,y∖n)] _q_θ(y x) [ _n^N p_φ(y_n x,y_N n) ] (pseudolikelihood) (pseudolikelihood) =qθ(y∣x)[∑nNlogpϕ(yn∣xn,y∖n)] =E_q_θ(y x) [ _n^N p_φ(y_n x_n,y_N n) ] =∑n=1Nlogpϕ(yn∣on,y∖n) = _n=1^N p_φ(y_n o_n,y_N n) =∑n=1Nlogpϕ(yn∣on,yNB(n)). = _n=1^N p_φ(y_n o_n,y_NB(n)). (7) The details of equation estimated as pseudo likelihood in that form is because that for pϕ(y∣x)p_φ(y x), there exist node and edge potentials ψn(⋅∣x),ψnm(⋅,⋅∣x) _n(· x), _nm(·,· x) such that pϕ(y∣x)=1Zϕ(x) p_φ(y x)= 1Z_φ(x) ∏n=1Nψn(yn∣xd) _n=1^N _n (y_n x_d ) ×∏(n,m)∈Eψnm(yn,ym∣xn). × _(n,m)∈ E _nm (y_n,y_m x_n ). We also define that for all the nodes x=(x1,…,xN)x= (x_1,…,x_N ) with labels y=(y1,…,yN)y= (y_1,…,y_N ) on a graph with neighbor sets NB(n)NB(n), then by the graphical Markov property we have: yn⟂y\(n∪NB(n))∣(x,yNB(n)),n=1,…,N.y_n y_N (\n\ (n)) (x,y_NB(n) ),n=1,…,N. Given the features x and the labels of n-th node’s neighbors yNB(n)y_NB(n), the label yny_n is statistically independent of all other non-neighbor labels. Thus, we have: pϕ(y∣x)≈∏n=1Npϕ(yn∣xn,y∖n) p_φ(y x)≈ _n=1^Np_φ (y_n x_n,y_N n ) =∏n=1Npϕ(yn,y∖n∣xn)∑n′pϕ(yn′,y∖n∣xn) = _n=1^N p_φ (y_n,y_N n x_n ) _n p_φ (y_n ,y_N n x_n ) =∏n=1Nψn(yn∣xn)∏m∈NB(n)ψnm(yn,ym∣xn)∑n′ψn(yn′∣xn′)∏m∈NB(n′)ψnm(yn′,ym∣xn′) = 372.91345pt$ _n=1^N _n (y_n x_n )Π _m (n) _nm (y_n,y_m x_n )Σ _n _n (y_n x_n )Π _m (n ) _nm (y_n ,y_m x_n )$ =∏n=1Npϕ(yn∣xn,yNB(n)). = _n=1^Np_φ (y_n x_n,y_NB(n) ). (8) The first line follows from the definition of the pseudo-likelihood. The transition from the second line to the third line applies the definition of the full conditional in fractional form, where we cancel the global normalizing constant Zϕ(x)Z_φ(x). The transition from the third line to the fourth line uses the Graphical Markov property. For the expectation qθ(y∣x)[⋅]E_q_θ(y x)[·], we approximate it using Monte-Carlo Sampling with a single sample, replacing the discrete input xnx_n with the LM-encoded embedding ono_n. As a result, pϕ(yn∣on,y∖n)p_φ(y_n o_n,y_N n) represents that GNN makes prediction given the LM-encoded embedding ono_n and all surrounding nodes y∖ny_N n. As a result, M-step will train a GNN using the following supervised objective on the labelled nodes: O(p)=∑n=1Vlogp(yn∣on,y∖n). O(p)= _n=1^V p(y_n o_n,y_N n). (9) A.2 Optimal Mean-Field Variational Update Theorem 4.1 Given the fixed pϕ(yn∣xn)p_φ(y_n x_n), the local optimum upate qθ∗(yn∣xn)q_θ^*(y_n x_n) satisfies: logqθ∗(yn∣xn) q_θ^*(y_n x_n) =qθ(y∣x)[logpϕ(yn∣xn,yNB(n))]. =E_q_θ(y x) [ p_φ(y_n x_n,y_NB(n)) ]. (10) Proof The goal of E-step is to optimize qθ(y∣x)q_θ(y x) by minimizing KL[qθ(y∣x)||pϕ(y∣x)]KL [q_θ(y x)||p_φ(y x) ]. Therefore, the objective function for qθ(y∣x)q_θ(y x) could be formulated as follows: ℒ(qθ(yn∣xn)) (q_θ(y_n x_n)) =−KL[qθ(y∣x)∥pϕ(y∣x)] =-KL [q_θ(y x)\|p_φ(y x) ] =qθ(y∣x)[logpϕ(y∣x)−logqθ(y∣x)] =E_q_θ(y x) [ p_φ(y x)- q_θ(y x) ] =qθ(y∣x)[logpϕ(y∣x) =E_q_θ(y x) [ p_φ(y x) −∑n′=1Nlogqθ(yn′∣xn′)] - _n =1^N q_θ(y_n x_n ) ] =qθ(y∣x)[logpϕ(y∣x) =E_q_θ(y x) [ p_φ(y x) −∑n′≠nNlogqθ(yn′∣xn′) - _n ≠ n^N q_θ(y_n x_n ) −logqθ(yn∣xn)] - q_θ(y_n x_n) ] =qθ(y∣x)[logpϕ(y∣x) =E_q_θ(y x) [ p_φ(y x) −logqθ(yn∣xn)]+const - q_θ(y_n x_n) ]+const =qθ(y∣x)[logpϕ(yn∣x,y∖n) =E_q_θ(y x) [ p_φ(y_n x,y_N n) +logpϕ(y∖n∣x) + p_φ(y_N n x) −logqθ(yn∣xn)]+const - q_θ(y_n x_n) ]+const =qθ(y∣x)[logpϕ(yn∣x,y∖n) =E_q_θ(y x) [ p_φ(y_n x,y_N n) −logqθ(yn∣xn)]+const - q_θ(y_n x_n) ]+const =qθ(y∣x)[logpϕ(yn∣xn,yNB(n)) =E_q_θ(y x) [ p_φ(y_n x_n,y_NB(n)) −logqθ(yn∣xn)]+const - q_θ(y_n x_n) ]+const =−KL[qθ(yn∣xn)∥pϕ(yn∣xn,yNB(n))] =-KL [q_θ(y_n x_n)\|p_φ(y_n x_n,y_NB(n)) ] +const. +const. (11) In GNNs, we leverage the Markov property that node n is conditionally independent of all non-neighbors given its neighbors, reducing y∖ny_N n to yNB(n)y_NB(n). As a result, the objective O(qθ(yn,∣xn))O(q_θ(y_n, x_n)) becomes the is the KL divergence between qθ(yn∣xn)q_θ(y_n x_n) and pϕ(yn∣xn,yNB(n))p_φ(y_n x_n,y_NB(n)). The local optimum satisfies qθ∗(yn∣xn)∝qθ(y∣x)[pϕ(yn∣xn,yNB(n))]q_θ^*(y_n x_n) _q_θ(y x) [p_φ(y_n x_n,y_NB(n)) ]. With Monte Carlo sampling, this local optimum becomes qθ∗(yn∣xn)∝pϕ(yn∣on,yNB(n))q_θ^*(y_n x_n) p_φ(y_n o_n,y_NB(n)), where ono_n is the LM‑encoded embedding for node n. This indicates that optimal LM learning corresponds to using the GNN’s predictions for yny_n based on LM‑encoded embeddings ono_n and neighbor labels yNB(n)y_NB(n), integrating local textual information and global patient relationships into clinical prediction. A.3 Comparison with Classical EM and GLEM Under Full Supervision The classical EM algorithm maximizes an observed-data likelihood by introducing an associated complete-data model. In the classical EM formulation (Dempster et al., 1977), the observed data are fixed, whereas the complete data include components that are not directly observed. The E-step computes the conditional expectation of the complete-data log-likelihood given the observed data and the current parameter estimate, and the M-step maximizes this expected complete-data log-likelihood. Thus, classical EM is tied to an observed-data likelihood and a particular complete-data augmentation. Our proposed VEM algorithm is related to this EM algorithm only at the level of its alternating inference-prediction structure. It is not an incomplete-label EM algorithm in the classical sense: all training labels are observed, and we do not impute missing labels or recover unobserved components of the training data. Instead, PLM introduces an auxiliary predictive distribution to couple an LM-based semantic encoder with a GNN-based structured predictor. This distinction is important because the variational object in PLM is not a posterior distribution over missing ground-truth labels. Degeneracy of GLEM under full supervision. GLEM is formulated for semi-supervised node classification on text-attributed graphs. Let V=L∪UV=L∪ U denote the set of all nodes, where L is the labelled set and U is the unlabelled set. The observed quantities are the node texts or semantic features sVs_V, the graph structure A, and the labels yLy_L of the labelled nodes. The labels yUy_U of the unlabelled nodes are treated as latent variables. A generic GLEM-style semi-supervised variational objective can be written as ℒGLEM(θ,ϕ) _GLEM(θ,φ) =qθ(yU∣sU)[logpϕ(yL,yU∣sV,A) =E_q_θ(y_U s_U) [ p_φ(y_L,y_U s_V,A) −logqθ(yU∣sU)], - q_θ(y_U s_U) ], where qθ(yU∣sU)q_θ(y_U s_U) is the LM-induced variational distribution over the missing labels and pϕ(yL,yU∣sV,A)p_φ(y_L,y_U s_V,A) is the GNN-induced graph-structured predictive model. Proposition A.1 (Collapse of GLEM’s missing-label variational objective). Suppose all nodes are labelled, so that U=∅U= and L=VL=V. Then the GLEM-style variational objective reduces to ℒGLEM(θ,ϕ)=logpϕ(yV∣sV,A).L_GLEM(θ,φ)= p_φ(y_V s_V,A). Consequently, the variational distribution over missing labels (i.e., latent variables) disappears, and the objective contains no nontrivial missing-label inference problem. Proof. If all nodes are labelled, then U=∅U= and hence yU=∅y_U= . The variational distribution over the empty label vector is the degenerate distribution qθ(yU∣sU)=qθ(∅∣s∅)=1.q_θ(y_U s_U)=q_θ( s_ )=1. Therefore, −qθ(yU∣sU)logqθ(yU∣sU)=−log1=0.-E_q_θ(y_U s_U) q_θ(y_U s_U)=- 1=0. Moreover, the expectation with respect to qθ(yU∣sU)q_θ(y_U s_U) is an expectation over a point mass, so qθ(yU∣sU)[logpϕ(yL,yU∣sV,A)] _q_θ(y_U s_U) [ p_φ(y_L,y_U s_V,A) ] =logpϕ(yL∣sV,A). = p_φ(y_L s_V,A). Since L=VL=V, this becomes ℒGLEM(θ,ϕ)=logpϕ(yV∣sV,A).L_GLEM(θ,φ)= p_φ(y_V s_V,A). Thus, the missing-label latent-variable component of the GLEM objective vanishes under full supervision. ∎ Proposition A.1 shows that the variational interpretation of GLEM relies on the existence of unlabelled-node labels. When U=∅U= , the latent space is |U|=0=∅,Y^|U|=Y^0=\ \, a singleton. Hence, there is no nontrivial E-step left for inferring missing labels. Any remaining alternating update is then an implementation choice for training LM and GNN components, rather than a non-degenerate missing-label variational EM procedure. Non-degeneracy of PLM under full supervision. PLM differs from GLEM in the role of the variational distribution. In supervised prediction, PLM treats all ground-truth labels as latent variables, while the latent variables in GLEM are vanished. We denote the ground-truth labels by ygt=(y1gt,…,yNgt).y^gt=(y^gt_1,…,y^gt_N). These labels enter the supervised losses. Separately, PLM introduces an auxiliary predictive label vector y~=(y~1,…,y~N), y=( y_1,…, y_N), which is used only to couple the LM and GNN components. The LM parameterizes a factorized auxiliary predictive distribution qθ(y~∣x)=∏n=1Nqθ(y~n∣xn),q_θ( y x)= _n=1^Nq_θ( y_n x_n), where x=xnn=1Nx=\x_n\_n=1^N denotes the patient sequences. The GNN parameterizes a graph-structured predictive distribution pϕ(y~∣x,A),p_φ( y x,A), where A is the patient graph. A supervised PLM objective may be written schematically as PLM(θ,ϕ) _PLM(θ,φ) =ℒsupLM(θ;ygt)+ℒsupGNN(ϕ;ygt) =L_sup^LM(θ;y^gt)+L_sup^GNN(φ;y^gt) +γℒcoup(θ,ϕ), +γ\,L_coup(θ,φ), where γ≥0γ≥ 0 controls the strength of the LM–GNN coupling. The coupling term is ℒcoup(θ,ϕ) _coup(θ,φ) =qθ(y~∣x)[logpϕ(y~∣x,A)−logqθ(y~∣x)]. =E_q_θ( y x) [ p_φ( y x,A)- q_θ( y x) ]. Thus, ygty^gt denotes the ground-truth labels, while y~ y denotes auxiliary predictive labels used for variational coupling. The latter should not be interpreted as missing ground-truth. When exact maximization of the graph-structured likelihood is intractable, PLM may replace the graph likelihood with a pseudo-likelihood approximation, leading to ℒcoupPLM(θ,ϕ) _coup^PLM(θ,φ) =qθ(y~∣x)[∑n=1Nlogpϕ(y~n∣xn,y~(n),A) =E_q_θ( y x) [ _n=1^N p_φ ( y_n x_n, y_N(n),A ) −∑n=1Nlogqθ(y~n∣xn)], - _n=1^N q_θ( y_n x_n) ], where (n)N(n) denotes the neighbours of node n in the patient graph. Proposition A.2 (Non-collapse of the PLM auxiliary variational space). Assume that qθ(y~∣x)q_θ( y x) is an auxiliary predictive distribution used to couple the LM and GNN components, rather than a posterior distribution over unobserved ground-truth labels. Then full supervision alone does not make the PLM variational space collapse. In particular, the PLM coupling objective reduces to ordinary supervised likelihood only if qθ(y~∣x)q_θ( y x) is additionally constrained to be a point mass at the observed labels. Proof. In GLEM, the variational distribution is defined only over the missing labels yUy_U of unlabelled nodes. Thus, when U=∅U= , the variational space collapses to a singleton. In PLM, by contrast, PLM defines auxiliary predictive distribution over y~=(y~1,…,y~N)∈N. y=( y_1,…, y_N) ^N. This space remains nontrivial whenever N>0N>0 and ||>1|Y|>1. The fact that the ground-truth labels yobsy^obs are observed does not by itself imply that the auxiliary predictive distribution qθ(y~∣x)q_θ( y x) is degenerate. Consequently, the entropy term −qθ(y~∣x)logqθ(y~∣x)-E_q_θ( y x) q_θ( y x) is generally nonzero, and the coupling term qθ(y~∣x)logpϕ(y~∣x,A)E_q_θ( y x) p_φ( y x,A) continues to depend jointly on the LM-induced predictive distribution and the GNN-induced structured predictive distribution. Hence, PLM’s objective contains a nontrivial LM-GNN alignment term even under full supervision. Only under the additional degenerate constraint qθ(y~∣x)=δygt(y~)q_θ( y x)= _y^gt( y) does the coupling expectation collapse to ℒcoup(θ,ϕ)=logpϕ(ygt∣x,A),L_coup(θ,φ)= p_φ(y^gt x,A), up to the zero entropy of the point mass. Therefore, PLM reduces to an ordinary supervised graph-predictive likelihood only after imposing this additional degeneracy condition on the auxiliary predictive distribution, not merely because all labels are observed. ∎ The distinction can therefore be summarized as follows. In GLEM, the variational space is |U|Y^|U|, which collapses to the singleton ∅\ \ under full supervision. In PLM, the auxiliary predictive space is NY^N, which remains nontrivial under full supervision. Thus, full supervision removes the missing-label latent variables in GLEM, but it does not remove the auxiliary predictive distribution used by PLM to align the LM and GNN. Accordingly, PLM should be viewed as a supervised VEM-inspired coupling framework rather than as a classical incomplete-data EM algorithm. Its novelty is not the use of EM itself, but the construction of a supervised variational coupling objective that remains meaningful when all training labels are observed. Table 3: Comparison between Patients-like-me and GLEM. Aspect Patients-like-me GLEM Objective Supervised prediction of all labels y Semi-supervised completion of missing labels yUy_U Observed variables Sequences x Node text features sVs_V and labels on labeled nodes yLy_L Latent variables Labels y Missing labels yUy_U ELBO qθ(y∣x)[logpϕ(x,y)−logqθ(y∣x)]E_q_θ(y x) [ p_φ(x,y)- q_θ(y x) ] qθ(yU∣sU)[logpϕ(yL,yU∣sV)−logqθ(yU∣sU)]E_q_θ(y_U s_U) [ p_φ(y_L,y_U s_V)- q_θ(y_U s_U) ] E-step Generate patient embedding via LM qθ(y∣x)q_θ(y x) guided by graph structural Fit LM qθ(yU∣sU)q_θ(y_U s_U) given GNN-imputed missing labels M-step Predict targets via GNN pϕ(y∣x)p_φ(y x) using LM embeddings Fit GNN pϕp_φ using LM-imputed missing labels Under full supervision EM maintained via latent targets y EM collapse to joint training due to empty latent variables Appendix B Details of Patients-Like-Me B.1 Prompt Template Prompt Template for Readmission Prediction Demographics: Insurance: <INSURANCE>; Language: <LANGUAGE>; Religion: <RELIGION>; Marital status: <MARITAL_STATUS>; Ethnicity: <ETHNICITY>. The patient has <VISIT_NUM> ICU visits. In visit 1, diagnoses: <DIAG_NAME>, …; procedures: <PROC_NAME>, …; medications: <MED_NAME>, …. … Current visit: diagnoses: <DIAG_NAME>, …; procedures: <PROC_NAME>, medications: <MED_NAME> …. Task: Predict whether the patient will be readmitted within 15 days after the current visit. Choose exactly one class from 0,1, where 0 = no readmission within 15 days and 1 = readmission within 15 days. Output format: A single class label in 0,1. Prompt Template for Length-of-Stay Prediction Demographics: Insurance: <INSURANCE>; Language: <LANGUAGE>; Religion: <RELIGION>; Marital status: <MARITAL_STATUS>; Ethnicity: <ETHNICITY>. The patient has <VISIT_NUM> ICU visits. In visit 1, diagnoses: <DIAG_NAME>, …; procedures: <PROC_NAME>, …; medications: <MED_NAME>, …. … Current visit: diagnoses: <DIAG_NAME>, …; procedures: <PROC_NAME>, …; medications: <MED_NAME>, …. Task: Predict the ICU length-of-stay (LOS) category for the current visit. Choose exactly one class from 0,1,2,3,4,5,6,7,8,9, where: 0: LOS << 1 day; 1: LOS = 1 day; 2: LOS = 2 days; 3: LOS = 3 days; 4: LOS = 4 days; 5: LOS = 5 days; 6: LOS = 6 days; 7: LOS = 7 days; 8: LOS = 8-14 days; 9: LOS ≥ 15 days. Output format: A single class label in 0,1,2,3,4,5,6,7,8,9. Prompt Template for Drug Recommendation Demographics: Insurance: <INSURANCE>; Language: <LANGUAGE>; Religion: <RELIGION>; Marital status: <MARITAL_STATUS>; Ethnicity: <ETHNICITY>. The patient has <VISIT_NUM> ICU visits. In visit 1, diagnoses: <DIAG_NAME>, …; procedures: <PROC_NAME>, …; medications: <MED_NAME>, …. … Current visit: diagnoses: <DIAG_NAME>, …; procedures: <PROC_NAME>, …. Task: Predict the medications to prescribe for the current visit. Output a multi-label set of 3-level ATC codes. Choose any number of codes from the medication vocabulary. Note that a patient may have multiple medications. Output format: A multi-label set of ATC level-3 codes. Appendix C Additional Experimental Details C.1 Datasets and Cohort Statistics We used two large-scale, de-identified EHR datasets MIMIC-I and MIMIC-IV from Beth Israel Deaconess Medical Center (BIDMC), released via PhysioNet under a credentialed data use agreement (Johnson et al., 2016, 2023). Both datasets contain inpatient EHR records, including diagnoses, procedures, medication prescriptions, and admission/discharge timestamps. Patient identifiers are removed and all timestamps are shifted by a patient-specific offset. This preserves relative temporal ordering within each patient sequence while preventing alignment of absolute calendar time across patients. We construct visit-based patient sequences by treating each inpatient admission as a visit. Table 4 summarizes the basic cohort statistics after preprocessing. Table 4: Basic statistics of the preprocessed MIMIC-I and MIMIC-IV datasets. “#” denotes counts and “/patient” denotes the average number of records per-patient. Dataset #patients #visits #visits/patient #conditions/patient #procedures/patient #drugs/patient MIMIC-I 35,707 44,399 1.24 12.89 4.54 33.71 MIMIC-IV 123,488 232,263 1.88 21.74 4.70 43.89 C.2 Preprocessing and Concept Normalization We applied the same preprocessing pipeline to both MIMIC-I and MIMIC-IV, including cohort selection, code normalization, and data splitting, to obtain consistent visit-based patient EHR sequences (Jiang et al., 2024). We used the same normalization procedure across datasets to ensure a fair comparison. Cohort construction and splitting. We treat each inpatient admission as a visit and represent each patient record as an ordered sequence of visits. We perform patient-level train/validation/test splits so that all visits from the same patient appear in only one split, avoiding information leakage across splits. We randomly divide the dataset into training, validation and test sets in a 60%: 20%: 20% ratio. After splitting patients, we independently construct GtrainG_train, GvalG_val, and GtestG_test using patients within each split. Diagnosis and procedure codes. For each visit, we extract diagnoses and procedures from the corresponding MIMIC tables and represent them as unordered sets of codes. We retain the original ICD coding systems (ICD-9-CM and ICD-10-CM for diagnoses; ICD-9 and ICD-10-PCS for procedures) and preserve the recorded ICD version to avoid conflating codes across systems. We map ICD-9/10 diagnosis and procedure codes to CCS categories, and de-duplicate repeated codes within each visit. Specifically, 11,736 ICD-9-CM and 72,446 ICD-10-CM diagnosis codes are mapped to 285 CCS condition categories, while 4,670 ICD-9-PROC and 79,758 ICD-10-PCS procedure codes are mapped to 231 CCS procedure categories, yielding 285 condition and 231 procedure concepts in total. Medication mapping and normalization. Medications are extracted from the prescriptions tables, where drugs are originally recorded as National Drug Codes (NDC). Following GraphCare (Jiang et al., 2024), we standardize medications using a two-stage mapping pipeline: NDC → RxNorm (RxCUI) → ATC. We keep ATC level-3 codes to preserve fine-grained medication semantics. Unmappable NDC entries (e.g., non-drug supplies, institution-specific packaging codes, or missing references) are excluded from the medication representation. This yields 269 unique drugs represented as ATC level-3 codes. Visit-level feature construction. Each visit xix_i is represented by three unordered sets: diagnosis codes, procedure codes, and mapped ATC-3 medication codes. These per-visit sets are then composed into a visit-based patient history x1,…,xt\x_1,…,x_t\ by admission time. Basic cohort statistics after preprocessing are reported in Appendix C.1. Task instance generation. Given a patient sequence x1,…,xt\x_1,…,x_t\, we generate training instances from visit prefixes (x1:i)(x_1:i) and apply supervision on the subsequent visit, consistent with the task definitions in Appendix C.3. For each task instance, to avoid future information leakage, graph edges are computed using only diagnosis codes available within the corresponding input prefix. C.3 Task Definitions and Metric Computation We follow the task formulation and evaluation protocol as prior work (Jiang et al., 2024). Given a multi-visit patient sequence x1,…,xt\x_1,…,x_t\, we treat each visit xix_i as the current visit and construct a multi-visit inputs (x1),(x1,x2),…,(x1,…,xt)\(x_1),(x_1,x_2),…,(x_1,…,x_t)\ for prediction. Unless otherwise stated, we report the average metrics and standard deviation over 10 independent runs, as shown in Table 1. Readmission prediction. For each visit xix_i with i≤t−1i≤ t-1, τ(xi)τ(x_i) denote the admission time of visit xix_i (in the de-identified timeline, preserving within-patient time intervals). For multi-visit patient sequence (x1:i)(x_1:i) with i≤t−1i≤ t-1, we define the readmission label as yireadm=[τ(xi+1)−τ(xi)≤σ],y readm_i=I [τ(x_i+1)-τ(x_i)≤σ ], (12) where we set σ=15σ=15 days in our experiments, and the label indicates whether a subsequent visit occurs within the readmission window. We evaluate readmission using AUPRC and AUROC on both MIMIC-I and MIMIC-IV datasets. LOS prediction. LOS prediction is formulated as a C-category classification problem with C=10C=10 classes, corresponding to stays of <1<1 day (0), 11-77 days (1-7), 88-1414 days (8), and ≥15≥ 15 days (9). For each visit xix_i, the goal is to predict the LOS category of the current visit xix_i. We report AUPRC and (macro-)F1 scores on both MIMIC-I and MIMIC-IV datasets. Drug recommendation. Drug recommendation is formulated as a multi-label classification task. Let drugV_ drug denote the medication vocabulary and Si⊆drugS_i _ drug denote the medication set prescribed at the current visit xix_i. The target is represented as a multi-hot vector yidrug∈0,1|drug|y drug_i∈\0,1\^|V_ drug|. For the current visit xix_i, the model input consists of all EHR codes from prior visits x1,…,xi−1\x_1,…,x_i-1\, including their diagnosis, procedure, and medication codes, together with the diagnosis and procedure codes observed at the current visit xix_i.The medication set prescribed at the current visit, denoted by SiS_i, are excluded from the input and serve as the prediction target. We report AUPRC, F1, and Jaccard on both datasets. For F1 and Jaccard, we convert predicted probabilities to a binary set using a fixed threshold of 0.20.2. Jaccard is computed per instance as |Si∩S^i||Si∪S^i| |S_i∩ S_i||S_i∪ S_i| and averaged over the test set. Training objectives. We use binary cross-entropy (BCE) with sigmoid for readmission and drug recommendation, and cross-entropy (CE) with softmax for the multi-class LOS task. C.4 Baselines and Reproducibility Notes We describe the baseline implementations and the evaluation protocol used to ensure fair comparison. Baselines. We compare against representative patient-level sequence models, namely Deepr (Nguyen et al., 2017), RETAIN (Choi et al., 2016b), GRAM (Choi et al., 2017), StageNet (Gao et al., 2020),AdaCare (Ma et al., 2020), and GRASP (Zhang et al., 2021). We further compare with LM-graph models that incorporate external medical knowledge or graph, including G-BERT (Shang et al., 2019a), LEADER (Liu et al., 2024), GLEM (Zhao et al., 2023), GraphCare (Jiang et al., 2024), KARE (Jiang et al., 2025), ColaCare(Wang et al., 2025). For drug recommendation, we also include task-specific models SafeDrug (Yang et al., 2021b), MICRON (Yang et al., 2021a), GAMENet (Shang et al., 2019b), MoleRec (Yang et al., 2023), and UDC (Zhao et al., 2025). Implementation sources and training protocol. Where official implementations are available, we follow the authors’ released code and default settings. Otherwise, we re-implement baselines based on the original papers. All methods are trained and evaluated under the same preprocessing pipeline, patient-level data split, and task instance construction (Appendix C.2-C.3). Hyperparameters for each baseline are tuned on the validation set and the selected configuration is used for final testing. Unless stated otherwise, we use early stopping based on the validation metric corresponding to each task and report results over multiple runs (Table 1). Input alignment across methods. For sequence-based baselines such as Deepr, RETAIN, GRAM, and StageNet, each visit is represented with the same set-valued clinical concepts used by our method, and visits are ordered chronologically. LM/graph-enhanced baselines such as G-BERT, LEADER, GLEM, and GraphCare are evaluated using the same normalized concept vocabularies and the same patient-level split to ensure comparability. Applicability constraints. AdaCare and GRASP are designed for binary clinical classification and do not directly apply to the multi-class LOS or multi-label drug recommendation settings under our formulation; therefore, we report their results only on the readmission task, consistent with Table 1. Drug-specific knowledge. For drug recommendation baselines that rely on drug-drug interaction (DDI) knowledge (e.g., SafeDrug and GAMENet), we follow the original papers and corresponding implementations to construct the required drug graphs/knowledge, while keeping the same ATC-level medication normalization used in our preprocessing (Appendix C.2). All methods are evaluated on the same target medication space. Reproducibility. We fix random seeds for each run and report mean and standard deviation across 10 runs. Complete hyperparameters, training budgets, and hardware details are provided in Appendix C.5. C.5 Implementation Details and Hyperparameters This section summarizes model configurations, optimization settings, and reproducibility details. Backbone language models. We instantiate our framework with two biomedical language models: (i) an encoder-only BioBERT-large model (biobert-large-cased-v1.1) and (i) a decoder-only BioMistral model. For BioBERT-large, we fine-tune the last six Transformer layers together with a task-specific prediction head, while freezing all earlier layers. For BioMistral, we extract sequence-level representations using LLM2Vec (BehnamGhader et al., 2024) and fine-tune the model via parameter-efficient LoRA adapters, freezing all non-adapted parameters. LLM2Vec configuration. We follow the standard LLM2Vec setup to obtain fixed-length sequence embeddings from the decoder-only backbone. Specifically, we use the last-layer hidden states and apply mean pooling over the token dimension, followed by a linear projection to match the hidden dimension of the downstream graph model. This configuration is shared across all tasks and datasets. LoRA configuration (BioMistral). LoRA adapters are applied to the attention projection layers of BioMistral. Unless otherwise stated, we use a uniform configuration across tasks: rank r=8r=8, scaling factor α=16α=16, and LoRA dropout 0.10.1. Adapters are inserted into the q_proj,k_proj,v_proj,o_proj\q\_proj,k\_proj,v\_proj,o\_proj\ modules of each attention block. GNN. We use a standard three-layer GCN with a hidden dimension of 128 and ReLU activation, followed by a softmax classifier (Kipf and Welling, 2017b). The GCN layer updates node representations based on a weighted patient graph with adjacency A, where Aij≥0A_ij≥ 0 encodes patient similarity: H(ℓ+1)=σ(D~−12A~D~−12H(ℓ)W(ℓ)),A~=A+I, splitH^( +1)&=σ\! ( D^- 12 A D^- 12H^( )W^( ) ),\\ A&=A+I, split (13) where D~ D is the degree matrix of A~ A. Equivalently, for node i (ignoring σ), hi(ℓ+1)=∑j∈(i)∪iA~ijd~id~jhj(ℓ)W(ℓ).h_i^( +1)= _j (i)∪\i\ A_ij d_i\, d_j\;h_j^( )W^( ). (14) This formulation naturally supports weighted patient graphs: the edge weight AijA_ij enters directly through A~ij A_ij and modulates the strength of message passing between patients. Since the LM produces higher-dimensional embeddings (e.g., 768 for BioBERT and 4096 for BioMistral), we use an MLP projection to map them to 128 as the GCN input. Optimization and training. All models are optimized using AdamW. We use separate learning rates for the language model and the graph/prediction components: ηlm=1×10−5 _ lm=1× 10^-5 and ηgnn=1×10−3 _ gnn=1× 10^-3. Weight decay is set to 1×10−21× 10^-2. We train with batch size 3232 and use gradient accumulation of 22 steps when necessary. The maximum input sequence length is capped at 512512 tokens. We adopt a linear learning-rate scheduler with a warmup ratio of 10%10\% of the total training steps. Early stopping is applied based on the primary validation metric of each task with a patience of 55 epochs. Training budget and hardware. All experiments are conducted on NVIDIA Dgx Spark GPU with 128GB of memory. We train models for up to 1010 epochs with early stopping. Mixed-precision training (FP16) is used to improve training efficiency. Random seeds and reporting. We use 10 random seeds for 10 runs, with one seed assigned to each run. The same seeds are used across all methods and baselines, and we report the mean and standard deviation over the 10 runs. C.6 Details of the Weighted Patient Graph We provide additional details on the construction of the weighted patient graph. Following prior work (Lu and Uddin, 2021), we first represent the cohort as a bipartite graph between patients and diagnosis codes. A bipartite graph is a special class of graph consisting of two disjoint sets of vertices. In this study, we use an undirected bipartite graph to represent relationships between patients and diseases. Here, we project the bipartite graph into the patient side, to construct a weighted patient network: two patients are connected if they share at least one diagnosed disease, and the edge weight equals the number of diseases they have in common (i.e., their number of shared neighbors). To avoid counting repeated occurrences, we compute edge weights using unique diagnosis codes per patient. As illustrated in Fig. 1b, patients 11 and 22 have one EHR code, resulting in a weight of 1 between them; patients 33 and 44 have a weight of 2 due to two shared EHR codes. In the projected graph, each patient node retains its original properties, and two patients are connected if they share at least one neighbor (i.e., a diagnosis code) in the original bipartite graph. Prior studies suggest that diseases that co-occur can reflect shared biological mechanisms (e.g., via disease–gene interactions) and that patients with the same chronic conditions often share common risk factors (e.g., smoking history, obesity, and insufficient physical activity) (Lu and Uddin, 2021). Motivated by these observations, we use the weighted patient network to capture latent patient–patient relationships induced by shared clinical profiles. We deliberately construct the patient graph from diagnosis overlap only. This design is leakage-safe across our tasks: incorporating current-visit medications would directly leak targets for drug recommendation, while incorporating procedures or other current-visit signals could partially reveal the LOS target. We therefore use diagnosis overlap as a simple, task-agnostic cohort prior that remains compatible with all three prediction settings. Richer graph signals, such as temporality, disease severity, or demographics, may further improve graph quality, but they are not straightforward to use with the MIMIC datasets. Figure 4: Edge-weight distributions in the weighted patient graphs for the training cohort (6,000 patients) and test cohort (2,000 patients). Edge weights are defined as the number of shared diagnosis codes between two patients and are binned as 1-10 and 10+. The two cohorts exhibit highly similar distributions with Jensen–Shannon Divergence of 0.0097 and Total Variation Distance of 0.0052, suggesting minimal practical distribution shift. Fig. 4 shows the edge-weight distributions of the weighted patient graphs constructed from 6,000 training patients and 2,000 test patients sampled from MIMIC-I. The histogram is heavily skewed toward small weights (primarily 1–3 shared diagnosis codes), and the frequency drops rapidly as the weight increases. This long-tailed pattern indicates a sparse graph in which most patient pairs are only weakly connected, with a small minority exhibiting strong similarity through many shared diagnoses.We compare the edge-weight distributions of the training and test patient graphs using a G-test on binned weights. The distributions are highly similar, with Jensen–Shannon Divergence of 0.0097 and Total Variation Distance of 0.0052, indicating minimal distribution shift. This indicates similar graph sparsity and edge-weight patterns across the training and test cohorts. To keep the resulting patient graph sparse, we further sparsify the projected graph by retaining an edge only when its weight exceeds a minimum threshold, i.e., Aij≥τA_ij≥τ. Based on the distribution in Fig. 4, we set τ=8τ=8 to filter out the large number of weak connections. This choice yields a substantially sparser patient graph, with an average node degree on the order of ∼ 8 in the sampled cohort, while preserving clinically meaningful links among genuinely similar patients. The training, validation, and test graphs are constructed independently. The model is trained only on GtrainG_train. During testing, all test patients are encoded by the fixed LM and jointly processed by the fixed GNN on GtestG_test. Test–test edges therefore enable representation-level message passing, while test labels are never used during inference. C.7 Additional Details on Interpretability Analysis Figure 5: Comparison of two explanation signals for the same target patient (57139). (a) Patients-like-me ranks reference patients and assigns a patient-level importance score (Imp) indicating each reference patient’s contribution to the target prediction. (b) Baseline visualization using the weighted patient graph, where the edge weight AijA_ij equals the number of shared diagnosis codes between the target patient and each reference patient. Figure 5 contrasts two explanations for the same target patient 57139 in the insulin drug recommendation task. We visualize an ego-network comprising the target node and its 10 reference patients retrieved by Patients-like-me. Each reference node is labeled with its importance score (Imp). The higher Imp is encoded by a shorter edge length (closer to the target) and a darker blue node color, enabling qualitative comparison across references. In panel (a), the target patient is placed at the center and the 10 reference patients are arranged around it for clarity. Each reference node is annotated with an importance score that quantifies its influence on the target prediction. To improve interpretability, we additionally annotate a the patient-level clinical evidence using CCS codes: (i) direct diabetes CCS evidence that is strongly associated with insulin use (e.g., 50: Diabetes mellitus with complications), and (i) indirectly related CCS categories capturing endocrine conditions or diabetes-related complications (e.g., chronic kidney disease). These annotations are extracted from the mapped CCS codes in each patient’s structured record and shown near each reference node. Panel (b) provides a structural baseline explanation derived from the weighted patient graph used by our model. We use the same target patient and the same set of 10 reference patients, but replace the Patients-like-me importance score with the edge weight AijA_ij between patients i and j. We visualize the same ego graph, encoding the edge weight AijA_ij by edge length (larger AijA_ij yields a shorter edge, indicating stronger structural proximity), and using node color to reflect the same weight for consistency. This structural baseline reflects a common assumption in cohort graphs: patients are considered more related if they share more structured codes. In panel (a), Patients-like-me assigns the highest importance to patient 305 with an importance score 0.41, which is clinically aligned with the presence of direct diabetes CCS code 50. Moreover, patients 107 and 384 also receive high importance scores (0.33 and 0.29), consistent with multiple indirectly related endocrine conditions and diabetes complications annotated near those nodes. Other reference patients have fewer indirectly related CCS codes and therefore receive lower importance scores. In contrast, the structural baseline in panel (b) ranks references purely by code overlap, , which can yield less informative explanations. For example, patient 305, although the most influential reference in panel (a), has only an edge weight of 9, only slightly above the sparsification threshold. This indicates that edge weights alone only capture structural shared-code overlap and can under-emphasize clinically salient evidence that drives the model’s prediction. Overall, Fig. 5 shows that Patients-like-me provides explanations by retrieving reference patients that are not only structurally similar in the cohort graph but also clinically and semantically relevant to the target prediction. This comparison highlights that our patient-level importance scores are more informative than explanations based only on graph connectivity. C.8 Monte Carlo Sampling Analysis We next examine sensitivity to the number of Monte Carlo samples used in the E-step. The expectation in Eq. 4 is taken over latent labels y, not over high-dimensional continuous LM embeddings; the LM embedding ono_n and GNN representation hn(L)h_n^(L) are deterministic forward-pass outputs. Table 5 shows that increasing the estimator from one to five samples on MIMIC-IV readmission yields only modest gains, supporting our default single-sample setting. Table 5: Sensitivity to the number of Monte Carlo samples in the E-step. Results are shown on the MIMIC-IV readmission task with the BioBERT backbone. Number of MC samples AUPRC (%) AUROC (%) Single sample (default) 47.8 (± 0.4) 80.1 (± 0.4) Five samples 49.0 (± 0.2) 80.8 (± 0.3) C.9 Backbone Generalization To assess whether the gains from VEM extend across different LM backbones, we report MIMIC-IV results on readmission, LOS, and drug recommendation for four backbones spanning both encoder-only models (BCMBERT and BioBERT) and decoder-only models (Meerkat and BioMistral). Table 6 shows that VEM consistently improves performance over the corresponding LM-only baseline across all tasks and all reported metrics. The gains are especially pronounced for BCMBERT and Meerkat, where VEM improves readmission AUPRC by 5.0–5.8 points, LOS AUROC by 5.8–6.2 points, and drug recommendation Jaccard by 3.4–4.6 points. Even for the stronger BioBERT and BioMistral baselines, VEM remains consistently beneficial; for example, it improves BioBERT readmission AUPRC from 45.4% to 47.8% and BioMistral drug recommendation F1 from 58.3% to 62.4%. Overall, these results indicate that Patients-like-me is not tied to a specific LM architecture and generalizes well across both encoder-only and decoder-only biomedical backbones. Table 6: Generalization of PLM across LM backbones on MIMIC-IV. We evaluate PLM with the VEM algorithm across four backbones on readmission, LOS, and drug recommendation. VEM consistently improves performance over the corresponding LM-only baseline, showing that PLM generalizes across backbones. Backbone Training Readmission LOS Drug Recommendation AUPRC AUROC AUROC F1 AUPRC F1 Jaccard BCMBERT-396M VEM 48.6 80.6 82.8 35.0 75.4 64.8 48.9 LM-only 43.6 74.5 77.0 30.7 70.2 61.4 45.5 Meerkat-8B VEM 49.7 80.3 85.3 35.6 77.1 65.7 50.4 LM-only 43.9 75.4 79.1 30.9 71.3 60.7 45.8 BioBERT-300M VEM 47.8 80.1 80.5 33.8 75.1 64.2 48.5 LM-only 45.4 77.4 78.1 33.7 73.1 61.7 46.8 BioMistral-7B VEM 48.3 78.5 84.8 34.5 74.7 62.4 47.5 LM-only 42.9 74.7 82.1 31.3 71.9 58.3 45.2 C.10 Efficiency Analysis Table 7: Efficiency on MIMIC-IV with Meerkat-8B. Time is hours per epoch, memory is peak GPU memory (GB), and metrics are AUPRC/F1/Jaccard (%). Training Time (h) Memory (GB) AUPRC F1 Jaccard VEM 2.55 80.4 77.1 65.7 50.4 Alternating 2.52 80.8 72.4 60.9 47.1 E2E 2.61 91.1 70.4 60.2 45.4 2-stage 2.63 80.8 73.2 63.6 47.9 LM-only 2.45 76.8 71.3 60.7 45.8 Table 7 reports runtime and peak GPU memory on MIMIC-IV with the Meerkat-8B backbone. Relative to LM-only, VEM increases training time from 2.45 to 2.55 hours per epoch and peak memory from 76.8 to 80.4 GB, i.e., only a 4.1% runtime increase and a 4.7% memory increase. It also has lower runtime than E2E and 2-stage training, while maintaining comparable memory usage, showing that VEM delivers consistent performance gains with only modest computational overhead. Compared with alternating optimization, VEM requires nearly identical training time (2.55 vs. 2.52 hours per epoch) and slightly less GPU memory (80.4 vs. 80.8 GB), while improving AUPRC, F1, and Jaccard by 4.7, 4.8, and 3.3 %, respectively.