Prediction of bullying at work: A data-driven analysis of the Finnish public sector cohort study

Aim: To determine the extent to which change in (i


Introduction
Bullying at work refers to constant, repeated isolation of a member of the working community, belittling one's work effort, threats, talking behind one's back or other forms of pressure (Nielsen and Einarsen, 2018).Reports of workplace bullying are slightly more common among women than men.In 2020, for example, 10% of women and 9% of men working in the Finnish public sector reported workplace bullying.
There are at least two broad categories of factors that may be assumed to predict bullying at work.The work environment hypothesis suggests that bullying is a consequence of job design and social environment, whereas the individual-disposition hypothesis suggests that individual characteristics predispose some individuals to bullying or to being a bully (Nielsen and Einarsen, 2018).Supporting the first hypothesis, workplace bullying often occurs in combination with other psychosocial risk factors at work, such as role conflicts, problems in the organization of work and work tasks, high job demands, excess workload, insecurity, poor team climate, dissatisfaction with leadership and towards organization.There is also evidence that supports the latter hypothesis by linking personality traits, such as high neuroticism and low extraversion, to increased likelihood of being a target of bullying (Nielsen and Einarsen, 2018).Thus, both work environment and individual characteristics have been associated with risk of bullying.However, less is known about factors that predict change in bullying.This information would be particularly useful for organizations in efforts to prevent bullying at workplaces.
There are standard tools for assessing risk of health outcomes, including risk scores for cardiovascular disease (Alaa et al., 2019;Damen et al., 2016;Sahle et al., 2017), diabetes (Abbasi et al., 2016), various cancers (Gray et al., 2016;McCarthy et al., 2020;Stark et al., 2019), and those for long-term sickness absence and disability pension (Airaksinen et al., 2017(Airaksinen et al., , 2018)).However, a risk score for workplace bullying is not available.Bullying is more difficult to measure and predict than the onset of an illness due to its subjective nature (people may perceive bullying differently), and because risk factors are likely to be related to work context, i.e., work characteristics and work unit psychosocial factors, in addition to individual-level factors.A workplace survey measuring a wide range of both individual and work contextual factors may provide a basis for development of such prediction model.
Harrell's C-index (also known as the concordance index) is a widely used measure of predictive performance (Harrell, 2022).It provides a goodness of fit value for risk scores.Values above 0.7 indicate that the score is good at determining which of two randomly selected participants will have the outcome whereas value 0.5 means the risk score prediction is no better than a coin flip.Other useful indices in prediction are a) the detection rate, which denotes the proportion of test positive individuals among people who experience the outcome at follow-up, and b) false positive rate, that is, the proportion of individuals with a positive test result among those who did not experience the outcome.A further useful indicator is the ratio of true to false positives, also referred to as the odds of being affected given a positive result.This can be calculated by using the detection and false positive rate as well as information on frequency of the outcome over a specified time.
Using work unit and individual-level data on all sociodemographic and work-related variables of the Finnish Public Sector Study, the aim of this data-driven study was to determine the extent to which changes in workplace bullying can be predicted by employee responses to standard workplace surveys.

Study design and participants
The Finnish Public Sector (FPS) study is a prospective cohort study of employees in 11 towns, three hospital districts, and two other health care organizations in Finland (Airaksinen et al., 2019;Ervasti et al., 2022).The study has received ethical approval from the Ethical Committee of the Helsinki and Uusimaa hospital district (HUS/1210/2016).Data from two survey waves were used in the present analyses: 2017-2018 (T1) and 2019-2020 (T2).
We developed and examined four prediction models.The selection of the samples for each model is described in Fig. 1.In the Model 1, the aim was to identify individual-level factors predicting the start of bullying.Employees who were bullied at work at T1 were excluded.The outcome variable was reporting being bullied at T2.The analytical sample was 43,635 participants.
In Model 2, the aim was to identify factors predicting the start of bullying with predictor variables aggregated to work unit level and assigned to T1 survey respondents.If the work unit had less than 5 respondents, the work unit aggregate was calculated from the next level unit up in organization hierarchy.We retained the exclusion criteria from the first model and the final analytical sample was similar to first model (n = 43,635 participants from 5777 work units).
In Models 3 and 4 (presented in supplementary material), the aim was to identify factors predicting the end of bullying among employees who reported bullying at T1 (n = 4902).The predictor variables were from individual-level data in Model 3 and from work unit level in Model 4.

Measures
Workplace bullying was measured with the following question: "Psychological violence or bullying at work refers to the constant, repeated isolation of a member of the working community, belittling one's work effort, threats, talking behind one's back or other forms of pressure.Have you been the target of such bullying in the past 12 months?(yes/no).
In addition, participants were asked a total of 87 questions on their sociodemographic characteristics, work characteristics, psychosocial work environment factors, leadership, turnover, and retirement intentions.Detailed description is provided in Web appendix 1.
Sociodemographic characteristics included respondents' sex, age, type of job contract, occupational title (ISCO-coded), job tenure, working time (full or part-time; day work or shift work; years in shift work).
For work characteristics, the survey included 5 items on job demands and 9 items on job control (Karasek, 1979;Karasek and Theorell, 1990).Effort at work was measured with one question and rewards were assessed with three questions (Kivimäki et al., 2007).Worktime control was measured using 7 items (Ala-Mursula et al., 2004;Vahtera et al., 2010).Job insecurities were measured with 5 items, and changes at work with 2 items.
Procedural justice was measured with 7 items, and relational justice with 6 items (Moorman, 1991).Supervisor support was measured with 4 items, and support from the work unit to supervisor also with 4 items.Performance appraisals/career development discussions were measured with 2 items: having had such a discussion within the last 12 months, and whether the discussion was perceived useful.
Team climate was measured with 14 items on four dimensions: participation safety, support for innovation, vision, and task orientation (Anderson and West, 1996;Kivimaki and Elovainio, 1999).Discrimination at work was measured with a single item: Is there discrimination due to age, gender, education, opinion, status, origins, language, religion, believes/convictions, political activity, trade union activity, health, disability, sexual orientation, or gender identity/gender expression?
Job satisfaction was measured with 5 items, retirement intentions and turnover intentions with a single item each.
The response format in most of the items in the original survey was a five-point scale from 1 = strongly agree to 5 = strongly disagree.For ease of interpretability, we reverse-coded the items so that a high score always indicated a stronger agreement.
For prediction models 1 and 3 we used responses to these questions as individual-level variables.If individual-level information was missing, it was replaced with the corresponding work unit level aggregate (totaling 0.7% of imputed observations).In sensitivity analysis including health variables, missing data was imputed using predictive mean matching.
For prediction models 2 and 4, questionnaire responses and organizational records on employee turnover and work unit size were aggregated to work unit level resulting in 89 predictor variables.

Statistical analysis
To develop prediction models, we split the data with a 75/25 split, to training and test datasets stratifying for the outcome so that the proportion of those who had been bullied remained the same in both datasets.We then used bootstrap enhanced least-absolute-shrinkageand-selection-operator (lasso) with logistic regression and crossvalidation to find the optimal lambda value to select our predictors (Bach, 2008).
As a sensitivity analysis for the best model, instead of selecting the optimal lambda from the cross validation, we chose a lambda value that was one standard error from the optimal value.This commonly used approach allowed us to get a more parsimonious prediction model while still retaining reasonable accuracy compared to the model using the optimal value (Chen and Yang, 2021;Hastie et al., 2009).
Health and lifestyle are sensitive issues not often included in personnel surveys as employees may hesitate to give this information to the employer.Thus, as the aim was to develop a tool for prediction of bullying that could be used by employers, we did not include health and lifestyle -related survey questions in the main analysis.However, we performed a sensitivity analysis to determine whether including marital status, psychological distress (Goldberg et al., 1997), sleep problems (Juhola et al., 2021), perceived health (Robine et al., 2003), perceived work ability (Ahlstrom et al., 2010;Ilmarinen et al., 1997;Tuomi et al., 1997), body mass index (Halonen et al., 2014), alcohol use (Ervasti et al., 2018), smoking (Heikkilä et al., 2012), and physical activity (Fransson et al., 2012) would improve the predictive performance of the best model.The measures are described in detail in Web appendix 1.
All predictors were standardized for lasso.The predictors selected by regular lasso may vary depending on the sample and how strongly the predictors are correlated.With bootstrap enhanced lasso the final predictors are selected as those predictors that are present in a set proportion of the bootstrap replications.For the current study, we used 100 bootstrap replications and set the threshold for predictor selection to 95%.
We then fit a model using the predictors retained from the bootstrap enhanced lasso model to the test dataset.Using that model we compared the predictions against the observed bullying cases, plotted a ROC curve, and computed the area under curve (AUC).We also plotted the predicted probabilities for bullying to start (and bullying to stop) for ease of risk comparison between those who were bullied and those who were not.

Results
Of the 102,399 employees in the eligible population at T1, 72,968 (71%) participated in the T1 survey.Of them, 13,373 (18%) had left the organization and were no longer eligible at T2.Leaving the organization was more common among those who reported being bullied at T1 (OR = 1.33, 95% CI 1.25-1.40).After excluding those not eligible and those with missing data on work unit aggregates at T1, 59,346 were eligible at T2.The response rate at T2 was 83% resulting in 49,456 participants.After excluding those with missing data on bullying at T1 or T2, the final analytic samples included 43,635 employees who were not bullied at T1 and 4902 employees who were bullied at T1 (Fig. 1).
The descriptive characteristics of the study population by prediction model are shown in Table 1.Participants included in models 1 and 2 predicting the start of bullying were more often with higher socioeconomic status than participants included in models 3 and 4 predicting the end of bullying while any other differences were small.

Model 1: Individual-level factors predicting the start of bullying
Statistical significance of all unadjusted bivariate associations between individual-level factors and start of bullying is shown in Fig. 2. Most of the statistically significant associations were related to work unit/team climate.Some variables related to leadership and management also reached statistically significance.
Using a 95% threshold in the bootstrap enhanced lasso models, we obtained a model that included 25 of the 87 original predictors.The individual-level predictors for future bullying, and their odds ratios with 95% confidence intervals are shown in Table 2.The probability of bullying to start was higher with increasing discrimination (OR = 1.68, 95% CI 1.49-1.89),but also with other factors, including increasing number of women at work unit, when employees felt they had to work unreasonable amounts of work, work required well-developed skills, there was a threat that some work tasks would be terminated, when rewards were gained through personal satisfaction, with increasing amount of investment of abilities and resources on work, work included lots of repeated tasks, changes at workplace were evaluated as mostly negative, and with increasing time in shift work (range of ORs between 1.01 and 1.21).
The probability of bullying to start was lower with managers, senior officials and professionals, when the employees felt they were understood and accepted, the supervisor could be trusted, rewards were gained through recognition and respect, supervisor's personal biases did not influence decisions, work was perceived as challenging enough, members of the work unit considered suggestions for improvement made by others members of the work unit and understand the goals of the work unit, employees were satisfied with the level of challenges at work, and with longer time with the current employer and in the same position (range of ORs 0.80 to 0.99, Table 2.) For four variables ('rewards through personal satisfaction', 'possibilities to take breaks during working day', 'goals of the work unit are achievable', and 'supervisor considers subordinates' opinions in important matters'), the direction of bivariate association changed to the opposite in Lasso model incorporating multiple predictors (Table 2).Thus, to get a more parsimonious prediction model, we performed a new iteration of Lasso analysis increasing the lambda value as one standard error from the optimal value.This resulted in a model with the following six predictors associated with start of workplace bullying: unreasonably high workload, threat that some work tasks will be terminated, discrimination at work unit, and reversely with belonging to a work unit where everyone feels they are understood and accepted, having a supervisor who can be trusted, and a longer time in current position (presented in bold in Table 2).

