Two-vehicle crashes have been dominating all types of traffic accidents, wherein the vehicle drivers have been sustaining the highest risk of injury among all vehicle occupants. To understand the critical factors to the drivers’ injury severity of two-vehicle crashes, we employed the random parameters multinomial logit model as a data analyzing tool. To capture the unobserved heterogeneity and potential temporal instability, we combined two strategies: Bayesian random parameter logit and explicitly correlated outcomes. The random parameter logit models were validated with a nine-year large-scale dataset compiled by combining the Crash Report Sampling System (CRSS) and General Estimates Sampling (GES) databases. The results underscore the importance of explicit modeling of inter-outcome correlation, which captured the potential transition probability between adjacent levels of injury severity and improved the model’s predictability. Our model also highlighted substantial per-case and per-driver heterogeneity, which respectively explained 22.8% and 29.4% of the total variance (minor injury) and 25.4% and 24.9% of the variance (severe injury). We found that the female drivers, old ( years) drivers, unbuckled drivers, speeding drivers sustained a higher injury risk in their corresponding groups. Drivers in lighter and older vehicles suffer higher injury risks. Several other factors also considerably affect the injury severity outcomes, such as the road’s speed limit and variables that are proxies of traffic volume (intersection type, whether at the peak hours). Regarding Bayesian modeling, we observed that using weakly informative prior distribution has little effect on the parameter estimates. We also pointed to the directions to further improve the proposed modeling framework.