Computational head injury models are promising tools for understanding and predicting traumatic brain injuries. However, most available head injury models are "average" models that employ a single set of head geometry (e.g., 50th-percentile U.S. male) without considering variability in these parameters across the human population. A significant variability of head shapes exists in U.S. Army soldiers, evident from the Anthropometric Survey of U.S. Army Personnel (ANSUR II). The objective of this study is to elucidate the effects of head shape on the predicted risk of traumatic brain injury from computational head injury models. Magnetic resonance imaging scans of 25 human subjects are collected. These images are registered to the standard MNI152 brain atlas, and the resulting transformation matrix components (called head shape parameters) are used to quantify head shapes of the subjects. A generative machine learning model is used to generate 25 additional head shape parameter datasets to augment our database. Head injury models are developed for these head shapes, and a rapid injurious head rotation event is simulated to obtain several brain injury predictor variables (BIPVs): Peak cumulative maximum principal strain (CMPS), average CMPS, and the volume fraction of brain exceeding an injurious CMPS threshold. A Gaussian process regression model is trained between head shape parameters and BIPVs, which is then used to study the relative sensitivity of the various BIPVs on individual head shape parameters. We distinguish head shape parameters into 2 types: Scaling components ${T_{xx}}$, ${T_{yy}}$, and ${T_{zz}}$ that capture the breadth, length, and height of the head, respectively, and shearing components (${T_{xy}},{T_{xz}},{T_{yx}},{T_{yz}},{T_{zx}}$, and ${T_{zy}}$) that capture the relative skewness of the head shape. An overall positive correlation is evident between scaling components and BIPVs. Notably, a very high, positive correlation is seen between the BIPVs and the head volume. As an example, a 57% increase in peak CMPS was noted between the smallest and the largest investigated head volume parameters. The variation in shearing components ${T_{xy}},{T_{xz}},{T_{yx}},{T_{yz}},{T_{zx}}$, and ${T_{zy}}$ on average does not cause notable changes in the BIPVs. From the Gaussian process regression model, all 3 BIPVs showed an increasing trend with each of the 3 scaling components, but the BIPVs are found to be most sensitive to the height dimension of the head. From the Sobol sensitivity analysis, the ${T_{zz}}$ scaling parameter contributes nearly 60% to the total variance in peak and average CMPS; ${T_{yy}}$ contributes approximately 20%, whereas ${T_{xx}}$ contributes less than 5%. The remaining contribution is from the 6 shearing components. Unlike peak and average CMPS, the VF-CMPS BIPV is associated with relatively evenly distributed Sobol indices across the 3 scaling parameters. Furthermore, the contribution of shearing components on the total variance in this case is negligible. Head shape has a considerable influence on the injury predictions of computational head injury models. Available "average" head injury models based on a 50th-percentile U.S. male are likely associated with considerable uncertainty. In general, larger head sizes correspond to greater BIPV magnitudes, which point to potentially a greater injury risk under rapid neck rotation for people with larger heads.