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