It depends on the values of $p$ and $N$ - sometimes there will be counterexamples and sometimes there will not. For instance, any function $f \in W^{1,q}$ for $q>N$ is classically differentiable almost everywhere, and hence the above result is true for all $1\leq p < \infty$. When $1\leq q < N$, then the result holds for all $1\leq p \leq q^\ast$, and this is optimal.

To see these claims, check Evan's proof of differentiability of $W^{1,q}$ for $q>N$ and observe that it is a clever use of Morrey's embedding to show that

$\frac{|f(y)-f(x)-\nabla f(x)(y-x)|}{|y-x|} \leq C \frac{1}{r^N}\int_{B(x,r)}|\nabla f(y)-\nabla f(x)|^q\;dy$,

and the right hand side goes to zero by the classical Lebesgue differentiation theorem for $\nabla f \in L^q$.

A similar proof allows one to use the Sobolev-Gagliardo-Nirenberg embedding to let $f \in W^{1,q}$ and consider the function $u(y)=f(y)-f(x)-\nabla f(x)(y-x)$ which has average value on a ball $\frac{1}{\alpha(N)r^N} \int_{B(x,r)} u(y)\;dy= \frac{1}{\alpha(N)r^N} \int_{B(x,r)} f(y)\;dy - f(x)$ to deduce that

\begin{align*}
\lim_{r\to 0}\frac{1}{r^{N+q^\ast}}\int_{B(x,r)} |f(y)-\frac{1}{\alpha(N)r^N}\int_{B(x,r)}f(y)dy-\nabla f(x)(y-x)|^{q^\ast}\;dy =0,
\end{align*}
and thus, writing
\begin{align*}
\frac{1}{r^{N+q^\ast}}&\int_{B(x,r)} |f(y)-f(x)-\nabla f(x)(y-x)|^{q^\ast}\;dy \newline
&\leq \frac{C}{r^{N+q^\ast}}\int_{B(x,r)} |f(y)-\frac{1}{\alpha(N)r^N}\int_{B(x,r)}f(y)dy-\nabla f(x)(y-x)|^{q^\ast}\;dy \newline
&\;\;+\left[\frac{\tilde{C}}{r^{N+1}}\int_{B(x,r)}|f(y)- f(x)-\nabla f(x)(y-x)|\;dy \right]^{q^\ast}
\end{align*}

and applying the standard theorem on $L^q$ differentiability the second term tends to zero. Thus, there is a possible improvement of differentiability theorems for Sobolev functions. To see that this is optimal, just take a function $f \in W^{1,q}$ and $f \notin L^p$ of any open set (for $p>q^\ast$. Then the above integral is plus infinity. For example, consider the function

$g_\epsilon(x) = |x|^{1-\frac{N}{q}+\epsilon}$.

Then $\int_{B(0,r)} |g_\epsilon|^p\;dx = |S^{N-1}|\int_0^1 r^{p(1-\frac{N}{q}+\epsilon)+N-1}\;dr$,

and $|g|^q$ will be not be integrable if

$p(1-\frac{N}{q}+\epsilon)+N-1 \leq -1$,

which is equivalent to

$p \geq \frac{Nq}{N-q-\epsilon q}$.

Now for any $p>\frac{Nq}{N-q}$, we can find an $\epsilon>0$ small such that

$p \geq \frac{Nq}{N-q-\epsilon q}$

and necessarily $g_\epsilon \in W^{1,q}$ by construction, since

\begin{align*}
\int_{B(0,r)} |\nabla g_\epsilon|^q\;dx \leq C \int_0^1 r^{-N+\epsilon q+N-1}\;dr <\infty
\end{align*}

Now, let $f_\epsilon(x):= \sum_n \frac{1}{2^n} g_\epsilon(x-x_n)$

for $\{x_n\}$ be a dense sequence of some bounded open set, and then we have

$\int_{B(x,r)}|f_\epsilon(x+h)-f_\epsilon(x)-\nabla f_\epsilon(x)h|^p\;dh = +\infty$,

for every $x$ in this set, since $f \notin L^p$ and the result is demonstrated.