Hi Sereina.
Why lambda=0.5? Normally investigators tend to compare a model where
lambda is estimated to one in which it is fixed at 1 which corresponds
to Brownian evolution; or 0 which corresponds to no phylogenetic
correlation in the residual error of the model.
We can compare two fitted models, one in which lambda is estimated and
another in which lambda is fixed (at 0, 1, or some other arbitrary
value) using a likelihood ratio test. This is what that would look like:
## (Note that this uses gls in nlme instead of pgls
## in caper)
## fit models
## DATA is data frame containing x & y
## tree is object of class "phylo"
library(nlme)
fitBM<-gls(y~x,data=DATA,correlation=corBrownian(1,tree),
method="ML")
fitLambda<-gls(y~x,data=DATA,correlation=corPagel(1,tree,
fixed=FALSE),method="ML")
## function for likelihood ratio test
lrtest<-function(model1,model2){
lik1<-logLik(model1)
lik2<-logLik(model2)
LR<--2*(lik1-lik2)
degf<-attr(lik2,"df")-attr(lik1,"df")
P<-pchisq(LR,df=degf,lower.tail=FALSE)
cat(paste("Likelihood ratio = ",
signif(LR,5),"(df=",degf,") P =",
signif(P,4),"\n",sep=" "))
invisible(list(likelihood.ratio=LR,p=P))
}
## run likelihood-ratio test
lrtest(fitBM,fitLambda)
For small trees, we may want to generate a null distribution for the LR
using simulation under the null rather than the parametric distribution.
All the best, Liam
Liam J. Revell, Assistant Professor of Biology
University of Massachusetts Boston
web: http://faculty.umb.edu/liam.revell/
email: liam.rev...@umb.edu
blog: http://blog.phytools.org
On 8/20/2014 9:32 AM, Sereina Graber wrote:
Dear all,
I have a question conserning the pgls regression in package caper. The
function allows to estimate or fix three branch length transformations.
I wanna figure out which transformation gives me the best model fit by
comparing for example a model with lambda estimated (lambda="ML") to a
model where I fix lambda at 0.5 (lambda=0.5). However, if I use
anova(model1, model2) to compare the two models I get the error message
that the two models are run with different branch length
transformations, which makes sense. But is there any other possbility to
compare models with different branch length transformations? How do I
know which model is better?
For any help I am very grateful.
Best,
Sereina
_______________________________________________
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/
_______________________________________________
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/
