This study involved the development and external validation of a clinical prognostic model. The reporting follows the Transparent Reporting of a multivariable prediction model for Individual Prognosis Or Diagnosis (TRIPOD) statement (Collins et al., 2015). Data for model development and internal validation were sourced from the prospective PPMI cohort. PPMI is a large, multicenter, longitudinal observational study initiated by The Michael J. Fox Foundation to identify PD progression biomarkers by enrolling drug-naive, early PD participants ¹ (Parkinson Progression Marker Initiative, 2011). The PPMI study was approved by the institutional review board (IRB) at each participating site, and all participants provided written informed consent. Detailed inclusion and exclusion criteria for PPMI are described elsewhere (Parkinson Progression Marker Initiative, 2011). PPMI data (RRID:SCR_006431) were accessed by us in April 2025. For this study, we included de novo PD participants from the PPMI cohort who met the following criteria: (1) no history of falls at their baseline visit, based on fall questionnaire responses; (2) Hoehn and Yahr (H&Y) stage 1–2 at baseline; and (3) at least 36 months of follow-up data or follow-up until the occurrence of a first fall within this period. Participants with missing data for key variables or with comorbid conditions significantly affecting gait and balance (e.g., major stroke sequelae, peripheral neuropathy, myopathy, vertigo, visual impairment) were excluded from this analysis.

External validation data were sourced from an ongoing, single-center, prospective PD cohort initiated in 2022 at the Department of Neurology, Yueyang Hospital of Integrated Traditional Chinese and Western Medicine, Shanghai University of Traditional Chinese Medicine. Inclusion criteria were: (1) diagnosis of idiopathic PD according to the Movement Disorder Society (MDS) clinical diagnostic criteria (Postuma et al., 2015); (2) age 40–80 years; (3) H&Y stage 1–2 at enrollment; and (4) no reported PD-related falls before enrollment. Exclusion criteria were: (1) secondary parkinsonism (e.g., vascular, drug-induced, traumatic, normal pressure hydrocephalus, tumor-related) or parkinsonism-plus syndromes (e.g., progressive supranuclear palsy, multiple system atrophy); (2) comorbid conditions significantly affecting gait and balance; (3) prior deep brain stimulation therapy; and (4) severe dementia or psychiatric disorders impairing assessment or follow-up. This local study was approved by the Ethics Committee of Yueyang Hospital of Integrated Traditional Chinese and Western Medicine, Shanghai University of Traditional Chinese Medicine (IRB No. 2022-106), and all procedures were conducted in accordance with the principles of the Declaration of Helsinki. All participants provided written informed consent.

Baseline clinical data were collected at enrollment for both cohorts. For the PPMI cohort, demographic and medical history information was collected at baseline, including age, sex, years of education, body mass index (BMI), family history of parkinsonism, dominant symptom side, disease duration, and H&Y stage. (1) Motor symptom assessment: The severity of the four cardinal motor features of PD was assessed by summing relevant sub-items from the MDS-Unified Parkinson’s Disease Rating Scale (MDS-UPDRS) Parts II and III during the “off” medication state: bradykinesia score (items 3.2, 3.4–3.9, 3.14), tremor score (items 2.10, 3.15–3.18), rigidity score (item 3.3), and postural instability and gait disorder (PIGD) score (items 2.12–2.13, 3.10–3.12) (Goetz et al., 2008; Wolke et al., 2023). Patients were classified as tremor-dominant (TD) or non-TD based on the TD/PIGD mean ratio (Stebbins et al., 2013). (2) Non-motor symptom assessment: The 13 sub-items of MDS-UPDRS Part I were used to provide a general assessment of overall PD non-motor symptoms (Goetz et al., 2008). Given the very low scores for non-motor symptoms in early PD, which could lead to floor effects if raw scores were used directly, each sub-item was operationalized as a binary variable (score ≥ 1 indicating symptom presence, 0 indicating absence). The total score of the Scale for Outcomes in Parkinson’s disease for Autonomic Symptoms (SCOPA-AUT) was used to evaluate the level of autonomic dysfunction (Visser et al., 2004). Due to the skewed distribution of Geriatric Depression Scale-15 (GDS-15) scores in early PD, a score > 5 was considered indicative of at least mild depressive mood (Shin et al., 2019). Anxiety levels were assessed using the State-Trait Anxiety Inventory (STAI) and its subscales (Oei et al., 1990). Global cognitive function was evaluated using the adjusted Montreal Cognitive Assessment (MoCA) score (Nasreddine et al., 2005; Lu et al., 2011). The presence of hyposmia was assessed using the revised University of Pennsylvania Smell Identification Test (UPSIT), with a UPSIT percentile ≤ 15% considered indicative of hyposmia (Brumm et al., 2023). Sleep-wake dysfunction was assessed using the Epworth Sleepiness Scale (ESS) and the REM Sleep Behavior Disorder Screening Questionnaire (RBDSQ) (Kumar, 2003; Stiasny-Kolster et al., 2007). Asymptomatic orthostatic hypotension (OH) was defined as a drop in systolic blood pressure of ≥ 20 mmHg or diastolic blood pressure of ≥ 10 mmHg within 1–3 min of standing relative to supine values, provided that the MDS-UPDRS Part I Item 1.12 (“Lightheadedness on Standing”) score was 0. This combination confirmed the presence of orthostatic hypotension in the absence of typical symptoms (Beach and McKay, 2024). Further details on these variables have been described in previous publications (Parkinson Progression Marker Initiative, 2011).

