The proof below uses only the chain rule (to change integration variables). It's what I, as a physics student, would write.
But first, I believe that you have an error in your original posting. I believe that the 1D identity should be:
\begin{align}
\int \Big( \sum_{\substack{\textrm{roots $a_i$}\\ \textrm{ of $g$}}} \frac{1}{g'(a)} \Big)^{-1} \delta(g(a)) da = 1.
\end{align}
I know that you asked about the multidimensional case. However, let's review the one-dimensional case first. They are extremely similar.
One-dimensional Case
\begin{align}
\int_{\substack{\textrm{small $a$-region} \\ \textrm{ containing root $a_i$}}} \delta (g(a)) da = &\int_{\substack{\textrm{small $g$-region} \\ \textrm{ corresponding to} \\ \textrm{ root $a_i$}}} \delta(g) \frac{1}{g'(a)} dg
\end{align}
We knew $g$ as a function $a$ (i.e., $g(a)$), but we can also think of $a$ as a function of $g$. This works because in the small region near the root, $g(a)$ is invertible. So we can think of $\frac{1}{g'(a)}$ as a function of $g$. The delta function extracts the value of this function near $g=0$. But in the small region we care about, $g=0$ means $a=a_i$. Thus:
\begin{align}
\int_{\substack{\textrm{small $a$-region} \\ \textrm{ containing root $a_i$}}} \delta (g(a)) da = \frac{1}{g'(a_i)}
\end{align}
Integrating over all space instead of a small region gives:
\begin{align}
\int \delta(g(a)) da = \sum_{\textrm{roots $a_i$ of $g$}} \frac{1}{g'(a_i)}
\end{align}
which is a special case of an identity you can find at the Wikipedia article on delta functions
We can rewrite this to obtain:
\begin{align}
\int \Big( \sum_{\substack{\textrm{roots $a_i$}\\ \textrm{ of $g$}}} \frac{1}{g'(a)} \Big)^{-1} \delta(g(a)) da = 1.
\end{align}
Which is (except for the inversion) the 1D identity that you had. You can bring the factor $\left( \sum_{\textrm{roots $a_i$ of $g$}} \right)^{-1}$ inside the sum because it's just a constant.
Multidimensional case
The only difference in the multidimensional case is that when you change variables you need to use the Jacobian determinant (which Peskin and Schroeder call $\frac{\partial g_i }{\partial a_j}$). So, instead of $da = \frac{1}{g'(a)} dg$, we have
\begin{align}
\prod_k d a_k = \frac{1}{\det{ \frac{\partial g_i }{\partial a_j} } } \prod d g_k
\end{align}
Note that $\frac{\partial g_i }{\partial a_j}$ really stands for ``the matrix whose $ij$ entry is $\frac{\partial g_i }{\partial a_j}$''.
Hence the final identity from before becomes:
\begin{align}
\left( \prod_k \int d a_k \right) \Big( \sum_{\substack{\textrm{roots $\vec{a}_i$}\\ \textrm{ of $g$}}} \frac{1}{\det{ \frac{\partial g_i }{\partial a_j} } } \Big)^{-1} \delta^{(n)}(\vec{g}(\vec{a})) = 1.
\end{align}
where $\vec{a}_i$ are the roots of $\vec{g}(\vec{a})$ and each determinant is evaluated at the root $\vec{a}_i$.
If there's only one root
, then the sum only has one term, and so inverting the terms is undone by inverting the sum:
\begin{align}
\left( \prod_k \int d a_k \right) \det{ \frac{\partial g_i }{\partial a_j} } \delta^{(n)}(\vec{g}(\vec{a})) = 1.
\end{align}
[Which is the identity on p295 of Peskin and Schroeder.] The fact that these two inversions cancel is nice, but I think it is what led you to the error in your original post.