A Cardiovascular Risk Score for Use in Occupational Medicine

Cardiovascular disease is one of the most frequent causes of long-term sickness absence from work. The study aims to develop and validate a score to assess the 10-year risk of unsuitability for work accounting for the cardiovascular risk. The score can be considered as a prevention tool that would improve the cardiovascular risk assessment during health surveillance visits under the assumption that a high cardiovascular risk might also translate into high risk of unsuitability for work. A total of 11,079 Italian workers were examined, as part of their scheduled occupational health surveillance. Cox proportional hazards regression models were employed to derive risk equations for assessing the 10-year risk of a diagnosis of unsuitability for work. Two scores were developed: the CROMA score (Cardiovascular Risk in Occupational Medicine) included age, sex, smoking status, blood pressure (systolic and diastolic), body mass index, height, diagnosis of hypertension, diabetes, ischemic heart disease, mental disorders and prescription of antidiabetic and antihypertensive medications. The CROMB score was the same as CROMA score except for the inclusion of only variables statistically significant at the 0.05 level. For both scores, the expected risk of unsuitability for work was higher for workers in the highest risk class, as compared with the lowest. Moreover results showed a positive association between most of cardiovascular risk factors and the risk of unsuitability for work. The CROMA score demonstrated better calibration than the CROMB score (11.624 (p-value: 0.235)). Moreover, the CROMA score, in comparison with existing CVD risk scores, showed the best goodness of fit and discrimination.


Introduction
Cardiovascular disease (CVD) is the most common cause of death globally, with 17.9 million deaths estimated each year. Of these deaths, 85% are due to cerebrovascular diseases, and one third of these deaths occur in people under 70 years of age [1,2]. CVD is a major cause of disability for all age-groups worldwide, including working-age individuals [3]. Moreover, several occupational risk factors including work-related stress, long working hours or manual handling of heavy loads are associated with increased risk of CVD [4,5]. The hypothesized mechanism between long working hours and cardiovascular health is related to reduced time available for other activities besides work. Employees working long hours spend more time at the workplace and thus they might be increasingly exposed to psycho-social, chemical and/or physical occupational risk factors which could adversely impact the cardiovascular system [6] or, on the other hand, workers engaged in performing sedentary jobs (i.e., office work) have reduced physical activity and increased CVD risk [7]. Obviously, these associations might be stronger for specific job categories, which involve exposure to occupational risk factors whose damaging action can mainly affect the cardiovascular system but, to certain extent, they apply for all types of jobs [4,5]. Furthermore, CVD may negatively impact productivity and quality of work [8]. CVD is one of the most frequent causes of long-term sickness absence from work. Additionally, the illness consequences may cause a radical change of work activity [9], which prompts the question of whether a thorough cardiovascular risk assessment and appropriate preventive strategies are appropriately integrated as part of the occupational health surveillance. The majority of existing cardiovascular risk scores in the literature have been adopted in medical settings different from occupational medicine [10][11][12][13][14][15][16]. These scores have been created to evaluate the risk of diagnosis cardiovascular diseases or risk of mortality due to CVD either in a general population or in a certain occupational cohort (i.e., male industrial workers or government officials) [7,17,18]. Interestingly, previous research recalibrated the Framingham risk score to evaluated the risk of unsuitability for work in a cohort of Italian workers without, however, obtaining clinically meaningful results [3]. Therefore, in this context, this study aims is to derive, calibrate, and validate a new score to predict the risk of unsuitability for work, accounting for the cardiovascular risk factors. Specifically, the study aims to provide a prevention tool to improve the cardiovascular risk assessment during health surveillance visits under the assumption that a high cardiovascular risk might also translate into a high risk of unsuitability for work. The score, named Cardiovascular Risk in Occupational Medicine (CROM), was also compared with CVD risk scores already present in the literature.

Study Design and Study Population
Th study is a retrospective cohort study using data from a cohort of 11,079 workers. Workers were employed by the Naples municipalities and were examined as part of their scheduled occupational health surveillance, at the occupational health outpatient clinic of the Department of Public Health of the University "Federico II" of Naples between January 2006 and December 2016. All clinical assessments were part of clinical practice in a university setting and data were fully anonymized. All subjects signed the general informed consent form, authorizing the use of observational clinical data for research purposes.

