: An accurate creep deformation model is needed for detailed description of creep behavior of high temperature structural materials, e.g., Grade 91 steels used in boiler tubes of thermal power plants. Two types of creep constitutive equations are known, as follows: the one, e.g., modified theta method, assumes the existence of a steady state; and the other, e.g., theta method, does not. So far, both types have been selected on a case by case basis and there is no consensus on whether or not the steady state should be assumed even if limited in the Grade 91 steels. In this study, we examine which creep constitutive equation, the modified theta method (steady-state type) or the theta method (non-steady-state type), is more appropriate in the Grade 91 steel based on data science. Here we calculated Bayesian free energy for each model with a linear regression framework by introducing hyperparameters for the exponential terms in the models. The modified theta method had lower Bayesian free energy than the theta method and exhibited extremely high posterior probabilities of nearly 1 for all test conditions. This means that the modified theta method is supported by the data and selected as a more appropriate model that should generated the data when the two models are compared. Note that the steady-state region accounted for more than a half of the creep rupture time for all test conditions. This shows that the steady state played a significant role in the creep deformation at the present test condition for the Grade 91 steel.