Risk factors for the time to development of retinopathy of prematurity in premature infants in Iran: a machine learning approach

Abstract

Background

Retinopathy of prematurity (ROP), is a preventable leading cause of blindness in infants and is a condition in which the immature retina experiences abnormal blood vessel growth. The development of ROP is multifactorial; nevertheless, the risk factors are controversial. This study aimed to identify risk factors of time to development of ROP in Iran.

Methods

This historical cohort study utilized data from the hospital records of all newborns referred to the ROP department of Farabi Hospital (from 2017 to 2021) and the NICU records of infants referred from Mahdieh Hospital to Farabi Hospital. Preterm infants with birth weight (BW) ≤ 2000 g or gestational age (GA) < 34 wk, as well as selected infants with an unstable clinical course, as determined by their pediatricians or neonatologists, with BW > 2000 g or GA ≥ 34 wk. The outcome variable was the time to development of ROP (in weeks). Random survival forest was used to analyze the data.

Results

A total of 338 cases, including 676 eyes, were evaluated. The mean GA and BW of the study group were 31.59 ± 2.39 weeks and 1656.72 ± 453.80 g, respectively. According to the criteria of minimal depth and variable importance, the most significant predictors of the time to development of ROP were duration of ventilation, GA, duration of oxygen supplementation, bilirubin levels, duration of antibiotic administration, duration of Total Parenteral Nutrition (TPN), mother age, birth order, number of surfactant administration, and on time screening. The concordance index for predicting survival of the fitted model was 0.878.

Conclusion

Our findings indicated that the duration of ventilation, GA, duration of oxygen supplementation, bilirubin levels, duration of antibiotic administration, duration of TPN, mother age, birth order, number of surfactant administrations, and on time screening are potential risk factors of prognosis of ROP. The associations between identified risk factors were mostly nonlinear. Therefore, it is recommended to consider the nature of these relationships in managing treatment and designing early interventions.

Introduction

Retinopathy of prematurity (ROP) is a condition that impacts preterm infants especially those with extremely low birth weights, leading to diminished vision and, in severe cases, blindness [1]. ROP arises due to the incomplete development of retinal blood vessels and exclusively impacts premature infants [1]. It is a prevalent and preventable condition that can cause childhood blindness. In 2010, ROP was responsible for an estimated 32,300 cases of blindness and visual impairment globally, with the highest incidence occurring in middle-income countries where neonatal care coverage was expanding [2]. More recently, in 2022, Zhang et al. [3]., used cause-specific vision loss data from the Global Health Data Exchange to estimate that in 2019, ROP led to moderate vision impairment in around 49,100 cases, severe vision impairment in 27,500 cases, and blindness in 25,000 cases.

The enhancement of neonatal healthcare in developing countries, including Iran, has led to improved survival rates among premature infants, consequently resulting in a rise in the occurrence of ROP [4]. Numerous studies worldwide have focused on the incidence of ROP and have identified various risk factors, including GA, BW, oxygen therapy concentration [5], bilirubin levels, sex [6], multiple gestation [7] and intraventricular hemorrhage (IVH) [8]. However, the major risk factors remain controversial.

Time to developing ROP is an outcome that considers interval from birth to developing ROP. Nonetheless, some infants that are loss to follow-up, and their developing status are unknown, resulting in censoring in observations. Investigating the risk factors for the time to developing ROP may provide the clinicians with valuable information. This outcome requires statistical analysis under censoring, known as survival analysis. The typical model for analyzing survival data is the Cox proportional hazards model [9]. However, the fundamental premise of this model relies on the proportional hazards assumption as its key determinant. In real-world scenarios, the explanatory variables may not meet this assumption, or they might exhibit strong collinearity and non-linear or complex relationships that cannot be captured by the traditional statistical models [10]. Detecting interactions, especially those involving multiple variables, poses challenges. This process typically involves exhaustive exploration (examining all possible two-way and three-way interactions, for example) or depends on subjective expertise to narrow down the search [10]. In the few past decades, machine learning (ML) models, including random forests (RF), have been used extensively. ML, is a subset of artificial intelligence, involves developing algorithms and statistical models that enable computers to learn from data and make predictions or decisions. This learning process involves training algorithms on large datasets to identify patterns, relationships, and trends, which can then be used to make predictions or decisions with new data. One variant of ML is the supervised learning, where the algorithm is trained on labeled data, and the input data is paired with the correct output. The algorithm learns to map the input data to the correct output by minimizing the error between its predictions and the true labels.

