Our study proposes an explainable AI-based exploratory AI recommender (recommender henceforth) designed to systematically uncover inherent data relationships and optimize the design of predictive cohort studies. The methodology follows three stages (Figure 1): (1) establishing a baseline predictive model, representing the initial analytic stage of a research team; (2) deploying our recommender to generate actionable insights and recommendations for refining study design; and (3) evaluating the recommender’s utility by quantifying the performance gains of models augmented with these recommendations. We apply this methodology to a clinical case study: a prognostic model for the 1-year risk of FRI, developed using DataLoch data.

The analysis was conducted on a high-dimensional clinical dataset provided by Lothian DataLoch, a secure, regional data repository linking primary and secondary care data in South East Scotland [31]. Because DataLoch holds the routinely collected records of all patients registered with contributing general practices, the cohort is population-based rather than a selected sample. It represents adults aged 50 years or older in one Scottish health board region. This study included all individuals aged 50 years or older, permanently registered for at least 1 year with a Lothian general medical practice (GP) that contributed to the DataLoch dataset between January 1, 2015, and February 28, 2020. This pre–COVID-19 period was specifically chosen to mitigate potential biases on FRI that might arise from pandemic-related changes in health care–seeking behavior. The index date was defined as the latest of January 1, 2015, the patient’s 50th birthday, or the date on which 1 year of continuous registration was completed. Each patient entered the cohort once. Follow-up began on the index date and ended at the earliest of the first subsequent FRI, death, GP deregistration, 1 year after the index date, or February 28, 2020. Individuals were excluded if data for age, sex, or the Scottish Index of Multiple Deprivation (SIMD) were missing. The dataset was randomly split into 80% training (176,272 rows) and 20% held-out test data (44,069 rows). This study follows the Transparent Reporting of a Multivariable Prediction Model for Individual Prognosis or Diagnosis (TRIPOD+AI) reporting guidelines [32] (checklist in Multimedia Appendix 1). No study protocol was published, and the study was not registered.

The primary outcome was the time from cohort entry to the first occurrence of FRI, with related injuries defined as low-energy fractures of hip, wrist, and proximal humerus. Diagnosis of FRI was defined by clinical codes from primary care (Read version 2), hospital admissions (International Statistical Classification of Diseases and Related Health Problems 10th Revision [ICD-10]), and emergency department records (Table S1 in Multimedia Appendix 2). Death before FRI was treated as a censoring event, so the analysis modeled the cause-specific hazard of FRI.

The primary analysis used 59 features selected from 409 available features based on established literature following the knowledge-driven approach. These include 4 patient demographic characteristics, 3 lifestyle risk factors, 28 comorbidities, the eFI [5], and exposure to 23 Fall-Risk-Increasing Drugs (list of features in Table S4 in Multimedia Appendix 2). Comorbidities were identified from all available primary care and hospital records dated before the index date. Exposure to each medication class was defined as at least 1 recorded dispensing during the 84 days before the index date. The same ascertainment window was applied to all medication classes. BMI, eFI, smoking status, and alcohol consumption were obtained from the most recent record before the index date. Age was calculated at the index date. The SIMD quintile was derived from the patient’s residential postcode. For patients with more than 1 SIMD record, the most frequent quintile was used. Age and eFI were continuous variables, and BMI category and SIMD quintile were coded using ordered integer codes. Ethnicity, smoking status, and alcohol consumption were one-hot encoded, with missing values treated as separate categories. We chose this approach because missingness in routine records can itself be informative, and a model applied to routine data will encounter the same missing values [33]. The remaining predictors were binary. Full predictor definitions and preprocessing are provided in Note S1 in Multimedia Appendix 3.

First, we established a baseline predictive model using the standard multivariable CPH model. This model represents the starting point of a conventional, literature-driven predictive modeling study [28]. The baseline model was fitted on the training data using 59 selected candidate predictors, deliberately omitting advanced techniques such as feature selection, nonlinear terms, or interaction terms. Its performance on the held-out test set served as the baseline for evaluating recommender-driven refinements.

