In Hörmander, Analysis of Linear Partial Differential Operators (page 85), it states that with $x_+^n=x^n$, for $x>0$ and $x_+^n=0$, for $x\leq 0$,
$x_+^{a-1} * x_+^{b-1}=B(a,b)x_+^{a+b-1}$,
where $B(a,b)=\int_{0}^1 t^{a-1}(1-t)^{b-1} \,dt$ and $a>0$,$b>0$.
When trying to solve this directly, I get $x_+^{a-1} * x_+^{b-1}=\int_{\mathbb{R}_{\geq 0}} y^{a-1}(x-y)^{b-1} \,dy$. I know I should be pulling a term of $y^{a+b-2}$ out an then integrating it, but I don't see how to proceed here.
There is an error in your upper integration bound which should be:
$$\int_{y=0}^{y=\color{red}{x}}$$
Nevertheless, I think that, as it is often the case with convolution, coming back to its definition as an integral is not the best solution. It is better to use some of its properties.
Here, as we have functions that are $0$ for $x<0$, it is natural to think to Laplace Transform (L.T.). Here is how we can use it.
Relationship:
$$x_+^{a-1} * x_+^{b-1}=B(a,b)x_+^{a+b-1}\tag{1}$$
will be established if we prove that the LHS and RHS of (1) have the same L.T.
Using the fact that L.T. transforms a convolution into a product and that
$$t_+^{a-1} \xrightarrow{L.T.} \tfrac{\Gamma(a)}{s^{a}}$$
(valid for any real $a>0$), we see that indeed the LHS and RHS of (1) have the same L.T.:
$$\frac{\Gamma(a)}{s^{a}}\frac{\Gamma(b)}{s^{b}}=\frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}\frac{\Gamma(a+b)}{s^{a+b}}$$
Remark: another very different proof would use the property $$(f \star g)' = f' \star g = f \star g'$$
where the "primes" denote derivation.