Development and validation of predictive models for the hemodialysis adequacy index (Kt/V) using multiple linear regression and machine learning algorithms

1 Introduction

Chronic kidney disease (CKD) has an increasing prevalence in China and throughout the world, and a recent study estimated the global prevalence rate was 14.3% (1). Thus, it is necessary to provide better treatments for CKD to increase the survival time and improve the quality of life in these patients (2). CKD is most common in patients with diabetes or high blood pressure, and it can progress at different rates to end-stage renal disease (ESRD) (3). Renal replacement therapy (RRT), including hemodialysis (HD), peritoneal dialysis, and kidney transplantation, is the only effective treatment for ESRD. Due to its convenience, safety, and effectiveness, HD has become the main method of RRT for patients with ESRD (4).

Adequate HD can be broadly defined (5) as achieving and maintaining a good clinical status through dialysis treatment (based on improved blood pressure, intravascular volume, cardiac function, hemoglobin level, appetite, and quality-of-life), and more narrowly defined as a satisfactory clearance rate of small-molecule solutes. Kt/V (known as the “urea clearance index” or the “hemodialysis adequacy index”), is the most commonly used quantitative indicator of the clearance rate of small-molecule solutes (6). Inadequate HD is associated with clinical complications, poor quality of life, and poor survival rates (7, 8). Therefore, it is important for clinicians to assess dialysis adequacy in patients receiving maintenance hemodialysis (MHD), and Kt/V is widely used as an indicator to formulate and adjust dialysis treatment plans.

Currently, there are two internationally recognized methods for determining Kt/V. The first method is based on the urea kinetic model proposed by Gotch and Sargent, which calculates Kt/V by measuring the pre- and post-dialysis plasma urea nitrogen levels from blood samples (9). However, this method requires collection of multiple blood samples, which can exacerbate anemia and patient discomfort, and can have an unacceptable economic burden. The results from a sample after a single dialysis session can be inadequate for formulation of a dialysis plan using this method. The second method, online clearance monitoring, uses parameters from the dialysis machine that are related to urea clearance (10). However, many HD machines cannot record these data, and no large-scale clinical studies have validated the accuracy of machine-based online detection. Therefore, the clinical practice guidelines for CKD (NKF-K/DOQI, 2006) do not recommend use of components built into dialysis machines for detection of dialysis adequacy. For these reasons, both of these methods can be problematic in clinical settings.

Clinical prediction models, defined as statistical or machine learning models that are based on the phenotypic characteristics of a disease, can predict the probability of the future occurrence of an outcome or a specific numerical value or indicator in a population with certain clinical features (11). Studies of other diseases utilized similar models to predict the risk of disease onset and progression, and treatment outcomes (11, 12). The purpose of the present study was to develop and validate predictive models for HD adequacy (Kt/V). We employed multiple linear regression (MLR), seven machine learning algorithms for point prediction (XGBoost, LightGBM, Random Forest, Gradient Boosting, KNN, SVM, and BPNN), an ensemble model combining three of these algorithms, and two probabilistic prediction algorithms (LightGBM Quantile Regression and NGBoost). The predictive performance of all models was then compared using internal and external validation sets. Our aim was to develop a reliable model for prediction of Kt/V that can be used by practicing clinicians, by helping them to formulate HD treatment plans in a more timely and scientifically informed manner, and thereby improve patient quality-of-life and long-term survival.

2 Data and methods

2.1 Data preparation

2.1.1 Study design and study population

This was a multicenter, prospective observational study conducted across seven hemodialysis (HD) centers in Hohhot City, Inner Mongolia, China, from January 2022 to September 2022. The participating centers included the Affiliated Hospital of Inner Mongolia Medical University, Hohhot First Hospital, Inner Mongolia People's Hospital, Togtoh County Hospital, Inner Mongolia Armed Police General Hospital, Inner Mongolia Traditional Chinese Medicine Hospital, and the First Hospital of Inner Mongolia Autonomous Region Prison Administration.

The study was approved by the ethics committee of the Affiliated Hospital of Inner Mongolia Medical University (S.2021086) and was performed in accordance with the ethical standards in the 1964 Declaration of Helsinki and its later amendments. Written informed consent was obtained from all participants or their legal representatives prior to enrollment.

A total of 1,200 consecutive maintenance hemodialysis (MHD) patients were initially screened. Inclusion criteria were: (1) age between 18 and 80 years; (2) diagnosis of chronic kidney disease (CKD) requiring MHD; (3) clear consciousness and ability to communicate normally; (4) willingness to participate and provision of informed consent. Exclusion criteria were: (1) presence of language, visual, auditory, or intellectual impairment preventing participation; (2) age <18 or >80 years; (3) New York Heart Association Functional Class IV; (4) frequent recent angina attacks; (5) resting oxygen saturation <90%; (6) any other serious illnesses; or (7) refusal to participate. Patients with missing or incomplete data were excluded from the final analysis. After applying the inclusion and exclusion criteria, 1,017 eligible patients were enrolled.

2.1.2 Data collection and outcome definition