To ensure comparability between cohorts, the local cohort followed operational procedures consistent with PPMI for the collection of demographic information and primary clinical data. However, a few differences existed. Firstly, for olfactory function assessment, the local cohort used subjective reporting rather than the objective UPSIT used in PPMI. Secondly, unlike the PPMI cohort where participants were drug-naive at baseline, most patients in the local cohort had initiated levodopa replacement therapy at enrollment; therefore, their baseline levodopa equivalent daily dose (LEDD) was recorded, calculated according to established methods (Tomlinson et al., 2010; Jost et al., 2023). All motor assessments in both cohorts were conducted during a practically defined “off” medication state to minimize the influence of anti-PD medications on assessment results.

The primary outcome of this study was the time (in months) from baseline visit to the first fall experienced by PD patients within a 36-month follow-up period. A fall was defined according to widely accepted criteria: an unexpected event in which the participant comes to rest on the ground, floor, or other lower level (van der Marck et al., 2014; Montero-Odasso et al., 2022). This study focused on the first fall event, excluding near falls and not differentiating subsequent recurrent falls. In the PPMI cohort, fall information was collected using a standardized questionnaire at each regular follow-up visit (typically every 6 months) (Parkinson Progression Marker Initiative, 2011). In the local cohort, fall data were collected prospectively through structured interviews during routine clinical follow-ups (typically every 3 months). The time origin for the study was the date of the baseline visit for the PPMI cohort and the date of the first completed assessment for the local cohort. All participants included in the final analysis were successfully followed for 36 months or until their first fall, with no deaths or loss to follow-up for other reasons. Therefore, standard survival analysis methods were employed, and specific models for competing risks were not required. For participants who did not experience a first fall during the follow-up period, their data were right-censored at 36 months, and their follow-up time was recorded as 36 months.

To address potential bias from missing data in the PPMI cohort, multiple imputation was performed using the “mice” R package (version 3.17.0). Approximately 86% of participants had complete data for all candidate predictors. For variables with missing values, logistic regression imputation (“logreg”) was used for dominant side (missingness ∼12%) and UPSIT (∼3%), while predictive mean matching (“pmm”) was used for disease duration (∼11%). Five complete imputed datasets (m = 5) were generated with 10 iterations each. Iteration trace plots and density plots of imputed data showed good convergence and plausible distributions of imputed values, indicating satisfactory imputation quality (Supplementary Figures 1A–C). Subsequent statistical analyses were performed on these five imputed datasets, and results were pooled using Rubin’s rules.

Continuous variables were described as mean ± standard deviation (SD) or median [interquartile range (IQR)] based on their distribution, while categorical variables were presented as counts (n) and percentages (%). Between-group comparisons were performed using independent samples t-tests, Mann-Whitney U tests, chi-squared tests, or Fisher’s exact tests, as appropriate. Each of the five imputed PPMI datasets was randomly split into a training set (70%) and an internal validation set (30%). In each imputed training set, univariable Cox proportional hazards regression analysis was performed for all candidate predictors. Variables with a p < 0.05 in at least three of the five imputed datasets were included in the candidate pool for the multivariable model. These selected variables were then entered into a multivariable Cox proportional hazards regression model. A backward stepwise selection strategy based on the Akaike Information Criterion (AIC) was used to identify the most parsimonious and statistically significant combination of independent predictors. The Bayesian Information Criterion (BIC) was also calculated to assist in model comparison. Hazard ratios (HRs) and their 95% confidence intervals (CIs) were reported for each factor. A nomogram was constructed based on the final multivariable Cox model to serve as a visual tool for clinicians to estimate the 36-month probability of a first fall in fall-naive, early PD patients. The model’s overall performance was assessed in the training set, internal validation set, and external validation cohort based on discrimination, calibration, and clinical utility at 18 and 36 months. Overall discrimination was evaluated using Harrell’s C-index. Time-dependent C-indices and time-dependent areas under the receiver operating characteristic curve (AUCs) were calculated at 18 and 36 months to assess predictive accuracy at specific time points. An optimism-corrected C-index was calculated for internal validation via bootstrapping (500 resamples). Calibration was assessed by generating calibration curves at 18 and 36 months. Decision curve analysis (DCA) was used to evaluate the net benefit of the model across a range of threshold probabilities at 18 and 36 months. To further assess the model’s risk stratification capacity, patients were categorized into high-risk and low-risk groups based on the median of their predicted risk scores. Kaplan-Meier survival curves for time to first fall were plotted and compared between these groups using the log-rank test. All statistical analyses were performed using R software (version 4.5.0; The R Foundation for Statistical Computing, Vienna, Austria)² and relevant packages. A two-sided p < 0.05 was considered statistically significant.

