Let $\{e_k\}_{1\le k\le 16}$ denote the standard basis of $\mathbb{R}^{16}$, and  $x:=(x_1,\dots,x_{16})$. Let's consider the difference operator in the $k$-th variable, $\delta_k:\mathbb{R}[x]\to\mathbb{R}[x ]$, that is $\delta_kp(x):=p(x+e_k)-p(x)$. So $$\delta_{13}\delta_9 \delta_5 \delta_1p(x)=\sum_{\epsilon } (-1)^{|\epsilon|_1}p(x+\epsilon_i ),$$
the sum being extended over all $\epsilon\in\{0,1\}^{16}$ with support in the set $S:=\{1,5,9,13\}$: it vanishes if and only if $p$ is of the form $p=\sum_{k\in S} p_k$ for some $p_k\in\mathbb{R}[x]$ with  $\deg_k p_k=0 $,  for any $k\in S$.
Now if $g\in \mathbb{R}[x]$ has $g^{-1}(0)\cap\{0,1\}^{16}=\mathcal{Z},$ we have $\delta_{13}\delta_9 \delta_5 \delta_1g(0)=g(e_1+e_5+e_9+e_{13})\neq0$, proving that $g$ contains a monomial of positive degree in all  variables $x_1, x_5, x_9,$ and $x_{13}$ (and for the same reason, it must also contain  any term of the expansion of $f$).