The following demographic and clinical information was collected: gender, age, height, ethnicity, pre-dialysis weight, post-dialysis weight, HD vintage, dialysis duration per session, dialysis frequency, dialyzer membrane surface area and material, dialysate sodium concentration, blood flow rate, dialysate flow rate, access type, anticoagulant type, primary disease, and laboratory parameters including pre-dialysis and post-dialysis serum urea, hemoglobin (Hb), albumin (ALB), triglycerides (TG), total cholesterol (TC), C-reactive protein (CRP), platelets, parathyroid hormone (PTH), serum phosphorus, serum calcium, and serum potassium.

Serum urea concentration was measured before and after dialysis and the data were used to measure Kt/V. Before dialysis, blood was collected from the arterial end of the dialysis circuit without heparin or saline; after dialysis, ultrafiltration was stopped and the blood was collected from the arterial end of the dialysis circuit after the blood flow rate decreased to 100 mL/min for 15 s.

The primary outcome variable was Kt/V, calculated using the second generation Daugirdas formula (Equation 1):where R is the ratio of pre-dialysis to post-dialysis serum urea, t is the dialysis duration (hours), UF is the ultrafiltration volume (liters), and W is the post-dialysis body weight (kg).

Kt/V=−ln(R−0.008×t)+(4−3.5×R)×UF/W(1)

2.2 Dataset partitioning and data preprocessing

2.2.1 Dataset partitioning

Eligible patients were divided into three independent datasets following a stratified approach: (1) Training set (n = 677): used for feature selection, hyperparameter tuning, and final model training; (2) Internal validation set (n = 220): strictly withheld from all training and tuning procedures; used only once after model finalization for initial unbiased performance assessment; (3) External validation set (n = 120): sourced from three different institutions not represented in the training set; used as the ultimate independent test set to evaluate generalizability.

This strict three-way split ensured that the internal and external validation sets were never exposed to any model development step, providing an unbiased evaluation of model performance.

2.2.2 Missing data imputation

Missing values in continuous predictor variables were imputed using the median value calculated from the training set only. Missing values in categorical variables were imputed using the mode from the training set. The same imputation values derived from the training set were applied to the internal and external validation sets to prevent data leakage.

2.2.3 Feature standardization

To ensure compatibility with distance-based and gradient-based algorithms (SVM, KNN, BPNN), all continuous predictor variables were standardized using Z-score normalization.

Specifically, the mean (μ) and standard deviation (σ) for each continuous feature were calculated exclusively from the training set. These training-set parameters were then applied to transform the training, internal validation, and external validation sets (Equation 2):This approach simulates a real-world deployment scenario, where new patients are standardized using fixed parameters derived from the development data, and prevents any information leakage from the validation sets.

x_standardized=(x−μ_train)/σ_train(2)

2.3 Feature engineering

2.3.1 Several constructed derived features

To capture clinically meaningful information beyond the raw measurements, we constructed several derived features from the original variables: All original and constructed features were considered as candidate predictors for subsequent screening.

• Circulating blood volume (mL) = actual blood flow rate (mL/min) × actual dialysis duration (h) × 60;
• Standardized blood flow (mL/min/m²) = actual blood flow rate (mL/min)/body surface area (m²);
• Ultrafiltration rate (L/h) = actual ultrafiltration volume (L)/actual dialysis duration (h);
• Pulse pressure (mmHg) = pre-dialysis systolic BP—pre-dialysis diastolic BP;
• Mean arterial pressure (mmHg) = diastolic BP + (systolic BP—diastolic BP)/3;
• Lipid accumulation index = total cholesterol (mmol/L) × triglycerides (mmol/L);
• Nutritional index = albumin (g/L)/hemoglobin (g/L);
• Weight to BSA ratio = pre-dialysis weight (kg)/body surface area (m²).

2.3.2 Feature selection

We employed a multi-algorithm fusion strategy to identify the most robust set of predictors, integrating eight complementary feature selection methods that span different methodological paradigms: For each candidate variable, the ranks from the first seven algorithms were averaged to generate a consensus importance ranking. The Boruta algorithm categorizes each candidate variable into one of three classes, confirmed, rejected, or tentative, based on its importance relative to shadow features, allowing the investigator to decide whether a variable should be retained for subsequent model development. This ensemble-based approach ensures that a variable is retained only if it is consistently important across fundamentally different methodological paradigms, linear vs. non-linear, filter vs. embedded vs. wrapper, and tree-based vs. model-agnostic.

• F test (filter method, for linear associations);
• Mutual Information (filter method, for non-linear associations);
• Random Forest importance (embedded method, tree-based with interaction capture);
• XGBoost importance (embedded method, gradient-boosting with gain-based importance);
• LASSO coefficients (embedded method, L1-regularized linear model for sparse selection);
• Permutation Importance (model-agnostic diagnostic, assessing drop in performance after shuffling);
• RFE (Recursive Feature Elimination) (wrapper method, iteratively removing the least important features);
• Boruta (wrapper method, random forest-based with shadow feature comparison) (13).

2.3.3 Outlier removal

