The hypothesis that $F$ is closed is not necessary. However, note that the distance from a set is the same as the distance from its closure. By passing from $F$ to its closure we lose no generality, and in fact make a stronger statement because the estimate is claimed to hold a.e. on a larger set.
One proof I like: the function $\delta$ is Lipschitz, therefore differentiable a.e. (it has bounded variation on finite intervals). In particular, it is differentiable at a.e. point of $F$. But every point of $F$ is a point of minimum for $\delta$. So the derivative of $\delta$ is $0$ whenever it exists.
If you want to argue in terms of density, I think it's easier to do the contrapositive. Namely, assume that there is a sequence $y_k\to 0$ such that $\delta(x+y_k)\ge \epsilon |y_k|$, with fixed $\epsilon>0$, and show that $F$ does not have density $1$ at $x$. Reason: there's a hole of size $2\epsilon |y_k|$ at $x+y_k$. So, the intersection of $F$ with $[x-2|y_k|,x+2|y_k|]$ has measure at most $2|y_k|(1-\epsilon)$, which is not enough to have density $1$.
First, by the Lebesgue density theorem, almost every point $x\in A$ has density
$$
\begin{align}
d(x)=\lim_{\epsilon\to0}\frac{\mu(B_\epsilon(x)\cap A)}{B_\epsilon(x)}&&{\rm(1)}
\end{align}
$$
equal to 1. Similarly, almost every point in $\mathbb{R}\setminus A$ has Lebesgue density $d(x)=0$. Every other point is known as an an exceptional point. By what has just been said, the set of exceptional points of a measurable $A\subseteq\mathbb{R}$ has measure $0$.
On the other hand, it is known that for any measurable $A$ such that both $A$ and $\mathbb{R}\setminus A$ have positive measure then there is always at least one exceptional point. By definition, at such points, one of the following holds.
- The limit (1) exists and is strictly between $0$ and $1$.
- The limit (1) does not exist.
While it is well known that such exceptional points do exist, the question asks whether there always exists exceptional points of case 1. That is, the density exists but is strictly between 0 and 1. For this, the answer is no. There are such sets $A$ for which the Lebesgue density does not exist at any exceptional point.
First, the follwoing argument shows that exceptional points exist. If $x$ is a density point of $A$ and $y$ is a density point of $\mathbb{R}\setminus A$ then, assuming (wlog) that $x < y$ we define $f\colon[x,y]\to\mathbb{R}$ by $f(t)=\mu((x,t)\cap A)-t/2$. We have $f^\prime(x)=1/2$, $f^\prime(y)=-1/2$. So, the maximum of $f$ occurs in the interior $(x,y)$. If $t^*$ maximises $f(t^*)$ then,
$$
\frac{\mu(B_\epsilon(t^*))}{\mu(B_\epsilon(x))}=\frac{f(t^*+\epsilon)-f(t^*-\epsilon)}{2\epsilon}+\frac12=\frac{f(t^*)-f(t^*-\epsilon)}{2\epsilon}-\frac{f(t^*)-f(t^*+\epsilon)}{2\epsilon}+\frac12
$$
As long as $0\le t^*-\epsilon \le t^*+\epsilon\le1$, each of the first two fractions on the right hand side is nonnegative and bounded by 1/4, showing that the Lebesgue density, if it exists, is in the range $[1/4,3/4]$, so $t^*$ is an exceptional point.
However, the question asks whether there are exceptional points where the density $d(x)$ exists (and is, by definition, strictly between $0$ and $1$). The answer to this is no, not necessarily. I will now construct a set $A$ for which $d(x)$ does not exist at any exceptional point.
I'll follow the construction given in the thesis by Jack Grahl, originally due to Szenes and Kurka. Let $U$ be a finite union of open intervals in the unit interval $[0,1]$ and $C=(-\infty,0)\cup U$. For each of the finite set of points $x\in\partial C=\bar C\setminus C$ suppose that there exists positive reals $r_1(x),r_2(x)$ with $\mu(B_{r_1(x)}\cap C)/\mu(B_{r_1(x)})\not=\mu(B_{r_2(x)}\cap C)/\mu(B_{r_2(x)})$. For example this is true for $U=(u,v)$ with $0 < u < v\le1$ and $u\not=v/2$. Then choose $\delta > 0$ with $\mu(B_{r_1(x)}\cap C)/\mu(B_{r_1(x)})-\mu(B_{r_2(x)}\cap C)/\mu(B_{r_2(x)}) > \delta$ for all $x\in\partial C$.
For any open $V$ I'll write $\partial_L V$ (resp. $\partial_R V$) for the set of limit points of $V$ which are left (resp. right) end points of intervals in $V$. For a fixed $\epsilon > 0$ set
$$
\begin{align}
&A_0=(0,1),\\
&A_{n+1}=A_n\cup\bigcup_{x\in\partial_L A_n}(x-\epsilon^{n+1}U)\cup\bigcup_{x\in\partial_R A_n}(x+\epsilon^{n+1}U),\\
&A=\bigcup_{n=0}^\infty A_n.
\end{align}
$$
For $\epsilon$ small enough, the Lebesgue density will not exist at any point of $\partial A$. I give a proof of this below, follow the arguments of 1. It is a little messy though...
Any $x\in\partial A$ is within a distance $\epsilon^{n+1}+\epsilon^{n+2}+\cdots=\epsilon^{n+1}/(1-\epsilon)\le 2\epsilon^{n+1}$ of a point $y$ of $\partial A_{n}$ (assuming $\epsilon\le1/2$). It follows that
$$
\begin{align}
\left\lvert\mu(B_r(x)\cap A)-\mu(B_r(y)\cap A)\right\rvert\le 2\epsilon^{n+1}.&&{\rm(2)}
\end{align}
$$
Next, letting $L$ be the minimum length of the connected components of $U$, the connected components of $A_{n-1}$ have length at least $\epsilon^{n-1}L$, so none are strictly contained in the intervals $(y-\epsilon,y)$ or $(y,y+\epsilon)$ so long as
$$
\begin{align}
r\le\epsilon^{n-1}L.&&{\rm(3)}
\end{align}
$$
Letting $M$ be the minimum length of a connected component of $\mathbb{R}\setminus C$ then the minimum length $M_n$ of connected components of $\mathbb{R}\setminus A_n$ satisfies $M_1=\epsilon M$ and $M_{n+1}\ge\min(\epsilon^{n+1}M,M_n-2\epsilon^{n+1})$. So long as
$$
\begin{align}
(M+2)\epsilon\le M&&{\rm(4)}
\end{align}
$$
this gives $M_{n+1}\ge\epsilon^{n+1}M$. We assume that $\epsilon$ is small enough for (4) to hold.
Now, as $y$ is a limit point of $A_n$, there is a $z\in\partial_R A_{n-1}$ (resp. $\partial_L A_{n-1}$) such that $y=z$ or $y\in z+\epsilon^n\partial U$ (resp. $y\in z-\epsilon^n\partial U$). Assume, wlog, that $z\in\partial_R A_{n-1}$. Then, $y$ is in the boundary of $z+\epsilon^nC$. If (3) holds then
$$
B_r(y)\cap(-\infty,z)=B_r(y)\cap(-\infty,z)\cap A_{n-1}=B_r(y)\cap(-\infty,z)\cap A_{n}.
$$
Next, if
$$
\begin{align}
r+\epsilon^n\le\epsilon^{n-1}M&&{\rm(5)}
\end{align}
$$
(so $r+\epsilon^n\le M_n$) then $(y,y+r+\epsilon^n)$ does not strictly contain any interval of $\mathbb{R}\setminus A_{n-1}$. This means that $B_r(y)$ does not intersect $z^\prime-\epsilon^n U$ for any $z^\prime > z$ in $A_{n-1}$. Hence,
$$
B_r(y)\cap[z,\infty)\cap A_n=B_r(y)\cap(z+\epsilon^nU).
$$
So, if (3) and (5) are satisfied, we have
$$
\begin{align}
&B_r(y)\cap A_n= B_r(y)\cap(z+\epsilon^n C),\\
&\frac{\mu(B_r(y)\cap A_n)}{\mu(B_r(y)}=\frac{\mu(B_{\epsilon^{-n}r}(y^\prime)\cap C)}{\mu(B_{\epsilon^{-n}r}(y^\prime)}&&{\rm(6)}
\end{align}
$$
where $y^\prime=\epsilon^{-n}y-z$. Letting $r^*$ be the maximum of $r_1(x)$ and $r_2(x)$ over $\partial C$ we can plug in $r=\epsilon^n r_1(y^\prime)$ and $r=\epsilon^n r_2(y^\prime)$ respectively and (2),(5) will be satisfied as long as $\epsilon$ is small enough that
$$
\begin{align}
&\epsilon r^*\le L\\
&\epsilon(r^*+1)\le M.
\end{align}
$$
Also, in this case, $(A_{m+1}\setminus A_m)\cap B_r(y)$ contains at most $N^{m-n+1}$ copies of $U$ scaled by $\epsilon^{m+1}$ (where $N$ is the size of $\partial C$). So, for small $\epsilon$,
$$
\mu(B_r(y)\cap A)-\mu(B_r(y)\cap A_n)\le\sum_{m=n}^\infty N^{m-n+1}\epsilon^{m+1}=N\epsilon^{n+1}/(1-\epsilon)\le 2N\epsilon^{n+1}.
$$
by at most $2\epsilon^{n+1}$ and this differs from $\mu(B_r(y)\cap A_{n})$ by at most $(3\epsilon)^{n+1}$. So,
$$
\left\lvert\frac{\mu(B_r(x)\cap A)}{\mu(B_r(x))}-\frac{\mu(B_r(y)\cap A_{n})}{\mu(B_r(y))}\right\rvert\le r^{-1}(4\epsilon)^{n+1}.
$$
Combining this with (2) and (6) gives
$$
\left\lvert\frac{\mu(B_{s_1}(x)\cap A)}{\mu(B_{s_1}(x)}-\frac{\mu(B_{s_2}(x)\cap A)}{\mu(B_{s_2}(x)}\right\rvert\ge\delta-\frac{2\epsilon}{\tilde r}-\frac{2N\epsilon}{\tilde r}
$$
where $s_1,s_2\le \epsilon^n r^*$ and $\tilde r$ is the minimum of $r_1(x)$ and $r_2(x)$ over $\partial C$. Choosing $\epsilon$ small enough that the right hand side is strictly positive, the left hand side does not vanish as $s_1,s_2\to0$, and the Lebesgue density does not exist at any point of $\partial A$.
Best Answer
I think we only need to assume that almost every $x \in E$ is a point of Lebesgue density of $E,$ which does not require that $E$ is a Lebesgue measurable set. [Indeed, if we can assume that $E$ is a Lebesgue measurable set, then we can conclude much more than what you said, namely that at almost every $x \in E^c$ the Lebesgue density of $E$ at $x$ is equal to $0$ (instead of "is less than $1$"). But I don't see that we need to know anything about the Lebesgue density of $E$ at points not in $E.$]
The result you want follows from the following more precise result.
Theorem: Let $E \subseteq {\mathbb R}^N$, $x \in {\mathbb R}^N$, and assume that $x$ is a point of Lebesgue density of $E.$ Then $$\lim_{|y| \rightarrow 0}\frac{\delta(x+y)}{|y|} = 0$$ Proof: If not, then $\;\limsup_{|y| \rightarrow 0}\frac{\delta(x+y)}{|y|} > 0,\;$ from which it follows that there exists $\epsilon > 0$ and a sequence $\{y_n\}$ with $y_n \rightarrow x$ such that for each $n$ we have $\delta(x+y_n) \geq \epsilon |y_n|.$ The key idea is to notice that this last inequality tells us there is a sequence of arbitrarily small balls centered at $x,$ each of which contains a "sufficiently big" subball having empty intersection with $E,$ where "sufficiently big" is big enough to contradict the assumption that the Lebesgue density of $E$ at $x$ is equal to $1.$ Specifically, the balls centered at $x$ have radii $|y_n| + \epsilon |y_n|$ and the corresponding subballs are centered at $x+y_n$ and have radii $\epsilon |y_n|.$ Note that the ratio of each subball radius to the corresponding ball radius is ${\epsilon}/({1 + \epsilon}),$ and so the ratio of the volume of each subball to the corresponding volume of the ball is $\left(\frac{\epsilon}{1 + \epsilon}\right)^N.$ Because this ratio of volumes is a positive constant (that depends only on $F,$ $x,$ and $N),$ and no points of $E$ belong to any of the subballs, we cannot have the Lebesgue density of $E$ at $x$ equal to $1,$ which gives a contradiction.
To Summarize: From the assumption "the Lebesgue density of $E$ at $x$ equals $1,$" it follows that, as we shrink down to $x$ using balls centered at $x,$ the ratios of the measures of $E^c$ intersect the balls to the measures of the balls must approach $0$ (because the ratios of the measures of $E$ intersect the balls to the measures of the balls must approach $1$). However, $\delta(x+y_n) \geq \epsilon |y_n|$ provides us with a sequence of balls shrinking down to $x$ in which the ratios of the measures of $E^c$ intersect the balls to the measures of the balls is bounded above $0.$ [Note that the measure of $E^c$ intersect any of the balls is greater than or equal to the measure of the corresponding subball.]
(ADDED NEXT DAY) Something I overlooked: If we don't assume that $E$ is Lebesgue measurable then, in various places above, "measure" and "density" need to be replaced with "outer measure" and "outer density".