Bayesian machine learning enables discovery of risk factors for hepatosplenic multimorbidity related to schistosomiasis

Domain: chemistry
Published: 2026-03-03
Editorial outcome: major_revision

Abstract
One in 25 deaths worldwide is related to liver disease, and often with multiple hepatosplenic conditions. Yet, little is understood of the risk factors for hepatosplenic multimorbidity, especially in the context of chronic infections. We present a novel Bayesian multitask learning framework to jointly model 45 hepatosplenic conditions assessed using point-of-care B-mode ultrasound for 3155 individuals aged 5-91 years within the SchistoTrack cohort across rural Uganda where chronic intestinal schistosomiasis is endemic. We identify distinct and shared biomedical, socioeconomic, and spatial risk factors for individual conditions and hepatosplenic multimorbidity, and introduce methods for measuring condition dependencies as risk factors. Notably, for gastro-oesophageal varices, we discover key risk factors of older age, lower hemoglobin concentration, and schistosomal periportal fibrosis. Our findings provide a compendium of risk factors to inform surveillance, triage, and follow-up, while our model enables improved prediction of hepatosplenic multimorbidity, and if validated on other anatomical systems, general multimorbidity.

Introduction
Multimorbidity is the co-occurrence of two or more chronic conditions and is of growing global concern in existing health systems due to the lack of clinical guidelines, strains on resources and staff, and potential mismanagement or misdiagnosis of patients [1–3]. Individuals with multimorbidity are among the most socially and economically disadvantaged with greater years lived with disability and higher mortality compared to those without multimorbidity [2]. In low and middle-income countries (LMICs), over one third of individuals are affected by complex multimorbidity due to both infectious and non-communicable causes [4, 5]. Yet, there is little understanding of multimorbidity due to chronic infections in LMICs with most multimorbidity studies focusing on aging or cardiovascular diseases in high-income countries [6–8].

Schistosomiasis, caused by the intestinal species Schistosoma mansoni, is a major cause of hepatosplenic multimorbidity in sub-Saharan Africa. The most common hepatosplenic conditions that occur in schistosomiasis-endemic populations include periportal fibrosis, early cirrhosis-like or hepatitis-like livers, splenomegaly, moderately shrunken or enlarged livers, and thickened gall bladder walls–all of which are often accompanied by clinical symptoms such as diarrhea, abdominal pain, anaemia, jaundice, and melena [9–11]. Hepatosplenic multimorbidity due to chronic schistosome infections encompasses two or more interacting or co-occurring conditions ranging from combinations of less severe forms of periportal fibrosis and changes in spleen or liver organometry to inclusion of life-threatening complications of gastro-oesophageal varices and liver cirrhosis [10]. More than one in every two individuals aged five years and older living in areas endemic with S. mansoni have been estimated to have hepatosplenic multimorbidity [10]. Yet, there is a limited understanding of the risk factors of hepatosplenic multimorbidity or, more generally, diverse hepatosplenic outcomes beyond periportal fibrosis [9]. Currently, only simple composite outcomes are studied where diverse conditions are collapsed into one indicator (periportal fibrosis [9]) or less severe conditions (diarrhoea, etc. [12]) are assessed as independent outcomes. From these singular outcome studies, current schistosome infection has been shown to be irrelevant for periportal fibrosis prediction [9], despite still being used as proxy indicator for morbidity by the World Health Organization (WHO) [13–15]. In contrast, co-infection with human immunodeficiency virus (HIV) and hepatitis B (HBV) infections have been shown to confer over two times higher odds of periportal fibrosis coupled with the past history of schistosome exposure [9]. The complex interplay of pathogenesis provides support for a renewed focus on multimorbidity and interacting conditions. It is unknown whether co-infections result in biologically-related conditions that interact to exacerbate schistosomiasis-specific conditions or cause independent co-occurring conditions due to shared socioeconomic and spatial risk factors. Life-threatening conditions related to schistosomiasis, such as gastro-oesophageal varices that are indicative of portal hypertension, remain still to be studied outside of health facilities and beyond basic demographic risk factors [11] and little is known about how multimorbidity influences schistosomiasis-related condition progression.

For schistosomiasis, the common approach of focusing narrowly on one condition at a time fails to account for all conditions currently relevant to a patient. Consequently, there is a limited understanding of risk factors shared across conditions and, critically, condition inter-dependencies, to be able to identify early indicators of severe conditions and disease progression. This oversimplified approach to multimorbidity modelling is not unique to schistosomiasis and applies more generally to recent attempts to study multimorbidity. Many studies simply identify multimorbidity as a binary indicator of two or more conditions from an arbitrary list of predefined conditions regardless of the actual condition type or interaction [2, 5, 16, 17].

Given multimorbidity arises from shared risk factors including condition inter-dependencies, we developed a Bayesian multitask learning architecture [18, 19], allowing each condition to learn from the risk factors of other conditions [20]. Graphs (or networks) were utilized to encode inter-dependencies between conditions, and parameterized to homogenize the predictions based on node connectivity. We jointly predicted 45 hepatosplenic conditions assessed through point-of-care B-mode ultrasound (with the WHO Niamey Protocol for defining different periportal fibrosis patterns [21]) for 3155 individuals aged 5-91 years within the community-based cohort, SchistoTrack, in rural Uganda. The aim of this study was to employ Bayesian machine learning and multitask modelling to discover shared risk factors from a wide range of biomedical, sociodemographic, and spatial variables for hepatosplenic multimorbidity, in particular life-threatening severe hepatosplenic conditions, and to measure the influence of condition dependencies.

Results
Risk factors for hepatosplenic conditions

An illustration of how the results were obtained can be found in Fig. 1. In our study population, 82% (2578/3155) of individuals had at least one condition, and 54% (1708/3155) had two or more conditions. When considering the maximum count of multimorbidity across all study participants including those without multimorbidity, the top 5.9% (185/3155) of all study participants had six or more conditions (maximum number of conditions in any one individual was 13); the top 5.9% across all study participants corresponded to the top 10.8% of multimorbid participants (185/1708). Mild liver fibrosis (Niamey protocol patterns C1 and C2), liver shrunkenness, and spleen enlargement were the most commonly observed conditions (see Table 1 for prevalence of all conditions). By examining a wide set of participant covariates (Table 2), we first identified significant risk factors for each of the individual 45 conditions (Table 3).

Risk factors significant for at least one condition included three clinical measurements, seven sociodemographic variables, and three spatial factors. Hemoglobin (Hb) concentration was the most informative measurement and was negatively related to 24.44% (11/45) of conditions (natural log-transformed average median odds ratios (OR) 0.20, range 0.03 - 0.54); among associated conditions, four were severe, four moderate, and three mild with respect to portal hypertension risk. Notably, among these conditions, significant variation in WHO categories of anaemia status were observed for chronic hepatitis or early cirrhosis, gastro-oesophageal varices, splenorenal shunts, ascites, and severely enlarged spleens (Supplementary Table S1). Malaria infection and S. mansoni eggs per gram (EPG) were associated with few non-severe conditions (three and two, respectively). Age was relevant for most conditions (64.44%, 29/45) irrespective of severity (10 severe, four moderate, 12 mild) with average median OR of 1.03 (range 0.97 - 1.06), and almost all significant conditions were positively related to older age (89.66%, 26/29). Both gender (female) and being a fisherman were related to six conditions (13.33%). All associations with being female were negative (average median OR 0.64, range 0.53 - 0.72) including two severe conditions, and all associations with being a fisherman were positive (average median OR 1.91, range 1.51 - 2.13) including three severe conditions. Access to an improved drinking water source had two negative (median OR range 0.74 - 0.75) and one positive (median OR 1.57) association with only non-severe conditions. The number of individuals in the household (positive), belonging to the majority tribe (negative), and belonging to a home that was owned (i.e. not rented) (positive) were related to one non-severe condition each (OR range 0.66 - 1.53). Spatial risk factors included living in the Western region of Uganda and distance to open freshwater sites. Pakwach district was associated to 31.11% (14/45, 12/14 positive) of conditions, and Buliisa district to 17.78% (8/45, 6/8 positive) when compared to Mayuge district, of which only three and one condition, respectively, were severe. The minimum distance (km) to an open freshwater site was positively associated with only one mild condition. The OR for every condition is presented in Supplementary Fig. S1, presented by the mean to show the directional effect.

Gastro-oesophageal varices were the primary focus of our analysis because they indicate severe portal hypertension. All covariate significance values for the condition are shown in Supplementary Table S2; the corresponding posterior ORs are visualized in Supplementary Fig. S2. Age and Hb concentration were the only significant variables. Each additional year of age was associated with a 4% increase in the odds of gastro-oesophageal varices (median OR 1.04, 95% credible interval (CI) 1.02 - 1.08). A 10% increase in the log-transformed Hb concentration corresponded to a 22% decrease in the likelihood of gastro-oesophageal varices (OR 0.07, 95% CI 0.01 - 0.37). From the covariate densities, a number of clear moderate to large non-zero effects of other risk factors were observed, but were borderline insignificant from uncorrected p-values. HIV (OR 2.52, CI 0.97 - 10.77), HBV (OR 1.40, CI 0.81 - 3.72), female (OR 0.67, CI 0.26 - 1.04), and fisherman (OR 1.79, CI 0.91 - 4.44) had credible intervals narrowly including 1. The non-zero effect was visible through the heavy tails of skewed posterior densities, while insignificance was due to the modal OR remaining close to one. In sensitivity analyses, we investigated the models with alternative priors as discussed in Methods. The posterior ORs for gastro-oesophageal varices from the alternative priors is shown in Fig. S3. The inference was robust to prior choice as long as they remained theoretically appropriate for the model components, as the OR distributions were extremely similar to Fig. S2.

Group risk factors were computed to profile periportal fibrosis associated to schistosomiasis through combined likelihood ratios of the five liver fibrosis patterns (Niamey protocol patterns C - F as defined by the features adjacent to portal vasculature) (Supplementary Table S2). Direction of association of each covariate was determined by the mean from combining the five posterior ORs densities. Periportal fibrosis shared many risk factors with the five patterns individually, with the addition of HIV despite being insignificant alone for singular liver fibrosis patterns. HIV exhibited a positive mean effect, with infection leading to a 1.77 higher odds of developing periportal fibrosis (Niamey protocol patterns C-F). Other significant associations were found with older age (mean OR 1.04), being female (mean OR 0.62), being a fisherman (mean OR 2.00), the number of individuals in the household (mean OR 1.09), access to improved drinking water source (mean OR 0.95), and the Western districts (Buliisa mean OR 1.31, Pakwach mean OR 1.99).