To prevent extreme values from unduly influencing model training, we performed outlier removal based on the Kt/V distribution in the training set. Specifically, samples with |Z| > 3 (i.e., beyond 3 standard deviations from the mean) in the training set were identified as extreme outliers and excluded prior to model training. This step was performed after feature selection to ensure that the model learned the dominant physiological patterns rather than fitting to extreme cases.

2.4 Model development for point prediction

2.4.1 Candidate models

We developed and compared eight different predictive models using the five selected features: All models were implemented in Python 3.10 using scikit-learn, XGBoost, and LightGBM libraries.

• MLR (Multiple Linear Regression), served as the baseline linear model;
• XGBoost, an extreme gradient boosting algorithm based on decision trees;
• LightGBM, a lightweight gradient boosting framework with leaf wise tree growth;
• RandomForest, a classical bagging ensemble of decision trees;
• GradientBoosting, a sequential boosting algorithm with decision trees as base learners;
• KNN (K Nearest Neighbors), a distance-based non parametric regression method;
• SVM (Support Vector Machine), with a radial basis function (RBF) kernel for non linear regression;
• BPNN (Back Propagation Neural Network), a multi layer perceptron with two hidden layers.

2.4.2 Hyperparameter tuning

Hyperparameter optimization was performed using GridSearchCV with 5-fold cross-validation on the training set. The objective function was to maximize the average R² across the 5 folds.

For each model, we pre-defined a search space of candidate hyperparameters. The full hyperparameters search spaces are provided in Supplementary Table S1.

Key examples of the tuned hyperparameters included: After the optimal hyperparameters were identified via GridSearchCV, the final model for each algorithm was retrained on the full training set using these optimal parameters.

• XGBoost: n_estimators, max_depth, learning_rate, subsample, colsample_bytree, gamma, reg_lambda, min_child_weight
• SVM: C, gamma, epsilon (RBF kernel)
• BPNN: hidden_layer_sizes, activation, alpha, learning_rate_init
• KNN: n_neighbors, weights, p

For boosting-based models (XGBoost, LightGBM, GradientBoosting), early stopping was applied during training to prevent overfitting: training was terminated if the validation loss did not improve for 30 consecutive rounds.

2.4.3 Ensemble model

In addition to the individual models, we constructed a simple average ensemble combining the predictions of the top performing models (SVM, GradientBoosting, and XGBoost), which were selected based on their superior test set RMSE and MAE.

The ensemble prediction was calculated as (Equation 3):This ensemble approach leverages the complementary strengths of geometrically based (SVM) and tree based (GradientBoosting, XGBoost) models to achieve more robust and stable predictions.

y^_ensemble=(y^_SVM+y^_GradientBoosting+y^_XGBoost)/3(3)

2.4.4 Model calibration assessment

To assess whether the predictions from the three selected individual models (SVM, GradientBoosting, and XGBoost) and their simple average ensemble exhibited systematic bias (i.e., consistent over- or under-estimation), we performed a calibration diagnosis on the internal validation set.

A linear calibration model was fitted for each candidate model (Equation 4):where ŷ_raw is the original model prediction, a is the calibration slope, and b is the calibration intercept. A perfectly calibrated model would have a slope of 1 and an intercept of 0. The residual distribution (residual = true Kt/V−predicted Kt/V) was examined using: (1) residual-vs.-predicted scatter plots; (2) histograms of residuals; (3) true-vs.-predicted scatter plots; and (4) Q-Q plots for normality.

y^_calibrated=a×y^_raw+b(4)

To quantify the potential benefit of post-hoc calibration, we calculated the improvement in MAE (mean absolute error) and the change in R² after applying the linear calibration transformation. A model was considered well-calibrated and not requiring post-hoc adjustment if the MAE improvement was less than 0.005, as this threshold was deemed clinically negligible for Kt/V prediction.

2.4.5 Model interpretability (SHAP analysis)

To explain the contribution of each feature to the model predictions, we performed SHAP (SHapley Additive exPlanations) analysis for the top performing models (SVM and GradientBoosting).

SHAP values quantify the marginal contribution of each feature to the prediction for each individual patient, providing both global feature importance rankings and local explanations for individual predictions.

2.5 Performance evaluation for point prediction models

Model performance was assessed using the following metrics: (1) R² (Coefficient of Determination), the proportion of variance in Kt/V explained by the model; (2) RMSE (Root Mean Squared Error), the square root of the average squared prediction error, sensitive to large errors; (3) MAE (Mean Absolute Error), the average absolute prediction error, interpretable in the original Kt/V units; (4) MAPE (Mean Absolute Percentage Error), the average percentage error. All metrics were computed separately for the internal validation set (initial unbiased assessment) and the external validation set (final generalizability assessment).

To evaluate overfitting, we compared performance between the training set (from cross validation) and the validation sets. A small performance gap indicates that the model does not suffer from substantial overfitting.

2.6 Development and performance evaluation for probabilistic prediction models

