## Solution to Challenge Problem 7: Exploring multi-regime models for continuous traits

The challenge problem was as follows.

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?

``````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")
``````

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
``````

``````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.