Study Variables
Data collection was performed during a medical examination for scheduled occupational health surveillance. According to the Italian occupational medicine legislation and above all considering the different levels and types of occupational and/or job strain, workers were classified into four groups, which translated into a different periodicity of their scheduled occupational health surveillance visit (every one/two/three/five years) [19]. A detailed list of each job included in this analysis can be found in the Appendix A. At the end of each visit for health surveillance, according to Italian occupational medicine guidelines, a fitness for work judgment was issued by the occupational medicine physician and later confirmed by a senior occupational medicine consultant. The fitness for work judgment has three possible outcomes: (i) suitability; (ii) suitability with limitations or prescriptions (with a consequent reduction or remodeling of the job strain for the worker); (iii) total unsuitability, with a radical change of activities within the job. Incident diagnosis of unsuitability for work (both partial with limitations or prescriptions and total) over a 10-year follow-up period was considered as study outcome. As part of the medical examination, socio-demographic and clinical history data, as well as laboratory and instrumental data, were collected. Study covariates included age, sex, smoking status, blood pressure (systolic and diastolic), body mass index (BMI), height, diagnosis of hypertension, diagnosis of diabetes, Ischemic heart disease, mental disorders (bipolar disorder and moderate/severe depression) and prescription of drugs, especially anti-diabetic medications (biguanides, sulfonylureas, insulin, others), anti-hypertensive (ACE inhibitor (ACEi) or angiotensin receptor blocker (ARB); others). Additionally, taking into account the considerable heterogeneity of the work carried out by the study participants, in order to consider the different exposures (and the relative extent) to occupational risk factors and job strain impacting on the cardiovascular system, we have introduced a new variable "job risk class" (highest risk class, medium-high risk class, medium-low risk class and lowest risk class) based on the health surveillance protocol adopted which, in turn, is established on the basis of the type and level of exposure to occupational risk factors [3]. On the contrary, a categorical variable describing each job was not included due to sample size constraints. Analysis of variables was conducted at baseline, considering only the first visit for each individual. For each covariate in the study, baseline values were considered. To reduce the amount of missing data at baseline for study covariates, we used the first available clinical recording within five years for each individual for individuals [20]. Cigarette smoking status was ascertained by self-reporting. The individuals with missing data on smoking were classified as non-smokers. Similarly, the absence of recording for CVD was considered as the absence of the condition [3,20,21].

Statistical Analysis
Missing data were present for the following variables: body mass index (3.25% missing), systolic blood pressure (2.67% missing), diastolic blood pressure (2.66% missing), height (3.25% missing). We used multiple imputations by chained equations to replace missing values for body mass index, systolic blood pressure, diastolic blood pressure and height [22,23]. Predictor variables to include in the imputation models were identified employing a missing-pattern analysis that used a multivariable logistic regression model to assess factors associated with missing data. We included sex and diagnosis of hypertension variables in the imputation model as likely to be associated with the recording of risk factor. We carried out five imputations, as this has relatively high efficiency and was a pragmatic approach accounting for the size of the data set and capacity of the available servers and software [13]. We used the Rubin rule to combine the results across the imputed data set. Dir references in covariate distribution between occupational risk groups were explored employing ANOVA and χ-squared, as appropriate. Variables such as sex and age were considered, but they were not statistically significant [10]. We implemented Cox proportional hazards regressions to predict the risk of a diagnosis of unsuitability over a 10-year follow-up period (Appendix C). To build the multivariate Cox proportional hazard regression model, for the model derivation a model-building strategy was employed. Similarly to previous research [20], the approach was based on evaluating both the statistical significance and clinical relevance of the variables to ensure those candidate variables were likely to be clinically important and to reduce the over-fitting and optimism of the model. First, a multivariate Cox proportional hazard regression model including all candidate independent variables was run. In line with previous research [3,19,23], variables were included in the multivariate model if they had a hazard ratio of less than 0.90 or more than 1.10 (for binary variables) and were statistically significant at the 0.01 level in the univariate analyses. Although not significant, age and sex were considered because they were likely to be associated with the outcome [10]. Interactions between predictor variables and age at baseline were also examined and then significant interactions were included in the final models. All the continuous variables were naturally logarithmically transformed to improve discrimination and calibration of the models and to minimize the influence of extreme observations [10]. In addition, a multivariate Cox proportional hazard regression model using data where cases with missing data are excluded was run (Appendix B). Model discrimination was evaluated through an interval validation employing a k Fold Cross-Validation procedure and the C-statistic's measure [24]. Harrell's C statistic is a goodness of fit measure for models which produce risk scores and values approaching 1 indicate a good model.
The calibration of the validated risk score, a measure of agreement between observed and predicted events, was evaluated using the Gronnesby and Borgan Test based on martingale residuals [25]. The output returns a chi-square value (chi-squared) and a p-value (e.g., Pr > ChiSq). Small p-values mean that the model is a poor fit.
Using the person's absolute risk [26], empirical research on the risk threshold was undertaken, as previously suggested [27]. Finally, the validated score was compared with the widely used CVD risk scores. Three CVD risk scores already present in the literature (the Framingham score, European SCORE and Q-risk score) were recalibrated employing Cox proportional hazard regression models. Specifically, one model included risk factors present in the Framingham score (age, sex, smoking status, systolic blood pressure and antihypertensive), a second model included the risk factors present in the European SCORE (age, sex, smoking status, systolic blood pressure, body mass index (BMI)) and a third model was developed using risk factors present in the Q-risk score (age, sex, smoking status, systolic blood pressure, body mass index (BMI), diagnosis of diabetes, diagnosis of Ischemic heart disease, anti-hypertensive, diagnosis of mental disorders) [12,28].
Since the cohort used in the study is a healthy working age cohort, the coefficients of unavailable covariates were kept constant and equal to 1 under the assumption that no recording meant absence of the condition.