To address the editor's request for prediction uncertainty estimation, we implemented two probabilistic prediction methods capable of generating prediction intervals alongside point estimates: The performance of probabilistic models was evaluated using: (1) PICP (Prediction Interval Coverage Probability), the proportion of true values falling within the predicted interval (ideal: close to the nominal confidence level); (2) MPIW (Mean Prediction Interval Width), the average width of the prediction intervals (narrower is better); (3) CWC (Coverage Width Criterion), a composite metric balancing PICP and MPIW (lower is better); (4) CRPS (Continuous Ranked Probability Score), a strictly proper scoring rule evaluating the entire predictive distribution (lower is better); (5) Pinball Loss, for quantile regression models; (6) NLL (Negative Log Likelihood), for NGBoost.

• LightGBM Quantile Regression, direct modeling of conditional quantiles (0.05, 0.50, 0.95) using the quantile regression objective. This approach makes no distributional assumptions and can handle asymmetric error distributions.
• NGBoost (Natural Gradient Boosting), a probabilistic gradient boosting framework that models the parameters of a log normal distribution for log transformed Kt/V. It outputs a full probability distribution from which prediction intervals and probability density curves can be derived.

2.7 Bridging post-hoc development and pre-dialysis application

A critical point warrants explicit clarification regarding the intended clinical use of our model. During model development, we used post-dialysis parameters (actual dialysis duration, actual blood flow rate, and actual ultrafiltration volume) because these were the only objectively documented values available in our retrospective dataset.

In routine clinical practice, however, the dialysis prescription is determined by the clinician before the session begins, i.e., the planned dialysis duration, planned blood flow rate, and planned ultrafiltration volume are entered into the dialysis machine prior to treatment initiation. Importantly, due to standardized treatment protocols and stable machine performance, the prescribed parameters are virtually identical to the delivered actual parameters in the vast majority of treatments. In our dataset, the mean absolute difference between prescribed and actual values was negligible (<2% for all three parameters), confirming this consistency.

Therefore, although the model was developed using post-hoc actual values, its intended real-world application is pre-dialysis prediction, i.e., clinicians can input the planned prescription parameters (which are known before dialysis starts) into the identical model to obtain a reliable prediction of the Kt/V that would be achieved if that prescription is followed. This enables real-time, pre-emptive adjustment of the dialysis prescription, for example, increasing blood flow rate or extending dialysis duration, before the session actually begins, a functionality that the Daugirdas formula (which requires post-dialysis urea measurements by blood sampling and testing) cannot provide.

2.8 Development of a clinical application interface

To facilitate clinical adoption, we developed a local standalone executable application that implements the final ensemble model and the NGBoost probabilistic predictor. The application was packaged using PyInstaller and is distributed as a self-contained file package that can be executed offline without network connection. The interface accepts five input parameters (sex, height, pre-dialysis weight, blood flow rate, dialysis duration) and returns point predictions, prediction intervals, and probability density curves.

2.9 Statistical analysis and software

All statistical analyses and model development were performed using Python 3.10 (Python Software Foundation) with the following core libraries: pandas 2.3.3, numpy 2.2.6, scikit learn 1.7.2, xgboost 3.2.0, lightgbm 4.6.0, boruta 0.4.3, matplotlib 3.10.9, seaborn 0.13.2, shap 0.49.1, and properscoring 0.1.

Descriptive statistics are presented as mean ± standard deviation for normally distributed continuous variables, median (Q1, Q3) for non-normally distributed variables, and frequencies (percentages) for categorical variables. The Shapiro–Wilk test was used to assess normality.

All code used for data preprocessing, feature selection, model training, hyperparameter tuning, and evaluation are available. The final trained models are provided from the corresponding author upon reasonable request.

3 Results

3.1 Study population and baseline characteristics

A total of 1,017 eligible patients were enrolled from seven hemodialysis centers. The training set comprised 677 patients (66.6%), the internal validation set 220 patients (21.6%), and the external validation set 120 patients (11.8%). Figure 1 shows the flow chart for disposition of patients.

The median Kt/V values were 1.33 (Q1–Q3: 1.14–1.56) in the training set, 1.31 (Q1–Q3: 1.11–1.51) in the internal validation set, and 1.31 (Q1–Q3: 1.11–1.51) in the external validation set. The training set had a median age of 55.0 years (Q1–Q3: 44.0–65.0), with 396 (58.5%) male patients. The median dialysis vintage was 42.0 months (Q1–Q3: 15.0–87.8). The mean body surface area (BSA) was 1.81 ± 0.19 m² in the training set. Detailed baseline characteristics of the three datasets are presented in Table 1.

Table 1. Baseline characteristics of participants.

Normally distributed continuous data are presented as mean ± standard deviation, non-normally distributed continuous data as median (Q1, Q3), and categorical data as n (%).

P < 0.05 for comparison with the training set.

3.2 Feature engineering and selection

From the original 26 candidate variables, we constructed eight clinically meaningful derived features: circulating blood volume, standardized blood flow, ultrafiltration rate, pulse pressure, mean arterial pressure, lipid accumulation index, nutritional index, and weight to BSA ratio.

