Re: [R-sig-phylo] Question regarding PGLS in caper

[Alejandro Gonzalez Voyer]

Hello Oliver,

I completely agree with Julien. I’d add however that when running a PGLS with 
different branch length transformations you are in fact comparing among 
different models of trait evolution, which make different assumptions about how 
traits change along the phylogeny. Thus, I would add that if you do go down the 
route of comparing such models, you need to think carefully first whether the 
models actually make sense for the traits you are comparing. 

Cheers

Alejandro
___
Dr Alejandro Gonzalez Voyer
Investigador Titular B

Coordinación de Docencia y Formación de Recursos Humanos
Instituto de Ecología

Laboratorio de Conducta Animal
Instituto de Ecología
Circuito Exterior S/N
Ciudad Universitaria
Universidad Nacional Autónoma de México
México, D.F.
04510
México

Tel: +52 55 5622 9044
E-mail: alejandro.gonza...@iecologia.unam.mx

Web Site: http://alejandrogonzalezvoyer.com 


> El 19 oct 2021, a las 9:22, Julien Clavel  escribió:
> 
> Hi Oliver,
> 
> Your model comparison tells you that the "lambda" is better than "kappa", and 
> hence that it should provides a better phylogenetic structure for the 
> corresponding regression test.
> 
> Making a model comparison to choose the structure is probably best practice, 
> but it also depends on the practicability. If I would have to stuck with a 
> single model from these three one (kind of omnibus transform) I would 
> probably choose "lambda" then, because at one extreme it reduces to the 
> classical OLS (lambda=0) and intermediate values (strength of "phylogenetic 
> signal") correspond to using a mixed model with BM random effects. This is 
> easier to interpret and directly relate to the signal to noise ratio of 
> phylogenetic structure in your model residuals.
> 
> Cheers,
> 
> Julien
> 
> 
> De : R-sig-phylo  > de la part de Oliver Betz 
> mailto:oliver.b...@uni-tuebingen.de>>
> Envoy� : mardi 19 octobre 2021 01:02
> � : r-sig-phylo@r-project.org  
> mailto:r-sig-phylo@r-project.org>>
> Objet : [R-sig-phylo] Question regarding PGLS in caper
> 
> 
> Dear all:
> 
> I am applying the pgls function in caper to test for the relationship
> between a dependent and an independent variable. Usually, the
> regressions can be well interpreted when applying the lambda branch
> lenghth transformation. However, in one case, the regression became
> only significant when I applied the kappa transformation, but remained
> non-significant when I applied the lambda transformation. Does this
> automatically mean that I should choose the kappa transformation in
> this case? Or should I still stick to the AICc values and always
> choose the model with the lowest AICc and the highest model weight? In
> the described case, AICc for the (non-significant) lambda model was
> lowest, i.e. it was lower (-36.7) than the (significant) kappa model
> (-32.8).
> 
> This brings me to a more general question: when using pgls in caper,
> should I always compare (via AICc) bewteen the three branch length
> transformation models (lambda, kappa, delta), or would it be
> sufficient to stick to lambda to get an idea of the phylognetic
> non-indepedence of the residual errors of the regression model?
> 
> Thank you for any suggestions,
> 
> Oliver Betz
> 
> University of T�bingen, Germany
> Institute of Evolution and Ecology
> 
> ___
> 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/ 
> 
> 
>   [[alternative HTML version deleted]]
> 
> ___
> 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/ 
> 

[[alternative HTML version deleted]]

___
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/

Re: [R-sig-phylo] Question regarding PGLS in caper

[Julien Clavel]

Hi Oliver,

Your model comparison tells you that the "lambda" is better than "kappa", and 
hence that it should provides a better phylogenetic structure for the 
corresponding regression test.

Making a model comparison to choose the structure is probably best practice, 
but it also depends on the practicability. If I would have to stuck with a 
single model from these three one (kind of omnibus transform) I would probably 
choose "lambda" then, because at one extreme it reduces to the classical OLS 
(lambda=0) and intermediate values (strength of "phylogenetic signal") 
correspond to using a mixed model with BM random effects. This is easier to 
interpret and directly relate to the signal to noise ratio of phylogenetic 
structure in your model residuals.

Cheers,

Julien


De : R-sig-phylo  de la part de Oliver Betz 

Envoy� : mardi 19 octobre 2021 01:02
� : r-sig-phylo@r-project.org 
Objet : [R-sig-phylo] Question regarding PGLS in caper


Dear all:

I am applying the pgls function in caper to test for the relationship
between a dependent and an independent variable. Usually, the
regressions can be well interpreted when applying the lambda branch
lenghth transformation. However, in one case, the regression became
only significant when I applied the kappa transformation, but remained
non-significant when I applied the lambda transformation. Does this
automatically mean that I should choose the kappa transformation in
this case? Or should I still stick to the AICc values and always
choose the model with the lowest AICc and the highest model weight? In
the described case, AICc for the (non-significant) lambda model was
lowest, i.e. it was lower (-36.7) than the (significant) kappa model
(-32.8).

This brings me to a more general question: when using pgls in caper,
should I always compare (via AICc) bewteen the three branch length
transformation models (lambda, kappa, delta), or would it be
sufficient to stick to lambda to get an idea of the phylognetic
non-indepedence of the residual errors of the regression model?

Thank you for any suggestions,

Oliver Betz

University of T�bingen, Germany
Institute of Evolution and Ecology

___
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/

[[alternative HTML version deleted]]

___
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/

Re: [R-sig-phylo] question regarding PGLS

[Julien Clavel]

Just to illustrate it:

values <- seq(0,to=1,by=0.01) # grid of lambda values
var_t = diag(vcv(SteninaeTree.cortado))  # heights from the root to the tips
ll <- sapply(values, function(x){
  cp<-corPagel(x, SteninaeTree.cortado, form =~species, fixed = TRUE)
  model_ml<-gls(force ~ anzahl, data=SteninaeData, correlation=cp, 
weights=varFixed(~var_t), method = "ML")
  model_reml<-gls(force ~ anzahl,  data=SteninaeData, correlation=cp, 
weights=varFixed(~var_t), method = "REML")
  c(logLik(model_ml), logLik(model_reml))
})
par(mfrow=c(1,2))
plot(ll[1,]~values, main="ML", type="l", xlab="lambda", ylab="Log-lik")
plot(ll[2,]~values, main="REML", type="l", xlab="lambda", ylab="Log-lik")
values[apply(ll,1,which.max)] # best ML and REML values

You can try on simulated data (e.g. fit with a BM model, use the empirical 
sigma^2 to add BM noise to the predicted values) to appreciate how both 
approaches will behave on that tree. You will likely find that some lambda 
estimates will often hit the zero bound under ML if your tree is rather small 
(or/and that a lot of branching events are deep).
Hope this helps

Julien

De : R-sig-phylo  de la part de Simone 
Blomberg 
Envoyé : lundi 31 mai 2021 03:39
À : r-sig-phylo@r-project.org 
Objet : Re: [R-sig-phylo] question regarding PGLS 
 
Another issue is that the tree is not ultrametric. If you use nlme::gls 
to fit the model and you have a non-ultrametric tree, you need to use 
the weights argument to pass the tip heights. Otherwise gls() assumes 
the tree is ultrametric. This could be part of the problem. I wrote a 
post about this in this forum several years ago, but I can no longer 
find it in the archive. The REML/ML distinction might also be important 
here.

Cheers,

Simone.

On 31/5/21 10:42 am, Julien Clavel wrote:
> Hi Oliver,
>
> The "gls" from nlme uses REML by default. I think that caper or phylolm use 
> ML instead. On small sample size ML estimates (e.g., "lambda") are known to 
> be biased and you may sometimes not have enough power to estimate them. With 
> increasing sample size (bigger trees) this is less an issue as the bias 
> should vanishes and ML will likely converge to similar values to REML.
>
> Also, “gls” estimates a correlation rather than a covariance structure. On 
> non-ultrametric trees (such as yours) this will lead to different results.
>
> Regards,
>
> Julien
>
>
>
>
> De : R-sig-phylo  de la part de Oliver 
> Betz 
> Envoyé : dimanche 30 mai 2021 22:49
> À : r-sig-phylo@r-project.org 
> Objet : [R-sig-phylo] question regarding PGLS
>   
>
> Dear list members:
>
>
> I tried various R packages to calcuate a PGLS with the data set (csv
> and nwk) I have attached to this email. I would like to use the Pagels
> lambda model to attain an index that measures whether data exhibit
> phylogenetic dependence or not.
>
> While doing so, I came up with the following problem. Using ape - for
> example, according to the script of the Latin American Macroevolution
> Workshop and others - I attained results (regarding intercept, slope
> and lambda) that differ from the results I received from caper or
> phylolm. For example, with ape, lambda amounts to 0.57, whereas
> according to caper or phylolm lambda amounts to 0.
>
> Please find the R script that I have used attached (first block is
> just the data reading and organisational part; 2nd block: PGLS Pagels
> Lambda according to ape, 3rd block: PGLS Pagels Lambda according to
> caper).
>
> Do you have any idea why these packages produce these different
> results and which software I should follow?
>
> Thanks for any suggestions.
>
> Oliver Betz
> University of Tübingen
> Tübingen
> Germany
> ___
> 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/

-- 
Simone Blomberg, BSc (Hons), PhD, MAppStat.
Senior Lecturer and Consultant Statistician
School of Biological Sciences
The University of Queensland
St. Lucia Queensland 4072
Australia
T: +61 7 3365 2506
email: S.Blomberg1_at_uq.edu.au
Twitter: @simoneb66
UQ ALLY Supporting the diversity of sexuality and gender at UQ.

Policies:
1.  I will NOT analyse your data for you.
2.  Your deadline is your problem.



Basically, I'm not interested in doing research
and I never have been. I'm interested in
understanding, which is quite a different thing.
And often to understand something you have to
work it out for yourself because no one else
has done it. - David Blackwell

___
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https:

Re: [R-sig-phylo] question regarding PGLS

[Simone Blomberg]

Another issue is that the tree is not ultrametric. If you use nlme::gls 
to fit the model and you have a non-ultrametric tree, you need to use 
the weights argument to pass the tip heights. Otherwise gls() assumes 
the tree is ultrametric. This could be part of the problem. I wrote a 
post about this in this forum several years ago, but I can no longer 
find it in the archive. The REML/ML distinction might also be important 
here.


Cheers,

Simone.

On 31/5/21 10:42 am, Julien Clavel wrote:

Hi Oliver,

The "gls" from nlme uses REML by default. I think that caper or phylolm use ML instead. 
On small sample size ML estimates (e.g., "lambda") are known to be biased and you may 
sometimes not have enough power to estimate them. With increasing sample size (bigger trees) this 
is less an issue as the bias should vanishes and ML will likely converge to similar values to REML.

Also, “gls” estimates a correlation rather than a covariance structure. On 
non-ultrametric trees (such as yours) this will lead to different results.

Regards,

Julien




De : R-sig-phylo  de la part de Oliver Betz 

Envoyé : dimanche 30 mai 2021 22:49
À : r-sig-phylo@r-project.org 
Objet : [R-sig-phylo] question regarding PGLS
  


Dear list members:


I tried various R packages to calcuate a PGLS with the data set (csv
and nwk) I have attached to this email. I would like to use the Pagels
lambda model to attain an index that measures whether data exhibit
phylogenetic dependence or not.

While doing so, I came up with the following problem. Using ape - for
example, according to the script of the Latin American Macroevolution
Workshop and others - I attained results (regarding intercept, slope
and lambda) that differ from the results I received from caper or
phylolm. For example, with ape, lambda amounts to 0.57, whereas
according to caper or phylolm lambda amounts to 0.

Please find the R script that I have used attached (first block is
just the data reading and organisational part; 2nd block: PGLS Pagels
Lambda according to ape, 3rd block: PGLS Pagels Lambda according to
caper).

Do you have any idea why these packages produce these different
results and which software I should follow?

Thanks for any suggestions.

Oliver Betz
University of Tübingen
Tübingen
Germany
___
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/


--
Simone Blomberg, BSc (Hons), PhD, MAppStat.
Senior Lecturer and Consultant Statistician
School of Biological Sciences
The University of Queensland
St. Lucia Queensland 4072
Australia
T: +61 7 3365 2506
email: S.Blomberg1_at_uq.edu.au
Twitter: @simoneb66
UQ ALLY Supporting the diversity of sexuality and gender at UQ.

Policies:
1.  I will NOT analyse your data for you.
2.  Your deadline is your problem.



Basically, I'm not interested in doing research
and I never have been. I'm interested in
understanding, which is quite a different thing.
And often to understand something you have to
work it out for yourself because no one else
has done it. - David Blackwell

___
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/

Re: [R-sig-phylo] question regarding PGLS

[Liam J. Revell]

Dear Oliver & Julien.

> Also, “gls” estimates a correlation rather than a covariance
> structure. On non-ultrametric trees (such as yours) this will lead to
> different results.

This is a great point. If your tree is non-ultrametric you can do 
something like:


w<-diag(vcv.phylo(tree))
spp<-names(w)
corLambda<-corPagel(1,phy=tree,form=~spp)
fit<-gls(...,correlation=corLambda,weights=varFixed(~w))

(I'm not sure the exact code, but it should be something like that.)

All the best, Liam

Liam J. Revell
University of Massachusetts Boston [Assoc. Prof.]
Universidad Católica de la Ssma Concepción [Adj. Res.]

Web & phytools:
http://faculty.umb.edu/liam.revell/, http://www.phytools.org, 
http://blog.phytools.org


Academic Director UMass Boston Chile Abroad:
https://www.umb.edu/academics/caps/international/biology_chile

U.S. COVID-19 explorer web application:
https://covid19-explorer.org/

On 5/30/2021 8:42 PM, Julien Clavel wrote:

EXTERNAL SENDER

Hi Oliver,

The "gls" from nlme uses REML by default. I think that caper or phylolm use ML instead. 
On small sample size ML estimates (e.g., "lambda") are known to be biased and you may 
sometimes not have enough power to estimate them. With increasing sample size (bigger trees) this 
is less an issue as the bias should vanishes and ML will likely converge to similar values to REML.

Also, “gls” estimates a correlation rather than a covariance structure. On 
non-ultrametric trees (such as yours) this will lead to different results.

Regards,

Julien




De : R-sig-phylo  de la part de Oliver Betz 

Envoyé : dimanche 30 mai 2021 22:49
À : r-sig-phylo@r-project.org 
Objet : [R-sig-phylo] question regarding PGLS


Dear list members:


I tried various R packages to calcuate a PGLS with the data set (csv
and nwk) I have attached to this email. I would like to use the Pagels
lambda model to attain an index that measures whether data exhibit
phylogenetic dependence or not.

While doing so, I came up with the following problem. Using ape - for
example, according to the script of the Latin American Macroevolution
Workshop and others - I attained results (regarding intercept, slope
and lambda) that differ from the results I received from caper or
phylolm. For example, with ape, lambda amounts to 0.57, whereas
according to caper or phylolm lambda amounts to 0.

Please find the R script that I have used attached (first block is
just the data reading and organisational part; 2nd block: PGLS Pagels
Lambda according to ape, 3rd block: PGLS Pagels Lambda according to
caper).

Do you have any idea why these packages produce these different
results and which software I should follow?

Thanks for any suggestions.

Oliver Betz
University of Tübingen
Tübingen
Germany
___
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-phylodata=04%7C01%7Cliam.revell%40umb.edu%7Ca0076f879f9041059c2008d923ccf9bf%7Cb97188711ee94425953c1ace1373eb38%7C0%7C1%7C637580185553969207%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000sdata=3pSEeLqSLv2l7K0fgnUCu0hgUiCQYPpOuu0t%2FpQSvJc%3Dreserved=0
Searchable archive at 
https://nam10.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.mail-archive.com%2Fr-sig-phylo%40r-project.org%2Fdata=04%7C01%7Cliam.revell%40umb.edu%7Ca0076f879f9041059c2008d923ccf9bf%7Cb97188711ee94425953c1ace1373eb38%7C0%7C1%7C637580185553969207%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000sdata=lqErdhi2GLfaN5%2F01McrOWg2ePmWeuKnmFDB%2B%2BStolI%3Dreserved=0



___
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/

Re: [R-sig-phylo] question regarding PGLS

[Julien Clavel]

Hi Oliver,

The "gls" from nlme uses REML by default. I think that caper or phylolm use ML 
instead. On small sample size ML estimates (e.g., "lambda") are known to be 
biased and you may sometimes not have enough power to estimate them. With 
increasing sample size (bigger trees) this is less an issue as the bias should 
vanishes and ML will likely converge to similar values to REML.

Also, “gls” estimates a correlation rather than a covariance structure. On 
non-ultrametric trees (such as yours) this will lead to different results.

Regards,

Julien




De : R-sig-phylo  de la part de Oliver Betz 

Envoyé : dimanche 30 mai 2021 22:49
À : r-sig-phylo@r-project.org 
Objet : [R-sig-phylo] question regarding PGLS 
 

Dear list members:


I tried various R packages to calcuate a PGLS with the data set (csv  
and nwk) I have attached to this email. I would like to use the Pagels  
lambda model to attain an index that measures whether data exhibit  
phylogenetic dependence or not.

While doing so, I came up with the following problem. Using ape - for  
example, according to the script of the Latin American Macroevolution  
Workshop and others - I attained results (regarding intercept, slope  
and lambda) that differ from the results I received from caper or  
phylolm. For example, with ape, lambda amounts to 0.57, whereas  
according to caper or phylolm lambda amounts to 0.

Please find the R script that I have used attached (first block is  
just the data reading and organisational part; 2nd block: PGLS Pagels  
Lambda according to ape, 3rd block: PGLS Pagels Lambda according to  
caper).

Do you have any idea why these packages produce these different  
results and which software I should follow?

Thanks for any suggestions.

Oliver Betz
University of Tübingen
Tübingen
Germany
___
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/