We numerically study the eigenvalue statistics along the thermal to many-body localization (MBL) transition in random spin systems, and compare them to two random matrix (RM) models, namely the mixed ensemble and Gaussian β ensemble. We show both RM models are capable of capturing the lowest-order level correlations during the transition, while the deviations become non-negligible when fitting higher-order ones. Specifically, the mixed ensemble will underestimate the longer-range level correlations, while the opposite is true for β ensemble. Strikingly, a simple average of these two models gives nearly perfect description of the eigenvalue statistics at all disorder strengths, even around the critical region, which indicates the interaction range and strength between eigenvalue levels are responsible for the phase transition.