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 <[email protected]> on behalf of Antonella
Soro <[email protected]>
Sent: Wednesday, January 6, 2016 12:39 PM
To: [email protected]
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 - [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]/
