jc69.rtmcmc: MCMC sampling of r and t under JC69

In dosreislab/simplephy: A simple package for Bayesian phylogenetic analysis

Description

MCMC sampling of r and t under JC69

Usage

jc69.rtmcmc(init.r, init.t, N, x, n, shape.r, rate.r, shape.t, rate.t, w.r,
  w.t)

Examples

data(humanorang)
n <- 100
r <- seq(from=0, to=6, len=n) * 1e-3
t <- seq(from=0, to=50, len=n)
rt <- expand.grid(t=t, r=r)
shape.r <- 2; rate.r <- .5 * 1e3 # diffuse prior with mean 4e-3 substitutions per My
shape.t <- 16.2; rate.t <- .72   # informative prior age of root with mean 22.5 Ma
# MLE of molecular distance, note d here is distance from tip to root
d.mle <- jc69.2smle(humanorang$x, humanorang$n) / 2

# Joint posterior on r and t:
z3 <- exp(jc69.rtlnP(rt$r, rt$t, x=humanorang$x, n=humanorang$n, shape.r, rate.r, shape.t, rate.t))
z3 <- matrix(z3, ncol=n)
image(t, r, z3, xlab="Root age (Ma)", ylab="Molecular rate (s/s/My)", main="Posterior", col=rev(heat.colors(12)), las=1)
contour(t, r, z3, add=TRUE, drawlabels=FALSE)
curve(d.mle/x/1e-3, col="blue", lwd=2, add=TRUE)

# MCMC sampling
rt.mcmc <- jc69.rtmcmc(init.r=1e-3, init.t=10, N=1e4, x=humanorang$x, n=humanorang$n,
  shape.r, rate.r, shape.t, rate.t, w.r=2.5e-3, w.t=20)
# Add trace to contour plot:
points(rt.mcmc$t, rt.mcmc$r, col="blue", pch='.')

dosreislab/simplephy documentation built on May 11, 2019, 7:27 p.m.