This retrospective cohort study used data from couples who underwent IVF or ICSI at the Reproductive Medical Center of Shanghai First Maternity and Infant Hospital between January 2019 and January 2025. This study was approved by the Research Ethics Committee of Shanghai First Maternity and Infant Hospital (KS25468). The inclusion criteria were: (1) a primary diagnosis of male factor infertility, including oligozoospermia, asthenozoospermia, or teratozoospermia, defined according to the WHO 5th edition criteria (22); (2) treatment with conventional IVF or ICSI; (3) complete follow-up records for pregnancy outcomes; (4) only the first IVF/ICSI treatment cycle was included for each couple. The exclusion criteria were: (1) severe uterine malformations or intrauterine adhesions in the female partner; (2) chromosomal karyotype abnormalities in either partner; (3) cycles involving donor sperm or oocytes; and (4) missing values in non-imputable administrative or eligibility variables (none in the final analytic cohort). Ultimately, 2,565 couples were included. To enhance the methodological rigor and better reflect real-world clinical application, the cohort was partitioned strictly according to the date of embryo transfer (23). The earliest 70% of cases (n = 1,797) were assigned to the training set, and the most recent 30% (n = 768) were reserved as an internal validation set. The training set was used exclusively for feature selection, model development, and hyperparameter tuning, whereas the validation set was used only for final model evaluation.

Variables were extracted from the electronic medical record (EMR) system. Demographic characteristics included female age, male age, female body mass index (BMI), male BMI, infertility duration, infertility type, female education, and male education. Clinical characteristics included menstrual regularity; basal follicle-stimulating hormone (FSH), luteinizing hormone (LH), estradiol (E2), progesterone (P), testosterone (T), prolactin (PRL), and anti-Müllerian hormone (AMH). The outcome was clinical pregnancy, defined as the presence of a gestational sac with fetal cardiac activity in the uterine cavity on transvaginal ultrasound 28–35 days after embryo transfer (24). Absence of a gestational sac or biochemical pregnancy was classified as non-clinical pregnancy.

For missing data, multiple imputation was performed using the mice package in R (25). In this dataset, all candidate predictors exhibited low levels of missingness (<5% in both the training and validation sets; Supplementary Table S1 and Supplementary Figure S1) (25). No participants were excluded due to missing non-imputable administrative or eligibility-defining variables (e.g., outcome follow-up), and therefore the final analytic cohort contained complete information on all eligibility-defining variables. Given the low proportion of missingness across predictors, multiple imputation was applied to all predictors to avoid unnecessary case deletion while minimizing potential instability associated with highly incomplete variables. Five imputed datasets were generated (m = 5, seed = 123). To avoid information leakage, imputation was conducted after the temporal split and performed separately within the training set and the temporal validation set. In the training set, the imputation model included all candidate predictors and the outcome variable in order to preserve predictor–outcome associations for model development. In the validation set, the imputation model included only predictors and explicitly excluded the outcome variable, thereby preventing outcome-informed imputation during model evaluation. Predictive mean matching was used for continuous variables and logistic regression for binary variables, with 20 iterations per imputation. The predictor matrix followed the default mice setting in which predictors were allowed to inform each other, except where structurally inappropriate; specifically, the outcome variable was excluded from all validation-set imputation models, and administrative or eligibility-defining variables were not imputed. The full imputation methods and predictor matrix for the training set are reported in Supplementary Table S2 and Supplementary Figure S2. Feature selection was performed using the least absolute shrinkage and selection operator (LASSO) within the training set only (26). LASSO was implemented in R using the glmnet package with internal standardization (standardize = TRUE). LASSO was fitted separately within each of the five imputed training datasets, and predictors were retained in the final feature set if they showed stable selection at λ₁SE across imputations, operationalized as non-zero coefficients in at least four of the five imputed datasets (27). The penalty parameter λ was selected using 10-fold cross-validation, and the minimum deviance occurred at λ_min = 0.00406, while the 1-standard-error criterion selected λ_1SE = 0.01801. For downstream machine learning analyses, continuous predictors were standardized to Z-scores using the StandardScaler function in the Python scikit-learn library (28), with scaling parameters learned from the training set and then applied to the validation set. Model development proceeded separately within each imputed training dataset. Each fitted model was then applied to each imputed validation dataset to generate predicted probabilities; for each individual, predicted probabilities were averaged across imputations to obtain pooled predictions. All performance metrics were calculated using these pooled predicted probabilities.

