Retinopathy of prematurity: A comprehensive risk analysis for prevention and prediction of disease

Background Retinopathy of prematurity (ROP) is a blinding morbidity of preterm infants. Our current screening criteria have remained unchanged since their inception and lack the ability to identify those at greatest risk. Objectives We sought to comprehensively analyze numerous proposed maternal, infant, and environmental ROP risk variables in a robustly phenotyped population using logistic regression to determine the most predictive model for ROP development and severity. We further sought to determine the statistical interaction between significant ROP risk variables, which has not previously been done in the field of ROP. We hypothesize that our comprehensive analysis will allow for better identification of risk variables that independently correlate with ROP disease. Going forward, this may allow for improved infant risk stratification along a time continuum from prenatal to postnatal development, making prevention more feasible. Methods We performed a retrospective cohort analysis of preterm infants referred for ROP screening in one neonatal intensive care unit from 2010–2015. The primary outcome measure was presence of ROP. Secondary outcome measures were ROP requiring treatment and severe ROP not clearly meeting current treatment criteria. Univariate, stepwise regression and statistical interaction analyses of 57 proposed ROP risk variables was performed to identify variables which were significantly associated with each outcome measure. Results We identified 457 infants meeting our inclusion criteria. Within this cohort, numerous factors showed a significant individual association with our ROP outcome measures; however, stepwise regression analysis found the most predictive model for overall ROP risk included estimated gestational age, birth weight, the need for any surgery, and maternal magnesium prophylaxis. The corresponding Area Under the Curve (AUC) for this model was 0.8641, while the traditional model of gestational age and birth weight predicted ROP disease less well with an AUC of 0.8489. Development of severe ROP was best predicted by estimated gestational age (week), the need for any surgery and increased probability of death or moderate-severe BPD at 7 days. Finally, the model most predictive for type 1 ROP included estimated gestational age (week) and the presence of severe chronic lung disease. No significant statistical interaction was found between variables. Conclusions Our work is unique as we report comprehensive analysis of the greatest number of proposed ROP risk variables to date in a robustly phenotyped population. We describe novel risk models for our ROP outcome measures and demonstrate independence of these variables using statistical modeling not previously applied to ROP. This may better allow for individual infant risk stratification and importantly mitigation of future risk.


Objectives
We sought to comprehensively analyze numerous proposed maternal, infant, and environmental ROP risk variables in a robustly phenotyped population using logistic regression to determine the most predictive model for ROP development and severity. We further sought to determine the statistical interaction between significant ROP risk variables, which has not previously been done in the field of ROP. We hypothesize that our comprehensive analysis will allow for better identification of risk variables that independently correlate with ROP disease. Going forward, this may allow for improved infant risk stratification along a time continuum from prenatal to postnatal development, making prevention more feasible.

Methods
We performed a retrospective cohort analysis of preterm infants referred for ROP screening in one neonatal intensive care unit from 2010-2015. The primary outcome measure was presence of ROP. Secondary outcome measures were ROP requiring treatment and severe ROP not clearly meeting current treatment criteria. Univariate, stepwise regression and statistical interaction analyses of 57 proposed ROP risk variables was performed to identify variables which were significantly associated with each outcome measure.

Results
We identified 457 infants meeting our inclusion criteria. Within this cohort, numerous factors showed a significant individual association with our ROP outcome measures; however, stepwise regression analysis found the most predictive model for overall ROP risk included PLOS