Significance against overall hepatosplenic multimorbidity also was computed on the 45 conditions (Supplementary Table S2). Predicting overall multimorbidity represented the expected morbidity counts in each participant, scenarios that lead to the highest overall multimorbidity would represent the expected maximum morbidity count. Six risk factors previously observed to be significant for individual conditions remained. The factors leading to higher expected morbidity counts included lower Hb concentration, older age, being male, working as a fisherman, and living in the Western districts (with the strongest effect in Pakwach district). To identify conditions with the highest probabilities for overall hepatosplenic multimorbidity, an extreme scenario was simulated using the top 99% percentile or the bottom 1% percentile of the significant covariate values (depending on the direction of association) while keeping insignificant covariates at median value (Fig. S23). The expected number of conditions the patient will have in this extreme scenario is 7.86, while the most likely conditions to develop were severe enlargement of the spleen (probability 0.69), and any form of enlargement of the portal vein (probability 0.69 for moderate enlargement and 0.64 for severe enlargement). Findings for individual conditions, periportal fibrosis, and overall hepatosplenic multimorbidity remained robust when stratified by children (aged 5-17 years) and adults (aged 18+ years), notably revealing most risk factors associations were found in adults (Supplementary Tables S3 - S6).

Graph convolution and condition inter-dependencies

Condition dependencies for gastro-oesophageal varices using the graph are shown in Fig. 2 where node weightings (condition dependency measure) for the condition were controlled by parameter α that was chosen as the median value from its posterior MCMC samples (Supplementary Fig. S4). A set of probabilities between all conditions can be found in Supplementary Fig. S5. The self-influence for gastro-oesophageal varices indicated a probability of 0.397, measuring how much was predicted from its own model and risk factors/covariates. From the other condition models, the highest probability contributions to gastro-oesophageal varices came from related biological complications of splenic varices and splenorenal shunting, conferring additional probabilities of 0.102 and 0.081, respectively. The next most relevant conditions encompassed potential fibrotic precursors to portal hypertension and included moderate and severe forms of schistosomal periportal fibrosis. The ruff liver pattern (Niamey protocol pattern D) where there was extensive fibrosis at the point of portal bifurcation increased the probability by 0.045. Liver patches (Niamey protocol pattern E), which were indicative of occluded vessels and fibrosis blocking the main portal vein, conferred an increased probability of 0.059. Other condition dependencies included a severely enlarged spleen (0.049 probability), and severely enlarged (mean) portal vein diameter (0.046 probability). All other conditions had less than 0.030 probability influence, with the next most relevant being bird’s claw, another severe schistosomal periportal fibrosis pattern that represents a blocked main portal vein, occluded vessels, and extensive fibrosis extending to the liver capsule (Niamey protocol pattern F) with 0.027 probability.

The graph convolution probability between splenic and gastro-oesophageal varices was among the highest donated probabilities between conditions (Table S9), and hence among the strongest connections between any two conditions. The robustness of the probabilities around gastro-oesophageal varices under perturbation are visualized in Fig. S20, S21, and S22 for 5%, 10%, and 20% perturbations to edges, respectively. When both addition and removal of edges were applied, the weights of reported dependencies of gastro-oesophageal varices consistently remained within the inter-quartile range of the perturbed convolution weights for all percentages. When edges were added only, the weights, being normalized, generally became reduced in proportion, but the upper-quartiles continued to contain the reported dependency weights.

For periportal fibrosis (Niamey protocol patterns C to F), all patterns had self-influence probabilities greater than 0.4 (Supplementary Fig. S14 - S18). The patterns were highly inter-dependent and featured in their dependencies. The D grade, being the intermediate grade between assumed fibrosis progression from C1/C2 to E was in the top five for all three dependencies with influence probabilities from 0.023 - 0.038 and 0.030, respectively. C1 and C2 are two different cross-sections of the liver that capture similar levels of vessel-related fibrosis, and as such exhibited strong dependence (0.275 - 0.278 probabilities). The thickening and presentation of many second order branches (spider thickening fibrosis, Niamey protocol pattern B2, a non-periportal fibrosis associated liver pattern), also appeared in the top five for the two C patterns (0.018 - 0.032 probabilities). The most severe grades of schistosomal periportal fibrosis, which were anatomically nested (Niamey protocol patterns E and F) also were strongly connected (> 0.175 probabilities), and strikingly both E and F were highly dependent on a gross undulating liver surface (0.041 - 0.069 probabilities) typically indicative of liver cirrhosis.

To give an alternative view on condition groupings using the graph convolution network, we show the clustering of conditions found through community detection on the graph (Fig. S24). Three main clusters were found. One cluster mainly contained measurements of liver enlargement, fatty livers, and other abnormalities not specifically attributable to infections and potentially due to alcohol use. Another cluster that included cirrhosis conditions that represent features, complications, of consequences of cirrhosis that may be attributable to alcohol use, HIV, or HBV such as rounded liver edges, shrunken liver lobes, and gross liver surface irregularities. The third cluster encompassed mostly schistosomal periportal fibrosis Niamey patterns, its features such as changes in the main portal vein, and portal hypertension consequences including gastro-oesophageal varices and splenorenal shunts.

Multimorbidity prediction

The predictive performances based on 50:50 training-testing split of the data were compared against widely used machine learning models featuring various frameworks for multiple outcomes. Models used for comparison included logistic regression, XGBoost, multi-output random forest, and shallow neural networks with single tasks and multitask architectures. These results can be found in Table 4, with some additional metrics also presented in Supplementary Table S7. Our Bayesian multitask model using both the graph and covariate learning (BMT-G&CL), which was the focus of all aforementioned results, was superior even against models with far more complex structures, producing the highest AUC, AUC-precision-recall (AUC-PR), F1 score, precision, and positive predicted values (PPV). The reduced version consisting of only the covariate learning multitask element (BMT-CL) performed nearly as well on all metrics and also was among the top performing models. There was a notable difference between models that were multitask and single task (i.e. treating outcomes separately). Our Bayesian multitask models (full model and without either the graph or covariate learning component) along with the multitask neural network (MT-NN) had markedly higher AUC, F1, and PPV, while also performing marginally better with AUC-PR and precision due to the low prevalence leading to low values across all models.

We show the predictions from the BMT-G&CL for individual conditions using the AUC-PR in Supplementary Fig. S6, with the prevalence of each condition marked for reference. The AUC-PR alone indicated that the best predicted conditions were generally of higher prevalence, but we also assessed AUC-PR accounting for the prevalence, which is equivalent to comparing against random classification (which will have expected AUC-PR equal to the prevalence). This comparison was represented by the AUC-PR-to-prevalence ratio, and can be found in Fig. 3. The ratios indicated that our model generally had greater proportional improvement on rarer and often severe conditions compared to random classification. In particular, the model was over 20 times better at predicting the existence of ascites, pancreaticoduodenal varices, and recanalized paraumbilical vein, and most important, over 15 times better at predicting gastro-oesophageal varices and non-specified portosystemic collaterals when compared to random classification.

When comparing predictions on individual conditions, the BMT-G&CL did not show any consistent differences to the MT-NN (best performing alternative model) on any particular metric (AAUC Fig. S8, AUC-PR Fig, S10, and PPV Fig. S12). Against logistic regression (the standard epidemiological approach), BMT-G&CL had superior AUC (Fig. S7) in a number of conditions, but few significant differences were observed in AUC-PR (Fig. S9) and PPV (Fig. S11). AUC and AUC-PR curves can be found in Fig. S13.

To understand district generalisability, a predictive experiment holding out one district at a time is presented in Supplementary Table S8. District-based prediction depended on the metric and the method used to select probability thresholds for converting to binary, but Mayuge district was the best predicted in terms of AUC-PR despite having the lowest overall prevalence in its population.

Methods
Ethics approvals

Our study complies with all relevant ethical regulations. Data collection and use were reviewed and approved by Oxford Tropical Research Ethics Committee (OxTREC 509-21), Vector Control Division Research Ethics Committee of the Uganda Ministry of Health (VCDREC146), and Uganda National Council for Science and Technology (UNCST HS 1664ES). Written informed consent was obtained for adult participants, with adults consenting on behalf of the children after receiving their informed assent. Ethics approval was obtained from all local ethics committees. Participants were compensated in-kind for their time with a large (1kg) bar of soap.

Participant sampling

This cross-sectional study was conducted as a nested design within the SchistoTrack cohort [41] during the first annual follow-up between 17 January and 16 February 2023. Participants were selected from 1952 randomly sampled households from 52 villages across Buliisa, Pakwach, and Mayuge Districts of Uganda; 38 of the villages were sampled in the baseline of 2022 [9]. One child aged 5-17 years and one adult aged 18-91 were selected by the household head or spouse and invited for clinical assessments. 3224 individuals were examined with point-of-care ultrasound in 2023. There were 3186 individuals out of 32 24 with non-missing ultrasound data, and were featured in learning the graph [10] utilized for this study. All participants had household demographic covariate data, but 3155 had non-missing data from blood and stool assessments and were considered for analysis in this work. A breakdown of missing data can be found in the participants flow diagram in Supplementary Fig. S25.

Hepatosplenic outcomes

The outcomes were hepatosplenic conditions obtained by point-of-care B-mode ultrasound in abdominal view. Philips Lumify C5-2 curved linear array transducers were used with the Philips Lumify Ultrasound Application v3.0 on Lenovo 8505-F tablets with Android 9 Pie. Lossless DICOM images and videos were saved for quality assurance [9]. We assessed a total of 45 hepatosplenic conditions listed in Table 1. Sonographers were blind to participant covariates except for age and gender. The conditions were measured as binary indicators [10], including focal and diffuse liver fibrosis patterns following the WHO Niamey protocol [21], liver surface irregularities, caudal liver edge assessments, fatty and cirrhotic livers, liver and spleen organometry, portal vein dilation or restriction, portosystemic collaterals, ascites, gall bladder obstruction, among others. For the left and right liver lobes, spleen, and portal vein diameter, we considered abnormality as one or two standard deviations above or below an internal healthy reference population from the same year split by height, as described in [10]. Importantly, all conditions were assigned a severity rating in relation to schistosomal portal hypertension as shown in the final column, these were determined from expert opinion by a clinical epidemiologist (GFC), sonographer (VA), and gastroenterologist (CKO) as detailed in [10].

Participant covariates

