I just posted a function that I wrote a while ago (and updated today) for phylogenetic reduced major axis (RMA) regression. The function is available on my R-phylogenetics page (direct link here).
The calculation of the RMA regression slope is straightforward. It is described in Ives et al. (2007). Basically, to get the slope we just compute the phylogenetic variance of y & x and then take the square root of their ratio.
I.e., for phylogeny tree, and vectors x & y in order tree$tip.label:
> ax<-sum(solve(C,x))/sum(solve(C)) # phylogenetic mean for x
> ay<-sum(solve(C,y))/sum(solve(C)) # same for y
To get the intercept, we just plot through the phylogenetic means:
My function also (optionally) optimizes Pagel's lambda & computes the residuals in the y dimension. When I get a chance I will add hypothesis testing of beta0 & beta1 following Ives et al. (2007).