Prediction of Lithium Therapeutic Target Attainment in Adolescents Using an XGBoost Model: A Single-Center Retrospective Study

Introduction

Lithium carbonate remains a cornerstone pharmacotherapy for the maintenance treatment of severe mood dysregulation in adolescents.¹ Importantly, affective disorders are increasingly recognized not as entirely discrete diagnostic categories, but as a spectrum with shared biological underpinnings.² Although lithium is widely used in the treatment of these disorders, its clinical management is complicated by a narrow therapeutic index and substantial inter-individual pharmacokinetic variability.³ Lithium carbonate is uniquely characterized by negligible plasma protein binding, lack of hepatic metabolism, and almost exclusive renal elimination unchanged via glomerular filtration and proximal tubular reabsorption. Following oral administration, it is rapidly and completely absorbed, with peak serum concentrations generally reached within 1 to 4 hours.⁴ Within the therapeutic range, lithium exhibits linear pharmacokinetics, and its apparent volume of distribution averages 0.7–0.9 L/kg. Crucially, the elimination half-life is highly variable; while it averages 18–24 hours in adults, it is often shorter in pediatric and adolescent populations due to age-related differences in renal clearance capacity.⁵ This distinct pharmacokinetic profile renders it highly susceptible to rapid physiological changes in adolescents.⁶

According to the Association for Neuropsychopharmacology and Pharmacopsychiatry (AGNP) Consensus Guidelines, the recommended therapeutic reference range is 0.5–1.2 mmol/L, which is often stratified into a maintenance target of 0.5–0.8 mmol/L and an acute treatment limit up to 1.2 mmol/L.⁷ Within this range, 0.5 mmol/L is generally regarded as the critical lower threshold for efficacy; concentrations below this level are associated with an increased risk of relapse and insufficient therapeutic response.⁸ In clinical practice, a key question is whether the patient has attained this minimum effective threshold. Identifying patients at risk of subtherapeutic exposure (<0.5 mmol/L) may therefore facilitate earlier recognition of inadequate exposure and support timely dose adjustment.

Adolescents should not be regarded simply as “smaller adults”, as differences in body size, body composition, and organ maturation result in pharmacokinetic characteristics that differ significantly from those of adults.⁹ In general, pediatric patients undergo rapid changes in their organs during maturation, which can lead to differences in physiological conditions and drug disposition.¹⁰ Therefore, these developmental variations must be taken into account when evaluating lithium exposure and therapeutic target attainment specifically in adolescent patients.

Currently, therapeutic drug monitoring (TDM) is the standard of care for managing lithium therapy.¹¹ However, TDM is inherently reactive, often requiring patients to reach steady-state before dose adjustments can be made, which may expose them to periods of subtherapeutic or supratherapeutic levels. While population pharmacokinetic (PopPK) modeling has been employed to guide initial dosing, these models typically rely on predefined compartmental assumptions.¹² They often struggle to fully capture the complex, non-linear interactions between heterogeneous clinical covariates, particularly in pediatric and adolescent cohorts where data is frequently sparse.

To address these challenges, machine learning (ML) offers a valuable complementary approach. Studies based on real-world evidence extracted from medical records often involve large amounts of clinical data and therefore require appropriate methods for data processing and analysis. As a data-driven method, ML can integrate multidimensional information, identify key influencing factors, and build prediction models.^13–15 Recent studies based on real-world evidence have developed ML models to predict blood concentrations for neuropsychiatric drugs (eg., levetiracetam and quetiapine), thereby demonstrating potential in providing insights to optimize dosage regimens.^16,17 However, studies applying ML to predict lithium therapeutic target attainment in adolescents remain scarce.

Therefore, this study aimed to identify factors associated with lithium therapeutic target attainment in adolescents through machine-learning methods and to develop and internally validate a prediction model to assess whether lithium concentrations reach the therapeutic target range.

Methods

Patients and Data Collection

This single-center, retrospective study was conducted at Xi’an Mental Health Center and approved by the institutional Ethics Committee. Adolescent inpatients who received lithium carbonate therapy at Xi’an Mental Health Center between January 1, 2024, and May 31, 2025, were systematically identified through the hospital electronic medical record system. Patients were included if they met the following criteria: (1) age < 18 years; (2) received a constant daily dose of lithium for >7 consecutive days under nursing supervision; and (3) availability of at least one steady-state TDM measurement. Patients were excluded if they had lithium plasma concentrations below the lower limit of quantification (LLOQ).

