The REHIPED registry2 is a Spanish, voluntary, multicenter registry conducted under the aegis of the Spanish Society of Pediatric Cardiology and Congenital Heart Disease (SECPCC) and includes patients older than two months and younger than 18 years with PH. The registry started in 2010 and includes prospective patients (recruited in 2010 and thereafter) and retrospective patients (diagnosed between 1998 and 2009). The registry includes 24 Spanish hospitals and, since 2016, 4 Colombian centers (5% of the registry patients).

The inclusion criterion was an established diagnosis of PH. The diagnosis of PH was generally obtained by right heart catheterization, requiring a mean pulmonary artery pressure ≥ 25mmHg (before the World Symposium on Pulmonary Hypertension of 2024) and> 20mmHg and indexed pulmonary vascular resistance ≥ 3 WUm2 (after 2024). When catheterization was not performed (eg, in an unstable patient), the patient was included based on transthoracic echocardiography or a pathology sample, after evaluation by the registry committee. Written parental informed consent was obtained in all cases. The protocol was approved by institutional review boards and ethics committees in all participating centers.

For the implementation of this risk score, patients in groups 1 (PAH) and 3 (PH associated with lung disease or hypoxia) were included, as these 2 etiopathogenic groups represent most PH cases in pediatrics. A patient inclusion flowchart is shown in figure 1.

Figure 1.

Patient inclusion flowchart.

(0.2MB).

The risk score model was trained with the information provided by the REVEAL Lite 2 variables (heart rate, systolic blood pressure, functional class, renal function, 6-minute walk test distance, natriuretic peptides), with the addition of sex, age at diagnosis, weight, and weight-for-height. The model was trained using variables collected at the time of diagnosis; consequently, all patients were treatment naïve.

All clinical variables (except functional class and age at diagnosis) were normalized according to each patient's age, sex, weight, and height, using reference values from healthy populations.10 The variable “natriuretic peptides” corresponds to the z score for N-terminal pro-B-type natriuretic peptide or B-type natriuretic peptide, normalized to the respective age-specific normal value,11,12 as one or the other was available depending on the center.

Renal function was evaluated using creatinine because glomerular filtration rate was not included in the REHIPED registry. It was also normalized according to age and weight.13,14

Sex, age at diagnosis, weight, and weight-for-height were included in the model because they are required for normalization of the REVEAL variables. As these variables must be collected for normalization purposes, their inclusion was not expected to affect the future usability of the score. We therefore decided to incorporate the information they provide into the model.

For each of the patients, a survival time Ts relative to the diagnosis date, tdiagnosis, was defined as follows. The event was either death of lung transplantation. Ts was truncated at 7 years to avoid noise from learning longer-term survival.

• •

  Ts∈[7years, +∞), right-censored, if tlast known−tdiagnosis>7years, where tlast known was the last known follow up or event date.

• •

  Ts=tevent−tdiagnosis, for patients with a known event at date tevent.

• •

  Ts ∈[tlast follow up− tdiagnosis,+∞), right censored, for patients with no events, using the last known follow up date tlast follow  up .

A gradient boosting training algorithm was implemented to estimate each patient's survival time. To address the limited amount of data, Monte Carlo cross-validation was performed to tune the model and estimate performance on unseen data (figure 2). Once the model was tuned, a final training on the whole population generated the score. Risk stratification provided the cutpoints to classify new patients into 3 risk strata.

Figure 2.

Machine learning process overview. ROC AUC, receiver operating characteristic area under the curve.

(0.58MB).

To obtain a model for the survival time, Ts, a gradient boosting machine learning process15 was implemented using the XGBoost library, with the accelerated failure time log-likelihood as objective. To generate an interpretable model, the base learners were restricted to depth-1 decision trees (gradient-boosted stumps: single decisions with a cutoff on 1 of the available variables), whose effects can be added after the model is computed, to obtain a short decision list.

Only the number of boosting iterations B was tuned. The learning objective, learning rate parameters, and tree depth were fixed.

As the available registry data were scarce, a 250-fold Monte Carlo cross-validation (MCCV) was used to estimate the generalization error and optimize the number of boosting iterations, B, which is the only tuned model parameter. MCCV16 was previously found to be superior to leave-one-out cross-validation17 and has been used in studies with scarce data.18

The area under the receiver operating characteristic curve was used as the validation metric, evaluating the capability to predict 1-year event-free survival (1-year ROC AUC), using only the validation dataset.

A repeated random subsampling validation was conducted. On each repetition: a) a new random split of the full available population was performed, with 80% in the training set and 20% in the validation set; b) a new model was trained using the training data set; and c) the model was evaluated every 100 boosting iterations against unseen data in the validation dataset. The 1-year ROC AUC on unseen data AUCb was obtained as a function of the number of boosting iterations b.

