The challenge problem was as follows.
Download the following two data files.
Use stochastic character mapping to generate 10 stochastic character maps of feeding
mode on the tree of elopomorphs, then use brownie.lite and
OUwie to fit a multi-rate & a multi-regime OU model to the continuous
character, transformed to a log-scale. What do you rind?
First, let's read the data:
library(phytools)
## Loading required package: ape
## Loading required package: maps
##
## # maps v3.1: updated 'world': all lakes moved to separate new #
## # 'lakes' database. Type '?world' or 'news(package="maps")'. #
tree<-read.tree("elopomorph.tre")
X<-read.csv("elopomorph.csv",row.names=1)
Next, let's generate 10 stochastic mapped trees for feeding mode on the tree using
the phytools function make.simmap:
fmode<-setNames(X[,1],rownames(X))
trees<-make.simmap(tree,fmode,model="ARD",nsim=10)
## make.simmap is sampling character histories conditioned on the transition matrix
##
## Q =
## bite suction
## bite -0.01560942 0.01560942
## suction 0.01749038 -0.01749038
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## bite suction
## 0.5 0.5
## Done.
trees
## 10 phylogenetic trees with mapped discrete characters
## visualize:
obj<-densityMap(trees,states=c("suction","bite"),plot=FALSE)
## sorry - this might take a while; please be patient
plot(obj,fsize=c(0.6,0.9),outline=TRUE)
To start, let's fit the multi-rate Brownian model to each of these trees in which the rate of our continuous trait depends on the state:
tl<-setNames(log(X[,2]),rownames(X))
fitBM<-list()
for(i in 1:length(trees)) fitBM[[i]]<-brownie.lite(trees[[i]],tl)
## finally, let's get the parameters out from our list of model fits
rates<-matrix(NA,length(trees),length(fitBM[[1]]$sig2.multiple))
for(i in 1:length(trees)) rates[i,]<-fitBM[[i]]$sig2.multiple[levels(fmode)]
colnames(rates)<-paste("sig2(",levels(fmode),")",sep="")
rates
## sig2(bite) sig2(suction)
## [1,] 0.013231637 0.01512746
## [2,] 0.008662904 0.02393120
## [3,] 0.014550128 0.01443698
## [4,] 0.011770278 0.01688145
## [5,] 0.010935225 0.01879003
## [6,] 0.013160133 0.01550023
## [7,] 0.013247143 0.01528550
## [8,] 0.011522100 0.01906208
## [9,] 0.009866175 0.01902837
## [10,] 0.012003249 0.01593864
Obviously, in an empirical study we would pull many other things out of the fitted model, such as its overall fit.
Now, we can move to the multi-optimum OU analysis.
First, let's build out OUwie data frame:
OUdata<-data.frame(Genus_species=rownames(X),Reg=fmode,X=tl)
OUdata
## Genus_species Reg
## Albula_vulpes Albula_vulpes suction
## Anguilla_anguilla Anguilla_anguilla suction
## Anguilla_bicolor Anguilla_bicolor suction
## Anguilla_japonica Anguilla_japonica suction
## Anguilla_rostrata Anguilla_rostrata suction
## Ariosoma_anago Ariosoma_anago suction
## Ariosoma_balearicum Ariosoma_balearicum suction
## Ariosoma_shiroanago Ariosoma_shiroanago suction
## Bathyuroconger_vicinus Bathyuroconger_vicinus bite
## Brachysomophis_crocodilinus Brachysomophis_crocodilinus bite
## Conger_japonicus Conger_japonicus suction
## Conger_myriaster Conger_myriaster suction
## Conger_verreaxi Conger_verreaxi suction
## Conger_wilsoni Conger_wilsoni suction
## Congresox_talabonoides Congresox_talabonoides suction
## Cynoponticus_ferox Cynoponticus_ferox bite
## Dysomma_anguillare Dysomma_anguillare bite
## Elops_saurus Elops_saurus suction
## Facciolella_gilbertii Facciolella_gilbertii bite
## Gavialiceps_taeniola Gavialiceps_taeniola bite
## Gnathophis_longicauda Gnathophis_longicauda suction
## Gorgasia_taiwanensis Gorgasia_taiwanensis suction
## Gymnothorax_castaneus Gymnothorax_castaneus bite
## Gymnothorax_flavimarginatus Gymnothorax_flavimarginatus bite
## Gymnothorax_kidako Gymnothorax_kidako bite
## Gymnothorax_moringa Gymnothorax_moringa bite
## Gymnothorax_pseudothyrsoideus Gymnothorax_pseudothyrsoideus bite
## Gymnothorax_reticularis Gymnothorax_reticularis bite
## Heteroconger_hassi Heteroconger_hassi suction
## Ichthyapus_ophioneus Ichthyapus_ophioneus bite
## Kaupichthys_hyoproroides Kaupichthys_hyoproroides bite
## Kaupichthys_nuchalis Kaupichthys_nuchalis bite
## Megalops_cyprinoides Megalops_cyprinoides suction
## Moringua_edwardsi Moringua_edwardsi bite
## Moringua_javanica Moringua_javanica bite
## Muraenesox_bagio Muraenesox_bagio bite
## Muraenesox_cinereus Muraenesox_cinereus bite
## Myrichthys_breviceps Myrichthys_breviceps suction
## Myrichthys_maculosus Myrichthys_maculosus suction
## Myrichthys_magnificus Myrichthys_magnificus suction
## Myrophis_vafer Myrophis_vafer bite
## Nemichthys_scolopaceus Nemichthys_scolopaceus bite
## Nettastoma_melanurum Nettastoma_melanurum bite
## Ophichthus_serpentinus Ophichthus_serpentinus suction
## Ophichthus_zophochir Ophichthus_zophochir suction
## Oxyconger_leptognathus Oxyconger_leptognathus bite
## Parabathymyrus_macrophthalmus Parabathymyrus_macrophthalmus suction
## Paraconger_notialis Paraconger_notialis suction
## Pisodonophis_cancrivorus Pisodonophis_cancrivorus bite
## Poeciloconger_kapala Poeciloconger_kapala suction
## Rhinomuraena_quaesita Rhinomuraena_quaesita bite
## Rhynchoconger_flavus Rhynchoconger_flavus suction
## Saurenchelys_fierasfer Saurenchelys_fierasfer bite
## Scolecenchelys_breviceps Scolecenchelys_breviceps suction
## Scuticaria_tigrina Scuticaria_tigrina bite
## Serrivomer_beanii Serrivomer_beanii bite
## Serrivomer_sector Serrivomer_sector bite
## Simenchelys_parasitica Simenchelys_parasitica suction
## Uroconger_lepturus Uroconger_lepturus suction
## Uropterygius_micropterus Uropterygius_micropterus bite
## Venefica_proboscidea Venefica_proboscidea bite
## X
## Albula_vulpes 4.644391
## Anguilla_anguilla 3.912023
## Anguilla_bicolor 4.787492
## Anguilla_japonica 5.010635
## Anguilla_rostrata 5.023881
## Ariosoma_anago 4.094345
## Ariosoma_balearicum 3.555348
## Ariosoma_shiroanago 3.688879
## Bathyuroconger_vicinus 4.477337
## Brachysomophis_crocodilinus 4.787492
## Conger_japonicus 4.941642
## Conger_myriaster 4.605170
## Conger_verreaxi 5.298317
## Conger_wilsoni 5.010635
## Congresox_talabonoides 5.521461
## Cynoponticus_ferox 5.298317
## Dysomma_anguillare 3.951244
## Elops_saurus 4.605170
## Facciolella_gilbertii 4.110874
## Gavialiceps_taeniola 4.169761
## Gnathophis_longicauda 3.653252
## Gorgasia_taiwanensis 4.305416
## Gymnothorax_castaneus 5.010635
## Gymnothorax_flavimarginatus 5.480639
## Gymnothorax_kidako 4.516339
## Gymnothorax_moringa 5.298317
## Gymnothorax_pseudothyrsoideus 4.382027
## Gymnothorax_reticularis 4.094345
## Heteroconger_hassi 3.688879
## Ichthyapus_ophioneus 3.806662
## Kaupichthys_hyoproroides 3.401197
## Kaupichthys_nuchalis 2.772589
## Megalops_cyprinoides 5.010635
## Moringua_edwardsi 2.708050
## Moringua_javanica 4.787492
## Muraenesox_bagio 5.298317
## Muraenesox_cinereus 5.393628
## Myrichthys_breviceps 4.624973
## Myrichthys_maculosus 4.605170
## Myrichthys_magnificus 4.356709
## Myrophis_vafer 3.828641
## Nemichthys_scolopaceus 4.867534
## Nettastoma_melanurum 4.347694
## Ophichthus_serpentinus 4.219508
## Ophichthus_zophochir 4.477337
## Oxyconger_leptognathus 4.094345
## Parabathymyrus_macrophthalmus 3.850148
## Paraconger_notialis 4.138361
## Pisodonophis_cancrivorus 4.682131
## Poeciloconger_kapala 4.094345
## Rhinomuraena_quaesita 4.867534
## Rhynchoconger_flavus 5.010635
## Saurenchelys_fierasfer 3.912023
## Scolecenchelys_breviceps 4.094345
## Scuticaria_tigrina 4.787492
## Serrivomer_beanii 4.356709
## Serrivomer_sector 4.330733
## Simenchelys_parasitica 4.110874
## Uroconger_lepturus 3.951244
## Uropterygius_micropterus 3.401197
## Venefica_proboscidea 4.605170
We are already ready to run our analysis:
library(OUwie)
## Warning: package 'OUwie' was built under R version 3.3.1
## Loading required package: nloptr
## Warning: package 'nloptr' was built under R version 3.3.1
## Loading required package: lattice
fitOU<-list()
for(i in 1:length(trees))
fitOU[[i]]<-OUwie(trees[[i]],OUdata,model="OUM",simmap.tree=TRUE)
## Initializing...
## Finished. Begin thorough search...
## Finished. Summarizing results.
## Initializing...
## Finished. Begin thorough search...
## Finished. Summarizing results.
## Initializing...
## Finished. Begin thorough search...
## Finished. Summarizing results.
## Initializing...
## Finished. Begin thorough search...
## Finished. Summarizing results.
## Initializing...
## Finished. Begin thorough search...
## Finished. Summarizing results.
## Initializing...
## Finished. Begin thorough search...
## Finished. Summarizing results.
## Initializing...
## Finished. Begin thorough search...
## Finished. Summarizing results.
## Initializing...
## Finished. Begin thorough search...
## Finished. Summarizing results.
## Initializing...
## Finished. Begin thorough search...
## Finished. Summarizing results.
## Initializing...
## Finished. Begin thorough search...
## Finished. Summarizing results.
theta<-matrix(NA,length(trees),length(fitOU[[1]]$tot.states))
for(i in 1:length(trees))
theta[i,]<-fitOU[[i]]$theta[order(fitOU[[i]]$tot.states),1]
colnames(theta)<-sort(fitOU[[1]]$tot.states)
theta
## bite suction
## [1,] 4.475331 4.337566
## [2,] 4.345699 4.415305
## [3,] 4.589966 4.274586
## [4,] 4.240477 4.512560
## [5,] 4.277375 4.481104
## [6,] 4.404827 4.323435
## [7,] 4.645394 4.255197
## [8,] 4.351069 4.384847
## [9,] 4.465252 4.218027
## [10,] 4.448947 4.308154
Once again, in a “real” study we would extract much more from each fitted model than simply the MLEs of θ.
Written by Liam J. Revell. Last updated 30 Jun. 2016.