Development of a parsimonious predictive model for height gain outcome in children with central precocious puberty following gonadotropin-releasing hormone agonist therapy

Introduction

Background

Central precocious puberty (CPP) is a pediatric endocrine condition defined by premature activation of the hypothalamic-pituitary-gonadal (HPG) axis, leading to the appearance of secondary sexual characteristics before the age of 8 years in girls and 9 years in boys (1,2). Although the condition can occur in both sexes, it is substantially more common in females, with male cases accounting for fewer than 10% of all diagnoses. Epidemiological surveys in East Asian populations have documented a rising incidence over recent decades, with one large Korean registry-based study estimating a prevalence of approximately 55.9 per 100,000 among girls (3). A more recent Korean population-based analysis covering 2008 to 2020 reported that CPP incidence increased by 15.9 times in girls and 83.3 times in boys over the study period (4). A 2025 systematic review and meta-analysis synthesizing global data further confirmed these upward trends, reporting a pooled prevalence of 1.26% in girls, with notable geographic heterogeneity across regions (5). Similar upward trends have been reported in European and North American cohorts, prompting growing clinical and public health attention. It should be acknowledged, however, that these apparent increases may partly reflect a broader secular trend toward earlier timing of pubertal onset in the general population, together with heightened clinical awareness and more frequent referral and testing, rather than a true rise in pathological precocious puberty alone. This distinction matters clinically: some children who meet diagnostic criteria may not require treatment, especially when the expected height benefit is limited.

The clinical consequences of CPP extend well beyond the visible emergence of pubertal features at an inappropriately young age. Psychosocially, affected children may experience significant emotional distress, body image difficulties, and social challenges arising from physical development that is markedly discordant with their chronological age and cognitive maturity (4). From a somatic standpoint, one of the primary concerns of premature HPG axis activation is accelerated skeletal maturation. Supraphysiological concentrations of sex steroids stimulate chondrocyte proliferation in the growth plates, transiently increasing linear growth velocity but simultaneously advancing epiphyseal fusion. The clinical impact of this process varies with the age of onset and the tempo of pubertal progression. In severe cases, particularly those with earlier onset and more rapidly progressive puberty, this sequence can lead to early cessation of longitudinal bone growth and a final adult height that falls substantially below genetic potential. In contrast, children with later or more slowly progressive presentations (for example, girls in whom puberty begins closer to 7.5–8 years of age) may experience a less pronounced effect on final height (6,7).

Gonadotropin-releasing hormone agonists (GnRHa) have been the mainstay of CPP pharmacotherapy since their introduction in the early 1980s, and their efficacy and safety profile are now supported by more than four decades of clinical experience (2,7,8). The pharmacological rationale rests on the paradoxical effect of continuous, non-pulsatile GnRH receptor stimulation: after an initial brief flare of gonadotropin release, sustained receptor occupation results in desensitization and downregulation of pituitary GnRH receptors, effectively suppressing the premature secretion of luteinizing hormone (LH) and follicle-stimulating hormone (FSH). When treatment is initiated sufficiently early, while substantial growth potential remains, the resulting reduction in gonadal steroid output decelerates skeletal maturation, slows epiphyseal fusion, and preserves the window for continued linear growth; conversely, when therapy begins late, with little residual growth potential, the achievable benefit is correspondingly limited (8). Multiple studies have documented that appropriately timed GnRHa treatment can augment adult height by 1.1 to 10.5 cm relative to pretreatment predicted adult height (PAH) (9-11), although the magnitude of benefit varies widely among individual patients.

Rationale and knowledge gap

Despite the well-established efficacy of GnRHa therapy at a population level, a recurrent observation across published cohorts is the substantial heterogeneity in individual height outcomes (12). Some children achieve or even exceed their mid-parental target height, while others gain considerably less than expected during treatment. This variability has prompted extensive investigation into potential prognostic factors. Prior studies have linked final or near-adult height to a range of baseline variables, including chronological age at treatment initiation, bone age, height standard deviation score (SDS), body mass index (BMI), stimulated LH peak, pubertal (Tanner) stage at treatment initiation (11), and genetic target height (12-14). However, existing predictive frameworks have generally incorporated numerous covariates, and the resulting models, while statistically robust, can be cumbersome to apply in routine outpatient settings where rapid, point-of-care risk assessment is most needed.

An additional gap in the literature concerns the role of insulin-like growth factor-1 (IGF-1) as a prognostic biomarker in this context. IGF-1 is the principal downstream effector of growth hormone action on the skeletal growth plate, mediating both the endocrine and paracrine signals that drive chondrocyte proliferation and hypertrophy (13). Recent reviews have further underscored that IGF-1 produced locally from chondrocytes is especially important for chondrocyte differentiation, proliferation, and hypertrophy, and that it stimulates extracellular matrix production and ossification within the growth plate (15). GnRHa administration, by suppressing gonadal steroid production, can secondarily influence IGF-1 secretion and bioavailability, potentially modulating the pace of linear growth during treatment (8,16). Although individual studies have reported associations between IGF-1 dynamics and growth velocity changes on GnRHa, the prognostic value of baseline IGF-1 SDS as a component of a parsimonious predictive tool has not been systematically examined.

Objective

