Multi-regime challenge question


Multi-optimum OU evolution

Multi-optimum OU evolution

## load OUwie & phytools
library(OUwie)
## Loading required package: ape
## Loading required package: nloptr
library(phytools)
## Loading required package: maps
## read anole data from file
X<-read.csv("Anolis.csv",row.names=1,header=TRUE)
## read anole tree
tree<-read.tree("Anolis.tre")
## make analysis input data.frame
ecomorph<-as.matrix(X)[,"ecomorph"]
svl<-as.matrix(X)[,"svl"]
data<-data.frame(Genus_species=rownames(X),Reg=ecomorph,X=as.numeric(svl))
## perform stochastic mapping
trees<-make.simmap(tree,ecomorph,nsim=10)
## make.simmap is sampling character histories conditioned on the transition matrix
## Q =
##            CG         GB          TC         TG          Tr         Tw
## CG -0.3933357  0.0000000  0.19321489  0.0000000  0.00000000  0.2001209
## GB  0.0000000 -0.5024038  0.00000000  0.1959991  0.00000000  0.3064047
## TC  0.1932149  0.0000000 -0.74116782  0.0000000  0.09187982  0.4560731
## TG  0.0000000  0.1959991  0.00000000 -0.4428234  0.00000000  0.2468243
## Tr  0.0000000  0.0000000  0.09187982  0.0000000 -0.23442525  0.1425454
## Tw  0.2001209  0.3064047  0.45607311  0.2468243  0.14254542 -1.3519684
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##        CG        GB        TC        TG        Tr        Tw 
## 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667
## Done.
## hack
trees<-lapply(trees,function(x){ x$node.label<-getStates(x,"nodes"); x })
class(trees)<-"multiPhylo"
## fit OU
fitOU<-lapply(trees,OUwie,data=data,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.
## summarize the results across analyses
theta<-t(sapply(fitOU,function(x) x$theta[,1]))
colnames(theta)<-fitOU[[1]]$tot.states
theta
##             CG       GB       TC       TG       Tr       Tw
##  [1,] 5.195380 3.704090 4.138414 4.118540 3.817658 3.908048
##  [2,] 5.686878 3.614559 4.049840 4.133195 3.819009 3.979025
##  [3,] 5.263409 3.695316 4.152514 4.107550 3.835559 3.879574
##  [4,] 5.179145 3.663865 4.118413 4.117061 3.787933 3.935703
##  [5,] 5.812586 3.647810 4.082441 4.110418 3.816441 3.929778
##  [6,] 5.336393 3.679318 4.234161 4.110181 3.804458 3.969039
##  [7,] 5.232107 3.668695 4.144867 4.108663 3.760949 3.896545
##  [8,] 5.230976 3.657498 4.131510 4.105782 3.892274 3.867889
##  [9,] 5.583534 3.652226 4.101999 4.144991 3.821046 3.908917
## [10,] 5.539585 3.624411 4.165444 4.119165 3.816741 3.929854
colMeans(theta)
##       CG       GB       TC       TG       Tr       Tw 
## 5.405999 3.660779 4.131960 4.117555 3.817207 3.920437
apply(theta,2,var)
##           CG           GB           TC           TG           Tr 
## 0.0528268638 0.0007994982 0.0024971193 0.0001582489 0.0011305799 
##           Tw 
## 0.0012768616

Multivariate Brownian evolution

## load phytools
library(phytools)
## read centrarchid tree
tree<-read.tree("Centrarchidae.tre")
## read centrarchid data
X<-read.csv("Centrarchidae.csv",header=TRUE,row.names=1)
X
##                         feeding.mode gape.width buccal.length
## Acantharchus_pomotis            pisc      0.114        -0.009
## Lepomis_gibbosus                 non     -0.133        -0.009
## Lepomis_microlophus              non     -0.151         0.012
## Lepomis_punctatus                non     -0.103        -0.019
## Lepomis_miniatus                 non     -0.134         0.001
## Lepomis_auritus                  non     -0.222        -0.039
## Lepomis_marginatus               non     -0.187        -0.075
## Lepomis_megalotis                non     -0.073        -0.049
## Lepomis_humilis                  non      0.024        -0.027
## Lepomis_macrochirus              non     -0.191         0.002
## Lepomis_gulosus                 pisc      0.131         0.122
## Lepomis_symmetricus              non      0.013        -0.025
## Lepomis_cyanellus               pisc     -0.002        -0.009
## Micropterus_coosae              pisc      0.045        -0.009
## Micropterus_notius              pisc      0.097        -0.009
## Micropterus_treculi             pisc      0.056         0.001
## Micropterus_salmoides           pisc      0.056        -0.059
## Micropterus_floridanus          pisc      0.096         0.051
## Micropterus_punctulatus         pisc      0.092         0.002
## Micropterus_dolomieu            pisc      0.035        -0.069
## Centrarchus_macropterus          non     -0.007        -0.055
## Enneacantus_obesus               non      0.016        -0.005
## Pomoxis_annularis               pisc     -0.004        -0.019
## Pomoxis_nigromaculatus          pisc      0.105         0.041
## Archolites_interruptus          pisc     -0.024         0.051
## Ambloplites_ariommus            pisc      0.135         0.123
## Ambloplites_rupestris           pisc      0.086         0.041
## Ambloplites_cavifrons           pisc      0.130         0.040
## stochastic mapping of feeding mode on the tree (we would normally do x100)
fmode<-as.matrix(X)[,1]
trees<-make.simmap(tree,fmode,nsim=100)
## make.simmap is sampling character histories conditioned on the transition matrix
## Q =
##            non      pisc
## non  -4.228856  4.228856
## pisc  4.228856 -4.228856
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##  non pisc 
##  0.5  0.5
## Done.
## data
X<-as.matrix(X[,2:3])
## fit model
fit<-lapply(trees,function(x,y) evol.vcv(x,y),y=X)
r<-t(sapply(fit,function(x) 
c(cov2cor(x$R.multiple[,,"non"])[1,2],cov2cor(x$R.multiple[,,"pisc"])[1,2])))
colnames(r)<-c("r(non)","r(pisc)")
r
##              r(non)   r(pisc)
##   [1,] -0.003329626 0.7570405
##   [2,] -0.190872207 0.7557320
##   [3,] -0.086503148 0.7761596
##   [4,] -0.142029527 0.7606595
##   [5,] -0.056392012 0.7778228
##   [6,] -0.087586017 0.7687934
##   [7,] -0.156202349 0.7582008
##   [8,] -0.050055310 0.8057935
##   [9,] -0.105622360 0.7979087
##  [10,]  0.204348775 0.7161280
##  [11,]  0.127016343 0.7598835
##  [12,]  0.175557385 0.7263260
##  [13,] -0.125129232 0.7581984
##  [14,] -0.123772382 0.7608000
##  [15,] -0.128673602 0.7655795
##  [16,] -0.081680488 0.7608447
##  [17,] -0.034678116 0.7625138
##  [18,] -0.109736134 0.7656513
##  [19,]  0.264963264 0.7140560
##  [20,] -0.042744788 0.8350070
##  [21,] -0.100281745 0.7642891
##  [22,] -0.102982383 0.7674248
##  [23,] -0.026692827 0.8071230
##  [24,] -0.083104479 0.7735995
##  [25,] -0.103250588 0.7666826
##  [26,] -0.155412969 0.7558549
##  [27,]  0.082215179 0.7854529
##  [28,] -0.133140749 0.7645241
##  [29,] -0.053748832 0.7772855
##  [30,] -0.064307392 0.7595214
##  [31,] -0.014544488 0.7803895
##  [32,] -0.082179908 0.7682772
##  [33,]  0.302742593 0.6846339
##  [34,] -0.100570235 0.7598663
##  [35,]  0.190266696 0.8094365
##  [36,] -0.064397365 0.7656073
##  [37,]  0.143464311 0.7475610
##  [38,] -0.094360147 0.7593885
##  [39,]  0.066206612 0.7607164
##  [40,] -0.114894182 0.7616189
##  [41,]  0.090723594 0.7987334
##  [42,] -0.094498800 0.7552973
##  [43,] -0.134277649 0.7534743
##  [44,] -0.117715158 0.7772095
##  [45,] -0.115284787 0.7668147
##  [46,]  0.264354271 0.7698039
##  [47,] -0.093793693 0.7809012
##  [48,]  0.297341767 0.6766665
##  [49,] -0.103536031 0.7609834
##  [50,] -0.081802967 0.7619563
##  [51,] -0.034187281 0.7595237
##  [52,] -0.088139788 0.7608838
##  [53,] -0.043452787 0.7704884
##  [54,] -0.059796429 0.7641235
##  [55,]  0.298766506 0.6974202
##  [56,]  0.039348127 0.8755706
##  [57,] -0.150255539 0.7545902
##  [58,] -0.121486033 0.7888943
##  [59,] -0.172333681 0.7633727
##  [60,] -0.101178274 0.7699019
##  [61,] -0.091594307 0.7783453
##  [62,]  0.261420238 0.7080909
##  [63,] -0.158991460 0.7504172
##  [64,] -0.040127768 0.7582792
##  [65,] -0.115613907 0.8366310
##  [66,] -0.086048867 0.7942098
##  [67,] -0.190773421 0.8176592
##  [68,] -0.070912153 0.7908942
##  [69,] -0.125485187 0.7794628
##  [70,] -0.108208560 0.8116884
##  [71,] -0.083306025 0.7686787
##  [72,] -0.148114867 0.7561150
##  [73,] -0.224080980 0.7318766
##  [74,] -0.082629730 0.7635915
##  [75,] -0.066203839 0.7800408
##  [76,] -0.052850847 0.7695926
##  [77,]  0.259644203 0.7063583
##  [78,]  0.275135622 0.7926476
##  [79,] -0.059735530 0.7671309
##  [80,] -0.339474619 0.8503920
##  [81,] -0.056491863 0.7635030
##  [82,] -0.059277307 0.8378233
##  [83,] -0.113294212 0.7646433
##  [84,] -0.051694932 0.7589224
##  [85,]  0.196545842 0.7256409
##  [86,] -0.013303623 0.8583480
##  [87,] -0.063802277 0.7681354
##  [88,] -0.218783311 0.7698223
##  [89,] -0.098294084 0.8236113
##  [90,]  0.046482958 0.8078842
##  [91,] -0.134307758 0.7556345
##  [92,] -0.206182723 0.7443465
##  [93,] -0.158014739 0.7636064
##  [94,] -0.070572659 0.7930967
##  [95,] -0.110957671 0.7653692
##  [96,] -0.018472734 0.7479191
##  [97,] -0.119710489 0.7776119
##  [98,]  0.009240777 0.7770134
##  [99,]  0.260074614 0.7113830
## [100,] -0.104631559 0.7652208
colMeans(r)
##      r(non)     r(pisc) 
## -0.04116699  0.76900600
apply(r,2,var)
##      r(non)     r(pisc) 
## 0.017788245 0.001071565

Written by Liam J. Revell. Last updated July 4, 2015