Abstract Aim This study aimed to identify the risk factors for gingival invagination during orthodontic treatment after premolar extraction. Materials and Methods The medical records of 135 patients who had undergone interdental space closure after premolar extraction were collected, and cone beam computed tomography was performed to determine the presence of gingival invagination. The risk factors were examined using mixed‐effects models and generalized propensity score weighting (GPSW) to develop a predictive model. Results Univariate analysis revealed that the extraction site, buccal bone thickness 4 mm apical to the cemento‐enamel junction (MB1), mid‐root buccal bone thickness (MB2) and vertical skeletal relationships were related to gingival invagination ( p < .05). Furthermore, a subsequent multivariable mixed‐effects model analysis indicated a significantly increased risk of gingival invagination at MB1 < 1 mm ( p < .001; odds ratio [OR MB1≤0.5mm ] = 29.304; 95% confidence interval [CI]: 8.986–93.807; OR 0.5<MB1<1mm = 22.309; 95% CI: 9.890–50.320). Row mixing model analysis performed after balancing covariates using GPSW showed that MB1 and vertical skeletal relationships were associated with gingival invagination ( p < .05) with an increase in ORs. Therefore, in addition to MB1, the risk of gingival invagination is higher in hyperdivergent and normodivergent patients. During the establishment and internal validation of the predictive model, the area under the curve for all three models exceeded 0.7. Conclusions The risk of gingival invagination is higher in patients with MB1 < 1 mm and in normodivergent or hyperdivergent patients.