Model 2: Work unit -level factors predicting the start of bullying
The unadjusted bivariate associations between work unit level factors and start of bullying are shown in Fig. 3.Most of the statistically significant associations were related to work unit/team climate, work characteristics and leadership.Of the work unit sociodemographic characteristics, only higher proportion of managers, senior officials and professionals was statistically significantly associated with start of bullying.
The bootstrapped enhanced lasso model found that the start of workplace bullying was best predicted using 11 variables (Table 3).Of all the 11 predictors, only the item referring to critical appraisal of weaknesses was not statistically significant when the model was fitted using the test portion of the data.The probability of start of bullying was higher if the work unit was characterized by items, such as lots of repeated tasks, discrimination, high proportion of temporary employees, high turnover intentions, and employees reporting longer duration of shift work.The probability of start of bullying was lower if the work unit was characterized by consistent decision making, hard work, implementing performance appraisals/career development discussions, long work careers with the same employer, and high proportion of managers, senior officials, and professionals.A further analysis with lambda set at one standard error from the optimal value resulted in a model with a single predictor, discrimination at work unit.

Model performance: start of bullying
The ROC-curves and density plots for individual and work unit level models are shown in Fig. 4. The curves of those being bullied and those not being bullied were largely overlapping and the predictive performance was modest for both models, with AUC scores being 0.68 and 0.62, respectively.
To obtain a more parsimonious model, we performed another sensitivity analysis with lambda value set as 1 standard error from the optimal value.The resulting C-statistics were 0.66 (individual) and 0.60 (work unit).The corresponding ROC-curves and density plots are shown in Web Fig. 2).
Table 4 shows detection rate, false positive rate, and ratio of true to false positive for the two models using various risk thresholds for test positive result.For a 5% cut-off for positive test results (i.e., being bullied), the detection rate was 79-84%, false positive rate was 56-71%, and ratio of true to false positive test result was 1 to 11-12.That is, while many cases were detected, this came at a price of many false positives, i. e., cases where bullying does not start despite our prediction.With 15% cut-off as the threshold, detection rate decreased to 3-17%, false positive rate to 1-5%, and ratio of true to false positives was 1-5.Increasing the cut-off point from 5% to 15% thus resulted in better accuracy (less false positives), but also to a very poor detection rate.