A total of 403 patients from the PPMI cohort and 150 patients from the local cohort were included in this study; the screening process is detailed in Figure 1. Patients in the PPMI cohort were randomly divided into a training set (n = 283) and an internal validation set (n = 120). No statistically significant differences were observed between these two sets for any baseline characteristics (all P > 0.05) (Supplementary Table 1). Significant differences in multiple baseline characteristics were observed between the total PPMI cohort and the local external validation cohort (Table 1). Compared to the PPMI cohort, PD patients in the local cohort were older, had fewer years of education, lower BMI, and a lower proportion of males. Regarding clinical features, patients in the local cohort had a longer disease duration, a higher proportion of H&Y stage 2 and non-tremor-dominant phenotype, and a lower proportion with a family history of parkinsonism. Motor function assessments revealed that patients in the local cohort had higher MDS-UPDRS Part II and Part III scores. For non-motor symptoms, patients in the local cohort had higher MDS-UPDRS Part I scores and SCOPA-AUT scores, but lower MoCA scores. Concurrently, the prevalence of GDS score > 5 and asymptomatic OH was also higher in the local cohort.

During the 36-month follow-up period, 68 (16.9%) patients in the PPMI cohort experienced a first fall. The Kaplan-Meier estimated cumulative incidence of first falls within 3 years in the PPMI cohort was 16.9% (95% CI: 13.1–20.5%), with a median time to first fall of 23.0 months (IQR: 14.0–31.0 months). Similarly, within the 36-month follow-up period, 31 (20.7%) patients in the local cohort experienced a first fall. Their 3-year cumulative incidence of first falls was 20.7% (95% CI: 13.9–26.9%), with a median time to first fall of 19.0 months (IQR: 13.5–28.0 months). Although the 3-year cumulative fall incidence was slightly higher in the local cohort than in the PPMI cohort, this difference did not reach statistical significance (log-rank P = 0.27) (Figure 2).

In the PPMI training set, univariable Cox proportional hazards regression analysis was initially performed on all candidate demographic and clinical variables to screen for potential predictors (Supplementary Table 2). A total of 13 variables were significantly associated with the risk of a first fall in the univariable Cox analysis and were subsequently considered for the multivariable Cox proportional hazards regression model. The 13 variables found to be significant in univariable Cox analysis were subsequently entered into a multivariable Cox proportional hazards regression model. After initial fitting and assessment of this multivariable model, 7 candidate predictors were retained for further selection procedures (Table 2). “Disease duration” did not reach the predefined statistical significance level (P = 0.164) and was thus excluded first. Among the remaining 6 statistically significant predictors, “Hallucinations and psychosis” although statistically associated with fall risk, had an extremely low prevalence in early PD patients, accounting for only 3.18% in our training cohort. We compared the overall performance of Model 1 (including this variable) and Model 2 (excluding this variable), with results presented in Supplementary Table 3. Model 1 showed slightly better AIC, Bayesian Information Criterion (BIC), and C-index in the training set, but Model 2 had a slightly higher C-index in the internal validation set (0.768 vs. 0.763), suggesting that Model 1 might be slightly overfitted. Due to the limited improvement in predictive performance and its potential to diminish model robustness, “Hallucinations and psychosis,” was removed from the model. Ultimately, our clinical prediction model incorporated 5 independent predictors: lower BMI, presence of asymptomatic OH, lower MoCA score, GDS score > 5, and higher PIGD score (Figure 3). The variance inflation factor for all variables in the final model was less than 5, indicating no multicollinearity. Based on the Schoenfeld residuals test, all covariates in the model met the proportional hazards assumption (all P > 0.05), and the global test for the proportional hazards assumption yielded P = 0.068. For ease of clinical application, the final model was visualized as a nomogram (Figure 4). By locating the corresponding values on their respective axes to obtain points, summing these points to get a total score, which allows for the estimation of an individual’s 18 and 36-month risk of a first fall from the corresponding risk prediction axes.

