This retrospective cohort study included patients with prostate cancer treated at the Brunei Cancer Centre (TBCC), Jerudong Park Medical Centre (JPMC) between January 2018 and December 2024. The study was conducted in accordance with the TRIPOD-AI (Transparent Reporting of a multivariable prediction model for Individual Prognosis or Diagnosis adapted for Artificial Intelligence) reporting guidelines to enhance transparency, reproducibility, and clinical applicability (Multimedia Appendix 1). Eligible participants were male patients with a histologically confirmed diagnosis of prostate cancer who received at least 1 treatment modality at TBCC. Patients with secondary metastatic prostate cancer or multiple primary malignancies were excluded following clinical review by oncologists. In total, 212 eligible patients contributed 478 time-indexed observations, capturing longitudinal treatment pathways and disease progression, as illustrated in Figure 1. Demographic, clinical, and treatment data were extracted from multiple sources, including the Brunei Darussalam Healthcare Information and Management System, the JPMC laboratory database, and both electronic and manually documented clinical records.

Data were structured across 2 hierarchical levels: the individual case level and the longitudinal cohort level. At the case level, collected variables encompassed demographic characteristics (including age, ethnicity, and marital status), medical history (eg, comorbidities and family history of malignancy), lifestyle factors (particularly smoking status), anthropometric measures (height, weight, and BMI), tumor-related features (Gleason score, tumor stage, and ECOG performance status), and survival outcomes. At the cohort level, comprehensive treatment information and 23 longitudinal laboratory measures were assembled, including baseline and 1-year posttreatment PSA values, full blood count parameters, liver and renal function indicators, erythrocyte sedimentation rate, calcium, and phosphate levels. Absolute and relative changes in PSA were derived, and tumor staging, Gleason grade groups, and National Comprehensive Cancer Network prognostic risk categories were harmonized to address inconsistencies in coding practices [24]. Except for PSA, all laboratory parameters were baseline measurements obtained before each treatment modification, enabling clinicians to assess laboratory results when determining subsequent therapeutic strategies. Data integrity was ensured through systematic cross-validation across multiple data sources and clinical verification by oncologists before analysis. The primary outcome was PFS.

PFS was calculated on a treatment-episode basis from treatment initiation to the earliest occurrence of PSA-defined progression, which was determined using PCWG-consistent criteria based on treatment-specific baseline PSA, posttreatment PSA nadir, and confirmed PSA rises (≥25% and ≥2.0 ng/mL, with timing and confirmation rules applied) or death from any cause, with patients censored at last follow-up in the absence of an event. A PFS event value of 0 indicates no biochemical progression or survival without progression at the last follow-up (censored), whereas a value of 1 indicates the occurrence of biochemical PSA progression or death from any cause. To conduct a survival analysis, PFS events, which denote the occurrence of progression or death, were used rather than PFS, which indicates the time a patient remains progression-free.

Statistical analyses were performed using RStudio (version 4.3.3; R Foundation for Statistical Computing), Python (version 3.10; Python Software Foundation) implemented via Google Colab, and Visual Studio Code (version 1.106.0; Microsoft). Two analytic datasets were created: a case-level dataset consisting of 212 patients and a longitudinal cohort-level dataset comprising 478 time-indexed observations. Data preprocessing procedures included data cleaning, harmonization of tumor staging and Gleason grade groups, elimination of duplicate records, and construction of patient-specific longitudinal sequences. At the case level, continuous variables were compared using the Mann-Whitney U test, while categorical variables were assessed using the Fisher exact test, with a significance threshold set at P<.05. For the cohort-level analysis, mixed-effects linear and logistic regression models were used to evaluate differences between survival groups. The differences between groups, including those stratified by comorbidity status, were assessed using the log-rank test.

The overall ML model development pipeline is illustrated in Figure 2A. To minimize information leakage, data were split at the individual case level into training (80%) and testing (20%) sets using a random group-splitting strategy [25]. Missing values in both static and longitudinal variables were addressed using extremely randomized trees–based iterative imputation, a nonparametric approach derived from the random forest algorithm that accommodates mixed data types and captures nonlinear relationships [26]. Categorical variables were transformed using one-hot encoding, and continuous variables were standardized using z score normalization to ensure comparable feature scales. All preprocessing steps, including imputation, scaling, and encoder fitting, were performed solely on the training dataset and subsequently applied to the testing dataset.

