We develop an R package MLRMPA for fitting a pool of models between response variable and chemical descriptors. It is an embedded method combining feature selection with model building. The feature selection procedure is a cluster sampling method and different from model population analysis (MPA) that was implemented in a previously published study. The modeling process performs multiple stepwise regression analysis using the sampled features from the clustered group. This paper provides the algorithm and method implemented in the R package, which includes VarCor feature selection, cluster sampling, model building and model checking. This package is applied to establish an optimal linear model to predict the response and detect outliers from sub-optimal models.