A dataset, including lithium administration details, demographic information, laboratory parameters, and concomitant therapies, was retrieved from the hospital’s Electronic Medical Record System. Under steady-state conditions, blood samples were drawn during the 30-minute window preceding the scheduled morning administration. Serum lithium concentrations were measured as part of routine clinical care at Xi’an Mental Health Center using standard validated methods.

Data Preparation

The primary outcome was defined as attainment of the minimum therapeutic target concentration, binary-coded as 1 if the serum lithium concentration was ≥0.5 mmol/L and 0 otherwise.

The overall workflow of this study is illustrated in Figure 1. Initially, variables with a missing data rate exceeding 20% were excluded. The entire dataset was divided into three subsets at the patient level to prevent data leakage caused by multiple observations from the same individual. Specifically, the dataset was partitioned into a Training Set (56%), an Internal Validation Set (24%), and a Held-out Testing Set (20%).

Figure 1 The workflow of model development and validation. Abbreviations: GBM, Gradient Boosting Machine; SVM, Support Vector Machine; LR, Logistic Regression; RF, Random Forest; Ada, AdaBoost; KNN, K-Nearest Neighbors; GNB, Gaussian Naïve Bayes; XGB, XGBoost; LGBM, LightGBM; Cat, CatBoost.

The partitioning process followed a two-stage stratified strategy: First, 20% of the unique patients were randomly held out as the Held-out Testing Set for the final performance evaluation. The remaining 80% of patients constituted the development set, which was further split into a Training Set (70% of development data) and an Internal Validation Set (30% of development data). Crucially, we employed a Stratified Group K-Fold cross-validation strategy (using the StratifiedGroupKFold algorithm). This approach ensures that all observations from the same patient are assigned to the same set, thereby providing an unbiased evaluation of the model’s generalization ability to new subjects. Missing values in the remaining variables were imputed using the median value of the respective feature. Feature scaling using z-score normalization (StandardScaler) was applied only to algorithms sensitive to feature scale, whereas tree-based and boosting-based models were trained using the original unscaled features.

Feature Selection

To identify relevant clinical predictors and reduce data dimensionality, we employed the Boruta algorithm, an “all-relevant” feature selection wrapper method built around the Random Forest classifier.¹⁸ Unlike traditional methods that seek a minimal optimal feature subset, Boruta aims to identify all features relevant to the outcome variable. It achieves this by creating randomized “shadow” features and comparing the importance (Z-score) of real features against the maximum Z-score among these shadow features in iterative Random Forest runs. To ensure reproducibility and prevent overfitting, the algorithm was configured with a maximum of 100 iterations, and the underlying Random Forest estimator was limited to a maximum depth of 5. To mitigate potential bias during the feature selection phase, the Boruta algorithm was configured with balanced class weights. Features that consistently outperformed the best shadow feature were deemed “confirmed” and selected. To validate the efficacy of this feature selection strategy, we conducted a comparative analysis by developing baseline models using both the full spectrum of clinical variables (Full Feature Set) and the Boruta-identified subset (Selected Feature Set). This rigorously filtered feature subset was used for all subsequent model development steps to enhance model stability and interpretability.¹⁹ Specifically, Boruta feature selection was performed once on the training set only, and the selected feature subset was then carried forward into all subsequent modeling procedures.

Model Development and Optimization

Ten supervised machine learning algorithms were initially screened: Gradient Boosting Decision Tree (GBDT), Support Vector Machine (SVM), Logistic Regression (LR), CatBoost, Random Forest (RF), XGBoost, LightGBM, AdaBoost, K-Nearest Neighbors (KNN), and Gaussian Naïve Bayes (GNB).^20,21 Hyperparameter optimization was conducted for each algorithm on the training Set using grid search with 10-fold cross-validation, utilizing the Area Under the Receiver Operating Characteristic Curve (AUC) as the optimization metric.²² The top three performing algorithms based on their AUC scores in the internal validation set were selected as candidate base learners. To address the issue of class imbalance and achieve an optimal balance between sensitivity and specificity, the optimal classification threshold for the top-performing base learners and ensemble models was determined using the Youden Index (J = Sensitivity + Specificity − 1) on the Internal Validation Set, rather than using a default threshold of 0.5.²³

