The study obtained approval from the ethics committee of the First Medical Center of Chinese PLA General Hospital (No. 2025–071). As a single-center retrospective study, it waived the requirement for informed consent forms.

Imaging data were retrospectively collected from patients who underwent cranial and facial MRI scans for meningioma at the First Medical Center of Chinese PLA General Hospital between January 2014 and December 2024. All cases were stratified by age group (youth group: ≥ 18 and < 35 years; middle-aged group: ≥ 35 and < 60 years; senior group: ≥ 60 years) and randomly divided into training and test sets in a 7:3 ratio using stratified random sampling. A total of 237 patients were included (Figure 1), with 27 in the youth group, 59 in the middle-aged group, and 151 in the senior group (Table 1, Table 2).

All patients met the following inclusion criteria: (1) age ≥18 years; (2) no history of head or facial surgery, with primary disease excluding periorbital fat involvement; (3) absence of malignant tumors or immune system disorders; (4) no prior radiotherapy, chemotherapy, or glucocorticoid therapy; (5) no history of craniofacial dysplasia; (6) cranial and facial MRI scans with 1-mm slice thickness, covering a range from the skull vertex superiorly to the upper incisor plane inferiorly.

Patients were excluded based on the following criteria: (1) body mass index (BMI) < 18.5 or ≥ 28.0; (2) MRI images with poor resolution or indistinct boundaries, precluding accurate regions of interest (ROI) annotation; (3) incomplete coverage of target regions in cranial and facial MRI scans; (4) absence of T1-weighted sequences in cranial and facial MRI protocols; (5) history of facial or periorbital deformities and prior surgeries.

Clinical information and MRI images were collected retrospectively from the hospital’s electronic medical record system and imaging database. Patients diagnosed with meningioma by the neurosurgery department underwent preoperative high-resolution head and facial T1-weighted imaging, primarily with a slice thickness of 1 mm, meeting the study’s requirements. Identical MRI equipment and unchanging imaging parameters ensured data stability and reliability for this investigation.

Preoperative patients routinely underwent cranial and facial MRI scanning using high-resolution T1-weighted imaging. The scanning was performed with the following parameters: slice thickness of 1 mm, interslice gap of 1 mm, matrix size of 260 × 260, and an average of 2 signal acquisitions. The imaging was conducted on a 1.5 T high-field superconducting magnet (Siemens Espree, Erlangen, Germany) equipped with a 32-channel phased-array body coil.

Based on reviewing previous research data (19), three anatomical regions of periorbital fat were identified (Figure 2A, Figure 3): retro-orbicularis oculi fat (ROOF), sub-orbicularis oculi fat (SOOF) and deep medial cheek fat (DMCF). Two senior plastic surgeons (each with over 10 years of clinical experience) independently annotated the ROIs on all imaging data using 3D Slicer software (version 5.7.0). For complex or ambiguous images, annotations were guided by one expert plastic surgeon (with over 20 years of clinical experience). To ensure reproducibility, two senior plastic surgeons concurrently ambiguous imaging data from 50 randomly selected patients; an intraclass correlation coefficient (ICC) evaluation was subsequently performed (20). Metrics with ICC values below the reliability threshold (defined as ICC ≤ 0.75) were excluded as unstable indicators.

All craniofacial T1-weighted MRI images underwent N4 bias field correction algorithm using the Python (version v3.9.23) programming language to mitigate potential artifacts arising from local magnetic field in homogeneities. Subsequently, coordinate system standardization was performed on all imaging data to ensure accurate feature extraction.

The open-source library “PyRadiomics” (version v3.1.0) was employed to extract radiomics feature data from segmented ROIs using the Python programming language. Prior to feature extraction, all images were preprocessed: images underwent resampled to a uniform voxel size of 1 × 1 × 1 mm³, and grey-level data were discretized into 25 bins using nearest-neighbor interpolation. Ultimately, a total of 1,688 radiomics features were extracted, including 14 shape-based features, 324 first-order statistical features, 432 gray level co-occurrence matrix (GLCM) features, 288 gray level run length matrix (GLRLM) features, 288 gray level size zone matrix (GLSZM) features, 252 gray level dependence matrix (GLDM) features, and 90 neighboring gray tone difference matrix (NGTDM) features.