Seven commonly used machine learning algorithms were developed to predict pregnancy outcomes: logistic regression (LR), decision tree (DT), random forest (RF), support vector machine (SVM), artificial neural network (ANN), XGBoost, and LightGBM. Hyperparameters for each model were optimized using 5-fold cross-validation with grid search within the training set. To ensure comparability and avoid optimistic bias under multiple imputation, hyperparameter tuning was performed only once using the first imputed training dataset, with a fixed random seed (seed = 123) (29). The resulting optimal hyperparameters were then held constant and applied to all five imputed training datasets for model fitting. This strategy ensured that model complexity and tuning degrees of freedom were consistent across imputations while allowing uncertainty due to imputation to be reflected in model estimation.

Model performance was evaluated in the validation set using the area under the receiver operating characteristic curve (AUC), accuracy, sensitivity, specificity, positive predictive value (PPV), negative predictive value (NPV) and F1 score. For all classification metrics, predicted probabilities were dichotomized using a fixed operating threshold of 0.5, which was applied consistently across all models and datasets to ensure comparability.

Calibration was evaluated using the Brier score (30) and calibration curves constructed by grouping predicted probabilities into deciles and plotting observed versus predicted outcome probabilities. Clinical utility was examined using decision curve analysis (DCA) over a clinically plausible threshold probability range of 20–80%, reflecting the range in which clinicians may reasonably consider counseling or intervention in the context of IVF/ICSI outcome prediction.

Models were fitted separately within each of the five imputed datasets. For each individual and each model, predicted probabilities were averaged across imputations to obtain a single pooled predicted probability. All discrimination metrics (including AUC), classification metrics, calibration analyses, and decision curve analyses were computed based on these pooled predicted probabilities.

The 95% confidence intervals for AUC were estimated using nonparametric bootstrap resampling (1,000 replications) applied to the pooled predicted probabilities within each dataset. Pairwise comparisons of AUCs between models were conducted using DeLong’s test for correlated receiver operating characteristic curves based on ROC curves constructed from the pooled predicted probabilities. Specifically, the AUCs of each model were compared against those of LightGBM and logistic regression, respectively, and the corresponding P values were reported.

The representative model was then interpreted using SHAP to quantify the marginal contribution of each feature (31). SHAP summary plots, beeswarm plots, and dependence plots were generated to illustrate the model’s decision patterns. SHAP values were computed using the standard interventional SHAP implementation provided by the SHAP Python package. We acknowledge that when predictors are correlated, feature attributions may be influenced by feature dependence, and therefore SHAP results should be interpreted as model-based associations rather than causal effects.

All analyses were performed using R (version 4.2.0) and Python (version 3.9.0). In R, data preprocessing, multiple imputation, descriptive analyses, dataset partitioning, and LASSO feature selection were conducted using standard statistical packages, including mice, caret, tableone, and glmnet. All machine learning model development, hyperparameter tuning, and performance evaluation were conducted in Python. Graphical analyses, including ROC curves, calibration curves, and decision curve analysis, were generated using matplotlib and ggplot2.

Baseline characteristics of the total cohort, training set, and validation set are presented in Table 1. The distributions of outcomes, demographic variables, and clinical characteristics were highly comparable between the training and validation sets. All standardized mean differences (SMDs) were below 0.25, with most variables showing SMDs < 0.10, indicating good balance between the two datasets. Slight imbalances were observed for female age (SMD = 0.218) and AMH (SMD = 0.113), but the overall clinical characteristics remained broadly similar across the two subsets, supporting the appropriateness of the temporal split for model development and validation.

