[R-sig-phylo] Fwd: Error message when running gls with estimated lambda

Zenil-Ferguson,Rosana

2016-01-06 19:49:09 UTC

Hi Antonella,
You should try optimizations starting with multiple lambdas and see if you reach the same results in the optimization. That is a sign that there is a global maximum likelihood in your dataset.
You can look at the profile likelihood of Pagel's lambda, it might not be nice as Emmanuel Paradis suggested and that is giving you trouble when optimizing.
To look at the (relative) profile likelihood of lambda you can use the following code

library(MASS)
library(ape)
library(nlme)

lambda0<-seq(0, 1, 0.005)
long1<-length(lambda0)
log.likes<-rep(0,long1)
for(i in 1:long1){
model.i<-gls(y~x,correlation=corPagel(lambda0[i], phy.tree,fixed=TRUE),method="ML")
log.likes[i]<-logLik(model.i)[1] # y~x is your model, phy.tree your phylogeny
}
rel.loglike<-exp(log.likes-(max(log.likes)))

plot(lambda0,rel.loglike, type="l", xlab=paste("Pagel's", expression(lambda),sep=" "),ylab="Relative profile likelihood")

It is possible that the likelihood of lambda looks flat or multimodal that tells you that you don't have sufficient information to determine phylogenetic signal, and that might be messing all your estimates

Best,
-Rosana Zenil-Ferguson
PhD Candidate. University of Florida

________________________________________
From: R-sig-phylo <r-sig-phylo-***@r-project.org> on behalf of Antonella Soro <***@zoologie.uni-halle.de>
Sent: Wednesday, January 6, 2016 12:39 PM
To: r-sig-***@r-project.org
Subject: [R-sig-phylo] Fwd: Error message when running gls with estimated lambda

dear all,

a couple of month ago i have sent this message to [R-sig-phylo]:

>>> Le 02/11/2015 18:29, Antonella Soro a écrit :

>>>> dear all,
>>>> like kari allen in june 2011, i am trying to "to run Pagels GLS in R and I have consistently encountered
>>>> an error message (see script and error message below).
>>>> "
>>>>> pagel.data<-corPagel(1,phy=gls.phylo1,fixed=FALSE)
>>>>> pagel.gls8<-gls(Gst~log.Body_size+Diet+Behavior+Max_Geo+He,correlation=pagel.data,data=data,na.action=na.exclude, method="ML")
>>>> Error in corFactor.corStruct(object) :
>>>> NA/NaN/Inf in foreign function call (arg 1)

>>>

i have now found out that the code above works for my data if i change the (initial) value of the parameter lambda from “1” to “0”. given this positive result, i am very much tempted to change that value and carry on with my analyses, but i wonder how entitled i am in making that change. what is the meaning of giving “1” as initial value of the parameter lambda, rather than any other value between “0” and “1"? and why would the outcome of the analysis depend so drastically (error vs not error) on that initial value?

emmanuel paradies’s response to the same questions on the legitimacy of changing the initial value of the parameter lambda was:

“Hi Antonella,

I see no problem in changing the initial value of lambda as long as it gives an estimate. I do not know if there is a reason to expect a "regular" shape of the likelihood function for this kind of problem; so adjusting the initial values of the parameters to "help" the optimization algorithm is not a bad thing. You can send this question on r-sig-phylo; maybe others will have another opinion"

any comment or a different opinion from the one of emmanuel is very much appreciated

thank you in advance

antonella

_______________________________________________
R-sig-phylo mailing list - R-sig-***@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-***@r-project.org/

_______________________________________________
R-sig-phylo mailing list - R-sig-***@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-***@r-project.org/