Early childhood trajectories of domain-specific developmental delay and gestational age at birth: An analysis of the All Our Families cohort

Objective To describe developmental domain-specific trajectories from ages 1 through 5 years and to estimate the association of trajectory group membership with gestational age for children born between ≥34 and <41 weeks gestation. Methods Using data from the All Our Families cohort, trajectories of the domain-specific Ages & Stages Questionnaire scores were identified and described using group-based trajectory modeling for children born ≥34 and <41 weeks of gestation (n = 2664). The trajectory groups association with gestational age was estimated using multinomial logistic regression. Results Across the five domains, 4–5 trajectory groups were identified, and most children experienced changing levels of risk for delay over time. Decreasing gestational age increases the Relative risk of delays in fine motor (emerging high risk: 1.46, 95% CI: 1.19–1.80; resolving moderate risk: 1.11, 95% CI: 1.03–1.21) and gross motor (resolving high risk: 1.21, 95% CI: 1.04–1.42; and consistent high risk: 1.64, 95% CI: 1.20–2.24) and problem solving (consistent high risk: 1.58 (1.09–2.28) trajectory groups compared to the consistent low risk trajectory groups. Conclusion This study highlights the importance of longitudinal analysis in understanding developmental processes; most children experienced changing levels of risk of domain-specific delay over time instead of having a consistent low risk pattern. Gestational age had differential effects on the individual developmental domains after adjustment for social, demographic and health factors, indicating a potential role of these factors on trajectory group membership.


Introduction
Gestational age exists along a continuum, and the risk of developmental delays has an inverse dose-response association with gestational age [1,2].At a population level, the gestational age distribution at birth has shifted over the last 40 years, with the mean gestational age decreasing from 40 to 39 weeks [3,4].Though adverse developmental outcomes are more severe and best established with infants born extremely preterm, the increasing numbers of infants born between �34 and <39 weeks may result in a higher total number of children at risk in the population [5].A growing body of literature demonstrates this potential gradient of risk for developmental delays, neurodevelopmental disabilities, cognitive and academic difficulties, behavioral and emotional challenges across the gestational age continuum [6][7][8][9][10].
The most receptive stage of development is during early childhood; during this period, there is increased potential to permanently alter developmental trajectories and minimize disadvantages from accumulating through intervention [11].Many interventions have sought to improve outcomes for infants born before 34 weeks gestation [12], but few interventions have been developed that target or include infants >34 weeks [8].Understanding the short and long-term risk associated along the continuum of gestational age is vital to developing early childhood interventions for developmental delay [8].However, few studies have examined longitudinal development patterns or how those patterns of development differ across domains [1,2,[13][14][15].The objectives of the present study were to 1) describe developmental domain-specific trajectories as measured by the Ages & Stages Questionnaire 3 rd edition (ASQ-3) from 1 to 5 years of age, and 2) to estimate the association of trajectory group membership with gestational age for children born between �34 and <41 weeks gestation.

Study design
This observational study uses data from the All Our Families (AOF) study, a communitybased prospective pregnancy cohort that recruited 3387 pregnant women via prenatal serology testing, primary healthcare clinics, and community posters in Calgary, Alberta, from May 2008 until May 2011 [16].Details of AOF cohort methods, including study design, setting, and participant selection/eligibility/follow-up, have been reported elsewhere, and persons were eligible if they were �18 years of age, accessing prenatal care in Calgary, could complete questionnaires using English, and their pregnancy was �25 completed weeks gestation [17].Demographically, the AOF cohort is broadly representative of the pregnant and parenting population of urban centers in Canada [16].
Participants were initially invited to complete three questionnaires at �25 weeks pregnant, 34-36 weeks pregnant, and 4 months post-partum and consented to linkage to labor and delivery health records [16,18].Annual follow-up questionnaires were sent out on the child's birthdate to those participants who consented to be contacted for future research [16].Some follow-up questionnaires were delayed to secure funding and ethical approval, resulting in participants becoming ineligible as their children aged out of the questionnaire's administrative window (Fig 1) [16,19].Additional inclusion criteria for the current study (n = 2644) were participant dyads �18 years of maternal age who speak English with singleton pregnancies and children born �34 and <41 weeks of gestation with no congenital anomalies (Fig 1).

