Simple OU

• Simple Ornstein-Uhlenbeck Models

rb                                         # or maybe ./rb depending on your system
getwd()
setwd("/mnt/c/Users/oschwery/Documents/Teaching/RevBayes_Workshop_UI2021/2-Traits/")

Data

dataset = 11
T = readTrees("data/trees.nex")[dataset]
data = readContinuousCharacterData("data/traits.nex")[dataset]

The OU Model

Tree

tree <- T

Rate Parameter

L <- 1e-3
U <- 1
sigma2 ~ dnLoguniform(L, U)

Adaptation Parameter

alpha ~ dnExponential(10)

Optimum

A <- -10
B <- 10
theta ~ dnUniform(A, B)

OU Model

We’re using dnPhyloOrnsteinUhlenbeckREML, the REML refers to the used algorithm to compute the likelihood (“reduced likelihood”, as described in Felsenstein J. 1985. Phylogenies and the comparative method. The American Naturalist.:1–15. 10.1086/284325). Alternatives would be dnPhyloOrnsteinUhlenbeckMVN (a phylogenetic covariance matrix approach) or dnPhyloOrnsteinUhlenbeckThreePoint (the fast three-point algorithm of Ho LS, Ané C. Asymptotic theory with hierarchical autocorrelation: Ornstein–Uhlenbeck tree models. The Annals of Statistics. 2013 Apr;41(2):957-81.).

X ~ dnPhyloOrnsteinUhlenbeckREML(tree, alpha, theta, sigma2^0.5, rootStates=theta)
X.clamp(data)

mymodel = model(theta)

Running an MCMC

moves    = VectorMoves()
monitors = VectorMonitors()

moves.append( mvScale(sigma2, weight=1.0) )
moves.append( mvScale(alpha, weight=1.0) )
moves.append( mvSlide(theta, weight=1.0) )

monitors.append( mnModel(filename="output/simple_OU.log", printgen=10) )
monitors.append( mnScreen(printgen=1000, sigma2, alpha, theta) )

mymcmc = mcmc(mymodel, monitors, moves, nruns=2, combine="mixed")
mymcmc.run(generations=50000)