We present a general diagrammatic theory for determining consistent electromagnetic response functions in strongly correlated fermionic superfluids. The general treatment of correlations beyond BCS theory requires a new theoretical formalism not contained in the current literature. Among concrete examples are a rather extensive class of theoretical models which incorporate BCS-BEC crossover as applied to the ultra cold Fermi gases, along with theories specifically associated with the high-$T_c$ cuprates. The challenge is to maintain gauge invariance, while simultaneously incorporating additional self-energy terms arising from strong correlation effects. Central to our approach is the application of the Ward-Takahashi identity, which introduces collective mode contributions in the response functions and guarantees that the $f$-sum rule is satisfied. We outline a powerful and very general method to determine these collective modes in a manner compatible with gauge invariance. Finally, as an alternative approach, we contrast with the path integral formalism. Here, the calculation of gauge invariant response appears more straightforward. However, the collective modes introduced are essentially those of strict BCS theory, with no modification from correlation effects. Since the path integral scheme simultaneously addresses electrodynamics and thermodynamics, we emphasize that it should be subjected to a consistency test beyond gauge invariance, namely that of the compressibility sum-rule. We show how this sum-rule fails in the conventional path integral approach.