Purpose <p>This study aimed to identify and quantify the contribution of various features and their interactions to the prediction of return to work (RTW) following medical rehabilitation due to musculoskeletal disorders, using a machine learning approach.</p> Methods <p>This retrospective cohort study used data from the Rehabilitation Statistics Database (RSD) of the German Pension Insurance. A total of 685,890 individuals who received multimodal medical rehabilitation for musculoskeletal disorders between January 2018 and December 2021 were included. A Light Gradient Boosting Machine (LightGBM) model was trained on 75% of the sample and tested on the remaining 25%. Model predictions were explained using Shapley values to determine the contributions of individual features and their interactions.</p> Results <p>Overall, 65.0% of individuals achieved stable RTW. On the test dataset (threshold 0.5), the model achieved an accuracy of 0.815, precision of 0.830, recall of 0.901, and an AUC-ROC of 0.867. The most important predictors of RTW were days in employment, income, incapacity for work at admission, occupational status, duration of incapacity for work, and age. The model also identified interaction effects, including an interaction between postacute rehabilitation and incapacity for work at admission, and an interaction between age and income.</p> Conclusion <p>The machine learning model demonstrated good predictive performance for RTW after medical rehabilitation for musculoskeletal disorders and confirmed the relevance of established predictors, such as prior employment and incapacity for work. In addition, the analysis identified interactions that are difficult to detect with traditional regression approaches.</p>