The present study was therefore designed with two principal objectives. First, we sought to identify independent predictors of suboptimal height gain during GnRHa therapy by analyzing a comprehensive panel of baseline clinical and biochemical variables in a well-characterized cohort of children with idiopathic CPP. Second, we aimed to construct and internally validate a parsimonious predictive model that retains clinically meaningful discriminative power while minimizing the number of input variables required, thereby enhancing practical applicability at the bedside. We hypothesized that a streamlined model incorporating IGF-1 SDS and bone age advancement would achieve predictive accuracy comparable to more complex multi-variable models, and that such a tool could facilitate early identification of patients who may more frequent on-treatment growth monitoring (for example, more closely spaced auxological and growth-velocity assessments during GnRHa therapy) or consideration of adjunctive growth-promoting interventions such as recombinant human growth hormone (rhGH). We present this article in accordance with the TRIPOD reporting checklist (available at https://tp.amegroups.com/article/view/10.21037/tp-2026-0254/rc).

Methods

Study design and patient population

The screening, eligibility assessment, cohort allocation, and data analysis are summarized in Figure 1. This retrospective cohort study enrolled children diagnosed with idiopathic CPP who received triptorelin therapy at the Pediatric Department of Yueqing Maternal and Child Health Hospital from February 2020 to November 2025. The diagnostic criteria for CPP were based on the 2009 international consensus (2) and included: (I) onset of breast development (Tanner stage B2 or above) before 8 years of age in girls; (II) GnRH stimulation test demonstrating a peak LH concentration ≥5 IU/L; and (III) bone age advancement of at least 1 year relative to chronological age. At our institution, the GnRH stimulation test is performed routinely as part of the standardized diagnostic work-up for all children presenting with suspected CPP; consequently, a stimulated LH peak was available for every patient in this cohort. We recognize that the 2009 international consensus also accepts a diagnosis based on a sufficiently elevated basal LH concentration without a stimulation test, and that many contemporary practices adopt this approach. Because the GnRH-stimulated LH peak was uniformly available in our cohort, we used it as the hormonal diagnostic criterion to ensure diagnostic consistency across all enrolled patients. Inclusion criteria comprised: (I) confirmed diagnosis of idiopathic CPP meeting all three diagnostic criteria; (II) completion of a standardized triptorelin treatment regimen (intramuscular injection of triptorelin at standard doses administered every 28 days); and (III) availability of complete baseline and follow-up clinical data. Exclusion criteria included: (I) CPP attributable to organic etiologies such as central nervous system tumors, hamartomas, or other intracranial lesions as determined by cranial MRI; (II) concurrent endocrinopathies known to affect linear growth, including hypothyroidism, growth hormone deficiency, Turner syndrome, or adrenal disorders; (III) prior administration of rhGH or other growth-modifying therapies before enrollment; and (IV) incomplete medical records precluding accurate outcome ascertainment.

Figure 1 Flow diagram of the study. Children with idiopathic CPP treated with triptorelin (GnRHa) were screened for eligibility. Inclusion criteria required confirmed CPP diagnosis, completed treatment course, and availability of full follow-up data. Exclusions included organic CPP, other endocrinopathies, and prior GH use. A total of 163 patients were enrolled. Baseline variables included age, height SDS, BMI SDS, Tanner stage, bone age advancement (BA-CA), LH peak, E2, IGF-1 SDS, and MPH. Height gain was dichotomized at the cohort median into favorable (n=113, 69.3%) and unfavorable (n=50, 30.7%) outcomes. Univariate screening (P<0.10) followed by multivariable logistic regression identified BMI SDS (OR =2.33), IGF-1 SDS (OR =1.94), and BA-CA (OR =21.49) as independent predictors. Three models were developed: a comprehensive model (M0; AUC =0.838), a univariate model (M1; AUC =0.664), and a parsimonious model (M2; AUC =0.815). Model performance was evaluated using ROC analysis, calibration, DCA, and NRI/IDI with bootstrap validation (B=1,000). The parsimonious model M2 (IGF-1 SDS + BA-CA) demonstrated comparable performance to M0 (dAUC =−0.023, P=0.24; calibration slope =0.957). AUC, area under the curve; BA-CA, bone age minus chronological age; BMI, body mass index; CPP, central precocious puberty; dAUC; DCA, decision curve analysis; E2, estradiol; GH, growth hormone; GnRHa, gonadotropin-releasing hormone agonist; IDI, integrated discrimination improvement; IGF-1, insulin-like growth factor-1; LH, luteinizing hormone; MPH, mid-parental height; NRI, net reclassification improvement; OR, odds ratio; ROC, receiver operating characteristic; SDS, standard deviation score.

This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. The study was approved by the Ethics Committee of Yueqing Maternal and Child Health Hospital (approval No. LY2025006). Given the retrospective nature of the study, the requirement for individual informed consent was waived by the Ethics Committee.

Data acquisition

Baseline clinical and laboratory parameters were systematically extracted from medical records at the time of initial presentation and prior to commencement of GnRHa therapy. The following data were collected: (I) anthropometric measurements including chronological age (years), standing height (cm), and body weight (kg); (II) pubertal staging according to the Tanner classification system, performed by experienced pediatric endocrinologists; (III) bone age assessment from a left hand-wrist radiograph interpreted by the Greulich-Pyle atlas method, with all radiographs independently read by experienced radiologists at our institution; (IV) hormonal evaluations comprising the GnRH stimulation test with measurement of LH peak, basal serum estradiol (E2), and serum IGF-1 concentrations; and (V) genetic parameters including parental heights and calculated mid-parental height (MPH), derived as (father’s height + mother’s height – 13)/2 for girls. BMI, calculated as body weight (kg) divided by the square of standing height (m), and the corresponding height SDS, BMI SDS, and IGF-1 SDS were derived d using age- and sex-specific national reference standards. Bone age advancement [bone age minus chronological age (BA-CA)] was defined as the arithmetic difference between bone age (years) and chronological age (years).

Outcome definition and patient stratification

The primary outcome measure was height gain during the treatment period, operationally defined as the difference between height at the conclusion of GnRHa therapy (or at the last available follow-up visit during active treatment) and baseline height at treatment initiation. Patients were dichotomized into two groups using the cohort median height gain as the threshold: the favorable outcome cohort (y=0, n=113) comprised patients whose height gain was at or above the median, while the unfavorable outcome cohort (y=1, n=50) comprised patients whose height gain fell below the median value. The cohort median was selected as the dichotomization threshold because no universally accepted, validated cut-off for “suboptimal” treatment-period height gain currently exists for children with CPP; a data-driven, distribution-based threshold therefore provides an objective and reproducible partition that does not depend on arbitrary external values. The resulting groups were not exactly equal (favorable n=113, unfavorable n=50) because a cluster of patients shared height-gain values at or immediately adjacent to the median, and these tied values were assigned to the at-or-above (favorable) group; this produced an approximately 2:1 split rather than a strict 50:50 division. We recognize that dichotomization at the median is inherently a pragmatic simplification and that the choice of threshold influences group sizes and effect estimates; the rationale for, and limitations of, this approach—together with alternative outcome definitions such as the change in PAH—are discussed further in the Discussion (Strengths and limitations).

Statistical analysis

All statistical analyses were performed using Python (version 3.9) with the scikit-learn, statsmodels, and scipy packages. Continuous variables were summarized as mean ± standard deviation (SD) and compared between groups using the independent-samples t-test for normally distributed data or the Mann-Whitney U test for non-normally distributed data, as determined by the Shapiro-Wilk test. Categorical variables were expressed as frequencies and percentages and compared using Pearson’s Chi-squared test or Fisher’s exact test as appropriate.

Candidate predictors were initially screened through univariate logistic regression. Variables achieving a significance level of P<0.10 in univariate analysis were entered into multivariable logistic regression using a forward stepwise approach. The candidate variables were not selected by an automated algorithm alone; they were pre-specified on the basis of established clinical and biological plausibility as recognized determinants of growth response in CPP (chronological age, BMI SDS, IGF-1 SDS, bone age advancement, and Tanner stage), and the data-driven screening served only to confirm and rank their contribution. Because stepwise selection can be unstable in modest datasets, we limited the candidate predictors and assessed optimism by bootstrap validation. To mitigate these concerns, the number of candidate predictors was deliberately kept small relative to the number of events, the comprehensive model (M0) was pre-specified rather than allowed to expand indefinitely, and the optimism attributable to model development was quantified directly through bootstrap internal validation (see below). Penalized regression approaches [for example, least absolute shrinkage and selection operator (LASSO)] represent a complementary strategy and are identified as a direction for future work in larger, externally validated cohorts. Three predictive models were constructed for systematic comparison:

• M0 (comprehensive model): y ~ IGF-1 SDS + BMI SDS + BA-CA + Tanner stage + age;
• M1 (univariate model): y ~ IGF-1 SDS;
• M2 (parsimonious model): y ~ IGF-1 SDS + BA-CA.

The rationale for these three model structures was as follows. M0 served as the reference comprehensive model, incorporating all clinically pre-specified candidate predictors that survived univariate screening. M1 was constructed with IGF-1 SDS as the single predictor because IGF-1 SDS was the novel biomarker of primary interest in this study; M1 therefore represents the minimal “biomarker-only” benchmark against which the incremental value of adding skeletal maturation could be judged. M2 was designed to be the most parsimonious model that retained strong discriminative performance: among the independent predictors identified in M0, IGF-1 SDS and bone age advancement (BA-CA) are the two that are routinely available at the initial diagnostic work-up of every child with CPP and that together carry the majority of the prognostic signal (BA-CA being the strongest single predictor and IGF-1 SDS the biomarker of interest). In the model-building process we examined alternative two-variable combinations using the same candidate set (for example, BA-CA paired with BMI SDS, age, or Tanner stage, and IGF-1 SDS paired with BMI SDS or age). Combinations built around BA-CA performed similarly to one another in discrimination, reflecting the dominant contribution of skeletal maturation; however, the IGF-1 SDS + BA-CA pairing was selected as M2 because it combined the strongest skeletal-maturation predictor with the study’s biomarker of interest while requiring no investigation beyond the standard CPP work-up, and because adding further variables to this pair yielded negligible improvement in discrimination at the cost of greater complexity.

Model discrimination was quantified using the area under the receiver operating characteristic curve (AUC), with 95% confidence intervals (CIs) derived via DeLong’s method (17). Pairwise comparison of AUC values was performed using the DeLong test. Model calibration was evaluated by the Hosmer-Lemeshow goodness-of-fit test and by calibration plots comparing predicted probabilities against observed event rates across deciles of predicted risk. Reclassification performance was assessed using the continuous net reclassification improvement (NRI) and integrated discrimination improvement (IDI) (18). Decision curve analysis (DCA) was employed to evaluate clinical utility by estimating net benefit across a range of clinically relevant risk thresholds (1–50%) (19). Internal validation was performed using bootstrap resampling with 1000 iterations to compute optimism-corrected AUC, calibration slope, and Brier score. All statistical tests were two-tailed, and P<0.05 was considered statistically significant.

Results

Patient characteristics

A total of 163 children with idiopathic CPP who met all eligibility criteria were enrolled in this study. Of these, 113 patients (69.3%) were classified into the favorable outcome cohort and 50 patients (30.7%) into the unfavorable outcome cohort. Comparative analysis of baseline clinical and biochemical characteristics between the two groups is presented in Table 1. Children in the unfavorable outcome group were significantly older at diagnosis (7.64±0.74 vs. 7.37±0.74 years, P=0.03), had higher BMI SDS (1.21±0.90 vs. 0.63±0.58, P=0.003), higher IGF-1 SDS (1.94±0.92 vs. 1.43±0.90, P=0.001), and greater bone age advancement (2.76±0.34 vs. 2.46±0.25 years, P<0.001) compared with those in the favorable outcome group. Baseline height SDS, stimulated LH peak, serum estradiol concentrations, and Tanner stage distribution did not differ significantly between groups.

Univariate and multivariable logistic regression analysis

Univariate logistic regression analysis identified chronological age, BMI SDS, IGF-1 SDS, and bone age advancement (BA-CA) as significantly associated with unfavorable height gain outcome (all P<0.05). Stimulated LH peak, baseline height SDS, serum estradiol concentration, and MPH were not significantly associated with outcome in univariate analysis.

Multivariable logistic regression incorporating these candidate predictors (Model M0) revealed three independently significant predictors of unfavorable outcome: BMI SDS (OR =2.33, 95% CI: 1.14–4.76, P=0.02), IGF-1 SDS (OR =1.94, 95% CI: 1.24–3.02, P=0.004), and bone age advancement (OR =21.49, 95% CI: 5.41–85.30, P<0.001). Chronological age and Tanner stage were not independently associated with outcome after covariate adjustment (both P>0.05). The regression coefficients for all three models are presented in Tables 2,3.

Model performance comparison

Discriminative performance of the three models is illustrated in Figure 2 and summarized in Table 4. The comprehensive model (M0) achieved an AUC of 0.838 (95% CI: 0.774–0.903), the univariate model (M1) an AUC of 0.664 (95% CI: 0.573–0.754), and the parsimonious model (M2) an AUC of 0.815 (95% CI: 0.748–0.883). DeLong test demonstrated no statistically significant difference in discriminative performance between M2 and M0 (ΔAUC =–0.023, P=0.24). Hosmer-Lemeshow goodness-of-fit testing indicated adequate calibration for all three models (M0: χ²=13.78, P=0.09; M1: χ²=11.34, P=0.18; M2: χ²=11.13, P=0.19).

Figure 2 Receiver operating characteristic curves for the predictive models. ROC curves illustrating discriminative performance of three predictive models for unfavorable height outcome. M0: comprehensive model incorporating five baseline predictors (IGF-1 SDS, BMI SDS, bone age advancement, Tanner stage, and chronological age), AUC =0.838 (95% CI: 0.774–0.903); M1: univariate model with IGF-1 SDS exclusively, AUC =0.664 (95% CI: 0.573–0.754); M2: parsimonious model comprising IGF-1 SDS and bone age advancement, AUC =0.815 (95% CI: 0.748–0.883). Diagonal reference line denotes random classification (AUC =0.5). DeLong test demonstrated no statistically significant difference between M2 and M0 (ΔAUC =−0.023, P=0.24). AUC, area under the curve; BMI, body mass index; CI, confidence interval; IGF-1, insulin-like growth factor-1; ROC, receiver operating characteristic; SDS, standard deviation score.

Reclassification analysis showed that, relative to M0, M2 yielded a continuous NRI of –0.439 (95% CI: –0.764 to –0.114, P=0.008) and an IDI of –0.058 (95% CI: –0.165 to 0.050, P=0.29). The negative NRI indicates that although overall model discrimination is comparable, the parsimonious model may provide somewhat less refined individual-level risk reclassification. DCA (Figure 3) revealed comparable clinical net benefit for M0 and M2 across the 1–50% risk threshold range, with both models substantially outperforming M1 and the default strategies of treating all patients or treating none.

Figure 3 Decision curve analysis of the three predictive models. The curves depict the net clinical benefit of each model (M0, comprehensive; M1, univariate; M2, parsimonious) across a range of threshold probabilities for predicting an unfavorable height-gain outcome. The “Treat all” line represents the default strategy of intervening in all patients, and the “Treat none” line the strategy of intervening in none. Across the clinically relevant 1–50% threshold range, M0 and M2 yielded comparable net benefit, both exceeding M1 and the two default strategies, supporting the clinical utility of the parsimonious model M2 for risk-stratified decision-making.

Bootstrap internal validation

Internal validation via bootstrap resampling (B=1,000) was conducted for all three models, and the results are summarized in Table 5. Optimism-corrected AUC values were 0.810, 0.661, and 0.807 for M0, M1, and M2, respectively. The parsimonious model M2 exhibited the lowest optimism coefficient (0.008), suggesting minimal susceptibility to overfitting relative to M0 (optimism =0.028). Assessment of calibration slopes revealed that M2 (0.957, 95% CI: 0.642–1.389) most closely approximated the ideal value of unity, outperforming M0 (0.834) and M1 (1.121). Brier scores were 0.160 for M0, 0.201 for M1, and 0.167 for M2, indicating comparable overall predictive accuracy between the comprehensive and parsimonious models. Calibration plots are presented in Figure 4.

Figure 4 Calibration plots derived from bootstrap internal validation (B=1,000 iterations) for M0 (comprehensive model), M1 (univariate model), and M2 (parsimonious model). The diagonal reference line represents ideal calibration (predicted probability equals observed probability). Apparent calibration curve (solid line) reflects model fit to derivation data; bias-corrected calibration curve (dashed line) estimates expected performance in new patient populations. M2 demonstrated optimal calibration slope (0.957, 95% CI: 0.642–1.389) closest to the ideal value of unity, superior to M0 (calibration slope =0.834) and M1 (calibration slope =1.121), indicating minimal overfitting and robust agreement between predicted and observed event probabilities. CI, confidence interval.

Taken together, the parsimonious model M2 incorporating IGF-1 SDS and BA-CA demonstrates comparable discrimination, calibration, and clinical decision utility relative to the comprehensive five-predictor model, while exhibiting greater parsimony and practical clinical applicability.

Discussion

Key findings

The present study retrospectively analyzed clinical data from 163 children with idiopathic CPP treated with triptorelin to identify prognostic determinants of treatment-period height gain and to develop a streamlined predictive model suitable for clinical application. The central finding is that a parsimonious two-variable model (M2) incorporating IGF-1 SDS and bone age advancement achieves discriminative performance essentially equivalent to the comprehensive five-predictor model (M0), with AUC values of 0.815 vs. 0.838 (P=0.24 by DeLong test). Bootstrap internal validation further confirmed the robustness of M2, demonstrating the lowest optimism coefficient among the three candidate models (0.008 vs. 0.028 for M0), the calibration slope closest to unity (0.957 vs. 0.834), and a Brier score comparable to M0 (0.167 vs. 0.160). DCA revealed that M2 provides net clinical benefit similar to M0 across the range of risk thresholds most relevant to clinical practice (1–50%), reinforcing its potential utility as a pragmatic bedside prediction tool.

Strengths and limitations

Several strengths of this study merit consideration. First, the systematic head-to-head comparison of models ranging from a single predictor to a five-variable comprehensive model allows transparent evaluation of the trade-off between model complexity and predictive performance. Second, the use of multiple complementary performance metrics—including discrimination (AUC), calibration (Hosmer-Lemeshow test, calibration plots and slopes), reclassification (NRI, IDI), clinical utility (DCA), and internal validation (bootstrap)—provides a thorough assessment that goes beyond simple AUC comparison. Third, the two predictors retained in M2 are routinely measured during the initial diagnostic evaluation of every child with suspected CPP, meaning no additional laboratory or imaging investigations are required to operationalize the model.

However, several methodological limitations must be acknowledged. First, the retrospective single-center design introduces inherent risks of selection bias and information bias, as data were extracted from medical records not originally collected for research purposes. Second, the sample size (n=163) is modest, particularly for the unfavorable outcome group (n=50), which may limit the stability of parameter estimates and reduce the power to detect smaller effect sizes. Third, while bootstrap internal validation provides a reasonable estimate of model generalizability within the source population, it cannot substitute for external validation in geographically and demographically distinct cohorts; multicenter prospective validation remains an essential next step, and the model should be regarded as developmental and not yet ready for clinical deployment until such validation is completed. Fourth, the primary outcome was height gain during the treatment period rather than final adult height, which is the end point of greatest clinical relevance. Height gain during therapy is a practical surrogate, but its relationship to final adult stature may be influenced by subsequent growth after GnRHa discontinuation, the timing of puberty resumption, and other factors not captured in this analysis. Fifth, an important limitation of absolute treatment-period height gain as the outcome is that the residual height still achievable depends strongly on pubertal stage and bone age at treatment initiation: a child with more advanced maturation has, by definition, less height left to gain, so the outcome may partly reflect baseline maturational status rather than the true response to GnRHa therapy. We chose this endpoint because it was uniformly and objectively recorded for all patients and is clinically familiar. More informative endpoints would be the change in PAH over the treatment interval, a height change adjusted for chronological age, pubertal stage, and bone age advancement, or—definitively—final (near-)adult height after completion of growth; these were not uniformly available in our retrospective dataset but should be the basis of future prospective work. This same caveat bears on the very large odds ratio (OR) observed for bone age advancement (OR =21.49 in M0; OR =26.65 in M2): although biologically consistent with the dominant role of skeletal maturity, the wide confidence interval likely also reflects the modest sample size and the partial dependence of the outcome on baseline maturation, and the estimate should therefore be interpreted with caution pending external validation. Sixth, the study enrolled predominantly girls with idiopathic CPP; findings may not be directly generalizable to boys or to patients with organic etiologies.

Comparison with similar research

The predictive significance of bone age advancement in determining height outcomes after GnRHa therapy is well supported by the existing literature. Trujillo and colleagues (12) analyzed 100 girls with idiopathic CPP and identified bone age at treatment initiation as a key determinant of final height outcome, consistent with the present observation that bone age advancement was the strongest single predictor (OR =21.49). Tseng et al. (9) similarly found that PAH at baseline—itself heavily influenced by bone age—was strongly associated with final height, and concluded that patients with the greatest bone age advancement and the lowest predicted heights derived the most benefit from prolonged GnRHa therapy. Knific and colleagues (6), in a systematic review and meta-analysis of final adult height data, confirmed that degree of skeletal maturation at diagnosis consistently predicts height outcomes across diverse study populations and treatment protocols. Mondkar et al. (11), in a longitudinal assessment of Indian girls with idiopathic CPP, further demonstrated that lesser skeletal maturity and lower baseline estradiol at treatment initiation translated into better adult height outcomes. Lerman et al. (20) likewise reported, in a cohort of girls with CPP and early fast puberty, that baseline auxological and skeletal-maturation parameters were important determinants of the height response to GnRHa therapy, reinforcing the prognostic relevance of skeletal maturity. , The role of IGF-1 as a prognostic variable in CPP has received comparatively less systematic attention. Park et al. (16) reported a significant positive correlation between changes in IGF-1 SDS and height SDS during GnRHa therapy (r=0.405, P<0.001), providing evidence that the IGF-1 axis is mechanistically relevant to growth kinetics on treatment. However, that study focused on longitudinal IGF-1 changes rather than the predictive value of baseline IGF-1 concentrations. The present investigation extends this work by demonstrating that a single baseline IGF-1 SDS measurement, when combined with bone age advancement, contributes sufficient prognostic information to match the discriminative performance of models incorporating three to five additional covariates.

It is also instructive to compare the AUC performance of our models with those reported in other CPP prediction studies. Although direct comparisons are complicated by differences in outcome definitions, patient populations, and statistical methodologies, the AUC of 0.815 for M2 compares favorably with published models in pediatric endocrinology, which typically report AUC values in the range of 0.70–0.85 for clinical prediction tools. The added value of M2 lies not in achieving the highest possible discrimination but in delivering strong performance with the smallest number of readily available inputs.

Explanations of findings

The biological plausibility of the two retained predictors deserves careful consideration. Bone age advancement serves as an integrated surrogate marker of cumulative sex steroid exposure on the skeleton. A greater discrepancy between bone age and chronological age at the time of CPP diagnosis reflects more prolonged or intense premature HPG axis activation, which translates into a more advanced state of epiphyseal maturation and, consequently, less residual growth potential (9,13). Even when GnRHa therapy effectively suppresses further pubertal progression, the irreversible component of prior skeletal maturation constrains the absolute height gain achievable during treatment. This principle is consistent with the large OR observed for BA-CA in our multivariable model (OR =21.49 in M0; OR =26.65 in M2), underscoring the outsized importance of skeletal maturity status at the time of treatment initiation.

IGF-1 SDS, the second component of M2, provides mechanistically complementary information. IGF-1 is the principal downstream mediator of growth hormone action on the growth plate (15). Elevated baseline IGF-1 SDS in the context of CPP likely reflects augmented somatotropic axis activity, which may be driven by sex steroid stimulation of GH secretion and hepatic IGF-1 production. When GnRHa therapy suppresses gonadal steroid output, the associated decline in this stimulatory signal may lead to a proportionally greater reduction in IGF-1 levels and growth velocity in patients who started with the highest baseline concentrations (10,18). Thus, paradoxically, children with higher baseline IGF-1 SDS may be at increased risk of deceleration in growth rate during treatment—an observation that aligns with our finding that elevated IGF-1 SDS was independently associated with unfavorable outcome. We emphasize that growth velocity was not directly measured or modelled in the present study; the link between IGF-1 dynamics and growth velocity is invoked here only as a biological mechanism to interpret our findings, drawing on prior longitudinal work (16). Because our retrospective dataset did not contain serial measurements at uniform intervals suitable for reliable growth-velocity computation, we deliberately refrained from reporting growth velocity to avoid presenting unstable estimates. The absence of directly assessed growth velocity is therefore acknowledged as a limitation, and prospective evaluation of on-treatment growth-velocity trajectories in relation to baseline IGF-1 SDS is identified as an important avenue for future research.

An additional consideration is the exclusion of BMI SDS from the parsimonious model despite its significance in the M0 multivariable analysis (OR =2.33, P=0.02). This apparent discrepancy is likely explained by the well-established biological interconnections among obesity, bone maturation, and IGF-1 status in childhood. Pediatric obesity is associated with accelerated skeletal maturation and advanced bone age (21), and obese children tend to exhibit higher circulating IGF-1 concentrations driven by hyperinsulinemia and enhanced hepatic GH receptor sensitivity (22). Consequently, simultaneous inclusion of IGF-1 SDS and BA-CA in M2 effectively captures much of the prognostic information carried by BMI SDS, rendering its explicit inclusion statistically redundant without meaningful loss of predictive accuracy. These same biological interconnections suggest that, in clinical settings where IGF-1 is not routinely measured, a model combining BMI SDS and bone age advancement could serve as a practical surrogate, since BMI SDS appears to carry part of the prognostic information conveyed by IGF-1 SDS. Because BMI and bone age are universally available at minimal cost, such a model would be more broadly applicable than one requiring IGF-1 measurement, and we regard the formal development and head-to-head comparison of an IGF-1-free (BMI + BA-CA) model as a valuable extension of the present work. We have refrained from reporting performance figures for this alternative model here, as a rigorously developed and internally validated BMI + BA-CA model warrants dedicated analysis in an adequately powered, ideally prospective and externally validated cohort to avoid the data-dredging and overfitting that post hoc reporting of an additional model in the present modest sample could introduce.

Implications and actions needed

From a clinical perspective, the principal value of the M2 model lies in its simplicity and immediate applicability. Both IGF-1 SDS and bone age advancement are routinely ascertained during the standard diagnostic workup for CPP, at our institution, meaning that the model can be applied at the time of initial diagnosis without any additional investigations or costs in this setting. We acknowledge, however, that serum IGF-1 measurement is not a routine component of the CPP evaluation in all centers or health systems, and that applying the IGF-1-based model elsewhere would entail additional, potentially costly laboratory testing. This limits the immediate generalizability of M2 to centers that do not already measure IGF-1, and reinforces the value of developing an IGF-1-free alternative (for example, the BMI + BA-CA combination discussed above) for broader applicability. Clinicians could use the predicted probability generated by M2 to stratify patients into risk categories at the point of treatment initiation. Such risk stratification could inform the treatment decision in more than one direction. In addition to flagging high-risk children for closer monitoring, the model could contribute to the shared decision of whether GnRHa therapy is worthwhile at all: because not every child derives a clinically meaningful height benefit from treatment, families and clinicians might reasonably elect to defer or forgo therapy—and avoid its associated cost and burden—when the predicted height response is poor and the primary indication is height preservation rather than psychosocial concern. We note, however, that this study was not designed to evaluate treatment-versus-no-treatment decisions, and any such use would require validation against final-height outcomes in untreated and treated comparison groups. Children identified as high risk for suboptimal height gain could then be targeted for more frequent growth monitoring during GnRHa therapy, with e-discussion of adjunctive interventions in carefully selected cases where clinically indicated. We are careful not to overstate this point: the present data do not establish that adding rhGH on the basis of this model improves final adult height, and our findings should not, on their own, be taken as justification for routinely intensifying therapy. Any role for adjunctive rhGH guided by such risk stratification would need to be tested directly in prospective trials with adult-height endpoints (10).

It is important to emphasize that the M2 model is intended as a clinical decision-support tool rather than a deterministic rule. The negative NRI observed for M2 relative to M0 (–0.439, P=0.008) indicates that the parsimonious model, while comparable in overall discrimination and calibration, may provide less refined individual-level risk stratification in certain probability strata. Clinicians should therefore integrate model-derived predictions with clinical judgment, patient and family preferences, and evolving growth trajectory data obtained during follow-up.

Looking ahead, several research priorities emerge. External validation in independent, multi-institutional cohorts of children with CPP from diverse geographic and ethnic backgrounds is the most immediate need. Prospective studies incorporating longitudinal follow-up to final adult height would permit assessment of whether treatment-period height gain—the outcome modeled here—serves as a reliable surrogate for the ultimate clinical endpoint. Additionally, investigation of whether dynamic changes in IGF-1 SDS during therapy can further refine risk prediction beyond what is achievable with baseline values alone represents a promising avenue for model enhancement. Finally, the development of user-friendly digital tools (e.g., nomograms, web-based calculators, or mobile applications) integrating the M2 model could facilitate its adoption in routine clinical practice. In the interim, and to make the model immediately usable at the bedside without specialized statistical software, the M2 model can be applied as an explicit formula derived from the regression coefficients in Table 2. The predicted probability of an unfavorable height-gain outcome is given by P=1 / (1 + e−L), where L = −10.249 + 0.573 × (IGF-1 SDS) + 3.283 × (BA-CA in years). As a worked example, for a child with an IGF-1 SDS of 2.0 and a bone age advancement of 2.8 years, L = −10.249 + (0.573 × 2.0) + (3.283 × 2.8) ≈0.087, yielding a predicted probability of an unfavorable outcome of approximately 1 / (1 + e−0.087) ≈0.52, i.e., roughly a 52% predicted risk. This calculation requires only a basic calculator and the two routinely available inputs, allowing the model to be operationalized in clinic before dedicated digital tools become available.

Conclusions

This study developed and internally validated a parsimonious predictive model incorporating IGF-1 SDS and bone age advancement for prognostication of treatment-period height gain outcome in children with CPP receiving GnRHa therapy. The two-variable model achieved discriminative performance comparable to a comprehensive five-predictor model (AUC 0.815 vs. 0.838, P=0.24), with favorable calibration characteristics (calibration slope 0.957) and minimal overfitting susceptibility (optimism coefficient 0.008) on bootstrap internal validation. This model enables risk stratification at initial presentation to identify patients with potential suboptimal height gain during therapy, thereby supporting more intensive on-treatment growth monitoring in selected patients. The present data do not, however, establish that this model can guide decisions about adjunctive interventions such as rhGH. Multicenter external validation and prospective studies with final adult height endpoints are warranted to confirm the generalizability and long-term prognostic value of this model before widespread clinical implementation.

Acknowledgments

None.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tp.amegroups.com/article/view/10.21037/tp-2026-0254/rc

Data Sharing Statement: Available at https://tp.amegroups.com/article/view/10.21037/tp-2026-0254/dss

Peer Review File: Available at https://tp.amegroups.com/article/view/10.21037/tp-2026-0254/prf

Funding: The research was supported by The Wenzhou Medical and Health Research Project (No. 2025027).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tp.amegroups.com/article/view/10.21037/tp-2026-0254/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. The study was approved by the Ethics Committee of Yueqing Maternal and Child Health Hospital (approval No. LY2025006). Given the retrospective nature of the study, the requirement for individual informed consent was waived by the Ethics Committee.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Srilanchakon K, Supornsilchai V, Wacharasindhu S, et al. Precocious puberty: a comprehensive review of diagnosis and clinical presentation, etiology, and treatment. Asian Biomed (Res Rev News) 2025;19:69-77. [Crossref] [PubMed]
2. Carel JC, Eugster EA, Rogol A, et al. Consensus statement on the use of gonadotropin-releasing hormone analogs in children. Pediatrics 2009;123:e752-62. [Crossref] [PubMed]
3. Kim YJ, Kwon A, Jung MK, et al. Incidence and Prevalence of Central Precocious Puberty in Korea: An Epidemiologic Study Based on a National Database. J Pediatr 2019;208:221-8. [Crossref] [PubMed]
4. Kang S, Park MJ, Kim JM, et al. Ongoing increasing trends in central precocious puberty incidence among Korean boys and girls from 2008 to 2020. PLoS One 2023;18:e0283510. [Crossref] [PubMed]
5. Zhang X, Xu Y, Yan L, et al. Global prevalence and incidence of precocious puberty: a systematic review and meta-analysis. BMC Public Health 2025;26:99. [Crossref] [PubMed]
6. Knific T, Lazarevič M, Žibert J, et al. Final adult height in children with central precocious puberty - a retrospective study. Front Endocrinol (Lausanne) 2022;13:1008474. [Crossref] [PubMed]
7. Zevin EL, Eugster EA. Central precocious puberty: a review of diagnosis, treatment, and outcomes. Lancet Child Adolesc Health 2023;7:886-96. [Crossref] [PubMed]
8. Bangalore Krishna K, Fuqua JS, Rogol AD, et al. Use of Gonadotropin-Releasing Hormone Analogs in Children: Update by an International Consortium. Horm Res Paediatr 2019;91:357-72. [Crossref] [PubMed]
9. Tseng CH, Lee YJ, Huang CY, et al. The effects of gonadotropin-releasing hormone agonist on final adult height among girls with early and fast puberty. Front Endocrinol (Lausanne) 2023;14:1271395. [Crossref] [PubMed]
10. Yang C, Liu Z, Zeng L, et al. Pharmacotherapy for children with central precocious puberty or early puberty: A systematic review and meta-analysis. Medicine (Baltimore) 2025;104:e41936. [Crossref] [PubMed]
11. Mondkar S, Khadilkar V, Yewale S, et al. Longitudinal assessment of auxological parameters, adult height outcome and its determinants in leuprolide-treated Indian girls with idiopathic central precocious puberty. J Pediatr Endocrinol Metab 2024;37:62-8. [Crossref] [PubMed]
12. Trujillo MV, Lee PA, Reifschneider K, et al. Using change in predicted adult height during GnRH agonist treatment for individualized treatment decisions in girls with central precocious puberty. J Pediatr Endocrinol Metab 2023;36:299-308. [Crossref] [PubMed]
13. Lee PA, Neely EK, Fuqua J, et al. Efficacy of Leuprolide Acetate 1-Month Depot for Central Precocious Puberty (CPP): Growth Outcomes During a Prospective, Longitudinal Study. Int J Pediatr Endocrinol 2011;2011:7. [Crossref] [PubMed]
14. Park JH, Hwang IT, Yang S. Relationship between growth velocity and change of levels of insulin-like growth factor-1, insulin-like growth factor binding protein-3 and, IGFBP-3 promoter polymorphism during GnRH agonist treatment. Ann Pediatr Endocrinol Metab 2020;25:234-9. [Crossref] [PubMed]
15. Cho JH, Jung HW, Shim KS. Growth plate closure and therapeutic interventions. Clin Exp Pediatr 2024;67:553-9. [Crossref] [PubMed]
16. White CA, Serrat MA. Dysregulation of insulin-like growth factor-1 signaling in postnatal bone elongation. Biochem Cell Biol 2023;101:388-93. [Crossref] [PubMed]
17. DeLong ER, DeLong DM, Clarke-Pearson DL. Comparing the areas under two or more correlated receiver operating characteristic curves: a nonparametric approach. Biometrics 1988;44:837-45.
18. Pencina MJ, D'Agostino RB Sr, D'Agostino RB Jr, et al. Evaluating the added predictive ability of a new marker: from area under the ROC curve to reclassification and beyond. Stat Med 2008;27:157-72; discussion 207-12. [Crossref] [PubMed]
19. Vickers AJ, Elkin EB. Decision curve analysis: a novel method for evaluating prediction models. Med Decis Making 2006;26:565-74. [Crossref] [PubMed]
20. Lerman L, Yackobovitch-Gavan M, Phillip M, et al. Gonadotropin-releasing hormone analogs treatment in girls with central precocious puberty and early fast puberty. Pediatr Res 2024;95:1051-9. [Crossref] [PubMed]
21. Kythreotis AK, Nicolaou M, Mitsinga E, et al. The Interplay Between Body Weight and the Onset of Puberty. Children (Basel) 2025;12:679. [Crossref] [PubMed]
22. Reinehr T, Roth CL. Is there a causal relationship between obesity and puberty? Lancet Child Adolesc Health 2019;3:44-54. [Crossref] [PubMed]