Consistent with conventional radiomics, the Python programming language was employed to extract machine learning features from segmented medical images. A 3D ResNet18 backbone network was utilized, optionally embedding squeeze-and-excitation (SE) modules after residual blocks. The SE mechanism compressed spatial information through global average pooling, generated channel-wise weights via fully-connected layers (compression ratio: 16), and recalibrated features using sigmoid activation to enhance discriminative channel responses. Input MRI scans (including images and masks) were uniformly resampled to 1 × 1 × 1 mm³ voxels. ROIs were extracted and cropped to minimum bounding boxes, with intensity values normalized to the 0–1 range using min-max scaling. ROI volumes were padded or cropped to 32 × 32 × 32 tensors (zero-padded for undersized volumes), retaining only masked regions. Features were extracted through global average pooling and fully-connected layers, yielding a 512-dimensional output vector without preserved spatial dimensions.

Using CR and ML features extracted from the training set, we conducted separate training procedures for nine different ML models (Figure 2B). Given the sample size and distribution characteristics of the study cohort, five-fold cross-validation with five repeats was employed to enhance model stability during validation. Procedure: (1) Features with ICC > 0.75 were retained due to high stability. All datasets were standardized using the Z-score method. (2) Analysis of variance (ANOVA) screened features exhibiting statistically significant differences (p < 0.05) across youth, middle-aged, and senior groups. (3) Pearson and Spearman correlation analyses were applied to normally and non-normally distributed features, respectively, to remove redundant features with correlation coefficients (r) > 0.9. (4) Synthetic minority over-sampling technique (SMOTE) was employed to balance sample sizes in the two non-senior subgroups, aligning them with the senior group (n = 105). (5) Random forest (RF) feature selection was implemented to retain only features with importance scores meeting or exceeding the mean value. (6) Least absolute shrinkage and selection operator (LASSO) regression was used for feature screening, with the penalty coefficient (λ) determined by minimum mean square error (MSE). (7) Principal component analysis (PCA) was synchronously applied to reduce dimensionality of selected features in the training, internal validation, and test sets. (8) nine ML models were trained on the training set, with performance rigorously evaluated on an independent test set.: multiclass logistic regression (MLR), neural network (NN), support vector machine (SVM), multilayer perceptron (MLP), random forest (RF), gradient boosting machine (GBM), light gradient boosting machine (Light GBM), naïve Bayes, and extreme gradient boosting (XGBoost).

MRI feature data from each patient were extracted from three ROIs (Figure 2C): Roof, Soof, and DMCF. These regions generated three distinct sets of feature data. Initially, the filtered CR features from these groups were fused through direct concatenation to construct a CR model (identical methodology was applied to the ML model). Subsequently, the screened features from CR and ML were fused through direct concatenation to construct a CR-ML fusion model. Finally, the predicted probabilities generated by the CR-ML fusion model were utilized as training data to build a stacking ensemble learning framework, with a logistic regression classifier employed as the meta-learner.

All data analyses were conducted using open-source libraries in Python. The ICC was employed to evaluate feature reproducibility, with ICC > 0.75 indicating good consistency. Model discriminative performance was assessed using the internal validation set and testing set through the following metrics: area under the curve (AUC), AUC macro-average (AUC-macro), 95% confidence interval (CI), accuracy (Acc), positive predictive value (PPV), sensitivity (Sen), F1-score, and confusion matrices. Owing to the imbalanced class distribution in the study sample, AUC-macro (calculated via interpolation) was selected as the primary evaluation metric due to its minimal susceptibility to sample distribution bias and superior stability (21). Based on the macro-average ROC curve, the AUC-macro (Macro_ AUC) is calculated using the trapezoidal rule:

f j and f j − 1 are adjacent false positive rate (FPR) grid points. macro _ tpr ( f j ) and macro _ tpr ( f j − 1 ) are the macro-average true positive rate (TPR) values at the corresponding FPR points. M is the total number of points in the FPR grid.