In the PPMI training set, the C-index for the prediction model was 0.847. Internal validation using 500 bootstrap resamples estimated an optimism of 0.003, resulting in an optimism-corrected C-index of 0.844 (95% CI: 0.786–0.907). In the PPMI internal validation set, the model’s C-index was 0.768 (95% CI: 0.656–0.880). In the local external validation cohort, the C-index was 0.825 (95% CI: 0.768–0.882). Overall and time-specific discriminative performance metrics, including C-indices and AUCs for all datasets, are summarized in Table 3. To further evaluate the model’s discrimination at specific time points, time-dependent ROC curves were plotted. At 18 months, the AUCs for the model in the PPMI training set (Figure 5A), internal validation set (Figure 5B), and external validation set (Figure 5C) were 0.870 (95% CI: 0.806–0.934), 0.794 (95% CI: 0.659–0.929), and 0.787 (95% CI: 0.699–0.876), respectively. At 36 months, the corresponding AUCs were 0.932 (95% CI: 0.881–0.983), 0.909 (95% CI: 0.820–0.998), and 0.884 (95% CI: 0.829–0.939) in the PPMI training set, internal validation set, and external validation set, respectively. These results indicate that the model demonstrated good to excellent discrimination across the training, internal validation, and external validation cohorts, particularly for predicting longer-term first fall risk.

Calibration curves for predicting 18 and 36-month first fall risk were generated for the training set, internal validation set, and external validation cohort, with confidence intervals estimated using 500 bootstrap resamples. For 18-month first fall risk prediction, the calibration curve for the PPMI training set (Figure 6A) showed good agreement across most risk strata, with slight deviation only at very high predicted risks ( > 0.75). The internal validation set (Figure 6B) showed acceptable calibration in the high-risk range. The 18-month calibration curve for the external validation cohort (Figure 6C) exhibited some fluctuation at predicted risks above 0.7. For 36-month first fall risk prediction, the calibration curve in the PPMI training set (Figure 6A) closely followed the ideal diagonal, demonstrating excellent calibration. In the PPMI internal validation set (Figure 6B) and the local external validation cohort (Figure 6C), the calibration curves also showed good agreement with the ideal line in the medium-to-high risk strata, although observed fall risk was slightly lower than predicted in the moderate-risk range. Overall, the calibration curves suggest that the model achieved acceptable calibration for predicting 36-month first fall risk across the three datasets, as predicted probabilities generally aligned with observed fall frequencies. However, for 18-month interim risk prediction, the model’s calibration was less consistent across datasets, with some deviation between predicted and observed values in certain risk strata, particularly in the external validation cohort.

DCA was employed to evaluate the net benefit of using the model to guide clinical decisions across different threshold probabilities. For 18-month first fall risk prediction, the model demonstrated some clinical net benefit in all three datasets. In the PPMI training set (Figure 7A), internal validation set (Figure 7B), and external validation cohort (Figure 7C), using the model was superior to the “treat all” or “treat none” strategies within threshold probability ranges of approximately 0.05–0.65, 0.05–0.45, and 0.05–0.35, respectively, with corresponding maximum net benefits of approximately 0.07, 0.08, and 0.06. Compared to 18-month predictions, the model showed superior and more robust clinical utility for predicting 36-month first fall risk. In the PPMI training set (Figure 7D), the model provided a net benefit over the two extreme strategies across a very wide range of threshold probabilities (0.05–0.95), with a maximum net benefit of about 0.16. In the PPMI internal validation set (Figure 7E), the model was beneficial within a threshold range of approximately 0.05–0.74, achieving a maximum net benefit of about 0.18. In the local external validation cohort (Figure 7F), the model demonstrated clinical net benefit within a threshold range of approximately 0.05–0.75, also with a maximum value of about 0.18. Overall, the DCA results indicate that this prediction model, particularly for 36-month first fall risk, provides a net benefit across a wide range of clinical decision thresholds in all three datasets, suggesting good clinical utility. However, for 18-month interim predictions, the range and magnitude of net benefit provided by the model were relatively limited.

To further evaluate the model’s risk stratification ability, patients in each dataset were divided into high-risk and low-risk groups based on the median of the 36-month predicted first fall risk scores generated by the model. Kaplan-Meier curves were then used to compare the cumulative fall-free rates between the two groups. In the PPMI training set (Figure 8A), the 36-month fall-free survival probability was significantly lower in the high-risk group compared to the low-risk group (P < 0.001). In the PPMI internal validation set (Figure 8B), the model also effectively distinguished patients into different risk strata, with the fall-free survival curve for the high-risk group being significantly lower than that for the low-risk group (P = 0.005). In the local external validation cohort (Figure 8C), the model’s risk stratification capacity was further confirmed, as patients in the high-risk group had a significantly higher risk of experiencing a first fall than those in the low-risk group (P < 0.001). These results demonstrate that the fall-free survival probability was consistently and significantly lower in the high-risk group across the training, internal validation, and external validation sets.