Bayesian hierarchical models (BHM) are being routinely used for astronomy data. However, with the recent advent of computing power, although a lot of complex statistical models can be fitted using Monte Carlo methods, it has largely remained illusive how to validate these complex models when the data are observed with heterogeneous large measurement errors. We illustrate the methodology using a non-trivial extension of the M–R relation by including the incident flux as an additional variable. By using BHM that leverages the flexibility of finite mixture models, a probabilistic mass–radius–flux relationship (M–R–F relation) is obtained based on a sample of 319 exoplanets. We find that the flux has non-negligible impact on the M–R relation, while such impact is strongest for hot Jupiters. On the population level, the planets with higher level of flux tend to be denser, and high flux could trigger significant mass loss for plants with larger radii. We present two novel methods to examine model assumptions, which can be used not only for the M-R-F models but can also be adapted for other statistical models.