Since the [OP asked for a discussion of features](https://mathoverflow.net/questions/440271/higher-integrability-for-sobolev-functions#comment1135653_440271), I provide one by way of an explanation of Christian Remling's counterexample:

Holder's inequality states that 
$$ | \int fg | \leq \| f\|_p \|g\|_q $$
if $p^{-1} + q^{-1} = 1$, and the dual characterization of $L^p$ states that 
$$ f \in L^p \iff \exists M \forall g\in L^q, \int fg \leq M \|g\|_q $$

Your question can be interpreted as "does the dual characterization of $L^p$ work if, instead of all $g$ in $L^q$, we only look at those $g$s which are characteristic functions of balls (here we use that for $g$ the characteristic function of $B(x,r)$ in $\mathbb{R}^n$, its $L^q$ norm is $\approx r^{n/q}$). 

While one may expect that the set of characteristic functions of balls form a nice "basis" for $L^q$, the problem however is that we don't see some convexity properties: knowing that the bounds hold for $g_1$ and $g_2$ only gives

$$ \int f (g_1 + g_2) \leq M (\|g_1\|_q  + \|g_2 \|_q) $$

the right hand side of which is not $\leq \|g_1 + g_2\|_q$. In Christian Remling's answer his $f$ is built out of functions 
$$ \chi_{B(0, 2^{-k} + 2^{-k^2 - k})} - \chi_{B(0,2^{-k})} $$
We have for this function, in $\mathbb{R}^n$
$$ \|\chi_{B(0, 2^{-k} + 2^{-k^2 - k})} - \chi_{B(0,2^{-k})}\|_q \approx 2^{-nk/q - k^2/q} \ll 2^{-nk/q} \approx \|\chi_{B(0, 2^{-k} + 2^{-k^2 - k})}\|_q + \|\chi_{B(0,2^{-k})}\|_q $$
These functions are, in some sense, "orthogonal" to the characteristic functions of balls, and so their largeness are not detectable by integration against balls.