Ensemble Strategy Construction

To investigate whether combining base learners could further enhance predictive performance, two distinct ensemble strategies were developed and evaluated. First, an optimized weighted-average ensemble was constructed. The integration weights for each base learner were determined using the Sequential Least Squares Programming (SLSQP) optimization algorithm. The objective function was defined to minimize the Logarithmic Loss (LogLoss) on the internal validation set, ensuring calibrated probability outputs and optimal weight allocation.¹³ Second, a stacking ensemble model was developed.²⁴ Using a 10-fold cross-validation scheme to prevent data leakage, the base learners generated probability predictions (meta-features), which were recalibrated using Platt scaling (sigmoid calibration) and then utilized by a Logistic Regression meta-learner to produce the final calibrated prediction.²⁵ Logistic Regression was explicitly selected for this role due to its interpretability and effectiveness in probability calibration, offering a robust linear combination of base model outputs while mitigating the risk of meta-level overfitting.

Model Evaluation and Interpretation

Prior to the final evaluation, the identified top-performing base learners and the ensemble models were retrained on the entire development dataset (combining both the Training and Internal Validation Sets) using the optimized hyperparameters.

A comparative evaluation was then conducted on the held-out Testing Set to determine the final model for deployment. Model performance was comprehensively assessed using discrimination metrics including AUC, accuracy, precision, recall, and F1-score. Calibration curves based on patient-level out-of-fold predictions and DCA on the independent testing set were used to determine the goodness-of-fit and clinical utility of the model, respectively. Additionally, to address the “black-box” nature of the model, we employed SHapley Additive exPlanations (SHAP) to provide transparent model interpretation. SHAP summary plots, dependence plots, and force plots were generated to elucidate the global importance of features and to explain individual predictions, ensuring plausibility.²⁶

Sensitivity Analysis

To examine whether the main findings were sensitive to the choice of therapeutic threshold, we conducted an additional sensitivity analysis using an alternative outcome definition. Specifically, the binary endpoint was redefined as serum lithium concentration >0.8 mmol/L, corresponding to the upper therapeutic reference range for maintenance treatment. Importantly, the identical patient-level data partition, the predictor set identified by Boruta, and the same modeling strategy were retained. No additional feature selection, hyperparameter tuning, or model redevelopment was performed. Model performance was assessed using the same discrimination and calibration metrics as in the primary analysis.

Development of an Exploratory Web-Based Interface

To demonstrate the proof-of-concept and facilitate future external validation, a lightweight web application was developed. The system was built using the Flask web framework (v2.3.3) within a Python 3.12.2 environment. To ensure strict consistency between offline model development and online inference, the deployment environment enforced identical library dependencies, specifically scikit-learn (1.5.2), xgboost (2.1.3), and joblib (1.4.2).

Powered by the final optimized model, the web-based interface accepts the identified key clinical features to generate real-time probability estimates of therapeutic target attainment, serving as an exploratory platform for future validation and research use. Furthermore, to enhance model transparency, a dynamic SHAP waterfall plot module was integrated into the interface to visually break down the specific contribution of each input variable to the final prediction.

Statistical Analysis

Continuous variables were expressed as median [interquartile range (IQR)] and compared using the Mann–Whitney U-test. Categorical variables were described as frequencies (percentages) [n (%)] and analyzed using the Chi-square test. All data processing, modeling, and statistical analyses were performed using Python (version 3.12.2). To account for intra-patient correlation arising from multiple observations per patient, 95% confidence intervals (CIs) for the reported performance metrics were estimated using patient-level cluster bootstrapping with 2,000 resamples.

Results

Baseline Characteristics

