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://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
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
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/
[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
(((D.scabricollis:0.0281398198,D.punctiventris:0.0522260568):0.0165304832,(D.chetri:0.0324715501,D.srivichaii:0.036577528):0.0042592543):0.016112973,(D.obliquenotatus:0.0347332056,D.siwalikensis:0.0240496956):0.0356443411):0.0076514838,(D.fornicifrons:0.0673765515,D.nitidicollis:0.1004560337):0.0145728323):0.0121603904,(((D.karen:0.0504328722,D.pseudacutus:0.0471609249):2.7666e-06,D.andrewesi:0.0480975507):0.0069669242,D.tortuosus:0.0382059675):0.010811741):0.0020684997,(D.luteolunatus:0.0562242167,D.siamensis:0.1021470027):0.010950687):0.009366087,(D.betzi:0.0402231372,D.meo:0.0430754102):0.0414707113):0.0231370948,D.tubericollis:0.0883745395):0.0084178945,(S.amoenus:0.0018564107,S.notaculipennis:0.0010542338):0.0781434518,(S.pustulatus:0.0884092627,S.puthzianus:0.0623078578):0.0266474282):0.0198965576,((S.angusticollis:0.0131904559,S.thanonensis:0.0197882377):0.0669117612,S.humeralis:0.0679819232):0.0160558121):0.0160685651,S.rorellus:0.0968319412):0.0325738113,S.biplagiatus:0.0146882231,S.virgula:2.7327e-06):0.0181242578,S.luteomaculatus:0.041652629):0.0501030856,S.subthoracicus:0.0702948875):0.0192588442,((S.explanipennis:0.0669204403,S.liangtangi:0.0837475091):0.0191896571,S.stigmatias:0.0859050079):0.027375944):0.0340151206):0.0058362379):0.010800078,(((S.bivulneratus:0.0636351848,(S.diversiventris:2.4147e-06,S.gestroi:0.0641563764):0.0797494813):0.0277285246,S.subsimilis:0.0642784929):0.0167015008,S.feae:0.0607533321):0.0472946346):0.034510427,(S.diffidens:0.0021005521,S.megacephalus:0.0034191894):0.1022193644);Species;LOG_rel_Breite_4.Vordertarsus;LOG_rel_Anzahl_Hafthaare_4.Vordertarsus;LOG_rel_Attachment_force_smooth
D.andrewesi;0,033620906;0,897109673;1,206825876
D.betzi;0,024434016;0,936221491;2,01494035
D.chetri;0,042472404;0,700346373;1,01709
D.fornicifrons;0,046418895;0,952039522;1,149219113
D.karen;0,036448989;0,935587553;1,874481818
D.luteolunatus;0,03592642;0,815651322;1,086359831
D.meo;0,034042043;0,917882193;1,644438589
D.nitidicollis;0,046209628;0,697899073;0,968482949
D.obliquenotatus;0,034006252;0,781513275;1,686636269
D.pseudacutus;0,037919024;0,914995734;2,027757205
D.scabricollis;0,030655227;0,834123208;1,222716471
D.siamensis;0,063980241;1,042864906;2,201943063
D.siwalikensis;0,029994301;0,730857451;1,117271296
D.srivichaii;0,04035635;0,879907419;1,64246452
D.tortuosus;0,031486925;0,824229931;1,408239965
D.tubericollis;0,025100961;0,840167747;1,857332496
S.amoenus;0,042810582;1,107036321;1,40654018
S.diversiventris;0,037265779;1,087117381;1,269512944
S.humeralis;0,05243937;1,290835955;1,7355989
S.luteomaculatus;0,03659563;1,057227453;1,460897843
S.pustulatus;0,041894643;1,128142516;1,752816431
S.puthzianus;0,039201731;1,039932755;1,503790683
S.thanonensis;0,046789187;1,181896504;1,7355989
#PGLS
#Data management to prepare correct data and tree file
library(ape)
library(nlme)# regression modelling
library(geiger)
library(phytools)
library(OUwie)
library(picante)
library(slouch)
library (phylolm)
library (caper)
#shows curent working directory
getwd()
# my folder is at the top level of my user directory
setwd("C:/temp")
#Reding the data
SteninaeData<-read.csv2("Steninae1_PGLS_attachment force smooth versus number
tenent hairs.csv",header=TRUE)
names(SteninaeData) <- c("species", "breite", "anzahl", "force") # rename
colums with shorter names, save dataframe as "dat"
rownames(SteninaeData) <- SteninaeData$species # use species names to
set rownames
head(SteninaeData)
#reads tree
SteninaeTree<-read.tree("DataR/Steninae_reduced.nwk")
#Plot Tree
plot(SteninaeTree,no.margin=TRUE,edge.width=2)
#compares csv and tree for consistency
obj<-name.check(SteninaeTree,SteninaeData)
obj
#prune taxa without data from tree
SteninaeTree.cortado<-drop.tip(SteninaeTree, obj$tree_not_data)
#check again
name.check(SteninaeTree.cortado,SteninaeData)