Participant consent and ethics statement
Cohort participants were informed of and agreed to (via signature of consent forms) the known research intent, risks, benefits, and privacy considerations during initial recruitment at each follow-up questionnaire for themselves and their child [16].Participants included in the current study also agreed, via the same informed consent processes, for their de-identified research data to be stored and used for secondary data analysis subsequent to ethical review.The present study, using de-identified participant data, was approved by the Conjoint Health Research and Ethics Board at the University of Calgary (REB21-0249).

Exposure
The exposure of interest, gestational age at birth, was derived from early ultrasound dating within electronic medical records (EMR) or maternal reported last menstrual period if EMR data were unavailable.EMR gestational age was collected for 85% of participants in this study.Previous analysis using this cohort data showed high validity for the maternal report of gestational age at birth and labor and delivery outcomes [18].This study included only those born �34 and <41 weeks gestational age, and gestational age was reverse coded with 40 weeks as a baseline for modeling.

Outcome
Child developmental outcomes were measured using the ASQ-3, a caregiver-completed screening tool [20].The ASQ-3 was administered at 1, 2, 3, and 5 years of age [16,20].The ASQ-3 contains 30 items divided into five domains of development: communication, problem solving, personal social, fine motor, and gross motor.Each item is assigned a point value to generate a domain-specific score out of 60, with higher scores indicating a lower risk of delay within the specific domain.Scores 1 standard deviation (SD) below the United States (US) normalized mean suggest that the child is at risk of delay and should be monitored, and scores �2 SD below indicate that the child is at risk of delay and should be further assessed [20].

Potential confounding variables
Various perinatal sample characteristics and factors influencing child development were assessed as confounding variables.Maternal age at delivery (years), maternal education (�high school or >high school), maternal self-identified ethnicity (dichotomised into White or BIPOC [Black, Indigenous, People of Color] due to sample size limitations), parity (primiparous or multiparous) are self-reported variables.Antenatal anxiety (�40 on the Spielberg State Anxiety scale) [21], antenatal depression (�10 on the Edinburgh Postpartum Depression Scale) [22], and antenatal stress (�19 on the Perceived Stress Scale, which is 1 SD above the mean [23]) variables were collected via validated maternal self-report screening tools during the second trimester.Binge drinking (defined as consuming �5 alcoholic drinks), smoking (cigarettes), and drug use (defined as broadly as street drugs, with no differentiation for cannabis) were self-reported variables assessed for the entirety of pregnancy, including before pregnancy recognition.Child sex (male/female), small for gestational age (SGA; birthweight below 10 th percentile of sex-specific birthweight as a proxy for intrauterine growth restriction [24]), and neonatal intensive care unit admission (NICU; a proxy for birthing related morbidities [25]) variables originated from the newborn discharge abstract of the EMR.

Statistical analysis
All analyses used Stata 16 IC software [26].Descriptive analysis presents the distribution of sample characteristics and factors influencing child development for the full sample, and across the range of gestational ages via stratification by gestational age categories.These gestational age categories are: late preterm (LPT: �34 and <37 weeks gestation), early term (ET: �37 and <39 weeks gestation), and full term (FT: >39 and <41 weeks gestation).
Trajectories on ASQ-3 domain-specific scores between 1 and 5 years of age were examined using group-based trajectory modeling (GBTM).GBTM applies finite mixture modeling using a maximum likelihood estimation, combining single-group models within a multiple-group model structure to identify clusters of individuals with similar trajectories [27,28].This method estimates the probability of trajectory group membership used in subsequent regression [28].We considered each child's age at ASQ-3 administration at 1, 2, 3, and 5 years as time 0, 1, 2, and 3. Domain-specific ASQ-3 scores were modeled separately within the TRAJ STATA plug-in using a censored normal model to account for clustering at either end of the scale [29].Model selection used both the Bayesian Information Criterion and the Akaike Information Criterion for base specification [30], followed by altering the polynomial orders of the trajectories using the Wald z-test and visualization of the trajectories to determine the most parsimonious number of groups [28].The number of trajectory groups was also visually assessed via graphing the mean outcome values and confidence intervals throughout the model selection process to ensure that distinct trajectory groups were not combined or separated based solely on Bayesian Information Criterion and the Akaike Information Criterion values.Model adequacy was determined using: 1) average posterior probability of assignment for all groups, 2) odds of correct classification, and 3) probability of group membership [28].Proportions within each domain-specific trajectory group that scored 1 SD and 2 SD below the normalized mean at least once from 1 to 5 years of age were used to describe the level of risk of developmental delay among trajectory groups [20].Proportions within each domainspecific trajectory group that scored 1 SD below the normalized mean at each assessment point were used to characterize the trajectory groups' transitions between risk levels."Emerging" described an increase in risk over time, "resolving" described a decrease in risk over time, "consistent" described relatively constant risk over time, and "transient" described changing risk, which may be both increasing or decreasing over time but was not consistent.
Once the most parsimonious GBTM model had been fit for each domain-specific outcome, a multinomial logistic regression within the TRAJ plug-in estimated the association between gestational age and posterior probabilities of trajectory group membership using starting values [28].Incorporation of multinomial logistic regression within the trajectory model makes the probability of group membership a function of the covariates, allowing accurate estimation of the association of gestational age with group membership probability [28].Risk ratios (RR) were generated by appending observations with no outcome data to generate group membership probabilities based on the inputted TRAJ risk factor settings, and 95% confidence intervals (CI) were generated through parametric bootstrapping [28].Informed by crosstabulation of potential confounding factors and trajectory groups, the inclusion of confounding variables was determined through backward elimination.