Second, we deployed the recommender to systematically explore underlying data relationships and generate concrete recommendations for refining the above knowledge-driven predictive study design. This began with fitting an advanced exploratory AI model to the training data to capture complex, nonlinear relationships and potential interactions between all features and outcomes. For this purpose, we used random survival forest (RSF) [34] as the exploratory AI model. Unlike typical machine learning applications, the RSF was used purely to capture high-dimensional data structures rather than for final clinical prediction. RSF was fitted with 300 trees, with min samples leaf as 100, and random seeds as 20 (Table S5 in Multimedia Appendix 2). The RSF returns a risk score, calculated as the sum of the predicted cumulative hazards over the distinct event times in the training data. We used the logarithm of this score as the model output to be explained. Next, the RSF model was interpreted using explainable AI techniques, such as SHAP [20], to translate “black-box” patterns into actionable insights. SHAP quantifies each feature’s contribution (termed FAs) to the RSF’s prediction for each patient (see Note S2 in Multimedia Appendix 3 for further details on FAs). FAs are expressed on the log scale of the risk score. On this scale, an additive effect of a feature yields an FA that depends only on that feature, and multiplying the score by a constant adds the same constant to every prediction. The FAs therefore do not depend on the arbitrary scale of the score. Standard FA analysis often relies on population-level averaging, which can obscure risk factors that are clinically significant only within specific subgroups. We therefore analyzed FAs in 2 subcohorts targeted at different levels of predicted risk: a low-risk and high-risk subcohort. Each subcohort had a reference patient: the mean feature vector of all training patients without an FRI (low-risk subcohort) or with an FRI (high-risk subcohort). The low-risk subcohort included training patients whose RSF log risk score differed from that of the low-risk reference patient by less than a margin, measured in SDs of the log risk score among training patients without an FRI. The high-risk subcohort was defined in the same way, relative to the high-risk reference patient and training patients with an FRI. The margin was chosen on the training data from a prespecified grid (0.025 to 1.0 SD). We selected the smallest value at which the recommendations did not change when the margin was widened further, because a narrower margin keeps each subcohort close to its reference patient (Note S4 in Multimedia Appendix 3). FAs in each subcohort were computed relative to the same reference patient and therefore describe deviations from a typical patient of that group (Note S2 in Multimedia Appendix 3). The FAs from these cohorts were then processed to derive modeling recommendations. While computing FAs might be more intensive on average than training a baseline Cox model, this exploratory phase is a one-time cost that precedes the deployment of the more efficient augmented model.

Our recommender offers the following 3 types of recommendations (Figure 2): (1) Feature exclusion (Figure 2E). A feature was recommended for exclusion if the upper bound of the 95% CI of its mean absolute FA was below a data-driven significance threshold. The threshold was defined as 5% of the standard deviation of mean absolute FA values across all candidate features, a heuristic to distinguish features with negligible contributions. (2) Nonlinear relationships (Figure 2D). We screened continuous and ordinal features for potential nonlinear effects. If a feature has a linear effect in the Cox model, its FAs are a linear function of its values on the log-risk scale. For each feature, we therefore regressed its FAs on its values using both a straight line and a flexible function. For ordinal features, the flexible function estimated each category’s mean FA by using 0/1 indicator variables, allowing changes between successive levels to differ. For continuous features, it was a natural cubic spline with 4 degrees of freedom. The gain in the coefficient of determination (R²) from the straight line to the flexible function, that is, the additional variation in the FAs explained by allowing a nonlinear shape, was tested with a nested F test. A feature was recommended if its adjusted P value (see Statistical Analysis) was below .05 and its gain in R² exceeded 0.02 [35]. The validation of the nonlinearity rule is through a simulation considering different cases of linear and nonlinear data generation processes (Table S8 in Multimedia Appendix 2); (3) Feature interactions (Figure 2; A-C). Interactions were screened in 4 steps using training data only (Note S3 in Multimedia Appendix 3). First, features whose FAs were not explained by their own values (lower limit of the 95% bootstrap interval of the normalized residual SD ≥0.10) were selected as focal features. Second, each focal feature was paired with every other candidate feature, and we tested whether the second feature's within-stratum FA contrast (binary features) or product-term coefficient (continuous or ordinal features) varied with the focal feature, using Wald tests with robust standard errors. Under an additive model with a fixed reference, these quantities do not change with differences between strata, so a pair is flagged only if the model is nonadditive on the explained scale. Pairs with an adjusted P<.05 (Bonferroni within pair, then Benjamini-Hochberg across pairs) and an effect of at least 10% of the mean absolute FA advanced. The interaction screening rule is validated through simulation on only additive null interaction, weak interaction (0.1), and strong enough interaction with type I error below 10%, and interaction detection rate above 90% in both interaction cases (Table S9 in Multimedia Appendix 2). Third, pairs were retained if a Cox score test against a main-effects reference model had a Benjamini-Hochberg-adjusted P<.05. Fourth, the number of added parameters was capped using the sample-size criterion of Riley et al [35].

