# R/treeOperators.R In coalescentMCMC: MCMC Algorithms for the Coalescent

#### Documented in NeighborhoodRearrangementTipInterchange

```## treeOperators.R (2019-01-17)

##   Trees Operators for Running MCMC

## This file is part of the R-package `coalescentMCMC'.
## See the file ../COPYING for licensing issues.

.branching_times <- function(phy)
{
phy <- reorder(phy)
.Call(branchingTimesCall, phy\$edge,  phy\$edge.length)
}

NeighborhoodRearrangement <- function(phy, n, THETA, brtimes)
{
e1 <- phy\$edge[, 1L] # local copy
e2 <- phy\$edge[, 2L]

## sample a node excluding the root:
target <- sample.int(n - 2L, 1L) + n + 1L

### i1, i2, and i3 are indices of some edges
### target, anc, and sister are indices of some nodes

i1 <- which(e2 == target)
anc <- e1[i1] # the ancestor of 'target'
i2 <- which(e1 == target) # the 2 edges where 'target' is basal
i3 <- which(e1 == anc) # this includes i1, so:
i3 <- i3[i3 != i1]
sister <- e2[i3] # the sister-node of 'target'
sel <- sample.int(2L, 1L)
i2.move <- i2[sel]
i2.stay <- i2[-sel]
phy\$edge[i3, 2L] <- child2move <- e2[i2.move]
child2stay <- e2[i2.stay]
phy\$edge[i2.move, 2L] <- sister

bt.anc <- brtimes[anc - n]
bt.child2move <- if (child2move > n) brtimes[child2move - n] else 0
bt.sister <- if (sister > n) brtimes[sister - n] else 0
bt.child2stay <- if (child2stay > n) brtimes[child2stay - n] else 0

## adjust the branch length that was subtending 'sister':
phy\$edge.length[i3] <- bt.anc - bt.child2move

## random age for 'target' between the ones of 'sister' and 'anc':
## newage <- runif(1, max(bt.sister, bt.child2stay), bt.anc)

## random coalescent time:
## (gives a slightly better acceptance rate than the previous one)
agemin <- max(bt.sister, bt.child2stay)
pmax <- 1 - exp(-THETA * (bt.anc - agemin))
p <- runif(1, 0, pmax)
newage <- agemin - log(1 - p) / THETA

phy\$edge.length[i1] <- bt.anc - newage
phy\$edge.length[i2.move] <- newage - bt.sister
## adjust the branch length below the child that has not been moved:
phy\$edge.length[i2.stay] <- newage - bt.child2stay

attr(phy, "order") <- NULL
phy <- reorder(phy)
newNb <- integer(2*n - 1L)
newNb[n + 1L] <- n + 1L
sndcol <- phy\$edge[, 2L] > n
phy\$edge[sndcol, 2L] <- newNb[phy\$edge[sndcol, 2L]] <- 1:(n - 2L) + n + 1L
phy\$edge[, 1L] <- newNb[phy\$edge[, 1L]]
phy
}

TipInterchange <- function(phy, n)
{
e1 <- phy\$edge[, 1L]
e2 <- phy\$edge[, 2L]
repeat {
ij <- sample.int(n, size = 2L)
k <- match(ij, e2)
## check that the two tips in 'k' are not sisters
if (e1[k[1]] != e1[k[2]]) break
}
phy\$edge[k, 2L] <- ij[2:1]
phy
}

###EdgeLengthJittering <- function(phy)
###### all edge lengths are added a random value on U[-MIN, MAX]
###### (the ultrametric nature of the tree is kept)
###{
###    el <- phy\$edge.length
###    MIN <- min(el)
###    MAX <- max(el)
###    phy\$edge.length <- el + runif(1, -MIN, MAX)
###    phy
###}

## moves from Drummond et al. (2002, Genetics):

ScalingMove <- function(phy)
{
## values decreased compared to Drummond et al.'s
phy\$edge.length <- runif(1, 0.95, 1.05) * phy\$edge.length
phy
}