Missing data
The dataset has no missing data for the exposure (gestational age), and GBTM accommodates missing outcome data due to intermittent or missed assessments.However, missing data within potential confounding variables does have the potential to bias the estimate of the association between gestational age and the probability of group membership.To address missing data, the percentage of missing data for all potential confounding variables was computed during the descriptive analysis.Providing that the fitted models did not include variables with greater than 5% data missingness [31], a complete cases approach is used to handle missing data.

Descriptive analysis
Overall, 2,644 children were included in this analysis, of which 6.1% were LPT, 29.7% were ET, and 64.2% were FT.Children born LPT were significantly more likely to be admitted to the NICU (neonatal intensive care unit) than infants born either ET or FT, and their mothers were significantly more likely to have experienced depressive symptoms during pregnancy when compared to infants born FT.No other significant differences (via exclusivity of presented confidence intervals) were noted in the distribution of sample characteristics between LPT, ET, and FT born infants (Table 1).

Description of trajectory groups by model
A 4-group model emerged as the best fit for the communication, personal social, and fine motor domains, while a 5-group model was used for gross motor and problem solving.Further details on model adequacy are presented in Table 2.These identified trajectory groups are graphically represented by line graphs of the of mean ASQ-3 scores from ages 1 though 5 in Fig 2.
Description of the sample sizes, polynomial order, slope direction(s), and levels and patterns of risk for these trajectory groups are detailed in Table 3.Three levels of risk were defined when considering the proportion of children within each trajectory group who score 1 SD below the corresponding domain-specific US normalized mean at least once before 5 years of age: low (below 10%), moderate (10-49%) and high (50-100%).Patterns of risk, characterized using the proportions of scores 1 SD below the normalized mean at each assessment point, describe how the trajectory groups move between levels of risk, are shown within Table 4.
Table 3 shows that within the communication domain, most of the sample (64.6%; 95% CI:62.8-66.4%) is identified as belonging within trajectory group 4. The proportion of children who score 1 SD below the domain-specific US normalized mean is 0.8% (0.4-1.5%); therefore, within the communication domain, trajectory group 4 is described as having a low level of risk of developmental delay (Table 3).Table 4 shows that 0.2% (0.0-1.3), of the children within trajectory group 4 score 1 SD below the domain-specific US normalized mean at 1 year, and 0.0% (0.0-0.4), 0.9% (0.4-1.9), and 0.0% (0.0-0.4) at ages 2, 3, and 4 respectively.As these proportions are relatively consistent over time, the pattern of risk for trajectory group 4 within the communication domain is described as consistent (Table 4), and trajectory group 4 is described as having a low consistent pattern of risk of developmental delay (Table 3).Within all other domains, the majority of participants are classified as moderate risk with either a transient (personal social) or resolving (fine motor, gross motor, problem solving) pattern.

