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 <r-sig-phylo-boun...@r-project.org> de la part de Oliver Betz
<oliver.b...@uni-tuebingen.de>
Envoyé : dimanche 30 mai 2021 22:49
À : r-sig-phylo@r-project.org <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-phylo&data=04%7C01%7Cliam.revell%40umb.edu%7Ca0076f879f9041059c2008d923ccf9bf%7Cb97188711ee94425953c1ace1373eb38%7C0%7C1%7C637580185553969207%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=3pSEeLqSLv2l7K0fgnUCu0hgUiCQYPpOuu0t%2FpQSvJc%3D&reserved=0
Searchable archive at
https://nam10.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.mail-archive.com%2Fr-sig-phylo%40r-project.org%2F&data=04%7C01%7Cliam.revell%40umb.edu%7Ca0076f879f9041059c2008d923ccf9bf%7Cb97188711ee94425953c1ace1373eb38%7C0%7C1%7C637580185553969207%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=lqErdhi2GLfaN5%2F01McrOWg2ePmWeuKnmFDB%2B%2BStolI%3D&reserved=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/