Introduction
The reason for problematic incorporation of epidemiologic data into our screening paradigm is likely due to the variability in analysis approach, relatively small simultaneous factor analysis, changes in neonatology practice over time with regard to proposed risk factors such as oxygen supplementation, and small cohort size. For example, the majority of studies perform only univariate analysis, which does not account for the confounding effects of individual variable significance, chiefly known significant factors such as early gestational age. For those studies that do use a logistic regression or multivariate approach, most analyze only modest numbers of risk variables. [26], [31] Thus the analysis is not comprehensive with regard to suggested ROP risk variables. Finally, given the particularly low incidence of severe ROP as noted previously, sample size is often small for epidemiologic ROP studies. However, data compiled over a long duration, while it may increase sample size, includes infants treated very differently per neonatal intensive care unit (NICU) technology, target oxygen saturations, available interventions and medications. Thus, this also can introduce confounding.
Therefore, a knowledge gap remains, as despite the number of epidemiologic studies attempting to define ROP risk variables, we still cannot adequately predict infants at greatest risk for disease. [6], [32] More specific screening criteria are desperately needed however as our current screening exams are stressful for infants and caregivers and further, the majority of infants screened do not develop vision threatening disease. [20], [33] Thus, going forward we need to preserve our ability to identify all infants at risk for blinding disease, while better targeting our exams to appropriately exclude those at low risk. In our study we sought to address this knowledge gap through simultaneous analysis of the greatest number of proposed ROP risk variables to date in a roubstly phenotyped cohort to determine the most predictive models for ROP development and severity. We further sought to determine if there was a statistically significant interaction between significant variables, an analysis technique that has not previously been done in the field of ROP. Finally, we analyzed the predictive value of our models as compared with current screening criteria and a recently cited model. We hypothesize that by including a greater number of proposed ROP risk variables in a simulataneous anlaysis and within a well-phenotyped population receiving consistent NICU care, we would be better able to define individually singificant ROP risk factors and determine the most predictive risk models for ROP development and severity.

Materials and methods
Study Cohort: We performed a retrospective cohort study of preterm infants at one neonatal intensive care unit between 2010 and 2015. Infants met inclusion criteria if they were deemed appropriate for ROP screening per the clinical NICU guidelines. For the NICU in our study this included infants born less than 30 weeks gestational age and 1250 grams, or infants who displayed an unstable medical course as determined by the attending neonatologist. The ladder is common practice thorughout the US to ensure all infants with any ROP risk receive appropriate screening; neonatologists making this clinical judgement were not involved in study analysis therefore study enrollment was not a confounding factor in their decision making. Infants were excluded from our study if they did not meet our NICU clinical criteria for ROP screening. We did not set an a priori study size but rather designed our study to include years during which neonatology care practice and importantly oxygen saturation standards were uniform to remove these potentially confounding factors. All data collection was approved by the Institiutional Human Subjects Committees at the University of Utah, adhered to the Declariation of Helsinki, and was compliant with the Health Insurance Portability and accountability Act (HIPAA). Data were collected with waiver of consent for de-identified data. No authors have a conflict of interest. Data Collection: Infant and environmental risk factor data were collected for the duration of infant NICU stay. This duration is sufficient to document all collected indices. Maternal data were collected retrospectively at the time of infant enrollment. Data for a total of 57 indices were collected and analyzed. (Table 1) These risk factors were chosen through collaboration with pediatric ophthalmology and neonatology physicians who specialize in the diagnosis and treatment of ROP. In brief, factors were chosen based on the preponderance of literature suggesting their role in ROP and specialist recommendation. [23], [24], [26], [28] The probability of death or moderate-severe bronchopulmonary dysplasia at 7 and 14 days is a clinically calculated indices commonly assessed by neonatologists, as referenced. [34] Outcome Measure Phenotype: All infant ROP phenotypes were assessed by pediatric ophthalmologists during the course of normal clinical care. The primary outcome measure was development of ROP. Secondary outcome measures were 1) development of ROP disease requiring treatment termed type 1 ROP as defined by the Early Treatment of ROP study. [7] and 2) development of severe ROP which we defined as zone 1 or 2 stage 3 ROP or greater in the worse eye. We used the universally accepted Revised International Classification of ROP guidelines to define the location and extent of disease within the retina.
[35], [36] These criteria represent consensus guidelines for describing where the ROP disease occurs within the retina, whether it is close to the area of central vision and therefore more concerning (zone 1) or in a more peripheral portion of the retina (zone 3). These criterion also allow for standardized discussion of disease severity by stage. For example, stage 1 is the least severe form of ROP, correlating with a demarcation line between vascular and avascular retina and stage 5 is the most severe, correlating with total retinal detachment. Notably for our study, stage three as included in our "severe ROP" outcome measure, is the first stage where aberrant neovascularization is found. Our delineation of severe ROP therefore includes infants with eye disease that is worrisome by clinical standards, but does not have an absolute indication for treatment by current standards. For infants requiring multiple exams, ROP zone and stage for analysis was from the highest severity examination in the worse eye as is standard reporting practice for ROP studies. Importantly, these data were collected seperately from the time of clinical exam, therefore screening clinicians were not influenced by patient study inclusion.
Statistical Analysis: Proposed ROP risk variables were tested for association with each ROP outcome measure using logistic regression in SAS (v9.3, SAS Institute Inc., Cary, NC). Nominal significance was considered a p-value of <0.05 and estimates were reported as odds ratios with 95% confidence intervals. Subsequently, stepwise logistic regression was performed in order to determine the most predictive model for ROP. Specifically, those variables showing statistical significance at p <.05 in the single factor analysis were included in the stepwise regression model. Subsequently, each variable is added one by one to the model at a significance level of p .1. After the addition of each variable, variables already in the model are removed if they do not remain significant at p < .05.
To test the effectiveness of our model, we created receiver operating characteristic (ROC) curves using SAS to compare the stepwise regression model from the current study to both the traditional ROP risk model which informs our current screening criteria, namely gestational age less than 30 weeks and birth weight less than 1250 grams, and the model recently proposed by Slidsborg et al. to be a standard model. [26] Additionally, we performed interaction analysis to test for possible interactions within the current dataset. In order to control for potential confounders that may interfere with the interpretation of the interaction model, tests for interactions between the risk factors were performed using interaction terms as well as main effects in the logistic regression model, following the methodology proposed by Keller. [37] Specifically, factors that were shown to be