Using a multi-algorithm fusion strategy that integrated rankings from seven different feature selection methods (F test, Mutual Information, Random Forest importance, XGBoost importance, LASSO coefficients, Permutation Importance and Recursive Feature Elimination), we generated a consensus importance ranking of all candidate features. Besides, the Boruta algorithm provided the importance classification for each candidate variable (confirmed, rejected, or tentative), based on its importance relative to shadow features. Table 2 shows the ranking of feature selection for the top ten features and the ranking of feature selection for all the features are shown in Supplementary Table S2. Based on the average rank and clinical expert review, five features were selected as the final input variables for model development: (1) Body surface area (BSA, m²); (2) Pre-dialysis weight (kg); (3) Standardized blood flow (mL/min/m²); (4) Sex; (5) Circulating blood volume (mL).

Table 2. Feature selection ranking across eight algorithms for the top ten features.

F, F test; MI, Mutual information; RF, Random forest importance; XGB, XGBoost importance; LASSO, LASSO coefficients; Perm, permutation importance; RFE, recursive feature elimination.

The feature “weight-to-BSA ratio” was excluded due to high collinearity with BSA and pre-dialysis weight, and because its importance rank was lower than that of the selected variables.

For the multiple linear regression (MLR) model, the final fitted Equation 5 was

Kt/V=1.305−0.276×BSA(m2)+0.010×predialysisweight(kg)−0.240×standardizedbloodflow(ml/min)+0.110×sex(1forfemale,0formale)+0.276×circulatingbloodvolume(ml)(5)

3.3 Hyperparameter tuning outcomes

Hyperparameter optimization was performed using GridSearchCV with 5-fold cross validation on the training set, with the objective of maximizing the average R² across folds. The optimal hyperparameters for each model are summarized in Table 3. The full hyperparameters search spaces are provided in Supplementary Table S1.

Table 3. The optimal hyperparameters for each model.

For XGBoost, the optimal configuration was n_estimators=100, max_depth=5, learning_rate=0.05, subsample=0.8, colsample_bytree=0.8, gamma=0.1, reg_lambda=1.0, and min_child_weight=1. For SVM with RBF kernel, the best parameters were C=50, gamma=0.01, and epsilon=0.2. For the BPNN, the optimal architecture consisted of two hidden layers with 5 and 3 nodes respectively (hidden_layer_sizes=(5,3)), activation='tanh', alpha=0.05, and learning_rate_init=0.005. For KNN, the optimal n_neighbors was 9 with weights='uniform' and p=2.

3.4 Model performance and model selection

3.4.1 Model performance on internal validation

All eight individual models and the ensemble model were evaluated on the internal validation set (n = 220), which was strictly withheld from all training and tuning procedures. Performance metrics on internal validation are presented in Table 4.

Table 4. Performance metrics of point prediction on internal and external validation sets.

Val refers to internal validation set and Test refers to external validation set. Ensemble (Avg) refers to a simple average ensemble combining the predictions of the top performing models (SVM, GradientBoosting, and XGBoost), which are selected based on their superior test set RMSE and MAE.

On the internal validation set, the SVM model achieved the highest R² (0.4633), the lowest RMSE (0.2258), and the lowest MAE (0.1512) among individual models. The GradientBoosting model achieved an R² of 0.4126, RMSE of 0.2362, and MAE of 0.1648. The XGBoost model achieved an R² of 0.4293, RMSE of 0.2329, and MAE of 0.1628. The MLR baseline model achieved an R² of 0.3337, RMSE of 0.2516, and MAE of 0.1574.

The simple average ensemble (combining SVM, GradientBoosting, and XGBoost) outperformed all individual models on the internal validation set, achieving an R² of 0.4478, RMSE of 0.2291, and MAE of 0.1582. The BPNN, LightGBM, RandomForest, and KNN models demonstrated intermediate performance, with R² values ranging from 0.3873 to 0.4197 on the internal validation set.

3.4.2 Model performance on external validation

To assess generalizability to completely unseen patient populations from different clinical centers, all models were evaluated on the external validation set (n = 120). Performance metrics on external validation are presented in Table 4. On the external validation set, all models demonstrated improved performance compared to the internal validation set, which may be attributed to the smaller standard deviation of Kt/V in the external validation cohort (0.257 vs. 0.309 in the internal validation set).

The SVM model continued to perform best among individual models, achieving an R² of 0.5931, RMSE of 0.1633, and MAE of 0.1214. The GradientBoosting model achieved an R² of 0.5748, RMSE of 0.1670, and MAE of 0.1276. The XGBoost model achieved an R² of 0.5681, RMSE of 0.1683, and MAE of 0.1288.

The simple average ensemble achieved the highest overall performance on the external validation set, with an R² of 0.5911, RMSE of 0.1637, and MAE of 0.1240, demonstrating robust generalizability across different clinical settings. Notably, the MLR baseline model also performed competitively on the external validation set (R²=0.5389, RMSE=0.1739, MAE=0.1278), suggesting that the selected linear features capture substantial predictive information.

3.4.3 Model selection for clinical translation