Multinomial regression results
Table 5 shows the relative risk for the association of gestational age at birth with the probability of group membership to trajectory groups per domain-specific model.Gestational age, infant sex, maternal age, maternal ethnicity, and prenatal anxiety and depression were included as covariates in all models (see Table 6 for sample characteristics by trajectory group).Gestational age was not consistently associated with membership in high-risk trajectories.While membership in high-risk trajectory groups was often inversely associated with gestational age in crude models, these associations were attenuated following adjustment for confounders.In fully adjusted models, gestational age was inversely associated with membership in the high-risk group for problem solving, gross motor, and fine motor domains.Decreasing gestational age was not associated with an elevated risk of delay in the communication or personal social domains after adjustment for confounding.

Discussion
Using a large community-based prospective pregnancy cohort, this study sought to describe ASQ-3 domain-specific trajectories of child development from ages 1 through 5 years and estimate the association of trajectory group membership with gestational age for children born between �34 and <41 weeks, controlling for covariates.The five child development domains measured within the ASQ-3 were independently modeled to determine the number and shape of potential trajectory groups within each domain, and the trajectories were described according to the pattern and level of risk for developmental delay.This study found that the probability of group membership within high-risk groups for fine motor, gross motor, and problem solving domains of development increases with decreasing gestational age.
Our findings are consistent with previous studies showing a dose-response association between gestational age at birth and developmental delays assessed via the ASQ [1,2,14,15].In 2015, Schonhaut et al. found an inverse dose-response association between gestational age and developmental delay using the ASQ-3, looking cross-sectionally at infants born LPT, ET, and FT [2].In 2017, Hornman et al. found that children born preterm had increased odds of scoring below 2 SD on any domain compared to FT when using the Dutch version of the Ages & Stages Questionnaire [32].The literature supports the idea that gestational age for those born between �34 and <41 weeks is inversely related to an increased risk of developmental delay.In examining children longitudinally from 4 months through 3 years of age, Hochstedler et al. so found gestational age was inversely associated with the odds of scoring below 2 SD on any ASQ domain [14].However, literature examining the patterns of development from birth to 5  years of age is more heterogeneous and limited.Hornman and colleagues' (2017) paper found that the pattern of developmental delay between ages 4 and 5 years was different between EPT and LPT groups and between the differing domains [32].And, a 2018 paper, assessing ASQ-3 at 8-, 18-, and 24-months scores from children born between �34 and <41 weeks, found that decreasing gestational age at birth was a risk factor for boys who had scored 2 SD below the norm in 1 or 2 of the ASQ domains during that time [33].The current study shows an inverse dose-response association between gestational age and developmental delay over the entirety of early childhood, and that the association differs by developmental domain.
To optimize child developmental outcomes through the provision of early intervention, those who fall within the consistent or emerging high risk trajectories in a specific domain may benefit from further assessment.Those within the consistent low risk trajectories describe what we might consider typical development, reaching those domain-specific developmental milestones at the appropriate chronological age.However, in all domains (except communication) most children within this study experienced changing levels of risk of developmental delay, which can be described as either resolving or transient.In a cross-sectional analysis, as with a one-time developmental screening program, these children might even appear typical.However, many of these children encounter domain-specific difficulties at one point before 5 years of age.The purpose of screening tools such as the ASQ-3 is to identify children with disabilities or difficulties and those at risk for further challenges.These at-risk children may not qualify for formal intervention programs, but that does not mean they do not require or would not benefit from intervention.However, further research on factors concerning these resolving and transient groups, aside from gestational age at birth, is necessary to determine who may benefit most and from which interventions.

