This is not a direct answer, but it may be useful.

Define, for any subset of sites $A$ $$ \sigma_A \equiv \prod_{i\in A} \sigma_i $$

For systems with only ferromagnetic couplings (ferromagnetic is defined as the couplings being positive as defined below, regardless of whether you preserve global spin-flip symmetry. For example, a magnetic field that prefers $\sigma_i = 1$ is allowed.) (note that there is no disorder averaging), $$ H = -\sum_{A} J_{A} \sigma_A ,\quad J_{A} \geq 0 $$ there are the Griffiths-Kelly-Sherman inequalities and its extension, the Ginibre inequality. One of its statements is under these conditions, $\langle \sigma_A\sigma_B\rangle -\langle \sigma_A\rangle \langle \sigma_B\rangle \geq 0 $. This is useful for us because if you consider \begin{align} \frac{\partial \langle \sigma_A \rangle}{\partial J_B} &= \frac{\partial}{\partial J_B}\Bigg(\frac{1}{\sum_{\sigma} e^{-H[\sigma]}} \sum_{\sigma}\sigma_A e^{-H[\sigma]}\Bigg) \\ &= \langle \sigma_A\sigma_B\rangle -\langle \sigma_A\rangle \langle \sigma_B\rangle \geq 0 \end{align} using the above theorem. It tells us that any correlation function is non-decreasing in any of the "ferromagnetic" couplings.

This is not a direct answer, but it may be useful.

Define, for any subset of sites $A$ $$ \sigma_A \equiv \prod_{i\in A} \sigma_i $$

For systems with only ferromagnetic couplings (ferromagnetic is defined as the couplings being positive as defined below, regardless of whether you preserve global spin-flip symmetry. For example, a magnetic field that prefers $\sigma_i = 1$ is allowed.) (note that there is no disorder averaging), $$ H = -\sum_{A} J_{A} \sigma_A ,\quad J_{A} \geq 0 $$ there are the Griffiths-Kelly-Sherman inequalities and its extension, the Ginibre inequality. One of its statements is under these conditions, $\langle \sigma_A\sigma_B\rangle -\langle \sigma_A\rangle \langle \sigma_B\rangle \geq 0 $. This is useful for us because if you consider \begin{align} \frac{\partial \langle \sigma_A \rangle}{\partial J_B} &= \frac{\partial}{\partial J_B}\Bigg(\frac{1}{\sum_{\sigma} e^{-H[\sigma]}} \sum_{\sigma}\sigma_A e^{-H[\sigma]}\Bigg) \\ &= \langle \sigma_A\sigma_B\rangle -\langle \sigma_A\rangle \langle \sigma_B\rangle \geq 0 \end{align} using the above theorem. It tells us that any correlation function is non-decreasing in any of the "ferromagnetic" couplings. In particular, if your $c_{ij}\geq 0$, you can show $$ \frac{\partial \langle \sigma_k \sigma_l \rangle}{\partial c_{ij}} \geq0 \quad \text{for all } i,j,k,l. $$ I am not sure if there are additional subtleties with spontaneous symmetry breaking, but I don't think so since this correlation function is invariant under the global symmetry.