Interaction selection for linear regression models with categorical predictors is useful in many fields of modern science, yet very challenging when the number of predictors is large. Existing interaction selection methods focus on finding one optimal model. While attractive properties such as consistency and oracle property have been well established for such methods, they actually may perform poorly in terms of stability for high-dimensional data, and they do not deal with categorical predictors. In this paper, we introduce a variable importance based interaction modeling (VIBIM) procedure for learning interactions in a linear regression model with both continuous and categorical predictors. It delivers multiple strong candidate models with high stability and interpretability.  We apply the VIBIM procedure to a Corona Virus Disease 2019 (COVID-19) data used in Tian et al. (2020) and measure the effects of relevant factors, including transmission control measures on the spread of COVID-19. We show that the VIBIM approach leads to better models in terms of interpretability, stability, reliability and prediction, compared to the models in Tian et al. (2020) and by some group variable selection methods.