The AUCb was then averaged over 250 repetitions of the validation procedure, obtaining the average 1-year ROC AUC on unseen data AUCb¯ as a function of the number of boosting iterations b. The optimal number of boosting iterations B* is the one that maximizes the AUCb¯: B*= argmaxb(AUCb¯).

For the optimal number of boosting iterations, B*, the distribution of AUCB* values obtained in the 250 repetitions represents an estimate of the model performance on unseen data. From this distribution, the 1-year ROC AUC on unseen data and its 95% confidence interval (95%CI) were estimated.

The full available population was finally used for training with the optimal number of boosting iterations, B*. In the resulting model, the effects of coincident decision points were added to shorten the list of decisions, obtaining the final model.

This model represents an estimate of the survival time in years, truncated at 7 years:

A risk stratification analysis was performed, including only the population with known survival or an event within the first year (n=360), with 3 risk cutoffs.

The risk was estimated among patients in a neighborhood of the score value, using Gaussian kernel filters (σ  =  0.18 score points). The risk was stratified in our population according to European guidelines19: a) low risk: 0% to 5% 1-year mortality or transplantation; b) intermediate risk: 5% to 20% 1-year mortality or transplantation; and c) high risk:> 20% 1-year mortality or transplantation.

Statistical analysis was used to describe the baseline characteristics of the population and to analyze survival in the 3 risk groups. Survival analysis was performed using Kaplan-Meier curves, and significant differences were assessed using the log-rank test and Cox proportional hazards regression. All statistical analyses were performed using Stata 15.

A total of 420 children and adolescents were analyzed: 309 patients had group 1 (PAH), and 111 had group 3 PH (associated with lung diseases and/or hypoxia). Their baseline characteristics at diagnosis are shown in table 1.

Among group 1 patients, 204 (66%) had PAH associated with congenital heart disease, 87 (28%) had idiopathic PAH, and 18 (6%) had PAH due to other causes. In addition, 18% of group 1 patients had significant respiratory comorbidities, including bronchopulmonary dysplasia and other developmental lung disorders.

Within the subgroup of patients with PAH associated with congenital heart disease, 15% had Eisenmenger syndrome, 38% had a nonrestrictive shunt, 17% had a restrictive shunt, and 30% had post-repair PH.

Among group 3 patients, 51 (46%) had bronchopulmonary dysplasia, 33 (30%) had other developmental lung disorders, 20 (19%) had interstitial lung disease, and 7 (6%) had PH due to other causes.

The estimated 1-year and 3-year survival rates, free from transplantation, of the entire cohort were 84% (95%CI, 80-97) and 76% (95%CI, 71-80), respectively. Lung transplantation in the first 3 years of follow-up occurred in 12 patients (5%).

After we performed the MCCV analysis with 250 random splits, the average 1-year ROC AUC curve obtained is shown in figure 4, as a function of the number of boosting iterations. The optimal number of boosting iterations was B*= argmaxb(AUCb¯) = 3700, maximizing the average 1-year ROC AUC.

Figure 4.

Average 1-year ROC AUC of the validation set in the 250 random splits used to identify the optimal number of boosting iterations (vertical red line). ROC AUC, receiver operating characteristic area under the curve.

(0.18MB).

The final risk score was obtained after training on the full population using the optimal number of boosting iterations B*=3700.

The risk score model produces a score that estimates survival time in years. This means that the patient is estimated to survive longer for larger score values and thus is inversely proportional to PH risk.

The contribution of each variable to the risk score function is listed in table 2: negative values increase risk, and larger absolute values have a higher impact on the prediction. Table 2 also shows the assigned score when the value is missing.

As can be observed in table 1, the 6-minute walk test was available in only 17% of the patients, mainly because of age limitations. B-type natriuretic peptide or N-terminal pro-B-type natriuretic peptide was available in 48% of the patients. Regarding the rest of the model variables, 8 patients (2%) had exactly 1 missing value, 29 patients (7%) had 2 missing values, and 2 patients (< 1%) had 3 or more missing values.

The variable with the largest impact was functional class. Age at diagnosis (with worse prognosis for younger children), heart rate (which conferred worse prognosis in bradycardic patients), low weight for height, low blood pressure, and natriuretic peptides also had an important impact on survival.

Sex and the normalized (z score) values of creatinine were not assigned a factor at any cutoff, indicating that their impact on survival was lower than that of the other variables and thus they were never selected in the training process.

The C-statistic of the risk score function was 0.82 for 1-year survival when evaluated in the full available population. The ROC curve is shown in figure 5.

Figure 5.

