Suppose you have a system with physical field $\sigma(\vec{r})$, which is constrained by some law $L[\sigma](\vec{r})=0$. Naively, to do statistical mechanics, one has to carry around the constraint $L[\sigma]=0$ in all computations.
In some cases, in particular when $L$ is linear and built out of derivatives, the constraint can be explicitly solved by writing $\sigma = K[\phi]$, in terms of a new field $\phi$. This allows us to work directly on the constrained manifold. The price of this representation is that this new field is, in general, a gauge field. This gives a general mechanism for gauge theory in statistical mechanics.
For example, if $\sigma=\sigma_{ij}$ is the stress tensor then the condition of mechanical equilibrium is $\partial_i \sigma_{ij}=0$. This can be solved (in 3 dimensions) by writing $\sigma_{ij} = \epsilon_{ikl} \partial_k \phi_{lj}$. The new field $\phi$ cannot directly be a physical field, since we can add any function $\Delta \phi_{lj}$ for which $\epsilon_{ikl} \partial_k \Delta \phi_{lj}=0$, such as constants. We are thus led to consider the gauge theory for $\phi$. The original physical Lagrangian, once written in terms of $\phi$, will have Wilson-loop-like terms, after using Stokes' theorem. For example, the surface integral $\int_{A} dS n_i \sigma_{ij} = \int_{\partial A} dl_i \phi_{ij}$, where $A$ is some surface with boundary $\partial A$.
The example here is taken from the physics of amorphous solids at low temperature, for which the states of mechanical equilibrium control the vibrational and rheological properties of the solid. In this case one needs to also impose torque balance, so $\phi$ is further rewritten in terms of a more primitive field, but the result is the same: gauge theory.