Commit c743615e authored by Nicolas Salamin's avatar Nicolas Salamin
Browse files

Solved an issue with the passing of ancestral states to the likelihood...

Solved an issue with the passing of ancestral states to the likelihood function. It now solves the overestimation of the sigma2 under MCMC...
parent b68ad308
......@@ -37,30 +37,6 @@ fN2<-function(xa, xb, vpa, vpb, sigma2, anc) {
return(-(((xa - anc)^2 / vpa) + ((xb - anc)^2 / vpb)) / (2 * sigma2) - log(sqrt(2 * pi * sigma2 * (vpa + vpb))));
}
condLike<-function(sigma2=1.0, condvec, kids, anc) {
#following the paper by Felsenstein 1985 + Novembre and Slatkin 2008
#given two descendants {a, b} and an ancestor c:
#v'_a = v_a + S^2_a
#v'_b = v_b + S^2_b
#L_c = N(x_a - x_b, (v'_a + v'_b)*sigma^2)
#S^2_c = v'_a * v'_b / (v'_a) + v'_b)
#x_c = (v'_b * x_a) / (v'_a + v'_b) + (v'_a * x_b) / (v'_a + /v'_b): This is not used anymore as I'm keeping a vector of ancestral values as parameters
a <- kids[1];
b <- kids[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))
lnL <- fN2(condvec[a, 2], condvec[b, 2], vpa, vpb, sigma2, anc) + condvec[a, 1] + condvec[b, 1];# the log likelihood for the node
S <- (vpa * vpb) / (vpa + vpb);
return(c(lnL, anc, S));
}
#set the condL for the tips based on the data and the internal nodes
condvecInit<-function(tree, data, ntips, tnodes) {
condvec<-matrix(nrow=tnodes, ncol=4); #four columns: the first for the logLn, the second for the state, the third for extra BM variance and the fourth is the brlen subtending the node
......@@ -89,10 +65,40 @@ condvecInit<-function(tree, data, ntips, tnodes) {
return(condvec);
}
lnLike<-function(tree, data, anc, sigma2) {
condLike<-function(sigma2=1.0, condvec, kids, ancestor, useContrast=FALSE) {
#following the paper by Felsenstein 1985 + Novembre and Slatkin 2008
#given two descendants {a, b} and an ancestor c:
#v'_a = v_a + S^2_a
#v'_b = v_b + S^2_b
#L_c = N(x_a - x_b, (v'_a + v'_b)*sigma^2)
#S^2_c = v'_a * v'_b / (v'_a) + v'_b)
#x_c = (v'_b * x_a) / (v'_a + v'_b) + (v'_a * x_b) / (v'_a + /v'_b): This is not used anymore as I'm keeping a vector of ancestral values as parameters
a <- kids[1];
b <- kids[2];
vpa <- condvec[a, 4] + condvec[a, 3];
vpb <- condvec[b, 4] + condvec[b, 3];
#if ancestral states are not to be estimated, use the weighted mean (ie independent contrasts)
if(useContrast) {
anc <- ((condvec[a, 2] * vpb) / (vpa + vpb)) + ((condvec[b, 2] * vpa) / (vpa + vpb))
}
else {
anc <- ancestor;
}
lnL <- fN2(condvec[a, 2], condvec[b, 2], vpa, vpb, sigma2, anc) + condvec[a, 1] + condvec[b, 1];# the log likelihood for the node
S <- (vpa * vpb) / (vpa + vpb);
return(c(lnL, anc, S));
}
lnLike<-function(tree, data, anc, sigma2, useContrast=FALSE, getAnc=FALSE) {
ntips<-length(tree$tip.label);
nnodes<-tree$Nnode;
tnodes<-ntips+nnodes;
nodeOrder<-numeric(length(anc));
#conditional likelihood vectors
condvec<-condvecInit(tree=tree, data=data, ntips=ntips, tnodes=tnodes);
......@@ -102,29 +108,42 @@ lnLike<-function(tree, data, anc, sigma2) {
for(i in tnodes:(ntips+2)) {
row<-which(tree$edge[,2] == i);
kids<-which(tree$edge[,1] == i);
condvec[row, 1:3]<-condLike(sigma2=sigma2, condvec=condvec, kids, anc[i - ntips]); #the vector anc has the same order than the internal nodes of the tree. The root is thus at position 1, etc...
condvec[row, 1:3]<-condLike(sigma2=sigma2, condvec=condvec, kids, anc[i - ntips], useContrast); #the vector anc has the same order than the internal nodes of the tree. The root is thus at position 1, etc...
nodeOrder[i - ntips] <- row
}
#lnL<-tmp;
#finish with the root
kids<-which(tree$edge[,1]==(ntips+1));
condvec[tnodes, 1:3]<-condLike(sigma2=sigma2, condvec=condvec, kids, anc[1])
condvec[tnodes, 1:3]<-condLike(sigma2=sigma2, condvec=condvec, kids, anc[1], useContrast);
nodeOrder[1] <- tnodes;
return(condvec[tnodes, 1] + fN(anc[1], anc[1], condvec[tnodes, 3], 0, sigma2));
if(!getAnc) {
return(condvec[tnodes, 1] + fN(anc[1], anc[1], condvec[tnodes, 3], 0, sigma2));
}
else {
return(condvec[nodeOrder, 2]);
}
}
fitBM <- function(tree, x) {
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) {
lik <- function(par, tree, x, sigma) lnLike(tree, x, par[-1], sigma); #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, sigma=sigma, method = "L-BFGS-B", lower = rep(-Inf, tree$Nnode), upper = rep(Inf, tree$Nnode), control = list(fnscale = -1));
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 ###
#########################################
......@@ -175,21 +194,23 @@ draw_N2<-function(sigma2=1.0, condvec, kids, anc) {
#if estimating with weighted average:
#anc <- ((condvec[a, 2] * vpb) / (vpa + vpb)) + ((condvec[b, 2] * vpa) / (vpa + vpb))
#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))
#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);
#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));
#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));
#this is not possible because any move will necessarily decrease the lnL as for a given sigma, new_m is the best solution
new_anc <- rnorm(1, mean=new_m, sd=sqrt(sigma2*new_s));
#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));
#therefore, simply return it, which means not optimizing the ancestral states
#new_anc <- new_m
#to return the contrast only, uncomment the following:
new_anc <- new_m
return(c(new_anc, new_s));
}
......@@ -205,59 +226,63 @@ runGibbs <- function(tree, data, sigma2, anc) {
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[i - ntips]); #the vector anc has the same order than the internal nodes of the tree. The root is thus at position 1, etc...
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[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) {
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=rgamma(1, shape=shapeGamma, scale=scaleGamma), ancestor=runif(nodes, min=minUnif, max=maxUnif));#JO ????#
if(sum(fossil,na.rm=T)>1){
ml$anc[which(!is.na(fossil))]<-fossil[which(!is.na(fossil))]
}
current <- list(gen=1, lnL=0, priorS=0, priorA=rep(0, nodes), post=-100000, hastings=0, acceptance=0, sigma2=ml$sig2, ancestor=ml$anc);
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);
print(current$lnL);
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)}
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))
#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))]<- 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);
}
else {
current$priorA <- sapply(current$ancestor, dunif, min=minUnif, max=maxUnif, log=TRUE);
}
#if dunif pour NA et dnorm pour valeur fossil
......@@ -273,13 +298,17 @@ myMCMC <- function(tree, data, fossil,winsize=2, ngen=5000, samplefreq=100, minU
new <- current; #initiate the new generation with the values of the old one
if(runif(1) > 0.9) { #update ancestral state in 10% of the time
#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;
new$ancestor <- runGibbs(tree, data, current$sigma2, current$ancestor);
new$hastings <- 10000000; #always accept it
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
}
}
else { #update sigma2
prop <- multiplier(current$sigma2, winsize, runif(1));
......@@ -292,7 +321,7 @@ myMCMC <- function(tree, data, fossil,winsize=2, ngen=5000, samplefreq=100, minU
if(sum(fossil,na.rm=T)>1){
#input vecteur avec des NA valeurs fossiles.
new$priorA<-rep(NA,length(current$ancestor))
#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]
......@@ -301,20 +330,15 @@ myMCMC <- function(tree, data, fossil,winsize=2, ngen=5000, samplefreq=100, minU
#}
#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)
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);
}else{
}
else {
new$priorA <- sapply(new$ancestor, dunif, min=minUnif, max=maxUnif, log=TRUE);
}
#to modify if we have fossils.
#new$priorA <- sapply(new$ancestor, dnorm, mean=ancValue, sd=sqrt(new$sigma));
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
......@@ -334,9 +358,18 @@ myMCMC <- function(tree, data, fossil,winsize=2, ngen=5000, samplefreq=100, minU
}
require(phytools);
require(geiger);
tree<-pbtree(n=20, scale=1);
data<-fastBM(tree, sig2=0.2, a=0);
ml<-fitBM(tree, data)
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))
......@@ -351,8 +384,13 @@ plot(x[-1], out[-1])
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=5000, samplefreq=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)
\ No newline at end of file
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
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