How to calculate coupling strength from a VNA?

I have a VNA (M980x) which I have used to calculate the S-parameters at many different frequencies for a resonator. As such, I have \[S_{11}\], \[S_{12}\], \[S_{21}\] and \[S_{22}\] data, and I would like to calculate the coupling strength at each port for a certain resonance mode. However, due to impedance mismatch, only \[S_{12}\] and \[S_{21}\] match up, as expected. However, \[S_{11}\], \[S_{21}\] and \[S_{22}\] and all have different resonant frequencies.

I was using \[k_i = S^2_{21}/(1-S^2_{ii} - S^2_{21})\] to calculate the coupling at the \[i^{th}\] port, however, this gives me nonsensical couplings (negative, no trend), primarily due to the fact that I do not which peak I should choose. Both the \[S_{21}\] and \[S_{ii}\] resonant frequencies give me nonsensical couplings.

I could not find an equation that takes impedance mismatch into account either. Am I missing something obvious? Even the formula above was used in a thesis without any citation.

I am sorry, I don't know how to fix the LaTeX here.