Random survival forest (RSF) is a supervised machine learning technique introduced by Hemant Ishwaran and Udaya B. Kogalur to analyze survival data [10]. In this method, a bivariate vector consisting of the time to a specific event and the status of individuals is considered as the label. This model handles all the mentioned difficulties of the traditional models automatically using forests.

So far, limited attention have been given to the time to developing ROP, with the majority of studies focusing on the incidence of ROP. Several studies have been aimed to determine the risk factors of developing ROP [11,12,13,14]. However, these studies have treated ROP development as a binary outcome [15,16,17]. Moreover, few studies have examined the time to developing ROP, and those that have typically used traditional statistical methods like log-rank test and Cox regression, which assume only linear relationships. To the best of our knowledge, no study has identified risk factors for the time to developing ROP using machine learning for survival data, which can account for nonlinear and complex relationships among many covariates. Moreover, screening criteria for ROP can vary based on advancements in perinatal care, geographic location, and diverse population groups [18,19,20]. This highlights the need for determining risk factors of ROP in every Neonatal Intensive Care Unit (NICU). Therefore, this study aimed to determine risk factors for the time to developing ROP based on the utilization of an advanced ML algorithm of RSF. The method used can capture non-linear relationships and interactions among features more effectively. This sophisticated approach allows for a more nuanced analysis of the data and may uncover novel risk factors or patterns that were previously overlooked by conventional statistical methods.

Material and methods

Data collection

In this historical cohort study, data from the hospital records of all newborns referred to the ROP department of Farabi Hospital (from 2017 to 2021) and the NICU information of infants referred from Mahdieh Hospital to Farabi Hospital were extracted based on a checklist created from patient records and in accordance with the opinions of specialists. In this regard, preterm infants with BW ≤ 2000 g or GA < 34 wk as well as selected infants with an unstable clinical course according to discretion of pediatricians or neonatologists with BW > 2000 g or GA ≥ 34 wk were included. The infants were examined first between 4 and 6 wk after birth or within the 31st and 33rd weeks postmenstrual age, whichever occurred later.

The predictors at checklist included:

• Maternal Characteristics: Maternal age, GA, multiple gestations, birth order, corticosteroid use during pregnancy, delivery method, infertility history, chorioamnionitis, preeclampsia, antibiotic use pre-delivery, abortion history, intrauterine death history, assisted reproduction methods, diabetes status, surfactant administration, Granulocyte Colony-Stimulating Factor GCSF use.

• Infant Characteristics: Gender, sildenafil administration, positive blood culture, resuscitation at birth, resuscitation procedures, clinical sepsis, thrombocytopenia, Intraventricular Hemorrhage (IVH), anticoagulant use (Enoxaparin), milrinone administration, inotropic drug, diuretic, Apgar scores at 1 and 5 min, BW, duration of invasive/non-invasive ventilation or oxygen therapy, number of surfactant administrations, number of tuberculosis pack injections, duration of antibiotic therapy in infants, duration of Total Parenteral Nutrition (TPN) administration, treatment for Patent Ductus Arteriosus (PDA), and highest recorded bilirubin level.

The outcome variable was considered as time from birth to developing ROP (in weeks). Developing ROP is diagnosed through eye examinations by four experienced retina specialists according to the International Classification of Retinopathy of Prematurity (ICROP 2) [21]. Pupils were fully dilated using Mydrax < 0.5% and phenylephrine < 2.5% eye drops. After topical anesthesia, each infant was examined with an indirect ophthalmoscope using a lid speculum, a scleral depressor and a + 20 D or + 30 D lens.

Random survival forest

