Let's prove that $a_3=a_1$. Note that the recurrence for $R$ translates to the following PDE for the generating $F(x,y):=\sum_{n,m\geq0} R(n,m,p,q) \frac{x^n}{n!}\frac{y^m}{m!}$: $$\frac{\partial}{\partial x} F(x,y) = \frac{\partial}{\partial y} \bigg(F(x,y) + \frac{p}q(e^{qy} - 1)F(x,y)\bigg)$$ with the boundary condition $F(0,y) = e^y$. This PDE is well solvable in CAS like Maple, it can be easily verified that $F(x,0)=\sum_{n\geq0} R(n,0,p,q) \frac{x^n}{n!}$ does indeed coincide with the derivative of the e.g.f for $a_1$ given in Ira's answer. QED
Equality $a_1 = a_4 = a_5$ follows from the Salas and Sokal paper, see formula (3.10) and Proposition 3.3 for $(w,y,u,v)=(1,p,q,p)$.
ADDED. Here are the requested details of computation. It turns out that while Maple solves our PDE, it is not so good at simplifying symbolic radicals (or I'm not familiar with the best practice here) in the solution. So, I will use a combination of Maple and Sage for solving the PDE and simplifying the result, respectively:
sage: S = maple('simplify(subs(y=0, rhs( pdsolve( [ diff(F(x,y),x)=diff(F(x,y) + p/q * (exp(q*y)-1)*F(x,y), y), F(0,y) = exp(y) ], F(x,y) ) ) ))').sage()
sage: p,q,x = S.variables()
sage: ascii_art(S.simplify_real().canonicalize_radical().simplify_real())
-1
q + 1 ---
----- q
q / q*x p*x\ x*(p + 1)
(-p + q) *\- p*e + q*e / *e
-----------------------------------------------
q p - 2*q
----- + -------
q*x p - q p - q p*x
- p*e + q *e
sage: ascii_art( diff(((p-q)/(p-q*exp(x*(p-q))))^(1/q),x).canonicalize_radical() )
-1
q + 1 ---
----- q
q / q*x p*x\ x*(p + 1)
(-p + q) *\- p*e + q*e / *e
-----------------------------------------------
q*x p*x
- p*e + q*e
Sage just left us with noticing that $\frac{q}{p-q}+\frac{p-2q}{p-q}=1$ to conclude equality of the two expressions.