LASSO regression was applied to reduce multicollinearity and identify key predictors. As shown in Figures 1A, B, the penalty parameter was selected using 10-fold cross-validation. The minimum cross-validated deviance occurred at λ_min = 0.00406, while the 1-standard-error criterion selected λ_1SE = 0.01801, yielding a parsimonious set of five predictors with non-zero coefficients. Across the five imputed training datasets, five predictors were consistently selected by LASSO at λ₁SE in at least four imputations: female BMI, male BMI, basal FSH, AMH, and female age. Variables such as menstrual regularity, infertility type, and education level were not stably selected and were therefore excluded. This stability-based selection suggests that spousal BMI and ovarian reserve markers constituted the most robust predictors within the available feature set.

The predictive performance of the seven models in the training and validation sets is summarized in Tables 2 and Supplementary Table S3, with ROC curves, calibration plots, and decision curve analyses shown in Figure 2, and confusion matrices presented in Figure 3. In the training set, ensemble models (Random Forest, XGBoost, and LightGBM) achieved higher AUCs than other algorithms (AUCs 0.903–0.923), whereas model performance became more comparable in the temporal validation set, with AUCs ranging narrowly from 0.840 to 0.857 across all models (Table 2, Figure 2B). LightGBM, XGBoost, and Random Forest demonstrated almost identical discriminative ability (all AUC = 0.857). DeLong tests confirmed that there were no statistically significant differences in AUC between LightGBM and Random Forest (P = 0.918), XGBoost (P = 0.985), or Logistic Regression (P = 0.067).

Beyond discrimination, LightGBM showed a balanced performance profile in the validation set, with an accuracy of 0.775, high specificity of 0.909 and moderate sensitivity (0.596). Calibration analysis suggested reasonable agreement between predicted and observed risks (Brier score = 0.145; Figure 2F), and decision curve analysis indicated net clinical benefit across a range of clinically plausible threshold probabilities (Figure 2D). Taken together, LightGBM was selected as the representative model for subsequent interpretability analyses due to its overall stability across discrimination, calibration, and clinical utility, rather than on statistically superior AUC alone.

Analysis of confusion matrices (Figure 3) further illustrated model behavior. LightGBM achieved a high true-negative rate (92.20%), reflecting strong specificity, while maintaining a sensitivity comparable to other models. From a clinical perspective, this tendency to limit false-positive predictions may be advantageous for avoiding overly optimistic prognostic assessments in couples with low likelihood of pregnancy.

To enhance interpretability of the selected model, SHAP was applied to visualize the contribution of individual predictors to model output (Figure 4). The SHAP bar plot (Figure 4A), based on mean absolute SHAP values, indicated that male BMI and female BMI showed the highest average contributions within the fitted model and the available feature set, followed by basal FSH, AMH, and female age. The SHAP beeswarm plot (Figure 4B) further illustrated the direction of these associations: higher BMI values in either partner (red points) were predominantly located on the negative side of the x-axis, suggesting that higher BMI was associated with lower predicted probability of clinical pregnancy in the model output.

The SHAP force plot (Figure 4C) presents an illustrative individual case, showing how each feature contributed to shifting the prediction from the baseline toward a lower probability of pregnancy. In this example, elevated male BMI and older female age exerted negative contributions that outweighed the positive contribution of AMH. This visualization demonstrates how the model integrates multiple features to generate a personalized prediction, while reflecting model behavior rather than biological causation.

Finally, the SHAP dependence plots (Figure 4D) suggested nonlinear relationships between predictors and model output. Both female BMI and male BMI showed a threshold-like pattern: SHAP values remained relatively neutral within the lower range but declined sharply once BMI exceeded approximately the upper-normal range. AMH demonstrated a modest positive association at low-to-moderate levels, while higher FSH and increasing female age were associated with progressively negative SHAP values. These patterns reflect how the fitted model utilizes these predictors and should be interpreted as model-based associations rather than evidence of specific biological thresholds.