Participants covariates were collected through a combination of household surveys, and diagnoses by laboratory technicians and nurses. Information was recorded using Open Data Kit version. 2022.4. These are listed in Table 2, detailed definitions can be found in [9, 22]. Clinical measurements included S. mansoni status as measured by microscopy and eggs per gram, malaria rapid diagnostic test outcomes, HBV rapid diagnostic test outcomes, self-reported HIV, and point-of-care Hb concentration. Although not a model covariate, anaemia status was investigated by constructing indicators based on WHO definitions where thresholds for mild, moderate, and severe anaemia considered age, gender, and pregnancy [42]. Socio-demographic information included age, gender, belonging to the majority tribe of their district, belonging to the majority religion of their district, years of educational attainment, occupations (farmer, fisherman, and fishmonger, with unemployed and other occupations as reference). Home quality score was calculated based on the quality of materials used to build the house. Individuals belonged to a household with social status if any adult member held or previously held a position in the local council, village health team, religious or clan leadership, or influential beach management committees. Additional household variables included number of individuals in the household, number of years household lived in the village, if they owned the home they lived in, number of rooms in the household, current alcohol use based on self-reported consumption within the year preceding recruitment (based on the WHO STEPwise surveys [43]), improved water source defined by if they obtained water from any of protected well or spring, borehole, village tap, and rainwater, number of water activities, and the year of recruitment (2022 or 2023, with the baseline recruitment year of 2022 as reference). Notably, for alcohol use, we also asked about past drinking following the WHO STEPwise approach, but did not consider this variable further as only 1.43% (45/3155) of all study participants had previously drank and stopped drinking. We did not focus on the frequency or amount of alcohol consumed as most individuals in our study were not regular drinkers [44]. Finally, locational variables measured the distance (in kilometres) to nearest freshwater site, and government health centre, and the respective districts (Buliisa, Pakwach, with Mayuge as reference). The data was standardized pre-modelling to make priors easier to choose as detailed in later sections.

Model specification

We adopted a Bayesian multitask model for predicting the 45 hepatosplenic conditions diagnosed in this study. Multitask models have often been presented as neural networks (e.g. [45–48]), but we proposed an architecture and inference framework that was arguably more general and suitable for clinical studies. Importantly, existing work in Bayesian multitask regression [49–51] lacked our level of interpretability. The model was inferred in a Bayesian manner such that the possible relationships between covariates and outcomes were incorporated more comprehensively as distributions [52], while we were able to maintain interpretations of significance and strength of associations. A Bayesian procedure was implemented to also improve the handling of imbalanced outcomes, as was the case for many severe conditions. Having such uncertainty in the model had been shown to be robust against class imbalance [52], mitigating the need for resampling methods such as SMOTE [53]. Importantly, we designed a new model architecture following the properties of multitasking, introducing ways for parameters to be exposed to multiple outcomes. All factors incorporated led to a final model with far superior predictive performance. Each step of model interpretation and utilization is summarized in Fig. 1.

Models were developed to predict M = 45 possible hepatosplenic conditions within an individual, participant data were indexed by i = 1,...,N. The outcome conditions were represented by indicator vectors y = (y ,...,y )⊤, such that y = 1 if person i has condition j, and 0 otherwise. Each participant was represented by (x ,y ), where we predicted the set of conditions y given the individual covariates x . We accounted for dependencies between the conditions through a thresholded multimorbidity graph G representing positive inter-relationships between the 45 hepatosplenic conditions. The graph was learned based on graphical lasso [54], considering the positive edges only with a threshold of 50% validated following the pipeline from previous study [10].

To progress towards multitask modelling, we started with single outcome logistic regression defined as

y = g(w +w x +···+w x )+ϵ = g(w⊤x ). (1)
ij 0 1 i1 iP P Ei i

In the case when M outcomes (M > 1) existed for analysis, we built M models for person i by

y = g(w⊤x )
i1 1 i
y = g(w⊤x )
i2 2 i
.
.
y = g(w⊤x ). (2)
iM M i

For more compact notation, we combined the regression weights/coefficients (w ) as

  
← w⊤ →
1
← w⊤ → 
y = g . 2 x  (3)
i  .  i
.
  
← w →
M
= g(W⊤x ). (4)
i

Model (4) was the basis of multi-output modelling [55], where it was important to note that y ∈ y continued to be predicted by w⊤x only.
ij i j i

We extended model (4) to a multitask framework by proposing the introduction of two new terms:

y = g(f (G)W⊤diag(Γ)x ). (5)
i α i

f is a function defined from the graph, and Γ = (γ ,...,γ )⊤ is a vector of probabilities for the P possible covariates. The roles of the two terms were to introduce multitask model structures as explained in the following sections.

Multitask graph convolution

We introduced the graphical function f (G) ∈ RM×M that convolved the models on each outcome. The aim was to replicate the known inter-relatedness of biological interactions of hepatosplenic conditions by introducing inherent relational inter-dependence between the models on each condition. Without this term, relations between the models would not be incorporated in the model. We chose the graphical function from signal processing [56, 57] to act as filtering functions to smooth node data based on the graph

f (G) = (I+αL)−1 (6)
α

where L was the regularized graph Laplacian defined by

L = D−1 (D−A)D−1 (7)
2 2

for adjacency matrix A and diagonal degree matrix D based on the graph, and α is an additional positive parameter controlling the influence of the graph.

By multiplying the multi-output model with f (G) ∈ RM×M, each outcome would be predicted by not only the corresponding model, but also the models on proximate nodes based on the graph structure. The choice of f allowed the construction of continuous weightings that accounted for both edge weights and graph distances. The weights assigned to the nodes were further controlled by the parameter α and optimized as part of the model inference. The ith rows of f (G) provided the weighting across the graph with the ith node as centre, ensuring model i primarily focused on condition i but drawing influence from models on proximate nodes. The size of the weights exhibited similar scale to probabilities, and were therefore normalized to one so that they could be interpreted as such. Condition i would be predicted by a weighted average of models, with the weighting on proximate nodes indicating the additional likelihood of developing condition i contributed from the other models. The value on node i itself would be interpreted as self-influence, reflecting solely the model developed for that condition and therefore the influence of the covariates on the singular condition. Isolated nodes in the graph represented unrelated conditions and would be treated as independent heterogeneous tasks. The graph filtering function (6) provided weightings for all nodes, incorporating more information than existing studies on multimorbidity networks which restricted analysis to 1-hop neighbourhoods. As the conditions came in groups based on different diagnoses of sonography (see first column of Table 1 for grouping), some conditions naturally co-occurred, therefore, our analysis also differentiated the two types of influences, with conditions from different groups classified as risk factors, while same group connections were interpreted as co-occurring conditions.

Multitask covariate learning

The term diag(Γ) ∈ R(P+1)×(P+1) performed the role of multitask covariate learning, with each γ acting as global inclusion probabilities of the covariates across all models. When combined with the covariates, diag(Γ)x = (γ ,γ x ,...,γ x )⊤ (assuming the first term in the covariate vector was the intercept). Each γ was inferred to represent the utility of covariate i for all outcomes as a joint distribution, allowing for covariate information sharing between models. For categorical variables, the categories were encoded as indicators, but we enforced the same inclusion probability across all indicators of the same variable. For instance, if x , x , and x were possible outcomes of x , then diag(Γ)x = (γ ,...,γ x ,γ x ,γ x ,...)⊤. Doing this provided a measure of relevance of x as a whole instead of each indicator independently.

While the inclusion probabilities acted as global multipliers shared across all models, inferring them in a Bayesian manner expanded their function as a parameterization of all the possible associations between the covariates against multiple outcomes. Both favorable (high probabilities) and uninformative (low probabilities) associations would be reflected in a single posterior distribution, which were fully incorporated for prediction. The shared nature of Γ would also influence the covariates effect size across tasks, covariates with many favourable associations would have their effect size increased on other conditions, and vice versa.

The role of inclusion probabilities could also be interpreted in the context of Bayesian variable selection, drawing similarities to the spike-and-slab prior when coupled with the regression coefficients. Each γ could potentially induce shrinkage effects on the covariates, and as a result, covariates were standardized before model fitting to ensure equal, comparable shrinkage, as consistently practiced with common Bayesian variable selection priors [58].

Prior selection

The multitask model took the following form:

      
y ← w⊤ → γ 0 1
1   N1 0
y   ← w⊤ → γ x 
 

. .
.
2  

= g 

 f α (G)   

I. .
.
2  

 

1 ...  

 

. .
.
1  

 

. (8)
y ← w → 0 γ x
M M P P

We henceforth refer to f (G) as the graph convolution, W as the regression coefficient matrix, and Γ as the inclusion probabilities.

The model was inferred in a Bayesian manner with all parameters given priors, and predictions were found by marginalization [52, 59]. First, inclusion probabilities were assigned priors with support between [0,1]:

i.i.d.
γ ,...,γ ∼ ContinuousBernoulli(θ = 0.5). (9)
0 P

The continuous Bernoulli was chosen due to it similarity to the binary counterpart when θ was close to 0 or 1 [60, 61], making each γ easier to interpret as inclusion/exclusion; Beta distributions of various shapes were also tested and led to similar results.

As regression weights w could be positive or negative, we used the common Gaussian distribution. From testing, we found the predictive quality was sensitive to the distributional variance, and therefore assigned priors to allow a varied variance

w i. ∼ i.d. N(0,σ2) (10)
ij ij
σ2 i. ∼ i.d. Γ(2,2). (11)
ij

The positive priors were assigned to the variance of each Gaussian distribution to allow the regression weights to have non-fixed variance to adapt to the uncertainty of the covariates. The Gamma distribution was chosen because of its non-negativeness and flexible shape, its hyperparameters were chosen to have mode 1, made appropriate by the normalization of the covariates. With a non-fixed variance, we found Gaussian priors performed similarly to t-distributions with prior-assigned scale parameter (equivalent to the variance of the Gaussian) and could be interchanged.

Concerning the graph convolution with non-negative parameter α, we assigned the same positive prior as previously described to provide a suitable starting point during inference

α ∼ Γ(2,2). (12)

Other shapes of Gamma distributions were tested, and we found priors with modal values within the range [1,3] led to similar predictive performances, but larger modes outside of the range led to predictions that were overly smoothed and degraded performance.

Bayesian inference

Bayesian machine learning is characterized by the inference procedure involving integrating over the uncertainty over all parameters. This was the key advantage of the Bayesian approach, as predictions were based not on single point estimates of the parameters, but Bayesian averaged over all possible parameter configurations weighted by their posteriors. Here, we detail the steps to carry out Bayesian inference on the model.

Model parameters were inferred using Markov Chain (Hamiltonian) Monte Carlo (MCMC) [62, 63] using the NUTS kernel from probabilistic programming library Pyro [64]. Let Θ = {Γ,W,σ ,α} represent all parameters of the multitask model f(x ) = g(f (G)W⊤diag(Γ)x ), and data D = {x ,y }N . We were interested in obtaining the posterior distribution

