The short answer to the question is No. The consequence of Lorentz invariance is that it isn't possible to impose, *which* of the four components of the gauge fields is the ``temporal'' component and which the ``longitudinal'' component: the *only* invariant information is that *one* combination is temporal and one is longitudinal. Similarly, gauge invariance implies that it isn't possible to invariantly distinguish *which* components are gauge artifacts-only which are gauge-invariant. The former change as the gauge-fixing condition is changed, the latter don't. The only invariant information is that two are physical. The Ward-Takahashi (Slavnov-Taylor) idenities between the correlation functions describe what these statements mean in the quantum theory, beyond the equations of motion. This is well-known and can be found in any textbook on quantum field theory. If the gauge fixing condition is not covariant, then the fact that these identities hold implies that Lorentz invariance and gauge invariance are, in fact, restored, which means that the contribution of the gauge fixing terms cancels the contributions of the ghosts and only physical terms contribute to gauge-invariant processes. If the gauge fixing condiion is covariant and the regularization scheme is, as wsell, then Lorentz invariance is manifest and the identities constrain gauge invariance. No particular gauge fixing condition can be singled out. All gauge fixing conditions lead to the contribution of the physical degrees of freedom to gauge-invariant observables.

Gauge nvariance, in particular, implies Gauss' law on the gauge invariant configurations.

The quantization of non-abelian gauge theories is understood for more than forty years now. This includes, of course, an understanding of how gauge-variant quantities change, as one changes the gauge fixing condition. For a review, cf.

E. S. Abers and B. W. Lee,

``Gauge Theories,'', Phys.\ Rept.\ {\bf 9} (1973) 1.

and B. W. Lee's  course at Les Houches:

http://lss.fnal.gov/archive/1976/pub/Pub-76-034-T.pdf

This subject is a standard topic and exercise in quantum field theory and can be find in the textbooks these days.

The path integral formalism for compact gauge groups is well-defined through the lattice formulation-where the measure  is well-defined, the action is too. The discussion here makes sense for non-gravitational theories only, of course. It isn't known how to define the measure of the path integral, unambiguously, for non-compact gauge groups. The two issues don't have anything to do with each other and there isn't any point in mixing them up.

Well lattice field theory does give rise to physical quantities, that can be compared to experiment and is mathematically well-defined. The Haar measure on the lattice is on the group, not the algebra and thus the path integral of the pure gauge theory, with compact gauge group, which is of concern here, defined on the lattice is invariant under  large gauge transformations. This work is reviewed in the Physics Reports article by J. M. Drouffe and J. B. Zuber, for instance.

The operator formalism has been studied-cf., for instance, the analysis by Henneaux and Teitelboim.``Quantization of gauge systems,''

Princeton, USA: Univ. Pr. (1992) 520 p 

and of  P.C.Nelson and L.Alvarez-Gaume,``Hamiltonian Interpretation of Anomalies,''

Commun.\ Math.\ Phys.\ {\bf 99} (1985) 103.

The phase space of a field theory is infinite-dimensional and requires careful regularization-and the treatment of constraints even more care. These topics don't have *anything* to do with the subject under discussion here, however. 

No, meson masses are the *output* of lattice QCD simulations, not the *input*. And, more generally, the renormalized quantities are the output and correspond to the values calculated in the appropriate scaling limit. (Of course masses, as such, don't mean anything, absent a scale,  Lambda_QCD, for instance. So the mass ratios ar the physically intetesting quantities.) This is well-understood, theoretically and computationally; of course it is challenging to approach the scaling limit,  in practice, for certain observables-more so than for others and this is an active research area. But this doesn't have anything to do with the topic of the question, as such. 

And belief doesn't enter here, the point is that it is possible to perform well-defined calculations, within controlled approximations. Of course these approximations may not allow, yet,  the mathematical proof of certain statements. But the issue of the discussion isn't among them.

 The lattice QCD Lagrangian is determined by the (lattice) masses of the quarks and the lattice coupling constant (and the size and shape of the lattice).  From these it's possible to compute all correlation functions of  quark and gluon fields, by numerical techniques. The physical, i.e. renormalized, masses of the mesons, baryons etc. follow from such a calculation  and the mass ratios too: they are all defined as the values of appropriate combinations of  correlation functions of quark fields, evaluated in the scaling region. This is a difficult task in practice, but not in principle. The scaling limit amounts to tuning the lattice coupling constant to its critical value (for QCD  this is zero) and computing the renormalized values in its vicinity. In the scaling region the values of the lattice masses, as well as the properties of the lattice (spacing and shape) become, literally and technically, irrelevant and the mass ratios, in particular, calculated there are physical quantities.

They are,of course, subject to uncertainties, statistic and systematic and a lot of work is devoted to reducing them-it's not, just, a question of more computing power; a significant portion is devoted to imagining improved algorithms and better coding for them. Taking account of ``corrections to scaling'' is a vital part of working in this and related areas of computational physics.

The quenched approximation was, essentially, completed around the beginning of the 21st century and now a lot of work is devoted to ``flavor physics'', with ``dynamical'' quarks, in different contexts.   This is well known and reviewed in many places. But this doesn't have  anything to do with the topic discussed here. 

Once *one* mass scale is set, *all* mass ratios are determined by it, since no phase transition separates the different bound states.(Of course, at a given scale, it's not possible to observe *all* bound states simultaneously.)

And, once more, it's the ratios that are the physically relevant quantities. So, no, there isn't any data adjustment. The tuning of the quark masses means tuning the *lattice* quark masses to the scaling region, defined by the coupling constant.