How to estimate the order of this integral with parameter Some introduction:
Given a homogeneous structure called "dilation" in $R^n$: For $t\geq 0$
$$D_t: R^n\rightarrow R^n$$
$$D_t(x)=(t^{a_1}x_1,...,t^{a_n}x_n)$$
where $1=a_1\leq...\leq a_n$, and $a_i$ are all integers. And we call $Q=a_1+...+a_n$ the homogeneous dimension. In our problem, we only consider when $Q>n\geq  2$.
Now consider the integral:
$$J(r)=\int_{[0,1]^n}\frac{dx}{P(x,r)}=\int_{[0,1]^n}\frac{dx}{f_n(x)r^n+f_{n+1}(x)r^{n+1}+...+f_Q(x)r^Q}$$
where $f_k(x)$ satisfies:
(1) $f_k(D_t(x))=t^{Q-k}f_k(x)$ for all $x\in R^n$ and $t\geq0$
(2) $f_k(x)$ is the combination of some positive monomials. (Examples will be shown below)
(3) $f_Q(x)=Constant>0$. (This property follows from other theorems and propositions, but they are too many so I don't describe them here.)
Four examples are the followings:
(ex1)  In $R^2$, $D_t(x)=(tx_1,t^2x_2)$, so $Q=3$. And Let $P(x,r)=x_1r^2+r^3$.
(ex2) In $R^3$,  $D_t(x)=(tx_1,tx_2,t^2x)$, so $Q=4$. Let $P(x,r)=(x_1+x_2)r^3+r^4$
(ex3) In $R^3$, $D_t(x)=(t^{1}x_1,t^2x_2,t^{3}x_3)$, so $Q=6$. Let $P(x,r)= x_1^3r^3+(x_2+3x_1^2)r^4+5x_1r^5+3r^6$
(ex4)  In $R^3$, $D_t(x)=(t^{1}x_1,t^2x_2,t^{3}x_3)$, so $Q=6$. Let $P(x,r)= x_1x_2r^3+(x_2+2x_1^2)r^4+3x_1r^5+r^6$
(You will find that $x_n$ doesn't make effort. In my work $x_n$ do make no sense in the integral but this  follows from other theorems, and it doesn't matter here.  )
Problem: Find the order of $J(r)$ when $r$ goes to $0^+$. Like the following description.
Attempt and information: I guess $J(r)=\frac{1}{r^\alpha}I(r)$, where the $\alpha$ is the "critical value", that is:
(i) $\liminf_\limits{r\rightarrow0^+}I(r)>0$.
(ii) for any $\epsilon>0$, $\lim_\limits{x\rightarrow0^+}r^\epsilon I(r)=0$.
I will give the reason why I guess so in the below. I can show that $g_p(r)=r^p J(r)$, then there exists $p_0$ s.t. when $a<p_0$, $\lim_\limits{r\rightarrow0^+}g_a(r)>0$ and when $a>p_0$, $\lim_\limits{r\rightarrow0^+}g_a(r)=0$. But I can't show $\lim_\limits{r\rightarrow0^+}g_{p_0}(r)>0$, that is, I can't show the (i) above. (see https://math.stackexchange.com/questions/3769564/how-to-find-the-critical-index-a-of-xafx) One gave a counterexample for the proposition in that link. But its counterexample will not appear in this problem. Because this is a rational fractional integral. The $I(r)$ I guess will be like the combination of $\log$ and $\arctan$.
The four example have the order estimates:
(ex1) We can calculate directly:
$$J(r)=\frac{1}{r^2}\ln(1+\frac{1}{r})=\frac{1}{r^2}I(r)$$
where $ I(r)$ satisfies (i)(ii) above.
(ex2)
$$J(r)=\frac{1}{r^3}I(r)$$
where $I(r)$ can be calculate or one can use Dominate convergence theorem to estimate that $I(r)$ satisfies (i)(ii)
(ex3)
$$J(r)=\frac{1}{r^{3+2/3}}I(r)$$
see https://math.stackexchange.com/questions/3718932/estimate-a-integral-with-parameter
(ex4)
$$J(r)=\frac{1}{r^{3}}I(r)$$
First
$$J(r)=\frac{1}{r^3}\int_{[0,1]^2}\frac{dxdy}{xy+(y+2x^2)r+3xr^2+r^3}=\frac{1}{r^3}I(r)$$
we can show $I(r)$ satisfies (i)(ii):
(i) change variables:
$$I(r)=\int_{0}^{1/r^2}\int_{0}^{1/r}\frac{dxdy}{xy+(y+2x^2)+3x+1}$$
and then obviously.
(ii) for $3>\epsilon>0$ (the part $\epsilon\geq 3$ follows from the part $3>\epsilon>0$),
$$r^\epsilon I(r)=\int_{[0,1]^2}\frac{r^\epsilon}{xy+(y+2x^2)r+3xr^2+r^3}dxdy=\int_{[0,1]^2}h_r(x,y)dxdy=\int_{(0,1)^2}h_r(x,y)dxdy$$
Pointwisely $\lim_\limits{r\rightarrow0^+}h_r(x,y)=0$ in $(0,1)^2$. Now look for a dominating function in $(0,1)^2$:
$$\frac{1}{h_r(x,y)}\geq \frac{xy}{r^\epsilon}+r^{3-\epsilon}\geq C(xy)^{1-\frac{\epsilon}{3}}$$
So $h_r(x,y)\leq \frac{C}{(xy)^{1-\frac{\epsilon}{3}}}$ in $(0,1)^2$, which is integrable. By DCT, we have $I(r)$ satisfying (i)(ii). But this method doesn't work in other examples like (ex3).
Based on the four examples, I guess $$J(r)=\frac{1}{r^\alpha}I(r).$$ But I can't show how to find the critical value $\alpha$ and even it's difficult to show the existence of critical value
 A: It looks like you care only about the order of magnitude (i.e., an answer up to a constant factor), in which case it is fairly easy.
First, ignore all coefficients. Setting them to $1$ just changes the answer at most constant number of times. Now, suppose we have the denominator of the form $\sum_{(\alpha,\beta)} x^\alpha r^\beta$ where $\alpha$ is a multi-index with real entries and $\beta$ is a real number. The sum is assumed to be finite. Make the change of variable $x_j=e^{-y_j}$. Now, at each point, only the maximal term matters (up to a factor that is the total number of terms). In terms of $y$'s, the condition of maximality of $x^\alpha r^\beta$ is $y_j\ge 0$,
$$
\langle y,\alpha-\alpha'\rangle\le (\beta'-\beta)\log(1/r)
$$
for all $(\alpha',\beta')\ne(\alpha,\beta)$. This domain is just a fixed polyhedron $P_{\alpha,\beta}$ stretched $\log(1/r)$ times (we keep only those with non-empty interiors in what follows; also I call it a "polyhedron" though, technically, it can be unbounded). Thus,
$$
J(r)\asymp\sum_{(\alpha,\beta)}r^{-\beta}\int_{(log\frac 1r)P_{\alpha,\beta}}e^{\psi_{\alpha,\beta}(y)}\,dy
$$
where $\psi_{\alpha,\beta}(y)=\langle \alpha-e,y\rangle$, $e=(1,\dots,1)$.
Now the life becomes straightforward. All you need is to find the order of magnitude of each integral. I'll drop the indices $\alpha,\beta$ for brevity. Let $F$ be the face of $P$ on which $\psi$ attains its maximum $p$ and let $d$ be the dimension of $F$. If $\psi\equiv 0$ (i.e., $\alpha=e$), we just have $F=P$ and $\int_{(\log\frac 1r)P}e^{\psi}=V(P)\log^d(1/r)$. Consider now the non-trivial situation when $\psi$ is not $0$. Then we can rotate and shrink the coordinate system so that $-\psi(y)$ becomes a new variable $t$. Also we can shift $P$ along this coordinate so that the face $F$ lies on the corresponding coordinate hyperplane $\{t=0\}$. Then the integral in question is just
$$
e^{p\log(1/r)}(\log^{D-1}\frac 1r)\int_{0}^\infty e^{-t}S_P(\frac t{\log{1/r}})\,dt
$$
where $S_P(\tau)$ is the $D-1$-dimensional volume of the cross-section of $P$ by the hyperplane $\{t=\tau\}$. By the general convex geometry nonsense, for small $\tau$, $S_P(\tau)=v_d\tau^{D-1-d}+v_{d-1}\tau^{D-d}+\dots+v_0\tau^{D-1}$ where $v_d>0$ and then it becomes smaller (look up "mixed volumes" on Google if you are interested in the details), whence the leading term in the integral becomes $\log^d\frac 1r$ with some coefficient depending on $P$. Thus, the final answer for the integral we are interested in with the factor $r^{-\beta}$ is
$$
\asymp r^{-p_{\alpha,\beta}-\beta}\log^{d_{\alpha,\beta}}\frac 1r
$$
We have several competing terms like that, so the winning one is the one with largest $p+\beta$ and among those the one with the largest $d$.
In your last example $x_1x_2+x_1^2r+x_2r+x_1r^2+r^3$ (I ignore $r^3$ that can be carried out and all the coefficients), we have $5$ polyhedra and functionals (I drop the trivial restrictions $y_1,y_2\ge 0$):
$$
P_{1,1,0}=\{-y_1+y_2\le 1, y_1\le 1, y_2\le 2, y_1+y_2\le 3\},
\\
\psi_{1,1,0}(y)=0 
\\
P_{2,0,1}=\{y_1-y_2\le -1, 2y_1-y_2\le 0, y_1\le 1,2y_1\le 2\}, 
\\
\psi_{2,0,1}(y)=y_1-y_2
\\
et\ cetera.
$$
Here $P_{1,1,0}$ dominates and yields $\log^2\frac 1r$ but it may be instructive to find the contribution of $P_{2,0,1}$. In this case (just draw the picture) $p=-1$, $\beta=1$, $d=1$, so we get $\log\frac 1r$.
