Ttdfjt

When doing DFT calculations, some integrations are commonly done numerically on grids. [In fact, more than a single grid may be used at the same time for different integrals, e.g. approximations such as the resolution of identity (RIJCOSX, RIJK, etc., see e.g. J. Chem. Phys. 118, 9136 (2003)) use grid schemes too (I believe some programs call this approximation density fitting).]

There's the downside that those grid schemes introduce a source of error. In fact, I find that only by increasing the quality of the grid I can remove some imaginary frequencies. Can I rigorously compare energy values among calculations that used slightly different grid schemes?

I think yes and I reasoned as follows. Since I'm interested in energy differences, let's imagine two structures with true energies $E^*_1$ and $E^*_2$. Both energies are approximated by the calculated energies $E_i = E^*_i + epsilon_i$, where $epsilon_i$ is the error introduced by calculations (including grids, etc.). Now the energy difference $Delta E^* = E^*_2 - E^*_1$ is approximated by $Delta E = E_2 - E_1 = Delta E^* + Delta epsilon$, where $Delta epsilon = epsilon_2 - epsilon_1$. Since grid schemes introduce errors that are smaller than other sources of error (e.g. implicit solvation), everything should be fine as long as $Delta epsilon$ is acceptably small for the particular application at hand.

Does this seem reasonable? If not, what's wrong with the above?

Yes, your error analysis is valid for energy differences. However, I believe it is also valid for the absolute error $epsilon_i$ of any calculated quantity, not just the error $Deltaepsilon$ of an energy difference.

In any numerical computation, the key thing for situations like this is to ensure that these "structural" sources of error have been reduced to a magnitude that don't affect your results. In this case, this is achieved when you reach a grid quality $Q$ where the numerical results no longer change when you further increase to $Q+delta Q$.

To be clear: in this case elimination of the imaginary frequencies is not a good criterion for a sufficiently large $Q$. The proper stopping point for refinement of $Q$ is when the results no longer change appreciably.

Note that, as you allude to implicitly in the question, this situation is different than the one of trying to compare, e.g., absolute energies from computations run with different methods or at different levels of theory. In this case, there are fundamental/theoretical reasons why those comparisons cannot be meaningfully made.