Tuesday, March 31, 2020

Some other updates to mcmcMk for Bayesian MCMC with the Mk model

I have been continuing to work on mcmcMk in phytools. The function does Bayesian MCMC for the Mk model. The other day I described a load of different major updates I had made to the function. Since then, I have just been working on improvements around the edges.

By default, the method plots the trace of the likelihood as the MCMC is running. This actually slows down computation quite a bit, so for practical purposes probably should be avoided. I have not updated the S3 plot method for the object, though, so that it can be used to generate a static version of the trace plots after the MCMC is done running. Let's see:

library(phytools)
packageVersion("phytools") ## updated from GitHub
## [1] '0.7.24'
## our data:
tree
## 
## Phylogenetic tree with 200 tips and 199 internal nodes.
## 
## Tip labels:
##  t11, t12, t8, t89, t90, t53, ...
## 
## Rooted; includes branch lengths.
x
##  t11  t12   t8  t89  t90  t53 t121 t150 t151  t61  t32  t85 t104 t105  t55 t113 
##    1    1    1    1    1    1    1    1    1    0    1    0    1    1    1    1 
## t114  t16  t17 t195 t196 t183 t184   t3  t63  t70  t71  t33  t44  t45  t36  t37 
##    1    1    1    1    1    1    1    1    0    0    1    1    1    1    0    0 
##  t93  t94 t179 t180 t193 t194  t24 t197 t198  t41  t25 t191 t192 t187 t188 t115 
##    1    1    1    1    0    0    0    1    1    1    1    1    1    1    1    1 
##  t96  t69  t30  t42  t87  t88  t50  t51 t177 t178 t166  t14  t97  t98  t74  t75 
##    1    0    1    1    1    1    0    1    1    1    1    0    1    1    1    1 
##  t31  t60  t73  t81  t82  t57   t7 t185 t186 t158 t159  t72  t35  t66  t67  t52 
##    1    1    1    1    1    1    1    1    1    1    1    1    0    1    1    0 
## t181 t182  t65  t83  t84  t64 t122 t123  t46  t47 t142 t143  t43   t2   t4   t5 
##    1    1    1    1    0    1    0    0    0    1    1    1    1    1    1    1 
## t137 t138 t133 t199 t200 t107  t76 t146 t147  t95 t111 t112  t62 t119 t120   t9 
##    1    0    1    1    1    1    0    1    1    0    0    1    0    1    1    1 
## t117 t118 t101  t77  t78  t91  t92  t40 t110 t169 t170  t58  t59 t160 t161  t56 
##    1    0    1    1    1    1    0    1    1    1    1    1    0    1    1    1 
## t154 t155  t26  t27 t152 t153  t49 t108 t109 t116 t124 t125  t29 t173 t174  t86 
##    1    1    1    0    1    1    1    0    0    1    1    0    1    1    1    1 
## t135 t136 t134  t68 t144 t145  t38 t131 t132  t79  t80  t34  t18 t140 t141  t99 
##    1    1    1    1    1    1    1    1    1    1    1    1    0    1    1    0 
## t127 t128  t48 t175 t176 t129 t130  t13 t156 t157 t139 t148 t149   t6  t20 t102 
##    1    1    1    1    1    1    1    0    1    1    0    0    1    1    1    1 
## t103 t100  t54 t167 t168  t39 t189 t190 t162 t163 t171 t172 t126  t28  t10   t1 
##    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1    1 
## t164 t165 t106  t23  t19  t21  t22  t15 
##    1    1    1    1    1    1    1    1 
## Levels: 0 1
mcmc<-mcmcMk(tree,x,model="ARD",ngen=20000,plot=FALSE)
## Running MCMC....
## gen  [1,0]   [0,1]   logLik  accept
## 1    0.1357  0.1357  -105.1208
## 101  0.4063  1.735   -86.5779    0.18
## 201  0.3771  1.2137  -88.1899    0.27
## 301  0.5402  2.1238  -86.4929    0.25
## 401  1.5722  5.5137  -89.2698    0.41
## 501  0.6964  2.0926  -88.1108    0.4
## 601  0.2261  1.4338  -89.1161    0.37
## 701  0.366   1.6579  -86.7456    0.28
## 801  0.4479  1.9002  -86.4486    0.26
## 901  0.4064  1.7107  -86.5999    0.25
## 1001 0.3189  1.5964  -87.1588    0.33
## 1101 0.438   1.5434  -87.1182    0.3
## 1201 0.8423  3.1835  -87.2122    0.47
## 1301 0.7107  5.6161  -91.887 0.6
## 1401 1.0065  3.9128  -87.6407    0.58
## 1501 0.3402  1.3901  -87.1779    0.46
## 1601 0.2763  1.2481  -87.8487    0.21
## 1701 0.4288  1.799   -86.5162    0.28
## 1801 0.5119  1.6194  -87.6261    0.36
## 1901 0.6065  3.3359  -87.4174    0.49
## 2001 0.4448  1.5674  -87.0883    0.61
## 2101 0.336   1.2723  -87.5561    0.32
## 2201 0.4812  1.7179  -86.901 0.52
## 2301 0.5145  1.6413  -87.5549    0.45
## 2401 0.6562  2.3686  -86.8828    0.42
## 2501 0.5675  2.0535  -86.7683    0.4
## 2601 0.2433  1.0488  -88.8104    0.44
## 2701 0.3899  1.5582  -86.8204    0.37
## 2801 0.422   2.2021  -86.7757    0.54
## 2901 0.2799  1.1325  -88.1681    0.34
## 3001 0.4402  1.7646  -86.5723    0.42
## 3101 0.2665  1.45    -87.9473    0.42
## 3201 0.1926  1.2161  -90.0367    0.48
## 3301 0.5739  2.5963  -86.4933    0.53
## 3401 0.5441  1.4406  -89.3762    0.64
## 3501 0.5201  2.267   -86.3877    0.36
## 3601 0.5078  2.3586  -86.4477    0.65
## 3701 0.3458  2.5261  -89.2152    0.58
## 3801 0.5208  1.4209  -89.056 0.66
## 3901 0.8962  4.5994  -87.9658    0.69
## 4001 1.0288  3.048   -89.0988    0.72
## 4101 0.8663  2.7482  -88.0857    0.69
## 4201 0.4627  2.9743  -88.2027    0.57
## 4301 0.4744  2.1205  -86.3981    0.5
## 4401 0.3305  0.9427  -89.9773    0.36
## 4501 0.3477  1.9463  -87.2451    0.41
## 4601 0.3964  2.1607  -86.9865    0.42
## 4701 0.489   1.5179  -87.8358    0.39
## 4801 0.6379  2.0963  -87.3412    0.43
## 4901 0.2687  1.6137  -88.1847    0.46
## 5001 0.4152  1.3467  -87.8133    0.42
## 5101 0.6136  2.6642  -86.5224    0.52
## 5201 0.3952  1.7459  -86.594 0.45
## 5301 0.6625  2.2434  -87.1957    0.53
## 5401 0.2854  1.6087  -87.7582    0.49
## 5501 0.3056  2.1394  -88.8528    0.37
## 5601 0.5826  3.2777  -87.4935    0.67
## 5701 0.7897  3.3535  -86.9847    0.64
## 5801 0.5144  1.9652  -86.5741    0.66
## 5901 0.5183  2.0377  -86.4934    0.44
## 6001 0.4614  2.3375  -86.6564    0.52
## 6101 0.4624  2.6789  -87.3576    0.58
## 6201 0.4415  1.4296  -87.6741    0.42
## 6301 0.4091  1.4479  -87.2499    0.37
## 6401 0.7051  2.7397  -86.7828    0.53
## 6501 0.8938  3.4224  -87.3372    0.71
## 6601 0.8745  4.0751  -87.4503    0.69
## 6701 0.4535  1.8252  -86.5245    0.59
## 6801 0.4888  3.1156  -88.2105    0.53
## 6901 0.6753  3.9967  -88.2523    0.52
## 7001 0.3154  1.5064  -87.1777    0.49
## 7101 0.3618  1.7275  -86.7772    0.43
## 7201 0.4844  1.6148  -87.3014    0.31
## 7301 0.6027  2.062   -87.0624    0.36
## 7401 0.2523  1.1684  -88.3204    0.39
## 7501 0.377   1.6359  -86.7145    0.54
## 7601 0.4581  2.0792  -86.4215    0.56
## 7701 0.1747  1.3706  -91.4918    0.36
## 7801 0.4835  1.6534  -87.1384    0.44
## 7901 0.4069  2.3206  -87.2084    0.54
## 8001 0.5024  1.8284  -86.7707    0.57
## 8101 0.5539  2.8492  -86.8667    0.6
## 8201 0.3739  1.4012  -87.2044    0.39
## 8301 0.4593  2.0997  -86.4262    0.42
## 8401 0.3534  1.7777  -86.8922    0.54
## 8501 0.2838  1.2867  -87.6988    0.39
## 8601 0.3209  1.3456  -87.3285    0.47
## 8701 0.5749  2.2371  -86.538 0.53
## 8801 0.4912  1.9618  -86.4758    0.51
## 8901 0.2892  1.8215  -88.1947    0.37
## 9001 0.271   1.3481  -87.8153    0.44
## 9101 0.599   2.4966  -86.4827    0.41
## 9201 0.5998  1.3223  -92.0136    0.6
## 9301 0.3759  1.8464  -86.733 0.49
## 9401 0.677   3.0678  -86.7523    0.57
## 9501 0.3331  1.252   -87.6341    0.65
## 9601 0.5929  1.9877  -87.1746    0.37
## 9701 0.5204  2.0126  -86.5346    0.56
## 9801 0.639   2.1994  -87.0719    0.54
## 9901 0.7195  2.5608  -87.071 0.39
## 10001    0.5244  1.679   -87.5115    0.63
## 10101    0.4258  2.8371  -88.4371    0.44
## 10201    0.4905  1.7704  -86.8278    0.67
## 10301    0.6522  2.5721  -86.6365    0.64
## 10401    0.3869  2.1999  -87.2085    0.41
## 10501    0.2145  1.1373  -89.1884    0.63
## 10601    0.32    1.5494  -87.1279    0.32
## 10701    0.4971  1.4564  -88.3475    0.33
## 10801    0.3288  1.2583  -87.5975    0.33
## 10901    0.5582  1.9416  -86.9543    0.26
## 11001    0.298   1.5212  -87.4112    0.41
## 11101    0.3925  2.5797  -88.24  0.57
## 11201    0.3464  1.1514  -88.2792    0.53
## 11301    0.4453  1.1717  -89.7669    0.34
## 11401    0.6008  1.6247  -89.0949    0.36
## 11501    0.6686  2.8379  -86.6382    0.66
## 11601    0.4374  1.4569  -87.4838    0.59
## 11701    0.4941  2.0442  -86.4129    0.48
## 11801    0.3075  1.0109  -88.9884    0.33
## 11901    0.4577  2.0432  -86.4129    0.45
## 12001    0.739   1.9534  -89.5397    0.6
## 12101    0.216   1.1412  -89.1422    0.56
## 12201    0.2645  1.6228  -88.3276    0.39
## 12301    0.2987  1.6587  -87.5664    0.42
## 12401    0.2493  1.0974  -88.5565    0.5
## 12501    0.4823  2.6338  -87.0234    0.51
## 12601    0.6478  2.775   -86.5914    0.53
## 12701    0.9594  3.3289  -87.8335    0.6
## 12801    0.8183  4.056   -87.5294    0.69
## 12901    0.6067  2.0576  -87.1153    0.61
## 13001    0.6112  2.3358  -86.6256    0.34
## 13101    0.4766  1.5212  -87.6432    0.59
## 13201    0.5055  3.3602  -88.695 0.52
## 13301    0.6627  2.5107  -86.7396    0.6
## 13401    0.3916  1.8177  -86.6044    0.55
## 13501    0.4899  2.1087  -86.3851    0.24
## 13601    0.4012  1.1885  -88.7344    0.47
## 13701    0.4984  1.8929  -86.6079    0.35
## 13801    0.6393  2.5536  -86.5921    0.51
## 13901    0.4755  3.0014  -88.0751    0.46
## 14001    0.3013  1.2367  -87.6985    0.44
## 14101    0.6486  2.528   -86.6486    0.45
## 14201    0.484   3.3338  -88.9931    0.61
## 14301    0.4215  2.1503  -86.6971    0.6
## 14401    0.7127  2.8521  -86.7636    0.57
## 14501    0.5132  3.08    -87.7589    0.61
## 14601    0.6629  1.8395  -88.8478    0.54
## 14701    0.3666  1.7077  -86.7368    0.46
## 14801    0.5318  2.0888  -86.4938    0.41
## 14901    0.4312  1.5002  -87.2266    0.45
## 15001    0.6946  2.5976  -86.8436    0.39
## 15101    0.4205  1.8298  -86.5022    0.57
## 15201    0.5157  1.1729  -91.5686    0.34
## 15301    0.5417  2.1208  -86.5032    0.44
## 15401    0.7411  2.0828  -88.8443    0.52
## 15501    0.3283  1.3551  -87.2844    0.38
## 15601    0.4673  1.4274  -88.0499    0.4
## 15701    0.5516  2.1714  -86.495 0.49
## 15801    0.6027  2.8873  -86.7026    0.51
## 15901    0.8289  4.0787  -87.5287    0.55
## 16001    0.597   2.8553  -86.6824    0.7
## 16101    0.3006  1.6897  -87.5855    0.54
## 16201    0.5171  3.3293  -88.4142    0.48
## 16301    1.028   3.2646  -88.5387    0.61
## 16401    0.3654  1.5814  -86.7997    0.49
## 16501    0.6002  2.0449  -87.0842    0.42
## 16601    0.1872  0.9333  -90.4503    0.36
## 16701    0.2449  1.3405  -88.366 0.44
## 16801    0.4965  1.3603  -89.0685    0.46
## 16901    0.5584  2.7247  -86.6578    0.61
## 17001    0.5717  1.486   -89.5842    0.53
## 17101    0.5259  2.1595  -86.4159    0.69
## 17201    0.5337  1.8196  -87.0765    0.57
## 17301    0.365   1.7547  -86.7632    0.42
## 17401    0.6923  2.6825  -86.7562    0.5
## 17501    0.4162  2.4751  -87.4745    0.41
## 17601    0.6844  4.3256  -88.9225    0.58
## 17701    0.9668  2.9587  -88.6418    0.59
## 17801    0.866   2.0879  -91.107 0.49
## 17901    0.4766  1.4898  -87.8149    0.4
## 18001    0.2491  1.053   -88.7275    0.39
## 18101    0.4709  1.3427  -88.7136    0.37
## 18201    0.484   1.8785  -86.5585    0.54
## 18301    0.5785  2.8227  -86.705 0.52
## 18401    0.7235  3.4094  -87.0066    0.76
## 18501    0.8042  2.5049  -88.0541    0.65
## 18601    0.7403  2.6089  -87.1612    0.62
## 18701    0.7836  2.1585  -89.1572    0.47
## 18801    0.3354  1.7731  -87.1186    0.49
## 18901    0.4295  2.2385  -86.7654    0.4
## 19001    0.3407  1.5118  -86.9727    0.49
## 19101    0.4458  1.7744  -86.5735    0.5
## 19201    0.7451  4.1086  -87.937 0.63
## 19301    1.0418  3.4559  -88.2982    0.63
## 19401    0.9847  4.045   -87.5478    0.54
## 19501    0.6078  2.5901  -86.4995    0.69
## 19601    0.4425  2.0591  -86.4654    0.38
## 19701    0.3169  1.5819  -87.1755    0.45
## 19801    0.4329  1.5163  -87.1767    0.32
## 19901    0.6255  1.8414  -88.2058    0.45
## Done.
plot(mcmc)

