Renormalization group methods are well-established tools for the (numerical) investigation of the low-energy properties of correlated quantum many-body systems, allowing to capture their scale-dependent nature. The functional renormalization group (FRG) allows to continuously evolve a microscopic model action to an effective low-energy action as a function of decreasing energy scales via an exact functional flow equation, which is then approximated by some truncation scheme to facilitate computation. Here we report on our transcription of a recently developed multiloop truncation approach for electronic FRG calculations to the pseudo-fermion functional renormalization group (pf-FRG) for interacting quantum spin systems. We discuss in detail the conceptual intricacies of the flow equations generated by the multiloop truncation, as well as our essential refinements to the integration scheme for the resulting integro-differential equations. To benchmark our approach we analyze geometrically frustrated Heisenberg models on the simple cubic, face-centered cubic, and pyrochlore lattice, discussing the convergence of physical observables for higher-loop calculations and comparing with exact quantum Monte Carlo results where available. Combined, these methodological refinements systematically improve the pf-FRG approach to one of the numerical tools of choice when exploring frustrated quantum magnetism in higher spatial dimensions.