A little while ago I received the following query:
"To account for body size I used the phylogenetic regression of your phyl.resid function. I used the mean value of each trait per species. However, now I would like to plot the residuals against the ecological variables to show how they relate, including standard deviations for trait means. Is there a way to calculate phylogenetically corrected residuals with standard deviations? Or simply calculate residuals for more than one individual per species?"
Well, there are multiple ways we could go about doing this - with multiple possible solutions. The simplest thing, and the one I'd probably recommend for starters, is to fit the phylogenetic regression to x̄ and ȳ, the species means vectors; but then compute ŷ = [1 x]β and e = y - ŷ from the original data vectors, x & y.
Here's how that would look if we had input data vectors xi and yi, both with the same order but in which names(xi) and names(yi) contain only species names (& thus no unique identifiers):
## now fit PGLS model
## compute yhat based on individual data in xi
## compute individualized residuals
Note that the species-wise mean of ei will be exactly equal to our residuals from the fitted model using species means. E.g., compare: