Aim. It is essential to understand the extent to which job characteristics impact work-related musculoskeletal disorders (WMSDs), and to calculate the probability that an employee will suffer from a musculoskeletal disorder given their working conditions. The objective of this research is to identify the relationships between WMSDs and working characteristics, by developing a Bayesian network (BN) model to calculate the probability that an employee suffers from a musculoskeletal disorder. Methods. A conceptual model was constructed based on a BN. This was then statistically tested and corrected to establish a BN model. Results. Experiments verified that the BN model achieves a better diagnostic performance than artificial neural network, support vector machine and decision tree approaches, and is robust in diagnosing WMSDs given working characteristics. Conclusion. It was verified that working characteristics, such as working hours and pace, impact the incidence rate of WMSDs, and a BN model was developed to probabilistically diagnose WMSDs.