A total of 279 patients (406 observations) were included in the study. The primary diagnoses included bipolar disorder, manic episodes, and behavioral and emotional disorders with childhood or adolescent onset. Specifically, 186 patients had a single observation, while 67, 18, and 8 patients had two, three, and four observations, respectively. The dataset was divided into a development set (332 observations) and a held-out testing set (74 observations), with all observations from the same patient assigned to the same subset. The development set was further split into the training set (235 observations) and internal validation set (97 observations). Overall, 324 of the 406 observations (79.8%) achieved the therapeutic target concentration, including 192/235 (81.7%) in the training set, 75/97 (77.3%) in the internal validation set, 267/332 (80.4%) in the development set, and 57/74 (77.0%) in the held-out testing set. Tables 1 and 2 summarize the baseline characteristics of the 49 included variables. The training and testing cohorts were generally well-balanced; all variables showed no statistically significant differences (p ≥ 0.05).

Table 1 Baseline Characteristics of Patients in the Development and Held-Out Testing Sets

Table 2 Baseline Characteristics of the Internal Training and Internal Validation Sets

Feature Selection and Base Learner Screening

To identify the most informative clinical predictors, the Boruta algorithm was applied to the training Set. This rigorous selection process identified 5 important features from the original clinical variables: weight, height, daily dose, albumin (ALB) and alkaline phosphatase (ALP).

To validate the effectiveness of this feature subset and assess potential information loss, a comparative analysis was conducted. Ten candidate supervised machine learning algorithms were trained and optimized on the training Set using both the full feature set and the selected 5-feature subset. Hyperparameters for all models were tuned via 10-fold stratified group cross-validation. The comprehensive hyperparameter search spaces and the final optimal values for all evaluated algorithms are detailed in Supplementary Table S1. As presented in Table 3, the models utilizing the simplified feature subset achieved predictive performance comparable to, and in some cases superior to, those using the full feature set. This confirmed that the selected features captured the majority of the predictive information while significantly reducing model complexity.

Table 3 Performance Comparison of Ten Machine Learning Algorithms Using Full vs. Selected Features on the Internal Validation Set (n = 97 Observations from 67 Patients)

Consequently, the models optimized on the selected feature subset were evaluated on the held-out internal Validation Set to determine the final base learners. Based on the AUC scores in the internal validation set, the top three performing models were identified as candidates for the final ensemble strategy: Logistic Regression (AUC: 0.812, 95% CI: 0.690–0.916), AdaBoost (AUC: 0.794, 95% CI: 0.663–0.894), and XGBoost (AUC: 0.792, 95% CI: 0.674–0.897).

Ensemble Model Development and Final Evaluation

The generalization capability of the developed models, including the top base learners and two distinct ensemble frameworks (Weighted and Stacking), was comprehensively evaluated on the held-out Testing Set. As summarized in Table 4, the comparative analysis revealed that the single XGBoost model offered optimal predictive performance and interpretability compared to the complex ensemble strategies. While the Stacking ensemble achieved an AUC of 0.846 (95% CI: 0.723–0.942), comparable to XGBoost, it exhibited a conservative prediction tendency, resulting in a lower accuracy of 0.676 (95% CI: 0.563–0.790) and a recall of 0.614 (95% CI: 0.478–0.754). In contrast, the XGBoost model demonstrated a better balance of performance, achieving a higher accuracy of 0.770 (95% CI: 0.671–0.861) and a recall of 0.754 (95% CI: 0.630–0.865) while maintaining an identical AUC of 0.846 (95% CI: 0.742–0.928). This indicates that the XGBoost model effectively minimized false negatives without the computational complexity or sensitivity loss observed in the Stacking framework. Consequently, XGBoost was identified as the optimal performer in this study.

Table 4 Predictive Performance of Selected Base Learners and Ensemble Models on the Held-Out Testing Set (n = 74 Observations from 56 Patients)

Clinical Utility and Model Interpretation

The Brier score was 0.128. Decision curve analysis performed in the independent test set showed that the model provided a higher net benefit than both treat-all and treat-none strategies across a wide range of clinically relevant threshold probabilities (Figure 2A). Calibration analysis based on patient-level out-of-fold predicted probabilities demonstrated acceptable agreement between predicted and observed risks (Figure 2B).

Figure 2 Model calibration and clinical utility evaluation of the final XGBoost model. (A) DCA of the XGBoost model on the independent testing set. The x-axis indicates the threshold probability, and the y-axis indicates the clinical net benefit. The red solid line represents the net benefit of the model, the dashed grey line represents the assumption that all patients receive intervention (Treat All), and the dotted black line represents the assumption that no patients receive intervention (Treat None). (B) Calibration curve of the XGBoost model. The x-axis represents the predicted probability of target attainment, and the y-axis represents the observed actual probability. The diagonal dashed grey line indicates perfect prediction (ideal calibration), while the red solid line represents the apparent predictive performance of the model.