To select the most suitable models for clinical implementation, we ranked all individual models by their RMSE and MAE on the external validation set. The SVM, GradientBoosting, and XGBoost models consistently ranked among the top three on both metrics. These three algorithms were therefore chosen as the final model set for subsequent application development, as they demonstrated the best trade-off between predictive accuracy and generalizability across different clinical centers. Moreover, their simple average ensemble achieved even higher performance on the external validation set (R² = 0.5911, RMSE = 0.1637), and was thus also incorporated into the final application interface, alongside the three individual models.

3.5 Calibration assessment

To determine whether the predictions from the three selected individual models and their ensemble exhibited systematic bias, we performed a calibration diagnosis on the internal validation set. The calibration parameters, including the slope (a), intercept (b), residual mean, residual standard deviation, MAE improvement after calibration, and ΔR², for each model are summarized in Table 5.

Table 5. Calibration judgment results of each model.

For the SVM model, the calibration slope was 0.9333 and the intercept was 0.0958, with a residual mean of 0.0028 and a residual standard deviation of 0.1633. For the GradientBoosting model, the slope was 0.9832 and the intercept was 0.0294, with a residual mean of 0.0059 and a residual standard deviation of 0.1654. For the XGBoost model, the slope was 1.0473 and the intercept was −0.0499, with a residual mean of 0.0155 and a residual standard deviation of 0.1660. For the simple average ensemble, the calibration slope was 1.0038 and the intercept was 0.0029, with a residual mean of 0.0081 and a residual standard deviation of 0.1624.

The residuals for all four models were symmetrically distributed around zero with approximately normal distributions, as confirmed by the residual-vs.-predicted scatter plots, histograms, true-vs.-predicted scatter plots, and Q-Q plots (Figure 2 for the ensemble; Supplementary Figures S1–S3 for the three individual models). The true-vs.-predicted scatter plots showed points closely aligned with the 45° diagonal line, indicating no substantial systematic over- or under-estimation.

Consistent with these findings, post-hoc linear calibration resulted in only marginal improvements in MAE across all models: 0.0004 for SVM, 0.0006 for GradientBoosting, 0.0015 for XGBoost, and 0.0003 for the ensemble, all well below the pre-defined threshold of 0.005. The changes in R² after calibration were also negligible (+0.0032, +0.0007, +0.0049, and +0.0010, respectively). The negligible improvements confirmed that all four models were already well calibrated and did not require additional post-processing. The ensemble model, with its calibration slope closest to 1 and intercept closest to 0 among all models, demonstrated the best overall calibration and was therefore selected as the final model for clinical application.

3.6 Feature importance and model interpretability

To elucidate the contribution of each feature to the model predictions, we performed SHAP (SHapley Additive exPlanations) analysis for the SVM and GradientBoosting models (Figures 3,4).

For the SVM model, the SHAP beeswarm plot revealed that BSA was the most influential feature, followed by circulating blood volume, gender, standardized blood flow, and pre-dialysis weight. Higher BSA values were associated with negative SHAP values (i.e., lower Kt/V predictions), consistent with the clinical expectation that larger body size requires higher dialysis doses to achieve adequate clearance. For the GradientBoosting model, the SHAP summary bar plot showed a similar ranking: BSA had the highest mean SHAP value (0.095), followed by standardized blood flow (0.077), gender (0.042), circulating blood volume (0.017), and pre-dialysis weight (0.013).

These consistent rankings across fundamentally different model architectures (geometric SVM vs. tree based GradientBoosting) reinforce the robustness of the feature importance findings.

3.7 Probabilistic prediction performance

To quantify prediction uncertainty, we evaluated two probabilistic prediction methods: LightGBM quantile regression and NGBoost. Performance metrics of probabilistic prediction on internal and external validation sets are presented in Table 6.

Table 6. Performance metrics of probabilistic prediction on internal and external validation sets.

Validation set refers to internal validation set, and Test set refers to external validation set. CRPS and NLL are not applicable to quantile regression models.

On internal validation, for the 80% prediction interval, NGBoost achieved a PICP (Prediction Interval Coverage Probability) of 82.27%, which closely approximates the nominal 80% coverage, with an MPIW (Mean Prediction Interval Width) of 0.4924 and a CWC (Coverage Width Criterion) of 0.4924. The CRPS (Continuous Ranked Probability Score) was 0.1240, and the NLL (Negative Log Likelihood) was 0.0084 on internal validation.

On internal validation, LightGBM quantile regression achieved a PICP of 86.82% (overly conservative), with an MPIW of 0.6123 and a CWC of 0.6123. While LightGBM covered more true values on internal validation, this came at the cost of substantially wider prediction intervals (approximately 24% wider than NGBoost).

Based on its superior calibration (PICP closer to nominal 80%), narrower prediction intervals, and ability to provide complete probability density information on internal validation, NGBoost was selected as the preferred method for clinical uncertainty quantification.

3.8 Clinical application interface

To facilitate clinical translation, we packaged the final ensemble model and the NGBoost probabilistic predictor into a local standalone executable application. This application is distributed as a self-contained file package that can be copied directly to a local computer and executed offline, without requiring any network connection.

