Load distribution in ball screws is a representation of the ball contact stress, and it is fundamental to understanding the behavior of these machine elements. This work aims to conduct a multi-variable analysis of the load distribution in ball screws. For this purpose, a numerical tool is developed for the generation and calculation of ball screw finite element (FE) models, which has been validated against the state of the art. This tool is based on the combination of an analytical contact model and the use of high-order FE models for the analysis of the load distribution of ball screws and stands out for its accuracy (less than 1% error against high-order FE models), adaptability, versatility (models are generated with more than 20 design variables and they can be introduced as components in larger models) and efficiency (being the computational time 1.25% of that of a high-order FE models) with respect to other existing models. Many different design variables (number of start threads, pitch, contact angle, ball size, slenderness and load arrangement) are studied in order to obtain a general characterization of the morphology of the load distribution in ball screws. Among them, the most influential variables on the load distribution and therefore on the structural behavior of ball screws are, load arrangement (with ratio r variations of up to 25% on the same ball screw) and slenderness (with ratio variations of up to 13% on ball screws with two turns of difference). The two most characteristic features, the non-uniformity at a local and global level are identified, along with as the possible causes of their appearance and the consequences that they may cause.