RSF was designed specifically for dealing with survival data that contains censored information as a direct extension of random forest. RSF retains the fundamental principles of RF, inheriting all of its crucial characteristics [10]. RF, in general, are composed of numerous decision trees created from random samples with replacement. Each tree is comprised of decision nodes (variables) where classification or splitting occurs. The decision nodes represent questions or conditions about the features (independent variables) of the data. These questions are typically binary, meaning they result in a yes or no answer. For instance, in the present study, delivery method is a binary variable and results in a splitting the tree into two daughter nodes of natural delivery and cesarean type. To decide which variable should be considered at the top of the tree, splitting criteria are used. At each decision node, the algorithm selects the feature and the threshold that best splits the data into two subsets. The goal is to create subsets that are as homogeneous as possible with respect to the target variable (the variable we want to predict). In the context of survival analysis, the splitting of tree nodes is determined by maximizing the differences in survival rates between daughter nodes (new nodes). For each tree, the survival time and the patients' survival status are regarded as response variables. The process of selecting the best split and creating child nodes is performed recursively (called Recursive Partitioning) for each subset until a stopping criterion is met. This criterion could be a maximum tree depth, a minimum number of data points in each node, or other conditions aimed at preventing overfitting. Once the splitting process is completed, the terminal nodes of the tree are called leaf nodes or terminal nodes. These nodes represent the predicted value of the target variable for the observations that fall into that node. Subsequently, to predict the outcome for a new observation, the algorithm traverses the tree from the root node to a leaf node by following the path determined by the values of the features. The predicted value for the new observation is then taken as the value associated with the leaf node it reaches.

For survival forests, an ensemble estimate for the hazard function is generated by calculating the hazard for each sample within a dataset. Summing this ensemble over the observed survival times produces the predicted outcome. In this study, the RSF procedure was performed on the ROP dataset, involving the creation of 1000 trees using a log-rank splitting rule for each run.

Identifying risk factors of developing ROP was conducted using variable importance (VIMP) as a rapidly computable internal measure used for rank variables with greater values indicating a more predictive power of variables [22]. We considered a threshold value of 0.002 for variable selection. An imputation strategy was used to handle missing data based on RSF [22]. The extracted information was analyzed using R software (4.3.1) with randomforestSRC package. Hierarchical cluster analysis was conducted to identify high and low risk groups. Log-rank test was used to compare survival curves.

Results

Descriptive statistics

The information collected included data from the 676 eyes (338 infants) with a mean BW of 1656.72 g (SD: 453.8 g). According to Table 1, the majority of the infants were male (55.3%) and first-born (81.4%). Most pregnancies were single gestations (63%), and most infants were born at 32–35 weeks’ GA (45.3%). A significant proportion of births were by caesarean Sect. (79.9%). Moreover, most of infants had an Apgar score of ≥ 7 at 1 and 5 min (82% and 94% respectively). The mean and standard deviation of other variables were reported in Table 1. Additionally, the mean age of mothers were 29 years (standard deviation: 5.84 years). For more details of the characteristics of the infants see Table 1.

Table 1 Descriptive statistics of maternal and infant characteristics

Figure 1, shows the Kaplan–Meier estimate of the survival distribution function for 676 eyes related to 338 infants with 95% confidence intervals indicating the time to development of ROP along with at risk eyes and the number of events over time. The Y-axis in this figure represents the proportion of the eyes still not developing ROP and the X-axis represents the age of infant in weeks. As seen, the one-month survival rate was around 91% and the three-month survival rate was about 62%. The shape of the survival function also indicates that the short-term hazard of developing ROP is high early in the life of infants, but that risk decreases over time (here, it goes down after 30 weeks). Among all eyes included, 210 eyes developed ROP. Of them, 43 eyes (20.5%) were of stage I, 53 eyes (25.2%) were of stage II, 64 eyes (30.5%) were of stage III, and 50 eyes (23.8%) were of stage IV or higher. The median time to developing ROP was 5 weeks for those with stage I and II, and 6 weeks for those with stage ≥ III. Moreover, 90 eyes (13.31%) were Type 1 ROP (according to the ICROP2 guideline) and received treatment.

Fig. 1

Kaplan–Meier estimate of the survival distribution function for 676 eyes related to 338 infants with 95% confidence intervals indicating time to development of Retinopathy of prematurity (ROP)

Variable selection