The interpretability of the model is illustrated in Figure 3. Specifically, Figure 3A maps the SHAP values (x-axis) corresponding to individual feature values, where their summation determines the likelihood of achieving the lithium therapeutic threshold. As detailed in Figure 3B, features are ranked by their mean absolute SHAP values, identifying Daily dose and Weight as the dominant predictors.

Figure 3 SHAP-based feature importance analysis of the final XGBoost model. (A) SHAP summary plot. Each dot represents one observation. Color indicates the relative feature value (red, high; blue, low), and the x-axis shows the SHAP value, representing the impact of each feature on the model output. Positive SHAP values indicate a higher likelihood of therapeutic target attainment. (B) Mean absolute SHAP values ranking the importance of five variables. Abbreviations: SHAP, Shapley Additive Explanations; ALB, albumin; ALP, alkaline phosphatase.

Furthermore, the SHAP dependence plots (Figure 4) elucidate the direction of influence. Daily dose exhibited a clear positive relationship with the predicted probability of target attainment, confirming dose as the most direct driver of increased exposure. In contrast, physiological parameters such as Body weight, ALB, Height, and ALP generally exhibited negative correlations. Specifically, as body weight increased from 40 kg to 80 kg, the SHAP value progressively declined.

Figure 4 SHAP dependence plots for the five most influential predictors in the final XGBoost model. (A) Daily dose, (B) weight, (C) albumin (ALB), (D) height, and (E) alkaline phosphatase (ALP). Each dot represents one observation. The x-axis represents the actual value of each feature, while the y-axis represents the SHAP value, where positive values indicate a higher likelihood of therapeutic target attainment.

To demonstrate the model’s decision-making process at the individual patient level, we visualized specific prediction instances using SHAP force plots (Figure 5). These plots illustrate how each feature contributes to pushing the model output from the base value to the final prediction. Figure 5A exemplifies a case with a low predicted probability (non-attainment). Here, negative contributors (blue bars) dominated the prediction; specifically, a lower Daily dose (0.6 g) combined with a higher Weight (79.0 kg) and elevated ALB suppressed the output value to −1.52. Conversely, Figure 5B depicts a case with a high predicted probability (target attainment). In this scenario, the prediction was strongly driven by positive contributors (red bars); a higher Daily dose (0.9 g) and lower Weight (53.0 kg) significantly increased the model output to a score of 5.32.

Figure 5 Visualization of individual prediction explanations using SHAP force plots. The SHAP force plot illustrates how specific clinical features contribute to the model’s prediction relative to the base value (the mean prediction of the dataset). Features pushing the prediction higher are shown in red, whereas those pushing it lower are shown in blue. The bold number above each plot indicates the final model output, f(x), for the corresponding individual. (A) A representative low-exposure risk case, primarily driven by negative contributions. (B) A representative target attainment case, primarily driven by positive contributions. Abbreviations: ALB, albumin; ALP, alkaline phosphatase; SHAP, Shapley Additive Explanations.

Sensitivity Analysis for a Higher Therapeutic Threshold

Among the 406 steady-state TDM observations, 78 (19.2%) and 6 (1.48%) exceeded the 0.8 mmol/L and 1.2 mmol/L thresholds, respectively. To further assess the model under a higher maintenance threshold, a sensitivity analysis using an alternative outcome definition (>0.8 mmol/L) was performed. Under this alternative threshold, the model maintained good discrimination in the held-out testing set (AUC = 0.858; Supplementary Figure S1). Calibration also remained acceptable (Brier score = 0.135; Supplementary Figure S2). SHAP analysis showed that daily dose remained the dominant contributor, while the relative importance of the remaining predictors showed modest reordering under the alternative threshold (Supplementary Figure S3).

Web-Based Deployment and Verification