The interface accepts five input parameters (sex, height, pre-dialysis weight, blood flow rate, and dialysis duration) and returns: (1) Individual model point predictions (SVM, GradientBoosting and XGBoost); (2) Ensemble average prediction; (3) NGBoost derived median and 80% prediction interval; (4) Probability density curve visualization.

A screenshot of the application interface is provided in Figure 5. The complete executable package and source code are available from the corresponding author upon reasonable request.

4 Discussion

4.1 Study overview

This multicenter study developed and validated predictive models for hemodialysis adequacy (Kt/V) using multiple linear regression, seven machine learning algorithms for point prediction, an ensemble model, and two probabilistic prediction methods. The key findings are: (1) the SVM model consistently outperformed other individual models across both internal and external validation sets, achieving an R² of 0.5931 and RMSE of 0.1633 on external validation; (2) the simple average ensemble of SVM, GradientBoosting, and XGBoost achieved the highest overall performance (R² = 0.5911, RMSE = 0.1637); (3) NGBoost provided well-calibrated prediction intervals (PICP of 82.27% for nominal 80% coverage) with narrower intervals than quantile regression; and (4) BSA, standardized blood flow, and gender were consistently identified as the most influential predictors across different modeling frameworks. This study was strictly designed, implemented and reported according to the TRIPOD + AI guidelines for machine learning-based prediction models (14).

4.2 Methodological contributions and innovations

From a methodological perspective, this study makes several contributions to the field of clinical prediction modeling for continuous outcomes. First, we employed a multi-algorithm fusion strategy for feature selection, integrating eight fundamentally different methods (F test, Mutual Information, Random Forest importance, XGBoost importance, LASSO coefficients, Permutation Importance, Recursive Feature Elimination, and Boruta) to generate a consensus importance ranking. The Boruta algorithm, a random forest-based wrapper method that iteratively compares feature importance with randomly permuted “shadow” features, has been shown to offer high stability and effectiveness in high-dimensional clinical datasets. By averaging rankings across multiple algorithms, our approach reduces the risk of selecting features that are important only under a single methodological paradigm. This ensemble-based selection strategy represents a more robust alternative to the conventional stepwise regression approach, which has been criticized for its instability and tendency to capitalize on chance.

Second, we addressed the critical issue of prediction uncertainty quantification, a dimension often neglected in clinical prediction model studies. For continuous outcomes such as Kt/V, providing only a point estimate without corresponding uncertainty information limits the clinical utility of the model, as clinicians cannot gauge the reliability of individual predictions. We implemented NGBoost (Natural Gradient Boosting) (15), an algorithm that generalizes gradient boosting to probabilistic regression by treating the parameters of the conditional distribution as targets for a multiparameter boosting algorithm. Unlike conventional regression models that return only a point estimate, NGBoost outputs a full probability distribution over the outcome space, enabling predictive uncertainty estimation. The Natural Gradient correction ensures stable training dynamics for the multiparameter approach. In our validation, NGBoost achieved a PICP of 82.27% for the nominal 80% prediction interval, with an MPIW of 0.4924, substantially narrower than the 0.6123 width of LightGBM quantile regression. This superior calibration (narrowness trade-off) suggests that NGBoost provides clinically meaningful uncertainty estimates without being overly conservative. These findings are consistent with previous reports demonstrating NGBoost's competitive performance in probabilistic prediction tasks across various domains.

Third, we enhanced model interpretability using SHAP (SHapley Additive exPlanations) analysis. SHAP, a model agnostic framework based on coalition game theory (16), quantifies the marginal contribution of each feature to individual predictions. The consistent ranking of BSA as the most influential predictor across both the SVM (geometric) and GradientBoosting (tree-based) models reinforces the robustness of this finding. SHAP analysis has been increasingly recognized as a valuable tool for translating “black box” machine learning models into clinically actionable insights, facilitating clinician trust and fair model evaluation. By providing both global feature importance rankings and local explanations for individual predictions, our approach addresses the interpretability requirements emphasized in the TRIPOD + AI guidelines.

4.3 Clinical translation and practical utility

A key contribution of this study is the explicit bridging of post-hoc model development with pre-dialysis clinical application. Although the model was developed using actual post-dialysis parameters (the only objectively documented values available in the retrospective dataset), its intended clinical use is pre-dialysis prediction. In routine practice, dialysis prescriptions, including planned duration, blood flow rate, and ultrafiltration volume, are determined by clinicians before the session begins. The high concordance between prescribed and delivered parameters (mean absolute difference <2% in our dataset) ensures that planned values can be substituted into the model to obtain reliable predictions of expected Kt/V. This enables real-time, pre-emptive adjustment of dialysis prescriptions, allowing clinicians to modify parameters before the session starts if the predicted Kt/V falls below target.

This functionality fundamentally distinguishes our approach from the Daugirdas formula, which requires post-dialysis urea measurements and therefore can only provide retrospective adequacy assessment. Our model offers prospective decision support, potentially reducing the proportion of sessions with inadequate dialysis and its associated complications.

