Thank you guys for you helpful advice!
The probelm in my model is that it estimates lambda to be 1, however, my gut feeling would say this is rather unlikely. The model where lambda equals 1 looses significance completely, however, with lambda fixed at 0.85, it is highly significant. So what do I believe now? This is why I wanted to compare different models and see whats happening...Do you think it is justified to set lambda to a lower value and then compare the model to the one where lambda equals 1?
I also checked the log likelihood curve of lambda: the lower CI boundary is very close to 1, however, the upper CI boundary is not given (=NA). So I don`t know what happens with the loglik beyond 1 which makes it hard to interpret and be sure whether the estimate is reliable or not, right? (plot attached)
Thank you & all the best,
Sereina
Gesendet: Mittwoch, 20. August 2014 um 17:11 Uhr
Von: "Liam J. Revell" <liam.rev...@umb.edu>
An: "Sereina Graber" <sereina.gra...@gmx.ch>, r-sig-phylo@r-project.org
Betreff: Re: [R-sig-phylo] caper model comparison with different branch length transformations
Von: "Liam J. Revell" <liam.rev...@umb.edu>
An: "Sereina Graber" <sereina.gra...@gmx.ch>, r-sig-phylo@r-project.org
Betreff: Re: [R-sig-phylo] caper model comparison with different branch length transformations
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=""> method="ML")
fitLambda<-gls(y~x,data=""> 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/
>
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=""> method="ML")
fitLambda<-gls(y~x,data=""> 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/