1-year ROC curve of the risk score function in the full population. ROC AUC, receiver operating characteristic area under the curve.

(0.18MB).

The distribution of AUCB* values obtained in the validation set of each of the 250 splits in MCCV (figure 6) represents an estimate of the model performance on unseen data from the study population. From this distribution, the 1-year ROC AUC on unseen data was estimated, as the validation sample on each split was not part of the training population. The mean 1-year ROC AUC obtained was 0.73 (95%CI, 0.60-0.87). The final model was trained with the whole available population using the same number of boosting iterations.

Figure 6.

Distribution of the 1-year ROC AUC on unseen data in the validation set of 250 Monte Carlo cross-validation random splits, for the optimal number of boosting iterations B*=3700. This distribution is an estimate of the performance of the final model on new unseen data. ROC AUC, receiver operating characteristic area under the curve.

(0.2MB).

To calculate the risk for each score value, patients within a neighborhood of the score value were used, obtaining the risk as a function of the score value in figure 7.

Figure 7.

One-year mortality or transplantation risk as a function of the score. Horizontal red lines show the stratification risk limits. Vertical red lines show the score value limits separating the risk strata.

(0.13MB).

The risk stratification in table 3 was obtained by identifying the score values corresponding to the 5% and 20% risk limits. The survival rate for each group is shown in figure 8.

Figure 8.

Survival by risk stratification.

(0.31MB).

There was a significant difference in transplant-free survival between the 3 groups (log-rank P <.001). Hazard ratios for the 3 risk strata at 1-year are shown in table 4. Hazard ratios were obtained with respect to the low-risk category.

Subgroup analyses were performed to assess score performance by age (< 2 vs ≥ 2 years) and etiologic group (group 1 vs group 3). The score effectively stratified risk across all subgroups (log-rank P <.001; figure 9 and figure 10).

Figure 9.

Subgroup analysis by age.

(0.26MB).

Figure 10.

Subgroup analysis by etiologic group.

(0.26MB).

An additional subgroup analysis evaluated score performance in patients with (n=342) vs without (n=78) right heart catheterization. As shown in figure 11, the model failed to discriminate between low- and intermediate-risk categories in patients without right heart catheterization but successfully stratified low-intermediate vs high risk.

Figure 11.

Subgroup analysis by availability of right heart catheterization.

(0.27MB).

In comparison with previously published scores, our model is innovative, adapting the REVEAL Lite 2 adult score through variable normalization and allowing a similar scoring system to be used across all ages of the patient's lifespan. It is also valuable because the model is trained to be used both in groups 1 and 3, and it has been trained and tested in a population including both infants and children with a very wide age range. The exclusive use of noninvasive variables makes it applicable to all clinical settings, including intensive care or low-income countries.

Remarkably, only a small series (n=58) published by Haarman et al.24 exploring PH risk factors is close to ours in age (median age of 6.8, interquartile range [2.2–13.4]).

The other 2 published pediatric scores studied older children and adolescents, so there is less expected physiological variability in the non-invasive variables: Qian et al.25 (n=248) developed a risk score using 3 noninvasive variables at diagnosis, including weight z score>−2, functional class, and right ventricular enlargement, with a C index of 0.647. Griffiths et al.26 studied how the biomarker ST2 could improve prediction of prognosis in PH. The data were obtained from 2 different cohorts, and the REVEAL score was applied only to the cohort of older age (n=182, median 13 years [8-17]). The best outcome discrimination was found in the REVEAL score plus ST2, with a C-statistic of 0.78.

All these 3 previous works included only group 1 PH, which can be a problem given the common overlap of etiologies described in children.

The need for a pediatric risk score applicable to the most prevalent PH groups in pediatrics and applicable to infants, children, and adolescents was addressed in the last World Symposium on Pulmonary Hypertension Pediatric Taskforce report.21 The newly designed score addresses PH groups 1 and 3 and takes infants, children, and adolescents into account for training and testing.

This study has several limitations. First, it is based on data from a voluntary, multicenter registry. Although an effort was made to minimize missing data, this remains an important limitation. In addition, many patients were too young to perform a walking test or too young or unwell to undergo right heart catheterization.

Another limitation is that data for this study were collected exclusively at diagnosis, a decision motivated by the voluntary and multicenter nature of the registry, which does not include a standardized follow-up protocol.

The scarcity of the available data is also a limitation for the development and validation of the score. As a mitigation, MCCV was used to estimate the performance of the model on unseen data in the same population, with a mean AUC of 0.73 (95%CI, 0.60-0.87). However, it must be noted that MCCV does not substitute for external validation, which should ideally be performed in further studies.

An optimal stopping point was also computed to avoid overfitting, but this is a common pitfall of machine learning models that cannot be completely avoided.