The ghosts are not so much inserted, as they naturally arise. The path integral of a gauge theory naively defined will integrate over all fields, including those related by a gauge symmetry, which are seen by the theory as being equivalent.
The Faddeev-Popov procedure provides a means to split our integration over physically distinct configurations and those over gauge orbits. Consider the case of non-Abelian gauge theory, with,
$$\int \mathcal D[A] \exp \left[i \int d^4x \left( -\frac14 (F^a_{\mu\nu})^2\right) \right].$$
To integrate only over physically distinct configurations, we need to constrain the integral by a gauge-fixing procedure, $G(A) = 0$ in general. To fix $G(A) = 0$, we can use a delta function, but to do so, we need to take into account the appropriate Jacobian factor, so the identity is,
$$1= \int D[\alpha(x)] \delta(G(A^\alpha)) \det \frac{\delta G(A^\alpha)}{\delta \alpha}$$
where $A^\alpha$ is the field transformed, that is, $(A^\alpha)^a_\mu = A^a_\mu + g^{-1}D_\mu \alpha^a$. We then have the path integral,
$$\int \mathcal D[A] \, e^{iS[A]} = \left(\int \mathcal D[\alpha(x)] \right) \int \mathcal D[A] e^{iS[A]} \delta(G(A))\det \frac{\delta G(A^\alpha)}{\delta \alpha}.$$
We essentially factored it into integrations over the gauge orbits and physically distinct solutions. Now, for an $n \times n$ matrix, $M$, we can express the determinant as a Grassmann integral, namely,
$$\int e^{-\theta^T M \eta} d\theta d\eta = \det M$$
where we have vectors of Grassmann variables, $\theta$ and $\eta$. Going back to the path integral, the determinant is the determinant of a differential operator, and so we use an analogous formula to compute it. We then interpret the analogous $\theta$ and $\eta$ as being fields, or ghosts.
To put it yet another way, we essentially introduced dummy variables in order to express the determinant as an integral, and it turns out this integral when included has the same form as a Lagrangian for these variables, and so we can interpret them as fields.