The final XGBoost classifier was successfully deployed as an exploratory online tool (Figure 6). For intuitive interpretation, the interface outputs an exploratory binary risk label (“Adequate/High Exposure” vs. “Low Exposure Risk”) for risk stratification purposes. Taking a simulated patient as an example (Figure 6), the tool generates a “Low exposure risk” label. Concurrently, the SHAP waterfall plot elucidates this outcome, revealing that daily dose and ALB negatively drive the target attainment probability, while weight provides a positive contribution. The tool is hosted on a cloud application platform and is publicly accessible at https://blood-concentration-predictor.onrender.com/.

Figure 6 Interface of the exploratory web-based tool for predicting lithium therapeutic target attainment. Patient-specific covariates are entered into the input fields, and the model generates an exploratory risk label (“Adequate/High Exposure” or “Low Exposure Risk”). The waterfall panel visualizes the contribution of individual predictors to the model output, where positive SHAP values (green) increase the predicted probability of target attainment and negative SHAP values (red) decrease it.

Discussion

In this study, we successfully established an exploratory, interpretable machine learning model to predict lithium therapeutic target attainment in psychiatric inpatient adolescents. After a comprehensive evaluation of multiple algorithms, we identified XGBoost as the optimal classifier. Unlike traditional linear methods, the XGBoost framework effectively captured complex, non-linear interactions within the clinical data. The XGBoost algorithm demonstrated acceptable internal performance for this specific population, serving as an exploratory tool for sub-therapeutic risk stratification.

Traditionally, lithium dosing relies on a “reactive” paradigm, prescribing empirically and adjusting based on subsequent TDM. This exploratory model may provide a complementary approach for estimating lithium therapeutic target attainment, particularly in primary care institutions or resource-limited settings where TDM facilities are often unavailable. Unlike mechanism-based models that necessitate complex intermediate PK parameter estimation, it does not require prior estimation of pharmacokinetic parameters and was implemented in a web-based interface for exploratory use. Decision curve analysis suggested potential net benefit across a range of threshold probabilities in the testing set. In summary, this tool serves as an exploratory platform for estimating the probability of lithium therapeutic target attainment.

The SHAP analysis provides a granular interpretation of the model by quantifying both the importance hierarchy and the direction of impact for each feature. These results demonstrate that therapeutic target attainment can be effectively estimated using a combination of routinely available clinical indicators.

Daily dose was identified as the foremost predictor of therapeutic target attainment. This finding is underpinned by the established pharmacokinetic profile of lithium carbonate, which follows linear, first‑order elimination kinetics within the therapeutic range.²⁷ Unlike drugs with saturable metabolism (eg., valproic acid or phenytoin), lithium is not metabolized by the liver and exhibits low protein binding.¹ Consequently, its steady‑state serum concentration is directly proportional to the administered dose, making dosage the most immediate modifiable factor influencing target attainment.

Next, body weight emerged as a critical predictor. Since lithium is a water-soluble ion, its volume of distribution is closely tied to total body water.²⁸ The inverse relationship between weight and serum concentration is well-corroborated by previous PopPK studies. Research focusing on pediatric patients demonstrated that incorporating weight into allometric scaling functions for systemic clearance and volume of distribution significantly improved model fit, their simulation results explicitly indicated that increases in body weight were associated with downward trends in Cmax, AUC, and steady-state trough concentrations.²⁹ Extending this to the adolescent population, another study focusing on 8–18-year-olds with bipolar disorder concluded that once body size is accounted for, lithium pharmacokinetic parameters in adolescents are comparable to those in adults. Notably, this study advocated for weight-based dosing regimens (eg., 25 mg/kg total body weight/day) to optimize therapeutic outcomes.¹⁰ These findings align with our results and underscore the importance of weight monitoring in the adolescent population. Given that adolescence is characterized by rapid growth spurts and heterogeneous changes in body composition, weight serves as a more sensitive metric than age for guiding dosage.

Considering that lithium is primarily cleared by the kidneys, renal-related variables, including creatinine, estimated glomerular filtration rate (eGFR), sodium, and urea, were included in the original candidate feature set. Among these renal-related variables, creatinine and eGFR were each entered separately for feature selection. However, none of these markers was retained in the final model, which is consistent with previous studies reporting similar findings.^30,31 This may be because renal function was generally normal and relatively homogeneous in our adolescent cohort. As these renal indices showed a relatively compact distribution, they may have played a smaller role in the final model than other factors.