The data was analyzed using RSF model. A total number of 1000 trees were considered. One sample tree from pool of the 1000 trees of the created RSF model for the ROP data is shown in Fig. 2. The RSF was used for variable selection. Figures 3 and 4 show the minimal depth and variable importance plots of the selected risk factors of developing ROP, using RSF. The number of days receiving non-invasive ventilation, GA (weeks), number of days receiving free oxygen, bilirubin levels, duration of receiving antibiotics (days), duration of TPN use (days), mother age (years), birth order, number of days receiving invasive ventilation, number of times receiving surfactant, and timely screening were the most important variables correlated with the time to development of ROP using both criteria of minimal depth and VIMP.

Fig. 2

One sample tree of the 1000 trees of the created RSF for the Retinopathy of prematurity (ROP) data

Fig. 3

Minimal depth plot of the selected risk factors of developing Retinopathy of prematurity (ROP); smaller values of minimal depth indicate a more important variable

Fig. 4

Variable importance (VIMP) plot of the selected risk factors of developing Retinopathy of prematurity (ROP); smaller values of minimal depth indicate a more important variable

Comparison between RSF and Cox regression model

RSF provided a lower prediction error compared to the traditional Cox proportional hazards model (results not shown), with the integrated Brier score of the RSF being 0.11 compared to that of the Cox model (0.25). Furthermore, the C-index of the RSF was 0.878 indicating great concordance between predicted and observed survival times. The C-index value for the Cox model was 0.621.

Associations between selected variables and Time to ROP

The association between the selected variables and the survival time (time to development ROP) was also investigated using the RSF analysis. The resulted partial survival plots were displayed in Fig. 5 at 12 weeks, indicating estimated survivals of levels of each risk factor considering the effects of all other risk factors to be justified. According to the Fig. 4, the nonlinear associations between selected risk factors and the hazard/survival rate of developing ROP are evident. For example, the survival rate of developing ROP decreases nonlinearly as the number of days receiving non-invasive ventilation increases. The slope of the predicted line decreases after 10 days. Additionally, as the mother’s age increases up to 28 years, the survival rate is increasing, but then after this point, it decreases, so that the hazard rate of developing ROP increases for infants with mothers over 28 years old. For categorical selected risk factors, the relationships were provided by each level. For example, infants with a timely screening had a higher survival rate compared to others. Also, infants with a higher birth order had a lower survival rate.

Fig. 5

Partial 3-month predicted survival for selected influential covariates on time to developing Retinopathy of prematurity (ROP). Dashed red lines indicates ± 2 standard error bars

Figure 6 shows the RSF estimated 3-months survival as a function of BW, duration of receiving non-invasive ventilation and bilirubin levels. This plot displays the interaction between the three important variables. The infants with BW less than 1000 g with bilirubin values between 15 to 20 mg/dl had the worst survival rate (the first column, second row from above). Survival probability was the best for infants with BW over 2000 g and higher values of bilirubin (see last column) and further dependent on changes in the duration of ventilation.

Fig. 6

Random survival forest estimated 3-months survival as a function of birth weight, gestational age, and bilirubin. Smoothed curves are loess curves of the estimated survival

In this group, the 3-month predicted survival was over 90% for those infants who did not receive non-invasive ventilation. The survival probability reduces to 70% for the infants with bilirubin values over 10 (mg/dl) and to about 20% for the infants with bilirubin values less than 10 (mg/dl) as the duration of receiving non-invasive ventilation increase. It is important to note that these interactions and non-linear relationships were identified by the RSF, and not pre-specified by the analyst.

Identifying high risk infants

We also conducted cluster analysis based on the selected variables using RSF to identify high and low-risk groups for developing ROP. Figure 7 shows the survival curves for the two groups of high-risk (90 infants) and low-risk (586 infants). The Log-rank test revealed that there was a significant difference between two identified groups (Test statistics: 10.7; P < 0.001).

Fig. 7

Survival curves for the two identified risk groups based on selected variables in the RSF analysis

For further investigation, we categorized the survival groups into low and high risk as a binary response to evaluate the sensitivity and specificity of individual variables in predicting the time to the development of ROP. Given that the RSF method requires at least three variables, including at least one mandatory continuous variable, we employed the multi-layer perceptron method, which can utilize a single predictor due to its model structure.

