Bhatnagar-Gross-Krook (BGK) models are widely used to study rarefied gas dynamics. However, as simplified versions of the Boltzmann collision model, their performances are uncertain and need to be carefully investigated in highly nonequilibrium flows. In this study, several common BGK models, such as the Ellipsoidal Statistical BGK (ES-BGK) and Shakhov BGK (S-BGK) models, are theoretically analyzed using their moment equations. Then, numerical comparisons are performed between the Boltzmann collision model and BGK models based on various benchmarks, such as Fourier flow, Couette flow, and shock wave. The prediction performance of the ES-BGK model is better than that of the S-BGK model in Fourier flow, while prediction performance of the S-BGK model is better than that of the ES-BGK model in Couette flow and shock wave. However, with increasing Knudsen number or Mach number, the results of both ES-BGK and S-BGK deviate from the Boltzmann solutions. These phenomena are attributed to the incorrect governing equations of high-order moments of BGK models. To improve the performance of the current BGK models, the S-BGK model is extended by adding more high-order moments into the target distribution function of the original one. Our analytical and numerical results demonstrate that the extended S-BGK (S-BGK+) model provides the same relaxation coefficients as the Boltzmann collision model for the production terms of high-order moment equations. Compared to the other BGK models, the proposed S-BGK+ model exhibits better performance for various flow regimes.