Models 3 and 4: Individual and work unit -level factors predicting the end of bullying
In analysis of ending bullying among those exposed to bullying, most of the statistically significant associations were related to work unit/ team climate (Web Figs. 3 and 4).Using the bootstrapped enhanced

Table 2
Individual-level predictors of start of workplace bullying in bivariate analysis and Lasso-regression.The predictors that remained in the model where lambda value was set as one standard error from the optimal value, are in bold.* The latent construct/sum scale is indicated when the individual item is part of a scale.lasso, seven of the original 89 variables were retained in the final model: effort to share information within the work unit, opportunities to develop one's special skills, and time with the current employer were associated with higher probability of bullying to stop (range of ORs 1.01 to 1.25) and discrimination at the workplace, work required welldeveloped skills, threat of lay-offs, and changes at work (range of ORs 0.58 to 0.95) with lower probability of stopping bullying (Web Table 2.) The AUC of 0.62 indicated weak predictive performance (ROC-curve and density plot available in Web Fig. 5 Panels A and B).
With work unit level variables, only three factors were predictive of the end of bullying (Web Table 3).The only statistically significant predictor in Lasso was opportunities to influence and take breaks during the workday, and the direction of the estimate was reversed such that the probability for end of workplace bullying was lower for employees with influence on the breaks during working day.The performance of the model did not differ from flipping a coin (AUC 0.51, the ROC-curve and density plot in Web Fig. 5, Panels C and D, and detection rate, false positive rate, and ratio of true to false positive in Web Table 4).