ALB and ALP were identified as influential predictors in the model. Considering the negligible protein binding of lithium, the association of ALB with model output is unlikely to reflect a direct mechanistic effect on lithium pharmacokinetics. Instead, ALB may be more appropriately regarded as a surrogate marker of broader physiological status, such as nutritional condition, hydration, concomitant illness, or overall physiological reserve.³² Similarly, in adolescents, ALP may be more reasonably interpreted as a proxy for developmental status, including growth, puberty, or body composition, rather than as a direct determinant of lithium pharmacokinetics.³³ Within this machine learning framework, these markers are better understood as model-derived statistical associations rather than definitive causal determinants. Their inclusion may therefore reflect broader physiological variability without implying a direct mechanistic link to lithium pharmacokinetics.

Finally, height was included in the final model.³⁴ Given the hydrophilic nature of lithium, its apparent volume of distribution scales directly with total body water and mass. Since weight was already identified as a primary predictor, it captures the majority of the variance related to body size. Consequently, the independent impact of height is naturally diminished. However, it likely complements weight to implicitly account for variations in body habitus. In the adolescent population, height serves as a “fine-tuning” variable, enhancing the model’s precision by resolving residual variance in body shape that weight alone cannot fully address.³⁵

While our findings align with Hsu et al (2021) in identifying daily dose as the primary determinant of lithium levels, significant divergence exists in other critical predictors. In the adult-centered model of Hsu et al, Age emerged as a key feature, likely because renal function in adults declines with age.³⁶ Their model also required factors such as systolic blood pressure and concomitant medications to achieve robust performance, reflecting that adult patients more commonly present with complex drug-drug interactions and substance-related metabolic interferences. In contrast, our adolescent model prioritized body size and specific laboratory parameters.

It is worth noting that we also explored more complex ensemble strategies, including Stacking and Voting, to potentially boost performance. However, our comparative analysis revealed that these “meta-models” did not outperform the standalone XGBoost algorithm; in fact, they exhibited slightly inferior performance in key metrics such as accuracy and F1-score. First, XGBoost, as a gradient boosting framework, already possesses a strong inherent capacity to model complex non-linearities and interactions. Given our dataset’s characteristics, the XGBoost model likely approached the performance ceiling achievable with the available predictors, rendering the marginal gains from further model stacking negligible.³⁷ Second, and perhaps more importantly, is the trade-off between complexity and generalization. Advanced ensemble techniques like Stacking introduce additional structural complexity and parameters. In the context of a modest sample size, this added complexity risks fitting the model to noise rather than the underlying signal (overfitting), as evidenced by the drop in accuracy for the Stacking model.³⁸ Consequently, we prioritized the XGBoost model. It strikes the optimal balance by delivering maximum predictive power while maintaining a streamlined architecture that is easier to interpret and implement.

A sensitivity analysis redefining the outcome as >0.8 mmol/L (the upper maintenance limit) further supported the robustness of the overall modeling framework. Model performance remained stable under the alternative threshold, suggesting that the main findings were not overly dependent on a single threshold definition.

Several limitations of this study should be acknowledged. First, the limited sample size from a single center may constrain the model’s generalizability, necessitating future external multi-center validation. Second, the retrospective design inevitably introduces potential selection bias and confounding factors. Third, the primary outcome was modeled as a binary endpoint centered on the minimum therapeutic threshold (0.5 mmol/L). However, this approach does not fully capture the information contained in continuous lithium concentrations. Future studies with larger datasets should explore regression-based or multi-threshold models for a more granular assessment of both efficacy and safety. Fourth, due to the low incidence of markedly supratherapeutic concentrations (>1.2 mmol/L) in our cohort, developing a dedicated toxicity-prediction model was unfeasible. Future datasets enriched for high-concentration cases are needed to enable true toxicity risk stratification. Finally, some patients contributed repeated observations. Although patient-wise data splitting and patient-level cluster bootstrapping were used to reduce potential bias from intra-patient correlation, the models were trained on observation-level data rather than within a formal repeated-measures modeling framework.

Conclusion

In conclusion, we developed an interpretable machine learning model for predicting lithium therapeutic target attainment in adolescents. Using routinely available clinical variables, this model provides an exploratory approach for estimating the probability of target attainment. Furthermore, we implemented a web-based interface as an exploratory platform for future validation and research use.