You can use the Weierstrass preparation theorem and recursively apply the 1-D result. The theorem asserts that an analytic function $f(x,\mathbf{y})$, for which $f(0,0)=0$, can be written as a product $f(x,\mathbf{y}) = W(x,\mathbf{y}) g(x,\mathbf{y})$, where $W(0,0) \ne 0$ and
$$
g(x,\mathbf{y}) = x^k + g_{k-1}(\mathbf{y}) x^{k-1} + \cdots + g_0(\mathbf{y})
$$
is a Weierstrass polynomial of some degree $k$, with analytic coefficients such that $g_i(\mathbf{y}=0) = 0$.
By continuity, there is some $\delta > 0$ such that $0 < C_1 \le W(x,\mathbf{y}) \le C_2$ on $(x,\mathbf{y}) \in B_\delta$, so
$$
C_1 g(x,\mathbf{y}) \le f(x,\mathbf{y}) \le C_2 g(x,\mathbf{y})
$$
reduces the bounds to those on $g(x,\mathbf{y})$. If we have bounds on the coefficients $b_i(\mathbf{y}) \le g_i(\mathbf{y}) \le c_i(\mathbf{y})$ for $|\mathbf{y}| \le \delta$, then
$$
x^k + b_{k-1}(\mathbf{y}) x^{k-1} + \cdots + b_0(\mathbf{y})
\le g(x,\mathbf{y}) \le
x^k + c_{k-1}(\mathbf{y}) x^{k-1} + \cdots + c_0(\mathbf{y}) .
$$
To get the coefficient bounds, you can now apply the above steps recursively, possibly shrinking $\delta$ as necessary. Once you get to the 1-D case, you can terminate by using the bounds given by monomials. The result will be of the form
$$
b(x,\mathbf{y}) \le f(x,\mathbf{y}) \le c(x,\mathbf{y}) ,
$$
where $b$ and $c$ are polynomials. If your original function only had positive Taylor coefficients, I suspect that all the intermediate upper and lower bounds should also have positive coefficients, but I'm not sure. Perhaps some other trick could be used if a negative lower bound appears for one of the coefficients at an intermediate step.
The lower bound $b$ is in general not going to be a monomial, but if it does have only positive coefficients, you could drop from it as many terms as you would like. Applying this idea to your example
\begin{align*}
f(x,y) &= x^4 + 2 (y^2 + 1) x^2 + y^4 \\
&= W(x,y) \left(x^2 + \frac{1}{2} + 2y^2 - \sqrt{y^2 + \frac{1}{4}}\right) \\
&= W(x,y) (x^2 + y^4 - 2 y^6 + \cdots)
\end{align*}
gives the bounds
$$
B (x^2 + y^4) \le B_2 (x^2 + B_0 y^4) \le f(x,y) \le C_2 (x^2 + C_0 y^4) \le C (x^2 + y^4) .
$$
The outer bounds follow from the inner ones by setting for instance $B = \min(B_2, B_2 B_0)$ and $C = \max(C_2, C_2 C_0)$. Switching around the roles of $x$ and $y$, $f(x,y) = y^4 + 2 x^2 y^2 + (x^4 + x^2)$ is already in Weierstrass form and the resulting bounds are
$$
B (y^2 + x^2)^2 \le f(x,y) \le C (y^4 + x^2 y^2 + x^2) .
$$
The middle coefficient in the upper bound could also be any other positive constant.