A recurrent autoencoder (RAE) was trained on the training cohort to capture temporal patterns in laboratory and treatment histories (Figure 2B) [20]. PFS duration served as the primary temporal input, with all laboratory variables recorded as baseline measurements before each treatment decision. PSA measurements were used exclusively for PFS event determination according to biochemical progression criteria and were not incorporated as predictive features in model training. The PFS event and time were used only as indexing variables, not as features, to avoid learning in the autoencoder, thereby ensuring that the autoencoder learned purely unsupervised representations of longitudinal patient trajectories without exposure to the outcome. The RAE encoder, composed of stacked gated recurrent units (GRUs), transformed each patient’s sequence into a low-dimensional latent embedding, which the decoder then reconstructed to preserve the original temporal patterns. These latent embeddings were concatenated with case-level static features (eg, race and morbidities) to form a multimodal feature vector. To prevent data leakage, the RAE was trained exclusively on the training dataset, and the learned embeddings were fixed before application to the testing set. Hyperparameters included a hidden layer size of 128, an embedding dimension of 32, a learning rate of 0.01, 100 epochs, and a batch size of 32.

Ethics approval for this study was granted by the Medical and Health Research and Ethics Committee, Ministry of Health, Brunei Darussalam (MHREC/MOH/2024/22(1)), on January 25, 2025. All patient data were replaced with pseudonymous research IDs since we received the data from the hospital, and no personally identifiable information was accessed or retained within the analytic or modeling environment.

All analyses were performed in Python using established open-source libraries, including pandas and NumPy for data manipulation; scikit-survival for Cox proportional hazards modeling, RSF, GBS analysis, concordance index estimation, and time-dependent AUC; lifelines and pycox for survival modeling (including DeepSurv); PyTorch (torch, torchvision, and torchaudio) for deep learning implementation; fancyimpute for iterative imputation; scikit-learn for data splitting and preprocessing (GroupShuffleSplit and StandardScaler); and SHAP for model interpretability.

The study included 212 patients at the case level, contributing a total of 478 longitudinal observations over 6 years (detailed participant characteristics are presented in Multimedia Appendix 2). PFS events were reported in 45 observations during follow-up. At the case level (n=212), PFS outcomes differed according to tumor-related factors, including prostate cancer stage and Gleason grade group, whereas other comorbid conditions did not differ significantly between PFS outcome groups. Furthermore, at the cohort level (478 observations), significant differences between observations with and without a PFS event were observed in calcium (P<.001), chloride (P<.001), sodium (P<.001), red blood cell count (P<.01), hemoglobin (P<.01), neutrophils (P<.01), eosinophils (P<.01), baseline PSA (P<.001), and posttreatment PSA (P<.001).

The PFS analysis with significant variables is presented using the Kaplan-Meier curves in Figure 3. The overall PFS curve demonstrates a steady decline over time, with the probability of remaining progression-free markedly decreasing after approximately 2 years of follow-up. Stratified analyses revealed significant differences across merged clinical stages (log-rank test, P=.0002), with advanced stages showing poorer PFS than earlier stages. Similarly, significant variation in PFS was observed among different prostate cancer treatment groups (log-rank test, P=.0009), indicating heterogeneity in outcomes depending on the therapeutic approach. Furthermore, PFS differed significantly according to place of treatment (log-rank test, P=.0033), with the inpatient treatment group demonstrating consistently lower survival probability over time.

