We present a diagrammatic theory for determining consistent electromagnetic response functions in strongly correlated fermionic superfluids. While a gauge-invariant electromagnetic response is well understood at the BCS level, a treatment of correlations beyond BCS theory requires extending this theoretical formalism. The challenge in such systems 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 method, which determines these collective modes in the presence of correlation effects and in a manner compatible with gauge invariance. Since this method is based on fundamental aspects of quantum field theory, the underlying principles are broadly applicable to strongly correlated superfluids. As an illustration of the technique, we apply it to a simple class of theoretical models that contain a frequency-independent order parameter. These models include BCS-BEC crossover theories of the ultracold Fermi gases, along with models specifically associated with the high-Tc cuprates. 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 those of strict BCS theory, without any modification from additional correlations. As the path integral 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.