The features selected by our consensus ranking strategy, BSA, pre-dialysis weight, standardized blood flow, sex, and circulating blood volume, are all non-invasive and readily available at the point of care. This is consistent with recent findings that machine learning models based on non-invasively collected general and dialysis-related variables can maintain high prediction accuracy. The identification of vascular access, gender, BMI, ultrafiltration volume, and dialysis duration as significant predictors in previous studies further supports the clinical plausibility of our selected features.

4.4 Comparison with previous studies

Previous attempts to predict Kt/V using machine learning have been limited by single center designs, lack of external validation, and insufficient methodological detail. Alexander et al. (2003) used MLR, tree-based modeling, and artificial neural networks but achieved a correlation coefficient of only 0.745 and lacked internal and external validation (17). Barbieri C et al. (2019) developed a multi-endpoint ANN model but relied on a single center and did not describe predictive performance or conduct feature importance analysis (18). More recent studies have demonstrated the feasibility of machine learning for dialysis adequacy prediction using demographic, intradialytic, and laboratory variables (19). He et al. proposed a domain-adversarial neural network for Kt/V prediction (20), achieving high accuracy on data from 17 patients. However, these studies predominantly focused on point prediction without addressing uncertainty quantification or model interpretability.

Our study addresses these gaps through: (1) a large multicenter sample (1,017 patients from seven centers); (2) comprehensive internal and external validation; (3) probabilistic prediction with well-calibrated uncertainty intervals; (4) multi algorithm feature selection with consensus ranking; (5) SHAP-based interpretability; and (6) a deployable local standalone application for clinical use.

5 Limitations

Several limitations should be acknowledged. First, the training set size (n = 677) and external validation set (n = 120) are moderate. Although our sample size exceeds the “10 events per predictor” rule of thumb for the five selected features (21), future studies with larger and more diverse populations are needed to confirm generalizability. Second, the study population was drawn from a single geographic region (Inner Mongolia, China), which may limit applicability to other ethnic groups and healthcare settings. Third, while we demonstrated excellent performance on internal and external validation, prospective deployment studies are required to evaluate real-world clinical impact. Fourth, the NGBoost probabilistic predictions assume a log-normal distribution for log transformed Kt/V; although our residual diagnostics supported this assumption, future work could explore distribution-free uncertainty quantification methods such as conformal prediction. Fifth, the simple average ensemble, while effective, could potentially be improved through more sophisticated stacking or Bayesian model averaging approaches.

6 Recommendations

Future research should focus on: (1) prospective validation of the model in routine clinical practice, ideally with randomized controlled trial designs to evaluate impact on patient outcomes; (2) expansion to diverse geographic and ethnic populations; (3) integration of the prediction tool into electronic health record systems for seamless clinical workflow integration; (4) exploration of deep learning approaches that can capture complex temporal patterns in dialysis data; and (5) investigation of model fairness across different demographic subgroups, as recommended by the TRIPOD + AI framework.

7 Conclusion

This study demonstrates that machine learning models, particularly SVM, ensemble averaging, and NGBoost, can provide accurate, interpretable, and uncertainty-aware predictions of hemodialysis adequacy (Kt/V) using non-invasive, readily available parameters. The methodological innovations, including multi-algorithm feature selection, probabilistic prediction with NGBoost, and SHAP-based interpretability, collectively address the key limitations identified in previous studies and align with the updated TRIPOD + AI reporting standards. The development of a local standalone executable application further facilitates clinical translation, enabling pre-dialysis prediction and real-time prescription adjustment. If validated prospectively, this approach has the potential to improve dialysis quality, reduce complications, and support personalized treatment planning for patients receiving maintenance hemodialysis.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fdgth.2026.1771806/full#supplementary-material

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Footnotes

Edited by: Jiawen Deng, University of Toronto, Canada

Reviewed by: Solikhun Solikhun, STIKOM Tunas Bangsa, Indonesia

Myroslav Havryliuk, Lviv Polytechnic National University, Ukraine

Abbreviations ANN, artificial neural network; BPNN, back propagation neural network; CKD, chronic kidney disease; CRPS, continuous ranked probability score; CWC, coverage width criterion; DOQI, dialysis outcome quality initiative; ESRD, end stage renal disease; GB, gradient boosting; HD, hemodialysis; KNN, K nearest neighbor; Kt/V, urea clearance index; LASSO, least absolute shrinkage and selection operator; LightGBM, light gradient boosting machine; MAE, mean absolute error; MAPE, mean absolute percentage error; MHD, maintenance hemodialysis; MI, mutual information; MLR, multiple linear regression; MPIW, mean prediction interval width; MSE, mean squared error; NGBoost, natural gradient boosting; NLL, negative log likelihood; PICP, prediction interval coverage probability; PTH, parathyroid hormone; RF, random forest; RFE, recursive feature elimination; RMSE, root mean squared error; SHAP, SHapley additive exPlanations; spKt/V, single pool Kt/V; SVM, support vector machine; URR, urea reduction rate; XGBoost, extreme gradient boosting.