I realise I'm bumping into you again and already gave you an answer elsewhere after you posted this, but I thought I'd post my answer here for others to see. The answer is yes, $f$ has to be constant even if $\varepsilon = 0$. Here's the proof (where $\varepsilon = 0$):
Fix $x \in \mathbb{R}^n$. Then by Markov's inequality, we have:
$$\mu^n(\{y \in \mathbb{R}^n : |f(y) - f(x)| \geq s \} \cap B_r(x)) \leq \frac{1}{s} \int_{B_{r}(x)} |f(y) - f(x)| dy$$
Then by setting $s = cd \space r$ (where $c$ and $d$ are arbitrary positive constants) and dividing both sides by $r^n$, we have:
$$\frac{\mu^n(\{y \in \mathbb{R}^n : |f(y) - f(x)| \geq cd \space r \} \cap B_r(x))}{r^n} \leq \frac{1}{cd \space r^{n+1}} \int_{B_{r}(x)} |f(y) - f(x)| dy$$
The right-hand side tend to $0$ as $r \rightarrow 0^+$, therefore so does the left-hand side. In other words:
$$\lim_{r \rightarrow 0^+} \frac{\mu^n(\{y \in \mathbb{R}^n : |f(y) - f(x)| \geq cd \space r \} \cap B_r(x))}{r^n} = 0$$
Now (whenever $d < 1$) we can lower bound the above expression by:
$$\frac{\mu^n(\{y \in \mathbb{R}^n : |f(y) - f(x)| \geq c |y - x| \} \cap [B_r(x) \setminus B_{d \space r}(x)])}{r^n}$$
and note that $\frac{\mu^n(B_{d \space r}(x))}{r^n} = V \space d^n$, where $V := \mu^n(B_1(0))$. Therefore:
$$\limsup_{r \rightarrow 0^+} \frac{\mu^n(\{y \in \mathbb{R}^n : |f(y) - f(x)| \geq c |y - x| \} \cap B_r(x))}{r^n} \leq V \space d^n$$
However this limit holds for all $d \in \mathbb{R}^{>0}$, and so we can actually conclude that:
$$\lim_{r \rightarrow 0^+} \frac{\mu^n(\{y \in \mathbb{R}^n : |f(y) - f(x)| \geq c |y - x| \} \cap B_r(x))}{r^n} = 0$$
$$\lim_{r \rightarrow 0^+} \frac{\mu^n(\{y \in \mathbb{R}^n : \frac{|f(y) - f(x)|}{|y - x|} \geq c \} \cap B_r(x))}{r^n} = 0$$
Therefore $f$ is approximately differentiable at $x$ with $Df_{ap}(x) = 0$ (the definition of approximately differentiable is given here https://encyclopediaofmath.org/wiki/Approximate_differentiability).
Now I want to quote the fact that $f$ having an everywhere zero approximate derivative implies that $f$ is constant. Unfortunately I couldn't find a good source for this, expect in the case of $n=1$ where it follows from theorem 14.18 in The Integrals of Lebesgue, Denjoy, Perron, and Henstock by Russell A. Gordon.
So for arbitrary dimensions, here is my own personal proof for completeness, modified from the $n=1$ case. Feel free to let me know if there are any mistakes.
Assume, seeking a contradiction, that there exists $a, b \in \mathbb{R}^n$ such that $f(a) > f(b)$. Define $g(x) := f(x) + x \cdot \frac{f(a) - f(b)}{2|b - a|} \frac{b-a}{|b-a|}$.
Then first of all $g$ is approximately differentiable with $D_{ap}g(x)(h) = h \cdot \frac{f(a) - f(b)}{2|b - a|} \frac{b-a}{|b-a|}$ for all $x \in \mathbb{R}^n$.
Now without loss of generality, we can assume that $a = 0$ and $b = (1,0,...,0)$ just so I can talk about moving from $a$ to $b$ as moving from left to right. Note that $D_{ap}g(x)$ is positive in every direction which points rightward (by which I mean the direction has a positive dot product with $b$).
In what follows let $x <_1 y$ be defined to mean that the $1^{\text{st}}$ coordinate of $x$ is less than the $1^{\text{st}}$ coordinate of $y$. Same for similar relations like $=_1$.
Let $W \subseteq \mathbb{R}^n$ be the interior and boundary of some right-angled cone with vertex at $b$, and $a$ at the middle of its base.
Now define $A := \{x \in W : g(x) \geq g(a)\}$. Let $\mathcal{F}$ be a family of subsets of $A$, such that $S \in \mathcal{F}$ whenever all the following hold:
- If $x,y \in S$ are distinct, then $x \not =_1 y$
- If $x,y \in S$ with $x <_1 y$, then $\frac{\mu^n(A \cap B_{|x-y|}(y))}{\mu^n(B_{|x-y|}(y))} \geq \frac{1}{2}$
Note that $\mathcal{F}$ is non-empty as $\{a\} \in \mathcal{F}$, and is partially ordered by $\subseteq$ where every non-empty linearly ordered subset has an upper bound (the union). Therefore, by Zorn's lemma, $\mathcal{F}$ has a maximal element which we can denote by $Z$.
Assume, seeking a contradiction, that $Z$ does not contain an element with maximal $1^{\text{st}}$ coordinate. Then there exists $p$ which is a limit point of $Z$ and has supremum $1^{\text{st}}$ coordinate. Let $(z_k)_{k=1}^\infty$ be a sequence in $Z$ which converges to $p$.
Then for each $x \in Z$, we have:
$$\mu^n(A \cap B_{|x-p|}(p))$$
$$= \lim_{k \rightarrow \infty} \mu^n(A \cap B_{|x-z_k|}(z_k))$$
$$\geq \lim_{k \rightarrow \infty} \frac{1}{2} \mu^n(B_{|x-z_k|}(z_k))$$
$$= \frac{1}{2} \mu^n(B_{|x-p|}(p))$$
(The justification for the first and last equalities comes from the fact that $\lim_{k \rightarrow \infty} \mu^n(B_{|x-z_k|}(z_k) \Delta B_{|x-p|}(p)) = 0$, where $\Delta$ denotes the symmetric difference).
So $\frac{\mu^n(A \cap B_{|x-p|}(p))}{\mu^n(B_{|x-p|}(p))} \geq \frac{1}{2}$ for all $x \in Z$. Also note that $g$ is approximately continuous at $p$, so we must have that $p \in A$.
Then $Z \cup \{p\} \in \mathcal{F}$, which contradicts the maximality of $Z$. Therefore $Z$ contains an element with maximal $1^{\text{st}}$ coordinate. Denote the element as $q$.
Now we must have that $q =_1 b$, otherwise we could use the fact that $g$ is approximately differentiable at $q$ (with positive derivative in any rightward direction) to extend $Z$, again contradicting its maximality. But due to $q \in W$, we must have that $q = b$.
So $b = q \in Z \subseteq A$. Hence $g(a) \leq g(b)$, which contradicts the assumption that $f(a) > f(b)$ from the beginning.
Hence $f$ is constant, as $a$ and $b$ are arbitrary.