Third, we evaluated the effectiveness of the recommender. We constructed augmented Cox models that integrated these recommendations (excluding features, adding nonlinear terms and interactions) to predict the same 1-year risk of FRI. For nonlinear terms, continuous features entered the model with quadratic terms, and ordinal features as the indicator variables described above. All recommendations were derived in one pass from the FAs of the same exploratory model. Feature exclusion was applied first, and the nonlinearity and interaction rules were applied only to the retained features. To preserve the interaction hierarchy, an interaction was included only when both main effects were retained. Therefore, the order in which the recommendations are added does not change the final model. We then compared the predictive performance of these augmented models against the baseline on the held-out test set. We also compared the augmented Cox models with 2 approaches commonly used in routine prediction modeling. These were a Cox model with restricted cubic splines for continuous predictors (Cox-RCS) and a penalized Cox model including all main effects and all pairwise interactions (Coxnet; least absolute shrinkage and selection operator [LASSO], with the penalty selected by 10-fold cross-validation on the training data). Furthermore, the clinical plausibility of generated recommendations was discussed in light of the published literature (See Discussion).

To examine the Recommender in a very high-dimensional setting, we applied it to all 409 available features, without literature-based selection. The pipeline was unchanged: the RSF was fitted on all 409 features, FAs were computed in newly selected low-risk and high-risk subcohorts, and the 3 recommendation rules were applied.

To assess the applicability of the Recommender across datasets, outcomes, and clinical settings, we applied the Recommender on 3 additional, publicly available time-to-event datasets: German Breast Cancer Study Group 2 (GBSG2; breast cancer), AIDS Clinical Trial (ACT; human immunodeficiency virus), and PeakVO2 (systolic heart failure) [36-38]^.

Baseline characteristics of the study samples were described using frequencies and percentages for categorical variables and means (SDs) for continuous variables. Statistical analyses, model development, and testing were conducted in Python (version 3.12.5; Python Software Foundation), except for the proportional hazards tests, which were conducted in R (version 4.5.3) with the survival package (version 3.8.6; R Core Team). All preprocessing steps (encoding of categorical features, standardization of continuous features, and the centering used for quadratic terms) were estimated on the training set only and then applied unchanged to the test set. Discrimination was measured with Harrell concordance index (C-index) and Uno C (less sensitive to the censoring distribution, truncated at the 90th percentile of test-set follow-up times) [39]. Calibration was assessed with the conventional calibration slope. This is the coefficient of the complementary log-log transformation of the predicted survival probability at a fixed horizon in a Cox model fitted to the test set, with a Wald 95% CI. The horizons were 365 days for DataLoch, 1500 days for GBSG2, 320 days for ACT, and 1095 days for PeakVO2. Calibration plots compared the mean predicted survival with the Kaplan-Meier estimate in 5 quantile-based groups of predicted survival. Overall prediction error was summarized by the integrated Brier score (IBS). Unless stated otherwise, 95% CIs are percentile intervals from 1000 bootstrap resamples of the test set. Differences in C-index relative to the baseline model were computed on the same resamples. Intervals in the calibration plots were 95% percentile CIs from 500 bootstrap resamples of the test set.

P values in the recommendation rules were adjusted with the Benjamini-Hochberg procedure, with the false discovery rate (FDR) controlled at 5%. For the nonlinearity rule, the adjustment was across all candidate features within each subcohort. For the interaction screen, the repeated tests of each pair were first combined by Bonferroni adjustment, and then adjusted across all unique feature pairs. For the Cox score test, the adjustment was across the pairs that passed the screen. Proportional hazards were assessed for the augmented Cox model fitted to the training data, with global and per-term Schoenfeld residual tests. All tests were 2-tailed, and P<.05 was considered significant.

