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)
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)
Or we can show this as a matrix. (This is more useful for more than two states.)
plot(d,show.matrix=TRUE)
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)))
