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)

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)))

Major updates to mcmcMk for Bayesian MCMC with the Mk model

In the midst of all the craziness that is March of 2020, for some reason over the past few days I have been trying to improve the phytools function mcmcMk, just to see if I can possibly make it useful for some phytools users. mcmcMk does MCMC for the Mk model - but note that this can also been done using the package diversitree, and (in fact) this is what I would generally recommend.

Nonetheless, I have managed to make mcmcMk kind of cool - so it could be worth checking out.

The first thing I had to do was address some issues with the default prior distribution. I had been using an exponential distribution - but this turns out to be a really terrible idea. I have now switched to a Gamma distribution, which works much, much better.

I also (2) let the parsimony score for the trait inform the default prior. In particular, I set the default value of \(\alpha\)/\(\beta\) (the mean from the \(\Gamma\) distribution) to be the parsimony score divded by the total tree length. The default value of \(\alpha\) is (for now) 0.1, which corresponds to a pretty flat distribution (though you'll see that I will use an even lower value of \(\alpha\) here).

The other thing (3) I was able to accomplish was to speed up the likelihood calculation - which in turn speeds up the MCMC. I had been doing this using fitMk in phytools for a fixed value of the transition matrix Q. I realized that it should run faster if I just modified fitMk to export the likelihood function, and then used this likelihood function in the MCMC calculations.

In addition to this (4), I also tested geiger::fitDiscrete which already exported it's likelihood function. It turns out that fitDiscrete's likelihood calculation is significantly faster than fitMk. As such, I have now added this method as the default for computing the likelihood - but users can still use fitMk if they want. (They'll need to, for instance, if they don't have geiger installed or if they want to input their data for x as a probability matrix.)

Finally (5), I added the option (also the default) to auto-tune the proposal distribution to target a user-specified acceptance rate for the MCMC. This value defaults to 0.5, but can be adjusted up or down by the user. (For a one-parameter model, a good target acceptance rate is 0.44; for more complex models it should be lower.)

Let's test it out with some simple data:

library(phytools)
packageVersion("phytools")

## [1] '0.7.23'

## here's our data
tree

## 
## Phylogenetic tree with 100 tips and 99 internal nodes.
## 
## Tip labels:
##  t7, t44, t80, t81, t21, t18, ...
## 
## Rooted; includes branch lengths.

x

##   t7  t44  t80  t81  t21  t18  t13  t84  t85  t11  t12   t2   t4  t78  t79  t27 
##    0    1    1    1    1    1    0    0    0    1    1    0    1    1    1    1 
##  t36  t70  t71  t42  t43  t14  t15  t73  t74  t69   t5   t1  t52  t53  t97  t98 
##    1    1    1    0    1    1    1    0    1    0    0    1    0    1    1    1 
##  t47  t48  t75  t76  t72  t28  t35  t37  t93  t94  t46  t99 t100  t30  t56  t90 
##    1    1    1    1    1    1    0    1    0    0    0    0    0    1    0    0 
##  t91  t50  t31  t23  t24  t61  t62  t67  t68  t20  t88  t89  t77  t54  t55  t39 
##    1    0    1    1    0    0    1    1    1    0    0    0    0    0    1    1 
##  t63  t64  t57  t58  t45   t3  t34  t86  t87  t59  t60  t51  t25  t26  t82  t83 
##    1    1    1    1    0    1    1    1    1    0    0    0    1    1    1    1 
##  t32  t33   t9  t10  t19  t65  t66  t29   t6   t8  t22  t49  t95  t96  t92  t16 
##    0    1    1    1    0    1    1    1    0    0    0    1    1    1    1    1 
##  t17  t38  t40  t41 
##    1    0    1    0 
## Levels: 0 1

mcmc<-mcmcMk(tree,x,ngen=2000,model="ER",prior=list(alpha=0.001))

## Running MCMC....
## gen  [1,0]   logLik  accept
## 1    0.2956  -67.298
## 101  0.757   -62.0417    0.59
## 201  0.8819  -61.9181    0.6
## 301  0.6133  -62.529 0.5
## 401  1.5944  -62.9026    0.55
## 501  1.2939  -62.3513    0.56
## 601  1.2689  -62.3083    0.54
## 701  0.766   -62.0253    0.5
## 801  0.8051  -61.9692    0.49
## 901  0.7627  -62.0312    0.5
## 1001 1.4013  -62.544 0.57
## 1101 1.091   -62.0426    0.61
## 1201 1.1005  -62.0543    0.5
## 1301 2.497   -64.3112    0.44
## 1401 1.4033  -62.5477    0.46
## 1501 0.5925  -62.6448    0.53
## 1601 0.9704  -61.9339    0.54
## 1701 0.5553  -62.8899    0.44
## 1801 0.6921  -62.2035    0.47
## 1901 1.4791  -62.6882    0.47
## Done.

Now, let's compare it to our fitted ML model:

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 -1.037971  1.037971
## 1  1.037971 -1.037971
## 
## Median value of Q from the post burn-in posterior sample:
##            0          1
## 0 -0.9082315  0.9082315
## 1  0.9082315 -0.9082315
## 
## 95% HPD interval computed either from the post burn-in
## samples or using 'coda':
##           lower    upper
## [1,0] 0.3608725 1.910995

fitMk(tree,x,model="ER")

## Object of class "fitMk".
## 
## Fitted (or set) value of Q:
##           0         1
## 0 -0.904984  0.904984
## 1  0.904984 -0.904984
## 
## Fitted (or set) value of pi:
##   0   1 
## 0.5 0.5 
## 
## Log-likelihood: -61.915559 
## 
## Optimization method used was "nlminb"

Or, perhaps more meaningful, let's compare our posterior distribution on the transition rate, q, to the likelihood surface. Since our prior distribution was chosen to be quite uninformative, these should be similar:

par(mfrow=c(2,1))
plot(density(mcmc))

## Assuming 20% burn-in as no burn-in was specified....

## to compute our likelihood surface:
Mk.surf<-function(tree,x,range.q){
    m<-length(unique(x))
    Q<-matrix(1,m,m,dimnames=list(levels(x),levels(x)))
    diag(Q)<--rowSums(Q)+1
    index.matrix<-Q
    diag(index.matrix)<-rep(NA,m)
    lik.f<-fitMk(tree,x,model="ER")$lik
    q<-seq(range.q[1],range.q[2],length.out=100)
    logL<-sapply(q,function(q,Q) lik.f(q*Q),Q=Q)
    plot(q,exp(logL),type="l",bty="l",
    main="likelihood surface for q",font.main=1)
}
Mk.surf(tree,x,range(density(mcmc)$Density$x))

## Assuming 20% burn-in as no burn-in was specified....

That's pretty cool so far.

Now, for fun, let's try an “ARD” model:

mcmc<-mcmcMk(tree,x,ngen=2000,model="ARD",prior=list(alpha=0.1),
    auto.tune=0.3)

## Running MCMC....
## gen  [1,0]   [0,1]   logLik  accept
## 1    0.2956  0.2956  -67.298
## 101  0.6481  1.3574  -59.8223    0.4
## 201  0.636   1.5911  -60.5524    0.4
## 301  0.6464  0.7386  -61.3897    0.11
## 401  1.3343  3.2096  -61.7304    0.38
## 501  1.4642  2.3867  -60.42  0.49
## 601  1.6126  2.0809  -61.1757    0.35
## 701  0.607   0.7238  -61.3178    0.18
## 801  0.6791  1.2245  -59.5928    0.24
## 901  1.6073  3.1604  -61.0521    0.25
## 1001 0.8206  1.5465  -59.5629    0.31
## 1101 1.3371  2.1958  -60.1967    0.24
## 1201 1.1337  1.7848  -59.8579    0.3
## 1301 0.8941  1.0917  -60.4888    0.26
## 1401 1.1789  2.1283  -60.0104    0.3
## 1501 0.6855  1.4508  -59.8168    0.26
## 1601 1.0183  1.0874  -61.4304    0.41
## 1701 0.6388  1.8892  -61.7907    0.23
## 1801 0.4876  1.2122  -60.8139    0.22
## 1901 0.589   1.1416  -59.8566    0.29
## Done.

(This probably needs to run longer - but you get the idea.)

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 -1.4867913  1.4867913
## 1  0.8876883 -0.8876883
## 
## Median value of Q from the post burn-in posterior sample:
##           0         1
## 0 -1.325659  1.325659
## 1  0.777271 -0.777271
## 
## 95% HPD interval computed either from the post burn-in
## samples or using 'coda':
##           lower    upper
## [1,0] 0.2959638 1.805645
## [0,1] 0.4251765 3.080158

fitMk(tree,x,model="ARD")

## Object of class "fitMk".
## 
## Fitted (or set) value of Q:
##           0         1
## 0 -1.394135  1.394135
## 1  0.804726 -0.804726
## 
## Fitted (or set) value of pi:
##   0   1 
## 0.5 0.5 
## 
## Log-likelihood: -59.489243 
## 
## Optimization method used was "nlminb"

plot(density(mcmc))

## Assuming 20% burn-in as no burn-in was specified....

That's kind of cool, right?

The data for this exercise were simulated as follows:

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