Dear All,
Okay, so I tried the code that Dr. Blomberg provided, and it doesn't seem
to work. I rewrote the code slightly in order for it to try and predict the
body mass of the taxa used to calculate the line, rather than some new
taxon, because I need to do that in order to calculate prediction error and
be able to compare it to previous studies. I got a lambda of about 0.9,
which is approximately the same as what I got previously. So on the surface
it seems like the method should work.
The problem comes when I try to apply this to the data itself to calculate
mu. This is the modified code I used if that helps, though I don't think I
can send the data or tree to make it replicable because the message
wouldn't go through last time.
tree.trimmed <- drop.tip(trees[[1]], which(is.na(b$bm)))
b.trimmed <- b[complete.cases(b$bm),]
fit <- gls(log(bm)~I(log(ocw)^(2/3)), corPagel(0.5, phy=tree.trimmed,
form=~species),
data=b.trimmed)
(lambda <- as.numeric(fit$modelStruct))
model.matFULL <- model.matrix(~I(log(ocw)^(2/3)), data=b)
preds <- predict(fit,newdata=b)
mat <- vcv.phylo(trees[[1]])
C <- 1+(lambda-1)*(1-mat) ## Lambda correction
XnoMean <- scale(model.matFULL, scale=FALSE)
mu <- C %*% solve(C) %*% XnoMean ## Replaced CCompletedata and cihmat with
C, since both seem to be subsets of that data. This is the only place I
could think things might go wrong.
mu2<-(mu[match(b$species,row.names(mu)),])[,2] ## Have to re-sort mu based
on the data sheet. Mu is reported in phylogenetic order, whereas
predict.gls will always report the data in alphabetical order of the row
names even if a phylogeny is present.
exp(preds+mu2) ## "correct" prediction accounting for bias due to
correlations with other species. Exp() is because data were log-transformed.
So when I check the results of this equation the predicted values are
actually *less *accurate than the fitted values reported from gls without
correction. Normally the accuracy is about +/- 64% of the actual value, but
adding the correction here results in the data being +/- 84% of the actual
value. What's more, I took a look at the values of mu and they don't seem
to make sense when compared on the phylogeny. Namely, if mu were being
calculated right, one would expect monotremes to have a very negative value
compared to all other mammals given the deep divergence in trait states
between monotremes and all other taxa in the figure I previously attached.
But when I check the mu values, I get a value of 0.35 for *Zaglossus*,
-0.58 for *Tachyglossus*, and 0.59 for *Ornithorhynchus*. Which is not what
would be expected if mu were modelling a phylogenetic correction factor (if
nothing else, the range of mu should probably not span zero). So it
suggests to me I'm not doing something right. I know that PGLS doesn't
always result in better fits than OLS, but I'm surprised adding signal back
in made the results worse when the fit was so poor when it was fitted
through the ancestral node.
Sincerely,
Russell
[[alternative HTML version deleted]]
_______________________________________________
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