Results
Between January 2006 and December 2016, 11,079 workers were examined for health surveillance by trained physicians at the Occupational Medicine Outpatient Clinic of "Federico II" University Hospital. A total of 57.95% were men, the mean (SE) age was 52.35 (8.46), 6.31% were in the highest risk class, 53.06% were in the medium-high risk Class, 32.85% were in the medium-low risk Class and 7.79% were in the lowest risk class. The statistical description of the variables after multiple imputations is shown in Table 1. Tables 2 and 3 show the adjusted hazard ratios for the CROM A score and the CROM B score. The CROM A score included all the variables above; the CROM B score was the same as the CROM A score except that it included only variables that were statistically significant at the 0.05 level in the Multivariable Cox Model.
The number of events of unsuitability for work was sufficiently high to allow the risk score derivation. In fact, 852 events (i.e., 852 diagnoses of unsuitability for work) were recorded.
The results show a positive association between most of the cardiovascular risk factors and the risk of unsuitability for work. Moreover, the tables (Tables 2 and 3) show a positive proportional link between the worker risk classes and the risk of unsuitability for work. Moreover, the estimated regression coefficient increases with the level of exposure to occupational risk factors and so the risk of unsuitability for work gradually decreased with decreasing level of exposure to occupational risk factors.
It possible to observe this result also in Figure 1, where the hazard function value increases with both follow-up time and with the level of exposure to occupational risk factors (highest risk class, medium-high risk class, medium-low risk class and lowest risk class).   In addition to measures of discrimination of the model, for both scores CROM A and CROM B , we calculated Harrell's C statistics. The result for the CROM A score was 0.700 (95% CI: 0.698-0.702), and was 0.699 (95% CI: 0.696, 0.702) for the CROM B score.
For the calibration assessment, the CROM A score chi-square was 11.624 (p-value: 0.235) whilst was the chi-square for the CROM B score was 11.000 (p-value: 0.275); therefore, the CROM A score fitted better than the CROM B score to the set of observations. Comparing the two scores developed in the study, the CROM A score had better calibration and discrimination. Figure 2 shows the person's absolute risk of unsuitability for work [10,26] for the study cohort using the CROM A score. For over 60% of cases, the status of unsuitability for work occurred for a risk greater than or equal to 46.2%; therefore, the risk of 40% appears to be a good threshold for discrimination of high-risk workers.
We calculated Harrell's C statistics for widely used CVD risk scores. The result for the Framingham score was 0.581 (95% CI: 0.577, 0.585), for the European score it was 0.590 (95% CI: 0.588, 0.593) and for the Q-risk score it was 0.608 (95%CI: 0.609, 0.612). For the calibration assessment, the Framingham score's chi-square was 13.512 (p-value: 0.140), the European score's chi-square was 14.126 (p-value: 0.117) and the Q-risk score's chi-square was 18.592 (p-value: 0.028). The CROM A score, in comparison with the used CVD risk scores above, had the best goodness of fit and discrimination. Absolute risk: The histogram shows the person's absolute risk calculated for the study cohort using the CROM A score. About 67% of workers present a cumulative risk larger than or equal to 0.4; therefore, the risk of 40% appears to be a good threshold for the discrimination of high-risk workers.