2 i α i i i i=1

P(Θ|f,ID) ∝ P(D|Θ,f)P(Θ|f) (13)
= P(D|Θ,f)P(Γ)P(W)P(σ2)P(α). (14)
(cid:124) (cid:123)(cid:122) (cid:125)(cid:124) (cid:123)(cid:122) (cid:125)
Likelihood Priors

The priors above were defined by (9), (10), (11), and (12). For binary classification, the likelihood was the Bernoulli density:

N M
(cid:89)(cid:89)
P(D|Γ,W,α,σ2) = f(x ) yij(1−f(x )) 1−yij (15)
i (j) i (j)
i=1 j=1

over N participants and M conditions, where f(x ) denoted the jth element of f(x ) ∈ RM. We assumed conditional independence of the outcomes in the likelihood, but statistical dependence was induced through the multitask parameters capturing interactions among outcomes. This setup allowed the likelihood to reflect outcome dependencies while retaining a computationally tractable likelihood. This was the optimal choice compared to multivariate logistic [65] or probit [66] likelihoods, which were intractable with the number of outcomes, and the dependence specification was through random effects as covariance structures instead of fixed effect of the multitask parameters.

Risk factor association

We interpreted the strength and direction of risk factor association by converting model posteriors to odds ratios. Each element of the model - graph convolution, regression coefficients, and inclusion probabilities - had individual posteriors, but we used the joint distribution of the three combined to reflect the scale of the data. Odds ratios were outcome specific such that for covariate j, we take the entry [f (G)W⊤diag(Γ)] to find the association with condition i. Division by the standard deviation of each covariate was applied such that the odds could be interpreted in the original units (undoing the data standardization prior to model inference).

More precisely, let m represent the ij-th entry of [f (G)W⊤diag(Γ)], and {mk}P represent the MCMC samples, we examined the empirical distribution of the odds ratio by the points {exp(mk/σ )}P where σ was the standard deviation of covariate j. When considering grouped conditions, we combined the MCMC samples for each condition to provide the overall effect.

The posterior densities were generally skewed and asymmetric, thus, we reported the median when inspecting individual covariates, but used the mean for overall effects. Credible intervals to 95% were computed using the highest posterior density [69], which were better equipped to handle skewed distributions than equal-tailed intervals.

Graph convolution under perturbation

The graph convolutional matrix was tested under perturbation of the graph. We considered two scenarios, first involving random addition and removal of edges, and adding random edges only. We tested proportions of edges randomly perturbed at 5%, 10%, and 20%, which led to a balance between noise and informative variation on the weights of the convolution. For edge addition, pairs of nodes were randomly selected and given edge weights that were randomly bootstrap sampled from the list of existing edge weights. Random perturbation was ran 100 times to produce the overall effect on the edges, we examined this mainly in the context of dependencies of gastro-oesophageal varices.

Model prediction and evaluation

Prior analysis had been performed using 100% of the data, but to compare the predictive performance of multitask modelling against other frameworks, we now divided participants into training and testing sets of sizes 50% each (1593). The sets were randomly sampled, and we completed this ten times to obtain an overall evaluation of the predictions. This procedure was chosen instead of cross validation because some conditions occurred infrequently, and for these cases, random sampling would lead to more splits with at least one positive case appearing in the test set, allowing the model to be evaluated on infrequent conditions more often.

We ran various versions of our multitask modelling, where in addition to the main model introduced, we also tested multitask using only the graph convolution corresponding to model (8) without the covariate learning term, and multitask using only covariate learning corresponding to model (8) without the graph convolution matrix. The main model in this study featured both multitask learning elements (8). Additionally, a number of machine learning models were tested as baselines: logistic regression using AIC variable selection on each outcome independently, XGboost for each outcome separately, multi-output random forest that modelled all outcomes jointly, and single and multitask shallow neural networks (as depicted in [18]), where the single task trained separate neural networks for each outcome while the multitask trained one neural network to predict 45 outcomes jointly.

For our Bayesian models, the training stage involved sampling from the posterior distributions P(Θ|D,f) using (Hamiltonian ) MCMC, providing us with a set of posterior sample points for all parameters. For each split we ran 500 warm-up steps followed by 1000 posterior sampling steps. To make predictions given the test participants, the MCMC samples were used to compute the Monte Carlo estimator (19).

