Commit 5784622f authored by Nicolas Salamin's avatar Nicolas Salamin
Browse files

New version of the MCMC with the old lnLike function within Liam Revell code for the MCMC

parent c743615e
#################################
## likelihood calculation in R ##
#################################
require(ape);
# We assume in this function that the tree is fully dichotomous and that 2 descendants are provided
# for each internal node.
# The parameter sigma is the standard deviation of the rate of evolution of the trait.
# The branch lengths for internal nodes should be extended a bit because the variance of the internal
# nodes is a bit larger than what is observed on the tree
# see eqn 23.46 on page 405 in Joe's book.
# However, condvec[,3] stores the standard deviation of S, which is sigma^2 * (v1v2) / (v1 + v2) if
# you consider two descendants with branch lengths v1 and v2. We want here to recover only
# (v1v2) / (v1 + v2), so we need to divide it by sigma^2 and take the square of it.
# The value in condvec[,3] is 0 for terminal taxa, where the branch lengths are "correct".
# The ancestral value is simply the weighted mean of the values, using the modified branch lengths
# as the variance (modified because for internal nodes the variance is a bit larger).
# see eqn 23.43 on page 405 in Joe's book
# Extra variance due to the composite construction of the ancestor based on the two descendants
# given here is sigma^2 (v1 v2) / (v1 + v2)
# see eqn 23.46 on page 405 in Joe's book
# Here, we store the standard deviation and we thus take the square-root of the variance.
# The likelihood is then simply the density of the normal distribution for the contrast of the two values
fN<-function(xa, xb, vpa, vpb, sigma2) {
#the normal density. Sigma2 is the variance. Same as log(dnorm(xa-xb, 0, sqrt(sigma * sigma * (vpa + vpb))))
return((-(xa - xb)^2 / (2 * sigma2 * (vpa + vpb)) - log(sqrt(2 * pi * sigma2 * (vpa + vpb)))));
......@@ -126,271 +97,121 @@ lnLike<-function(tree, data, anc, sigma2, useContrast=FALSE, getAnc=FALSE) {
}
}
fitBM <- function(tree, x) { #optimize both sigma and ancestral nodes
lik <- function(par, tree, x) lnLike(tree, x, par[-1], par[1]); #if you want to optimize, the function MUST return the likelihood, whatever the wish of the user...
fit <- optim(c(1, rep(0, tree$Nnode)), lik, tree = tree, x = x, method = "L-BFGS-B", lower = c(1e-12, rep(-Inf, tree$Nnode)), upper = c(Inf, rep(Inf, tree$Nnode)), control = list(fnscale = -1));
return(list(sig2 = fit$par[1], anc = fit$par[-1], logL = fit$value, convergence = fit$convergence));
}
fitBM_anc <- function(tree, x, sigma) { #optimize only the ancestral nodes
lik <- function(par, tree, x, sigma) lnLike(tree, x, par, sigma); #if you want to optimize, the function MUST return the likelihood, whatever the wish of the user...
fit <- optim(rep(0, tree$Nnode), lik, tree = tree, x = x, sigma=sigma, method = "L-BFGS-B", lower = rep(-Inf, tree$Nnode), upper = rep(Inf, tree$Nnode), control = list(fnscale = -1));
return(list(anc = fit$par, logL = fit$value, convergence = fit$convergence));
}
fitBM_sigma <- function(tree, x) { #optimize only sigma and ancestral nodes are the independent contrasts
lik <- function(par, tree, x) lnLike(tree, x, rep(0, tree$Nnode), par, useContrast=TRUE);
fit <- optim(c(1), lik, tree = tree, x = x, method = "L-BFGS-B", lower = c(1e-12), upper = c(Inf), control = list(fnscale = -1));
return(list(sig2 = fit$par, logL = fit$value, convergence = fit$convergence));
}
#########################################
### putting it into an MCMC framework ###
#########################################
slidingWindow <- function(i, d, u=0) {
# Slidign window proporal unconstrained at maximum
# For details of the method see http://people.sc.fsu.edu/~pbeerli/BSC-5936/10-12-05/Lecture_13.pdf
#
# Args:
# i: current value
# d: window size
#
# Returns:
# Proposal value (integer).
ii <- i + (runif(length(i), 0, 1) - 0.5) * d #MrBayes trick
return(list(v=ii, lnHastingsRatio=0))
}
multiplier <- function(i, d, u) {
# Multiplier proporal
# For details of the method see http://people.sc.fsu.edu/~pbeerli/BSC-5936/10-12-05/Lecture_13.pdf
#
# Args:
# i: current value
# d: window size
# u: a random value from a uniform distribution [0,1]
#
# Returns:
# Proposal value (integer).
# TODO: if d = 1 then error -> needs fixing
lambda=2*log(d)
m=exp(lambda*(u-0.5))
ii=i*m
return(list(v=ii, lnHastingsRatio=log(m)))
}
draw_N2<-function(sigma2=1.0, condvec, kids, anc) {
a <- kids[1];
mua <- condvec[a, 2];
b <- kids[2];
mub <- condvec[b, 2];
vpa <- condvec[a, 4] + condvec[a, 3];
vpb <- condvec[b, 4] + condvec[b, 3];
#if estimating with weighted average:
#anc <- ((condvec[a, 2] * vpb) / (vpa + vpb)) + ((condvec[b, 2] * vpa) / (vpa + vpb))
#the ancestral value that is obtained by combining the densities of the two descendant lineages
new_m <- ((mua/(vpa * sigma2)) + (mub/(vpb * sigma2))) / ((1/(vpa * sigma2)) + (1/(vpb * sigma2)))
#new_m <- ((condvec[a, 2] * vpb) / (vpa + vpb)) + ((condvec[b, 2] * vpa) / (vpa + vpb))
new_s <- ((vpa * sigma2) * (vpb * sigma2)) / ((vpa * sigma2) + (vpb * sigma2));
#new_s <- (vpa * vpb) / (vpa + vpb);
#We should then multiply this by the prior on the ancestral state itself, which is N(0, 0.2)
priorM <- 0;
priorS <- 0.2;
new_m2 <- ((new_m/(new_s)) + (priorM/(priorS))) / ((1/(new_s)) + (1/(priorS)));
new_s2 <- ((new_s) * (priorS)) / ((new_s) + (priorS));
#new ancestor is then drawn from a normal density with new_m2 and new_s2
new_anc <- rnorm(1, mean=new_m2, sd=sqrt(sigma2*new_s2));
#to return the contrast only, uncomment the following:
new_anc <- new_m
return(c(new_anc, new_s));
}
runGibbs <- function(tree, data, sigma2, anc) {
ntips<-length(tree$tip.label);
nnodes<-tree$Nnode;
tnodes<-ntips+nnodes;
newanc <- anc;
condvec<-condvecInit(tree=tree, data=data, ntips=ntips, tnodes=tnodes);
for(i in tnodes:(ntips+2)) {
row<-which(tree$edge[,2] == i);
kids<-which(tree$edge[,1] == i);
condvec[row, 2:3]<-draw_N2(sigma2=sigma2, condvec=condvec, kids, anc[row - ntips]); #the vector anc has the same order than the internal nodes of the tree. The root is thus at position 1, etc...
newanc[i - ntips] <- condvec[row, 2];
}
#lnL<-tmp;
#finish with the root
kids<-which(tree$edge[,1]==(ntips+1));
condvec[tnodes, 2:3]<-draw_N2(sigma2=sigma2, condvec=condvec, kids, anc[tnodes - ntips]);
newanc[1] <- condvec[tnodes, 2];
return(newanc);
}
myMCMC <- function(tree, data, fossil, winsize=2, ngen=5000, samplefreq=100, minUnif=-20, maxUnif=20, shapeGamma=1.1, scaleGamma=5, ancValue=0, MH=TRUE) {
o <- 6; #number of columns before the start of the parameters
n <- tree$Nnode + 1; #number of parameters, which is the number of nodes + sigma
tips <- length(tree$tip.label);
nodes <- tree$Nnode;
m <- o + n;
current <- list(gen=1, lnL=0, priorS=0, priorA=rep(0, nodes), post=-100000, hastings=0, acceptance=0, sigma2=1, ancestor=rep(0, nodes));
#estimate by ML the initial sigma:
current$sigma2<-fitBM_sigma(tree, data)$sig2
current$ancestor<-lnLike(tree, data, rep(0, tree$Nnode), current$sigma2, useContrast=TRUE, getAnc=TRUE)
current$sigma2<-rgamma(1, shape=shapeGamma, scale=scaleGamma);
current$ancestor<-runif(tree$Nnode, minUnif, maxUnif);
if(sum(fossil,na.rm=T)>1){
current$ancestor[which(!is.na(fossil))]<-fossil[which(!is.na(fossil))]
}
cat(paste(c("gen", "lnL", "priorS", paste("priorNode_", seq(tips + 1, tips + nodes, by = 1), sep=""), "post", "hastings", "acceptance", "sigma2", paste("ancNode_", seq(tips + 1, tips + nodes, by = 1), sep="")), collapse="\t", sep="\t"), file="logs.txt", append=FALSE);
cat("\n", file="logs.txt", append=TRUE);
cat(paste(unlist(current, use.names=FALSE), collapse="\t"), "\n", file="logs.txt", append=TRUE);
current$lnL <- lnLike(tree, data, current$ancestor, current$sigma2);
current$priorS <- dgamma(current$sigma2, shape=shapeGamma, scale=scaleGamma, log=TRUE);
#current$ancestor vector of values
#ordre dans edge (1er c'est la racine)
#current$priorA <- sapply(current$ancestor, dunif, min=minUnif, max=maxUnif, log=TRUE); #to modify if we have fossils. #JO#
#fossil.prior.func<-function(ancestor,fossil.mean){sapply(ancestor, dnorm, mean=fossil.mean, sd=sqrt(current$sigma2))}
fossil.prior.func<-function(ancestor, fossil.mean) { sapply(ancestor, dnorm, mean=fossil.mean, sd=0.0000001); }
if(sum(fossil,na.rm=T)>1){
#input vecteur avec des NA valeurs fossiles.
#current$priorA<-rep(NA,length(current$ancestor)) #utile?!?
ancestor<-current$ancestor[which(!is.na(fossil))]
fossil.mean<-fossil[which(!is.na(fossil))]
current$priorA[which(!is.na(fossil))]<- mapply(fossil.prior.func, ancestor, fossil.mean)
current$priorA[which(is.na(fossil))] <- sapply(current$ancestor[which(is.na(fossil))], dunif, min=minUnif, max=maxUnif, log=TRUE);
}
else {
current$priorA <- sapply(current$ancestor, dunif, min=minUnif, max=maxUnif, log=TRUE);
}
#if dunif pour NA et dnorm pour valeur fossil
current$post <- current$lnL + current$priorS + sum(current$priorA);
for(i in 2:ngen) {
if((i %% 100) == 0) {
cat(".");
}
if((i %% 10000) == 0) {
cat(" [", i, "]\n");
}
new <- current; #initiate the new generation with the values of the old one
if(runif(1) > 0.5) { #update ancestral state in 10% of the time
if(MH) {
pos <- sample(1:nodes, 5); #sample 5 of the ancestral nodes to updage
prop <- slidingWindow(current$ancestor[pos], winsize);
new$ancestor[pos] <- prop$v;
new$hastings <- prop$lnHastingsRatio;
}
else {
new$ancestor <- runGibbs(tree, data, current$sigma2, current$ancestor);
new$hastings <- 10000000; #always accept it
}
x<-mcmc.revell(tree, data, 10000)
mcmc.revell <- function (tree, x, ngen = 10000, control = list())
{
temp <- phyl.vcv(as.matrix(x), vcv(tree), 1)
sig2 <- temp$R[1, 1]
a <- temp$alpha
y <- rep(a, tree$Nnode - 1)
pr.mean <- c(1000, rep(0, tree$Nnode))
pr.var <- c(pr.mean[1]^2, rep(1000, tree$Nnode))
prop <- rep(0.01 * max(temp$C) * sig2, tree$Nnode + 1)
con = list(sig2 = sig2, a = a, y = y, pr.mean = pr.mean,
pr.var = pr.var, prop = prop, sample = 100)
con[(namc <- names(control))] <- control
con <- con[!sapply(con, is.null)]
message("Control parameters (set by user or default):")
str(con)
likelihood.revell <- function(C, invC, detC, x, sig2, a, y) {
z <- c(x, y) - a
logLik <- -z %*% invC %*% z/(2 * sig2) - nrow(C) * log(2 *
pi)/2 - nrow(C) * log(sig2)/2 - detC/2
return(logLik)
}
else { #update sigma2
prop <- multiplier(current$sigma2, winsize, runif(1));
new$sigma2 <- prop$v;
new$hastings <- prop$lnHastingsRatio;
likelihood <- function(C, invC, detC, x, sig2, a, y) {
lnLike(tree, x, c(a, y), sig2)
}
new$lnL <- lnLike(tree, data, new$ancestor, new$sigma2);
new$priorS <- dgamma(new$sigma2, shape=shapeGamma, scale=scaleGamma, log=TRUE);
if(sum(fossil,na.rm=T)>1){
#input vecteur avec des NA valeurs fossiles.
#new$priorA<-rep(NA,length(current$ancestor)) #à nouveau, est-ce utile?!?
#for(i in 1:length(which(!is.na(fossil)))){
# ancestor<-new$ancestor[which(!is.na(fossil))][i]
# fossil.mean<-fossil[which(!is.na(fossil))][i]
# new$priorA[which(!is.na(fossil))[i]]<-fossil.prior.func(ancestor,fossil.mean)
#}
#new$priorA[which(!is.na(fossil))]<- sapply(new$ancestor[which(!is.na(fossil))], dnorm, mean=fossil[which(!is.na(fossil))], sd=sqrt(current$sigma2));
ancestor<-new$ancestor[which(!is.na(fossil))]
fossil.mean<-fossil[which(!is.na(fossil))]
new$priorA[which(!is.na(fossil))]<-mapply(fossil.prior.func, ancestor,fossil.mean)
new$priorA[which(is.na(fossil))]<- sapply(new$ancestor[which(is.na(fossil))], dunif, min=minUnif, max=maxUnif, log=TRUE);
log.prior <- function(pr.mean, pr.var, sig2, a, y) pp <- dexp(sig2,
rate = 1/pr.mean[1], log = TRUE) + sum(dnorm(c(a, y),
mean = pr.mean[2:length(pr.mean)], sd = sqrt(pr.var[1 +
1:tree$Nnode]), log = TRUE))
C <- vcvPhylo(tree)
if (any(tree$edge.length <= (10 * .Machine$double.eps)))
stop("some branch lengths are 0 or nearly zero")
invC <- solve(C)
detC <- determinant(C, logarithm = TRUE)$modulus[1]
sig2 <- con$sig2
a <- con$a
y <- con$y
x <- x[tree$tip.label]
if (is.null(names(y)))
names(y) <- length(tree$tip) + 2:tree$Nnode
else y[as.character(length(tree$tip) + 2:tree$Nnode)]
L <- likelihood(C, invC, detC, x, sig2, a, y)
Pr <- log.prior(con$pr.mean, con$pr.var, sig2, a, y)
X <- matrix(NA, ngen/con$sample + 1, tree$Nnode + 3, dimnames = list(NULL,
c("gen", "sig2", length(tree$tip) + 1:tree$Nnode, "logLik")))
X[1, ] <- c(0, sig2, a, y, L)
message("Starting MCMC...")
for (i in 1:ngen) {
j <- (i - 1)%%(tree$Nnode + 1)
if (j == 0) {
sig2.prime <- sig2 + rnorm(n = 1, sd = sqrt(con$prop[j +
1]))
if (sig2.prime < 0)
sig2.prime <- -sig2.prime
L.prime <- likelihood(C, invC, detC, x, sig2.prime,
a, y)
Pr.prime <- log.prior(con$pr.mean, con$pr.var, sig2.prime,
a, y)
post.odds <- min(1, exp(Pr.prime + L.prime - Pr -
L))
if (post.odds > runif(n = 1)) {
if (i%%con$sample == 0)
X[i/con$sample + 1, ] <- c(i, sig2.prime, a,
y, L.prime)
sig2 <- sig2.prime
L <- L.prime
Pr <- Pr.prime
}
else if (i%%con$sample == 0)
X[i/con$sample + 1, ] <- c(i, sig2, a, y, L)
}
else if (j == 1) {
a.prime <- a + rnorm(n = 1, sd = sqrt(con$prop[j +
1]))
L.prime <- likelihood(C, invC, detC, x, sig2, a.prime,
y)
Pr.prime <- log.prior(con$pr.mean, con$pr.var, sig2,
a.prime, y)
post.odds <- min(1, exp(Pr.prime + L.prime - Pr -
L))
if (post.odds > runif(n = 1)) {
if (i%%con$sample == 0)
X[i/con$sample + 1, ] <- c(i, sig2, a.prime,
y, L.prime)
a <- a.prime
L <- L.prime
Pr <- Pr.prime
}
else if (i%%con$sample == 0)
X[i/con$sample + 1, ] <- c(i, sig2, a, y, L)
}
else {
k <- j - 1
y.prime <- y
y.prime[k] <- y[k] + rnorm(n = 1, sd = sqrt(con$prop[j +
1]))
L.prime <- likelihood(C, invC, detC, x, sig2, a,
y.prime)
Pr.prime <- log.prior(con$pr.mean, con$pr.var, sig2,
a, y.prime)
post.odds <- min(1, exp(Pr.prime + L.prime - Pr -
L))
if (post.odds > runif(n = 1)) {
if (i%%con$sample == 0)
X[i/con$sample + 1, ] <- c(i, sig2, a, y.prime,
L.prime)
y <- y.prime
L <- L.prime
Pr <- Pr.prime
}
else if (i%%con$sample == 0)
X[i/con$sample + 1, ] <- c(i, sig2, a, y, L)
}
}
else {
new$priorA <- sapply(new$ancestor, dunif, min=minUnif, max=maxUnif, log=TRUE);
}
new$post <- new$lnL + new$priorS + sum(new$priorA);
if((new$post - current$post + new$hastings) >= log(runif(1))) { #the new proposition is better than the old one: save it to the file and replace the current generation
current <- new;
current$acceptance <- current$acceptance + 1;
}
if((i %% samplefreq) == 0) {
current$gen = i;
cat(paste(unlist(current, use.names=FALSE), collapse="\t"), "\n", file="logs.txt", append=TRUE);
}
}
if((i %% 10000) != 0) {
cat("[ ", i, " ]\n");
}
}
require(phytools);
require(geiger);
tree<-pbtree(n=20, scale=1);
data<-fastBM(tree, sig2=0.2, a=0);
ml<-fitBM_sigma(tree, data)
fitBM(tree, data)
fitBM_anc(tree, data, sigma=ml$sig2)
anc<-lnLike(tree, data, rep(0, tree$Nnode), ml$sig2, useContrast=TRUE, getAnc=TRUE)
lnLike(tree, data, ml$anc, ml$sig2, useContrast=TRUE, getAnc=TRUE)
lnLike(tree, data, anc, ml$sig2, useContrast=FALSE, getAnc=TRUE)
fitContinuous(tree, data, model="BM") #get same lnL as with fitBM_sigma
x<-seq(0.00001, 15, by=0.5)
out<-numeric(length(x))
for(i in 1:length(x)) {
out[i]<-fitBM_anc(tree, data, sigma=x[i])$logL
}
plot(x[-1], out[-1])
#tree.edge<-tree$edge[order(tree$edge[,2]),][-(1:length(tree$tip.label)),]
#tree.edge2<-rbind(c(NA,21),tree.edge)
fossil<-rep(NA,length(tree$tip.label)-1)
fossil[sample(1:(length(tree$tip.label)-1),10)]<-sample(1:(length(tree$tip.label)-1),10)
t<-myMCMC(tree, data, fossil, winsize=2, ngen=20000, samplefreq=10, MH=FALSE)
plot(tree)
nodelabels()
#phenogram
log<-read.table("logs.txt",h=T)
hist(log[,29], breaks=500)
plot(log[,1],log[,29])
mcmcAnc<-log[,27:45]
apply(mcmcAnc, 2, mean) - anc
\ No newline at end of file
message("Done MCMC.")
return(X)
}
\ No newline at end of file
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment