Trait evolution using mvMORPH
Christoph Liedtke
August, 2019
Intro
Recently, we wanted to test whether a specific trait evolves under different evolutionary models for different subclades of a tree. In our case genome size in amphibians. To fit Brownian motion and Ornstein-Uhlenbeck models, we used Julien Clavel’s mvMORPH R package.
Tree: pruned amphibian tree originally from Pyron 2014
Data: genome sizes (c-values) for amphibian species from Liedtke et al. 2018
Loading Data and Libraries
library(mvMORPH)
# install from CRAN or most recent release from github: https://github.com/JClavel/mvMORPH
# make sure dependencies are installed too!### load data:
amphibia_cval<-read.csv("./data/amphibia_cval.csv", row.names = 1, sep=";")
### load tree:
amphibia.tree<-read.tree("./data/amphibia.tre")
### match tree and data order
amphibia_cval<-amphibia_cval[amphibia.tree$tip.label,]
all(rownames(amphibia_cval)==amphibia.tree$tip.label)## [1] TRUE### log transform variable to be normally distributed
amphibia_logC<-log10(amphibia_cval$median_c); names(amphibia_logC)<-rownames(amphibia_cval) # log10 transform data
### plot phenogram to vizualize trait distribution
phenogram(amphibia.tree, amphibia_logC, ftype = "off", ylab="genome size [log10(c-value)]")
Setting up hypotheses
By looking at the phylogenetic distribution of the data, it seemed like salamanders were off doing their own thing with really large genomes (that isolated clade, all with large genomes) we wanted to test whether a single trait model for the whole tree is really best in this case, or whether fitting a multiple optimum/parameters model would be more just. This way we could also test the hypothesis of whether each of the three amphibian orders (Frogs, Salamanders and Caecilians) have evolved under different sets of parameters.
To set up these hypotheses, we created a simmap object using (phytools)[http://blog.phytools.org/], a dependency of mvMORPH, ‘painting’ (i.e. annotating) subclades of the tree that we want to model different processes for.
# make a colour scheme to use for vizualizing clades (in this case the three amphibian orders)
cols<-c("black","deepskyblue3","chartreuse3","gold"); names(cols)<-c(1,"caudata","anura","gymno")
# make a simmap object with two 'regimes', comparing salamanders to non-salamanders
two_regimes_tree<-paintSubTree(amphibia.tree,
node=getMRCA(amphibia.tree, rownames(amphibia_cval)[amphibia_cval$order=="Urodela"]),
state="caudata",
stem=T)
plotSimmap(two_regimes_tree,cols,lwd=2,pts=F, ftype="off")
# make a simmap object with three 'regimes' where each of the three amphibian orders is treated separately
three_regimes_tree<-paintSubTree(amphibia.tree,
node=getMRCA(amphibia.tree, rownames(amphibia_cval)[amphibia_cval$order=="Gymnophiona"]),
state="gymno",
stem=T,
anc.state = "gymno")
three_regimes_tree<-paintSubTree(three_regimes_tree,
node=getMRCA(amphibia.tree,rownames(amphibia_cval)[amphibia_cval$order=="Urodela"]),
state="caudata",
stem=T)
three_regimes_tree<-paintSubTree(three_regimes_tree,
node=getMRCA(amphibia.tree,rownames(amphibia_cval)[amphibia_cval$order=="Anura"]),
state="anura",
stem=T)
plotSimmap(three_regimes_tree,cols,lwd=2,pts=F, ftype="off")
Fitting models
With the clades of interest defined, we can now fit BM and OU models. note: as of late, mvMORPH also allows fitting early burst and shift models too!
##
## convergence of the optimizer has not been reached, try simpler model
## a reliable solution has been reached## successful convergence of the optimizer
## a reliable solution has been reached# two optima/parameter sets:
fit_bm2<-mvBM(two_regimes_tree,amphibia_logC,model="BMM", param=list(smean=F),echo = F)## successful convergence of the optimizer
## a reliable solution has been reached## successful convergence of the optimizer
## a reliable solution has been reached# three optima/parameter sets:
fit_bm3<-mvBM(three_regimes_tree,amphibia_logC,model="BMM", param=list(smean=F),echo = T)## successful convergence of the optimizer
## a reliable solution has been reached
##
## -- Summary results for multiple rate BMM model --
## LogLikelihood: 288.415
## AIC: -564.8299
## AICc: -564.6461
## 6 parameters
##
## Estimated rate matrix
## ______________________
## , , gymno
##
##
## 0.0003795316
##
## , , caudata
##
##
## 0.0002565381
##
## , , anura
##
##
## 0.0005345463
##
##
## Estimated root state
## ______________________
##
## gymno 0.6168402
## caudata 1.6304365
## anura 0.7741864## successful convergence of the optimizer
## a reliable solution has been reachedCompare model fit
Various things can be done here, for example we can calculate Akaike weights to compare models or we can perform likelihood ratio tests for each type of models (BM or OU)
## compare model fit using Akaike weights:
results<-list(fit_bm1,fit_ou1, fit_bm2, fit_ou2,fit_bm3,fit_ou3)
results<-aicw(results, aicc=TRUE)
results # note, the model names used can be checked in each of the model fit objects.## -- Akaike weights --
## Rank AIC diff wi AICw
## BMM default 3 1 -567 0.00 1.00e+00 8.03e-01
## BMM default 5 2 -565 2.81 2.46e-01 1.97e-01
## OUM 4 3 -547 20.75 3.13e-05 2.51e-05
## OUM 6 4 -546 21.85 1.80e-05 1.45e-05
## BM1 default 1 5 -530 37.23 0.00e+00 0.00e+00
## OU1 2 6 -528 39.26 0.00e+00 0.00e+00## Warning in LRT(fit_bm1, fit_bm2): The order of the models has been changed
## for the comparison## -- Log-likelihood Ratio Test --
## Model BMM versus BM1
## Number of degrees of freedom : 2
## LRT statistic: 41.29293 p-value: 1.079829e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Warning in LRT(fit_bm1, fit_bm3): The order of the models has been changed
## for the comparison## -- Log-likelihood Ratio Test --
## Model BMM versus BM1
## Number of degrees of freedom : 4
## LRT statistic: 42.58177 p-value: 1.263614e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Warning in LRT(fit_bm2, fit_bm3): The order of the models has been changed
## for the comparison## -- Log-likelihood Ratio Test --
## Model BMM versus BMM
## Number of degrees of freedom : 2
## LRT statistic: 1.288839 p-value: 0.5249671
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Warning in LRT(fit_ou1, fit_ou2): The order of the models has been changed
## for the comparison## -- Log-likelihood Ratio Test --
## Model OUM symmetric positive versus OU1 symmetric positive
## Number of degrees of freedom : 1
## LRT statistic: 20.55006 p-value: 5.809192e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Warning in LRT(fit_ou1, fit_ou3): The order of the models has been changed
## for the comparison## -- Log-likelihood Ratio Test --
## Model OUM symmetric positive versus OU1 symmetric positive
## Number of degrees of freedom : 2
## LRT statistic: 21.49263 p-value: 2.152453e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Warning in LRT(fit_ou2, fit_ou3): The order of the models has been changed
## for the comparison## -- Log-likelihood Ratio Test --
## Model OUM symmetric positive versus OUM symmetric positive
## Number of degrees of freedom : 1
## LRT statistic: 0.9425729 p-value: 0.331617
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The AIC comparison suggests that the best model was a Brownian motion model with multiple sets of parameters, one for salamanders and one for everything else (‘BMM default 3 model’ in the Akeike weights table), but the likelihood ratio test was less clear about this, suggesting that that the 2 and 3 parameter sets performed equally well (or at least not significantly better or worse).
Ancestral state reconstruction
Seeing as the 2 and 3 parameter sets performed best, we are going to reconstruct ancestral states with these. Although BM outperformed OU, just as an example, I will reconstruct ancestral states with both.
# BM (2 parameter model)
bm2.asr.fit<-mvBM(two_regimes_tree,amphibia_logC,model="BMM", param=list(smean=F)) ## successful convergence of the optimizer
## a reliable solution has been reached
##
## -- Summary results for multiple rate BMM model --
## LogLikelihood: 287.7705
## AIC: -567.5411
## AICc: -567.4539
## 4 parameters
##
## Estimated rate matrix
## ______________________
## , , 1
##
##
## 0.0005240248
##
## , , caudata
##
##
## 0.0002566532
##
##
## Estimated root state
## ______________________
##
## 1 0.6642491
## caudata 1.6354438bm2.asr.estim<-estim(two_regimes_tree, amphibia_logC, bm2.asr.fit, asr=TRUE)
# BM (3 parameter model)
bm3.asr.fit<-mvBM(three_regimes_tree,amphibia_logC,model="BMM", param=list(smean=F)) ## successful convergence of the optimizer
## a reliable solution has been reached
##
## -- Summary results for multiple rate BMM model --
## LogLikelihood: 288.415
## AIC: -564.8299
## AICc: -564.6461
## 6 parameters
##
## Estimated rate matrix
## ______________________
## , , gymno
##
##
## 0.0003795316
##
## , , caudata
##
##
## 0.0002565381
##
## , , anura
##
##
## 0.0005345463
##
##
## Estimated root state
## ______________________
##
## gymno 0.6168402
## caudata 1.6304365
## anura 0.7741864bm3.asr.estim<-estim(three_regimes_tree, amphibia_logC, bm3.asr.fit, asr=TRUE)
# OU (2 parameter model)
ou2.asr.fit<-mvOU(two_regimes_tree,amphibia_logC,model="OUM", param=list(root=F)) ## successful convergence of the optimizer
## a reliable solution has been reached
##
## -- Summary results --
## LogLikelihood: 277.3975
## AIC: -546.795
## AICc: -546.7078
## 4 parameters
##
## Estimated theta values
## ______________________
##
## 1 0.6399487
## caudata 2.2490503
##
## ML alpha values
## ______________________
##
## 0.003615548
##
## ML sigma values
## ______________________
##
## 0.000500618ou2.asr.estim<-estim(two_regimes_tree, amphibia_logC, ou2.asr.fit, asr=TRUE)
# OU (3 parameter model)
ou3.asr.fit<-mvOU(three_regimes_tree,amphibia_logC,model="OUM", param=list(root=F)) ## successful convergence of the optimizer
## a reliable solution has been reached
##
## -- Summary results --
## LogLikelihood: 277.8688
## AIC: -545.7375
## AICc: -545.6065
## 5 parameters
##
## Estimated theta values
## ______________________
##
## gymno 0.7781379
## caudata 2.0813825
## anura 0.4916251
##
## ML alpha values
## ______________________
##
## 0.003975036
##
## ML sigma values
## ______________________
##
## 0.0005099655We can now visualize these reconstructions using phenograms for example:
## BM
#2 means
bm2.asr.states<-10^c(amphibia_logC,bm2.asr.estim$estimates) # The key here is to paste the tip states together with the node states estimated by mvMORPH as one vector. Here I am converting log10() traits back to real space so that I can plot a y axis on a log scale. this is not necesary, but I think it is more easily interpreted this way.
names(bm2.asr.states)[names(bm2.asr.states)==""]<-which(names(bm2.asr.states)=="") # some nodes don't have node numbers, which confuses phenogram() so this is just to make sure they are all numbered nicely
par(mfrow=c(2,2))
par(mar=c(4,4,1,1))
phenogram(tree=two_regimes_tree,x=bm2.asr.states, ftype="off", ylim=c(1,120), log="y", colors=cols)## 1 2 3 4 5 6
## 4.471551 5.450384 11.798339 3.947989 6.607296 11.152378
## 7 8 9 10 11 12
## 4.233458 2.940315 9.667680 4.817965 5.467809 8.869577
## 13 14 15 16 17 18
## 11.027308 5.393683 8.342427 3.700000 7.894747 5.653214
## 19 20 21 22 23 24
## 7.356226 6.663574 5.933774 4.358541 55.000000 50.000000
## 25 26 27 28 29 30
## 46.490000 47.540000 24.620000 24.075000 21.650000 16.160000
## 31 32 33 34 35 36
## 20.450000 19.200000 16.880000 16.480000 48.000000 45.500000
## 37 38 39 40 41 42
## 54.000000 57.000000 56.760000 30.750000 42.000000 32.050000
## 43 44 45 46 47 48
## 29.500000 32.000000 26.965000 34.500000 25.800000 24.410000
## 49 50 51 52 53 54
## 26.170000 32.000000 27.600000 35.500000 20.600000 29.220000
## 55 56 57 58 59 60
## 39.300000 35.335000 32.550000 23.160000 27.300000 24.490000
## 61 62 63 64 65 66
## 38.000000 30.000000 35.100000 30.060000 28.790000 37.800000
## 67 68 69 70 71 72
## 43.700000 44.270000 47.250000 39.800000 43.090000 51.000000
## 73 74 75 76 77 78
## 50.480000 23.100000 26.340000 22.000000 30.845000 27.420000
## 79 80 81 82 83 84
## 28.890000 32.595000 25.670000 37.980000 34.970000 24.329944
## 85 86 87 88 89 90
## 26.960000 24.000000 27.240000 28.695000 24.875000 28.160000
## 91 92 93 94 95 96
## 30.955000 48.840000 120.565000 119.225000 82.550000 66.600000
## 97 98 99 100 101 102
## 79.745000 41.270000 50.000000 71.650000 42.630000 52.100000
## 103 104 105 106 107 108
## 43.190000 32.650000 42.900000 45.400000 42.400000 13.895000
## 109 110 111 112 113 114
## 15.060000 16.570000 15.150000 16.200000 15.130000 30.700000
## 115 116 117 118 119 120
## 48.550000 69.300000 38.800000 38.050000 33.630000 18.200000
## 121 122 123 124 125 126
## 22.640000 21.400000 20.650000 24.200000 22.600000 30.750000
## 127 128 129 130 131 132
## 33.700000 27.800000 28.000000 32.450000 23.250000 24.650000
## 133 134 135 136 137 138
## 22.230000 20.775000 27.350000 21.200000 25.000000 26.400000
## 139 140 141 142 143 144
## 24.500000 24.700000 26.200000 25.800000 26.500000 34.500000
## 145 146 147 148 149 150
## 33.200000 40.900000 28.080000 26.060000 34.930000 26.820000
## 151 152 153 154 155 156
## 29.680000 30.620000 31.520000 28.000000 36.900000 28.100000
## 157 158 159 160 161 162
## 40.900000 52.200000 44.900000 27.800000 28.500000 27.000000
## 163 164 165 166 167 168
## 20.800000 40.200000 46.600000 27.200000 30.500000 27.800000
## 169 170 171 172 173 174
## 40.300000 42.100000 36.900000 39.700000 52.400000 58.700000
## 175 176 177 178 179 180
## 48.700000 65.210000 53.900000 45.000000 42.300000 43.500000
## 181 182 183 184 185 186
## 42.000000 67.000000 62.700000 47.500000 47.800000 43.600000
## 187 188 189 190 191 192
## 51.000000 59.400000 52.700000 50.600000 50.200000 56.500000
## 193 194 195 196 197 198
## 4.160000 3.924195 8.630000 8.440000 8.400000 11.385000
## 199 200 201 202 203 204
## 9.600000 6.700000 9.780000 4.281789 11.495696 5.250000
## 205 206 207 208 209 210
## 5.300671 7.223399 4.790000 2.430000 4.012580 1.740000
## 211 212 213 214 215 216
## 3.685000 7.895000 3.560000 4.020000 6.095000 3.110000
## 217 218 219 220 221 222
## 3.150000 3.110000 1.644467 1.400000 1.235000 1.445000
## 223 224 225 226 227 228
## 2.270000 2.440000 3.960000 4.470000 3.956593 1.827083
## 229 230 231 232 233 234
## 2.570000 1.635449 3.900547 2.490000 7.008320 4.777544
## 235 236 237 238 239 240
## 4.708416 5.598801 3.097822 3.154025 2.722458 3.370000
## 241 242 243 244 245 246
## 4.470000 5.570000 4.066651 4.825000 4.250000 5.550000
## 247 248 249 250 251 252
## 3.380000 3.810000 4.735000 5.100000 1.540000 3.715000
## 253 254 255 256 257 258
## 4.870000 4.750000 5.450000 6.300000 5.750000 5.050000
## 259 260 261 262 263 264
## 2.000000 5.150000 3.350000 1.692412 1.433676 3.073715
## 265 266 267 268 269 270
## 1.088463 1.400000 2.610000 3.095000 5.020000 4.150000
## 271 272 273 274 275 276
## 3.085000 6.430000 6.750000 5.950000 6.000000 11.360000
## 277 278 279 280 281 282
## 8.150000 3.330000 5.350000 4.900000 6.340000 7.045000
## 283 284 285 286 287 288
## 4.250000 8.160000 6.520000 7.170000 5.830000 6.370000
## 289 290 291 292 293 294
## 6.430000 7.500000 9.195000 9.190000 7.910000 5.300000
## 295 296 297 298 299 300
## 5.250000 4.350000 5.270000 6.180000 4.750000 4.275000
## 301 302 303 304 305 306
## 5.340000 5.075000 4.490000 7.670000 5.320000 6.800000
## 307 308 309 310 311 312
## 6.055000 4.234224 5.660000 7.100618 7.660000 7.170000
## 313 314 315 316 317 318
## 3.250000 3.050000 3.500000 4.200000 3.900000 4.750000
## 319 320 321 322 323 324
## 5.700000 5.650000 4.800000 5.500000 4.800000 5.300000
## 325 326 327 328 329 330
## 5.250000 4.650000 4.150000 4.700000 4.800000 4.200000
## 331 332 333 334 335 336
## 4.550000 1.636627 2.675303 2.693029 2.770000 2.615000
## 337 338 339 340 341 342
## 0.950000 2.100000 2.400000 2.250000 3.790000 1.950000
## 343 344 345 346 347 348
## 3.250000 2.120000 8.750000 3.930000 4.710000 2.410000
## 349 350 351 352 353 354
## 3.570000 4.190000 3.300000 4.560000 2.310000 2.920000
## 355 356 357 358 359 360
## 3.140000 2.900000 3.301585 4.102954 3.350000 1.580000
## 361 362 363 364 365 366
## 1.622500 2.155000 1.580000 2.150000 2.420000 2.060633
## 367 368 369 370 371 372
## 2.948384 1.930000 2.589580 5.254850 6.750000 6.545000
## 373 374 375 376 377 378
## 4.790000 2.850000 4.530000 3.280000 3.760000 6.610000
## 379 380 381 382 383 384
## 2.630000 4.830000 2.730000 2.350000 3.710000 5.150000
## 385 386 387 388 389 390
## 2.280000 3.760000 3.000000 1.540000 2.880485 2.554085
## 391 392 393 394 395 396
## 2.850000 1.890000 3.350000 2.530000 1.950000 1.830000
## 397 398 399 400 401 402
## 2.400000 4.490000 4.115000 3.810000 3.820000 3.620000
## 403 404 405 406 407 408
## 4.600000 4.140000 4.050000 4.370000 4.787160 4.060000
## 409 410 411 412 413 414
## 4.760000 4.300000 6.450000 3.830000 4.875000 4.670000
## 415 416 417 418 419 420
## 4.640000 2.420000 2.560000 1.750000 1.860000 1.992604
## 421 422 423 424 425 426
## 1.510000 2.520000 3.775000 4.470000 2.650000 2.077421
## 427 428 429 430 431 432
## 2.374929 2.991761 2.980000 7.110000 8.490000 8.950000
## 433 434 435 436 437 438
## 4.688845 6.860000 3.150000 3.290000 4.760000 2.910000
## 439 440 441 442 443 444
## 5.300000 5.600000 5.300000 6.700000 3.320000 5.450000
## 445 446 447 448 449 450
## 5.290000 5.000000 5.650000 4.230000 4.600000 7.310000
## 451 452 453 454 455 456
## 3.310000 3.990000 5.700000 6.800000 6.600000 5.740000
## 457 458 459 460 461 462
## 4.260000 3.671180 4.830000 4.900000 3.610000 6.020000
## 463 464
## 7.220980 3.379117#3 means
bm3.asr.states<-10^c(amphibia_logC,bm3.asr.estim$estimates)
names(bm3.asr.states)[names(bm3.asr.states)==""]<-which(names(bm3.asr.states)=="")
par(mar=c(4,4,1,1))
phenogram(tree=three_regimes_tree,x=bm3.asr.states, ftype="off", ylim=c(1,120), log="y", colors=cols)## 1 2 3 4 5 6
## 4.471551 5.450384 11.798339 3.947989 6.607296 11.152378
## 7 8 9 10 11 12
## 4.233458 2.940315 9.667680 4.817965 5.467809 8.869577
## 13 14 15 16 17 18
## 11.027308 5.393683 8.342427 3.700000 7.894747 5.653214
## 19 20 21 22 23 24
## 7.356226 6.663574 5.933774 4.358541 55.000000 50.000000
## 25 26 27 28 29 30
## 46.490000 47.540000 24.620000 24.075000 21.650000 16.160000
## 31 32 33 34 35 36
## 20.450000 19.200000 16.880000 16.480000 48.000000 45.500000
## 37 38 39 40 41 42
## 54.000000 57.000000 56.760000 30.750000 42.000000 32.050000
## 43 44 45 46 47 48
## 29.500000 32.000000 26.965000 34.500000 25.800000 24.410000
## 49 50 51 52 53 54
## 26.170000 32.000000 27.600000 35.500000 20.600000 29.220000
## 55 56 57 58 59 60
## 39.300000 35.335000 32.550000 23.160000 27.300000 24.490000
## 61 62 63 64 65 66
## 38.000000 30.000000 35.100000 30.060000 28.790000 37.800000
## 67 68 69 70 71 72
## 43.700000 44.270000 47.250000 39.800000 43.090000 51.000000
## 73 74 75 76 77 78
## 50.480000 23.100000 26.340000 22.000000 30.845000 27.420000
## 79 80 81 82 83 84
## 28.890000 32.595000 25.670000 37.980000 34.970000 24.329944
## 85 86 87 88 89 90
## 26.960000 24.000000 27.240000 28.695000 24.875000 28.160000
## 91 92 93 94 95 96
## 30.955000 48.840000 120.565000 119.225000 82.550000 66.600000
## 97 98 99 100 101 102
## 79.745000 41.270000 50.000000 71.650000 42.630000 52.100000
## 103 104 105 106 107 108
## 43.190000 32.650000 42.900000 45.400000 42.400000 13.895000
## 109 110 111 112 113 114
## 15.060000 16.570000 15.150000 16.200000 15.130000 30.700000
## 115 116 117 118 119 120
## 48.550000 69.300000 38.800000 38.050000 33.630000 18.200000
## 121 122 123 124 125 126
## 22.640000 21.400000 20.650000 24.200000 22.600000 30.750000
## 127 128 129 130 131 132
## 33.700000 27.800000 28.000000 32.450000 23.250000 24.650000
## 133 134 135 136 137 138
## 22.230000 20.775000 27.350000 21.200000 25.000000 26.400000
## 139 140 141 142 143 144
## 24.500000 24.700000 26.200000 25.800000 26.500000 34.500000
## 145 146 147 148 149 150
## 33.200000 40.900000 28.080000 26.060000 34.930000 26.820000
## 151 152 153 154 155 156
## 29.680000 30.620000 31.520000 28.000000 36.900000 28.100000
## 157 158 159 160 161 162
## 40.900000 52.200000 44.900000 27.800000 28.500000 27.000000
## 163 164 165 166 167 168
## 20.800000 40.200000 46.600000 27.200000 30.500000 27.800000
## 169 170 171 172 173 174
## 40.300000 42.100000 36.900000 39.700000 52.400000 58.700000
## 175 176 177 178 179 180
## 48.700000 65.210000 53.900000 45.000000 42.300000 43.500000
## 181 182 183 184 185 186
## 42.000000 67.000000 62.700000 47.500000 47.800000 43.600000
## 187 188 189 190 191 192
## 51.000000 59.400000 52.700000 50.600000 50.200000 56.500000
## 193 194 195 196 197 198
## 4.160000 3.924195 8.630000 8.440000 8.400000 11.385000
## 199 200 201 202 203 204
## 9.600000 6.700000 9.780000 4.281789 11.495696 5.250000
## 205 206 207 208 209 210
## 5.300671 7.223399 4.790000 2.430000 4.012580 1.740000
## 211 212 213 214 215 216
## 3.685000 7.895000 3.560000 4.020000 6.095000 3.110000
## 217 218 219 220 221 222
## 3.150000 3.110000 1.644467 1.400000 1.235000 1.445000
## 223 224 225 226 227 228
## 2.270000 2.440000 3.960000 4.470000 3.956593 1.827083
## 229 230 231 232 233 234
## 2.570000 1.635449 3.900547 2.490000 7.008320 4.777544
## 235 236 237 238 239 240
## 4.708416 5.598801 3.097822 3.154025 2.722458 3.370000
## 241 242 243 244 245 246
## 4.470000 5.570000 4.066651 4.825000 4.250000 5.550000
## 247 248 249 250 251 252
## 3.380000 3.810000 4.735000 5.100000 1.540000 3.715000
## 253 254 255 256 257 258
## 4.870000 4.750000 5.450000 6.300000 5.750000 5.050000
## 259 260 261 262 263 264
## 2.000000 5.150000 3.350000 1.692412 1.433676 3.073715
## 265 266 267 268 269 270
## 1.088463 1.400000 2.610000 3.095000 5.020000 4.150000
## 271 272 273 274 275 276
## 3.085000 6.430000 6.750000 5.950000 6.000000 11.360000
## 277 278 279 280 281 282
## 8.150000 3.330000 5.350000 4.900000 6.340000 7.045000
## 283 284 285 286 287 288
## 4.250000 8.160000 6.520000 7.170000 5.830000 6.370000
## 289 290 291 292 293 294
## 6.430000 7.500000 9.195000 9.190000 7.910000 5.300000
## 295 296 297 298 299 300
## 5.250000 4.350000 5.270000 6.180000 4.750000 4.275000
## 301 302 303 304 305 306
## 5.340000 5.075000 4.490000 7.670000 5.320000 6.800000
## 307 308 309 310 311 312
## 6.055000 4.234224 5.660000 7.100618 7.660000 7.170000
## 313 314 315 316 317 318
## 3.250000 3.050000 3.500000 4.200000 3.900000 4.750000
## 319 320 321 322 323 324
## 5.700000 5.650000 4.800000 5.500000 4.800000 5.300000
## 325 326 327 328 329 330
## 5.250000 4.650000 4.150000 4.700000 4.800000 4.200000
## 331 332 333 334 335 336
## 4.550000 1.636627 2.675303 2.693029 2.770000 2.615000
## 337 338 339 340 341 342
## 0.950000 2.100000 2.400000 2.250000 3.790000 1.950000
## 343 344 345 346 347 348
## 3.250000 2.120000 8.750000 3.930000 4.710000 2.410000
## 349 350 351 352 353 354
## 3.570000 4.190000 3.300000 4.560000 2.310000 2.920000
## 355 356 357 358 359 360
## 3.140000 2.900000 3.301585 4.102954 3.350000 1.580000
## 361 362 363 364 365 366
## 1.622500 2.155000 1.580000 2.150000 2.420000 2.060633
## 367 368 369 370 371 372
## 2.948384 1.930000 2.589580 5.254850 6.750000 6.545000
## 373 374 375 376 377 378
## 4.790000 2.850000 4.530000 3.280000 3.760000 6.610000
## 379 380 381 382 383 384
## 2.630000 4.830000 2.730000 2.350000 3.710000 5.150000
## 385 386 387 388 389 390
## 2.280000 3.760000 3.000000 1.540000 2.880485 2.554085
## 391 392 393 394 395 396
## 2.850000 1.890000 3.350000 2.530000 1.950000 1.830000
## 397 398 399 400 401 402
## 2.400000 4.490000 4.115000 3.810000 3.820000 3.620000
## 403 404 405 406 407 408
## 4.600000 4.140000 4.050000 4.370000 4.787160 4.060000
## 409 410 411 412 413 414
## 4.760000 4.300000 6.450000 3.830000 4.875000 4.670000
## 415 416 417 418 419 420
## 4.640000 2.420000 2.560000 1.750000 1.860000 1.992604
## 421 422 423 424 425 426
## 1.510000 2.520000 3.775000 4.470000 2.650000 2.077421
## 427 428 429 430 431 432
## 2.374929 2.991761 2.980000 7.110000 8.490000 8.950000
## 433 434 435 436 437 438
## 4.688845 6.860000 3.150000 3.290000 4.760000 2.910000
## 439 440 441 442 443 444
## 5.300000 5.600000 5.300000 6.700000 3.320000 5.450000
## 445 446 447 448 449 450
## 5.290000 5.000000 5.650000 4.230000 4.600000 7.310000
## 451 452 453 454 455 456
## 3.310000 3.990000 5.700000 6.800000 6.600000 5.740000
## 457 458 459 460 461 462
## 4.260000 3.671180 4.830000 4.900000 3.610000 6.020000
## 463 464
## 7.220980 3.379117## OU
#2 optima
ou2.asr.states<-10^c(amphibia_logC,ou2.asr.estim$estimates)
names(ou2.asr.states)[names(ou2.asr.states)==""]<-which(names(ou2.asr.states)=="")
par(mar=c(4,4,1,1))
phenogram(tree=two_regimes_tree,x=ou2.asr.states, fsize=0.0001, spread.labels=F, ylim=c(1,120), log="y", colors=cols, ftype="reg")## 1 2 3 4 5 6
## 4.471551 5.450384 11.798339 3.947989 6.607296 11.152378
## 7 8 9 10 11 12
## 4.233458 2.940315 9.667680 4.817965 5.467809 8.869577
## 13 14 15 16 17 18
## 11.027308 5.393683 8.342427 3.700000 7.894747 5.653214
## 19 20 21 22 23 24
## 7.356226 6.663574 5.933774 4.358541 55.000000 50.000000
## 25 26 27 28 29 30
## 46.490000 47.540000 24.620000 24.075000 21.650000 16.160000
## 31 32 33 34 35 36
## 20.450000 19.200000 16.880000 16.480000 48.000000 45.500000
## 37 38 39 40 41 42
## 54.000000 57.000000 56.760000 30.750000 42.000000 32.050000
## 43 44 45 46 47 48
## 29.500000 32.000000 26.965000 34.500000 25.800000 24.410000
## 49 50 51 52 53 54
## 26.170000 32.000000 27.600000 35.500000 20.600000 29.220000
## 55 56 57 58 59 60
## 39.300000 35.335000 32.550000 23.160000 27.300000 24.490000
## 61 62 63 64 65 66
## 38.000000 30.000000 35.100000 30.060000 28.790000 37.800000
## 67 68 69 70 71 72
## 43.700000 44.270000 47.250000 39.800000 43.090000 51.000000
## 73 74 75 76 77 78
## 50.480000 23.100000 26.340000 22.000000 30.845000 27.420000
## 79 80 81 82 83 84
## 28.890000 32.595000 25.670000 37.980000 34.970000 24.329944
## 85 86 87 88 89 90
## 26.960000 24.000000 27.240000 28.695000 24.875000 28.160000
## 91 92 93 94 95 96
## 30.955000 48.840000 120.565000 119.225000 82.550000 66.600000
## 97 98 99 100 101 102
## 79.745000 41.270000 50.000000 71.650000 42.630000 52.100000
## 103 104 105 106 107 108
## 43.190000 32.650000 42.900000 45.400000 42.400000 13.895000
## 109 110 111 112 113 114
## 15.060000 16.570000 15.150000 16.200000 15.130000 30.700000
## 115 116 117 118 119 120
## 48.550000 69.300000 38.800000 38.050000 33.630000 18.200000
## 121 122 123 124 125 126
## 22.640000 21.400000 20.650000 24.200000 22.600000 30.750000
## 127 128 129 130 131 132
## 33.700000 27.800000 28.000000 32.450000 23.250000 24.650000
## 133 134 135 136 137 138
## 22.230000 20.775000 27.350000 21.200000 25.000000 26.400000
## 139 140 141 142 143 144
## 24.500000 24.700000 26.200000 25.800000 26.500000 34.500000
## 145 146 147 148 149 150
## 33.200000 40.900000 28.080000 26.060000 34.930000 26.820000
## 151 152 153 154 155 156
## 29.680000 30.620000 31.520000 28.000000 36.900000 28.100000
## 157 158 159 160 161 162
## 40.900000 52.200000 44.900000 27.800000 28.500000 27.000000
## 163 164 165 166 167 168
## 20.800000 40.200000 46.600000 27.200000 30.500000 27.800000
## 169 170 171 172 173 174
## 40.300000 42.100000 36.900000 39.700000 52.400000 58.700000
## 175 176 177 178 179 180
## 48.700000 65.210000 53.900000 45.000000 42.300000 43.500000
## 181 182 183 184 185 186
## 42.000000 67.000000 62.700000 47.500000 47.800000 43.600000
## 187 188 189 190 191 192
## 51.000000 59.400000 52.700000 50.600000 50.200000 56.500000
## 193 194 195 196 197 198
## 4.160000 3.924195 8.630000 8.440000 8.400000 11.385000
## 199 200 201 202 203 204
## 9.600000 6.700000 9.780000 4.281789 11.495696 5.250000
## 205 206 207 208 209 210
## 5.300671 7.223399 4.790000 2.430000 4.012580 1.740000
## 211 212 213 214 215 216
## 3.685000 7.895000 3.560000 4.020000 6.095000 3.110000
## 217 218 219 220 221 222
## 3.150000 3.110000 1.644467 1.400000 1.235000 1.445000
## 223 224 225 226 227 228
## 2.270000 2.440000 3.960000 4.470000 3.956593 1.827083
## 229 230 231 232 233 234
## 2.570000 1.635449 3.900547 2.490000 7.008320 4.777544
## 235 236 237 238 239 240
## 4.708416 5.598801 3.097822 3.154025 2.722458 3.370000
## 241 242 243 244 245 246
## 4.470000 5.570000 4.066651 4.825000 4.250000 5.550000
## 247 248 249 250 251 252
## 3.380000 3.810000 4.735000 5.100000 1.540000 3.715000
## 253 254 255 256 257 258
## 4.870000 4.750000 5.450000 6.300000 5.750000 5.050000
## 259 260 261 262 263 264
## 2.000000 5.150000 3.350000 1.692412 1.433676 3.073715
## 265 266 267 268 269 270
## 1.088463 1.400000 2.610000 3.095000 5.020000 4.150000
## 271 272 273 274 275 276
## 3.085000 6.430000 6.750000 5.950000 6.000000 11.360000
## 277 278 279 280 281 282
## 8.150000 3.330000 5.350000 4.900000 6.340000 7.045000
## 283 284 285 286 287 288
## 4.250000 8.160000 6.520000 7.170000 5.830000 6.370000
## 289 290 291 292 293 294
## 6.430000 7.500000 9.195000 9.190000 7.910000 5.300000
## 295 296 297 298 299 300
## 5.250000 4.350000 5.270000 6.180000 4.750000 4.275000
## 301 302 303 304 305 306
## 5.340000 5.075000 4.490000 7.670000 5.320000 6.800000
## 307 308 309 310 311 312
## 6.055000 4.234224 5.660000 7.100618 7.660000 7.170000
## 313 314 315 316 317 318
## 3.250000 3.050000 3.500000 4.200000 3.900000 4.750000
## 319 320 321 322 323 324
## 5.700000 5.650000 4.800000 5.500000 4.800000 5.300000
## 325 326 327 328 329 330
## 5.250000 4.650000 4.150000 4.700000 4.800000 4.200000
## 331 332 333 334 335 336
## 4.550000 1.636627 2.675303 2.693029 2.770000 2.615000
## 337 338 339 340 341 342
## 0.950000 2.100000 2.400000 2.250000 3.790000 1.950000
## 343 344 345 346 347 348
## 3.250000 2.120000 8.750000 3.930000 4.710000 2.410000
## 349 350 351 352 353 354
## 3.570000 4.190000 3.300000 4.560000 2.310000 2.920000
## 355 356 357 358 359 360
## 3.140000 2.900000 3.301585 4.102954 3.350000 1.580000
## 361 362 363 364 365 366
## 1.622500 2.155000 1.580000 2.150000 2.420000 2.060633
## 367 368 369 370 371 372
## 2.948384 1.930000 2.589580 5.254850 6.750000 6.545000
## 373 374 375 376 377 378
## 4.790000 2.850000 4.530000 3.280000 3.760000 6.610000
## 379 380 381 382 383 384
## 2.630000 4.830000 2.730000 2.350000 3.710000 5.150000
## 385 386 387 388 389 390
## 2.280000 3.760000 3.000000 1.540000 2.880485 2.554085
## 391 392 393 394 395 396
## 2.850000 1.890000 3.350000 2.530000 1.950000 1.830000
## 397 398 399 400 401 402
## 2.400000 4.490000 4.115000 3.810000 3.820000 3.620000
## 403 404 405 406 407 408
## 4.600000 4.140000 4.050000 4.370000 4.787160 4.060000
## 409 410 411 412 413 414
## 4.760000 4.300000 6.450000 3.830000 4.875000 4.670000
## 415 416 417 418 419 420
## 4.640000 2.420000 2.560000 1.750000 1.860000 1.992604
## 421 422 423 424 425 426
## 1.510000 2.520000 3.775000 4.470000 2.650000 2.077421
## 427 428 429 430 431 432
## 2.374929 2.991761 2.980000 7.110000 8.490000 8.950000
## 433 434 435 436 437 438
## 4.688845 6.860000 3.150000 3.290000 4.760000 2.910000
## 439 440 441 442 443 444
## 5.300000 5.600000 5.300000 6.700000 3.320000 5.450000
## 445 446 447 448 449 450
## 5.290000 5.000000 5.650000 4.230000 4.600000 7.310000
## 451 452 453 454 455 456
## 3.310000 3.990000 5.700000 6.800000 6.600000 5.740000
## 457 458 459 460 461 462
## 4.260000 3.671180 4.830000 4.900000 3.610000 6.020000
## 463 464
## 7.220980 3.379117#3 optima
ou3.asr.states<-10^c(amphibia_logC,ou3.asr.estim$estimates)
names(ou3.asr.states)[names(ou3.asr.states)==""]<-which(names(ou3.asr.states)=="")
par(mar=c(4,4,1,1))
phenogram(tree=three_regimes_tree,x=ou3.asr.states, fsize=0.0001, spread.labels=F, ylim=c(1,120), log="y", colors=cols, ftype="reg")
## 1 2 3 4 5 6
## 4.471551 5.450384 11.798339 3.947989 6.607296 11.152378
## 7 8 9 10 11 12
## 4.233458 2.940315 9.667680 4.817965 5.467809 8.869577
## 13 14 15 16 17 18
## 11.027308 5.393683 8.342427 3.700000 7.894747 5.653214
## 19 20 21 22 23 24
## 7.356226 6.663574 5.933774 4.358541 55.000000 50.000000
## 25 26 27 28 29 30
## 46.490000 47.540000 24.620000 24.075000 21.650000 16.160000
## 31 32 33 34 35 36
## 20.450000 19.200000 16.880000 16.480000 48.000000 45.500000
## 37 38 39 40 41 42
## 54.000000 57.000000 56.760000 30.750000 42.000000 32.050000
## 43 44 45 46 47 48
## 29.500000 32.000000 26.965000 34.500000 25.800000 24.410000
## 49 50 51 52 53 54
## 26.170000 32.000000 27.600000 35.500000 20.600000 29.220000
## 55 56 57 58 59 60
## 39.300000 35.335000 32.550000 23.160000 27.300000 24.490000
## 61 62 63 64 65 66
## 38.000000 30.000000 35.100000 30.060000 28.790000 37.800000
## 67 68 69 70 71 72
## 43.700000 44.270000 47.250000 39.800000 43.090000 51.000000
## 73 74 75 76 77 78
## 50.480000 23.100000 26.340000 22.000000 30.845000 27.420000
## 79 80 81 82 83 84
## 28.890000 32.595000 25.670000 37.980000 34.970000 24.329944
## 85 86 87 88 89 90
## 26.960000 24.000000 27.240000 28.695000 24.875000 28.160000
## 91 92 93 94 95 96
## 30.955000 48.840000 120.565000 119.225000 82.550000 66.600000
## 97 98 99 100 101 102
## 79.745000 41.270000 50.000000 71.650000 42.630000 52.100000
## 103 104 105 106 107 108
## 43.190000 32.650000 42.900000 45.400000 42.400000 13.895000
## 109 110 111 112 113 114
## 15.060000 16.570000 15.150000 16.200000 15.130000 30.700000
## 115 116 117 118 119 120
## 48.550000 69.300000 38.800000 38.050000 33.630000 18.200000
## 121 122 123 124 125 126
## 22.640000 21.400000 20.650000 24.200000 22.600000 30.750000
## 127 128 129 130 131 132
## 33.700000 27.800000 28.000000 32.450000 23.250000 24.650000
## 133 134 135 136 137 138
## 22.230000 20.775000 27.350000 21.200000 25.000000 26.400000
## 139 140 141 142 143 144
## 24.500000 24.700000 26.200000 25.800000 26.500000 34.500000
## 145 146 147 148 149 150
## 33.200000 40.900000 28.080000 26.060000 34.930000 26.820000
## 151 152 153 154 155 156
## 29.680000 30.620000 31.520000 28.000000 36.900000 28.100000
## 157 158 159 160 161 162
## 40.900000 52.200000 44.900000 27.800000 28.500000 27.000000
## 163 164 165 166 167 168
## 20.800000 40.200000 46.600000 27.200000 30.500000 27.800000
## 169 170 171 172 173 174
## 40.300000 42.100000 36.900000 39.700000 52.400000 58.700000
## 175 176 177 178 179 180
## 48.700000 65.210000 53.900000 45.000000 42.300000 43.500000
## 181 182 183 184 185 186
## 42.000000 67.000000 62.700000 47.500000 47.800000 43.600000
## 187 188 189 190 191 192
## 51.000000 59.400000 52.700000 50.600000 50.200000 56.500000
## 193 194 195 196 197 198
## 4.160000 3.924195 8.630000 8.440000 8.400000 11.385000
## 199 200 201 202 203 204
## 9.600000 6.700000 9.780000 4.281789 11.495696 5.250000
## 205 206 207 208 209 210
## 5.300671 7.223399 4.790000 2.430000 4.012580 1.740000
## 211 212 213 214 215 216
## 3.685000 7.895000 3.560000 4.020000 6.095000 3.110000
## 217 218 219 220 221 222
## 3.150000 3.110000 1.644467 1.400000 1.235000 1.445000
## 223 224 225 226 227 228
## 2.270000 2.440000 3.960000 4.470000 3.956593 1.827083
## 229 230 231 232 233 234
## 2.570000 1.635449 3.900547 2.490000 7.008320 4.777544
## 235 236 237 238 239 240
## 4.708416 5.598801 3.097822 3.154025 2.722458 3.370000
## 241 242 243 244 245 246
## 4.470000 5.570000 4.066651 4.825000 4.250000 5.550000
## 247 248 249 250 251 252
## 3.380000 3.810000 4.735000 5.100000 1.540000 3.715000
## 253 254 255 256 257 258
## 4.870000 4.750000 5.450000 6.300000 5.750000 5.050000
## 259 260 261 262 263 264
## 2.000000 5.150000 3.350000 1.692412 1.433676 3.073715
## 265 266 267 268 269 270
## 1.088463 1.400000 2.610000 3.095000 5.020000 4.150000
## 271 272 273 274 275 276
## 3.085000 6.430000 6.750000 5.950000 6.000000 11.360000
## 277 278 279 280 281 282
## 8.150000 3.330000 5.350000 4.900000 6.340000 7.045000
## 283 284 285 286 287 288
## 4.250000 8.160000 6.520000 7.170000 5.830000 6.370000
## 289 290 291 292 293 294
## 6.430000 7.500000 9.195000 9.190000 7.910000 5.300000
## 295 296 297 298 299 300
## 5.250000 4.350000 5.270000 6.180000 4.750000 4.275000
## 301 302 303 304 305 306
## 5.340000 5.075000 4.490000 7.670000 5.320000 6.800000
## 307 308 309 310 311 312
## 6.055000 4.234224 5.660000 7.100618 7.660000 7.170000
## 313 314 315 316 317 318
## 3.250000 3.050000 3.500000 4.200000 3.900000 4.750000
## 319 320 321 322 323 324
## 5.700000 5.650000 4.800000 5.500000 4.800000 5.300000
## 325 326 327 328 329 330
## 5.250000 4.650000 4.150000 4.700000 4.800000 4.200000
## 331 332 333 334 335 336
## 4.550000 1.636627 2.675303 2.693029 2.770000 2.615000
## 337 338 339 340 341 342
## 0.950000 2.100000 2.400000 2.250000 3.790000 1.950000
## 343 344 345 346 347 348
## 3.250000 2.120000 8.750000 3.930000 4.710000 2.410000
## 349 350 351 352 353 354
## 3.570000 4.190000 3.300000 4.560000 2.310000 2.920000
## 355 356 357 358 359 360
## 3.140000 2.900000 3.301585 4.102954 3.350000 1.580000
## 361 362 363 364 365 366
## 1.622500 2.155000 1.580000 2.150000 2.420000 2.060633
## 367 368 369 370 371 372
## 2.948384 1.930000 2.589580 5.254850 6.750000 6.545000
## 373 374 375 376 377 378
## 4.790000 2.850000 4.530000 3.280000 3.760000 6.610000
## 379 380 381 382 383 384
## 2.630000 4.830000 2.730000 2.350000 3.710000 5.150000
## 385 386 387 388 389 390
## 2.280000 3.760000 3.000000 1.540000 2.880485 2.554085
## 391 392 393 394 395 396
## 2.850000 1.890000 3.350000 2.530000 1.950000 1.830000
## 397 398 399 400 401 402
## 2.400000 4.490000 4.115000 3.810000 3.820000 3.620000
## 403 404 405 406 407 408
## 4.600000 4.140000 4.050000 4.370000 4.787160 4.060000
## 409 410 411 412 413 414
## 4.760000 4.300000 6.450000 3.830000 4.875000 4.670000
## 415 416 417 418 419 420
## 4.640000 2.420000 2.560000 1.750000 1.860000 1.992604
## 421 422 423 424 425 426
## 1.510000 2.520000 3.775000 4.470000 2.650000 2.077421
## 427 428 429 430 431 432
## 2.374929 2.991761 2.980000 7.110000 8.490000 8.950000
## 433 434 435 436 437 438
## 4.688845 6.860000 3.150000 3.290000 4.760000 2.910000
## 439 440 441 442 443 444
## 5.300000 5.600000 5.300000 6.700000 3.320000 5.450000
## 445 446 447 448 449 450
## 5.290000 5.000000 5.650000 4.230000 4.600000 7.310000
## 451 452 453 454 455 456
## 3.310000 3.990000 5.700000 6.800000 6.600000 5.740000
## 457 458 459 460 461 462
## 4.260000 3.671180 4.830000 4.900000 3.610000 6.020000
## 463 464
## 7.220980 3.379117Simulating traits
Now that we have established the best model(s), we might also want to see whether any species or clades in particular may be outliers, i.e. showing trait histories diverging from the expected. To do this, we can simulate trait histories under the given models to get a distribution of trait values per species.
# Simulate
simul_bm2<-mvSIM(two_regimes_tree,nsim=10000, model="BMM",param=fit_bm2)
# Find species that fall outside 90% confidence intervals
simul_bm2.ci<-t(apply(FUN=quantile, c(0.05, 0.95), MARGIN = 1, X = simul_bm2))
simul_bm2.outliers<-which(simul_bm2.ci[,1] > amphibia_logC | simul_bm2.ci[,2] < amphibia_logC)
simul_bm2.outliers## Necturus_lewisi Necturus_punctatus
## 93 94
## Desmognathus_wrighti Desmognathus_quadramaculatus
## 108 109
## Desmognathus_ochrophaeus Desmognathus_fuscus
## 111 113
## Petropedetes_cameronensis Platyplectrum_ornatum
## 265 337# plot expected distributions and true trait values
par(mfrow=c(4,2))
for(i in 1:length(simul_bm2.outliers)){
hist(simul_bm2[names(simul_bm2.outliers)[i],], main=names(simul_bm2.outliers)[i], cex=0.5, xlab="trait",las=1, border=NA, col="grey80")
abline(v=quantile(simul_bm2[names(simul_bm2.outliers)[i],], c(0.05, 0.95)), col="blue", lwd=1, lty=2) # quantile bands
abline(v=amphibia_logC[names(simul_bm2.outliers)[i]], col="red", lwd=2) # true value
}
Simulating ancestral states
I couldn’t find a way to get ancestral states for trait simulations (maybe there is a good reason for this), and so my work around has been to just perform ancestral state reconstructions on the simulated states. This seems somewhat circular to me, and so I want to put a big disclaimer here that maybe this isn’t the best way to go about doing things.
node.sim<-list()
## NOTE: here i am just using the first 100 simulations, but ideally you would run the loop for 1:ncol(simul_bm2)
for(i in 1:100){
node.sim[[i]]<-estim(two_regimes_tree, simul_bm2[,i], bm2.asr.fit, asr=TRUE)$estimates[,1]
}
node.sim<-as.data.frame(t(do.call(rbind, node.sim)))
# Now you have a table of node states for each simulated trait set (simulations as columns)
# the root note will inevitably be the same for all simulations as this is the nature of a BM model:
node.sim[1,]## V1 V2 V3 V4 V5 V6 V7
## 465 0.6642491 0.6642491 0.6642491 0.6642491 0.6642491 0.6642491 0.6642491
## V8 V9 V10 V11 V12 V13 V14
## 465 0.6642491 0.6642491 0.6642491 0.6642491 0.6642491 0.6642491 0.6642491
## V15 V16 V17 V18 V19 V20 V21
## 465 0.6642491 0.6642491 0.6642491 0.6642491 0.6642491 0.6642491 0.6642491
## V22 V23 V24 V25 V26 V27 V28
## 465 0.6642491 0.6642491 0.6642491 0.6642491 0.6642491 0.6642491 0.6642491
## V29 V30 V31 V32 V33 V34 V35
## 465 0.6642491 0.6642491 0.6642491 0.6642491 0.6642491 0.6642491 0.6642491
## V36 V37 V38 V39 V40 V41 V42
## 465 0.6642491 0.6642491 0.6642491 0.6642491 0.6642491 0.6642491 0.6642491
## V43 V44 V45 V46 V47 V48 V49
## 465 0.6642491 0.6642491 0.6642491 0.6642491 0.6642491 0.6642491 0.6642491
## V50 V51 V52 V53 V54 V55 V56
## 465 0.6642491 0.6642491 0.6642491 0.6642491 0.6642491 0.6642491 0.6642491
## V57 V58 V59 V60 V61 V62 V63
## 465 0.6642491 0.6642491 0.6642491 0.6642491 0.6642491 0.6642491 0.6642491
## V64 V65 V66 V67 V68 V69 V70
## 465 0.6642491 0.6642491 0.6642491 0.6642491 0.6642491 0.6642491 0.6642491
## V71 V72 V73 V74 V75 V76 V77
## 465 0.6642491 0.6642491 0.6642491 0.6642491 0.6642491 0.6642491 0.6642491
## V78 V79 V80 V81 V82 V83 V84
## 465 0.6642491 0.6642491 0.6642491 0.6642491 0.6642491 0.6642491 0.6642491
## V85 V86 V87 V88 V89 V90 V91
## 465 0.6642491 0.6642491 0.6642491 0.6642491 0.6642491 0.6642491 0.6642491
## V92 V93 V94 V95 V96 V97 V98
## 465 0.6642491 0.6642491 0.6642491 0.6642491 0.6642491 0.6642491 0.6642491
## V99 V100
## 465 0.6642491 0.6642491We could now look at the individual distributions per node if we are interested in a particular one, for example what would be the distribution of traits for the most recent common ancestor of salamanders:
Using the ggtree package, we could also plot these histograms directly on the nodes of the tree. To do this, I borrowed a bit of code from here.
## Registered S3 method overwritten by 'treeio':
## method from
## root.phylo ape## ggtree v1.17.4 For help: https://yulab-smu.github.io/treedata-book/
##
## If you use ggtree in published research, please cite the most appropriate paper(s):
##
## [36m-[39m Guangchuang Yu, Tommy Tsan-Yuk Lam, Huachen Zhu, Yi Guan. Two methods for mapping and visualizing associated data on phylogeny using ggtree. Molecular Biology and Evolution 2018, accepted. doi: 10.1093/molbev/msy194[36m-[39m Guangchuang Yu, David Smith, Huachen Zhu, Yi Guan, Tommy Tsan-Yuk Lam. ggtree: an R package for visualization and annotation of phylogenetic trees with their covariates and other associated data. Methods in Ecology and Evolution 2017, 8(1):28-36, doi:10.1111/2041-210X.12628##
## Attaching package: 'ggtree'## The following object is masked from 'package:ape':
##
## rotatelibrary(ggplot2)
# make ggtee object
gtree<-ggtree(amphibia.tree)
# define vector with node numbers as character strings for which we want to plot the histograms
nodes <- c("466","488","657","694")
# make the node annotation object
pd <- as.data.frame(node.sim[nodes,])
pdplots <- apply(pd, 1, function(y) {
ggplot(data.frame(y=y), aes(y, fill=..x..)) +
geom_histogram(binwidth = 0.1, alpha=0.75) +
xlim(range(unlist(pd))) +
scale_fill_gradient(low = "blue", high = "red") +
theme_inset() +
theme(axis.line.x = element_line(color="black", size = 0.5),
axis.text.x = element_text(size=8))
})
# plot tree and nodes
inset(gtree, pdplots, width=0.5, height=0.25, vjust=-1)## Loading required package: ggimage