Study design and data source

This study utilized a merged database from the National Health Insurance Service (NHIS) and the National Health Screening Program for Infants and Children (NHSPIC) in Korea. The NHIS is a single insurance system covering nearly the entire Korean population, making it a representative data source. The NHIS database provides baseline demographic characteristics, such as birth date, sex, insurance premium, and region of residence, as well as information on healthcare utilization, including the type of hospital visit, diagnosis codes (International Classification of Diseases 10th revision [ICD-10] codes), prescribed medication codes, and procedure codes. All children in Korea were eligible to undergo seven rounds of the NHSPIC, which were conducted at specific age intervals from four to 72 months of age. The rounds were scheduled as follows: 1st (4–6 months old), 2nd (9–12 months old), 3rd (18–24 months old), 4th (30–36 months old), 5th (42–48 months old), 6th (54–60 months old), and 7th (66–72 months old). The NHSPIC survey includes a general health questionnaire, the Korean Developmental Screening Test (K-DST), an anthropometric examination, and a physical examination [13].

The de-identified individual data were used only for research purposes. Patient consents were not required as this study was based on de-identified and publicly available data. The Institutional Review Board of the Korea National Institute for Bioethics Policy waived the need for informed consent. The study protocol was reviewed and approved by the Institutional Review Board of the Korea National Institute for Bioethics Policy (P01-201603-21-005). All methods were performed in accordance with the relevant guidelines and regulations.

Study population

The study population is shown in Fig. 1. We included out of the 2,395,966 Korean children born between 2008 and 2012 those children who received all rounds of the NHSPIC from the first to the fourth round, responded appropriately to the questionnaire (n = 408,077) and received the K-DST properly in 7th round (n = 714,364). In total, 180,563 children met the inclusion criteria. Subsequently, children who met the following criteria were excluded: (1) died (n = 52), (2) birth weight <2.5 kg (n = 8029) or >4 kg (n = 6193), (3) multiple births (n = 1899), (4) preterm birth (n = 6427), (5) diagnosis with disorders of newborn related to the length of gestation and fetal growth (n = 7139), convulsions of newborn/disturbances of cerebral status of newborn (n = 456), or congenital malformations/chromosomal abnormalities (n = 29,130), (6) admission in intensive care unit over 4 days before 1 year of age (n = 5458), and (7) who received general anesthesia before 1 year of age (n = 2615) and for >5 days before 2 years of age were excluded. Ultimately, we enrolled 133,243 eligible children.

Fig. 1Dietary patterns in young childhood

The information on dietary patterns from infancy to 3 years of age was provided from the NHSPIC questionnaire, spanning the first through fourth rounds. The details of the questionnaire were described in Supplementary Table 1. Specifically, the first round of the NHSPIC, conducted at the age of 4−6 months, includes questions about the types of milk infants consume. The second round, conducted at the age of 9−12 months, includes questions about the initiation time of the introduction complementary foods, the frequency of intake of complementary foods, and the ingredients included in complementary foods. The third round, conducted at the age of 18−24 months, included questions about the consumption frequency of fruit juice or sweetened beverages. Finally, the fourth survey, conducted at the age of 30−36 months, includes questions about the frequency of fruit juice or sweetened beverage consumption, meal frequency, and milk intake.

Clusters according to dietary patterns during young childhood

Polytomous Variable Latent Class Analysis (poLCA) was used to identify groups of similar cases within the manifest variables for dietary patterns during young childhood and determine whether they were statistically independent. We generated a series of models featuring a diverse range of latent clusters, spanning from two to ten. We evaluated the performance of each model to determine the optimal fit of the data and the greatest possible distinction between the identified clusters. We utilize several statistical measures to evaluate the quality of the model fit, including the maximum log-likelihood plot, which indicates the point at which the maximum log-likelihood ceases to increase significantly, and the elbow heuristic for the Bayesian Information Criterion (BIC) and Akaike Information Criterion (AIC), where the change in successive values becomes less noticeable. (Supplementary Table 3 and Supplementary Fig. 1) [14,15,16,17,18]. In addition, entropy values greater than 0.6 indicate good cluster separation [19, 20], and we considered the distribution of clusters acceptable when each cluster comprised more than 3% of the total participants. Based on the final model, four clusters were determined to provide the best fit.

Developmental status at preschool age

The preschoolers’ developmental status was assessed using the K-DST performed at the age of 66–72 months, which is a valid screening tool designed specifically for Korean children and is part of the NHSPIC inventory [21, 22]. The K-DST consists of six domains: gross motor, fine motor, cognition, language, sociality, and self-care. Each domain consisted of eight questions answered by a parent or legal guardian, and the results were interpreted in four stages. These stages were: advanced development (total score ≤1 standard deviation [SD] score), age-appropriate (total score ≥–1 SD score and <1 SD score), need for follow-up (total score ≥–2 SD score and <–1 SD score), and recommendations for further evaluation (total score <–2 SD score). Children whose results indicated the need for follow-up underwent retesting or further evaluation if the interviews indicated problems. If the results for any of the six domains indicated the need for follow-up or recommendations for further evaluation, the total K-DST score was considered to reflect the same. The outcome of interest was an unfavorable outcome of K-DST, defined as a result of a “need for follow-up” or “recommendation for further evaluation” in each domain or total score.

Covariates

Demographic variables such as sex, region at birth, economic status, and birth year were obtained from the NHIS database. The regions at birth were classified as Seoul, metropolitan, urban, or rural. Health insurance premiums were determined based on economic factors, including income level and assets. Therefore, the economic status was categorized into quintiles using health insurance premiums as the criteria for evaluation. Moreover, birth weight and head circumference at 4–6 months of age were considered baseline clinical variables and were obtained from the first round of the NHSPIC. In addition, diagnosed perinatal conditions, as baseline clinical variables, were observed using P-codes in ICD-10 codes. These condition included fetuses and newborns affected by maternal conditions, birth trauma, respiratory and cardiovascular disorders specific to the perinatal period, infections specific to the perinatal period, hemorrhagic and hematological disorders of the fetus and newborn, transitory endocrine and metabolic disorders specific to fetuses and newborns, digestive system disorders of the fetus and newborn, and conditions involving the integument and temperature regulation of the fetus and newborn. Furthermore, atopic dermatitis or food allergies, which can influence dietary habits, was assessed (details of disease definitions are provided in Supplementary Table 2).

Statistical analysis

Categorical variables are expressed as the total number (n) and percentage (%), and continuous variables are described as mean and SD. Categorical variables between clusters were compared using the chi-square test, and continuous variables were compared using the Student’s t test. A multivariate logistic regression model was used to calculate odds ratios (ORs) with 95% confidence intervals (CIs) to identify the associations between dietary patterns and unfavorable K-DST outcomes. In addition, interaction p value between ORs were calculated by comparing the logarithmic differences of the ORs. The standard errors of these differences were used to derive Z-scores, from which p values were obtained to assess statistical significance. All analyses were adjusted for sex, region at birth, economic status, calendar year at birth, birth weight, head circumference at 4–6 months of age, perinatal conditions, and comorbidities. All the analyses were performed using the poLCA package (ver. 1.6.0.1) of R package (ver. 4.1.3) and SAS version 9.4 (SAS Institute Inc, Cary, NC, USA). Two-sided P < 0.05 was considered statistically significant.