All data from DataLoch were linked and deidentified by the DataLoch service (Edinburgh, United Kingdom) and analyzed within a Secure Data Environment (reference: DL_2024_030). DataLoch enables access to deidentified extracts of health care data from the South-East Scotland region to approved applicants [40]. The study was reviewed and received approval under delegated authority from a regional National Health Service Research Ethics Committee (REC 22/NS/0093) and Caldicott Guardian. Individual patient consent was not required. Clinical trial number: not applicable. The public datasets are de-identified and openly available. Patients and the public were not involved in the design, conduct, reporting, or dissemination of this study.

The DataLoch cohort included 220,341 patients (mean age 64, SD 12 years), of whom 119,026 (54%) were females. Ethnicity was self-identified as White (Gypsy/Traveller, Irish, Other British, Other White, Polish, Scottish) in 187,254 (85%), Asian (Bangladeshi, Chinese, Indian, Other Asian, Pakistani) in 3093 (1.4%), and Black (African, Caribbean, Other Black) in 757 (0.3%) patients. The most prescribed medications were angiotensin-converting enzyme inhibitors (46,438 patients, 21.1%), beta-adrenoceptor blocking drugs (34,809 patients, 15.8%), and calcium-channel blockers (33,054 patients, 15%). The most prevalent comorbidities were hypertension (94,477 patients, 42.9%) and depression (58,260 patients, 26.4%). During follow-up, 8489 (3.85%) patients experienced an FRI. Median observed follow-up was 1 (IQR 1-1) year, with a total of 207,278 person-years of observation. The crude FRI incidence rate was 40.95 per 1000 person-years. The crude FRI incidence rate was 51.10 per 1000 person-years in women and 29.14 per 1000 person-years in men. In the time-to-FRI analysis, 211,852 patients (96.15%) were censored, including 4040 (1.83% of the cohort) who died before experiencing an FRI. The training set included 176,272 patients with 6791 FRIs, and the test set included 44,069 patients with 1698 FRIs. Table S4 in Multimedia Appendix 2 contains the patients’ characteristics of the study.

As shown in Table 1, the baseline CPH achieved a C-index of 0.715 (95% CI 0.702-0.727) and a conventional calibration slope at 365 days of 0.961 (95% CI 0.909-1.013). The exploratory RSF model achieved a C-index of 0.716 (95% CI 0.703-0.727). The forest discriminated similarly to the linear, additive Cox model, indicating that most of the predictive information in the candidate features was captured by their main effects.

Figure 3 details the process for conducting FA analysis and deriving recommendations. The subcohort margin, chosen on the training data, was 0.05 SD of the log risk score. It selected 8208 patients into the low-risk subcohort and 5608 patients into the high-risk subcohort from the training set (n=176,272). The recommendations did not change when the margin was widened to 0.1 SD. The FA analysis was conducted in the 2 subcohorts separately, and the resulting recommendations are summarized in Table S6 in Multimedia Appendix 2. First, it recommended 15 features for exclusion, containing 8 medication classes and 7 comorbidities. Second, among the remaining 44 features, 4 features were recommended for nonlinear terms. Age and the eFI entered with an added squared term, and BMI category and SIMD quintile entered as level indicators (Table S15 in Multimedia Appendix 2). Finally, the interaction screen performed 1807 tests assessing 946 unique pairs among the 44 remaining features. Of these pairs, 370 passed the FA screen (step 2), and 4 passed the Cox score test (step 3). All 4 were within the allowance of 540 parameters set by the sample-size criterion (step 4). The 4 recommended interactions were age with stroke, age with osteoporosis, age with loop diuretics, and age with the eFI.

We integrated the 3 types of recommendations into augmented Cox models independently and cumulatively to evaluate their individual and collective impact (Table 1 and Figure 4). In the independent assessment, an augmented Cox model incorporating only the feature exclusion recommendations (59 features reduced to 44) demonstrated no degradation in performance compared with the baseline Cox model (difference in C-index: 0.000, 95% CI –0.001 to 0.001; P=.87), while reducing the number of parameters from 73 to 58. Adding the nonlinear terms increased the C-index to 0.718 (95% CI 0.707-0.730; difference in C-index: 0.003, 95% CI 0.001-0.006; P=.02). Adding the 4 interactions did not improve discrimination (C-index 0.715, difference −0.001, 95% CI –0.002 to 0.000; P=.15). The model with all recommendations achieved the highest C-index, 0.719 (95% CI 0.708-0.731), an improvement of 0.004 (95% CI 0.001 to 0.007; P=.01) over the baseline model, with 70 parameters. Uno C followed the same pattern (0.713-0.718), and the IBS was almost identical across models (0.0228-0.0230), as expected with an event rate below 4%.