K is the total number of classes. tp r i ( f j ) is the true positive rate for the i-th class at FPR point f j ​, obtained through linear interpolation. f j is the j-th point in the FPR grid, where j = 0, 1, …, M , typically with f 0 = 0 and f M ​ = 1.

The F1-score, which integrates PPV and Sen, was designated a secondary metric as it may overestimate model performance in imbalanced data; other metrics served as supplementary indicators. To enable precise model comparison, all evaluation results were reported to three decimal places. For normally distributed data with homogeneity of variance, independent samples t tests and one-way ANOVA were applied; otherwise, non-parametric tests (Kruskal–Wallis H test) were utilized, with statistical significance defined as p < 0.05.

The baseline characteristics of the study participants are presented in Tables 1, 2. A total of 237 patients were enrolled: 27 in the young group (11 male), 59 in the middle-aged group (14 male), and 151 in the senior group (36 male). Through stratified sampling, the cohort was divided into a training set (including an internal validation set) of 165 cases: young group (n = 19, 28.5 ± 5.0 years, 7 male), middle-aged group (n = 41, 42.9 ± 4.7 years, 9 male), and senior group (n = 105, 60.0 ± 6.5 years, 26 male); along with a testing set of 72 cases: young group (n = 8, 28.6 ± 5.6 years, 4 male), middle-aged group (n = 18, 43.9 ± 4.1 years, 6 male), and senior group (n = 46, 58.8 ± 6.7 years, 10 male). Within identical age strata, no statistically significant differences in BMI were observed between the training and test sets (p > 0.05). Similarly, no statistically significant BMI differences existed across distinct age groups (p > 0.05).

Following feature selection and SMOTE oversampling for the Roof, Soof, and DMCF regions respectively, the feature dimensionality of the Roof and DMCF regions was reduced to 5 via PCA, while the Soof region was reduced to 8 features. Each patient thus contributed a total of 18 features. All models demonstrated stable performance on the training set (Table 3). In the test set, the Naive Bayes model exhibited optimal results with an AUC-macro of 0.757 (95% CI, 0.628–0.856) and an F1-score of 0.598.

Identical to the CR model, each patient contributed a total of 18 ML features. Within the training set evaluation metrics, all models exhibited stable performance, with results which were comparable to those of the CR model (Table 3). In the testing set, the optimal model was MLR, attaining an AUC-macro of 0.771 (95% CI, 0.651–0.862) and an F1-score of 0.504. The predictive capability of models developed using ML features was comparable to that of the CR model; the optimal model in the testing set showed an increase of 0.014 (1.8%) in AUC-macro but exhibited a decrease of 0.095 (−15.8%) in F1-score.

To further enhance model performance, we directly concatenated CR and ML feature data to develop a fusion model. In the evaluation results of the CR-ML fusion model (Table 3), the training set metrics remained stable. In the test set, the NN model performed optimally, achieving an AUC-macro of 0.819 (95% CI: 0.720–0.893) and an F1-score of 0.615. Following data fusion, the performance of most models (seven models) improved to varying degrees. Compared to the CR and ML models, the optimal model’s test set AUC-macro increased by 0.062 (8.2%) and 0.048 (6.2%), respectively, and the F1-score increased by 0.017 (2.8%) and 0.111 (22.1%), respectively.

A stacking ensemble learning model was developed using the probability data from both the training and test sets of the CR-ML fusion model, with a logistic regression model selected as the meta-learner. The base models produced nine sets of predicted probability data, and all possible combinations were exhaustively evaluated. As shown in Table 4, the evaluation metrics of the top 10 model combinations were compiled according to their test set AUC-macro ranking. In the training set, all models exhibited excellent performance; in the test set, the stacking model fusing three base models (GMB, Light GBM and NN) achieved optimal performance, with an AUC-macro of 0.833 (95% CI: 0.737–0.902), an F1-score of 0.614, an Acc of 0.597, and a positive predictive value (PPV) of 0.690. Compared to the optimal CR-ML fusion models (Table 5), the top-performing model demonstrated an increase in test set AUC-macro by 0.014 (1.7%), while other metrics remained relatively stable.