Discussion
Using data from a cohort of 11,079 workers with different jobs, this study derived and validated a new risk score to evaluate the 10-year risk of a diagnosis of unsuitability for work due to cardiovascular diseases. The equation incorporates canonical cardiovascular risk factors as well as new risk factors associated with the workplace. Our findings confirm the linkage between the CVD risk factor and the workplace [2,7,17]. The proportion of individuals with a diagnosis of unsuitability for work grows progressively from the lowest to the highest work risk classes. Many CVD risk scores have been derived and are currently used in clinical practice. However, those scores have been designed to be used on the general population and might show little accuracy when applied to specific populations like the working-age population.
In our opinion, the use of CROM score could have significant practical implications in the field of occupational medicine, especially in implementing the preventive potential of health surveillance. Indeed, occupational medicine is basically a preventive specialty and occupational physicians have a fundamental role in protecting and improving the health of employees in relation to their work and to ensure a continual improvement of working environment and preventive and/or protective measures [29,30]. Therefore, in this context, a high CROM score (predictive of a likely unsuitability of the worker in the medium-long term) should not be used to possibly anticipate such unsuitability (this would not be correct from an occupational medicine point of view and not even acceptable from an ethical perspective as, when the score is calculated, it is likely that the worker is completely fit to perform his work), but rather it should represent an alarm bell for the occupational physician but, more generally, the entire occupational safety and health management system should therefore undertake to apply specific and appropriate additional preventive measures to protect the worker with a high CROM score's health. Obviously, the preventive strategy that is decided upon can be put into practice by exploiting a wide range of tools whose use should be integrated and in any case designated on the basis of occupational risk factors (which in turn depend on the type of work performed) to which the worker is exposed to. In other words, a high CROM value, let's say for example, to a worker who carries out the activity of warehouseman and who, in order to carry out his job tasks, is exposed to manual handling of heavy loads, night work and working activities that involve energy cost superior to six metabolic equivalents of the task, could be suggested the opportunity to reshape workloads or implement the health surveillance protocol (carrying out additional blood tests). However, on the other hand, the same CROM value, determined in a completely different type of worker such as an office worker, exposed only to a video display unit, would have quite different practical implications. Indeed, in this case, considering the low impact of exposure to occupational risk factors on the cardiovascular system, the main task of the occupational physician should be to carry out an adequate counseling action (i.e., health promotion) to reduce overall CVD risk factors (eventually including work-related ones such as psychosocial factors). However, differences in accuracy might also be explained by the fact that other scores were derived using not only different populations but also considering different outcomes. In addition, the score showed good discrimination power (0.700176 95% CI:0.700, 0.780), considering that similar levels from different scores have already been accepted when implementing prevention tools in a similar setting [7,31]. The majority of cardiovascular risk scores were developed to assess the CVD risk using the whole population and including accepted cardiovascular risk factors associated with [7,10,12,27]. Previous research has only occupational stress or job strain with which to assess the link between CVD risk factors and workplace [7,32]. Therefore, using only psychological job demand and occupational stressors, the job conditions were not examined in their entirety. To our knowledge, this is the first study that aimed to assess the risk of unsuitability for work adjusted for cardiovascular risk factors using a heterogeneous cohort of workers and a long follow-up time. In addition, this is also the first time that the job conditions were weighed by assessing both occupational risks (chemical, physical and biological risks) and job strain. While most studies consider mental health to be the leading cause of unsuitability for work, we included several mental disorders in a single covariate, minimizing its weight [33][34][35]. Several caveats merit discussion. Based on data availability, it was not possible to classify diagnoses of unsuitability for work according to whether it was CVD-related. Whilst this might be considered as a main limitation, many studies have shown that diagnosis of unsuitability for work is often multifactorial and involves many risk factors, including CVD risk factors which are likely to have contributed to it [36]. However, future studies should be conducted to externally validate the CROM score using datasets including information on reason for diagnosis of unsuitability for work.
Additional study limitations include the presence of missing data for clinical variables such as body mass index, systolic blood pressure, diastolic blood pressure, and height. However, we overcame the latter issue by using multiple imputations by chained equations.
A further important limitation of the study is the absence in the literature of a score that can be comparable to the CROM score due to different outcomes considered in each score and the model derivation carried out using different populations and settings. However, given that the CROM score should not be used to possibly anticipate unsuitability for work but as support to apply specific preventive measures, we overcame this issue by recalibrating three CVD risk scores already present in the literature employing Cox proportional hazard regression models to make them comparable with the CROM score (the Framingham score, European SCORE and Q-risk score).
Another limitation of the study includes the use of a limited kind of job activity; in fact, only those available to us have been used and so are present in the dataset. The existing literature supports the ineffectiveness of models derived and validated on a specific ethnic group when translated to a different one [16]. Therefore, external validation should be conducted to estimate the model predictiveness and validity on a more heterogeneous cohort including other ethnic groups as well and a broader range of jobs.

