In the Ising Model with $+$ boundary condition in dimension $2$, but possibly one could ask about dimension $d$. Set $\Lambda_N:=[-N,N]^2$, let $\beta >0$ be the inverse of the temperatura and $h>0$ be the magnetic field. For $x,y \in \Lambda_N$, is it true that
$$ \langle \sigma_x \sigma_y \rangle^{+}_{\Lambda_N,\beta,h} \le \langle \sigma_x \rangle^{+}_{\Lambda_N,\beta,h}. $$
I am particularly interested in the case where $x$ and $y$ are nearest neighbors, but I suppose that this case is actually more difficult to prove.
No, the inequality $$ \langle\sigma_x\sigma_y\rangle_{\Lambda_N,\beta,h}^+ \leq \langle\sigma_x\rangle_{\Lambda_N,\beta,h}^+ \tag{1} $$ does not hold in general, in any dimension and for any pair $x\neq y$ in $\Lambda_N$. Indeed, assume that $\beta<\beta_c$. On the one hand, $$ \lim_{h\downarrow 0}\lim_{N\to\infty} \langle\sigma_x\rangle_{\Lambda_N,\beta,h}^+ = 0 , $$ so that the right-hand side of (1) can be made as small as you wish by taking $N$ large enough and $h$ small enough.
On the other hand, by the Griffiths inequalities, $$ \langle\sigma_x\sigma_y\rangle_{\Lambda_N,\beta,h}^+ \geq \langle\sigma_x\sigma_y\rangle_{\Lambda_N,\beta,0}^+ \geq (\tanh\beta)^{\|y-x\|_1} . $$ The last inequality is obtained by choosing a shortest path in $\Lambda_N$ connecting $x$ and $y$ and removing all coupling constants (that is, setting them to zero) not along this path; this yields the 2-point function for a one-dimensional Ising model, which is easily seen to equal the expression in the right-hand-side.
This shows that the left-hand side of (1) is bounded below away from zero, uniformly in $N$ and $h$. Combined with the observation above, we conclude that (1) does not hold in general.