Model development choices were informed by a previous review [8] and best-practice recommendations [15,16]. The model was developed using Cox proportional hazards regression and included time-varying covariates and predictor variables available within the EMR system. Performance was assessed in terms of discrimination and calibration using internal-external cross-validation, and the final model was fit using the full dataset.

After consultation with fall prevention experts and hospital geriatricians, we censored observations after 14 days of admission. After 14 days, hospital falls occurred at a lower rate and may be associated with a different clinical context for long-stay patients that may benefit from a more targeted modeling strategy. To ensure the data focused on the period in which fall interventions are initiated and most falls occur, we used EMR data for the first 14 days of admission from patients admitted to 5 Australian hospitals between 2018 and 2021 to develop our model. Our intention was for study findings to be representative of inpatient hospital wards broadly. There were no exclusion criteria applied to eliminate particular subgroups of hospital patients; patients from any ward and of any age were included for model development and validation. In each of our included hospitals, inpatient falls were recorded in the EMR using RiskMan (RiskMan International Pty Ltd), a proprietary patient safety system for tracking and managing adverse events [17]. RiskMan reports were linked to each inpatient admission record using a unique identifier, including date and time stamps. When the time of fall was not recorded, it was assumed to be at midnight on the morning of the date specified to ensure no predictions used data recorded after the fall.

To maximize the potential for generalizability and ease of potential transferability to external hospital systems while eliminating the need for clinical staff time for additional fall risk assessment and data entry, predictor variables used in the model were primarily derived from administrative information routinely recorded for all hospital admissions. This included the source of the admission (eg, whether they came from residential aged care or through the emergency department), type of admitting medical service (eg, neurosurgery or orthopedics), sex (coded as the binary variable “Female”), age in years, as well as years since 2018 to incorporate potential adjustment for longer-term trends in fall risk over the study period. We also included 2 time-varying covariates in the model—the number of previous inpatient falls during the current admission up until the time of prediction and the time since admission, in hours. The time since admission variable was updated every 12 hours throughout the admission.

Because coding for admission source and medical service categories with similar or overlapping meanings varied across hospital facilities, we aggregated these into overarching categories that could be readily applied to data from all hospitals based on clinical groupings, shown in Table 1. Using these overarching clinical groupings also mitigated the risk of categorical variables for admission source or medical service being exclusive to the training or validation datasets. The “Other” categories for these 2 variables were used to include missing data as well. There were no missing data for other covariates included in the model. Restricted cubic splines with 4 knots were used for patient age and the years since 2018, and a restricted cubic spline with 3 knots was used for the hours since admission variable. A spline function with n knots is made up of n+1 piecewise polynomial functions where the knot locations describe the place (value of the variable to which spline is being applied) where these 2 sequential polynomial functions meet. For example, in the full model, the restricted cubic spline applied to patient age has 4 knots, meaning that there are 5 polynomial functions that meet at the 4 locations (values of patient age; described in Multimedia Appendix 1). Further information on use and interpretation of spline terms are illustrated in Regression Modeling Strategies by Harrell [18]. The restricted cubic splines used in this study have an advantage for use in prediction models in that they are constrained to be linear in their 2 tails (above the last knot and before the first knot). This is particularly important considering the use of “years since 2018” as a predictor and the possible use of the model in years that exceed the range of the data used to fit it.

Following advice from Steyerberg and Harrell [16], we developed the model using the full dataset but evaluated model performance using internal-external cross-validation. Internal-external cross-validation is a process in which the folds used for cross-validation correspond to individual study sites, in this case, hospitals. This method is useful for reporting prediction performance, as allocation to validation folds is nonrandom, is reproducible, and provides better estimates of model generalizability [19]. It also means that every validation fold is an entirely separate hospital from those used for fitting the model, meaning that it provides a better estimate of generalizable performance to new settings than internal validation (where data from a single site are divided into a development dataset and an evaluation dataset). Table 2 shows a summary of how the internal-external cross-validation folds are constructed from the dataset comprising 5 hospitals. In this process, 5 models are fit—one for each hospital being used as a validation dataset, with the remaining 4 hospitals being used to fit the model. Table 2 shows a summary of the data used to fit each of the 5 models (fold 1 to 5) and the final model (fit using the data from all hospitals). Using the entire dataset for model development but internal-external cross-validation for reporting ensures that the model does not suffer from avoidable losses in performance by using a subset of the data and potentially omitting useful observations [16].

Model performance was evaluated in terms of discrimination and calibration. The area under the time-dependent receiver operating characteristic (ROC) curve was calculated for each fold and a “combined” curve was derived to show the moving average across folds. The time-dependent ROC curve was evaluated using prediction times of every 24 hours from 1 day post admission to 1 week after admission. The survivalROC R package [20] was used to estimate the time-dependent ROC curves. To estimate the parameter uncertainty of our performance metrics, we applied the approach described by LeDell and colleagues [21] by modifying the cvAUC R package [22] to combine area under the curve (AUC) values from multiple cross-validation folds and construct 95% CIs. For visualization purposes, the curves labeled “combined” on the ROC plots in this study were constructed using vertical averaging [23].

