This work provides a family of explicit phenomenological models both in the F−H and F−B variable space. These models are derived directly from an analytical implicit homogenization model for isotropic magnetorheological elastomers (MREs), which, in turn, is assessed via full-field numerical simulations. The proposed phenomenological models are constructed so that they recover the same purely mechanical, initial and saturation magnetization and initial magnetostriction response of the analytical homogenization model for all sets of material parameters, such as the particle volume fraction and the material properties of the constituents (e.g., the matrix shear modulus, the magnetic susceptibility and magnetization saturation of the particles). The functional form of the proposed phenomenological models is based on simple energy functions with small number of calibration parameters thus allowing for the description of magnetoelastic solids more generally such as anisotropic (with particle-chains) ones, polymers comprising ferrofluid particles or particle clusters. This, in turn, makes them suitable to probe a large set of experimental or numerical results. The models of the present study show that in isotropic MREs, the entire magnetization response is insensitive to the shear modulus of the matrix material even when the latter ranges between 0.003-0.3MPa, while the magnetostriction response is extremely sensitive to the mechanical properties of the matrix material