Elastography aims at estimating the mechanical properties of biological tissue by solving an inverse problem based on measured tissue deformations, enabling non-invasive disease diagnosis. In practice, this inference relies on computational models whose governing equations —including constitutive laws — are typically assumed to be correct. Such assumptions can give rise to a dangerous illusion of fit: inverse solvers may identify material parameters that match as best as possible observations even when the underlying constitutive model is structurally inadequate, without providing any indication of this mismatch. In clinical contexts, this unrecognized model error can compromise predictive reliability; for example, a stiff tumor may be inferred as substantially softer when a linear elastic isotropic model is incorrectly imposed. We present a Bayesian computational framework designed to explicitly address constitutive uncertainty within inverse biomechanical problems. The method separates well-established physical principles, such as conservation of linear momentum, from uncertain constitutive relationships. In contrast to classical model discrepancy approaches, which typically rely on global additive correction terms, our framework probabilistically identifies spatially localized regions where constitutive assumptions fail. This is achieved by introducing latent variables—including stress fields that are not directly observable—and enforcing governing equations through virtual likelihoods rather than hard constraints. In doing so, the approach maintains mechanical consistency while enabling systematic quantification of model-form uncertainty. The framework builds upon the Weak Neural Variational Inference (WNVI) methodology, retaining its forward-model-free formulation and computational efficiency. The proposed approach is demonstrated on synthetic, brain-like geometries that emulate key challenges encountered in biomechanical inverse problems, including heterogeneous material behavior, noisy measurements, sparse observation settings, and limited displacement components. These in silico studies show that the framework reconstructs quantities of interest together with associated uncertainty, while explicitly identifying regions where the assumed constitutive model is not supported by the data. In this way, the method prevents the illusion of fit and provides a computational testbed for assessing model adequacy under controlled conditions.