Table 2 presents the results. According to our findings, the duration of receiving TPN, the duration of receiving antibiotics in infants, the number of TB pack injections, the number of days on ventilation and oxygen, GA, and the administration of corticosteroids in infants exhibited the highest sensitivities in predicting ROP as a binary outcome. However, it is important to note that these variables were evaluated individually in the model. When considering all variables together, the sensitivity was found to be 0.910, while the specificity was 1.0, resulting in an AUC of 0.999. This remarkable performance underscores the intricate relationships and interactions among the predictors, highlighting the necessity of accounting for these complexities in the modeling process.

Table 2 Sensitivity and specificity of included variables in the model in predicting time-to-ROP based on the univariate analysis

Discussions

ROP is considered as the main preventable cause of blindness in infants in developing countries [20, 23, 24]. Despite several studies conducted to determine risk factors of developing ROP, little is known about the potential risk factors for the time to development of ROP, and few attempts have been made in this regard. In this study, we investigated the predictors of time to development of ROP in premature infants in Iran using machine learning models. The findings of this study can be used to improve the quality of services provided for the premature infants.

This study found that the duration of receiving ventilation (invasive/non-invasive), GA, duration of receiving oxygen supplementation, bilirubin levels, duration of receiving antibiotics, duration of TPN use, mother age, birth order, number of times receiving surfactant, and on time screening were the most important predictors of the time to development of ROP in premature infants in the present study. This aligns with the findings of other studies. In a study conducted by Shah et al., it was reported that BW (g), GA (weeks), Apgar score in first minute, duration of oxygen therapy, and mechanical ventilation for ROP are potential risk factors for the incidence of ROP [25]. Shah et al. reported risk factors including: BW, GA, Apgar in the first minute, intraventricular hemorrhage, duration of oxygen therapy and mechanical ventilation for ROP [25]. Similarly, Seiberth et al. reported that BW, GA, mechanical ventilation longer than 7 days, and surfactant administration are potential risk factors of the incidence of ROP [26]. Also, a number of studies reported on risk factors for the development of ROP. Recently, Owen et al. performed a retrospective cohort analysis of preterm infants referred for ROP screening [27]. They found that GA, BW, the need for surgery of any nature, and maternal magnesium prophylaxis were related to the development of ROP, irrespective of the stage of the disease.

In this study, it was found that the bilirubin level is a potential risk factor (nonlinear association) for developing ROP. According to the findings, although the survival probability increases as the bilirubin increases, it diminishes for the very high values of bilirubin. According to studies, the high serum bilirubin values can be harmful and lead to severe neurological damage; however, studies have shown the protective effect of bilirubin and biliverdin cycling against oxidative stress, which acts as an antioxidant [28]. There is no consensus on the protective effect of higher bilirubin levels for ROP risk in premature infants [29]; with both protective and harmful effects reported [30, 31]. “Furthermore, it is well known that extremely high levels of bilirubin can lead to severe neurological damage and needs to be treated. However, knowing the lower limit of bilirubin to maintain its antioxidant effect may be as important as knowing the upper limit to avoid its toxic side effects [14].” In our study, it was found that the survival probability tends to decrease for bilirubin values greater than 8 mg/dl. Bilirubin is a naturally occurring antioxidant that can protect against oxidative stress, which is thought to play a role in the development of ROP [29]. In premature infants, the retinal vascular development is stimulated by low tissue oxygenation in the fetus, but after extremely preterm birth, both room air and oxygen supplementation cause a state of relative hyperoxia, which can disrupt the delicate balance of retinal vascular development [29]. High bilirubin levels can counteract this hyperoxia and protect against ROP development. However, very high bilirubin levels can have adverse effects on the brain and other organs, leading to a diminished survival probability [32]. Therefore, it is crucial to maintain an optimal bilirubin level in premature infants to protect against ROP development while minimizing the risk of bilirubin-induced neurological damage.