Public Health and Occupational Medicine Implications
CVD is one of the most prevalent causes of long sickness absence from work and may involve the radical change of individual work activities. Therefore, the score developed in this work could be used at each visit of health surveillance as a clinical tool to predict the 10-year risk of unsuitability for work using the personal information and medical history of the worker. Considering that high cardiovascular risk is also associated with a higher risk of unsuitability for work and that an occupational medicine doctor might have limited knowledge of the worker's medical history, it might be useful to have a tool with which to conduct an accurate assessment of the risk of unsuitability for work due to cardiovascular factors and to adopt non-pharmaceutical and pharmaceutical interventions to reduce it [37].
A modern occupational physician is a leading expert on mitigating the impact of health conditions on workers' professional activity. Therefore, in order to achieve this important goal, occupational physicians should not only consider possible workplacerelated threats to workers' health but they should also take into account any diseases, health issues, disabilities or risk factors that might be an obstacle to the adequate and safe performance of job activities [38]. In this regard, the management of workers with high CVD risk scores is a rather challenging issue which necessarily requires evaluating the complex interplay between CVD risk factors and exposure to occupational risk factors impacting the cardiovascular system. In the health surveillance context, the use of the CROM score could be useful, since it is an integrated and comprehensive tool. This would provide occupational physician with interesting data to ensure adequate safe and healthy working conditions.

Conclusions
We derived and validated an accurate score to predict the 10-year risk of unsuitability for work in occupational medicine. This score could be used as a preventive strategy clinical tool for cardiovascular risk assessment during the scheduled medical examination for health surveillance to reduce cardiovascular risk and so reduce the probability of unsuitability for work associated with cardiovascular risk factors.
Author Contributions: R.P. and M.T. conceived and designed the study. G.A. analyzed the data. G.A., P.A., F.B.-A., L.F., R.P. and M.T. discussed the data analyses and interpreted the results. G.A. and R.P. wrote the first draft of the manuscript. All authors critically revised the manuscript and agreed to act as guarantors of the work. G.A. has full access to all data used in the study and takes responsibility for the integrity of the data and the accuracy of the data analysis. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.

Institutional Review Board Statement:
All procedures performed in this study were in accordance with the ethical standards of the institutional and/or national research committee and with the 1964 Helsinki declaration and its later amendments or comparable ethical standards. Ethical approval was obtained by the University "Federico II" of Naples (112/18).
Informed Consent Statement: ll subjects signed the general informed consent form, authorizing the use of observational clinical data for research purposes.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
Job risk factors, according to Italian legislation on occupational medicine, are classified into four groups. The groups reflect the workers' risks (physical, chemical and biological risks) levels and they translate into a different periodicity of the scheduled visit to monitor health at work (every one/two/three/five years). Table A1 shows a detailed list of each job included in this analysis where the occupational risk group is composed of jobs with a similar surveillance protocol and therefore a similar health surveillance visit scheduling.

Appendix B
In addition to analysing that described above, a multivariate Cox proportional hazard regression model using data where cases with missing data are excluded was run. The results of the two analyses turned out to be comparable (Tables A2 and A3). Calibration assessment of the CROM score was performed using the Gronnesby and Borgan Test based on martingale residuals. The CROM A and CROM B models present a good calibration, which was confirmed also by Figure A1, showing a good fit to the crude hazard estimate for both models.

Materials and Methods
Cox proportional hazards regression to predict the risk of a diagnosis of unsuitability for work during a follow-up period of 10 years was implemented: where h(t) is the hazard function and it is defined as the instantaneous risk that a worker is evaluated as unsuitability for work, within a very narrow time frame. h 0 (t) is the baseline or underlying hazard function and corresponds to the probability of unsuitability for work when all the explanatory variables are zero.
In order to choose covariates to include in Cox's multi-variables regression model, for each risk factor, the Hazard Ratio of Cox's uni-variable model was evaluated.
Variables that had a hazard ratio of less than 0.90 or more than 1.10 (for binary variables) and were statistically significant at the 0.01 level were included in the model.
An assessment test for multicollinearity in both CROM models was run. Moreover, we tested for interactions between each variable and age. However, these interactions were not included in the final model because they were not significant.