Hi Carl (also Brian),
This can effectively be accomplished in the geiger function
fitContinuous() by setting the upper & lower bounds for the parameter to
both be equal to the value for which you want to calculate the
likelihood. Unfortunately, these are not allowed to be exactly the
same; however, they can be arbitrarily close. So, for example, given
tree=tree, data=x, and a desired BM rate parameter of 10, you can compute:
>
results<-fitContinuous(tree,x,model="BM",bounds=list(beta=c(10.00000000,10.00000001)));
# for example
You can also just calculate the log-likelihood directly as follows:
> sig2<-desired.value; # set your desired value for the BM rate, say
sig2<-10; as before.
> C<-vcv.phylo(tree); C<-C[names(x),names(x)]; # compute VCV for tree,
and assure that the rows & columns are in the same order as x
> n<-nrow(C); ones<-matrix(1,n,1); # compute some preliminaries
>
a<-solve(t(ones)%*%solve(sig2*C)%*%ones)%*%t(ones)%*%solve(sig2*C)%*%x;
# estimate ancestor
>
logL<--(1/2)*(t(x-ones%*%a)%*%solve(sig2*C)%*%(x-ones%*%a))-(n/2)*log(2*pi)-(1/2)*log(det(sig2g*C));
# compute the log-likelihood
Hope this is helpful!
- Liam
Liam J. Revell
NESCent, Duke University
web: http://anolis.oeb.harvard.edu/~liam/
NEW email: lrev...@nescent.org
Carl Boettiger wrote:
Dear r-sig-phylo,
I'd like to be able to estimate the likelihood of a model with a given
Brownian rate parameter under a given data set. For OU models, this can be
accomplished through the ouch package with the flag fit=FALSE. However
there does not seem to be a corresponding option for Brownian motion in
ouch/geiger/ape algorithms for fitting a Brownian motion model? Perhaps
I've overlooked something?
-Carl
_______________________________________________
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