Our findings indicate that maternal age is a risk factor for the time to the development of ROP. We observed a nonlinear relationship between maternal age and survival/hazard, with a high risk of developing ROP being associated with both very young and very old maternal ages. The suggest that extremes of maternal age, may pose increased risks for ROP in premature infants. Maternal age has been regarded as a significant factor influencing the outcome of childbirth [33]. Studies have reported higher risks of preterm birth and very preterm birth for both younger and older mothers [34]. In younger ages, this may be due to inadequate prenatal care, poor nutrition, and a higher likelihood of pregnancy complications. For older mothers, the increased risk may be associated with other underlying maternal health conditions or pregnancy complications that can indirectly influence ROP risk in preterm infants [35]. Older maternal age is also reported to be related to subfertility, stillbirth, chromosomal abnormalities, multiple gestation, and abortion [36,37,38]. Some studies also reported an inverse relationship between maternal age and the incidence and progression of ROP [39], while some others failed to find an association [25, 40]. This suggests that the relationship between maternal age and developing ROP may be complex and influenced by other factors. One possible explanation for the discrepancies in the findings could be differences in study design, population, and sample size. Also, the inconsistency between the results of various studies may be due to the fact that all these studies have used traditional statistical models and consider a linear relationship. Therefore, further research with more advance statistical/machine learning methods is needed to clarify the relationship between maternal age and ROP, taking into account potential confounding factors and study design limitations.

According to our findings, there was a relatively high (based on the lower depth or greater VIMP) association between oxygen supplementation (and duration of invasive / non-invasive ventilation) and the survival/hazard of developing ROP among premature infants, so that longer duration of receiving oxygen is associated with lower survival probability. This is in agreement with the results of other studies [41], which confirmed the role of excessive oxygen supplementation in developing ROP in low and middle income countries [42]. Insufficient oxygen levels in the tissue encourage the growth of blood vessels in the fetus's retina. However, following a very premature birth, both exposure to normal air and the addition of oxygen lead to a condition of increased oxygen levels [43]. In conclusion, the association between oxygen supplementation, ventilation duration, and the risk of ROP in premature infants underscores the importance of individualized care that balances the benefits of oxygen therapy with the potential risks of ROP development. Close monitoring, appropriate oxygen management, and tailored ventilation strategies are key in reducing the incidence and severity of ROP in this vulnerable population. Therefore, these should be prioritized in clinical assessments or in designing preventive strategies, especially for those infants with low GA.

The findings of the present study revealed that the time to developing ROP was adversely associated with the duration of receiving antibiotics in infants. A study has shown that receiving antibiotic treatment for a period of 4 to 7 days is associated with increased risk of mortality or significant morbidities, including ROP. In the present study, a large number of infants received antibiotics treatment (96.7%) and 66.5% of them received antibiotics treatment for at least 6 days. According to the studies, each additional day of receiving antibiotics increases the odds of mortality and severe morbidity by 14% [44]. Antibiotics can disrupt the gut microbiota in infants, potentially leading to systemic effects on their health, including susceptibility to conditions like ROP. This is biologically plausible because the gut microbiota plays a crucial role in the development and maturation of the immune system, and antibiotics can alter the composition and diversity of the microbiota. A study found that early empiric antibiotic use in preterm infants was associated with lower bacterial diversity and higher relative abundance of pathogenic bacteria, such as Enterobacteriaceae, which is a risk factor for Necrotizing Enterocolitis (NEC) and sepsis [45]. The study also found that infants who received 5–7 days of empiric antimicrobial agents in the first week had an increased relative abundance of pathogenic bacteria, and they experienced more cases of NEC, sepsis, or death than those not exposed to antibiotics [45]. Another study found that prolonged initial empirical antibiotic treatment was associated with adverse outcomes in premature infants, including increased risk of NEC, sepsis, and death [46]. The study suggested that this practice is a potential target for antimicrobial stewardship. A review of studies found that antibiotic exposure in the neonatal period appeared to induce various potentially disease-promoting alterations in the gut microbiota, including reduced diversity and increased abundance of pathogenic bacteria, which can lead to antibiotic resistance development [47]. The review emphasized the need to reduce unnecessary antibiotic treatment in neonates, including improving preventive measures, stopping antibiotic therapy after 36–48 h if infection is only vaguely suspected and there is no growth in the blood culture, and restricting the empirical use of broad-spectrum antibiotic treatment. In summary, the findings suggest that antibiotics can disrupt the gut microbiota in infants, potentially leading to systemic effects on their health, including susceptibility to conditions like ROP. Antibiotic stewardship is crucial to reduce unnecessary antibiotic treatment in neonates and prevent adverse outcomes.