Compared with the 2 comparators, Cox-RCS achieved a C-index of 0.718 (95% CI 0.707-0.731) with 79 parameters, similar to the model with the recommended nonlinear terms. The model with all recommendations demonstrated slightly higher discrimination performance than Cox-RCS with fewer parameters (70 vs 79). Coxnet had lower discrimination performance than the baseline model (C-index 0.706, 95% CI 0.694-0.719). It retained 18 interactions, which included all 4 interactions recommended by the recommender.

Calibration at 365 days was close to ideal for all Cox models, with calibration slopes of 0.961 (95% CI 0.909-1.013) for the baseline model and 0.962 (95% CI 0.910-1.015) for the model with all recommendations (Figure 4). The final augmented model had 70 coefficients. The coefficient of each parameter in the augmented model is shown in Table S10 in Multimedia Appendix 2. With 6791 events in the training set, this corresponds to about 97 events per parameter. For the final augmented model, the global Schoenfeld residual test showed no evidence against the proportional hazards assumption (χ²₇₀=64.44; P=.67), and 4 of the 70 terms had a P value of <.05 (Table S11 in Multimedia Appendix 2).

When applied to all 409 available features, the recommender recommended 304 exclusions, 4 nonlinear terms, and 2 interactions (Table 2). The nonlinear terms were recommended for age, the eFI, BMI category, and SIMD quintile. The interactions were between age and the eFI, and lipid-regulating drugs and drugs used in megaloblastic anemia. As in the primary analysis, the model with all recommendations achieved the highest C-index, 0.722 (95% CI 0.710-0.734). It showed an improvement of 0.009 (95% CI 0.004-0.014; P<.001) over the baseline model, with much fewer parameters (129 vs 423). Exclusion improved discrimination by 0.006 (95% CI 0.002-0.010; P<.001) and nonlinear terms by 0.003 (95% CI 0.001-0.005; P=.02). Adding the 2 interactions did not change discrimination (difference in C-index: 0.000, 95% CI −0.001 to 0.001; P=.30), as expected with so few terms. Uno's C increased from 0.712 to 0.720. The calibration slope improved from 0.846 (95% CI 0.797-0.894) for the baseline model to 0.952 (95% CI 0.901-1.003) for the model with all recommendations (Figure 5). The model with all recommendations had a significantly higher C-index than Cox-RCS (0.715; difference in C-index 0.006, 95% CI 0.002-0.010; P<.001) with about one-third of the parameters (129 vs 429). Coxnet could not be fitted because all pairwise interactions of 409 features give more than 80,000 candidate terms. For the final augmented model, the global Schoenfeld residual test showed no evidence against the proportional hazards assumption (χ²₁₂₉=137.66; P=.28), and 10 of the 129 terms had a P value of <.05 (Table S12 in Multimedia Appendix 2).

Table 3 summarizes the results in the 3 public datasets. Detailed recommendations and model performances are provided in Table S13 and S14, respectively, in Multimedia Appendix 2. In each dataset, the augmented model with all recommendations had a higher C-index than the baseline Cox model (GBSG2: 0.701 vs 0.667; ACT: 0.749 vs 0.743; PeakVO2: 0.677 vs 0.670).

No interaction entered the model in any of the 3 datasets. The Recommender identified some candidate interactions, such as beta-blockers with resting heart rate in PeakVO2. However, in these smaller datasets, the main effects and nonlinear terms already used the parameter allowance of the sample-size criterion.

