S3 methods for evol.rate.mcmc

I just added a bunch of new S3 method functionality to the old phytools function evol.rate.mcmc.

evol.rate.mcmc implements a mostly obsolete Bayesian MCMC method that we published way back in 2011 (Revell et al. 2011) - though I promise that we submitted it for publication in 2010 when such a thing still seemed quite novel.

The method still works & is as valid as the day it was conceived - so, in the belief that those two attributes will never go completely out of style, I decided to add some improve functionality to the implementation through a few nice S3 methods to facilitate user analysis of the posterior sample.

The method uses Bayesian MCMC to simultaneously sample the position of a single change in the Brownian evolutionary rate, along with the rates of evolution before & after this shift (σ₁² and σ₂²) from their joint posterior probability distribution. Of course, this kind of approach has been largely surpassed by more sophisticated rjMCMC that don't require that we specify the number of rate regimes a priori.

Here's how the function works with its new methods:

library(phytools)
packageVersion("phytools")

## [1] '0.6.78'

## our data
contMap(tree,x,ftype="off")
nodelabels(bg="white")

A ten-fold rate change occurred somewhere on this tree. With the data projected onto the edges can you guess where?

Let's see how our MCMC method does:

mcmc<-evol.rate.mcmc(tree,x,ngen=100000,
    control=list(print=1000))

## Control parameters (set by user or default):

## List of 11
##  $ sig1      : num 3.85
##  $ sig2      : num 3.85
##  $ a         : num 0.124
##  $ sd1       : num 0.77
##  $ sd2       : num 0.77
##  $ sda       : num 0.0248
##  $ kloc      : num 0.575
##  $ sdlnr     : num 1
##  $ rand.shift: num 0.05
##  $ print     : num 1000
##  $ sample    : num 100

## Starting MCMC run....

