[Purposes] The paper aims to determine the optimal model for predicting the resilient modulus of subgrade soil and predict the resilient modulus of subgrade soil accurately. [Methods] The fully Bayesian Gaussian process regression (fB-GPR) approach is proposed in this study to establish the relationship between the resilient modulus of subgrade soil (MR) and the confining pressure (σ3), deviator pressure (σd), moisture content (w), dry density (ρd). The parameters of the GPR model are estimated accurately and the optimal model for the prediction of the resilient modulus of subgrade soil is determined objectively using the proposed approach in this study. The optimal model reaches a balance between the complexity and fitting degree. [Findings] The findings show that the resilient modulus of subgrade soil can be predicted accurately using the proposed approach in this study. The coefficient of determination (R2) and the mean absolute percentage error (RMAPE) of the optimal GPR model reach 0.99 and 1.51%, respectively, which is similar to the full model. In 100 experiments, the percentage of the optimal model being selected reaches 88%. [Conclusions] The proposed fB-GPR approach not only accurately predict the resilient modulus using the easily available indices of subgrade soil but also effectively eliminate the redundant input variables. The proposed approach reduces complexity without compromising the fitting degree of the GPR model, holding significance for its application and promotion in the future.