In GBSG2 (Figure 6), only nonlinear terms were recommended. The resulting model had a higher C-index than Cox-RCS (0.701 vs 0.687) with fewer parameters (14 vs 23), and the calibration slope moved from 1.24 to 0.94. Coxnet, which retained 2 main effects and 3 interactions, had a lower C-index than the baseline model (0.659). In ACT and PeakVO2 (see Figures 7 and 8, respectively), the augmented model performed similarly to the Cox-RCS and Coxnet. The Recommender recommended nonlinear terms for most continuous features in both datasets, so the augmented model was close to Cox-RCS. Coxnet retained few main effects (7 in ACT and 5 in PeakVO2) and many interactions (31 and 25). This suggests that some of its interaction terms may approximate nonlinear effects rather than true interactions. Coxnet had the best calibration in these 2 datasets (slopes 0.73 and 1.00), probably because penalization shrinks the coefficients and so moderates extreme predictions. The baseline, augmented, and Cox-RCS models had calibration slopes of 0.61 to 0.72 in these 2 datasets, which shows some overfitting. Among the comparators, Coxnet had the best calibration (slopes 0.73 and 1.00), probably because penalization shrinks the coefficients and so moderates extreme predictions. The LASSO-augmented model applies the same penalization to the recommended terms. It had C-indexes of 0.702, 0.761, and 0.677 and calibration slopes of 1.03, 1.02, and 0.89 in GBSG2, ACT, and PeakVO2, respectively.

Across the 3 datasets, the models built from the recommendations had similar or higher discrimination than the comparators. With penalization, they also had good calibration.

In this study, we developed and evaluated an exploratory AI recommender designed to support predictive study design in high-dimensional clinical datasets. Rather than using AI as an opaque final prediction model, the Recommender uses a flexible explainable-AI module as an exploratory engine to identify data-driven modeling insights that can then be translated into transparent statistical models. In a fall risk case study, the recommender recommended excluding 15 of 59 candidate features, adding nonlinear terms for 4 features, and adding 4 interactions from 946 candidate pairs. With all 409 available features, it recommended excluding 304 features, adding nonlinear terms for 4 features, and adding 2 interactions from 7021 candidate pairs. The augmented Cox model had a statistically significant gain in discrimination over the baseline model (C-index 0.719 vs 0.715), with 3 fewer parameters and similar calibration. With all 409 features, the gain was larger (0.722 vs 0.713), with about one-third of the parameters (129 vs 423), and the calibration slope moved closer to 1 (0.952 vs 0.846). The exploratory forest has similar discrimination as the baseline Cox model (0.716 vs 0.715), which indicates that the candidate features carry little nonlinear or interaction signal beyond their main effects. Large gains were therefore not expected. The main benefits of the Recommender are a simpler, transparent model and a short list of hypotheses for formal testing. Because AI is used only to explore the data, the final model keeps the auditability and transparency that are required in regulated clinical environments.

The Recommender’s outputs agree with published evidence. Age and frailty are established predictors of falls and related injuries [41], and their associations with fall risks are known to be nonlinear [42]. The recommended interaction between age and stroke also agrees with a population-based study, in which the relative risk of hip fracture after stroke was higher in younger than in older patients [43]. Beyond confirming existing knowledge, our Recommender also identified some novel hypotheses that deserve further exploration. For example, with all 409 features, it recommended an interaction between lipid-regulating drugs and drugs used in megaloblastic anemias. Statins alone are not consistently associated with falls in older adults [44], but they are associated with myopathy [45] and, less consistently, with peripheral neuropathy [46]. Vitamin B12 or folate deficiency can also cause neuropathy and impaired proprioception [47]. When these effects occur together, their combined effect on balance and gait may be larger than either effect alone. If this hypothesis is confirmed in other data, it may inform medication review in older adults who take both types of drugs. Crucially, the Recommender supports rather than replaces clinical judgment. The Recommender’s transparent design allows researchers to critically evaluate such data-driven recommendations.

This work extends previous research to predictive model designs in several ways. First, the Recommender is not simply a feature selection tool. Existing approaches, including robust variable selection and regularization-based methods for Cox models [48], provide rigorous strategies for deriving effective predictor sets, improving model stability, and reducing overfitting within prespecified modeling frameworks. However, these methods are primarily designed to determine which variables should be retained in a final model. In contrast, our Recommender combines a flexible nonlinear exploratory model with explainable AI to generate broader design-relevant insights, including candidate feature exclusions, potential nonlinear effects, interaction structures, and subgroup-specific predictor behavior. These outputs are intended not only to support model development, but also to help researchers understand why a feature may be informative, how its contribution may vary across patient subgroups, and which feature relationships may warrant formal statistical evaluation. In this sense, the Recommender complements rather than replaces established biostatistical methods by prioritizing clinically and statistically plausible signals from a large candidate space for targeted follow-up using conventional tools, including penalized Cox models, restricted cubic splines, interaction terms, and likelihood ratio tests. Closest to our work, Gosiewska et al [29] used a flexible surrogate model to derive feature transformations for an interpretable model, and Liu et al [30] used machine learning to identify new predictors that were then added to Cox models. However, each method addresses only one type of model change, and to our knowledge, neither screens for interactions or recommends feature exclusions. Our Recommender differs in proposing all 3 types of change from one exploratory model, with formal tests and a parameter limit, rather than adding or transforming features alone.