Discussion
With 43,635 Finnish public sector survey respondents at two time points, our aim was to build a prediction model for workplace bullying.The model with the best predictive ability for individualized assessment of bullying risk included 25 items related to sociodemographic characteristics of the work unit, type of work, work characteristics, and the characteristics of the team climate and the supervisor.The probability that a participant who experienced a start of bullying at workplace had a higher risk score than a participant who did not get bullied (the C-index) was 68%.We aimed for a more parsimonious model by using stricter lambda value in lasso modelling.That model included six predictors describing discrimination, job demands, job insecurities, team climate, relational justice, and job tenure, and performed nearly as well with a Cstatistic being 0.66.As sensitivity analyses, we included health and lifestyle variables; analyzed work unit level predictors; and used end of

Table 3
Work unit-level predictors of start of workplace bullying in bivariate analysis and Lasso-regression.The predictors that remained in the model where lambda value was set as one standard error from the optimal value, are in bold.bullying (instead of start of bullying) as the outcome variable, but the predictive performance did not improve or was even more modest in these models.While many work-related characteristics had a bivariate association with the start of bullying, the prediction model based on these variables did not reach a satisfactory predictive power (C-index >0.7) for individualized risk assessment.This may be due to the items being associated with a shared (latent) phenomenon, such as occupationspecific work characteristics (Kristensen, 1995).As bullying is more common within the health care sector (Zapf et al., 2003), the prediction model consists of characteristics that describe work in the health care sector, such as high proportion of women, repeated tasks, and shift work, whereas work characteristics linked with white-collar work, such as working very hard is seen as a protective factor.
Dichotomized prediction scores to test positive and test negative are useful in decision-making and can be used to inform decisions regarding which work units would benefit most from an intervention aimed at prevention of bullying at workplace.If the intervention has minimal negative side-effects and is inexpensive, it may be feasible to emphasize high detection rate and tolerate a high false positive rate.In contrast, if the intervention is expensive and new with uncertain profile of adverse side effects, the employer might prefer low false positive rates.In the model with the best predictive power (Model 1), a 5% rate of false positives detected only one in five employees at risk of getting bullied in the future and may therefore have limited practical value.
Employers in Finland are mandated to take steps to eliminate workplace bullying and ensure safety at workplaces.Thus, despite the relatively poor prediction for individualized risk assessment, this study provides useful information for employers tackling workplace bullying

