analysis of macroevolution with phylogenies
#save(sylvia.clus, taxa.sylvia, sylvia.eco, nj.est, file="D:\\OneDrive\\R packages\\Rcode\\sylvia.RData")
load(file="D:\\OneDrive\\R packages\\Rcode\\sylvia.RData")
library(ape)
nj.est$tip.label## Z73494 AJ534526 AJ534527
## "Sylvia_atricapilla" "Chamaea_fasciata" "Sylvia_nisoria"
## AJ534528 AJ534529 AJ534530
## "Sylvia_layardi" "Sylvia_subcaeruleum" "Sylvia_boehmi"
## AJ534531 AJ534532 AJ534533
## "Sylvia_buryi" "Sylvia_lugens" "Sylvia_leucomelaena"
## AJ534534 AJ534535 AJ534536
## "Sylvia_hortensis" "Sylvia_crassirostris" "Sylvia_curruca"
## AJ534537 AJ534538 AJ534539
## "Sylvia_nana" "Sylvia_communis" "Sylvia_conspicillata"
## AJ534540 AJ534541 AJ534542
## "Sylvia_deserticola" "Sylvia_balearica" "Sylvia_undata"
## AJ534543 AJ534544 AJ534545
## "Sylvia_cantillans" "Sylvia_melanocephala" "Sylvia_mystacea"
## AJ534546 AJ534547 AJ534548
## "Sylvia_melanothorax" "Sylvia_rueppelli" "Sylvia_abyssinica"
## AJ534549
## "Sylvia_borin"nj.est <- drop.tip(nj.est, "Chamaea_fasciata")
nj.est$tip.label## Z73494 AJ534527 AJ534528
## "Sylvia_atricapilla" "Sylvia_nisoria" "Sylvia_layardi"
## AJ534529 AJ534530 AJ534531
## "Sylvia_subcaeruleum" "Sylvia_boehmi" "Sylvia_buryi"
## AJ534532 AJ534533 AJ534534
## "Sylvia_lugens" "Sylvia_leucomelaena" "Sylvia_hortensis"
## AJ534535 AJ534536 AJ534537
## "Sylvia_crassirostris" "Sylvia_curruca" "Sylvia_nana"
## AJ534538 AJ534539 AJ534540
## "Sylvia_communis" "Sylvia_conspicillata" "Sylvia_deserticola"
## AJ534541 AJ534542 AJ534543
## "Sylvia_balearica" "Sylvia_undata" "Sylvia_cantillans"
## AJ534544 AJ534545 AJ534546
## "Sylvia_melanocephala" "Sylvia_mystacea" "Sylvia_melanothorax"
## AJ534547 AJ534548 AJ534549
## "Sylvia_rueppelli" "Sylvia_abyssinica" "Sylvia_borin"We drop the outgroup species(Chamaea fasciata) for which we have no ecological data.
sylvia.eco## mig.dist mig.behav geo.range
## Sylvia_abyssinica 0 resid trop
## Sylvia_atricapilla 5000 short temptrop
## Sylvia_borin 7500 long temptrop
## Sylvia_nisoria 5900 long temptrop
## Sylvia_curruca 5500 long temptrop
## Sylvia_hortensis 3400 long temptrop
## Sylvia_crassirostris 2600 long temptrop
## Sylvia_leucomelaena 0 resid trop
## Sylvia_buryi 0 resid trop
## Sylvia_lugens 0 resid trop
## Sylvia_layardi 0 resid trop
## Sylvia_subcaeruleum 0 resid trop
## Sylvia_boehmi 0 resid trop
## Sylvia_nana 2600 short temptrop
## Sylvia_deserti 0 resid temp
## Sylvia_communis 6500 long temptrop
## Sylvia_conspicillata 1500 short temptrop
## Sylvia_deserticola 500 resid temp
## Sylvia_undata 300 resid temp
## Sylvia_sarda 300 resid temp
## Sylvia_balearica 0 resid temp
## Sylvia_cantillans 3500 long temptrop
## Sylvia_mystacea 2000 long temptrop
## Sylvia_melanocephala 300 short temptrop
## Sylvia_rueppelli 2400 long temptrop
## Sylvia_melanothorax 1200 resid tempDF <- sylvia.eco[nj.est$tip.label,]
DF## mig.dist mig.behav geo.range
## Sylvia_atricapilla 5000 short temptrop
## Sylvia_nisoria 5900 long temptrop
## Sylvia_layardi 0 resid trop
## Sylvia_subcaeruleum 0 resid trop
## Sylvia_boehmi 0 resid trop
## Sylvia_buryi 0 resid trop
## Sylvia_lugens 0 resid trop
## Sylvia_leucomelaena 0 resid trop
## Sylvia_hortensis 3400 long temptrop
## Sylvia_crassirostris 2600 long temptrop
## Sylvia_curruca 5500 long temptrop
## Sylvia_nana 2600 short temptrop
## Sylvia_communis 6500 long temptrop
## Sylvia_conspicillata 1500 short temptrop
## Sylvia_deserticola 500 resid temp
## Sylvia_balearica 0 resid temp
## Sylvia_undata 300 resid temp
## Sylvia_cantillans 3500 long temptrop
## Sylvia_melanocephala 300 short temptrop
## Sylvia_mystacea 2000 long temptrop
## Sylvia_melanothorax 1200 resid temp
## Sylvia_rueppelli 2400 long temptrop
## Sylvia_abyssinica 0 resid trop
## Sylvia_borin 7500 long temptropMigratory behavior is tightly linked with geographical range.
table(DF$geo.range, DF$mig.behav)##
## long resid short
## temp 0 4 0
## temptrop 9 0 4
## trop 0 7 0syl.er <- ace(DF$geo.range, nj.est, type="d")
syl.er##
## Ancestral Character Estimation
##
## Call: ace(x = DF$geo.range, phy = nj.est, type = "d")
##
## Log-likelihood: -22.16276
##
## Rate index matrix:
## temp temptrop trop
## temp . 1 1
## temptrop 1 . 1
## trop 1 1 .
##
## Parameter estimates:
## rate index estimate std-err
## 1 5.4806 1.8448
##
## Scaled likelihoods at the root (type '...$lik.anc' to get them for all nodes):
## temp temptrop trop
## 0.01809245 0.90972802 0.07217953syl.sym <- ace(DF$geo.range, nj.est, type="d", model="SYM")## Warning in sqrt(diag(solve(h))): NaNs producedsyl.sym##
## Ancestral Character Estimation
##
## Call: ace(x = DF$geo.range, phy = nj.est, type = "d", model = "SYM")
##
## Log-likelihood: -20.38674
##
## Rate index matrix:
## temp temptrop trop
## temp . 1 2
## temptrop 1 . 3
## trop 2 3 .
##
## Parameter estimates:
## rate index estimate std-err
## 1 3.4911 1.9330
## 2 0.0000 NaN
## 3 10.1228 4.3518
##
## Scaled likelihoods at the root (type '...$lik.anc' to get them for all nodes):
## temp temptrop trop
## 0.001717209 0.777724754 0.220558038anova(syl.er, syl.sym)## Likelihood Ratio Test Table
## Log lik. Df Df change Resid. Dev Pr(>|Chi|)
## 1 -22.163 1
## 2 -20.387 3 2 3.552 0.1693co <- rep("grey", 24)
co[DF$geo.rang=="temp"] <- "black"
co[DF$geo.rang=="trop"] <- "white"plot(nj.est, "c" , FALSE, no.margin=TRUE, label.offset=1)
tiplabels(pch=22, bg=co, cex=2, adj=1)
nodelabels(thermo=syl.er$lik.anc, cex=0.8, piecol=c("black","grey","white"))
syl.er$lik.anc## temp temptrop trop
## [1,] 0.018092450 0.90972802 0.0721795346
## [2,] 0.002827070 0.94648874 0.0506841935
## [3,] 0.002947196 0.97695162 0.0201011864
## [4,] 0.005444586 0.98932612 0.0052292921
## [5,] 0.013922361 0.98475065 0.0013269915
## [6,] 0.012157331 0.98698747 0.0008551975
## [7,] 0.034829467 0.96014119 0.0050293437
## [8,] 0.002124209 0.99730644 0.0005693477
## [9,] 0.001710543 0.99723969 0.0010497684
## [10,] 0.090521686 0.90432522 0.0051530923
## [11,] 0.879786609 0.11406556 0.0061478306
## [12,] 0.985017601 0.01345206 0.0015303346
## [13,] 0.003522825 0.92142789 0.0750492821
## [14,] 0.003992921 0.61105246 0.3849546174
## [15,] 0.003342871 0.58904158 0.4076155530
## [16,] 0.010313012 0.48696744 0.5027195458
## [17,] 0.004415313 0.61121713 0.3843675593
## [18,] 0.002911273 0.56055168 0.4365370511
## [19,] 0.003647181 0.54939640 0.4469564223
## [20,] 0.006663643 0.09521884 0.8981175121
## [21,] 0.007375432 0.56153035 0.4310942181
## [22,] 0.026224038 0.87832504 0.0954509198
## [23,] 0.026870004 0.88059920 0.0925307967