The formula
$$r(z) \sim \int_0^{N_S(z)} N_T(z-N_S^{-1}(x)) dx$$
can be re-written in a more appealing way. With the change of variable $u=N_S^{-1}(x)$ it becomes
$$r(z) \sim \int_0^{z} N_T(z-u)N'_S(u) du,$$
where $N'_S(u)$ is the derivative of $N_S(u)$ with respect to $u$. With an additional change of variable $u=zv$ it becomes 
$$r(z) \sim z\int_0^{1} N_T(z(1-v))N'_S(zv) dv.$$
Likewise
$$t(z) \sim r'(z) = \frac{dr(z)}{dz} =z\int_0^{1} N'_T(z(1-v))N'_S(zv) dv .$$

An interesting case is when $S=T$ and
$$N_S(u) \sim  \frac{a u^b}{(\log u)^c}, \mbox{ with } 0<a, 0<b<1, \mbox{ and } c \geq 0.$$
This covers sums of two primes ($a=1, b=1, c=1$) and sums of two squares ($a=1, b=\frac{1}{2}, c=0$). We have:
$$r(z) \sim \frac{a^2b z^{2b}}{(\log z)^{2c}}\cdot \int_0^1 (1-v)^b v^{b-1}dv = 
\frac{a^2b z^{2b}}{(\log z)^{2c}}\cdot \frac{\Gamma(b)\Gamma(b+1)}{\Gamma(2b+1)}$$

$$r'(z) \sim \frac{2 a^2 b^2 z^{2b-1}}{(\log z)^{2c}}\cdot \int_0^1 (1-v)^b v^{b-1}dv = 
\frac{2a^2 b^2 z^{2b-1}}{(\log z)^{2c}}\cdot \frac{\Gamma(b)\Gamma(b+1)}{\Gamma(2b+1)}$$

**Notes**

 - Solutions such as $z=x+y$ and $z=y+x$ count as two solutions: $(x,y)$
   and $(y, x)$.
 - The asymptotic formula for $t(z) \sim r'(z)$, representing the number
   of solutions to $z=x+y$ with $x\in S, y\in T$ is true only *on
   average*, as $z$ becomes larger and larger. There may still be
   infinitely many integer $z$'s for which $t(z)=0$ even if $r'(z)
   \rightarrow\infty$ as $z\rightarrow\infty$.
 - We assume that the sets $S$ and $T$ are "well balanced", both for
   small and large values. For instance, if you remove the first
   $10^{5000}$ elements of $S$, the asymptotic formula for $N_S(u)$
   remains unchanged, but this is likely to cause many formulas to
   fail.
 - On some tests, I noticed that there are more solutions (on average)
   to $z=x+y$ with $x\in S, y\in T$ (here $x, y, z$ are integers), if
   $z$ is even.

 - If $S=T$ is the set of primes, some adjustments must be made because
   the primes are not "well balanced", they are less random than they
   seem (for instance the sum of two odd primes can not be an odd
   number, but there are also more subtle issues). This is best
   described in the [Wikipedia entry][1] about Goldbach's conjecture
   (see section about heuristics).

 - To generate a set like $S$, one way is as follows. Use a random
   number generator function $U$ returning independent uniform deviates on $[0,
   1]$. If $U(k) < N'_S(k)$ then add the integer $k$ to the set $S$,
   otherwise discard it. Do that for all integers.

 - For sums involving three terms, say $R+S+T$, you can proceed as follows: first work on $S'=R+S$ and derive all the asymptotics for $S'$ using the methodology proposed here. Then work on $S'+T$.

 - If there are singularities in the functions $N_S$ or $N_S'$, they
   must be handled properly in the integral formulas, unless the
   integrals are improper but converging.

**Generalization of the formula**

It also works if $S\neq T$. Say 

$$N_S(u) \sim  \frac{a_1 u^{b_1}}{(\log u)^{c_1}}, N_T(u) \sim  \frac{a_2 u^{b_2}}{(\log u)^{c_2}}$$
with  $0<a_1,a_2, 0<b_1, b_2 <1$, and  $c_1, c_2 \geq 0$. I am working on getting the formula right, I will publish here when double-checked. 

$$r(z) \sim  
\frac{a_1 a_2 z^{b_1 + b_2}}{(\log z)^{c_1+c_2}}\cdot \frac{\Gamma(b_1 +1)\Gamma(b_2+1)}{\Gamma(b_1 + b_2+1)}$$

$$r'(z) \sim  
\frac{a_1 a_2 z^{b_1 + b_2 -1}}{(\log z)^{c_1+c_2}}\cdot \frac{\Gamma(b_1 +1)\Gamma(b_2+1)}{\Gamma(b_1 + b_2)}$$


In particular, it applies to sums of a square and a prime, see [here][2].


  [1]: https://en.wikipedia.org/wiki/Goldbach%27s_conjecture
  [2]: https://math.stackexchange.com/questions/3710032/conjecture-all-but-21-non-square-integers-are-the-sum-of-a-square-and-a-prime/