Cross-frequency coupling is emerging as a crucial mechanism that coordinates the integration of spectrally and spatially distributed neuronal oscillations. Recently, phase-amplitude coupling, a form of cross-frequency coupling, where the phase of a slow oscillation modulates the amplitude of a fast oscillation, has gained attention. Existing phase-amplitude coupling measures are mostly confined to either coupling within a region or between pairs of brain regions. Given the availability of multi-channel electroencephalography recordings, a multivariate analysis of phase amplitude coupling is needed to accurately quantify the coupling across multiple frequencies and brain regions. In the present work, we propose a tensor based approach, i.e., higher order robust principal component analysis, to identify response-evoked phase-amplitude coupling across multiple frequency bands and brain regions. Our experiments on both simulated and electroencephalography data demonstrate that the proposed multivariate phase-amplitude coupling method can capture the spatial and spectral dynamics of phase-amplitude coupling more accurately compared to existing methods. Accordingly, we posit that the proposed higher order robust principal component analysis based approach filters out the background phase-amplitude coupling activity and predominantly captures the event-related phase-amplitude coupling dynamics to provide insight into the spatially distributed brain networks across different frequency bands.