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,017033339
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)
#_____________________________________________________________________________
#APE for PGLS to allow regression model taking account of phylogenetic
correlation
# Expected covariance structure that allows regression model taking account of
phylogenetic correlation among species under Pagel lambda model
library(nlme)
cp<-corPagel(1, SteninaeTree.cortado, form =~species)
cp
model<-gls(force ~ anzahl, data=SteninaeData, correlation=cp)#regress 2
variables under Pagel Lambda
summary(model)
#____________________________________________________________________________
#PGLS with caper
comp.data<-comparative.data(SteninaeTree.cortado, SteninaeData,
names.col="species", vcv.dim=2, warn.dropped=TRUE)
model<-pgls(force~anzahl, data=comp.data, lambda = "ML", param.CI = 0.95,
control = list(fnscale=-1), bounds = NULL)
summary(model)
_______________________________________________
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/
