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:
> C<-vcv.phylo(tree)
> ax<-sum(solve(C,x))/sum(solve(C)) # phylogenetic mean for x
> ay<-sum(solve(C,y))/sum(solve(C)) # same for y
> beta1<-sqrt(t(y-ay)%*%solve(C,y-ay)/(t(x-ax)%*%solve(C,x-ax)))
To get the intercept, we just plot through the phylogenetic means:
> beta0<-ay-beta1*ax
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).
See the following warning by Simon Blomberg:
ReplyDelete"I think it is important to point out, that while RMA may superficially be an attractive method, it relies on the ratio of error variances being unity. This is almost always incorrect. It usually results in a massive over-correction of the slope bias with respect to the OLS estimator. That is, the slope is made much too steep. I would not encourage anyone to use RMA for anything other than in the case where there is sufficient within-species replication to estimate the error variances with some precision, and then use an appropriate generalization of RMA that allows for the variance ratio to be other than unity. Fiddling around with 'phylogenetically-informed' RMA is like rearranging the deck chairs on the Titanic." (Great metaphor besides.)
Full message here.