Strengths & limitations
Strengths of this study include the cohort study design, analyzing the continuous ASQ-3 scores, and using group-based trajectory modeling (GBTM) to estimate the association between gestational age and domain-specific developmental delay.The AOF cohort is a longitudinal prospective pregnancy cohort that followed all 2664 mother-child dyads included in this study from pregnancy until their child was 5 years old, consistently measuring the risk of developmental delay with the ASQ-3 at four separate time points [16].Utilizing prospective data collection minimized temporal bias of confounding factors, gestational age at birth, and developmental delay outcomes [34].
The ASQ-3 scores within each domain range from 0 to 60, and most of the research using this tool to measure developmental delays in children use a cut-off of 2 standard deviations from the normalized mean to indicate increased risk for developmental delay [1,2,14,15].The purpose of a screening tool, such as the ASQ-3 is to capture all persons possibly at risk [35].A 2018 study found that using 1 SD rather than 2 SD enhanced the test characteristics as a screening tool for neurodevelopmental disabilities as well as developmental delay [36].
Trajectories of risk of developmental delay were modeled in the current paper using the 0 to 60 domain-specific scores and were described in relation to the normalized mean cut-offs using the 1 SD cut-off to convey the level of risk of developmental delay.In using the continuous rather than the dichotomized scores, this study determined that increased risk of delay within a domain is rarely consistent, and using continuous data provides a more nuanced understanding of the risk of developmental delay.
GBTM compliments the cohort study design in this analysis.The longitudinal data enabled the application of GBTM to analyze the outcome of interest over four time points; thus, enabling cubic polynomial orders during modeling for greater precision in the group membership probabilities of trajectory groups [28].Further, GBTM reduces the effect of attrition or non-participation bias inherent to longitudinal cohorts as it accommodates missing data due to intermittent or missed assessments [37].Using GBTM, rather than fixed/random/mixedeffects longitudinal regression models, marginal models, or traditional growth curve models for analysis, we identified resolving, emerging, and transient patterns of increased risk which would otherwise be missed by methods requiring a priori grouping [38].
Certain limitations have been identified within this study's design, including analytical power, use of the ASQ-3 to measure developmental delay, and generalization of the study's findings.As we cannot a priori determine sample sizes for trajectory groups using GBTM, there are relatively small sample sizes within some trajectory groups.While these smaller trajectory group sizes within the high emerging or high transient groups within the present study are aligned with those within prior research on stability of risk of developmental delay [32], this particular analytical method has not been employed in a larger cohort to our knowledge so we cannot discount the possibility that within a larger cohort alternative trajectory groups may be observed.Additionally, trajectory group sample may have contributed to a lack of sufficient power within regression models to determine a significant association between gestational age and trajectory group membership when adjusting for confounding factors.The sensitivity of the ASQ-3 can change with age of assessment [39][40][41]; with previous research identifying increased sensitivity of the screening tool at age three, and as the sensitivity increases, there will be more children identified as "at risk of delay" [42].Though this potentially could explain the increased risk at age 3 seen in the "transient high risk" trajectory group within the fine motor domain in this study, the limitations of psychometric consistency throughout each domain over time should be noted.And finally, though at recruitment participants within the AOF cohort represented the demographic characteristics of urban Alberta, Canada, the elevated median socioeconomic status of this cohort's population may limit the generalizability of this study's findings [16].As prior research has shown that lower socioeconomic status increases the risk of early childhood developmental delay [43], the cohort's elevated socioeconomic status may underestimate the association between gestational age and risk of developmental delay.

Table 4 . Frequency proportions of participants which score below 1 standard deviation within each ASQ-3 (ages & stages questionnaire 3 rd edition) domain specific trajectory group over time.
These proportions were used during description of the trajectory groups to gauge the transitions between risk levels.Emerging describes an increase in risk over time, resolving describes a decrease in risk over time, consistent describes relatively constant risk over time, and transient describes changing risk which may be both increasing or decreasing over time but is not consistent.

Table 5 . Risk ratios for the association of gestational age at birth (per week) with probability of group membership per individual ASQ-3 (ages & stages question- naire 3 rd edition) domain.
used to guide initial covariate selection when modeling the association between gestational age at birth and the probability of group membership within domain specific trajectory groups.Bolded values are considered significantly different from the non-bolded values marked with the corresponding superscript as determined by the exclusivity of the 95% confidence intervals.Bolded values are considered significantly different from the non-bolded values marked with the corresponding symbol as determined by the exclusivity of the 95% confidence intervals.https://doi.org/10.1371/journal.pone.0294522.t006 RR: Risk ratio, CI: Confidence interval, REF: Referent group.b The risk ratio can be interpreted as the rate of change in the probability of group membership per week before 40 a Risk ratios are adjusted for gestational age, maternal age, ethnicity, anxiety, depression, and infant sex.weeks relative to the referent group.c statistically significant at α = 0.05.SGA: Small for gestational age, NICU: Neonatal intensive care, BIPOC: Black, Indigenous, People of Colour, CI: Confidence interval.a This table was b c