Percent/ Range
No./ Ave. nominally significantly associated with each ROP subtype (p < .05) were included in the model along with their corresponding interaction terms.

Results
546 infants were initally identified for inclusion in our study; 66 infants died prior to the first eye examination and 23 did not meet screening criteria due to birth weight greater than 1250 grams; 457 infants met inclusion criteria. Data were collected for the duration of the infant's NICU stay which was on average 10 weeks from the time of first eye exam. Baseline characteristics relative to each proposed ROP risk variable for our cohort of 457 infants, as represented in Table 1, show a fairly equal male and female population with a cohort average gestational age of 27.12 weeks and weight of 987 grams. Maternal race analysis demonstrated a predominantly Caucasian population, consistent with the demographic of Utah. We found that the   overall ROP incidence proportion in our population was 47.5%, which decreased to 12% for severe ROP, and 7.2% for type 1 ROP. This is consistent with the published incidence of disease among preterm infants in the US. [1] Univariate analysis was performed for all 57 indices, which have been suggested in the literature to contribute to ROP risk, though have not been comprehensively analyzed. A nominally statistically significant (p< .05) relationship was found between a number of factors and our ROP outcome measures. (Table 2) Importantly, we found factors historically known to influence ROP risk such as gestational age and birth weight to demonstrate a significant association for all outcome measures in the univariate analysis. [6] In order to determine the most predictive model of ROP, we performed stepwise regression analysis. When comparing infants in our cohort who developed any form of ROP to those who did not, we found the most predictive model of overall ROP risk included estimated gestational age, birth weight, the need for any surgery, and maternal magnesium prophylaxis ( Table 3). The corresponding concordance index, c, (an estimate of the area under the receiver operating characteristic (ROC) curve) for this model was 0.870. Development of severe ROP was best predicted by estimated gestational age (week), the need for any surgery and increased probability of death or moderate-severe BPD at 7 days (c = .978). (Table 3) Finally, the model most predictive for type 1 ROP included estimated gestational age (week) and severe chronic lung disease when comparing to those infants without ROP(c = .990). (Table 3) ROC curves were created for the stepwise regression model from the current study as compared to the traditional ROP risk model and standard model used by Slidsborg et al. (Fig 1). The "Traditional Model" used included the variables estimated gestational age and birth weight less than 1250 grams as these are the classically accepted ROP risk variables currently utilized to inform infant ROP risk and screening. The "Current Study" included the results from our stepwise regression analysis of any ROP: estimated gestational age, birth weight, the need for any surgery, and maternal magnesium prophylaxis. The "Slidsborg et al" model included the traditional model used in the Slidsborg paper: Small for gestational age, estimated gestational age, gender, and multiple gestation. For the purposes of this study, the variable intrauterine growth restriction was used to represent small for gestational age as this is the corollary clinical parameter used in our NICU. When comparing the area under the curve (AUC) for each of these models, the stepwise regression model from the current study gave the greatest AUC, 0.8641, while the traditional model and the traditional model used by Slidsborg et al. provided the same AUC, 0.8489. Interestingly, the model that produced the greatest AUC was that using all variables included in each of the three models: estimated gestational age, birth weight, birth weight less than 1250 grams, intrauterine growth restriction, gender, multiple gestation, the need for any surgery, and maternal magnesium prophylaxis (AUC = 0.8658).
Finally, we sought to determine the statistical interaction between ROP risk variables for each disease phenotype. Keller et al. have demonstrated this technique as a way to determine interaction between disease risk variables. [37] However, no significant statistical interaction was seen between variables tested in our model.

