Explicating the interplay of maternal work stress, breastfeeding, mental health, and infant attachment using machine learning and statistical modeling

2.1 Data acquisition

This study used the datasets that employed a cross-sectional design to examine the effects of maternal work stress on breastfeeding behaviors and infant attachment. Publicly available datasets from the Open Science Framework (OSF) (https://osf.io/58dng) were utilized, alongside a master dataset containing 2,010 maternal and child health indicator observations (specifically the dataset curated by Braithwaite, Pearson, and Wright (updated May 6, 2025)). Data collection for the original dataset occurred between 2022 and 2024 through structured surveys administered to mothers recruited via community health centers and online parenting networks in the United Kingdom.

2.2 Data preprocessing

Comprehensive data pre-processing steps were applied to ensure analytical robustness and data integrity. Missing data were handled using median imputation for numerical variables and mode imputation for categorical ones (Gautam & Latifi, 2023). According to the inclusion criteria adopted for the original study, mothers were required to have at least one child under 5 years of age and be engaged in full-time or part-time work or be self-employed. The data of participants with missing maternal status or breastfeeding data were excluded to maintain the integrity and validity of the statistical and machine learning analyses. Categorical variables with more than 50 percent missing entries were eliminated to maintain statistical validity. Outliers for work stress levels and breastfeeding duration were detected using the interquartile range (IQR) method and handled accordingly (López-Fernández et al., 2023). To prevent the leaking of data, median imputation was done for the numerical variables and mode imputation for the categorical variables. These imputations were both trained on the training data and applied to the test data. Numerical features were scaled with StandardScaler, and categorical features were encoded with ​‍​‌‍​‍‌​‍​‌‍​‍‌OneHotEncoder. The Edinburgh Postnatal Depression Scale (EPDS) score refers to item 10 (suicidal ideation, range 0-3). Whereas, Generalized Anxiety Disorder 7-item (GAD-7) values reported are item-level means (0-3); total scores range from 0-21. Work stress was measured on a standardized 1-5 Likert scale. Earlier mentions of a 0-60 range reflected raw summed scores from item banks; herein, consistent reporting of the calibrated 1-5 scale across all analyses, figures, and text to ensure clarity and comparability was ensured. Biologically implausible values (e.g., negative ages, negative breastfeeding durations) were removed before inferential tests and machine learning analyses.

2.3 Statistical analysis

The exploratory data analysis (EDA) summarized the dataset’s features and identified patterns of interest. Descriptive statistics were computed for numerical variables such as maternal age, work hours per week, and breastfeeding duration. Frequency distributions were analyzed for categorical variables, including breastfeeding exclusivity, maternal employment type, and job stress levels. To assess the relationship between work stress and breastfeeding outcomes, correlation analysis was performed using Pearson and Spearman coefficients. Bivariate analyses, including chi-square tests, t-tests, and one-way Analysis of Variance (ANOVA), were used to examine variations in breastfeeding behaviors across employment and stress levels. Additionally, principal component analysis (PCA) was utilized to reduce dimensionality and uncover key patterns in high-dimensional data, allowing for a more refined understanding of factors influencing maternal work stress and breastfeeding behaviors. These statistical procedures were used to establish baseline associations and distributional properties, which informed variable selection and model specification for subsequent regression and machine learning analyses.

Regression analysis assessed the relationship between maternal work stress and breastfeeding behaviors. Breastfeeding duration was analyzed separately through descriptive statistics and classification models. Linear regression models estimated the impact of work hours and stress levels on breastfeeding duration. Bivariate regression analyzed interactions between maternal psychological well-being and infant attachment scores. Various regression models accounted for confounding variables that included socioeconomic status, education, and support networks. ANOVA was used to compare stress levels across different employment categories. The statistical significance was computed using p values (≤ 0.05), and residual diagnostics were conducted to assess the model’s validity.

2.4 Machine learning analysis

Supervised machine learning models were employed to predict breastfeeding behaviors and assess the impact of maternal work stress on infant attachment. The classification target was defined as breastfeeding continuation beyond 6 months (Yes/No). Classification models, comprising logistic regression, decision trees, and support vector machines (SVMs), were applied to categorize breastfeeding continuation beyond six months. Model selection relied on interpretability, robustness, and compatibility with the size of the dataset and the complexity of variables. Logistic Regression established a clear baseline, SVM effectively handled non-linear margins, and Random Forest excelled at handling feature interactions and missing values. eXtreme Gradient Boosting (XGBoost) was not included to ensure interpretability and computational efficiency, given the dataset’s modest size. Machine learning approaches were applied to capture (Fig. S2(a)) complex, non-linear interactions between maternal stress, sociodemographic factors, and breastfeeding outcomes that are not adequately modeled by conventional regression techniques.

The dataset was split into training (70%) and test (30%) sets using stratified random sampling, ensuring balanced representation of breastfeeding continuation categories. Furthermore, five-fold cross-validation was utilized on the training set to ensure model stability and reduce sampling bias throughout the classification process. A feature importance analysis was conducted to identify key predictors influencing breastfeeding and attachment outcomes. The Gini impurity criterion was used to assess feature importance in the Random Forest model, measuring each variable’s contribution to classification accuracy. To verify robustness, permutation importance was calculated to ensure reliable identification of predictors affecting breastfeeding continuation and maternal stress outcomes.

Hyperparameter tuning was applied through grid and random search over the following intervals: n_estimators (100-1000), max_depth (3-20), and min_samples_split (2-10). Five-fold cross-validation was applied to avoid overfitting and guarantee model generalizability. Average performance metrics were recorded for all folds to reflect consistency and reliability. Breastfeeding outcomes class imbalance (continuation vs. cessation) was handled by stratified sampling during data splitting and class-weight balancing in classification algorithms (Logistic Regression, SVM, Random Forest). These alterations guaranteed that both outcome categories received the same amount of input and minimized the model’s bias towards the majority class. Performance metrics were used for model evaluation, including accuracy, precision, recall, F1-score, and area under the receiver operating characteristic (ROC) curve. Moreover, unsupervised learning techniques, such as k-means clustering, were utilized to profile mothers’ data according to stress levels and breastfeeding behaviors. Statistical analyses were performed with Python (v3.10) and the libraries Pandas, Scikit-learn, NumPy, and Matplotlib. Seaborn and Plotly tools were used for data visualization. The Scikit-learn pipelines were used to implement and validate the machine learning models.

3.1 Exploratory data analysis

Data analysis explored maternal ages, with an average of 35.17 years, a median of 34, and a mode of 33. Symmetric distribution of child age was around 5 years, with an interquartile range of 3-7 years (Fig. S1a). Maternal education (1 = No formal education, 2 = Primary, 3 = Lower secondary, 4 = Upper secondary, 5 = Undergraduate, 6 = Postgraduate or Professional degree) is evenly distributed, with a spike at the highest level, suggesting broad socioeconomic representation (Fig. S1b). Gestational age clustered around 39-40 weeks with some preterm births (Fig. S1c), while delivery methods (1 = Normal vaginal, 2 = Assisted vaginal, 3 = Elective cesarean, 4 = Emergency cesarean, 5 = Other or Unspecified) were evenly divided (∼18-21% each), indicating heterogeneity in delivery experiences (Fig. S1d). Table 1 presents baseline characteristics by stress levels. Fig. 2(a) compares birth weight by sex, showing similar 2-5 kg ranges with slightly narrower distribution in one group. The pair-plot (Fig. 2b) shows near-normal distributions for key variables, with only gestational age and birth weight showing a mild positive association. Breastfeeding duration (Fig. 2c) peaks around 15-18 months. Correlation analysis (Fig. 3a) shows weak associations overall, with gestational age–birth weight correlation at r≈0.25. Breastfeeding duration fluctuates without a trend in Fig. 3(b). Whereas, Fig. 3(c) shows a noisy but positive gestational age-birth weight trend. Mental health impacts (Fig. 3d) are evenly distributed across three discrete levels. Sex (child) (Fig. 3e) is evenly distributed, and Fig. 3(f) again highlights child age anomalies. Joint distribution of gestational age and birth weight (Fig. 4a) clusters at term and average weight. Breastfeeding duration rises cumulatively to ∼40 months (Fig. 4b), peaking at 0-30 months (Fig. 4c; Fig. S2(a)). Maternal age peaks near 30 years (Fig. 4d). Delivery methods show clustering in modes 2 and 4 (Fig. 4e). Birth weights cluster between 2.8-4.0 kg (Fig. 4f). Mean birth weights differ little by delivery method Supplementary Fig. S2(b). Breastfeeding difficulty shows minor influence on mental health, Fig. S2(c), while radar plots indicate balanced mental health dimensions, Fig. S2(d). Breastfeeding duration weakly increases with maternal education, Fig. S2(e). Gestational age distributions vary slightly by delivery method, as depicted in Fig. S2(f).

Close

Fig. 2.

(a) Strip plot of birth weight (kg) by Sex (child). (b) Pairplot of maternal and child variables with histograms and scatter plots. (c) KDE plot of breastfeeding duration (months).

Export to PPT

Close

Fig. 3.

(a) Correlation heatmap of maternal, child, and perinatal variables. (b) Line plot of breastfeeding duration over the sample index. (c) Scatter plot of birth weight versus gestational age. (d) Strip plot of mental health impact score. (e) Bar plot of Sex (child) count. (f) Rug plot of child age in months.

Export to PPT

Close

Fig. 4.

(a) Joint plot of gestational age and birth weight. (b) Step plot of breastfeeding duration. (c) Histogram of breastfeeding duration. (d) Density plot of maternal age. (E) Donut chart of delivery method. (f) Box plot of birth weight.

Export to PPT

3.2 Inferential analysis

ANOVA tests highlight key factors affecting maternal and child well-being. Fig. 5(a), a dumbbell chart, shows child birth weights pre- and post-cleaning; most values fall between 2.5-4.0 kg, with reduced differences after preprocessing, underscoring data cleaning’s importance. Fig. 5(b), a stacked bar plot, shows that most mothers reported breastfeeding support, a potential confounder in infant feeding and maternal health analyses. Maternal age distribution (Fig. 5c), displayed as a circular bar chart, peaks between 30-35 years, confirming earlier histograms while enhancing symmetry visualization. The hexbin plot (Fig. 5d) of gestational age and birth weight is a refinement of Fig. 4(a), showing dense clustering near 38-40 weeks and 3.0–3.5 kg, consistent with term births. Fig. 5(e), a quantile-quantile (Q-Q)-like plot, indicates that the distribution of breastfeeding duration is close to a normal distribution with a slight skew in the tails. Stress levels (Fig. 5f) deviate from normality, displaying ordinal clustering, suggesting non-parametric or ordinal regression may be suitable. Boxplots in Fig. 6(a) reveal breastfeeding duration medians are consistent across stress levels, though variability exists. Fig. 6(b), a pie chart, shows that stress levels are evenly distributed (∼19-21%), reducing class imbalance. Fig. 6(c) groups breastfeeding durations, showing that extended breastfeeding (12-24 months, 40.8%) dominates. Scatterplots (Fig. 6d) and correlations (Fig. 6e, r≈0.015) confirm that the relationship between stress and breastfeeding is minimal. These findings suggest that while direct inferential associations are weak, more complex interaction effects may exist, warranting the use of multivariate and machine learning approaches to uncover hidden patterns. Overall, results suggest that maternal stress modestly influences breastfeeding, and supportive measures for mothers are necessary to address this issue.

Close

Fig. 5.

(a) Dumbbell chart of birth weight pre-/post-cleaning. (b) Stacked bar plot of access to breastfeeding support. (c) Circular bar chart of maternal age distribution. (d) Hexbin plot of gestational age vs. birth weight. (e) Q-Q plot of breastfeeding duration. (f) Q-Q plot of work stress level.

Export to PPT

Close

Fig. 6.

(a) Boxplot of breastfeeding duration by work stress levels; (b) Pie chart of work stress level distribution; (c) Breastfeeding duration category percentages; (d) Scatterplot of work stress vs. breastfeeding; (e) Heatmap showing weak correlation values.

Export to PPT

3.3 Regression analysis: Model evaluation and interpretation

This study compared Linear, Multiple, and Polynomial Regression models for forecasting birth weight. The Polynomial Regression model performed (Mean Squared Error (MSE): 183,037.66; R²: 0.2377) better than Linear and Multiple Regression (both MSE: 189,951.97; R²: 0.2089), reflecting the model’s improved capability to capture the non-linear relationships. However, the overall power of the model to explain the variation in the data was still quite low. Fig. 7(a) shows that MSE is lowest for Polynomial Regression, while Fig. 7(b) confirms higher R². Residual analysis revealed heteroscedasticity in Linear Regression (Fig. 7c) and slight quantile-quantile (Q-Q) plot abnormalities in Multiple Regression (Fig. 7d). Polynomial Regression showed improved fit, though performance weakened for high birth weights. The Multiple Regression model yielded a highly significant F-statistic (p < 0.0001). Maternal education emerged as a strong positive predictor, whereas maternal age showed no significant effect. Predicted vs. actual plots (Fig. 7e) demonstrated Polynomial Regression aligning most closely with observed values. Multicollinearity, identified by high variance inflation factors (VIFs) and near-zero eigenvalues, indicated coefficient instability, leading to dimension reduction (PCA) where VIF > 5. Overall, Polynomial Regression provided superior predictive accuracy but also made apparent the limitations of the data and the modeling approach.

Close

Fig. 7.

Regression model diagnostics for birth weight prediction: (a) MSE comparison, (b) R² values, (c) residuals for linear regression, (d) Q-Q plot for multiple regression residuals, (e) predicted vs. actual birth weight using polynomial regression.

Export to PPT

3.4 Machine learning model performance

Model comparison showed Polynomial Regression outperforming Linear Regression slightly, with lower MSE and higher R². Fig. 8(a) confirms that Polynomial Regression minimizes error more effectively, while Fig. 8(b) shows consistently higher R² values, indicating breastfeeding duration may follow non-linear patterns. Classification models trained on a 70/30 split further highlighted predictive differences. Logistic Regression (Accuracy = 75%, Precision = 70.2%, Recall = 73.5%, F1 = 71.8%), Support Vector Machine (Accuracy = 75%, Precision = 72.4%, Recall = 74.1%, F1 = 73.2%), and Random Forest achieved the best performance (85% accuracy; Precision 77.8%, Recall 87.5%, F1 82.3%), outperforming Logistic Regression and SVM (Fig. 8c). Predicted vs. Actual plots (Fig. 8d) show Polynomial Regression tracking closely to the 45° line, however, the dispersion indicates that there is some degree of individual variability. This performance indicates that breastfeeding cessation is driven by interacting risk factors rather than single predictor, supporting the clinical relevance of multivariable, data-driven risk stratification. The Confusion Matrix (Fig. 8e) demonstrates high sensitivity and specificity, with most cases correctly classified and minimal false positives/negatives. Unsupervised clustering (Fig. 8f) revealed three distinct subgroups based on maternal work stress and breastfeeding duration, patterns not captured by regression or classification alone. These clusters suggest the presence of hidden heterogeneity, indicating distinct maternal experiences which in turn imply the necessity of tailored intervention programs. Polynomial Regression and Random Forest provided superior predictive performance, while the clustering method aided in identifying the latent structures, and thus deepened our understanding of the complex, multi-faceted nature of the relationships between maternal stress and breastfeeding behavior. The cross-validated performance of the model for categorizing maternal breastfeeding outcomes has been presented in Fig. 9(a), which features Receiver Operating Characteristic (ROC) curves for 5-fold cross-validation, illustrating the performance of each fold (Fold 1-5) along with the average ROC. Area Under the Curve (AUC) values reflect consistently robust discrimination, with an average AUC of 0.87. Precision-Recall (PR) curves (Fig. 9b) for 5-fold cross-validation, illustrating individual fold curves alongside the average PR. The average AUC = 0.89, indicating strong accuracy at different recall thresholds. These metrics suggest that the Random Forest classifier surpasses both Logistic Regression and SVM, achieving a strong predictive performance while maintaining a balance between sensitivity and specificity. Cross-validation underscores model consistency across folds, reinforcing trustworthiness in forecasting breastfeeding duration based on maternal stress and associated factors.

Close

Fig. 8.

(a) Bar plot comparing MSE of linear and polynomial regression; (b) R² score comparison; (c) Classification accuracy of logistic, random forest, and SVM models; (d) Predicted vs. actual breastfeeding duration; (e) Confusion matrix; (f) Work stress–breastfeeding duration clustering.

Export to PPT

Close

Fig. 9.

Five-fold cross-validated ROC and PR curves showing consistent Random Forest model performance (mean ROC-AUC = 0.87; PR–AUC = 0.89).

Export to PPT