To evaluate model calibration, we followed advice from McLernon and colleagues [24] to obtain moderate calibration curves over the full-time range. As described in the supplement of their study, we used a Poisson model with the outcome regressed against the linear predictors, fit with a restricted cubic spline, and the log of the predicted cumulative hazard as an offset. Smoothed curves were applied to investigate whether observed fall event rates were equal to the predicted risk among patients with the same predicted risk. We created curves separately for each validation fold.

We considered a model with a possible 30 parameters, a fall rate of 0.0054, a conservative estimate of 0.143 Cox-Snell R² for sample size calculation, a selected prediction time of 24 hours, and mean follow-up of 37.5 hours (the mean length of stay within the truncated dataset). The Cox-Snell R² was obtained by using an estimated C-index of 0.7, the approximate performance of recently published inpatient fall prediction models [8], and transforming it into an estimate for the Cox-Snell R² using methods described in related works by Jinks et al [25] and Royston and Sauerbrei [26]. Using the pmsampsize R package [27], we estimated that the minimum sample size for fitting this Cox model would be 1619 patients, with 330 fall events, and 11.75 events per parameter. During our internal-external cross-validation, the least powered model that was fit included 2977 falls and 135 events per parameter (22 parameters; a summary of model parameters, patients, and events used for each model fit during cross-validation is illustrated in Multimedia Appendix 2). The final model was fit with 6004 falls and 231 events per parameter (26 parameters).

Approval was granted by the Metro South Human Research Ethics Committee (HREC/2020/QMS/64807). The requirement for informed consent was waived.

To protect the privacy and confidentiality of candidates, the data used in the study cannot be publicly shared. However, analysis code, including software and package versions, is available on GitHub [28].

Using only routinely recorded information in EMR systems from 5 hospitals of varying capacities and capabilities, we successfully developed a novel prognostic model for inpatient falls using a time-to-event modeling approach. Our approach has the potential to add value to inpatient falls prevention efforts by updating prognostic risk predictions with new information throughout a patient’s admission as it becomes available without requiring clinical staff to collect additional fall-risk screening or assessment data and manually estimate risks. It was encouraging to observe useful predictive performance that was not inferior to more labor- and data-intensive fall risk prediction approaches that featured in our previous review [8]. Studies reporting AUC within our previous review that included external validation and used data from adult inpatients typically reported AUC in the range of 0.7 to 0.85 [30-32]. By developing our model using only data from EMR, this labor-efficient approach to estimation of inpatient fall risk has the potential to be implemented in clinical decision support systems that aid in appropriate allocation of effective inpatient fall prevention interventions without dependency on physical functional performance assessment [33,34].

Overall, our model demonstrates similar or greater discriminatory performance than many other published models for inpatient falls, particularly those that have been externally validated [2,8]. Calibration was site dependent but generally poor. The variable model performance, in terms of both discrimination and calibration, and its association with hospitals may be explained by differences in patient cohorts, fall rates, and the nature of care provided in each hospital.

The poor calibration performance may suggest that the absolute model predicted risk should be interpreted with caution. However, predictions are likely to remain useful for if the models were implemented in a manner that prioritized higher-risk patients (based on rank of estimated risk) toward a fall prevention intervention rather than applying a selected cutoff point above which the patient receives the intervention.

This approach would better take advantage of the high discriminatory performance without relying heavily on the need for the model to be well calibrated and for the selected cutoff point to be used on potentially biased estimates of risk. The use of a cutoff point may not work effectively due to the observed miscalibration based on local clinical contexts, but model recalibration by hospital could be beneficial [35].

We observed that there was a high level of between-hospital variability in ROC curves, particularly for prediction times at day 4 and later. We suspect that this is due to there being relatively few patients and fallers who are still admitted at that time. For example, the particularly unstable ROC curve for hospital 4 observed in Figure 1 had a median length of stay of 0.4 (IQR 0.2-2.1) days, with likely relatively few admitted patients, and very few admitted patients with falls, remaining after day 3. We suspect that the ROC curve for this hospital for these later prediction times is likely influenced by relatively few data.

While some aspects of our model are directly comparable to existing inpatient fall prediction models, several are unique. We have not identified any time-to-event prediction models in the literature that evaluate performance in terms of discrimination and calibration [8]. As binary prediction models typically identify fallers as patients who fell over the remaining time in the admission without time-varying covariates, perhaps the closest comparison from our study would be our model’s discriminatory performance at day 1 before time-varying covariates had influenced risk predictions. For this scenario, the discrimination of our model demonstrates superior performance during internal-external cross-validation. In addition, the advantages of a time-to-event approach include the ability to update predictions as the admission progresses and new data become available, and better incorporation of current length of stay into estimated risk.

