Hi John & list.
Can I add to this by suggesting (in case it is not clear from Joe's
comment) that we should probably combine the variance among estimates
obtained from different trees in the posterior sample with the sampling
variance of any single estimate?
One fairly sensible way to do this is to compute the variance due to
phylogenetic uncertainty as the variance among estimates obtained from
the trees of the posterior sample; and then to compute the sampling
variance as the mean variance of the estimator from each tree; and then
add the two variances them. The standard error of our estimate (computed
as the mean across trees) is the square-root of this variance. To
conduct a hypothesis test on the regression coefficient, then, you would
compute the mean across trees and then the ratio of the parameter and
its standard error should have a t-distribution with n-2 degrees of
freedom for n taxa (not contrasts).
To do this from a practical perspective from trees in 'multiPhylo'
object (here "trees") and data in x & y for a simple bivariate
regression, we do the following:
# first define the following custom function
ff<-function(tree,x,y){
pic.x<-pic(x,tree)
pic.y<-pic(y,tree)
fit<-lm(pic.y~pic.x-1)
setNames(c(coef(fit),vcov(fit)),c("beta","var(beta)"))
}
# now apply to all trees in your sample
BB<-t(sapply(trees,ff,x,y))
# total variance in beta estimated by
varBeta<-var(BB[,"beta"])+mean(BB[,"var(beta)"])
t.beta<-mean(BB[,"beta"])/sqrt(varBeta)
P.beta<-2*pt(abs(t.beta),df=length(trees[[1]]$tip)-2,lower.tail=FALSE)
I think that's right.
All the best, Liam
Liam J. Revell, Assistant Professor of Biology
University of Massachusetts Boston
web: http://faculty.umb.edu/liam.revell/
email: [email protected]
blog: http://blog.phytools.org
On 4/1/2013 11:10 AM, Joe Felsenstein wrote:
John D asked:
Given that things have been quiet in the list lately, I think this
could be a good time for me to ask your opinion about this issue.
Imagine that I run a PIC analysis on two traits using 1000
post-burn-in trees. What would be the best way to summarize these
results? Average p-values across all analyses? Perhaps a specific
method to combine the resulting probabilities e.g. Fisher's test?
If you want to know about a regression coefficient, or other
parameter, the 1000 inferred values in the 1000 trees can be
taken as a posterior distribution if the trees come from a
Bayesian inference (I am assuming from your description that
they do).
If you want to test whether, say, the regression coefficient is
positive, you can take this posterior distribution and simply
ask whether the credible interval for the parameter includes
zero. In this case that interval would be the upper 95% of the
sample values.
Bringing in P values or Fisher's test just mixes Bayesian with
frequentist or likelihood methods, which will get two groups of
people annoyed at you, whereas using only one of these methods
has the advantage of irritating only one group of us.
I note that John's email address is "dobzhanski". Good to hear
from you again, Professor Dobzhansky. You may recall that you
gave me a tour of your lab at Rockefeller University in the
summer of 1963 when Dick Lewontin sent me from Rochester to
pick up some data paperwork for him for a joint project the two
of you were working on.
(Okay, the last part is "Poisson d'Avrile")
----
Joe Felsenstein, [email protected]
Dept. of Genome Sciences, Univ. of Washington
Box 355065, Seattle, WA 98195-5065 USA
_______________________________________________
R-sig-phylo mailing list - [email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/[email protected]/
_______________________________________________
R-sig-phylo mailing list - [email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/[email protected]/
