How do you derive Lagrange's equation of motion from a Routhian?

No. The coordinate $r$ stil follows the Euler-Lagrange equation, but $\phi$ and $p_\phi$ follow the Hamilton equations. But these are trivial, which is the whole point of the Routhian. The motivation is that the Routhian isn't really $R(r, \dot r, \phi, p_\phi)$ but just $R(r, \dot r)$ with a constant parameter $p_\phi$. $\phi$ isn't a coordinate because by definition, it was a cyclic coordinate in the Lagrangian and so it doesn't appear in the Routhian either. The momentum $p_\phi$, meanwhile, is conserved so it's really just some constant. With both of these conditions, we can simply take them both out, call $p_\phi$ a constant, and we wind up with a problem with 1 fewer effective dimensions. The same would be true using the complete Hamiltonian, but sometimes Lagrangians are easier to work with for the "hard" part of the problem (the non-cyclic coordinates), and the Routhian lets you stay "Lagrangian" for the hard part.

Here is a concrete example. Take the Lagrangian describing a harmonic potential in polar coordinates:

$$ \mathcal L = \frac{1}{2} m \left( \dot r^2 + r^2 \dot \phi^2 \right) - \frac{1}{2} k r^2 $$

Since $\phi$ does not appear in the Lagrangian, it is a cyclic coordinate. From the Euler-Lagrange equation

$$ \frac{d}{dt} \frac{\partial \mathcal L}{\partial \dot \phi} = \frac{\partial \mathcal L}{\partial \phi} = 0 $$

it follows that

$$ p_{\varphi} =\frac{\partial \mathcal L}{\partial \dot \phi} = m r^2 \, \dot \phi $$

is a conserved quantity. It would be nice to take $\dot\phi$ out of the set of coordinates and just replace it with the constant $p_\phi$. Then we'd essentially have a Lagrangian whose only coordinates are $r$ and $\dot r$, just with a constant parameter $p_\phi$. The way to do this correctly is to make the Routhian

$$ \mathcal R = \mathcal L - p_\phi \dot \phi = \frac{1}{2} m \dot r^2 - \frac{1}{2} k r^2- \frac{p_\phi^2}{2 mr^2} $$

(Aside: Try solving $\dot \phi$ in terms of $p_\phi$ and substitute into the Lagrangian. You'll get something very different that will lead to very wrong results. If you want to remove the cyclic coordinate, you must do a Legendre transform w.r.t. that coordinate.).

Since the Legendre transform has been done w.r.t. $\theta$ and $\phi$, these coordinates follow Hamilton's equations (with an appropriate change of sign). But these is trivial:

$$ \dot \phi = - \frac{\partial \mathcal R}{\partial p_\phi} = \frac{p_\phi}{mr^2}, \quad \dot p_\phi = \frac{\partial \mathcal R}{\partial \phi} = 0 $$

What we've achieved here is essentially separation of variables. $\phi$ has its own equation of motion, which involves $r$, but the equation of motion of $r$ is independent of $\phi$. We can solve the equation of motion of $r$, which is an effectively one-dimensional problem, and then go back to find out how $\phi$ evolves.

The effective one-dimensional problem has an effective potential

$$ V(r) = \frac{1}{2} k r^2 + \frac{p_\phi^2}{2mr^2} $$

This extra term, which blows up at small $r$, is called the centrifugal barrier, and accounts for the fact that angular momentum conservation makes it harder / impossible for the particle in question to reach the origin. The coordinates $r$ and $\dot r$ didn't have a Legendre transform applied, so they still follow the Euler-Lagrange equation

$$ \frac{d}{dt} \frac{\partial \mathcal R}{\partial \dot r} = \frac{\partial \mathcal R}{\partial r} \Rightarrow m \ddot r = -kr + \frac{p_\phi^2}{3mr^3} $$

A follow up question might be "Yes, its convenient to not stay with a full Lagrangian, but why not just go full Hamiltonian?" The answer here probably depends on the context of the specific problem. With the Routhian, we have a 2nd order differential equation to solve, and then an integral to find the e.o.m. of $\phi$. With a full Hamiltonian, we'd have a system of two 1st order differential equations, followed by the same integral. One might be easier to solve, or better computationally, etc.