The discriminative performance, as measured by the time-dependent AUC at 1, 3, and 5 years (Figure 4), and the Harrell C-index (Table 1), was used to develop and compare 4 survival models. Evaluation of model performance at 4 and 5 years indicated that ML-based models outperformed the traditional Cox proportional hazards models. The Cox model demonstrated limited discriminative ability, with time-dependent AUCs of 0.706 and 0.833 at 4 and 5 years, respectively, along with a C-index of 0.531. In contrast, the tree-based models, RSF and GBT-Survival, showed markedly superior discrimination. The RSF model demonstrated strong discriminative performance at 4 and 5 years (AUCs of 0.941 and 0.917, respectively; C-index 0.906), while the deep survival model (AUCs of 0.941 at 4 years and 0.917 at 5-year; C-index 0.719) and the GBS model (AUCs of 0.765 at 4 years and 0.833 at 5-year, C-index 0.844) demonstrated moderate discriminative performance.

Calibration performance of the PFS models was assessed using the IBS, where lower values indicate better agreement between predicted and observed outcomes (Table 1). Among the evaluated models, the Cox proportional hazards model demonstrated the best calibration accuracy with the lowest IBS of 0.0527. The GBS model demonstrated comparable calibration performance (IBS=0.0590), followed by the RSF model (IBS=0.0698). In contrast, the deep neural network survival model exhibited the highest IBS (IBS=0.0887), indicating comparatively poorer calibration despite its better discriminative performance.

Although similar time-dependent AUC values were observed between the RSF and deep neural network survival models at specific predefined time points, the limited number of events in the study may reduce the sensitivity of AUC to detect subtle performance differences between models. Therefore, we placed greater emphasis on the C-index and IBS, as these metrics provide a more comprehensive assessment of survival model performance across the entire follow-up period compared with isolated time-point AUC measurements. According to both discriminative and calibration performance, the RSF model demonstrated stable and well-calibrated risk prediction over time.

Model explainability for PFS prediction using the RSF was assessed using SHAP, providing both global and individual-level interpretability (Figure 5). At the global level, the SHAP feature importance heatmap (Figure 5A) indicated that biomedical and laboratory parameters, particularly baseline sodium, alanine aminotransferase (ALT), mean corpuscular hemoglobin (MCH) level, weight, and the treatment variable (Rx_C_D_Merged), were the most influential contributors to the model’s PFS predictions, while demographic variables showed comparatively lower impact. These features highlighted the dominant role of laboratory biomarkers in capturing progression risk in patients with prostate cancer. SHAP values reflect model attribution, not a causal effect.

At the individual level, with the sample patient, SHAP bar, force and waterfall plots (Figure 5B) illustrated how specific features contributed to the predicted PFS for a representative patient, demonstrating the direction and magnitude of each feature’s effect on expected survival time. ALT, sodium, MCH, platelet count, and treatment category increased the predicted risk of progression, whereas other features had minimal or neutral influence. These outputs offer clinical interpretability of how individual patient characteristics contribute to PFS predictions, thereby enhancing transparency, trust, and potential clinical utility in the prognosis of prostate cancer.

The study findings demonstrated that ML-based models showed improved internal discrimination compared with the traditional Cox proportional hazards model in predicting PFS among patients with prostate cancer, particularly at long-term horizons. While early PFS prediction at 1 to 3 years was not feasible due to insufficient progression events and limited comparable risk sets, robust discrimination emerged at 4 and 5 years. This reflects the delayed progression patterns and sluggish natural history of prostate cancer that are frequently seen in clinical practice [37-39]. Among the evaluated models, RSF showed the most balanced performance, achieving high time-dependent AUC values and stable calibration as measured by the IBS. These findings highlight the importance of aligning evaluation time points with disease biology and reinforce the limitations of early event prediction risk in prostate cancer cohorts [8,12,40].

Additionally, the study provided an explainability output from the RSF model using SHAP, which allows addressing the “black-box” barrier for clinical adoption of AI in oncology [19]. The SHAP-based global and individual explanations revealed that the laboratory biomarkers and treatment-related variables dominated the prediction of PFS. This aligns with emerging evidence that data-driven survival models capture complex, nonlinear associative predictors of risk that are often missed by conventional regression-based approaches [20,41]. Furthermore, the SHAP force and waterfall plots at the individual patient level enable clinicians to understand not only what the model predicts but how a particular prediction is made, which can further improve trust, auditability, and shared decision-making in clinical oncology [19,42,43].

