bayoudev

Development version tutorial

This tutorial guides the user through using some of the development features of bayou, including a new method of running mcmc objects and loading the data back into R, as well as increased capacity for including predictor variables. This is in constant development, so please bear with me as I increase functionality and options.

In this example, we will use bayou chelonia dataset to simulate data and estimate the simulated relationships. First, we need to install the development version of bayou.

require(devtools)
require(geiger)
install_github('uyedaj/bayou', ref = 'dev')

## Downloading github repo uyedaj/bayou@dev
## Installing bayou
## '/usr/lib/R/bin/R' --vanilla CMD INSTALL  \
##   '/tmp/RtmpZjl13c/devtools4fe078fddc37/uyedaj-bayou-dfe12a4'  \
##   --library='/home/josef/R/x86_64-pc-linux-gnu-library/3.1'  \
##   --install-tests 
## 
## Reloading installed bayou

require(bayou)

Now we need to simulate allometric relationships.

data(chelonia)
tree <- chelonia$phy
ntips <- length(chelonia$phy$tip.label)
pred <- cbind(chelonia$dat)
true.pars <- list(alpha=1, sig2=1, beta1=0.75, k=3, ntheta=4, theta=c(2,3,1,8))

Now we need to set where the shifts occur on the tree. We can do this manually, or interactively using bayou's "identifyBranches" function. This function will require you to click on 'k' branches to assign regimes.

shifts <- identifyBranches(tree, true.pars$k)

Now use the output to generate the simulated data. This will be more clean in the future, but for now this works.

true.pars$sb <- shifts$sb
true.pars$loc <- shifts$loc
true.pars$t2 <- 2:true.pars$ntheta
dat <- dataSim(true.pars, model = "OU", tree, SE=0)$dat

## [1] "Using mapped regimes from parameter list"

plot of chunk unnamed-chunk-5

cache <- bayou:::.prepare.ou.univariate(tree, dat,SE=0, pred=pred)
nn <- names(dat)
dat <- dat + true.pars$beta1*cache$pred
dat <- dat[,1]
names(dat) <-nn

Let's look at the simulated data. Shift locations:

tr <- pars2simmap(true.pars, tree)
plotSimmap(tr$tree, colors=tr$cols)

## no colors provided. using the following legend:
##        1        2        3        4 
##  "black"    "red" "green3"   "blue"

plot of chunk unnamed-chunk-6 Phenogram:

phenogram(tr$tree, dat, colors=tr$cols)

plot of chunk unnamed-chunk-7 Scatterplot against predictor variables:

plot(cache$pred, dat)

plot of chunk unnamed-chunk-8

Make the prior.

prior <- make.prior(tree, dists=list(dalpha="dhalfcauchy", dsig2="dhalfcauchy",
                                     dbeta1="dunif", dsb="dsb", dk="cdpois", dtheta="dunif"), 
                      param=list(dalpha=list(scale=1), dsig2=list(scale=1), 
                              dbeta1=list(min=-10, max=10), dk=list(lambda=15, kmax=200), dsb=list(bmax=1,prob=1), 
                                    dtheta=list(min=0, max=15)), model="ffancova")

plot of chunk unnamed-chunk-9 Make sure the prior function works:

prior(true.pars)

## [1] -41.31

First step is to create the MCMC object that will run our MCMC chain.

mymcmc <- bayou.makeMCMC(tree, dat, pred, SE=0, model="ffancova", prior=prior, startpar=list(alpha=0.1, sig2=3, beta1=1, k=1, ntheta=2, theta=c(4,4), sb=200, loc=0, t2=2), new.dir=TRUE, ticker.freq=1000)

plot of chunk unnamed-chunk-11

To execute the MCMC, specify the number of steps in parentheses. Additional calls using the same function will restart the chain where the previous run left off.

mymcmc$run(10000)