## gen  sig2(1) sig2(2) a       node    pos'n   logLik
## 0    3.8513  3.8513  0.1242  47  0.0136  -98.5199
## 1000 4.7472  1.1015  0.1901  52  0.7403  -91.6246
## 2000 6.5516  1.4274  0.5328  53  0.0617  -89.3125
## 3000 10.2626 2.0501  0.5239  53  0.7168  -90.0144
## 4000 7.1997  1.1218  0.3302  53  0.7341  -89.1354
## 5000 9.1095  2.5996  -0.1084 53  0.2925  -91.4485
## 6000 11.0535 0.7713  -0.6809 53  0.5378  -92.7336
## 7000 7.5563  1.7938  -1.053  53  0.5669  -89.9477
## 8000 8.9896  1.7616  -0.802  53  0.2126  -89.6181
## 9000 8.4694  1.7614  -0.5515 53  0.6622  -89.2766
## 10000    6.0513  1.0425  -0.7134 53  0.6643  -90.3772
## 11000    6.471   1.419   -0.6901 53  0.5803  -89.3697
## 12000    9.7293  1.8891  -0.7913 53  0.6397  -89.8248
## 13000    9.2759  1.4576  -0.3384 53  0.6261  -88.7554
## 14000    6.5533  2.8637  0.5887  53  0.1736  -92.6213
## 15000    7.2331  1.8722  0.0933  53  0.5766  -89.4507
## 16000    5.7497  1.2727  0.6667  53  0.0719  -89.9337
## 17000    6.8121  1.2405  0.1802  53  0.0513  -89.0958
## 18000    5.071   1.406   0.4616  53  0.3621  -90.1982
## 19000    9.6313  1.2834  0.4153  53  0.2423  -89.0269
## 20000    5.9635  1.5177  0.3161  53  0.226   -89.4581
## 21000    11.5237 1.129   -0.3834 53  0.0792  -89.6332
## 22000    10.5815 2.1737  -0.5135 53  0.23    -90.5686
## 23000    7.7775  1.5563  -0.0165 53  0.0858  -88.9187
## 24000    8.1525  1.0985  0.625   53  0.0985  -89.6473
## 25000    8.2706  1.6853  0.2743  53  0.4959  -89.0231
## 26000    8.0452  1.0981  0.4027  53  0.0352  -89.4198
## 27000    9.4267  1.6219  0.5374  53  0.1307  -89.2664
## 28000    5.6156  1.0878  0.7263  53  0.6033  -90.4557
## 29000    8.8053  1.4802  0.5739  53  0.6831  -88.8768
## 30000    9.7727  1.3725  0.1832  53  0.6082  -88.7225
## 31000    8.988   1.4811  0.4002  53  0.3966  -88.8803
## 32000    6.4915  2.1834  0.8295  53  0.6319  -90.7786
## 33000    10.254  1.1449  0.3363  53  0.6771  -89.1161
## 34000    7.0591  1.1797  -0.3883 53  0.5014  -89.1111
## 35000    6.1379  3.0045  -0.2445 53  0.3873  -93.0979
## 36000    8.3586  1.1711  -0.3547 53  0.3766  -88.9909
## 37000    7.0023  1.591   -0.4526 53  0.059   -89.2233
## 38000    6.4041  1.7832  -0.7375 53  0.4314  -89.9055
## 39000    1.2828  7.2478  -0.4113 52  0.2368  -91.8563
## 40000    1.8916  7.3795  -0.8626 53  0.3218  -89.9803
## 41000    2.265   7.7458  -1.6187 53  0.3123  -91.7966
## 42000    1.3243  4.7777  -0.9015 53  0.7689  -91.1631
## 43000    2.1923  9.6244  -0.789  52  0.6619  -94.6795
## 44000    3.5813  5.5468  -0.9939 53  0.7143  -95.4661
## 45000    1.371   6.1316  -0.7344 53  0.7284  -89.5926
## 46000    1.179   5.7775  -0.5685 53  0.4621  -89.8642
## 47000    1.5695  6.9089  -0.8996 53  0.7442  -89.5592
## 48000    1.6022  8.8566  -0.4645 53  0.4595  -89.0048
## 49000    1.3209  6.7781  -0.1912 53  0.6059  -88.8861
## 50000    1.6415  7.4061  0.1713  53  0.6908  -88.9324
## 51000    3.2032  11.4299 0.6192  56  0.0098  -94.7366
## 52000    3.9338  9.0551  0.3722  68  0.3228  -99.0375
## 53000    1.7368  7.026   -0.0803 53  0.3387  -89.2686
## 54000    1.4098  8.5177  0.2391  53  0.2528  -88.7882
## 55000    1.8195  6.7139  0.536   53  0.1726  -89.7451
## 56000    1.687   10.5018 1.2422  53  0.1475  -90.3855
## 57000    1.3256  10.1453 1.2338  53  0.7398  -89.9808
## 58000    9.0391  1.5572  0.6414  53  0.2098  -89.2143
## 59000    6.4992  1.2157  0.4783  53  0.5747  -89.2762
## 60000    9.8724  1.5831  0.4376  53  0.2202  -89.1691
## 61000    11.455  1.3943  1.023   53  0.7193  -89.7722
## 62000    9.5703  1.6223  0.9556  53  0.6436  -89.507
## 63000    10.3903 1.2642  1.0419  53  0.4744  -89.9045
## 64000    6.971   1.6686  0.3952  53  0.2802  -89.2752
## 65000    6.0203  1.6706  0.2386  52  0.5842  -91.9453
## 66000    6.2021  1.8017  0.691   53  0.1594  -90.0356
## 67000    10.4031 2.1176  0.324   53  0.3784  -90.2674
## 68000    6.2886  1.031   0.8026  53  0.2542  -90.5827
## 69000    3.7018  1.6842  0.4081  18  0.3528  -98.7487
## 70000    1.6217  4.6785  0.132   53  0.0197  -90.8257
## 71000    1.5233  6.8208  -0.1426 53  0.3395  -89.0002
## 72000    1.882   4.6169  -0.1161 53  0.0577  -91.4273
## 73000    1.3064  5.7539  0.3419  53  0.2928  -89.5444
## 74000    2.743   3.4316  -0.0436 81  0.5173  -97.9719
## 75000    2.4256  6.752   0.868   53  0.6803  -91.3509
## 76000    1.6759  6.0283  0.2312  53  0.5698  -89.5232
## 77000    1.2085  6.5887  0.1353  53  0.7581  -89.0178
## 78000    1.354   7.4218  0.3005  53  0.7033  -88.726
## 79000    1.067   8.0797  0.0765  53  0.3345  -89.2576
## 80000    2.1334  7.3914  0.0101  53  0.0781  -90.2052
## 81000    1.4198  7.3826  -0.0578 53  0.1697  -88.8282
## 82000    1.6638  12.4242 0.0161  53  0.7585  -89.3959
## 83000    4.8258  1.6303  -1.1191 53  0.27    -91.4818
## 84000    7.417   1.2518  0.0326  53  0.2785  -88.8508
## 85000    11.2467 1.0696  -0.7653 53  0.1236  -90.1938
## 86000    3.7294  4.4079  -0.8034 47  0.3295  -98.8226
## 87000    4.4821  1.4358  -0.2162 53  0.1034  -90.9554
## 88000    7.865   2.0663  -0.2737 53  0.2632  -89.9702
## 89000    6.7051  1.8765  -0.536  53  0.2218  -89.8579
## 90000    6.0155  1.7863  -0.5395 53  0.4477  -89.9448
## 91000    8.87    2.233   -0.7544 53  0.3387  -90.648
## 92000    7.7177  2.0816  -1.1313 53  0.372   -90.6617
## 93000    11.1098 2.114   -0.7646 53  0.1598  -90.6933
## 94000    6.9441  1.4662  -0.6712 52  0.3691  -92.0849
## 95000    7.9683  2.4851  -0.6059 53  0.1836  -91.2853
## 96000    7.1385  1.2423  -0.6553 53  0.4355  -89.2403
## 97000    8.3754  2.3841  0.0756  53  0.3514  -90.7713
## 98000    10.8048 2.3087  -0.237  53  0.5283  -90.7832
## 99000    6.4152  1.0962  0.1316  53  0.7499  -89.3923
## 1e+05    6.3585  1.0092  -0.1182 52  0.5094  -91.4805