branchSwapping <- function(phy, n, brtimes)
{
e1 <- phy\$edge[, 1L]
e2 <- phy\$edge[, 2L]
Nedge <- length(e1)
ROOT <- n + 1L
## sample directly the edges:
x <- sample.int(Nedge, 2L)
wi <- x[1L]
wj <- x[2L]
i <- e2[wi]
j <- e2[wj]
ip <- e1[wi]
if (ip == ROOT) return(phy)
jp <- e1[wj]
if (ip == j || ip == jp || i == jp) return(phy)
if (jp > n && i > n && brtimes[jp - n] < brtimes[i - n]) return(phy)
wip <- which(e2 == ip)
ipp <- e1[wip]
tmp <- which(e1 == ip)
witild <- tmp[tmp != wi]
itild <- e2[witild]
phy\$edge[witild, 2L] <- j
phy\$edge[wj, 2L] <- ip
phy\$edge[wip, 2L] <- itild
##if (j != ROOT) { # j cannot be the root if phy is rooted
phy\$edge.length[wip] <- phy\$edge.length[wip] + phy\$edge.length[witild]
bti <- if (i > n) brtimes[i - n] else 0
btj <- if (j > n) brtimes[j - n] else 0
btjp <- brtimes[jp - n]
new.age.ip <- runif(1, max(bti, btj), btjp)
phy\$edge.length[wj] <- btjp - new.age.ip
phy\$edge.length[wi] <- new.age.ip - bti
phy\$edge.length[witild] <- new.age.ip - btj
##}
reorder(phy, "postorder")
}

subtreeExchange <- function(phy, n, brtimes)
{
phybak <- phy
e1 <- phy\$edge[, 1L]
e2 <- phy\$edge[, 2L]
Nedge <- length(e1)
ROOT <- n + 1L
## sample directly an edge:
wi <- sample.int(Nedge, 1L)
i <- e2[wi]
if (i == ROOT) return(phy)
ip <- e1[wi]
if (ip == ROOT) return(phy)
wip <- which(e2 == ip)
jp <- e1[wip]
tmp <- which(e1 == jp)
wj <- tmp[tmp != wip]
j <- e2[wj]
if (j > n) # no need to test the following condition if j is a tip
if (brtimes[ip - n] < brtimes[j - n]) return(phy)
phy\$edge[wi, 2L] <- j
phy\$edge[wj, 2L] <- i
## no need to adjust edge lengths if both i and j are tips
if (i > n || j > n) {
a <- phy\$edge.length[wi]
b <- phy\$edge.length[wip]
c <- phy\$edge.length[wj]
ab <- a + b
phy\$edge.length[wj] <- ab
phy\$edge.length[wi] <- c * a/ab
phy\$edge.length[wip] <- c * b/ab
tmp <- which(e1 == ip)
wib <- tmp[tmp != wi]
phy\$edge.length[wib] <- phy\$edge.length[wib] - phy\$edge.length[wip] + b
}
if (any(phy\$edge.length <= 0)) phybak else reorder(phy, "postorder")
}

NodeAgeMove <- function(phy, n, brtimes)
{
phybak <- phy
e1 <- phy\$edge[, 1L]
e2 <- phy\$edge[, 2L]
Nedge <- length(e1)
ROOT <- n + 1L
## sample a node:
i <- sample.int(n - 1L, 1L) + n
wi <- which(e1 == i)
j <- e2[wi[1]]
k <- e2[wi[2]]
ti <- brtimes[i - n]
tj <- if (j > n) brtimes[j - n] else 0
tk <- if (k > n) brtimes[k - n] else 0
mx <- max(tj, tk)
if (i == ROOT) {
delta <- runif(1, 0.8333333, 1.2)
newti <- mx + delta * (ti - delta * mx)
dif <- newti - ti
} else {
wip <- which(e2 == i)
ip <- e1[wip]
newti <- runif(1, mx, brtimes[ip - n])
dif <- newti - ti
phy\$edge.length[wip] <- phy\$edge.length[wip] - dif
}
phy\$edge.length[wi] <- phy\$edge.length[wi] + dif
if (any(phy\$edge.length <= 0)) phybak else phy
}

randomWalkThetaMu <-
function(theta, mu, range.theta = c(0, Inf), range.mu = c(0, Inf))
{
wtheta <- 1e-3
wmu <- 1e-3
newtheta <- theta + runif(1, -wtheta, wtheta)
newmu <- mu + runif(1, -wmu, wmu)
if (newtheta < range.theta[1] || newtheta > range.theta[2]) newtheta <- theta
if (newmu < range.mu[1] || newmu > range.mu[2]) newmu <- mu
c(theta, mu)
}
```

## Try the coalescentMCMC package in your browser

Any scripts or data that you put into this service are public.

coalescentMCMC documentation built on April 22, 2022, 5:05 p.m.