This study was conducted to predict off-track development in infants between 0 and 6 months of age and to identify key predictors influencing this outcome. These insights are intended to inform early interventions aimed at improving developmental outcomes in low-resource settings. The following section outlines the specific steps and techniques used throughout the study.
Research environment and participants
The study was conducted in Kaloleni and Rabai, two rural counties located in Kilifi County on the coast of Kenya. These areas are among the most economically disadvantaged areas in the country.10 The region has an estimated population of 350,000 people and is home to approximately 47,000 households. Access to health care is limited, with only 40 health facilities serving a widely dispersed population. The region faces severe poverty, with approximately 70% of the population living below the national poverty line, and the majority, approximately 81%, relying on small-scale agriculture, handicrafts, informal labor and small-scale trade for their income. Furthermore, maternal and child health conditions in the region are significantly worse than the national average.11
The study included 1,995 singleton infants aged 0 to 6 months who were healthy at the time of data collection and had no reported symptoms such as fever or diarrhea. Participants were identified by community health extension workers (CHEWs) using the Kaloleni Rabai Health and Demographic Surveillance System (KRHDSS). The system is managed by the Department of Population and Health (DPH) at the Aga Khan University in Nairobi, Kenya.12 It is a population-based registry that monitors approximately 21,000 households in the region and has conducted data collection annually since 2017.12
Data collection and research tools
A cross-sectional study was conducted between March 2023 and March 2024 among 2085 children. This included 1995 singleton infants and 35 twin pairs with biological mothers and 20 singleton infants in the care of a nonbiological guardian, usually due to the absence or death of the biological mother. Eligible participants were identified at the household level with the support of CHEW and community health promoters (CHPs).
A data collection team comprised of trained researchers from the Aga Khan University Institute for Human Development, in collaboration with CHP, visited participants’ homes to collect data. Written informed consent was obtained for the participation of parents or mothers themselves and their children. Data were collected using tablets equipped with Open Data Kit software and through interviewer-administered questionnaires as described below.
This study utilized several validated tools to collect data. These include the World Health Organization Index of Infant and Young Child Development (WHO-IYCD), a caregiver-reported questionnaire that assesses motor, language, cognitive, social-emotional, and behavioral functioning.13 Patient health questionnaire regarding maternal depressive symptoms14;Generalized Anxiety Disorder Scale for Maternal Anxiety Symptoms.15 A home screening questionnaire was used to assess psychosocial stimulation in the home environment.16 Socioeconomic status was measured using a tool adapted from:17 and maternal and child health history was obtained through a yes/no parent report questionnaire. Finally, anthropometric measurements (height and weight) of both mothers and children were collected according to the procedures recommended by the WHO. Detailed information about these tools can be found in Supplementary Text S1.
data analysis
result variable
The outcome variable for this study was the child’s overall developmental status, classified as either “developmentally off track” or “typically developing.” Developmental status was assessed using the Developmental Z-Score for Age (DAZ), which was calculated from the WHO-IYCD score and the child’s age. DAZ values were generated using the Excel tool provided in the WHO IYCD Instrument Manual (version 1.2, Web Appendix VII).13 It automatically calculates a Z-score based on the child’s age and raw IYCD score. According to the WHO 2005 reference standard, a score < -2 is classified as “developmentally off track” and a score ≥ -2 is considered “typically developing”. The internal consistency of the WHO-IYCD scale was acceptable, with Cronbach's alpha of 0.91 (95% CI: 0.905 to 0.916).
Data preparation and exploratory data analysis
Data were analyzed using Python (version 3.11.3).18 The dataset contained 64 predictive features selected based on previous research and input from domain experts (Supplementary Table S1). Considering minimal missing data (Supplementary Table S1), imputation was performed using logical inference and statistical methods including median and mode to ensure consistency and precision. For numerical variables, missing values were imputed using the median to maintain consistency and reduce the effect of potential outliers. For categorical variables, missing entries were imputed using the most frequent category. Descriptive statistics were used to summarize maternal and child characteristics, including proportions for categorical variables and median interquartile range for continuous variables. Group comparisons were performed using the Pearson method. χ² test or Fisher’s exact test for categorical variables and the Kruskal-Wallis test for non-normally distributed continuous data. The TableOne package for Python was used. ‘htest_name = True’ Option to display the name of the applied statistical test. Fisher’s exact test is automatically applied to categorical variables when appropriate, and Pearson’s exact test is χ² If not, test it. All hypothesis tests were two-sided, and statistical significance was defined as: p < 0.05. To inform feature selection based on expert judgment and impact on model performance, Pearson correlation coefficients were computed and visualized pairwise for all continuous features considered important or tentative. For categorical variables, pairwise Cramér's V correlations were assessed and selections were also made based on domain expertise. To facilitate model interpretation, one-hot encoding was applied to categorical variables (Supplementary Table S2).
Model development and feature selection
This study evaluated three supervised ML models (classifiers): Ridge Logistic Regression (LR), Advanced Random Forest (RF), and Extreme Gradient Boosting (XGBoost). These were selected for their robust performance, wide application in health-related research, and ease of interpretation.19,20 The Ridge LR model served as a baseline to assess the consistency of findings and assess whether advanced models (RF and XGBoost) can improve predictive performance by capturing nonlinear relationships between features. The model was implemented using: sklearn.linear_model LR library, scran ensemble for RF, and xgboost Library for XGBoost. To maintain class balance in both training and testing subsets, the dataset was stratified by the outcome variable and split using an 80:20 ratio (80% for training, 20% for testing). Due to the exploratory nature of this predictive modeling study, no formal sample size calculations were performed. However, we followed recommended ML practices and used a training set above the generally accepted 10-20 EPV threshold to reduce overfitting, ensuring sufficient outcome events for each predictor variable (EPV).twenty one Feature scaling using standard scaler It was only applied to the LR model to ensure proper convergence and interpretation of the coefficients. No scaling was applied to the RF and XGBoost models, as these tree-based algorithms are not sensitive to feature scaling. To address class imbalance, the RF and LR models include class weight = “Balanced”The parameters assign weights that are inversely proportional to class frequencies, ensuring fairer model predictions.twenty two For XGBoost,Scale position weightThe parameters were set based on the ratio of positive class samples to negative class samples, which helped the model to properly focus on the minority class.twenty three Hyperparameter tuning was performed using a 10-fold cross-validation approach on the training data. In each iteration, 90% of the data was used for training and 10% for validation.twenty four A grid search strategy was applied to identify the optimal model configuration. Bootstrap with 1000 resamples (n= 1000) to estimate the stability and variability of the model’s performance metrics.twenty five The model performance was evaluated using the area under the receiver operating characteristic curve (AUC), which is an important metric for evaluating classifier performance in situations of class imbalance.26
Feature selection was informed by importance scores derived from SHapley Additive exPlanations (SHAP), based on the method proposed by Lundberg SM, Erion GG, and Lee SI.27 Unlike traditional approaches such as gain-based importance used in tree-based models, SHAP provides a more robust and consistent estimate of the influence of each feature. Based on game theory principles, SHAP attributes a value to each feature that quantifies its unique contribution to the model’s predictions, providing a fair and interpretable representation of the feature’s influence.
Evaluate and interpret model performance
Since the test data were from the same dataset, the final evaluation served as an internal validation. Once the best performing model was identified, its effectiveness was tested using 20% of the reserved data that was not included in training or validation. This testing phase evaluated the model’s ability to generalize to new and unknown data. Various performance metrics such as accuracy, sensitivity, and specificity were used to provide a comprehensive evaluation of the model.28 Youden index was applied to determine the optimal threshold for classification.29 Additionally, the AUC was calculated to evaluate the discriminative power of the model across different thresholds.30 Bootstrap with 1000 resamples (n= 1000) to estimate the stability and variability of the model’s performance metrics.
To improve model interpretability during internal validation, we used SHAP values to assess both feature importance and the contribution of individual features to model predictions. Two types of SHAP visualizations were utilized. Dependency plots, which show the influence of specific features on the model results, and summary plots, such as bar graphs and honeycomb plots, that provide an overview of feature importance. These visualizations were generated using a test dataset to evaluate the role and impact of each feature.