Second, we compared the Recommender with 2 standard approaches. Cox-RCS adds flexible nonlinear terms for all continuous features. In all datasets, the augmented model discriminated similarly to or better than Cox-RCS, with fewer parameters (eg, 70 vs 79 in DataLoch and 43 vs 78 in PeakVO2). With all 409 features, it had a higher C-index than Cox-RCS (difference 0.006, 95% CI 0.002-0.010) with 129 vs 429 parameters. Thus, the Recommender matched the discrimination of the more flexible Cox-RCS with a simpler model, which used squared terms instead of splines and fewer features. Coxnet selects terms from all main effects and all pairwise interactions by LASSO shrinkage. Its candidate set grows quadratically with the number of features (2624 terms with 59 features), so it could not be used with all 409 features. It also cannot separate nonlinear effects from interactions. In the public datasets, Coxnet retained few main effects but many interactions, so some of its interactions may approximate nonlinear main effects rather than reflecting true interactions. When the same LASSO penalty was applied to the recommended terms, the C-index was higher than that of Coxnet in all 3 datasets. The recommended terms may therefore be a better starting point for a penalized model than all pairwise interactions.

The public datasets also show how the Recommender behaves in smaller datasets with fewer events. In each dataset, the model before interactions already had more parameters than the sample-size criterion allowed, so the criterion correctly prevented additional interaction terms. The calibration slopes of all unpenalized Cox models in ACT and PeakVO2 were consistent with this risk of overfitting. In such small datasets, the recommended terms can be combined with penalization, which gave good calibration in all 3 datasets.

Third, the recommender provides a mechanism for exploring subgroup-specific predictor behavior that may be overlooked by conventional model development. Traditional models often emphasize high-risk populations, yet patients classified as low risk are not necessarily free from risk; rather, potential risk factors may be compensated by other protective factors [49]. Our subcohort analysis showed that nonsteroidal anti-inflammatory drugs (NSAIDs), although contributing minimally in the high-risk cohort, were a significant predictor specifically in the low-risk cohorts. Overall NSAID exposure was similar in patients with and without an FRI (6.5% vs 6.8%), so this pattern is not explained by a crude difference in prescribing. It may instead reflect how NSAIDs are prescribed or avoided in different groups of patients, that is, confounding by indication or contraindication. This provides a hypothesis for future causal studies of NSAID prescribing. Moreover, in the whole-cohort Cox model, NSAIDs were not associated with FRI (hazard ratio [HR] 1.06, 95% CI 0.96-1.17), which illustrates how a subgroup-specific pattern can be hidden in population-level estimates. The Recommender does not calculate fairness metrics, and we did not assess fairness. However, researchers could define strata by protected characteristics, such as sex, ethnicity, or socioeconomic deprivation, and compare FAs across them. By exposing these subgroup-specific deviations, the Recommender provides insights that might help researchers identify bias early in the study design.

Fourth, the recommender bridges the gap between AI power and clinical transparency. AI has demonstrated remarkable predictive capabilities with complex and large-scale biomedical data [50]^, but its opacity often prevents trust and clinical adoption [51]^. By using AI models purely as an exploratory tool to inform conventional statistical models, we address the accuracy-interpretability trade-off without the fidelity issues of post hoc explainable AI [52]^.

