In magnetic resonance elastography (MRE), displacement fields from shear waves are inverted to estimate underlying material properties. Modulus differences detected by MRE may be used to distinguish tumors or other localized pathology in tissue. The accuracy of modulus estimates depends on the choice of the assumed constitutive model, as well as on the inversion algorithm, image resolution, and signal-to-noise ratio. In particular, in simpler inversion methods such as direct inversion and three-dimensional local frequency estimation (3D-LFE) the constitutive model is minimal (linear, elastic or viscoelastic, and isotropic) and the simplifying assumption of local homogeneity is usually made. The assumption of local homogeneity is often inaccurate , since the shear wavelength is typically comparable to the size of the structures of interest. Notably, the residual error (in direct inversion) between the model and the experimental data increases sharply at the boundaries of inclusions, while the “certainty” of the 3D-LFE estimate decreases. These error metrics may be used to detect local stiffness heterogeneity, as well as indicate variations in appropriate constitutive models. The utility of model uncertainty is demonstrated in simulations and with MRE data from a heterogeneous gel phantom.