We used internal-external cross-validation because this is the most appropriate way to evaluate a clinical prediction model when (1) data are available from multiple sites and (2) there are enough data at each site to ensure that models developed from each cross-validation fold are adequately powered [16]. Many clinical prediction model development studies may only rely on data from 1 to 3 sites. In our previous review, we found that only 2 of 51 studies used data from more than 3 hospitals to develop their model (43 of which only used data from a single hospital) [8]. When only a few sites, it may be impractical to conduct internal-external cross-validation as there would only be 2- or 3-folds, and potentially have too small a portion of the entire dataset for those models to be adequately powered. Fortunately, in our study, we had data from 5 sites and, of the 5 models fit during cross-validation, the lowest number of events per parameter was 135, which was much higher than the minimum number estimated by our sample size calculations (11.75).

For evaluation of the discrimination, we used the area under the time-dependent ROC curve. To create a single summary measure, we used vertical averaging of the 5-folds [23]. A potential limitation of this is that hospitals are equally influential on the “combined” curve and the final estimate of discrimination, there is no additional weight given to hospitals with more data than others. However, in this study, the models with the lowest performance were typically those evaluated on hospitals with fewer admitted patients and events, rather than the largest hospital (“Fold 1” curve in Figure 1), which was consistently the best performing.

Most existing approaches to predicting inpatient falls use data collected upon admission to estimate risk for the duration of the entire admission. This may lead to confusing interpretations; if clinicians interpret predicted risks as a function of admission time, this can be misleading. For example, a user may interpret the patient’s daily fall risk as the total estimated risk divided by the expected length of stay in days, which inappropriately assumes falls risk is evenly distributed over the entire admission. A patient may have a relatively high risk of falling on any given day, but if they are admitted for a condition that is associated with shorter lengths of stay, then a binary prediction model may predict a deceptively low risk. Time-to-event models avoid this potentially erroneous conclusion by accounting for length of stay inherently. We anticipate that a model of this level of discriminatory performance may be useful when allocating patients for fall prevention intervention, but this would need to be evaluated prospectively or possibly within a simulation study [36].

The use of machine learning methods for clinical prediction models has been increasing in popularity, including neural networks or tree-based ensemble models such as random forests or extreme gradient boosting [37,38]. There are variants of these approaches capable of handling time-to-event models but none, to our knowledge, have been applied to develop a prognostic model for falls. For example, random survival forests have been used to predict colorectal cancer prognosis [39], and “pycox” [40], an implementation of a Cox model for a neural network using PyTorch, has been used to predict in-hospital events including mortality and discharge [41]. There has also been a recent increase in popularity of sequence models, including attention models and recurrent neural networks, which can incorporate a sequence of input data of variable lengths and of irregular intervals [42]. We expect that combining these approaches with a Cox layer, as is used in “pycox,” may lead to the benefits of both being able to use time-varying covariates and being able to appropriately predict falls as a time-to-event outcome. A recent study has implemented a similar approach to fit a time-dependent Cox survival neural network to generate dynamic predictions, with example datasets including age-related eye disease and time-to-liver transplant [43]. The machine learning methods included in our recent review [8] were limited to those predicting falls as either occurring any time in the admission or within a fixed time horizon [44,45], but none used a time-to-event outcome. Consequently, it is difficult to know whether falls may be better estimated by the additional complexity and ability to estimate interaction effects within these methods.

We expect that by incorporating more time-varying predictors in a Cox model to predict falls may improve performance but may require a larger dataset with more falls to avoid overfitting. Another limitation of using this model is that it requires computation to estimate risk and cannot be easily implemented without an EMR, as other commonly used analog falls risk assessment tools can [33,46]. Although a reduction in data collection burden for health professionals is a strength of our model, it is a potential barrier to implementation in settings without the required levels of digital maturity.

Conclusions

In this study, we have presented a prognostic model using time-to-event modeling. We used internal-external cross-validation with data from 5 hospitals with notable heterogeneity to estimate its generalizable performance and present the final model that incorporated data from all 1,107,556 inpatient admissions. We anticipate that the reported model estimates may be adjusted to improve calibration performance at each hospital, and that this model has potential use for improving clinical care in the context of prompting the initiation of fall prevention strategies. However, evaluating our model in a simulated clinical environment or prospective clinical trial is a requisite step in determining whether the model is clinically useful. These studies may guide the implementation of the presented model within a decision support system for falls prevention strategies in hospital inpatient environments. Given the use of time-to-event modeling approaches, we also anticipate that machine learning methods that account for time-varying predictors may lead to improved discrimination and calibration.