Furthermore, our recommender demonstrates some utility in a “hard-to-improve” scenario, as the improvement in C-index (0.004 with the 59 selected features and 0.009 with all 409 features) is statistically significant but numerically modest. Relatively simple models can perform reasonably well on high-quality, large-scale datasets, leading to great difficulty in achieving marginal gains. However, in large-scale population health, even marginal gains in discrimination and calibration can lead to significantly better resource allocation and risk stratification when applied to hundreds of thousands of patients. Our approach is scalable, as it rapidly explores thousands of pairwise interactions; efficient, as it avoids the computational demands of complex architectures like graph neural networks [53,54]; and adaptable, as it is applicable to different models (eg, DeepSurv [55] or diagnosis prediction) and clinical datasets. While SHAP is used here due to its theoretical foundations in game theory, the Recommender is attribution-agnostic. Researchers may consider other techniques like Integrated Gradients [56] or SurvSHAP(t) [57] for time-dependent explanations, depending on the computational budget and the specific exploratory model used.

The recommender offers practical implications for clinical research. It streamlines future data collection as it may reduce the number of variables that need to be collected without loss of predictive performance, can guide feature selection in further causal studies, and enables researchers to move beyond preconceived hypotheses [58]. More broadly, the Recommender encourages researchers to move beyond purely literature-driven or knowledge-based model specification by providing a systematic, data-driven mechanism for discovering important features, nonlinearities, interactions, and subgroup-specific effects. These insights may also inform subsequent causal studies, although they should not be interpreted as causal evidence without further causal design and validation.

The recommender is intended for analysts who develop prediction models, not for clinicians at the point of care. The released code automates all steps, from fitting the exploratory model and computing attributions to choosing the subcohort margin, applying the 3 recommendation rules, checking candidate interactions in a Cox model, and limiting the number of added parameters. The analyst only needs to specify the encoding of each feature (continuous, ordinal, nominal, or binary). Preselection of candidate features is optional, as the analysis of all 409 features shows, and the analyst can prespecify cut points for continuous features when a specific subgroup is of interest (eg, age 65 years according to a common older age threshold [59] ). The analyst should review each recommendation for clinical plausibility before adopting it. In this study, all recommendations were applied as generated so that the method itself could be evaluated.

This study has several limitations. First, our interaction terms were restricted to (already challenging) first-order pairs to ensure computational feasibility. Higher-order interactions were not explored, but are a potential focus for future research, where computational complexity is the main concern [60]. Second, the validity of the recommendations depends on the exploratory model's performance (eg, RSF). Suboptimal discrimination or calibration could introduce bias into FAs, necessitating careful model selection or ensemble approaches. Thus, we advise using the recommendations produced only with high-performing exploratory models (high C-index and stable calibration), and leave investigation into quality metrics of FAs to future work. Third, our analysis focused on large-scale populations rather than individual longitudinal trajectories, leaving temporal dependencies out of scope, though in principle applicable with extensions [57]. Fourth, we screened candidate interactions and checked them in a Cox model on the same training data, without a separate confirmation set. We chose this because a further split would reduce the events available to the exploratory forest and its subcohorts, and so reduce the power of the screen. Confirmation on separate data is a priority for future work in datasets with more events. Moreover, we adopted a random train/test split. The DataLoch extract has no practice identifier, so we could not test for clustering by practice. Patients from the same practice can be in both sets, and any practice-level effect on prescribing or care pathways would not be captured by this safeguard. A temporal split was also of limited use, because follow-up was 1 year and 78.7% of patients had the same index date (January 1, 2015). In our DataLoch analysis, death before an FRI was treated as a censoring event, so the model estimates the cause-specific hazard of FRI. Because competing deaths were uncommon (4040 patients, 1.83% of the cohort) and follow-up was limited to 1 year, we did not perform a Fine-Gray analysis. However, competing deaths mainly affect absolute risk estimates and are unlikely to change the comparison between the baseline and augmented models, which is the focus of this study. The cohort came from a single Scottish region, with few patients from ethnic minority groups (Table S4 in Multimedia Appendix 2) and higher fracture rates than some other United Kingdom regions [7]. The model, and the size of the gains from the recommender, may differ in other populations and health systems. Finally, we emphasize that FAs represent predictive associations, not causal relationships; they denote contributions to the model’s output rather than the clinical outcome itself. Any clinical interventions based on these findings would require subsequent causal discovery or randomized controlled trials.

We demonstrated that explainable AI can uncover complex feature relationships within high-dimensional clinical datasets. This data-driven exploratory approach showed the potential to improve the robustness of predictive cohort study design by augmenting traditional, expert-led methods. These findings also suggest that explainable AI-informed insights can generate novel hypotheses that may warrant further investigation.