$W^{s,p}(\mathbb{R}^{n})$ Is Not Closed Under Multiplication when $s\leq n/p$
Unless I'm missing something obvious, I believe I have a complete answer that demonstrates the existence of an $f\in W^{s,p}(\mathbb{R}^{n})$, for $0<s\leq n/p$ using a combination of hard analysis and soft analysis.
First, a small lemma:
Lemma 1. If $f\in W^{s,p}(\mathbb{R}^{n})$ is such that the operator defined by pointwise multiplication by $f$ is bounded on $W^{s,p}(\mathbb{R}^{n})$, then $f\in L^{\infty}(\mathbb{R}^{n})$.
Proof. Let $T$ denote the multiplier operator defined by $f$. By hypothesis $T$ is bounded $W^{s,p}\rightarrow W^{s,p}$. For all $g,h\in\mathcal{S}(\mathbb{R}^{n})$, we have that $$\int_{\mathbb{R}^{n}}T(g)\overline{h}=\int_{\mathbb{R}^{n}}fg\overline{h}=\int_{\mathbb{R}^{n}}g\overline{\overline{f}h}$$ By density, we see that the adjoint operator $T^{*}$ densely defined by multiplication by $\overline{f}$ is bounded $W^{-s,p'}\rightarrow W^{-s,p'}$, where $1/p+1/p'=1$. Taking complex conjugates, we see that the operator densely defined by multiplication by $f$ is bounded $W^{-s,p'}\rightarrow W^{-s,p'}$. Interpolating (see T. Tao's Notes 5 here), we see that multiplication by $f$ is bounded on $W^{s_{\theta},p_{\theta}}$, for all $$s_{\theta}:=\theta s-(1-\theta)s \quad p_{\theta}:=\dfrac{\theta}{p}+\dfrac{1-\theta}{p'}, \qquad\forall \enspace 0<\theta<1$$ In particular, for $\theta=1/2$, we obtain that the operator $T$ is bounded on $W^{0,2}=L^{2}$. Testing $T$ on simple functions, we see that $f$ is essentially bounded. $\Box$
Lemma 2. $W^{s,p}(\mathbb{R}^{n})\not\subset L^{\infty}(\mathbb{R}^{n})$, for $0<s\leq n/p$, $1<p<\infty$.
Assuming the lemma, take an $f\in W^{s,p}(\mathbb{R}^{n})\setminus L^{\infty}(\mathbb{R}^{n})$. By Lemma 1, there exists $g\in W^{s,p}(\mathbb{R}^{n})$ such that $fg\notin W^{s,p}(\mathbb{R}^{n})$. Define $h:=f+g\in W^{s,p}(\mathbb{R}^{n})$. If $f^{2}$ or $g^{2}\notin W^{s,p}(\mathbb{R}^{n})$, then we're done. Otherwise by the reverse triangle inequality, $$\|h^{2}\|_{W^{s,p}}\geq 2\|fg\|_{W^{s,p}}-\|f^{2}\|_{W^{s,p}}-\|g^{2}\|_{W^{s,p}}=\infty$$
We now prove Lemma 2. If $W^{s,p}$ embeds (not necessarily continuously) in $L^{\infty}$, then I claim it must do so boundedly. Indeed, Let $I$ denote the inclusion map. By the Closed Graph Theorem, it suffices to show that if $f_{n}\rightarrow f$ in $W^{s,p}$ and $If_{n}\rightarrow g$ in $L^{\infty}$, then $If=g$. Observe that $f_{n}\rightarrow f$ and $f_{n}\rightarrow g$ in distribution and therefore $f=g$ a.e. We arrive at a contradiction, as my analysis in the original post shows that that $W^{s,p}(\mathbb{R}^{n})$ does not continuously embed in $L^{\infty}(\mathbb{R}^{n})$.
What follows is not a solution to the original problem, but rather some length comments, which I thought best to post separately, so as to minimize the length of the original post.
First, it's not hard to see that for general $0<s\leq n/p$, it is insufficient to work at the $L^{p}$ level. Indeed, if $q$ is the Sobolev conjugate given by $1/q=1/p-s/n$ and $1/q\geq 1/2p$, then by Sobolev embedding we have that $\|f\|_{L^{2p}}\lesssim_{n,s,p}\|f\|_{W^{s,p}}$. In particular, this is true for the endpoint case $s=n/p$.
Second, it is sufficient to find an $f\in W^{s,p}(\mathbb{R}^{n})$ such that $$\|\sup_{k}|2^{ks}P_{k}f^{2}|\|_{L^{p}}=\infty$$ For those of you familiar with Besov spaces, this condition is equivalent to $\|f^{2}\|_{B_{p,\infty}^{s}}=\infty$. I mention this condition because I imagine that the $f\mapsto f^{2}$ is sharply bounded (i.e. the regularity parameter of the target space is sharp) from $W^{s,p}(\mathbb{R}^{n})\rightarrow W^{s-\epsilon,p}(\mathbb{R}^{n})$, for some $\epsilon>0$, and therefore we have the inclusions $$W^{s-\epsilon,p}(\mathbb{R}^{n})\subset B_{p,\infty}^{s-\delta}(\mathbb{R}^{n})\subset W^{s,p}(\mathbb{R}^{n}),\quad\forall \enspace 0<\delta<\epsilon$$
Third, when $0<s<n/p$, it is not too difficult to show that the nonlinear operator $f\mapsto f^{2}$ is unbounded on $W^{s,p}(\mathbb{R}^{n})$. Indeed, let $\psi$ be a nonnegative bump function adapted to the annulus $1\leq|\xi|\leq 2$ and consider $f_{j}=2^{-sj+nj/p}\widehat{\psi}(2^{j}\cdot)$, for positive integer $j$. Our analysis above shows that $\|f_{j}\|_{W^{s,p}}\lesssim_{n,s,p}1$. By dilation invariance, \begin{align*} \|2^{sk}P_{k}f_{k}^{2}\|_{L^{p}}^{p}&=2^{skp}\int_{\mathbb{R}^{n}}\left|\int_{\mathbb{R}^{n}}2^{kn}\varphi(2^{k}(x-y))2^{-2sk+2nk/p}\widehat{\psi}(2^{k}y)^{2}dy\right|^{p}dx\\ &=2^{-skp+2nk}\int_{\mathbb{R}^{n}}\left|\int_{\mathbb{R}^{n}}\varphi(2^{k}x-y)\widehat{\psi}(y)^{2}dy\right|^{p}dx\\ &=2^{-skp+nk}\int_{\mathbb{R}^{n}}\left|\int_{\mathbb{R}^{n}}\varphi(x-y)\widehat{\psi}(y)^{2}dy\right|^{p}dx \end{align*} So $\|f_{k}^{2}\|_{W^{s,p}}\gtrsim 2^{k(\frac{n}{p}-s)}\rightarrow\infty$, as $k\rightarrow\infty$. Maybe there's a soft analysis argument from which one can obtain the existence of the desired function $f$, but since the square operator is nonlinear, it's not clear to me what this argument would be.