Moreover, the application of a RAE to model the longitudinal cohort structure enabled effective compression of repeated observations into latent representations, providing a novel and efficient approach for managing complex temporal data while preserving clinically relevant information. The performance of RSF and GBS models in this cohort was consistent with earlier reports demonstrating the superiority of tree-based methods for cancer survival prediction [19,44]. Additionally, similar to previous studies on prostate cancer, RSF models outperformed Cox models when handling high-dimensional laboratory and longitudinal data, particularly in heterogeneous real-world cohorts [36,45]. The underperformance of Cox regression in this setting supports the rationale for exploring ML survival methods, particularly RSF, which can capture complex nonlinear interactions without proportional hazards assumptions. Although earlier studies reported strong short-term predictive performance, our findings suggest that PFS prediction in prostate cancer is inherently limited at early points [39,40]. The comparatively weaker calibration of the deep neural network survival model highlighted that higher discrimination does not necessarily translate into reliable absolute risk estimation.

According to the SHAP model, sodium, ALT, and hematological parameters such as MCH and platelet count emerged as influential predictors of PFS. Dysnatremia has been associated with cancer progression and systemic inflammation, while elevated ALT may reflect hepatic stress, metabolic dysfunction, or treatment-related toxicity, all of which have been linked to poorer oncological outcomes [46-48]. Moreover, studies have indicated that the aspartate aminotransferase to ALT ratio can be used as a prognostic biomarker for prostate tumor progression [49,50]. These findings support the growing recognition that routine laboratory biomarkers capture systemic physiological states relevant to cancer progression beyond tumor-centric variables alone [51]. However, the prominence of sodium and ALT in our SHAP analysis does not indicate that these markers are biologically superior to other prognostic factors; rather, within this dataset and modeling framework, they emerged as influential predictive contributors that improved the model’s ability to estimate risk.

This study has several limitations that should be acknowledged. First, the relatively small sample size, while sufficient for internal model development and evaluation, may restrict the generalizability of the results to other populations, particularly in light of regional differences in clinical practice and disease characteristics. This also limits the model’s ability to capture early PFS events. Performance estimates are likely optimistic due to limited event counts. Despite these limitations, our study offers important novel contributions that address key gaps in the literature. No prior study, to our knowledge, has collected the comprehensive range of variables included here or applied a sequential approach to define PFS, whereas most existing studies use a single PFS event per patient, limiting clinical applicability. Second, the retrospective nature of the study may introduce bias due to missing data, unmeasured confounding, and evolving treatment strategies over time. Furthermore, model performance was evaluated using a held-out testing dataset; although cross-validation was applied to reduce overfitting, the absence of external validation remains a limitation.

Third, although the analysis was retrospective and conducted at a single institution, the TBCC functions as the only comprehensive oncology facility in the country and is responsible for the management of nearly all prostate cancer cases nationwide. Nevertheless, the findings highlight the potential utility of ML models incorporating SHAP-based explainability for prognostic risk assessment in prostate cancer. Finally, radiographic progression according to PCWG criteria was not incorporated into the PFS analysis because of the limitations of the imaging data available for this cohort. Future research should include independent external validation and multi-institutional, multiomics datasets with larger and more diverse cohorts to strengthen the generalizability and robustness of these findings.

Our study demonstrates that the ML-based survival models, particularly tree-based ensemble methods, substantially outperform the traditional Cox proportional hazards model in predicting PFS among patients with prostate cancer, especially at longer follow-up horizons. Among the evaluated approaches, the RFS model achieved the most balanced performance, combining strong long-term discrimination with stable calibration, thereby offering a robust framework for prognostic modeling in real-world prostate cancer cohorts. The integration of SHAP-based explainability addressed a key barrier to clinical adoption by providing transparent global and individual-level insights into model predictions highlighting the dominant role of laboratory biomarkers and treatment-related variables in driving PFS risk predictions in this cohort. The application of a RAE to encode longitudinal cohort data further enhanced the model’s ability to capture the complex temporal patterns while reducing dimensionality. These findings support the clinical potential of explainable ML models for personalized risk stratification in prostate cancer, while also emphasizing the need for external validation, larger multi-institutional, multiomics datasets, and prospective evaluation prior to routine clinical implementation.