plot of chunk unnamed-chunk-1

OK, now let's compute a summary of the posterior sample. For fun, let's compare it to our ML estimates obtained using fitMk.

summary(mcmc)
## Assuming 20% burn-in as no burn-in was specified....
## 
## Mean value of Q from the post burn-in posterior sample:
##            0          1
## 0 -2.1532491  2.1532491
## 1  0.5070647 -0.5070647
## 
## Median value of Q from the post burn-in posterior sample:
##            0          1
## 0 -2.0228121  2.0228121
## 1  0.4725864 -0.4725864
## 
## 95% HPD interval computed either from the post burn-in
## samples or using 'coda':
##           lower    upper
## [1,0] 0.1864875 0.893992
## [0,1] 0.9809502 3.630833
fitMk(tree,x,model="ARD")
## Object of class "fitMk".
## 
## Fitted (or set) value of Q:
##           0         1
## 0 -2.174694  2.174694
## 1  0.501980 -0.501980
## 
## Fitted (or set) value of pi:
##   0   1 
## 0.5 0.5 
## 
## Log-likelihood: -86.381613 
## 
## Optimization method used was "nlminb"

Finally, let's see the other plotting methods that I've made for the mcmcMk object:

d<-density(mcmc)
## Assuming 20% burn-in as no burn-in was specified....
d
## 
##  Density [1,0] 
## 
##        x                 y            
##  Min.   :0.06512   Min.   :0.0000713  
##  1st Qu.:0.42849   1st Qu.:0.0330722  
##  Median :0.79185   Median :0.2792682  
##  Mean   :0.79185   Mean   :0.6873446  
##  3rd Qu.:1.15521   3rd Qu.:1.3244296  
##  Max.   :1.51857   Max.   :2.3310174  
## 
##  Density [1,0] 
## 
##        x                 y            
##  Min.   :0.06512   Min.   :0.0000713  
##  1st Qu.:0.42849   1st Qu.:0.0330722  
##  Median :0.79185   Median :0.2792682  
##  Mean   :0.79185   Mean   :0.6873446  
##  3rd Qu.:1.15521   3rd Qu.:1.3244296  
##  Max.   :1.51857   Max.   :2.3310174  
## 
## To plot enter plot('object_name') at the command line interface.
plot(d)

plot of chunk unnamed-chunk-3

Or we can show this as a matrix. (This is more useful for more than two states.)

plot(d,show.matrix=TRUE)

plot of chunk unnamed-chunk-4

That's pretty cool.

The tree & data were simulated as follows - so our estimates got pretty close to the truth, which is good too.

x<-sim.Mk(tree<-pbtree(n=200),Q<-matrix(c(-2,2,0.5,0.5),2,2,byrow=TRUE,
    dimnames=list(0:1,0:1)))