The HOMO-LUMO gap is only weakly correlated with the reactivity of the molecule. Many sources, however, only point out the correlation between the HOMO-LUMO gap and reactivity, but fail to point out that the correlation only holds in the statistical sense, and is by no means universal. In addition, among those people who ever attempted to compare the HOMO-LUMO gap against the reactivity, usually only those people who do find a correlation will publish their results, leading to considerable survivorship bias. As a result, more and more people try to use the HOMO-LUMO gap to rationalize their reactivity patterns, while not realizing that this does not work nearly as well as they thought.
The correlation between the HOMO-LUMO gap and reactivity stems from second-order perturbation theory. A molecule is said to be reactive, if either (1) it reacts exothermically and quickly with many kinds of molecules, or (2) it decomposes exothermically and quickly by itself. These two cases can be unified as: a reactive molecule is one that tends to significantly lower its energy in the presence of a perturbation, where the perturbation can be either external (from other molecules) or internal (from the structural change of the molecule itself).
Suppose that the ground state wavefunction of the molecule is $|\Psi_0\rangle$, the Hamiltonian of the molecules is $H$, and the perturbation is expressed as an operator $V$. Then from second-order perturbation theory, the (ground state) energy change caused by the perturbation is
\begin{align}
\Delta E & = \langle \Psi_0 | V | \Psi_0 \rangle - \sum_{I>0} \frac{|\langle \Psi_0 | V | \Psi_I \rangle|^2}{\omega_I} + O(V^3) \tag{1}
\end{align}
where $|\Psi_I\rangle, I=1,2,\ldots$ is the $I$-th excited state of the molecule, and $\omega_I$ is the excitation energy of the state $|\Psi_I\rangle$.
Now, the usual assertion "smaller HOMO-LUMO gap means higher reactivity" amounts to making the following 4 approximations to Eq. (1):
$\langle \Psi_0 | V | \Psi_0 \rangle$ is negligible. This means that $V$ cannot be too small, since the second term on the right hand side of Eq. (1) scales quadratically with respect to $V$, while the first term scales only linearly. But ironically, perturbation theory itself requires that $V$ cannot be too large, otherwise the theory does not hold. Therefore the assertion can only possibly hold for those reactivities where the molecule undergoes a not-too-large, but also not-too-small change.
In the sum $\sum_{I>0} \frac{|\langle \Psi_0 | V | \Psi_I \rangle|^2}{\omega_I}$, the $I=1$ term dominates. This can be partly justified by the fact that $\omega_1 < \omega_2 < \omega_3 < \ldots$, which means that the $I=1$ term has the smallest denominator. But there is no guarantee that it has the largest numerator! The numerator, $|\langle \Psi_0 | V | \Psi_1 \rangle|^2$, can even be zero, for example when $V$ is a totally symmetric operator but $|\Psi_0 \rangle$ and $|\Psi_1 \rangle$ belong to different irreducible representations. Only when the numerator $|\langle \Psi_0 | V | \Psi_1 \rangle|^2$ is not too small, can the $I=1$ term possibly dominate over the remaining terms of the series.
$\omega_1 \approx \epsilon_{LUMO} - \epsilon_{HOMO}$. In other words the first excitation energy is approximated as the HOMO-LUMO gap. This approximation is not a good one, with errors on the order of 1 eV or so, because the approximation neglects the exciton binding energy. But this may be the best approximation among all the four approximations, as the other approximations are even poorer.
$|\langle \Psi_0 | V | \Psi_1 \rangle|^2$ is approximately constant across the molecules that you want to compare. This can only be justified if your molecules are very similar (i.e. have similar $|\Psi_0 \rangle$ and $|\Psi_1 \rangle$), and possess similar reactivities (i.e. have similar $V$). However, similar molecules also tend to possess similar $\omega_1$, so the differences of $|\langle \Psi_0 | V | \Psi_1 \rangle|^2$ among the molecules may still overshadow the differences of $\omega_1$. You have to make sure that the molecules differ considerably in the HOMO-LUMO gap, but do not differ much in other aspects, including the shapes of HOMO and LUMO (which affect the wavefunctions $|\Psi_0 \rangle$ and $|\Psi_1 \rangle$).
Only with these four drastic approximations can we arrive at our result
\begin{align}
\Delta E & \propto \frac{1}{\epsilon_{LUMO} - \epsilon_{HOMO}} \tag{2}
\end{align}
i.e. the reactivity is inversely proportional to the HOMO-LUMO gap.
The obvious conclusion at this point is that, while the reactivity may have a negative statistical correlation with the HOMO-LUMO gap, the correlation can only be relied upon under heavily controlled conditions. The reactivity must not be too low, nor can it be too high; the matrix element $|\langle \Psi_0 | V | \Psi_1 \rangle|$ must be relatively large, at least it must not be zero; the exciton binding energy of the first excited state is either small or systematic, so that the trend of $\epsilon_{LUMO} - \epsilon_{HOMO}$ reflects the trend of $\omega_1$; and finally $|\langle \Psi_0 | V | \Psi_1 \rangle|$ must be roughly constant among the molecules to be studied. If any of the conditions are violated, then you can never explain your reactivity based on the HOMO-LUMO gap. I said never, because even if the HOMO-LUMO gap does correlate perfectly with the reactivity, you are achieving this correlation largely by luck, i.e. by exploiting an accidental cancellation of the contributions of $\langle \Psi_0 | V | \Psi_0 \rangle$, $|\langle \Psi_0 | V | \Psi_I \rangle|$, the $O(V^3)$ term, and the exciton binding energy.
Thus my suggestion is: observe if your HOMO-LUMO gap data do correlate with your reactivities. If yes, you can comment that the correlation exists, or that the trend of the HOMO-LUMO gaps contributes to the trend of the reactivities, but unless you provide evidence that the aforementioned four approximations hold, you should not say that you "explained" the reactivities by the HOMO-LUMO gaps. If no, then simply move on, and try other descriptors (atomic charges/spin populations/bond orders at transition states, for example) instead - there wasn't a strong reason for you to expect a good correlation in the first place.