Discussion
ROP is an important clinical problem that if undetected can result in permanent, life-long blindness. With improved preterm infant survival, not only in the developed but also the developing world, the scope of the problem is increasing. [1] Our present screening guidelines are largely unchanged since their inception in the 1990's and [32] most sources estimate that only 5-10% of infants screened under these guidelines will develop vision-threatening disease. [20], [33] Therefore, while the sensitivity is uniformly high for current screening guidelines, the specificity is low. [20] Work to model risk and better predict ROP has been complicated by the multifactorial nature of this disease. [32] Further, current statistical modeling of ROP risk has shown variable sensitivity across ethnic populations and does not account for broad sources of ROP risk. [38], [39]

Birth Position
Breech vs. Normal <0.01 <0.01->999.99 0.980 n/a n/a n/a n/a n/a n/a Transverse vs. Normal <0.01 <0.01->999.99 0.976 n/a n/a n/a n/a n/a n/a Unknown vs. Normal <0.01 <0.01->999.99 0.982 n/a n/a n/a n/a n/a n/a Vertex vs. Normal <0.01 <0.01->999.99 0.980 n/a n/a n/a n/a n/a n/a Other vs. Normal 1.00 <0.01->999.99 0.989 n/a n/a n/a n/a n/a n/a doi:10.1371/journal.pone.0171467.t002 Herein, we provide the most comprehensive analysis to date of proposed maternal, infant, and environmental ROP risk variables. Using this approach, we find that most of the factors assessed, while significant in a univariate analysis, fail to show significance in stepwise regression analysis. Consistent with the majority of literature, early gestational age demonstrates significance for all of our ROP outcome measures. [6], [22] Individual significance of early birth as compared with low birth weight is supported by other studies. For example, Woo et al., have shown that in twin gestations discordant for gestational weight, gestational age was a better predictor of ROP disease. [40] However, a number of factors postulated to be significant in current literature by univariate and multivariate analysis such as pre-eclampsia, patent ductus arteriosus, maternal age, blood transfusion, male gender, and multiple birth, do not show individual disease significance when corrected for confounding significance in our population. [26] The variability of our findings compared with reported work speaks to both the strengths and weaknesses of our study. Strengths of our study include our comprehensive analysis of the most numerous set of proposed ROP risk variables published to date. Therefore, when accounting for a larger number of variables we may be better able to more accurately control for confounding effects of individual variables and better determine those variables that are independently predictive of ROP disease. Further, we report a robustly phenotyped cohort with uniform clinical care over a 5-year time period to further strengthen our ability to identify individually significant ROP risk variables. As is noted by Slidsborg et al., variability in clinical NICU practice with regard to oxygen administration etc as well as clinical data collection confounds analysis of variable contribution to ROP risk and limits inclusion or analysis of continuous variables. [26] Several factors however, may limit the application of our findings and require further study. In an effort to include only infants with uniform neonatal care and thus minimize this important confounder, we limit our sample size. This limits the number of observations for some of our risk variables, most notably for the less common and most severe, Type 1 ROP outcome measure and may introduce a type 2 statistical error. However, increased sample size does not insure increased observations, particularly with respect to our more severe ROP outcomes. In fact, assessment of a more homogeneous population will counter this limitation and reduce the likelihood of spurious association. [39] The dichotomy between sample size and homogeneity of sample population is common in ROP literature, as often maximizing sample size requires inclusion of infants with dissimilar perinatal care. As these interventions comprise a significant proportion of our ROP risk variables, and indeed those currently suggested in the literature, we felt it most important to include only infants with uniform neonatal care. We are additionally limited by the demographic of our NICU setting, which is predominantly Caucasian. Racial differences have been demonstrated in the sensitivity and specificity of other proposed ROP risk models and possibly in ROP incidence and therefore, our findings should be replicated in a more diverse population. [41] Finally, the predictive utility of our model versus other models is an imperfect comparison as the sample populations are different in our cohort versus the comparison Danish population. Therefore, there may be elements of our population that better predispose to prediction using our model versus the proposed Danish model; this is not a parameter that can be fully controlled and should be noted. Overall however, we argue that our presentation of a smaller cohort with more uniform clinical care and robust phenotyping offers valuable insight into ROP risk. Certainly replication of our findings in a larger, more racially diverse, population will allow for clarity on this point.
In addition to our stepwise regression analyses, we also conducted a robust interaction analysis. This alternate statistical modeling system takes into account the effects of confounding variables in a slightly different manner in order to determine the degree of interaction between variables. This has not been done previously to characterize ROP risk. Using this methodology, we found no significant statistical interactions between the variables considered. Thus, the risk factors we identified are independent of one another, and by this analysis, confer independent risk.
In conclusion, better understanding of factors important for ROP development is critical for improved infant risk stratification and blindness prevention. Our approach is unique in the breadth and number of covariates queried and the degree of statistical modeling employed. [20], [28] While our findings reiterate the accepted importance of early birth and low birth weight to ROP risk, we also find that additional consideration of need for any surgery and maternal magnesium prophylaxis create the most predictive model for development of ROP. When considering the area under the curve in our ROC analysis our proposed model better predicts ROP development than the parameters used currently to stratify infant ROP risk.
Our work has also defined novel risk variables for severe ROP disease. Specifically, we describe risk of severe ROP disease with increased probability of death or moderate-severe BPD at 7 days. This is an emerging parameter used by neonatologists to assess primarily lung outcomes in preterm infants, though our work suggests utility for severe ROP prediction as well. [34] Future work will allow us to understand the predictive value of the models we describe, and also importantly, whether modulation of risk variables, when possible, can be exploited to alter the ROP development or severity. For example, infants identified as at risk for ROP development using our model may benefit from delay of elective surgical procedures until they reach a gestational age not associated with ROP development. [32] Certainly future work is necessary to clarify the validity and potential role of the models generated herein to ROP risk stratification.

Author Contributions
Conceptualization: LAO ROH BAY.