TPN was also among the selected risk factors of time to development of ROP in this study. Our findings revealed a nonlinear association of the duration of receiving TPN, so that as the duration of receiving TPN increases up to 15 days, the survival probability increases, or the hazard of developing ROP decreases; however, the survival probability tends to decrease sharply. A study conducted on over 11,000 infants in Sweden, revealed that TPN was significantly associated with the severity of ROP (Spearman correlation coefficient of r = 0.45; P < 0.001) [48]. They showed that the infants who received parenteral nutrition for ≥ 14 days had a higher odd of developing ROP (1.84; P < 0.001). This is in consistency with our findings. Interestingly, the obtained cutoff for TPN was similar to that study and others [49]. The observed nonlinear relationship suggests that there may be an optimal duration of TPN that is associated with a lower risk of developing ROP. The findings indicate that a moderate duration of TPN, up to 15 days, may have a protective effect against ROP development, but prolonged use beyond this threshold could potentially increase the risk of ROP. Our findings highlight the importance of carefully managing the duration of TPN in preterm infants to optimize their nutritional support and potentially reduce the risk of developing ROP. Balancing the duration of TPN to avoid both deficiencies and excesses is crucial in promoting better outcomes for these vulnerable infants.

Our study found that the duration of receiving surfactant was inversely associated with survival probability, meaning that a longer duration of surfactant therapy was associated with a lower time of developing ROP in preterm infants. Surfactant therapy is a common treatment for Respiratory Distress Syndrome in premature infants, which can disrupt gas exchange and lead to high oxygen saturation levels, damaging newly developed retinal capillaries. This damage can stimulate the overexpression of angiogenic factors, leading to vasoproliferation observed in ROP [50, 51]. Surfactant therapy is essential for extremely preterm infants, usually improving their respiratory condition. While single dose does not reduce mortality, it helps with respiratory distress syndrome and pneumothorax. However, some infants may require more doses due to worsening respiratory issues after initial improvement [52,53,54]. The results of other studies confirmed this finding and they have shown that receiving more surfactant is associated with adverse neonatal outcomes [55] including mortality. Our finding also revealed an association between BW and time to development of ROP. This finding is in line with other studies, where low birth weight [56, 57] was found to be a significant risk factor.

According to our findings, the RSF showed a greater C-index and lower IBS compared to the Cox regression model (the classical model for survival analysis). This finding was in agreement with those of other studies; in various diseases including progressing from HIV to AIDS [58], kidney graft failure [59], breast cancer [60, 61], survival of patients with hemodialysis [62], gastric cancer [63], colorectal cancer [64], etc. it has been shown that the RSF outperforms the Cox model.

The main limitation of this study was its retrospective design which makes it impossible to verify the quality control of the used data. Another limitation was that information on neonatal hyperglycemia; low levels of insulin-like growth factor, nutrition, ω-3 long-chain polyunsaturated fatty acids, maternal breast milk and other risk factors reported by some other studies were not available. An additional limitation of our study is that the data were collected from referral hospitals that frequently admit complex cases referred from other healthcare facilities. This aspect may restrict the generalizability of our findings. Despite the limitation, the strength of the study was using a state-of-the-art machine learning model to analyze the data. The used model handles the non-proportionality hazards of the variables (a mandatory assumption in classical models) as well as dependency among observations and multicollinearity. It is worth noting that ROP is a multifactorial condition and all covariates may affect simultaneously and have complex interactions with each other, or they may have nonlinear effects on the time to developing ROP. The used model can extract interactions between variables (i.e. linear, nonlinear) without prior knowledge and provide predictions for survival probabilities considering these interactions, the model can consider the relationship between one covariate on time to developing ROP controlling for other variables. The current study revealed important prognostic factors for time to ROP development in premature infants.

Conclusions

We found that the number of days of receiving non-invasive ventilation, BW, and number of days receiving free oxygen, bilirubin levels, duration of receiving antibiotics, GA, duration of receiving TPN, and mother age were the most important predictors of time to developing ROP in infants. Based on our findings, it is recommended that healthcare professionals closely monitor and manage the number of days infants receive ventilation and oxygen, according to the BW, oxygen, bilirubin levels, and duration of antibiotic use, GA, duration of total parenteral nutrition, and maternal age to help identify and intervene early in infants at risk of developing ROP.