## Done MCMC run.

Note that the rates are permitted to 'walk-past' each other, which explains why σ₂² can become σ₁² and vice versa. This will be taken care of when we summarize our posterior sample.

mcmc

## 
## Object of class "evol.rate.mcmc" containing the results from a
## the Bayesian MCMC analysis of a Brownian-motion rate-shift model.
## 
## MCMC was conducted for 1e+05 generations sampling every 100 
## generations.
## 
## The most commonly sampled rate shift(s) occurred on the edge(s)
## to node(s) 53.
## 
## Use the functions posterior.evolrate and minSplit for more detailed
## analysis of the posterior sample from this analysis.

obj<-summary(mcmc)

## 
## No burn-in specified. Excluding the first 20% by default.

## Loading required package: coda

obj

## 
## The shift with the minimum (squared) distance to all shifts in the
## posterior set is found 0.2178 along the branch leading to node 53.
## 
## Mean 'root-wise' rate from the posterior sample (sig1): 1.9069
##      95% HPD for sig1: [0.7315, 3.8458]
## Mean 'tip-wise' rate from the posterior sample (sig2): 7.3875
##      95% HPD for sig1: [3.1494, 12.528]
## 
## To plot the mean shift point run plot(...,type="min.split") on the object
## produced by this function.
## 
## To plot the probability of a shift by edge run plot(...,method="edge.prob")
## on the object produced by this function.

We can also plot our results in a couple of different ways:

plot(obj,method="edge.prob",piecol=c("blue","lightgrey"),
    ftype="off",mar=c(0.1,0.1,5.1,0.1))
title(main="posterior probability of a shift\nmapped along each edge of the tree",
    font.main=3)

plot(obj,ftype="off",
    mar=c(0.1,0.1,5.1,0.1))
title(main="position of the sampled rate-shift with the minimum\nsquared distance to all shifts in the sample",
    font.main=3)

In case you hadn't guessed, the shift occurred on the edge leading to node #53:

Here's how the data were simulated:

set.seed(99)
tree<-pbtree(n=50)
mtree<-paintSubTree(tree,53,"2",stem=0.5)
x<-sim.rates(mtree,setNames(c(1,10),1:2))