Table 4
Detection rate (%), false positive rate (%), and ratio of true to false positive for predictive probability using various risk thresholds.Job crafting means self-initiated changes that employees make to work characteristics to shape the working conditions (Tims et al., 2012).Job crafting can also be a collaborative and shared group activity (Tims et al., 2013), and could thus be used as a method of bullying prevention.Employees may together determine that there is a need to increase the feeling of being understood and accepted, and then determine the ways in which to alter and modify their work to attain this shared goal.Although collaborative job crafting may not be suitable or sufficient for all types of working conditions associated with increased risk of bullying (e.g., unreasonably high workload), it could be a useful tool for others, such as activities to seek and increase social support and togetherness (Harju et al., 2021;Seppälä et al., 2018;Tims et al., 2013).
Our findings confirmed earlier cross-sectional results (Dussault and Frenette, 2015;Nielsen, 2013) suggesting that favorable leadership behavior is associated with lower probability of bullying.Moreover, our finding on the association between threat of work task termination and increased probability of bullying supports the work environment hypothesis proposing that role stressors are important antecedents of bullying (Reknes et al., 2014).
Longitudinal studies have focused on the mental health outcomes of bullying and found that bullying is associated with increased risk of psychological distress (symptoms of depression and anxiety), but whether the association is reciprocal needs further investigation (Boudrias et al., 2021).Indeed, our findings suggest that psychological distress may precede bullying victimisation, although adding psychological distress and other health-related predictors did not improve the predictive performance of our risk models.
It seems that the end of bullying is even harder to predict than the start of bullying.This emphasizes preventive measures and zerotolerance rather than measures or interventions aimed at stopping bullying after it has begun.For example, good leadership may play a role in preventing bullying, but not in bullying to stop once started.
The strengths of this study include a large sample (>40,000) of participants who filled up a multifaceted questionnaire with 87 items, and a longitudinal (two-wave) study design which allowed us to examine the change (that is, the start and the end of) bullying at workplace.We performed our analysis with both individual-and work unit level predictor variables to be able to determine the model with the best predictive power.Despite the strengths, we were not able to achieve a model with satisfactory predictive power.
The study has several limitations.The data were from the public sector employees in Finland, so the generalizability should be investigated in other settings.While we had a large dataset, there were relatively few individuals among whom bullying had ended within the twoyear interval.Thus, the prediction of bullying to stop resulted models with higher uncertainty that those for prediction of initiation of bullying.It is also noteworthy that being bullied at T1 increased the odds of being non-eligible at T2, suggesting that those who were bullied left the organization before the next survey.The association between workplace bullying and subsequent turnover has been demonstrated in also earlier studies (Boudrias et al., 2021).Thus, our estimates may be underestimates of the true effect size.Moreover, further research is needed to assess the generalizability of our findings because predictive algorithms may vary by occupation, industry and settings.
In conclusion, while several predictors were identified at group level, reliable individualized detection of individuals at risk of becoming bullied at workplace was not successful.The predictive performance of the developed risk scores was suboptimal, and we do not recommend their use as an individual-level risk prediction tool.In contrast, they might be useful tool to inform decision-making when planning the contents of interventions to prevent bullying at an organizational level.

J
.Ervasti et al.

Fig. 2 .
Fig. 2. Bivariate associations between individual-level predictor items and workplace bullying to start.Colors indicate themes the individual items are grouped.Horizontal red line indicates Bonferroni corrected statistical significance (p = 0.05/87 = 0.0006).Predictors with the highest -log10 (p-values) are labeled.(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 3 .
Fig. 3. Bivariate associations between work unit level predictor items and workplace bullying to start.Colors indicate themes the individual items are grouped.Horizontal red line indicates Bonferroni corrected statistical significance (p = 0.05/89 = 0.0006).Predictors with the highest -log10 (p-values) are labeled.(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 4 .
Fig. 4. Model performance.Panel A: ROC-curve for work unit -level prediction of bullying to start.Panel B: Density plot for being bullied at T2 by predicted probability (logistic regression, work unit -level predictor variables).Panel C: ROC-curve for individual-level prediction of bullying to start.Panel D: Density plot for being bullied at T2 by predicted probability (logistic regression, individual-level predictor variables).

Table 1
Descriptive sociodemographic characteristics at T1 of the study participants by prediction model.

that some work tasks will be terminated 1.21 1.17, 1.25 1.12 1.08, 1.16
at organizational level.There were six characteristics of the individuals, work, and psychosocial work environment that were robustly associated with increased risk of bullying: discrimination at work unit, unreasonably high workload, threat that some work tasks will be terminated, lack of feeling of being understood and accepted, untrustful supervisors, and a short time in current position.This information characterizes specific factors that could be targeted with workplace interventions.