I’m looking for a method to account for and reduce phylogenetic autocorrelation in a linear model.

My approach has been to reduce residuals based on Moran's eigenvectors, and then rerun the ln model conditioned on these results. To do this, I convert my phylogeny into a phylogenetic distance matrix, and then to a listw object using mat2listw() in R's spdep. Then, choose a model that reduces the residuals using mem.select() from the adespatial package, using Moran's I (MIR). This returns a MEM matrix but will not select a model, and does not return an error message (but there is no MEM$MEM.select object). I’m not sure why I can’t seem to select a model, or that adapting a package designed for spatial autocorrelation is the best way to go about this.

Is there a preferred method for accounting for phylogenetic autocorrelation in trait data?