Modelling Nurses’ Intention to Leave the Hospital and the Profession using a Bivariate Additive Ordered Model via Penalized Likelihood with the R Package Pblm Article Swipe

IntroductionAddressing the shortage of healthcare workers requires a clear understanding of the factors associated with nurses’ intention to leave their current hospital (ITL1), or more critically, to leave the healthcare profession altogether (ITL2). Univariate models, which analyze these outcomes separately, often fail to account for their dependence, resulting in reduced estimation precision and a higher risk of Type II errors. Modelling their joint behavior is essential to improve estimation efficiency and reduce false negatives. Moreover, explicitly capturing the association structure enables the identification of discordant profiles, such as individuals intending to leave the hospital but not the profession, or vice versa. ObjectivesThis paper investigates the determinants of nurses’ intention to leave the hospital and the profession using a bivariate additive ordered logit model, emphasizing its ability to model complex dependence structures in ordinal categorical data. The analysis is implemented using the new R package pblm. MethodsThe data derive from the METEOR [1] cross-sectional survey conducted in 2022 in eight hospitals across Belgium, the Netherlands, Italy, and Poland. The METEOR Turnover Intention questionnaire (MTI), administered to nurses in these hospitals, was based on the Job Demands–Resources (JD-R) model. It included validated instruments measuring job satisfaction, work engagement, burnout, and turnover intentions (ITL1 and ITL2), both assessed on five-level Likert scales, along with individual and hospital-level covariates. A previous analysis [2] addressed these outcomes separately. In contrast, this study applies a bivariate additive ordered logit model [3] with an association structure governed by a penalty term [4], which constrains the association intercepts (log-global odds ratios, log-gOR) to follow a data-driven polynomial structure. P-splines are used to model non-linear effects of age in both marginal and association equations. The model is fitted using the R package pblm, soon to be released on CRAN. ResultsThe survey collected 1350 complete responses. In the marginal model for ITL1, significant factors included younger age (gOR = 0.95, p < 0.001), having experienced bullying (gOR = 1.31, p = 0.040), emotional exhaustion (gOR = 2.24, p < 0.001), low opportunities for professional development (gOR = 1.80, p = 0.022), low support from supervisors (gOR = 2.10, p < 0.001), low work prospects (gOR = 2.20, p < 0.001), poor physical working conditions (gOR = 1.30, p = 0.038), underuse of professional abilities (gOR = 1.64, p = 0.001), and low salary (gOR = 1.50, p < 0.001). The Netherlands showed the highest country effect (gOR = 1.84, p < 0.001), while Italy showed the lowest (gOR = 0.55, p = 0.005). In the marginal model for ITL2, significant factors included younger age (gOR = 0.96, p < 0.001), experiences of bullying (gOR = 1.31, p = 0.043), emotional exhaustion (gOR = 2.16, p < 0.001), depersonalization (gOR = 2.24, p < 0.001), work-life conflict (gOR = 1.56, p < 0.001), low professional development (gOR = 1.76, p = 0.026), low supervisor support (gOR = 1.56, p = 0.012), low work prospects (gOR = 1.86, p < 0.001), poor physical working conditions (gOR = 1.48, p = 0.002), underuse of professional abilities (gOR = 1.38, p = 0.030), low salary (gOR = 1.62, p < 0.001), and low overall job satisfaction (gOR = 1.56, p = 0.020). Again, the Netherlands exhibited the strongest country effect (gOR = 1.46, p = 0.001), while Italy the weakest (gOR = 0.28, p < 0.001). Regarding the association between ITL1 and ITL2, the strength of association increased significantly with age (relative gOR = 1.046, p < 0.001) and high working pace (relative gOR = 1.622, p < 0.001). Conversely, the association decreased significantly in the presence of health problems (relative gOR = 0.60, p = 0.021), low work prospects (relative gOR = 0.35, p < 0.002), underuse of professional abilities (relative gOR = 0.52, p = 0.025), and when working in the Netherlands (relative gOR = 0.50, p = 0.004). Figure 1 displays the observed and predicted log-gOR structures, along with the partial effects of age (centered at 22 years) estimated using P-splines. These indicate a non-linear increasing effect of age in both marginal and association equations. ConclusionsThe use of a bivariate additive ordered logit model allowed for the identification of a wider set of significant predictors for nurses’ intention to leave the hospital or the profession, compared to previous univariate analyses [2]. Specifically, five additional factors were identified for ITL1 (bullying, low professional development, low supervisor support, poor physical conditions, and low salary) and six for ITL2 (bullying, high working pace, work-life conflict, low supervisor support, poor physical conditions, and salary). These findings offer more comprehensive insights into the drivers of nurses’ turnover intentions and may support the development of more targeted retention strategies. Furthermore, the association model highlighted the presence of discordant profiles—nurses whose characteristics are associated with a weaker connection between ITL1 and ITL2. These insights may guide the design of future surveys or the refinement of questionnaire items. For instance, among more experienced nurses, the intention to leave the hospital and the profession may reflect a single underlying construct, whereas among younger nurses these outcomes appear more distinct.