The model outputs were probabilities, so evaluations were done firstly using area under the receiver operating curve (AUC), and AUC-precision-recall (AUC-PR), as well as binary-based metrics F1 (binary and macro), precision, recall, positive and negative predicted values (PPV and NPV), and balanced accuracy after conversion. For binary-based metrics, we considered choosing the thresholds on the predicted probabilities using two different clinical measures: 1. the threshold leading to the highest F1 score was used to compute the precision, recall, and balanced accuracy (both the binary and macro versions in sklearn were evaluated), 2. the threshold leading to 90% sensitivity was used to evaluate PPV and NPV.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Zhi and coworkers present a manuscript with an impressive and comprehensive approach to assess hepatosplenic multimorbidtiy in an ugandan population, including ultrasound-data. Overall, the manuscript is written very well, the analysis was conducted reliably and is offering new insights. The manuscript should be encougrad to be accepted after several revision: A major issues seems to be alochol induced hepatosplenic multimorbidity as outlined below Issues: -hepatosplenic multimorbidity should be explained in further detail. It is recommend, for the readers comfort, to include in the introduction the spectrum of diseases (that are prevalant in the patienst from this region) included in this defintion (e.g. liver cirrhosis and hepatosplenic multimorbidity) -How was liver fibrosis diagnosed/defined (line 103) -Line 308„burst gastro-oesophageal varices“. More commonly the term osophageal variceal bleeding is used. -Are any markers and scores for cirrhosis (e.g. FIB4-Score) and portal hypertension (e.g. Platelet-count) available? If possible, they should be included in the analysis. -The outcome and data on the specific disease entities in hepatlosplenic multimorbidity should be analysed on each entity (Schistomiasis, hepatitis B, liver cirrhosis etc.) in more detail and if possible, analysed separately. Quite surprisingls current alcohol use was not associated with any hepatosplenic conditions Even though, it seems the only information on alcohol consumption was current use yes or no (respectively „based on self- reported consumption within the year preceding recruitment“) How was current alcohol use defined (How many drinks per day over which time period? Is data availabe on past alcohol comsumption? Alcohl consumption yes or no does not seem to be sufficient enough. How prevalent is alcohol abuse in this population? Could further data be provided? Could further alcohol-specific analyses of this dataset provide further insight? More information is need on this matter. This might, if no further information is available be a strong limitation of this analysis (as it is already included in line 291/292) and should be explained in detail, how this confounder is not severely limiting the conclusions drawn from the analysis. -Are we sure what endpoints that are clinically relevant were assessed (hepatosplenic multimorbidity is a term not really helpfull for a clinician taking care for the patients, as it contains quiete different diseass entities)? This should be further discussed in the manuscript Can it be distinguished in the analysis which patients had hepatosplenic multimorbidity caused by cirrhosis (e.g. due to alcohol abuse or hepatitis B) and which portal hypertension due to schistosomiasis, as well as other causes? More comments on this matter should be included.", "reviewer_2": "Summary This manuscript presents an innovative Bayesian multitask learning framework for modeling 45 hepatosplenic conditions in a large schistosomiasis cohort. The integration of graph-informed condition dependencies, global covariate inclusion probabilities, and fully Bayesian inference is technically sophisticated and well motivated. The dataset is unique and clinically valuable, and the approach yields meaningful biological and epidemiological insights. The work has strong potential impact, but several aspects must be strengthened for publication. In particular, while the authors emphasize interpretability, which is appropriate for the clinical domain, the methodological evaluation needs improvement. Broader comparative analyses and more imbalance-aware performance metrics would make the contribution significantly more compelling. Major Comments 1. Performance comparison should go beyond interpretability-focused baselines The authors explicitly position their method as an interpretable alternative to opaque black-box models. This is a valid and important goal. However, the paper would be considerably more impactful if it also demonstrated where the proposed model stands relative to commonly used multitask or multilabel learning approaches. Even if the authors prefer simpler or more interpretable models, showing how their Bayesian architecture compares to at least a subset of widely used baselines, such as multi-output Random Forests, XGBoost, classifier chains, or a shallow shared-bottom neural network, would provide crucial context. These comparisons need not detract from the interpretability focus; rather, they would strengthen the overall contribution by showing either (i) competitive performance with substantially greater interpretability, or (ii) superior performance in the setting of rare outcomes. This addition would broaden the paper’s relevance and demonstrate the practical advantages of the proposed Bayesian framework beyond conceptual elegance. 2. AUC alone is not appropriate given the extreme class imbalance Many of the modeled conditions occur in less than 5% of individuals, and several in less than 1%. In these cases, ROC-AUC is known to be an insufficient metric, as it remains high even when precision for rare events is very low. Since the authors argue that their Bayesian model handles imbalance more effectively, this claim should be supported with metrics designed for imbalanced classification. The manuscript would benefit from reporting precision–recall curves, PR-AUC (or average precision), precision, recall, F1-score, and balanced accuracy. These metrics are standard in modern evaluation of imbalanced data and are essential for understanding model behavior in clinically relevant scenarios. 3. Threshold-based evaluation is needed for clinical relevance Because the authors position the model as a clinical decision-support tool, threshold-based evaluation is essential. Clinicians need to know how reliable a prediction is when the model actually issues it. Metrics such as Positive Predictive Value (PPV) (when the model predicts a patient is positive, how likely is it that the patient truly has the condition?) and Negative Predictive Value (NPV) (when the model predicts a patient is negative, how likely is it that the patient truly does not have the condition?) are crucial. These should be reported at clinically meaningful operating points, such as sensitivities fixed at 90% or 95%. This information determines whether the model’s improvements translate into practical utility, especially for rare conditions where AUC can be misleading despite poor PPV. To make this concrete, consider a condition with <1% prevalence. A model may achieve an AUC of 0.90 yet still have a PPV of only a few percent at clinically practical thresholds, meaning that the vast majority of “positive” predictions would be false alarms. Without PPV, NPV, the false-positive burden at fixed sensitivities, and appropriate precision–recall evaluation, it is impossible to determine whether the reported performance improvements or the conditional relationships inferred by the model translate into meaningful diagnostic value. 4. Graph-derived conditional dependencies must be explicitly validated and interpreted While the Discussion briefly relates some inferred condition relationships to known hepatosplenic pathways, this evaluation remains qualitative. Given that the graph-convolution component is central to the model’s interpretability claims, the manuscript should more explicitly and systematically analyze whether the conditional relationships inferred by the model align with established clinical mechanisms, such as known sequences of fibrosis progression, portal hypertension development, and splenic involvement. It would be valuable for the authors to indicate where the learned dependencies support existing pathophysiological understanding and where they diverge. Without such analysis, it is difficult to determine whether the smoothing imposed by the graph enhances biological realism or primarily reflects statistical correlation structure. 5. Evidence for generalizability would strengthen the contribution While external validation datasets for schistosomiasis are scarce, a geographically stratified validation (training on two districts and testing on the third) would provide valuable insight into model transportability. If this is not feasible, some robustness analysis (bootstrapping, district-weighted resampling, or sensitivity to graph structure) would still enhance the reliability of the findings. 6. Data availability and reproducibility While I understand the ethical and privacy constraints around participant-level data, it would be helpful for the authors to clarify whether any form of anonymized, de-identified, or synthetic dataset could be shared to facilitate reproducibility. If full individual-level data cannot be released due to ethics board restrictions or data governance rules, the authors might still consider providing a minimally de-identified feature matrix (removing geographic, demographic, and other potentially identifying variables) or a synthetic dataset generated to match the statistical properties of the original cohort. Either option would enable readers to reproduce the modeling pipeline and evaluate the proposed method more effectively. Minor Comment Figure 3, including its schematic and legend, should be more explicative. This figure introduces important components of the model but currently requires substantial inference from the reader. Additional labels or a clearer narrative description would improve clarity. Conclusion and Recommendation This manuscript addresses a highly relevant public health problem with an original, well-designed Bayesian modeling framework. The approach is conceptually strong, and the epidemiological insights are of genuine interest. Incorporating broader performance baselines, imbalance-aware metrics, clinically grounded threshold analysis, and explicit interpretation of graph-derived dependencies will significantly strengthen the manuscript. I therefore recommend Major Revision.", "reviewer_3": "This study deploys an interpretable Bayesian multitask learning framework to jointly model 45 hepatosplenic ultrasound outcomes in a rural Ugandan cohort. By combining clinical tests + sociodemographic + access/locational covariates with a graph-convolution over an outcome-dependency network and continuous inclusion probabilities for covariates, the authors both (i) improve discrimination—especially for rare, severe outcomes—and (ii) surface clinically plausible risk patterns (older age↑; low Hb↑; female sex protective; higher risk among fishers). The work is methodologically thoughtful and practically relevant to triage in resource-limited settings.So it is suggested to make some revisions with following comments: 1.Primary estimation & interpretability. Most substantive findings in the text derive from per-outcome logistic regressions (ORs), while the proposed multitask model is judged mainly by AUC. If interpretability is central, the authors should report posterior log-odds (converted to ORs) with 95% HPD intervals for key covariates within the multitask model (e.g., f_α (G) W^(⊤\"diag\" ) (Γ) mapped to ORs after reversing standardization). Please also explain how \\Gamma (continuous inclusion) relates to effect sizes and whether high-\\Gamma variables consistently exhibit larger posterior odds shifts across tasks. 2. Validation design (80/20 and beyond). Validation design should also be expanded and clearly described. If the data are split 80/20, please report repeated stratified 80/20 evaluation over many resamples; additionally, include a geographical transport test (leave-one-district-out) and, if feasible, a temporal validation to probe stability over time. For rare outcomes, present PR-AUC, F1-macro, and confusion matrices at clinically meaningful thresholds. Because clinical deployment hinges on threshold choices, decision-curve analysis would add practical value by quantifying net benefit across plausible risk thresholds. 3. Multimorbidity patterns (go beyond “better AUC”). Since multimorbidity modeling is a stated novelty, I encourage the authors to move beyond “better AUC” and explicitly characterize comorbidity patterns induced by the model. Identify covariate scenarios that maximize predicted multimorbidity counts (e.g., low Hb combined with older age and longer distance to facilities); present patient-level archetypes by clustering predicted outcome vectors; and summarize which outcomes “donate” the most information to others via the graph by reporting row-wise weights from the graph convolution. Interpreting these pathways in light of hepatosplenic pathophysiology would substantiate the claim that the model surfaces clinically actionable “danger signals.” Finally, robustness and reporting can be strengthened. Please conduct graph sensitivity analyses by perturbing edge weights/thresholds or re-learning the graph and assessing stability of performance and danger-signal rankings; perform prior sensitivity for the graph-strength and variance components; specify missing-data handling; and document prevalence and any class-imbalance strategies."}, "author_response": "Zhi and coworkers present a manuscript with an impressive and comprehen- sive approach to assess hepatosplenic multimorbidtiy in an ugandan popu- lation, including ultrasound-data. Overall, the manuscript is written very well, the analysis was conducted reliably and is offering new insights. The manuscript should be encougrad to be accepted after several revision: A major issues seems to be alochol induced hepatosplenic multimorbidity as outlined below Response We appreciate the feedback and recognition of the comprehensive nature of our approach. We have incorporated all suggestions in our revisions, which have improved the clarity from a clinical perspective of our manuscript. Comment: hepatosplenic multimorbidity should be explained in further detail. It is recommend, for the readers comfort, to include in the introduction the spec- trumofdiseases(thatareprevalentinthepatientsfromthisregion)included in this definition (e.g. liver cirrhosis and hepatosplenic multimorbidity) Response: Thank you for pointing out that a more direct definition would have been helpful in the introduction. We have added to the second paragraph the following (with citations in text) \"The most common hepatosplenic conditions that occur in schistosomiasis-endemic populations include periportal fibrosis, early cirrhosis-like or hepatitis-like livers, splenomegaly, moderately shrunken or enlarged livers, and thick- ened gall bladder walls–all of which are often accompanied by clinical symptoms such as diarrhea, abdominal pain, anaemia, jaundice, and melena. Hepatosplenic multi- morbidity due to chronic schistosome infections encompasses two or more interacting or co-occurring conditions ranging from combinations of less severe forms of peripor- tal fibrosis and changes in spleen or liver organometry to inclusion of life-threatening complications of gastro-oesophageal varices and liver cirrhosis.\" Comment: How was liver fibrosis diagnosed/defined (line 103) Response: In this particular sentence, liver fibrosis was defined as periportal fibrosis and specifically Niamey protocol pattern C for the diagnosis. We have clarified in text and ensured the Niamey protocol is cited in the Methods. We note in the last paragraph of the introduction (as well as in the methods) that diagnosis was by B-mode point- of-care ultrasound. To adhere to journal guidelines, we keep additional details to the Methods that comes after the results. Comment: Line 308 \"burst gastro-oesophageal varices\". More commonly the term osophageal variceal bleeding is used. Response: Corrected. 1 Comment: Are any markers and scores for cirrhosis (e.g. FIB4-Score) and portal hyper- tension (e.g. Platelet-count) available? If possible, they should be included in the analysis. Response: Unfortunately, as this study relied on point-of-care ultrasound, we did not have any biomarkers available for analysis. We now mention this single diagnostic as a limitation in the Discussion and that the inclusion of biomarkers such as platelets and liver enzymes (e.g. AST or ALT) may better enable severity grading. Comment: The outcome and data on the specific disease entities in hepatosplenic mul- timorbidity should be analysed on each entity (Schistomiasis, hepatitis B, liver cirrhosis etc.) in more detail and if possible, analysed separately. Response: Thank you for this comment. We provide models on each condition entity and their risk factors for any condition that has at least one significant risk factor in Table 3, we also have the breakdown for the individual condition risk factors separated by adults and children in Tables S3 an S4 respectively. There is a detailed presentation of these results in the second paragraph of the first Results section. We also now add results to indicate what combination of risk factors (e.g. infection, demographics, etc) produce the highest multimorbidity counts to give more detail on the multimorbidity. Comment: Quite surprisingly current alcohol use was not associated with any hep- atosplenic conditions Even though, it seems the only information on alcohol consumption was current use yes or no (respectively \"based on self-reported consumption within the year preceding recruitment\") How was current alco- hol use defined (How many drinks per day over which time period? Is data available on past alcohol consumption? Alcohol consumption yes or no does not seem to be sufficient enough. How prevalent is alcohol abuse in this population? Could further data be provided? Could further alcohol-specific analyses of this dataset provide further insight? More information is need on this matter. 2 Thismight, ifnofurtherinformationisavailablebeastronglimitationofthis analysis (as it is already included in line 291/292) and should be explained in detail, how this confounder is not severely limiting the conclusions drawn from the analysis. Response:We agree that this is an important area for future research, and as such have initiated a detailed investigation in the next annual timepoint for SchistoTrack to collect information on community consumption patterns, alcohol purchase behaviours, local definitions of alcohol use (both identified socially through the community and by clinicians handling these patients), and more information on type, quantity, and frequency. This information, though not at the individual level, would provide high level information to explain the alcohol use patterns. Comment: Are we sure what endpoints that are clinically relevant were assessed (hep- atosplenic multimorbidity is a term not really helpful for a clinician taking care for the patients, as it contains quite different disease entities)? This should be further discussed in the manuscript Response: We agree that the clinical endpoint of focus is important for patient man- agement. In this respect, we provide a breakdown of individual conditions and we focus throughout the manuscript on the example of gastro-oesophageal varices and its concurrent/interacting conditions. We add this further justification to the discussion paragraph we discuss the results of gastro-oesophageal varices. \"We focused on this example due to the need of having a concrete clinical endpoint for patient case man- agement in schistosomiasis-endemic areas. High mortality rates of 24% often are ob- served across low and high income contexts due to oesophageal variceal bleeding. Even higher mortality rates are assumed in low-income contexts where the burden of gastro- oesophageal bleeding has been poorly documented and resources for endoscopy, surgical management, beta blockers, and blood transfusion are severely limited.\" Comment: Can it be distinguished in the analysis which patients had hepatosplenic multimorbidity caused by cirrhosis (e.g. due to alcohol abuse or hepatitis B) and which portal hypertension due to schistosomiasis, as well as other causes? More comments on this matter should be included. Response: Thank you for this suggestion, which initiated additional exploration. We looked at grouping the conditions on the graph using traditional methods of clustering (community detection) and found that there were three distinct groups/clusters. We have added results to the Graph convolution section as follows. \"To give an alternative view on condition clusters using the graph convolution network, we show the grouping of conditions found through community detection on the graph (Fig. S24). Three main clusters were found. One cluster mainly contained measurements of liver enlargement, fatty livers, and other abnormalities not specifically attributable to infections and po- tentially due to alcohol use. Another cluster that included cirrhosis had conditions that represent features, complications, of consequences of cirrhosis that may be attributable to alcohol use, HIV, or HBV such as rounded liver edges, shrunken liver lobes, and gross 3 liver surface irregularities. The third cluster encompassed mostly schistosomal peripor- tal fibrosis Niamey patterns, its features such as changes in the main portal vein, and portal hypertension consequences including gastro-oesophageal varices and splenorenal shunts.\" We also added the following to the discussion. \"Community detection on the graph also found gastro-oesophageal varices belonged to the same cluster as all schistosomal periportal fibrosis patterns (Niamey protocol patterns C - F). Interestingly, this cluster appeared clearly separated from a cluster on cirrhosis and its features and complications as well as another cluster that appeared to entail fatty livers and gall bladder thickening. Future prospective work is needed to assess whether these separate clusters may repre- sent portal hypertension conditions more attributable to schistosomiasis versus HIV or HBV (cirrhosis cluster) vs alcohol use (fatty cluster).\" Though, we do not want to over-interpret these clusters as we think despite the dis- tinction being ideal, it is best left to prospective analyses where issues of temporality and sequence of exposures can be dealt with to avoid reverse causality. For the por- tal hypertension that may be due more specifically to schistosomiasis, we do try to make an argument in the Discussion showing how periportal fibrosis patterns (focus- ing on schistosomiasis-specific liver fibrosis patterns) are the most important condition risk factors for gastro-oesophageal varices (a clinical indication of portal hypertension). But, overall, we hypothesize (as has been shown elsewhere in Anjorin et al. 2024 Lancet Microbe) that we think the cirrhosis and portal hypertension actually may arise due to interactions between hepatitis b, HIV, and schistosomiasis. Again, these interactions would need to be untangled in prospective work. We mention this in a few places in the Discussion but now make it clear in the limitations. \"Clinically, the inferred dependen- cies between conditions suggest biological relatedness, but they do not imply causality without further prospective study, and therefore cannot be used to explain liver fibro- sis pathogenesis or the origin of portal hypertension being due solely to initial single infections or co-infection interactions over time.\" 4 Reviewer 2 Comment This manuscript presents an innovative Bayesian multitask learning frame- work for modeling 45 hepatosplenic conditions in a large schistosomiasis cohort. The integration of graph-informed condition dependencies, global covariate inclusion probabilities, and fully Bayesian inference is technically sophisticated and well motivated. The dataset is unique and clinically valu- able, and the approach yields meaningful biological and epidemiological in- sights. The work has strong potential impact, but several aspects must be strength- ened for publication. In particular, while the authors emphasize inter- pretability, which is appropriate for the clinical domain, the methodologi- cal evaluation needs improvement. Broader comparative analyses and more imbalance-aware performance metrics would make the contribution signifi- cantly more compelling. Response: We really appreciate the comments on the innovative aspects of the frame- work, and clinical value of the data and interpretive analysis. We are also thankful that the reviewer drew our attention to imbalance-aware metrics that will greatly improve the comparative analysis of the paper. All suggestions have been incorporated and have strengthened the manuscript. Comment: 1. Performance comparison should go beyond interpretability-focused base- lines The authors explicitly position their method as an interpretable alternative to opaque black-box models. This is a valid and important goal. However, the paper would be considerably more impactful if it also demonstrated where the proposed model stands relative to commonly used multitask or multilabel learning approaches. Even if the authors prefer simpler or more interpretable models, showing how their Bayesian architecture compares to at least a subset of widely used baselines, such as multi-output Random Forests, XGBoost, classifier chains, or a shallow shared-bottom neural network, would provide crucial context. These comparisons need not detract from the interpretability focus; rather, they would strengthen the overall contribution by showing either (i) competitive performance with substantially greater interpretability, or (ii) superior performance in the setting of rare outcomes. This addition would broaden the paper’s relevance and demonstrate the practical advantages of the proposed Bayesian framework beyond conceptual elegance. Response: We are thankful for the reviewer for emphasizing the broader context of performance comparison. We have added a number of model comparisons, which have further strengthened our results to show that our proposed model has very strong pre- dictive capabilities relative to other models. We introduced XGBoost and multi-output 5 random forest, and shallow deep learning based models of multi and single task learning as described in Caruana [18]; standard logistic regression remains another core compar- ison for our model given its standard use in epidemiology. We compared the models by inspecting the newly added Table 4 in the manuscript, which contains a number of relevant metrics accounting for class imbalance and relevant for clinical interpretation; some additional less informative metrics can be found in Table S7. These tables replaced the original Fig. S6 in the main paper and Results section \"Multimorbidity prediction\" has been updated with results as described below. TheperformanceofourBayesianmultitaskmodelanditsvariantsstillremainedamongst the best, with only the multitask neural network being competitive based on a number of the metrics including AUC, AUC-PR, binary F1 and precision, and PPV. On other metricssuchasMacroF1andnegativepredictedvalue(NPV),therewerelessnoticeable differences as all models performed similarly well. It is worth highlighting the compet- itive performance from the multitask neural network, which also supports our claim of the advantage of joint modelling. On individual conditions, we added model compar- isons on AUC-PR and PPV in addition to AUC. Fig. 3, which previously compared AUCs of individual conditions, has been changed to AUC-PRs-to-prevalence ratio and compared to the prevalence of the conditions. As there were many models, comparisons between models for individual conditions are left to Supplementary Results. To focus the reader’s attention, we chose to compare our main model against logistic regression (most interpretable) and multitask neural networks (best performing baseline) in the main text–evaluated using AUC (Fig. S7 - S8), AUC-PR (Fig. S9 - S10), and PPV (Fig. S11 - S12). We agree that AUC can be less informative on imbalanced data and are now mainly examining AUC-PR and PPV. Our model had higher AUC-PRs on 82.22% (37/45, three significant) of conditions than logistic regression, and 73.33% (33/45, one significant) of conditions against multitask neural network. Meanwhile, the PPVsofourmodelwashigherthanlogisticregressionon84.44%(38/45, foursignificant) of conditions, and higher on 66.67% (30/45, one significant) of conditions compared to multitaskneuralnetwork. Theseresultsindicatethatourmodelwasbetteratpredicting each individual condition by a small margin, rather than being significantly better on a select few. Even though the overall performance is very similar on imbalance aware metrics, logistic regression and any neural network cannot provide the rich interpretable elementsofourmodel. Nobaselinemodelscancommentonconditiondependence, while neural networks and other black-box models cannot identify risk factors, and strength and direction of association. Comment: 2. AUC alone is not appropriate given the extreme class imbalance Many of the modeled conditions occur in less than 5% of individuals, and several in less than 1%. In these cases, ROC-AUC is known to be an insuf- ficient metric, as it remains high even when precision for rare events is very low. Since the authors argue that their Bayesian model handles imbalance more effectively, this claim should be supported with metrics designed for imbalanced classification. The manuscript would benefit from reporting precision–recall curves, PR- AUC (or average precision), precision, recall, F1-score, and balanced accu- racy. These metrics are standard in modern evaluation of imbalanced data and are essential for understanding model behavior in clinically relevant sce- narios. 6 Response: We thank the reviewer for pointing this out as we agree that AUC alone does not uncover the full behaviour of our models. As noted in the response to the previous comment, we have added a number of model comparisons. Along with those models, we added metrics of AUC-PR, F1 (binary and macro), precision, recall, PPV, NPV, and balanced accuracy. For reference on the reported AUC-PR, the prevalence of all outcomes is 0.046. We found the AUC-PR achieved by our model to be a good level for such low prevalences, as this has been observed previously, such as the example in [1] (reference only in response document) where the authors worked with a disease of 4% prevalence and off-the-shelf models achieved similar AUC-PRs. Thus, we have chosen to base our main comparisons using AUC-PR instead of AUC. For threshold based metrics, a threshold was chosen based on the what has been considered clinically meaningful: precision and recall were computed based on the threshold that led to the highest F1 score for binary and macro respectively; F1 was chosen here because it provides a balanced measure of both metrics. PPV and NPV are determined using a second threshold selection method based on 90% sensitivity, which is further detailed in the next response. Comment: 3. Threshold-based evaluation is needed for clinical relevance Because the authors position the model as a clinical decision-support tool, threshold-based evaluation is essential. Clinicians need to know how reliable a prediction is when the model actually issues it. Metrics such as Positive Predictive Value (PPV) (when the model predicts a patient is positive, how likely is it that the patient truly has the condition?) and Negative Predictive Value (NPV) (when the model predicts a patient is negative, how likely is it that the patient truly does not have the condition?) are crucial. These shouldbereportedatclinicallymeaningfuloperatingpoints, suchassensitiv- ities fixed at 90% or 95%. This information determines whether the model’s improvements translate into practical utility, especially for rare conditions where AUC can be misleading despite poor PPV. To make this concrete, consider a condition with <1% prevalence. A model may achieve an AUC of 0.90 yet still have a PPV of only a few percent at clinically practical thresholds, meaning that the vast majority of “positive” predictions would be false alarms. Without PPV, NPV, the false-positive burdenatfixedsensitivities, andappropriateprecision–recallevaluation, itis impossible to determine whether the reported performance improvements or the conditional relationships inferred by the model translate into meaningful diagnostic value. Response: We have introduced PPV and NPV along with a number of metrics in Table 4, using a second threshold corresponding to 90% sensitivity in the predictions as the reviewer pointed out the practical utility in measuring false alarms. All models appear to handle negative outcomes particularly well, and mainly differ in the PPV. We noticed that models that involved some form of joint modelling produced better PPVs, as we can see in Table 4 that the multitask neural network and multi-output random forest produced significantly higher PPVs than any single-task/separate modelling. Our Bayesian multitask models still remain amongst the best performing by a small but statistically significant margin. For precision and recall, we describe what was carried out in the previous response: AUC-PR represents the average precision across all recall values and is evaluated on the raw probabilities; the AUC and AUC-PR curves are 7 presented in Fig. S13. For exact precision and recall, which require converting predicted probabilities to binary, the F1 score was used as the objective to select the threshold as it provides a balance between precision and recall. We then considered both the binary F1 focusing on the positives, and the macro F1 to provide a balance between the precision on the positive and negative outcomes. Comment: 4. Graph-derived conditional dependencies must be explicitly validated and interpreted While the Discussion briefly relates some inferred condition relationships to known hepatosplenic pathways, this evaluation remains qualitative. Given thatthegraph-convolutioncomponentiscentraltothemodel’sinterpretabil- ity claims, the manuscript should more explicitly and systematically analyze whether the conditional relationships inferred by the model align with es- tablished clinical mechanisms, such as known sequences of fibrosis progres- sion, portal hypertension development, and splenic involvement. It would be valuable for the authors to indicate where the learned dependencies support existing pathophysiological understanding and where they diverge. Without such analysis, it is difficult to determine whether the smoothing imposed by the graph enhances biological realism or primarily reflects statistical correlation structure. Response: Thank you for highlighting the importance of graph validation. We first added Table S9, showing pairs of nodes with the strongest probabilities in the graph convolution matrix, which also represent the strongest inter-dependencies in the graph. A paragraph of the clinical analysis on these edges has been added to the Discussion, where we have interpreted the biological correctness of Table S9, discussing a combina- tion of biological complications, disease progression, and data acquisition as compared to the existing biological knowledge in the literature though no relationships were ex- perimentally validated. The graph is learned from data and therefore inherently reflects the statistical relationships between the conditions. Efforts were made in using higher order distributional statistical measures to capture the most important relationship to best represent true relationships. The work in [10] does go into case analysis, visualiz- ing a number of subgraphs to allow close examination of individual conditions and their related morbidities. Comment: 5. Evidence for generalizability would strengthen the contribution While external validation datasets for schistosomiasis are scarce, a geograph- ically stratified validation (training on two districts and testing on the third) would provide valuable insight into model transportability. If this is not fea- sible, somerobustnessanalysis(bootstrapping, district-weightedresampling, or sensitivity to graph structure) would still enhance the reliability of the 8 findings. Response: Thank you for the suggestions for new experiments. We have conducted these new experimentsandaddedtothemtothepaper. First, weevaluatedourBayesianmultitask model on district-based data splitting, training on two district and predicting on the third, which we present in the new Table S8. Given that district was a significant covariate, we expect the population from each district to have strong shared risk factors andtheretobeobviousdistrictdifferences. ItalsoshouldbenotedthatSchistoTrackwas designedtoincludepurposelydifferentdistrictsintermsofriskfactorsandprevalenceso it is unclear whether a district validation shows generalisability or whether it is better to keep a model trained on diverse districts then let it see other settings that are contained in the diversity it is trained on. As expected, we find that splitting by district influences predictive performances. The result in Table S8 suggests mixed trends but generally Mayuge is the best predicted district with the highest average AUC-PR and PPV. But, this performance may be due to having most of the cases (highest condition prevalences) in the other two (Western districts) so for the metrics that considered imbalance there were most positive cases to be seen during training. These results are now included in the section \"Multimorbidity prediction\" and in the a Discussion paragraph. We have also conducted graph perturbation to measure the variability of the influence probabilities from the graph convolutional matrix. We tried random addition of edges, and addition and removal; when edges are added, pairs of nodes are randomly selected and given weights that are bootstrap sampled from the list of existing edges. We showed what happens to the convolution matrix weights around gastro-oesophageal varices with 5%, 10%, and 20% edge removal, when both adding and removing are performed and the total number of edges remains the same. This experiment was ran a total of 100 times, and we present these in newly added Fig. S20 - S22. The graph convolution was very robust to graph perturbation especially at low percentages (5% & 10%) with very little variability to the convolution weights, while the unperturbed weights nearly always remained within the inter-quartile range from perturbation even at 20%. Thus, we concluded the graph was robust to noise, and this point was added to the Discussion paragraphs on the validity of the graph. Comment: 6. Data availability and reproducibility While I understand the ethical and privacy constraints around participant- level data, it would be helpful for the authors to clarify whether any form of anonymized, de-identified, or synthetic dataset could be shared to facilitate reproducibility. If full individual-level data cannot be released due to ethics board restrictions or data governance rules, the authors might still consider providing a minimally de-identified feature matrix (removing geographic, de- mographic, andotherpotentiallyidentifyingvariables)orasyntheticdataset generated to match the statistical properties of the original cohort. Either option would enable readers to reproduce the modeling pipeline and evaluate the proposed method more effectively. Response: For testing purposes, we now provide generated data by random sampling. For continuous covariates, we sampled from normal distributions with means and stan- 9 dard deviations equal to each of the variables (the original covariates had log transforms applied if they became closer to normal). For binary covariates, we sample from a Bernoulli using the mean of the variables. We also did this for each of the outcomes, using the prevalence of each condition. We re-uploaded these source files to allow the code to run, however we are unable to reproduce the results of the paper because this would involve generating data in a way that mimics the original data, which would be challenging especially for multiple outcomes and beyond the scope of this study— whereas the main aim of providing the anonymised code is to enable other studies to test our models on new datasets. Comment: Figure 3, including its schematic and legend, should be more explicative. This figure introduces important components of the model but currently requires substantial inference from the reader. Additional labels or a clearer narrative description would improve clarity. Response: To make the pipeline figure clearer without reading the methods, we have addedmoreadetaileddescriptioninthecaptionexplainingwhateachpaneldoes. Given the number of interpretive elements to the model, we completely agree that this needed more explanation to somewhat stand alone without reading all of the methods. Remark on code availability: I examined the shared materials (README.md and pyro_mcmc.py). Al- thoughthecoreBayesianmultitaskmodelisimplemented, thecurrentrepos- itorydoesnotallowreaderstoreproducetheresultsreportedinthemanuscript. The main script depends on a non-included module (data_loader) and on participant-level data that are not provided, so the pipeline cannot be ex- ecuted end to end. The example at the bottom of the script does not run successfully either. It attempts to compute AUROC values but references a variable (test_id) that is undefined when no train–test split is created, meaning that the example cannot produce the metrics shown in the paper. Theunderlyingprobabilisticmodelappearscorrectlyimplemented,butwith- outaworkingexample,withoutanyformofanonymizedorsyntheticdataset, and without an evaluation pipeline that mirrors the manuscript, the results are not currently reproducible. To meet the reproducibility expectations of the journal, the authors should provide a functional example using either anonymized or synthetic data, ensure that AUROC values and other metrics can actually be computed and saved, and expand the evaluation to include imbalance-aware metrics. A more complete README and more robust co- variate handling would also improve clarity and usability. Response: The imported module “dataloader\" contained information about our inter- nal cloud structure and participants which was why it could not be shared. Its purpose was to load the condition outcomes, the rest of the “pyro_mcmc.py\" script showed load- ing of the covariates, which were then joined with the outcomes. We wanted to show that the variables in the code corresponded to the ones described in the paper. The variables “test_id\" and “train_id\" were global variables that populate the session once you run “run_mcmc\" or “run_mcmc_id\", we have changed this to explicitly defining them in the example usage. We now share the revised code that is much more reduced 10 with all data preparation removed, and will simply import the synthetic files. Note, the synthetic data does not replicate the actual data but is provided to allow the code to run smoothly for anyone learning to use our models. We cannot share the individual data or properties of it that can be relearned to identify communities due to our data protection impact assessment, ethics approvals, and ongoing nature of the cohort. References 1. Garriga, R. et al. Machine learning model to predict mental health crises from electronic health records. Nature medicine 28, 1240–1248 (2022). 11 Reviewer 3 This study deploys an interpretable Bayesian multitask learning framework to jointly model 45 hepatosplenic ultrasound outcomes in a rural Ugandan cohort. By combining clinical tests + sociodemographic + access/locational covariates with a graph-convolution over an outcome-dependency network andcontinuousinclusionprobabilitiesforcovariates, theauthorsboth(i)im- prove discrimination—especially for rare, severe outcomes—and (ii) surface clinically plausible risk patterns (older age↑; low Hb↑; female sex protective; higher risk among fishers). The work is methodologically thoughtful and practically relevant to triage in resource-limited settings. So it is suggested to make some revisions with following comments: Response: We appreciate the positive assessment of this work from the reviewer, and are particularly thankful for them recognizing the strength of the methodology and the practical application. We were happy to make clarifications and the revisions suggested by the reviewer. Comment: 1. Primary estimation & interpretability. Most substantive findings in the text derive from per-outcome logistic re- gressions (ORs), while the proposed multitask model is judged mainly by AUC. If interpretability is central, the authors should report posterior log- odds (converted to ORs) with 95% HPD intervals for key covariates within the multitask model (e.g., f (G)W⊤diag(Γ) mapped to ORs after reversing α standardization). Please also explain how Γ (continuous inclusion) relates to effect sizes and whether high-Γ variables consistently exhibit larger posterior odds shifts across tasks. Response: We agree with the reviewer for pointing out that more interpretive statistics should be highlighted and more attention drawn to when these materials are provided in the supplement. Firstly, we showed the posterior distributions of the odds ratios (posterior log-odds) for gastro-oesophageal varices, presented in supplementary Fig. S2. These are plotted on the log-scale and the 95% HPD intervals are also marked out for each covariate. All reported median and credible intervals in the results are based on these histograms. A set of 28 plots exist for each of the 45 conditions (so 28x45=1260 plots); such a large number of plots would have been difficult to include in a manuscript.As an example, we only displayed these plots for one condition (gastro-oesophageal varices) that we focused our analysis on. We also interpreted these plots in our results and discussion. We picked outthesignificantvariablesinageand(log)Hbconcentration, withthemedianposterior odds ratios indicating 4% increase per year of age, and 22% increase per 10% decrease in log-transformed Hb concentration. In Table 3, all odds ratios are also presented with 95%HPDintervals. Onedifferenceisthesignificanceofthecovariatesaredeterminedvia q-values from corrected p-values from likelihood ratio tests, instead of being determined by the HPD intervals containing odds ratio 1. The corrected q-values are a stricter measure to reduce false-positives and we can see that all significant variables had HPD intervals not including odds ratio 1. All odds ratios, both in the plots and tables, are 12 presented after reverse standardization and detailed this in the section \"Risk factor association\" in Methods. The inclusion probabilities Γ alone do not determine the effect size, but do exhibit strong correlations with the size of odds ratios across all conditions. Their function is to spread information about the utility of the covariates. If a covariate is a very good predictor for a number of conditions, then it would lead to high γ, which will increase the influence of that covariate on the rest of the outcomes. As a result, the effect size of the covariate could be increased on outcomes that previously did not have a strong relationship. However the effect size is dependent on all terms in f (G)W⊤diag(Γ), not α just Γ, and each covariates’ effect size can still be reduced by the W terms. Comment: 2. Validation design (80/20 and beyond). Validation design should also be expanded and clearly described. If the data are split 80/20, please report repeated stratified 80/20 evaluation over many resamples; additionally, include a geographical transport test (leave- one-district-out) and, if feasible, a temporal validation to probe stability over time. For rare outcomes, present PR-AUC, F1-macro, and confusion matrices at clinically meaningful thresholds. Because clinical deployment hingesonthresholdchoices,decision-curveanalysiswouldaddpracticalvalue by quantifying net benefit across plausible risk thresholds. Response: We apologize for the lack of clarity on the data splitting. The results are presented in two aspects, first, for explainability purposes, we obtained all statistics using the full 100% of the data which Table 3, S1 - S6, and Fig. 2, S1, - S4, and S6 - S11 were derived from; any report of medians and credible intervals were derived from 100% of the data. We now make this information clear (what data were used to obtain the results) in the caption of each table of figure. In he last section of the results, we reported testing the predictive ability of the model, and here we split the data into training and testing sets. We performed random splits of 50:50 for training and testing, and we resampled this way 10 times to produce standard deviations on the test AUCs. We chose 50:50 split because for rare conditions this will lead to more test sets that contain at least one positive and can be evaluated (note the model was able to learn from training data even if only one class existed). Fig. 3 and S5 show predictive performances and are therefore based on this method of splitting. We also clarify in the figure captions that this split was used. While we understand temporal analysis will add significant insights to hepatosplenic multimorbidity, we are not able to complete this here due to our study being cross- sectional and only considered one year of data. Modelling of temporal patterns are currently ongoing within the cohort. And, our discussion highlights the importance of prospective designs in future work (in the Limitations section). For the geographical prediction test, we have now evaluated our Bayesian multitask model on leave-one-out based on districts, we present this in the new Table S8. Given that district was a significant covariate, we expect the population from each district to have strong shared risk factors and there to be obvious district differences. It also should be noted that SchistoTrack was designed to include purposely different districts 13 in terms of risk factors and prevalence so it is unclear whether a district validation shows generalisability or whether it is better to keep a model trained on diverse districts then let it see other settings that are contained in the diversity it is trained on. As expected, we find that splitting by district influences predictive performances. The results suggest mixed trends but generally Mayuge is the best predicted district with the highest average AUC-PR and PPV. But, this performance may be due to having most of the cases (highest condition prevalences) in the other two (Western districts) so for the metrics that considered imbalance there were most positive cases to be seen during training. We understand the need for imbalance-aware metrics, and have since evaluated our modelsonalistofmetricsincludingAUC-PR,binaryandmacroF1, aswellasprecision, recall, and positive and negative predicted values (PPV & NPV). Importantly, we have changed our main predictive analysis to comparisons between AUC-PR and prevalence, replacing the original AUC comparisons. For binary based metrics, we have used two methods of choosing the threshold: first, the threshold leading to the highest F1 score, which we then used to evaluate the precision and recall of the model, as F1 provides a balance of the model precision and recall. We did this for both the binary F1 focusing on just the positive outcome, and the macro F1 which provides a balance between the positive and negative classes. For PPV and NPV, we set a second cut-off that leads to 90% sensitivity. Being a clinically informative metric, PPV then informs the likelihood of the patients truly having the condition having been predicted positive. Comment: 3. Multimorbidity patterns (go beyond “better AUC”). Sincemultimorbiditymodelingisastatednovelty, Iencouragetheauthorsto move beyond “better AUC” and explicitly characterize comorbidity patterns induced by the model. Identify covariate scenarios that maximize predicted multimorbidity counts (e.g., low Hb combined with older age and longer distance to facilities); present patient-level archetypes by clustering predicted outcomevectors; andsummarizewhichoutcomes“donate” themostinforma- tion to others via the graph by reporting row-wise weights from the graph convolution. Interpreting these pathways in light of hepatosplenic patho- physiology would substantiate the claim that the model surfaces clinically actionable “danger signals.” Response: We thank the reviewer for suggesting the many different ways to interpret multimorbidity patterns and have expanded this element of the paper. First, we have addedDiscussionpointsontheresultsonoverallmultimorbidityfromthelastparagraph of \"Risk factors for hepatosplenic conditions\" section. The risk factors of overall multi- morbidity were used to predict the probabilities of having all conditions, this means the sum of the probabilities represented the expected number of conditions, which informs the morbidity count. Therefore, we can determine the scenarios that would maximize expected morbidity count based on covariate significance and directions of the asso- ciations to \"Overall\" in Table S2; this includes low Hb concentration, older age, and four other factors. We also introduced Fig. S23 where we computed the probabilities of having every condition if a patient were to have this scenario of the covariates that maximizes morbidity counts, the setting is as follows: positively associated covariates (e.g. age) were assigned the 99th percentile value from our data, negatively associated 14 covariates (e.g. Hb concentration) were assigned the 1st percentile value, insignificant covariates were assigned the 50th percentile (median). To measure the pairwise interactions between the outcomes in the model, we introduced TableS9todocumentthehighestweightsbetweennodepairsfromthegraphconvolution matrix, representing the outcomes that donate the most information to each other. The interpretation of the biological correctness is based on the pairs highlighted in this table, and a paragraph is added in the Discussion which we would like to direct the reviewer to. Because of the nature of the problem tackled in this study, we like to clarify that cluster- ingcanbeperformedontheoutcomesortheparticipants,butthereisnotaclearconnec- tionbetweenthetwo. Therefore,wecannotcommentonwithin-patientarchetypesbased on clusters of conditions. We performed Louvain community detection on the graph to produce clusters of conditions and presented these in Fig. S24. We found a central cluster was particularly consistent with our analysis as it grouped gastro-oesophageal varices with all fibrosis patterns C - F, confirming its relation with periportal fibrosis. Other clusters also reflect biological patterns, with there existing a cluster containing all forms of enlargements on both liver lobes, while shrunkenness/restrictedness of both liver lobes were clustered together with various shrunkenness in the spleen and portal vein. Comment: Finally, robustnessandreportingcanbestrengthened. Pleaseconductgraph sensitivity analyses by perturbing edge weights/thresholds or re-learning the graphandassessingstabilityofperformanceanddanger-signalrankings; per- form prior sensitivity for the graph-strength and variance components; spec- ifymissing-datahandling; anddocumentprevalenceandanyclass-imbalance strategies. Response: We have firstly conducted graph perturbation and measured the variability of the influence probabilities from the graph convolutional matrix. We tried random addition of edges, and addition and removal; when edges are added, pairs of nodes are randomly selected and given weights that are bootstrap sampled from the list of existing edges. We set the proportions of edges to add/remove at 5%, 10%, and 20%, when both adding and removing are performed the total number of edges remain the same. This was run a total of 100 times. We measured the range of values from the graph convolution around gastro-oesophageal varices, which we presented in newly added Fig. S20 - S22. The graph convolution function we used was very robust to graph perturbation especially at low percentages (5% & 10%) there were very little variability to the convolution weights, while the original unperturbed convolution weights nearly always remained within the inter-quartile range even at 20%. This finding is integrated in the Discussion along with the interpretation on the strongest graph connections to support the validation of the multimorbidity graph. For prior selection, we have tested various options before settling on the priors specified in our paper. For the α parameter in the graph convolution, we tested Gamma priors Γ(α,β) with the mode values at 1, 3, and 5 (mean = α/β), the performance was similar for mode = 1 and 3, but worse for mode = 5, suggesting it should be a small value, and a mode of 1 was used in our final model. Additionally, the inclusion probabilities γ was the most robust to priors with variously shaped Beta distributions tested that 15 i all led to almost identical posteriors. For the regression coefficients, normal and t priors were tested and did not show noticeable difference, but the performance was sensitive to the variance/scale terms of the normal/t distribution. This led to us selecting the normal and putting a prior on the variance to allow it to vary. We have added these explanations to the \"Prior selection\" section in Methods, and displayed the posteriors using an alternative set of priors in Fig. S3 (see caption for priors used), to show that as long as the prior selection were reasonable for each model component the results are extremely robust. Participants are selected for this study only if they have all required data. Of the participants examined with point-of-care ultrasound, only 2.14% (69/3224) had missing ultrasound data or covariates, therefore we decided it was reasonable to exclude them from the dataset. This is detailed in Participants selection section of Methods, and a participants flowchart in Supplementary which details the number of missingness at each stage of the data collection. We have documented the prevalence of every outcome condition in Table 1, the number of participants having each condition is listed with the prevalence in bracket. There is a strong imbalance especially for a number of severe conditions, and this partly mo- tivated the model architecture and framework. Firstly, being a fully Bayesian model, the use of priors act as regularization on the model parameters. Bayesian inference also incorporates the distributions of the parameters in the prediction through marginaliza- tion and Bayesian averaging, instead of using any point estimate that could potentially lead to overfitting. Architecture-wise, the multitask framework induces parameter shar- ing between the models for each outcome, thereby all parameters are exposed to more outcomes, and in turn more positive outcomes. Furthermore, we have introduced a number of imbalance-aware metrics including AUC-PR, F1, precision, recall, and PPV to provide clearer evaluations of the predictions. 16"}, {"round": 2, "reviewer_comments": {"reviewer_1": "After the corrections, the manuscript can now be recommended for publication", "reviewer_2": "I thank the authors for their careful and thorough responses to the reviewers’ comments. The revisions have satisfactorily addressed the previously raised concerns, and the manuscript has been substantially improved in clarity, methodological transparency, and interpretation of the results. In its current form, the study is sound, well presented, and makes a meaningful contribution to the field. I have no further comments and recommend the manuscript for publication. Minor comment: In Tables 4 and S7 (Model performance comparison), the columns appear to lack sufficient spacing or margins, which makes some values difficult to read. I recommend adjusting the table formatting to improve readability. (Remarks on code availability) The authors have made the code publicly available and runnable, with adequate documentation. Due to privacy and governance constraints associated with the clinical data, the original dataset cannot be shared. To address this, the authors provide a synthetic dataset that allows users to run the pipeline and verify code functionality. While the published results are therefore not fully reproducible using the shared data, the implementation and workflow are transparent and reproducible, and the code constitutes a usable resource for the community.", "reviewer_3": "The revisions have significantly enhanced the clarity and rigor of the analysis.The authors have provided more detailed explanations, particularly regarding the definitions and diagnostic methods for hepatosplenic multimorbidity, which is greatly appreciated. The clarification on alcohol consumption, the inclusion of additional validation steps, and the expanded discussions on model performance metrics further strengthen the manuscript. I am satisfied with the revisions made and agree that the manuscript is now much clearer and better aligned with clinical relevance. All of the points raised in my previous review have been adequately addressed. (Remarks on code availability) Yes code provided can be run, in which the main implementation was instored in pryo_mcmc_synth.py file. The code provide a README file,with enough instructions，including: node and edge lists as csv, for constructing the multimorbidity graph required for the model, and synthetic condition outcomes and covariates csv, to allow the code to run."}, "author_response": ""}]

Decision letter