R/tess.analysis.efbd.R
In hoehna/TESS: Diversification Rate Estimation and Fast Simulation of Reconstructed Phylogenetic Trees under Tree-Wide Time-Heterogeneous Birth-Death Processes Including Mass-Extinction Events
Defines functions
tess.analysis.efbd
Documented in
tess.analysis.efbd
#' tess.analysis.efbd
#'
#' @description Running an analysis on a given tree and estimating the diversification rates including rate-shifts and mass-extinction events.
#'
#' @param tree the phylogeny
#' @param initialSpeciationRate The initial value of the speciation rate when the MCMC is started.
#' @param initialExtinctionRate The initial value of the extinction rate when the MCMC is started.
#' @param initialFossilizationRate The initial value of the fossilization rate when the MCMC is started.
#' @param empiricalHyperPriors .
#' @param empiricalHyperPriorInflation .
#' @param empiricalHyperPriorForm .
#' @param speciationRatePriorMean mean parameter for the lognormal prior on lambda
#' @param speciationRatePriorStDev standard deviation for the lognormal prior on lambda
#' @param extinctionRatePriorMean mean parameter for the lognormal prior on mu
#' @param extinctionRatePriorStDev standard deviation for the lognormal prior on mu
#' @param fossilizationRatePriorMean mean parameter for the lognormal prior on phi
#' @param fossilizationRatePriorStDev standard deviation for the lognormal prior on phi
#' @param initialSpeciationRateChangeTime The initial value of the time points when speciation rate-shifts occur.
#' @param initialExtinctionRateChangeTime The initial value of the time points when extinction rate-shifts occur.
#' @param initialFossilizationRateChangeTime The initial value of the time points when fossilization rate-shifts occur.
#' @param estimateNumberSpeciationRateChanges Do we estimate the number of rate shifts?
#' @param estimateNumberExtinctionRateChanges Do we estimate the number of rate shifts?
#' @param estimateNumberFossilizationRateChanges Do we estimate the number of rate shifts?
#' @param numExpectedRateChanges Expected number of rate changes which follow a Poisson process.
#' @param samplingProbability probability of uniform sampling at present
#' @param samplingStrategy Which strategy was used to obtain the samples (taxa). Options are: uniform|diversified|age
#' @param missingSpecies The number of species missed which originated in a given time interval (empirical taxon sampling)
#' @param timesMissingSpecies The times intervals of the missing species (empirical taxon sampling)
#' @param tInitialMassExtinction The initial value of the vector of times of the mass-extinction events.
#' @param pInitialMassExtinction The initial value of the vector of survival probabilities of the mass-extinction events.
#' @param pMassExtinctionPriorShape1 The alpha (first shape) parameter of the Beta prior distribution for the survival probability of a mass-extinction event.
#' @param pMassExtinctionPriorShape2 The beta (second shape) parameter of the Beta prior distribution for the survival probability of a mass-extinction event.
#' @param estimateMassExtinctionTimes Do we estimate the times of mass-extinction events?
#' @param numExpectedMassExtinctions Expected number of mass-extinction events which follow a Poisson process.
#' @param estimateNumberMassExtinctions Do we estimate the number of mass-extinction events?
#' @param MRCA does the tree start at the mrca?
#' @param CONDITION do we condition the process on nothing|survival|taxa?
#' @param BURNIN number of starting iterations to be dropped
#' @param MAX_ITERATIONS The maximum number of iteration of the MCMC.
#' @param THINNING The frequency how often samples are recorded during the MCMC.
#' @param OPTIMIZATION_FREQUENCY The frequency how often the MCMC moves are optimized
#' @param CONVERGENCE_FREQUENCY The frequency how often we check for convergence?
#' @param MAX_TIME The maximum time the MCMC is allowed to run in seconds.
#' @param MIN_ESS The minimum number of effective samples (ESS) to assume convergence.
#' @param ADAPTIVE Do we use auto-tuning of the MCMC moves?
#' @param dir The subdirectory in which the output will be stored.
#' @param priorOnly Do we sample from the prior only?
#' @param verbose Do we want to print the progress of the MCMC?
#'
#' @return nothing. outputs are written to the file system
#' @export
#'
#' @examples
#' # we load the conifers as the test data set
#' data(conifers)
#'
#' # for the conifers we know what the total number of species is
#' total <- 630
#' # thus, we can compute what the sampling fraction is
#' rho <- (conifers$Nnode+1)/total
#'
#'
#' # next, we specify the prior mean and standard deviation
#' # for the speciation and extinction rate
#' mu_lambda = 0.15
#' std_lambda = 0.02
#' mu_mu = 0.09
#' std_mu = 0.02
#' mu_phi = 0.10
#' std_phi = 0.02
#'
#'
#'
#' # now we can run the entire analysis.
#' # note that a full analyses should be run much longer
#' tess.analysis.efbd( tree=conifers,
#' initialSpeciationRate = exp(mu_lambda),
#' initialExtinctionRate = exp(mu_mu),
#' initialFossilizationRate = exp(mu_phi),
#' empiricalHyperPriors = FALSE,
#' speciationRatePriorMean = mu_lambda,
#' speciationRatePriorStDev = std_lambda,
#' extinctionRatePriorMean = mu_mu,
#' extinctionRatePriorStDev = std_mu,
#' fossilizationRatePriorMean = mu_phi,
#' fossilizationRatePriorStDev = std_phi,
#' numExpectedRateChanges = 2,
#' samplingProbability = rho,
#' numExpectedMassExtinctions = 2,
#' BURNIN = 100,
#' MAX_ITERATIONS = 200,
#' THINNING = 10,
#' dir = "analysis_conifer")
tess.analysis.efbd <- function( tree,
initialSpeciationRate,
initialExtinctionRate,
initialFossilizationRate,
empiricalHyperPriors = TRUE,
empiricalHyperPriorInflation = 10.0,
empiricalHyperPriorForm = c("lognormal","normal","gamma"),
speciationRatePriorMean = 0.0,
speciationRatePriorStDev = 1.0,
extinctionRatePriorMean = 0.0,
extinctionRatePriorStDev = 1.0,
fossilizationRatePriorMean = 0.0,
fossilizationRatePriorStDev = 1.0,
initialSpeciationRateChangeTime = c(),
initialExtinctionRateChangeTime = c(),
initialFossilizationRateChangeTime = c(),
estimateNumberSpeciationRateChanges = TRUE,
estimateNumberExtinctionRateChanges = TRUE,
estimateNumberFossilizationRateChanges = TRUE,
numExpectedRateChanges = 2,
samplingProbability = 1,
samplingStrategy = "uniform",
missingSpecies = c(),
timesMissingSpecies = c(),
tInitialMassExtinction = c(),
pInitialMassExtinction = c(),
pMassExtinctionPriorShape1 = 5,
pMassExtinctionPriorShape2 = 95,
estimateMassExtinctionTimes = TRUE,
numExpectedMassExtinctions = 2,
estimateNumberMassExtinctions = TRUE,
MRCA = TRUE,
CONDITION = "survival",
BURNIN = 10000,
MAX_ITERATIONS = 200000,
THINNING = 100,
OPTIMIZATION_FREQUENCY = 500,
CONVERGENCE_FREQUENCY = 1000,
MAX_TIME = Inf, MIN_ESS = 500,
ADAPTIVE = TRUE,
dir = "" ,
priorOnly = FALSE,
verbose = TRUE) {
USE_EBDSTP = FALSE
orgDir <- getwd()
if ( dir != "" ) {
if (file.exists(dir)){
setwd(file.path(dir))
} else {
dir.create(file.path(dir))
setwd(file.path(dir))
}
}
if ( BURNIN > 0 & BURNIN < 1 ) {
burn <- MAX_ITERATIONS * BURNIN
} else {
burn <- BURNIN
}
MAX_TRIES <- 1000
if ( class(tree) == "phylo" ) {
use_tree_sample = FALSE
times <- tess.branching.times(tree)
} else if ( class(tree) == "multiPhylo" ) {
use_tree_sample = TRUE
times <- list()
for (i in 1:length(tree)) {
times[[i]] <- tess.branching.times(tree[[i]])
}
NUM_SAMPLED_TREES <- length(tree)
}
if ( CONDITION != "time" && CONDITION != "survival" && CONDITION != "taxa" ) {
stop("Wrong choice of argument for \"CONDITION\". Possible option are time|survival|taxa.")
}
if ( samplingStrategy != "uniform" && samplingStrategy != "diversified") {
stop("Wrong choice of argument for \"samplingStrategy\". Possible option are uniform|diversified.")
}
write.nexus(tree,file="input_tree.tre")
if ( empiricalHyperPriors == TRUE ) {
if( verbose ) {
cat("Estimating empirical hyper-parameters.\n\n")
}
empirical.prior.likelihood <- function(params) {
b <- params[1]
d <- params[2]
f <- params[3]
if ( use_tree_sample == TRUE ) {
ln_probs <- c()
max_prob <- -Inf
for (i in 1:NUM_SAMPLED_TREES) {
if ( USE_EBDSTP == FALSE ) {
ln_probs[i] <- tess.likelihood.efbd(nodes=times[[i]],
lambda=b,
mu=d,
phi=f,
samplingStrategy = samplingStrategy,
samplingProbability = samplingProbability,
MRCA = MRCA,
CONDITION = CONDITION,
log=TRUE)
} else {
ln_probs[i] <- tess.likelihood.ebdstp(nodes=times[[i]],
lambda=b,
mu=d,
phi=f,
r=0.0,
samplingStrategy = samplingStrategy,
samplingProbability = samplingProbability,
MRCA = MRCA,
CONDITION = CONDITION,
log=TRUE)
}
if ( max_prob < ln_probs[i] ) max_prob <- ln_probs[i]
}
sum <- 0.0
for (i in 1:NUM_SAMPLED_TREES) {
sum <- sum + exp( ln_probs[i] - max_prob ) / NUM_SAMPLED_TREES
}
lnl <- log( sum ) + max_prob
} else {
if ( USE_EBDSTP == FALSE ) {
lnl <- tess.likelihood.efbd(nodes=times,
lambda=b,
mu=d,
phi=f,
samplingStrategy = samplingStrategy,
samplingProbability = samplingProbability,
MRCA = MRCA,
CONDITION = CONDITION,
log=TRUE)
} else {
lnl <- tess.likelihood.ebdstp(nodes=times,
lambda=b,
mu=d,
phi=f,
r=0.0,
samplingStrategy = samplingStrategy,
samplingProbability = samplingProbability,
MRCA = MRCA,
CONDITION = CONDITION,
log=TRUE)
}
}
return(lnl)
}
if ( use_tree_sample == TRUE ) {
prior.speciation <- function(x) { dexp(x,rate=max(times[[1]]$age)/log(length(times[[1]]$age)),log=TRUE) }
} else {
prior.speciation <- function(x) { dexp(x,rate=max(times$age)/log(length(times$age)),log=TRUE) }
}
prior.speciation <- function(x) { dexp(x,rate=1.0,log=TRUE) }
prior.extinction <- function(x) { dexp(x,rate=1.0,log=TRUE) }
prior.fossilization <- function(x) { dexp(x,rate=1.0,log=TRUE) }
priors <- c("speciation"=prior.speciation,"extinction"=prior.extinction,"fossilization"=prior.fossilization)
tries <- 1
repeat {
starting.values <- runif(3,0,1)
starting.lnl <- empirical.prior.likelihood(starting.values)
if ( is.finite(starting.lnl) || tries >= MAX_TRIES ) {
break
}
tries <- tries + 1
}
if ( tries >= MAX_TRIES ) {
stop("Could not find valid starting values after ",tries," attempts.\n")
}
samples <- tess.mcmc( empirical.prior.likelihood,
priors,
starting.values,
logTransforms=c(TRUE,TRUE,TRUE),
delta=c(0.1,0.1,0.1),
iterations=20000,
burnin=2000,
verbose=verbose)
samples.lambda <- samples[,"speciation"]
m.lambda <- mean(samples.lambda)
v.lambda <- empiricalHyperPriorInflation * var(samples.lambda)
samples.mu <- samples[,"extinction"]
m.mu <- mean(samples.mu)
v.mu <- empiricalHyperPriorInflation * var(samples.mu)
samples.phi <- samples[,"fossilization"]
m.phi <- mean(samples.phi)
v.phi <- empiricalHyperPriorInflation * var(samples.phi)
lambda.hyperprior.fits <- c("lognormal"=Inf,"normal"=Inf,"gamma"=Inf)
mu.hyperprior.fits <- c("lognormal"=Inf,"normal"=Inf,"gamma"=Inf)
phi.hyperprior.fits <- c("lognormal"=Inf,"normal"=Inf,"gamma"=Inf)
if ( "lognormal" %in% empiricalHyperPriorForm ) {
lambda.hyperprior.fits["lognormal"] <- suppressWarnings(ks.test(samples.lambda,'plnorm',meanlog=log((m.lambda^2)/sqrt(v.lambda+m.lambda^2)),sdlog=sqrt( log(1+v.lambda/(m.lambda^2)))))$statistic
mu.hyperprior.fits["lognormal"] <- suppressWarnings(ks.test(samples.mu,'plnorm',meanlog=log((m.mu^2)/sqrt(v.mu+m.mu^2)),sdlog=sqrt( log(1+v.mu/(m.mu^2)))))$statistic
phi.hyperprior.fits["lognormal"] <- suppressWarnings(ks.test(samples.phi,'plnorm',meanlog=log((m.phi^2)/sqrt(v.phi+m.phi^2)),sdlog=sqrt( log(1+v.phi/(m.phi^2)))))$statistic
}
if ( "normal" %in% empiricalHyperPriorForm ) {
lambda.hyperprior.fits["normal"] <- suppressWarnings(ks.test(samples.lambda,'pnorm',mean=m.lambda,sd=sqrt(v.lambda)))$statistic
mu.hyperprior.fits["normal"] <- suppressWarnings(ks.test(samples.mu,'pnorm',mean=m.mu,sd=sqrt(v.mu)))$statistic
phi.hyperprior.fits["normal"] <- suppressWarnings(ks.test(samples.phi,'pnorm',mean=m.phi,sd=sqrt(v.phi)))$statistic
}
if ( "gamma" %in% empiricalHyperPriorForm ) {
lambda.hyperprior.fits["gamma"] <- suppressWarnings(ks.test(samples.lambda,'pgamma',shape=(m.lambda^2)/v.lambda,rate=m.lambda/v.lambda))$statistic
mu.hyperprior.fits["gamma"] <- suppressWarnings(ks.test(samples.mu,'pgamma',shape=(m.mu^2)/v.mu,rate=m.mu/v.mu))$statistic
phi.hyperprior.fits["gamma"] <- suppressWarnings(ks.test(samples.phi,'pgamma',shape=(m.phi^2)/v.phi,rate=m.phi/v.phi))$statistic
}
lambda.hyper.prior.form <- names(which.min(lambda.hyperprior.fits[empiricalHyperPriorForm]))
mu.hyper.prior.form <- names(which.min(mu.hyperprior.fits[empiricalHyperPriorForm]))
phi.hyper.prior.form <- names(which.min(phi.hyperprior.fits[empiricalHyperPriorForm]))
if ( lambda.hyper.prior.form == "lognormal" ) {
speciationRatePriorMean <- mean( log(samples.lambda) )
speciationRatePriorStDev <- sqrt( empiricalHyperPriorInflation * mean( (log(samples.lambda)-speciationRatePriorMean)^2 ) )
} else if ( lambda.hyper.prior.form == "normal" ) {
speciationRatePriorMean <- m.lambda
speciationRatePriorStDev <- sqrt(v.lambda)
} else if ( lambda.hyper.prior.form == "gamma" ) {
speciationRatePriorMean <- (m.lambda^2)/v.lambda
speciationRatePriorStDev <- m.lambda/v.lambda
}
if ( mu.hyper.prior.form == "lognormal" ) {
extinctionRatePriorMean <- mean( log(samples.mu) )
extinctionRatePriorStDev <- sqrt( empiricalHyperPriorInflation * mean( (log(samples.mu)-extinctionRatePriorMean)^2 ) )
} else if ( mu.hyper.prior.form == "normal" ) {
extinctionRatePriorMean <- m.mu
extinctionRatePriorStDev <- sqrt(v.mu)
} else if ( mu.hyper.prior.form == "gamma" ) {
extinctionRatePriorMean <- (m.mu^2)/v.mu
extinctionRatePriorStDev <- m.mu/v.mu
}
if ( phi.hyper.prior.form == "lognormal" ) {
fossilizationRatePriorMean <- mean( log(samples.phi) )
fossilizationRatePriorStDev <- sqrt( empiricalHyperPriorInflation * mean( (log(samples.phi)-fossilizationRatePriorMean)^2 ) )
} else if ( phi.hyper.prior.form == "normal" ) {
fossilizationRatePriorMean <- m.phi
fossilizationRatePriorStDev <- sqrt(v.phi)
} else if ( phi.hyper.prior.form == "gamma" ) {
fossilizationRatePriorMean <- (m.phi^2)/v.phi
fossilizationRatePriorStDev <- m.phi/v.phi
}
initialSpeciationRate <- m.lambda
initialExtinctionRate <- m.mu
initialFossilizationRate <- m.phi
empirical.hyperprior.samples <- data.frame(1:nrow(samples),samples[,"likelihood"],samples[,"prior"],samples.lambda,samples.mu,samples.phi)
write.table(empirical.hyperprior.samples[-1,],file='samplesHyperprior.txt',row.names=FALSE,col.names=c('Iteration','Likelihood','Prior','speciation','extinction','fossilization'),sep='\t')
summary.file <- "empiricalHyperpriorSummary.txt"
cat("Summary of the empirical hyperprior analysis",file=summary.file)
cat("\n============================================\n\n",file=summary.file,append=TRUE)
cat("Specation rate hyperprior\n",file=summary.file,append=TRUE)
cat("-------------------------\n",file=summary.file,append=TRUE)
cat("Form:\t",lambda.hyper.prior.form,'\n',file=summary.file,append=TRUE)
if ( lambda.hyper.prior.form == "lognormal" ) {
cat("meanlog:\t",speciationRatePriorMean,"\n",file=summary.file,append=TRUE)
cat("sdlog:\t\t",speciationRatePriorStDev,"\n",file=summary.file,append=TRUE)
} else if ( lambda.hyper.prior.form == "normal" ) {
cat("mean:\t",speciationRatePriorMean,"\n",file=summary.file,append=TRUE)
cat("sd:\t\t",speciationRatePriorStDev,"\n",file=summary.file,append=TRUE)
} else if ( lambda.hyper.prior.form == "gamma" ) {
cat("shape:\t",speciationRatePriorMean,"\n",file=summary.file,append=TRUE)
cat("rate:\t\t",speciationRatePriorStDev,"\n",file=summary.file,append=TRUE)
}
cat("\nExtinction rate hyperprior\n",file=summary.file,append=TRUE)
cat("-------------------------\n",file=summary.file,append=TRUE)
cat("Form:\t",mu.hyper.prior.form,'\n',file=summary.file,append=TRUE)
if ( mu.hyper.prior.form == "lognormal" ) {
cat("meanlog:\t",extinctionRatePriorMean,"\n",file=summary.file,append=TRUE)
cat("sdlog:\t\t",extinctionRatePriorStDev,"\n",file=summary.file,append=TRUE)
} else if ( mu.hyper.prior.form == "normal" ) {
cat("mean:\t",extinctionRatePriorMean,"\n",file=summary.file,append=TRUE)
cat("sd:\t\t",extinctionRatePriorStDev,"\n",file=summary.file,append=TRUE)
} else if ( mu.hyper.prior.form == "gamma" ) {
cat("shape:\t",extinctionRatePriorMean,"\n",file=summary.file,append=TRUE)
cat("rate:\t\t",extinctionRatePriorStDev,"\n",file=summary.file,append=TRUE)
}
cat("\nFossilization rate hyperprior\n",file=summary.file,append=TRUE)
cat("-------------------------\n",file=summary.file,append=TRUE)
cat("Form:\t",phi.hyper.prior.form,'\n',file=summary.file,append=TRUE)
if ( phi.hyper.prior.form == "lognormal" ) {
cat("meanlog:\t",fossilizationRatePriorMean,"\n",file=summary.file,append=TRUE)
cat("sdlog:\t\t",fossilizationRatePriorStDev,"\n",file=summary.file,append=TRUE)
} else if ( phi.hyper.prior.form == "normal" ) {
cat("mean:\t",fossilizationRatePriorMean,"\n",file=summary.file,append=TRUE)
cat("sd:\t\t",fossilizationRatePriorStDev,"\n",file=summary.file,append=TRUE)
} else if ( phi.hyper.prior.form == "gamma" ) {
cat("shape:\t",fossilizationRatePriorMean,"\n",file=summary.file,append=TRUE)
cat("rate:\t\t",fossilizationRatePriorStDev,"\n",file=summary.file,append=TRUE)
}
} else { # no empirical hyper priors
lambda.hyper.prior.form <- mu.hyper.prior.form <- phi.hyper.prior.form <- empiricalHyperPriorForm[1]
}
likelihood <- function(lambda,lambdaTimes,mu,muTimes,phi,phiTimes,tMassExtinction,pSurvival) {
if ( priorOnly == FALSE ) {
if( any( lambda < 0 ) | any( mu < 0 ) ) {
lnl <- -Inf
} else {
if ( use_tree_sample == TRUE ) {
ln_probs <- c()
max_prob <- -Inf
for (i in 1:NUM_SAMPLED_TREES) {
if ( USE_EBDSTP == FALSE ) {
ln_probs[i] <- tess.likelihood.efbd(node = times[[i]],
lambda = lambda,
mu = mu,
phi = phi,
rateChangeTimesLambda = lambdaTimes,
rateChangeTimesMu = muTimes,
rateChangeTimesPhi = phiTimes,
massExtinctionTimes= tMassExtinction,
massExtinctionSurvivalProbabilities = pSurvival,
samplingStrategy = samplingStrategy,
samplingProbability = samplingProbability,
MRCA = MRCA,
CONDITION = CONDITION,
log = TRUE)
} else {
ln_probs[i] <- tess.likelihood.ebdstp(node = times[[i]],
lambda = lambda,
mu = mu,
phi = phi,
r = 0,
rateChangeTimesLambda = lambdaTimes,
rateChangeTimesMu = muTimes,
rateChangeTimesPhi = phiTimes,
massDeathTimes= tMassExtinction,
massDeathProbabilities = 1.0 - pSurvival,
samplingStrategy = samplingStrategy,
samplingProbability = samplingProbability,
MRCA = MRCA,
CONDITION = CONDITION,
log = TRUE)
}
if ( max_prob < ln_probs[i] ) max_prob <- ln_probs[i]
}
sum <- 0.0
for (i in 1:NUM_SAMPLED_TREES) {
sum <- sum + exp( ln_probs[i] - max_prob ) / NUM_SAMPLED_TREES
}
lnl <- log( sum ) + max_prob
} else {
if ( USE_EBDSTP == FALSE ) {
lnl <- tess.likelihood.efbd(node = times,
lambda = lambda,
mu = mu,
phi = phi,
rateChangeTimesLambda = lambdaTimes,
rateChangeTimesMu = muTimes,
rateChangeTimesPhi = phiTimes,
massExtinctionTimes= tMassExtinction,
massExtinctionSurvivalProbabilities = pSurvival,
samplingStrategy = samplingStrategy,
samplingProbability = samplingProbability,
MRCA = MRCA,
CONDITION = CONDITION,
log = TRUE)
} else {
lnl <- tess.likelihood.ebdstp(node = times,
lambda = lambda,
mu = mu,
phi = phi,
r = 0,
rateChangeTimesLambda = lambdaTimes,
rateChangeTimesMu = muTimes,
rateChangeTimesPhi = phiTimes,
massDeathTimes= tMassExtinction,
massDeathProbabilities = 1.0 - pSurvival,
samplingStrategy = samplingStrategy,
samplingProbability = samplingProbability,
MRCA = MRCA,
CONDITION = CONDITION,
log = TRUE)
}
}
}
} else {
lnl <- 0
}
return (lnl)
}
prior <- function(timesLambda,valuesLambda,timesMu,valuesMu,timesPhi,valuesPhi,timesMassExtinction,valuesMassExtinction) {
if( any( valuesLambda < 0 ) | any( valuesMu < 0 ) | any( valuesPhi < 0 ) ) {
return (-Inf)
}
#####################
## Speciation Rate ##
#####################
# prior on number of changes for lambda
k <- length(timesLambda)
lnp <- dpois(k,numExpectedRateChanges,log=TRUE)
if ( k > 0 ) {
# prior on change times
lnp <- lnp + sum( dunif(timesLambda,0,AGE,log=TRUE) )
}
# prior on change values
if ( lambda.hyper.prior.form == "lognormal" ) {
lnp <- lnp + sum( dlnorm(valuesLambda,meanlog=speciationRatePriorMean,speciationRatePriorStDev,log=TRUE) )
} else if ( lambda.hyper.prior.form == "normal" ) {
lnp <- lnp + sum( dnorm(valuesLambda,speciationRatePriorMean,speciationRatePriorStDev,log=TRUE) )
} else if ( lambda.hyper.prior.form == "gamma" ) {
lnp <- lnp + sum( dgamma(valuesLambda,speciationRatePriorMean,speciationRatePriorStDev,log=TRUE) )
}
#####################
## Extinction Rate ##
#####################
# prior on number of changes for mu
k <- length(timesMu)
lnp <- lnp + dpois(k,numExpectedRateChanges,log=TRUE)
if ( k > 0 ) {
# prior on change times
lnp <- lnp + sum( dunif(timesMu,0,AGE,log=TRUE) )
}
# prior on change values
if ( mu.hyper.prior.form == "lognormal" ) {
lnp <- lnp + sum( dlnorm(valuesMu,meanlog=extinctionRatePriorMean,extinctionRatePriorStDev,log=TRUE) )
} else if ( mu.hyper.prior.form == "normal" ) {
lnp <- lnp + sum( dnorm(valuesMu,extinctionRatePriorMean,extinctionRatePriorStDev,log=TRUE) )
} else if ( mu.hyper.prior.form == "gamma" ) {
lnp <- lnp + sum( dgamma(valuesMu,extinctionRatePriorMean,extinctionRatePriorStDev,log=TRUE) )
}
########################
## Fossilization Rate ##
########################
# prior on number of changes for phi
k <- length(timesPhi)
lnp <- lnp + dpois(k,numExpectedRateChanges,log=TRUE)
if ( k > 0 ) {
# prior on change times
lnp <- lnp + sum( dunif(timesPhi,0,AGE,log=TRUE) )
}
# prior on change values
if ( phi.hyper.prior.form == "lognormal" ) {
lnp <- lnp + sum( dlnorm(valuesPhi,fossilizationRatePriorMean,fossilizationRatePriorStDev,log=TRUE) )
} else if ( phi.hyper.prior.form == "normal" ) {
lnp <- lnp + sum( dnorm(valuesPhi,fossilizationRatePriorMean,fossilizationRatePriorStDev,log=TRUE) )
} else if ( phi.hyper.prior.form == "gamma" ) {
lnp <- lnp + sum( dgamma(valuesPhi,fossilizationRatePriorMean,fossilizationRatePriorStDev,log=TRUE) )
}
#####################
## Mass-Extinction ##
#####################
# prior on number of changes for mass-extinction events
k <- length(timesMassExtinction)
lnp <- lnp + dpois(k,numExpectedMassExtinctions,log=TRUE)
if ( k > 0 ) {
# prior on change times
lnp <- lnp + sum( dunif(timesMassExtinction,0,AGE,log=TRUE) )
# prior on change values
lnp <- lnp + sum( dbeta(valuesMassExtinction,pMassExtinctionPriorShape1,pMassExtinctionPriorShape2,log=TRUE) )
}
return (lnp)
}
if ( use_tree_sample == TRUE ) {
AGE <- max( as.numeric( tess.branching.times( tree[[1]] )$age ) )
if ( length(tree) > 1 ) {
for ( i in 2:length(tree)) AGE <- min( c(AGE, max(as.numeric( tess.branching.times( tree[[i]] )$age) ) ) )
}
} else {
AGE <- max( as.numeric( tess.branching.times( tree )$age ) )
}
# these are the parameters we sample
posterior <- c()
kLambda <- c()
kMu <- c()
kPhi <- c()
kMassExtinction <- c()
speciationRateValues <- list()
speciationRateChangeTimes <- list()
extinctionRateValues <- list()
extinctionRateChangeTimes <- list()
fossilizationRateValues <- list()
fossilizationRateChangeTimes <- list()
survivalProbability <- list()
massExtinctionTime <- list()
cat("SpeciationRates\n",file="SpeciationRates.txt")
cat("SpeciationRateChanges\n",file="SpeciationRateChanges.txt")
cat("ExtinctionRates\n",file="ExtinctionRates.txt")
cat("ExtinctionRateChangeTimes\n",file="ExtinctionRateChanges.txt")
cat("FossilizationRates\n",file="FossilizationRates.txt")
cat("FossilizationRateChangeTimes\n",file="FossilizationRateChanges.txt")
cat("SurvivalProbabilities\n",file="SurvivalProbabilities.txt")
cat("MassExtinctionTimes\n",file="MassExtinctionTimes.txt")
cat(paste("Iteration","posterior","NumSpeciation","numExtinction","numMassExtinctions",sep="\t"),"\n",file="samples_numCategories.txt")
# the "random" starting values
lambda <- initialSpeciationRate
mu <- initialExtinctionRate
phi <- initialFossilizationRate
lambdaChangeTimes <- initialSpeciationRateChangeTime
muChangeTimes <- initialExtinctionRateChangeTime
phiChangeTimes <- initialFossilizationRateChangeTime
pMassExtinction <- pInitialMassExtinction
tMassExtinction <- tInitialMassExtinction
initialPP <- prior(lambdaChangeTimes,lambda,muChangeTimes,mu,phiChangeTimes,phi,tMassExtinction,pMassExtinction) + likelihood(lambda,lambdaChangeTimes,mu,muChangeTimes,phi,phiChangeTimes,tMassExtinction,pMassExtinction)
currentPP <- initialPP
## check if the prior and posterior gives valid probabilities
MAX_TRIES <- 1000
for ( i in 1:MAX_TRIES ) {
tmpPP <- prior(lambdaChangeTimes,lambda,muChangeTimes,mu,phiChangeTimes,phi,tMassExtinction,pMassExtinction) + likelihood(lambda,lambdaChangeTimes,mu,muChangeTimes,phi,phiChangeTimes,tMassExtinction,pMassExtinction)
if ( is.finite(tmpPP) ) break
lambdaChangeTimes <- c()
lambda <- rlnorm(1,speciationRatePriorMean,speciationRatePriorStDev)
muChangeTimes <- c()
mu <- rlnorm(1,extinctionRatePriorMean,extinctionRatePriorStDev)
phiChangeTimes <- c()
phi <- rlnorm(1,fossilizationRatePriorMean,fossilizationRatePriorStDev)
tMassExtinction <- c()
pMassExtinction <- c()
}
## MCMC proposal sizes
deltaLambdaTime <- 0.1
lambdaTimeAccepted <- 0
lambdaTimeTried <- 0
deltaLambdaValue <- 0.1
lambdaValueAccepted <- 0
lambdaValueTried <- 0
deltaMuTime <- 0.1
muTimeAccepted <- 0
muTimeTried <- 0
deltaMuValue <- 0.1
muValueAccepted <- 0
muValueTried <- 0
deltaPhiTime <- 0.1
phiTimeAccepted <- 0
phiTimeTried <- 0
deltaPhiValue <- 0.1
phiValueAccepted <- 0
phiValueTried <- 0
deltaMassExtinctionTime <- 0.1
massExtinctionTimeAccepted <- 0
massExtinctionTimeTried <- 0
deltaMassExtinctionValue <- 0.1
massExtinctionValueAccepted <- 0
massExtinctionValueTried <- 0
finished <- FALSE
SAMPLE <- FALSE
min.ess <- 0
startTime <- Sys.time()
i <- 1
if ( verbose ) {
cat("\nPerforming CoMET analysis.\n\n")
if( burn > 0 ) {
cat("Burning-in the chain ...\n")
} else {
cat("Running the chain ... \n")
}
cat("0--------25--------50--------75--------100\n")
bar <- txtProgressBar(style=1,width=42)
}
while ( !finished ) {
# change the time of an event
if ( length(lambdaChangeTimes) > 0 ) {
# compute index of value that will be changed
idx <- floor( runif(1,0,length(lambdaChangeTimes)) ) + 1
# store current value
t <- lambdaChangeTimes[idx]
# compute current posterior
oldPP <- currentPP
# propose new value
tPrime <- rnorm(1,t,deltaLambdaTime)
lambdaChangeTimes[idx] <- tPrime
# compute new posterior
if ( tPrime > 0 && tPrime < AGE ) {
newPP <- prior(lambdaChangeTimes,lambda,muChangeTimes,mu,phiChangeTimes,phi,tMassExtinction,pMassExtinction) + likelihood(lambda,lambdaChangeTimes,mu,muChangeTimes,phi,phiChangeTimes,tMassExtinction,pMassExtinction)
} else {
newPP <- -Inf
}
# accept/reject
if ( is.finite(newPP - oldPP) == FALSE || log( runif(1,0,1) ) > ( newPP - oldPP ) ) {
# reject
lambdaChangeTimes[idx] <- t
} else {
currentPP <- newPP
lambdaTimeAccepted <- lambdaTimeAccepted + 1
}
lambdaTimeTried <- lambdaTimeTried + 1
}
# change the value of an event
# compute index of value that will be changed
idx <- floor( runif(1,0,length(lambda)) ) + 1
# store current value
v <- lambda[idx]
# compute current posterior
oldPP <- currentPP
# propose new value
u <- rnorm(1,0,deltaLambdaValue)
scalingFactor <- exp( u )
vPrime <- v * scalingFactor
lambda[idx] <- vPrime
# compute the Hastings ratio
lnHastingsratio <- log( scalingFactor )
# compute new posterior
newPP <- prior(lambdaChangeTimes,lambda,muChangeTimes,mu,phiChangeTimes,phi,tMassExtinction,pMassExtinction) + likelihood(lambda,lambdaChangeTimes,mu,muChangeTimes,phi,phiChangeTimes,tMassExtinction,pMassExtinction)
# accept/reject
if ( is.finite(newPP - oldPP + lnHastingsratio) == FALSE || log( runif(1,0,1) ) > ( newPP - oldPP + lnHastingsratio ) ) {
# reject
lambda[idx] <- v
} else {
currentPP <- newPP
lambdaValueAccepted <- lambdaValueAccepted + 1
}
lambdaValueTried <- lambdaValueTried + 1
if ( estimateNumberSpeciationRateChanges == TRUE ) {
# We need to randomly pick a birth or death move
# Otherwise we might give birth and die every time
u <- runif(1,0,1)
if ( u > 0.5 ) {
# birth move
# compute current posterior
oldPP <- currentPP
# randomly pick a new time
t <- runif(1,0,AGE)
# randomly pick a new value
if ( lambda.hyper.prior.form == "lognormal" ) {
v <- rlnorm(1,speciationRatePriorMean,speciationRatePriorStDev)
} else if ( lambda.hyper.prior.form == "normal" ) {
v <- rnorm(1,speciationRatePriorMean,speciationRatePriorStDev)
} else if ( lambda.hyper.prior.form == "gamma" ) {
v <- rgamma(1,speciationRatePriorMean,speciationRatePriorStDev)
}
# construct the new parameters
lambdaChangeTimes[length(lambdaChangeTimes)+1] <- t
lambda[length(lambda)+1] <- v
# compute new posterior
newPP <- prior(lambdaChangeTimes,lambda,muChangeTimes,mu,phiChangeTimes,phi,tMassExtinction,pMassExtinction) + likelihood(lambda,lambdaChangeTimes,mu,muChangeTimes,phi,phiChangeTimes,tMassExtinction,pMassExtinction)
# compute proposal ratio
kPrime <- length(lambdaChangeTimes)
if ( lambda.hyper.prior.form == "lognormal" ) {
pr <- 1 / ( dunif(t,0,AGE) * dlnorm(v,speciationRatePriorMean,speciationRatePriorStDev) )
} else if ( lambda.hyper.prior.form == "normal" ) {
pr <- 1 / ( dunif(t,0,AGE) * dnorm(v,speciationRatePriorMean,speciationRatePriorStDev) )
} else if ( lambda.hyper.prior.form == "gamma" ) {
pr <- 1 / ( dunif(t,0,AGE) * dgamma(v,speciationRatePriorMean,speciationRatePriorStDev) )
}
# accept/reject
if ( is.finite(newPP - oldPP + log(pr)) == FALSE || log( runif(1,0,1) ) > ( newPP - oldPP + log( pr ) ) ) {
# reject
lambdaChangeTimes <- lambdaChangeTimes[-length(lambdaChangeTimes)]
lambda <- lambda[-length(lambda)]
} else {
currentPP <- newPP
}
} else {
# death move
if ( length(lambdaChangeTimes) > 0 ) {
# compute current posterior
oldPP <- currentPP
# randomly pick an index
idx <- floor( runif(1,0,length(lambdaChangeTimes)) ) + 1
# store the current value
t <- lambdaChangeTimes
v <- lambda
# construct the new parameters
lambdaChangeTimes <- lambdaChangeTimes[-idx]
lambda <- lambda[-(idx+1)]
# compute new posterior
newPP <- prior(lambdaChangeTimes,lambda,muChangeTimes,mu,phiChangeTimes,phi,tMassExtinction,pMassExtinction) + likelihood(lambda,lambdaChangeTimes,mu,muChangeTimes,phi,phiChangeTimes,tMassExtinction,pMassExtinction)
# compute proposal ratio
kPrime <- length(lambdaChangeTimes) + 1
if ( lambda.hyper.prior.form == "lognormal" ) {
pr2 <- dunif(t[idx],0,AGE) * dlnorm(v[idx+1],speciationRatePriorMean,speciationRatePriorStDev)
} else if ( lambda.hyper.prior.form == "normal" ) {
pr2 <- dunif(t[idx],0,AGE) * dnorm(v[idx+1],speciationRatePriorMean,speciationRatePriorStDev)
} else if ( lambda.hyper.prior.form == "gamma" ) {
pr2 <- dunif(t[idx],0,AGE) * dgamma(v[idx+1],speciationRatePriorMean,speciationRatePriorStDev)
}
# accept/reject
if ( is.finite(newPP - oldPP + log(pr2)) == FALSE || log( runif(1,0,1) ) > ( newPP - oldPP + log( pr2 ) ) ) {
# reject
lambdaChangeTimes <- t
lambda <- v
} else {
currentPP <- newPP
}
}
}
}
#####################
## Extinction Rate ##
#####################
# change the time of an event
if ( length(muChangeTimes) > 0 ) {
# compute index of value that will be changed
idx <- floor( runif(1,0,length(muChangeTimes)) ) + 1
# store current value
t <- muChangeTimes[idx]
# compute current posterior
oldPP <- currentPP
# propose new value
tPrime <- rnorm(1, t, deltaMuTime)
muChangeTimes[idx] <- tPrime
# compute new posterior
if ( tPrime > 0 && tPrime < AGE ) {
newPP <- prior(lambdaChangeTimes,lambda,muChangeTimes,mu,phiChangeTimes,phi,tMassExtinction,pMassExtinction) + likelihood(lambda,lambdaChangeTimes,mu,muChangeTimes,phi,phiChangeTimes,tMassExtinction,pMassExtinction)
} else {
newPP <- -Inf
}
# accept/reject
if ( is.finite(newPP - oldPP) == FALSE || log( runif(1,0,1) ) > ( newPP - oldPP ) ) {
# reject
muChangeTimes[idx] <- t
} else {
currentPP <- newPP
muTimeAccepted <- muTimeAccepted + 1
}
muTimeTried <- muTimeTried + 1
}
# change the value of an event
# compute index of value that will be changed
idx <- floor( runif(1,0,length(mu)) ) + 1
# store current value
v <- mu[idx]
# compute current posterior
oldPP <- currentPP
# propose new value
u <- rnorm(1,0,deltaMuValue)
scalingFactor <- exp( u )
vPrime <- v * scalingFactor
mu[idx] <- vPrime
# compute the Hastings ratio
lnHastingsratio <- log( scalingFactor )
# compute new posterior
newPP <- prior(lambdaChangeTimes,lambda,muChangeTimes,mu,phiChangeTimes,phi,tMassExtinction,pMassExtinction) + likelihood(lambda,lambdaChangeTimes,mu,muChangeTimes,phi,phiChangeTimes,tMassExtinction,pMassExtinction)
# accept/reject
if ( is.finite(newPP - oldPP + lnHastingsratio) == FALSE || log( runif(1,0,1) ) > ( newPP - oldPP + lnHastingsratio ) ) {
# reject
mu[idx] <- v
} else {
currentPP <- newPP
muValueAccepted <- muValueAccepted + 1
}
muValueTried <- muValueTried + 1
if ( estimateNumberExtinctionRateChanges == TRUE ) {
# We need to randomly pick a birth or death move
# Otherwise we might give birth and die every time
u <- runif(1,0,1)
if ( u > 0.5 ) {
# birth move
# compute current posterior
oldPP <- currentPP
# randomly pick a new time
t <- runif(1,0,AGE)
# randomly pick a new value
if ( mu.hyper.prior.form == "lognormal" ) {
v <- rlnorm(1,extinctionRatePriorMean,extinctionRatePriorStDev)
} else if ( mu.hyper.prior.form == "normal" ) {
v <- rnorm(1,extinctionRatePriorMean,extinctionRatePriorStDev)
} else if ( mu.hyper.prior.form == "gamma" ) {
v <- rgamma(1,extinctionRatePriorMean,extinctionRatePriorStDev)
}
# construct the new parameters
muChangeTimes[length(muChangeTimes)+1] <- t
mu[length(mu)+1] <- v
# compute new posterior
newPP <- prior(lambdaChangeTimes,lambda,muChangeTimes,mu,phiChangeTimes,phi,tMassExtinction,pMassExtinction) + likelihood(lambda,lambdaChangeTimes,mu,muChangeTimes,phi,phiChangeTimes,tMassExtinction,pMassExtinction)
# compute proposal ratio
kPrime <- length(muChangeTimes)
if ( mu.hyper.prior.form == "lognormal" ) {
pr <- 1 / ( dunif(t,0,AGE) * dlnorm(v,extinctionRatePriorMean,extinctionRatePriorStDev) )
} else if ( mu.hyper.prior.form == "normal" ) {
pr <- 1 / ( dunif(t,0,AGE) * dnorm(v,extinctionRatePriorMean,extinctionRatePriorStDev) )
} else if ( mu.hyper.prior.form == "gamma" ) {
pr <- 1 / ( dunif(t,0,AGE) * dgamma(v,extinctionRatePriorMean,extinctionRatePriorStDev) )
}
# accept/reject
if ( is.finite(newPP - oldPP + log(pr)) == FALSE || log( runif(1,0,1) ) > ( newPP - oldPP + log( pr ) ) ) {
# reject
muChangeTimes <- muChangeTimes[-length(muChangeTimes)]
mu <- mu[-length(mu)]
} else {
currentPP <- newPP
}
} else {
# death move
if ( length(muChangeTimes) > 0 ) {
# compute current posterior
oldPP <- currentPP
# randomly pick an index
idx <- floor( runif(1,0,length(muChangeTimes)) ) + 1
# store the current value
t <- muChangeTimes
v <- mu
# construct the new parameters
muChangeTimes <- muChangeTimes[-idx]
mu <- mu[-(idx+1)]
# compute new posterior
newPP <- prior(lambdaChangeTimes,lambda,muChangeTimes,mu,phiChangeTimes,phi,tMassExtinction,pMassExtinction) + likelihood(lambda,lambdaChangeTimes,mu,muChangeTimes,phi,phiChangeTimes,tMassExtinction,pMassExtinction)
# compute proposal ratio
kPrime <- length(muChangeTimes) + 1
if ( mu.hyper.prior.form == "lognormal" ) {
pr2 <- dunif(t[idx],0,AGE) * dlnorm(v[idx+1],extinctionRatePriorMean,extinctionRatePriorStDev)
} else if ( mu.hyper.prior.form == "normal" ) {
pr2 <- dunif(t[idx],0,AGE) * dnorm(v[idx+1],extinctionRatePriorMean,extinctionRatePriorStDev)
} else if ( mu.hyper.prior.form == "gamma" ) {
pr2 <- dunif(t[idx],0,AGE) * dgamma(v[idx+1],extinctionRatePriorMean,extinctionRatePriorStDev)
}
# accept/reject
if ( is.finite(newPP - oldPP + log(pr2)) == FALSE || log( runif(1,0,1) ) > ( newPP - oldPP + log( pr2 ) ) ) {
# reject
muChangeTimes <- t
mu <- v
} else {
currentPP <- newPP
}
}
}
}
########################
## Fossilization Rate ##
########################
# change the time of an event
if ( length(phiChangeTimes) > 0 ) {
# compute index of value that will be changed
idx <- floor( runif(1,0,length(phiChangeTimes)) ) + 1
# store current value
t <- phiChangeTimes[idx]
# compute current posterior
oldPP <- currentPP
# propose new value
tPrime <- rnorm(1, t, deltaPhiTime)
phiChangeTimes[idx] <- tPrime
# compute new posterior
if ( tPrime > 0 && tPrime < AGE ) {
newPP <- prior(lambdaChangeTimes,lambda,muChangeTimes,mu,phiChangeTimes,phi,tMassExtinction,pMassExtinction) + likelihood(lambda,lambdaChangeTimes,mu,muChangeTimes,phi,phiChangeTimes,tMassExtinction,pMassExtinction)
} else {
newPP <- -Inf
}
# accept/reject
if ( is.finite(newPP - oldPP) == FALSE || log( runif(1,0,1) ) > ( newPP - oldPP ) ) {
# reject
phiChangeTimes[idx] <- t
} else {
currentPP <- newPP
phiTimeAccepted <- phiTimeAccepted + 1
}
phiTimeTried <- phiTimeTried + 1
}
# change the value of an event
# compute index of value that will be changed
idx <- floor( runif(1,0,length(phi)) ) + 1
# store current value
v <- phi[idx]
# compute current posterior
oldPP <- currentPP
# propose new value
u <- rnorm(1,0,deltaPhiValue)
scalingFactor <- exp( u )
vPrime <- v * scalingFactor
phi[idx] <- vPrime
# compute the Hastings ratio
lnHastingsratio <- log( scalingFactor )
# compute new posterior
newPP <- prior(lambdaChangeTimes,lambda,muChangeTimes,mu,phiChangeTimes,phi,tMassExtinction,pMassExtinction) + likelihood(lambda,lambdaChangeTimes,mu,muChangeTimes,phi,phiChangeTimes,tMassExtinction,pMassExtinction)
# accept/reject
if ( is.finite(newPP - oldPP + lnHastingsratio) == FALSE || log( runif(1,0,1) ) > ( newPP - oldPP + lnHastingsratio ) ) {
# reject
phi[idx] <- v
} else {
currentPP <- newPP
phiValueAccepted <- phiValueAccepted + 1
}
phiValueTried <- phiValueTried + 1
if ( estimateNumberFossilizationRateChanges == TRUE ) {
# We need to randomly pick a birth or death move
# Otherwise we might give birth and die every time
u <- runif(1,0,1)
if ( u > 0.5 ) {
# birth move
# compute current posterior
oldPP <- currentPP
# randomly pick a new time
t <- runif(1,0,AGE)
# randomly pick a new value
if ( phi.hyper.prior.form == "lognormal" ) {
v <- rlnorm(1,fossilizationRatePriorMean,fossilizationRatePriorStDev)
} else if ( phi.hyper.prior.form == "normal" ) {
v <- rnorm(1,fossilizationRatePriorMean,fossilizationRatePriorStDev)
} else if ( phi.hyper.prior.form == "gamma" ) {
v <- rgamma(1,fossilizationRatePriorMean,fossilizationRatePriorStDev)
}
# construct the new parameters
phiChangeTimes[length(phiChangeTimes)+1] <- t
phi[length(phi)+1] <- v
# compute new posterior
newPP <- prior(lambdaChangeTimes,lambda,muChangeTimes,mu,phiChangeTimes,phi,tMassExtinction,pMassExtinction) + likelihood(lambda,lambdaChangeTimes,mu,muChangeTimes,phi,phiChangeTimes,tMassExtinction,pMassExtinction)
# compute proposal ratio
kPrime <- length(phiChangeTimes)
if ( phi.hyper.prior.form == "lognormal" ) {
pr <- 1 / ( dunif(t,0,AGE) * dlnorm(v,fossilizationRatePriorMean,fossilizationRatePriorStDev) )
} else if ( phi.hyper.prior.form == "normal" ) {
pr <- 1 / ( dunif(t,0,AGE) * dnorm(v,fossilizationRatePriorMean,fossilizationRatePriorStDev) )
} else if ( phi.hyper.prior.form == "gamma" ) {
pr <- 1 / ( dunif(t,0,AGE) * dgamma(v,fossilizationRatePriorMean,fossilizationRatePriorStDev) )
}
# accept/reject
if ( is.finite(newPP - oldPP + log(pr)) == FALSE || log( runif(1,0,1) ) > ( newPP - oldPP + log( pr ) ) ) {
# reject
phiChangeTimes <- phiChangeTimes[-length(phiChangeTimes)]
phi <- phi[-length(phi)]
} else {
currentPP <- newPP
}
} else {
# death move
if ( length(phiChangeTimes) > 0 ) {
# compute current posterior
oldPP <- currentPP
# randomly pick an index
idx <- floor( runif(1,0,length(phiChangeTimes)) ) + 1
# store the current value
t <- phiChangeTimes
v <- phi
# construct the new parameters
phiChangeTimes <- phiChangeTimes[-idx]
phi <- phi[-(idx+1)]
# compute new posterior
newPP <- prior(lambdaChangeTimes,lambda,muChangeTimes,mu,phiChangeTimes,phi,tMassExtinction,pMassExtinction) + likelihood(lambda,lambdaChangeTimes,mu,muChangeTimes,phi,phiChangeTimes,tMassExtinction,pMassExtinction)
# compute proposal ratio
kPrime <- length(phiChangeTimes) + 1
if ( phi.hyper.prior.form == "lognormal" ) {
pr2 <- dunif(t[idx],0,AGE) * dlnorm(v[idx+1],fossilizationRatePriorMean,fossilizationRatePriorStDev)
} else if ( phi.hyper.prior.form == "normal" ) {
pr2 <- dunif(t[idx],0,AGE) * dnorm(v[idx+1],fossilizationRatePriorMean,fossilizationRatePriorStDev)
} else if ( phi.hyper.prior.form == "gamma" ) {
pr2 <- dunif(t[idx],0,AGE) * dgamma(v[idx+1],fossilizationRatePriorMean,fossilizationRatePriorStDev)
}
# accept/reject
if ( is.finite(newPP - oldPP + log(pr2)) == FALSE || log( runif(1,0,1) ) > ( newPP - oldPP + log( pr2 ) ) ) {
# reject
phiChangeTimes <- t
phi <- v
} else {
currentPP <- newPP
}
}
}
}
#####################
## Mass-Extinction ##
#####################
# change the time of an event
if ( estimateMassExtinctionTimes == TRUE && length(tMassExtinction) > 0 ) {
# compute index of value that will be changed
idx <- floor( runif(1,0,length(tMassExtinction)) ) + 1
# store current value
t <- tMassExtinction[idx]
# compute current posterior
oldPP <- currentPP
# propose new value
tPrime <- rnorm(1, t, deltaMassExtinctionTime)
tMassExtinction[idx] <- tPrime
# compute new posterior
if ( tPrime > 0 && tPrime < AGE ) {
newPP <- prior(lambdaChangeTimes,lambda,muChangeTimes,mu,phiChangeTimes,phi,tMassExtinction,pMassExtinction) + likelihood(lambda,lambdaChangeTimes,mu,muChangeTimes,phi,phiChangeTimes,tMassExtinction,pMassExtinction)
} else {
newPP <- -Inf
}
# accept/reject
if ( is.finite(newPP - oldPP) == FALSE || log( runif(1,0,1) ) > ( newPP - oldPP ) ) {
# reject
tMassExtinction[idx] <- t
} else {
currentPP <- newPP
massExtinctionTimeAccepted <- massExtinctionTimeAccepted + 1
}
massExtinctionTimeTried <- massExtinctionTimeTried + 1
}
# change the value of an event
if ( length(pMassExtinction) > 0 ) {
# compute index of value that will be changed
idx <- floor( runif(1,0,length(pMassExtinction)) ) + 1
# store current value
v <- pMassExtinction[idx]
# compute current posterior
oldPP <- currentPP
# propose new value
vPrime <- rnorm(1, v, deltaMassExtinctionValue)
pMassExtinction[idx] <- vPrime
# compute new posterior
if ( vPrime > 0 && vPrime < 1 ) {
newPP <- prior(lambdaChangeTimes,lambda,muChangeTimes,mu,phiChangeTimes,phi,tMassExtinction,pMassExtinction) + likelihood(lambda,lambdaChangeTimes,mu,muChangeTimes,phi,phiChangeTimes,tMassExtinction,pMassExtinction)
} else {
newPP <- -Inf
}
# accept/reject
if ( is.finite(newPP - oldPP) == FALSE || log( runif(1,0,1) ) > ( newPP - oldPP ) ) {
# reject
pMassExtinction[idx] <- v
} else {
currentPP <- newPP
massExtinctionValueAccepted <- massExtinctionValueAccepted + 1
}
massExtinctionValueTried <- massExtinctionValueTried + 1
}
if ( estimateNumberMassExtinctions == TRUE ) {
# We need to randomly pick a birth or death move
# Otherwise we might give birth and die every time
u <- runif(1,0,1)
if ( u > 0.5 ) {
# birth move
# compute current posterior
oldPP <- currentPP
# randomly pick a new time
t <- runif(1,0,AGE)
# randomly pick a new value
v <- rbeta(1,pMassExtinctionPriorShape1,pMassExtinctionPriorShape2)
# construct the new parameters
tMassExtinction[length(tMassExtinction)+1] <- t
pMassExtinction[length(pMassExtinction)+1] <- v
# compute new posterior
newPP <- prior(lambdaChangeTimes,lambda,muChangeTimes,mu,phiChangeTimes,phi,tMassExtinction,pMassExtinction) + likelihood(lambda,lambdaChangeTimes,mu,muChangeTimes,phi,phiChangeTimes,tMassExtinction,pMassExtinction)
# compute proposal ratio
kPrime <- length(pMassExtinction)
pr <- 1 / ( dunif(t,0,AGE) * dbeta(v,pMassExtinctionPriorShape1,pMassExtinctionPriorShape2) )
# accept/reject
if ( is.finite(newPP - oldPP + log(pr)) == FALSE || log( runif(1,0,1) ) > ( newPP - oldPP + log( pr ) ) ) {
# reject
tMassExtinction <- tMassExtinction[-length(tMassExtinction)]
pMassExtinction <- pMassExtinction[-length(pMassExtinction)]
} else {
currentPP <- newPP
}
} else {
# death move
if ( length(tMassExtinction) > 0 ) {
# compute current posterior
oldPP <- currentPP
# randomly pick an index
idx <- floor( runif(1,0,length(tMassExtinction)) ) + 1
# store the current value
t <- tMassExtinction
v <- pMassExtinction
# construct the new parameters
tMassExtinction <- tMassExtinction[-idx]
pMassExtinction <- pMassExtinction[-idx]
# compute new posterior
newPP <- prior(lambdaChangeTimes,lambda,muChangeTimes,mu,phiChangeTimes,phi,tMassExtinction,pMassExtinction) + likelihood(lambda,lambdaChangeTimes,mu,muChangeTimes,phi,phiChangeTimes,tMassExtinction,pMassExtinction)
# compute proposal ratio
kPrime <- length(tMassExtinction) + 1
pr2 <- dunif(t[idx],0,AGE) * dbeta(v[idx],pMassExtinctionPriorShape1,pMassExtinctionPriorShape2)
# accept/reject
if ( is.finite(newPP - oldPP + log(pr2)) == FALSE || log( runif(1,0,1) ) > ( newPP - oldPP + log( pr2 ) ) ) {
# reject
tMassExtinction <- t
pMassExtinction <- v
} else {
currentPP <- newPP
}
}
}
}
#################
## Auto-Tuning ##
#################
# if we use adaptive MCMC we might have to modify our parameters
if ( ADAPTIVE == TRUE & SAMPLE == FALSE ) {
if ( i %% OPTIMIZATION_FREQUENCY == 0 ) {
# delta for time of lambda change events
if ( lambdaTimeTried > 0 ) {
rate <- lambdaTimeAccepted / lambdaTimeTried
if ( rate > 0.44 ) {
deltaLambdaTime <- deltaLambdaTime * (1.0 + ((rate-0.44)/0.56) )
} else {
deltaLambdaTime <- deltaLambdaTime / (2.0 - rate/0.44 )
}
lambdaTimeTried <- 0
lambdaTimeAccepted <- 0
}
# delta for value of lambda
rate <- lambdaValueAccepted / lambdaValueTried
if ( rate > 0.44 ) {
deltaLambdaValue <- deltaLambdaValue * (1.0 + ((rate-0.44)/0.56) )
} else {
deltaLambdaValue <- deltaLambdaValue / (2.0 - rate/0.44 )
}
lambdaValueTried <- 0
lambdaValueAccepted <- 0
# delta for time of mu change events
if ( muTimeTried > 0 ) {
rate <- muTimeAccepted / muTimeTried
if ( rate > 0.44 ) {
deltaMuTime <- deltaMuTime * (1.0 + ((rate-0.44)/0.56) )
} else {
deltaMuTime <- deltaMuTime / (2.0 - rate/0.44 )
}
muTimeTried <- 0
muTimeAccepted <- 0
}
# delta for values of mu
rate <- muValueAccepted / muValueTried
if ( rate > 0.44 ) {
deltaMuValue <- deltaMuValue * (1.0 + ((rate-0.44)/0.56) )
} else {
deltaMuValue <- deltaMuValue / (2.0 - rate/0.44 )
}
muValueTried <- 0
muValueAccepted <- 0
# delta for time of phi change events
if ( phiTimeTried > 0 ) {
rate <- phiTimeAccepted / phiTimeTried
if ( rate > 0.44 ) {
deltaPhiTime <- deltaPhiTime * (1.0 + ((rate-0.44)/0.56) )
} else {
deltaPhiTime <- deltaPhiTime / (2.0 - rate/0.44 )
}
phiTimeTried <- 0
phiTimeAccepted <- 0
}
# delta for values of phi
rate <- phiValueAccepted / phiValueTried
if ( rate > 0.44 ) {
deltaPhiValue <- deltaPhiValue * (1.0 + ((rate-0.44)/0.56) )
} else {
deltaPhiValue <- deltaPhiValue / (2.0 - rate/0.44 )
}
phiValueTried <- 0
phiValueAccepted <- 0
# delta for time of massExtinction change events
if ( massExtinctionTimeTried > 0 ) {
rate <- massExtinctionTimeAccepted / massExtinctionTimeTried
if ( rate > 0.44 ) {
deltaMassExtinctionTime <- deltaMassExtinctionTime * (1.0 + ((rate-0.44)/0.56) )
} else {
deltaMassExtinctionTime <- deltaMassExtinctionTime / (2.0 - rate/0.44 )
}
massExtinctionTimeTried <- 0
massExtinctionTimeAccepted <- 0
}
# delta for values of mu
if ( massExtinctionValueTried > 0 ) {
rate <- massExtinctionValueAccepted / massExtinctionValueTried
if ( rate > 0.44 ) {
deltaMassExtinctionValue <- deltaMassExtinctionValue * (1.0 + ((rate-0.44)/0.56) )
} else {
deltaMassExtinctionValue <- deltaMassExtinctionValue / (2.0 - rate/0.44 )
}
massExtinctionValueTried <- 0
massExtinctionValueAccepted <- 0
}
}
}
# store the values
if ( i %% THINNING == 0 & SAMPLE ) {
sampleIndex <- length(kLambda)+1
kLambda[sampleIndex] <- length(lambda)
speciationRateValues[[sampleIndex]] <- lambda
speciationRateChangeTimes[[sampleIndex]] <- lambdaChangeTimes
kMu[sampleIndex] <- length(mu)
extinctionRateValues[[sampleIndex]] <- mu
extinctionRateChangeTimes[[sampleIndex]] <- muChangeTimes
kPhi[sampleIndex] <- length(phi)
fossilizationRateValues[[sampleIndex]] <- phi
fossilizationRateChangeTimes[[sampleIndex]] <- phiChangeTimes
kMassExtinction[sampleIndex] <- length(pMassExtinction)
survivalProbability[[sampleIndex]] <- pMassExtinction
massExtinctionTime[[sampleIndex]] <- tMassExtinction
posterior[sampleIndex] <- currentPP
# write.table(speciationRateValues,"SpeciationRates.txt",sep="\t")
cat(lambda,sep="\t",file="SpeciationRates.txt",append=TRUE)
cat("\n",file="SpeciationRates.txt",append=TRUE)
# write.table(speciationRateChangeTimes,"SpeciationRateChanges.txt",sep="\t")
cat(lambdaChangeTimes,sep="\t",file="SpeciationRateChanges.txt",append=TRUE)
cat("\n",file="SpeciationRateChanges.txt",append=TRUE)
# write.table(extinctionRateValues,"ExtinctionRates.txt",sep="\t")
cat(mu,sep="\t",file="ExtinctionRates.txt",append=TRUE)
cat("\n",file="ExtinctionRates.txt",append=TRUE)
# write.table(extinctionRateChangeTimes,"ExtinctionRateChanges.txt",sep="\t")
cat(muChangeTimes,sep="\t",file="ExtinctionRateChanges.txt",append=TRUE)
cat("\n",file="ExtinctionRateChanges.txt",append=TRUE)
# write.table(fossilizationRateValues,"FossilizationRates.txt",sep="\t")
cat(phi,sep="\t",file="FossilizationRates.txt",append=TRUE)
cat("\n",file="FossilizationRates.txt",append=TRUE)
# write.table(fossilizationRateChangeTimes,"FossilizationRateChanges.txt",sep="\t")
cat(phiChangeTimes,sep="\t",file="FossilizationRateChanges.txt",append=TRUE)
cat("\n",file="FossilizationRateChanges.txt",append=TRUE)
# write.table(survivalProbability,"SurvivalProbabilities.txt",sep="\t")
cat(pMassExtinction,sep="\t",file="SurvivalProbabilities.txt",append=TRUE)
cat("\n",file="SurvivalProbabilities.txt",append=TRUE)
# write.table(massExtinctionTime,"MassExtinctionTimes.txt",sep="\t")
cat(tMassExtinction,sep="\t",file="MassExtinctionTimes.txt",append=TRUE)
cat("\n",file="MassExtinctionTimes.txt",append=TRUE)
nSamples <- sampleIndex
write.table(list(Iteration=nSamples*THINNING,posterior=posterior[nSamples],NumSpeciation=kLambda[nSamples],numExtinction=kMu[nSamples],numFossilization=kPhi[nSamples],numMassExtinctions=kMassExtinction[nSamples]),"samples_numCategories.txt", sep="\t", row.names = FALSE, append=TRUE, quote = FALSE, col.names = FALSE)
}
# Progress by convergence
if ( i %% CONVERGENCE_FREQUENCY == 0 & SAMPLE & length(posterior) > 2 ) {
samples <- length(posterior)
if ( estimateNumberSpeciationRateChanges == TRUE ) {
spectralDensityLambda <- spectrum0.ar(kLambda)$spec
if ( spectralDensityLambda > 0 ) {
ess_kLambda <- (var(kLambda) * samples/spectralDensityLambda)
} else {
ess_kLambda <- Inf
}
} else {
ess_kLambda <- Inf
}
if ( estimateNumberExtinctionRateChanges == TRUE ) {
spectralDensityMu <- spectrum0.ar(kMu)$spec
if ( spectralDensityMu > 0 ) {
ess_kMu <- (var(kMu) * samples/spectralDensityMu)
} else {
ess_kMu <- Inf
}
} else {
ess_kMu <- Inf
}
if ( estimateNumberFossilizationRateChanges == TRUE ) {
spectralDensityPhi <- spectrum0.ar(kPhi)$spec
if ( spectralDensityPhi > 0 ) {
ess_kPhi <- (var(kPhi) * samples/spectralDensityPhi)
} else {
ess_kPhi <- Inf
}
} else {
ess_kPhi <- Inf
}
spec <- c()
for ( j in 1:length(speciationRateValues) ) {
spec[j] <- speciationRateValues[[j]][1]
}
spectralDensity <- spectrum0.ar(spec)$spec
if ( spectralDensity > 0 ) {
ess_spec <- (var(spec) * samples/spectralDensity)
} else {
ess_spec <- Inf
}
ext <- c()
for ( j in 1:length(extinctionRateValues) ) {
ext[j] <- extinctionRateValues[[j]][1]
}
spectralDensity <- spectrum0.ar(ext)$spec
if ( spectralDensity > 0 ) {
ess_ext <- (var(ext) * samples/spectralDensity)
} else {
ess_ext <- Inf
}
fos <- c()
for ( j in 1:length(fossilizationRateValues) ) {
fos[j] <- fossilizationRateValues[[j]][1]
}
spectralDensity <- spectrum0.ar(fos)$spec
if ( spectralDensity > 0 ) {
ess_fos <- (var(fos) * samples/spectralDensity)
} else {
ess_fos <- Inf
}
if ( estimateNumberMassExtinctions == TRUE ) {
ess_kMassExtinctions <- (var(kMassExtinction) * samples/spectrum0.ar(kMassExtinction)$spec)
} else {
ess_kMassExtinctions <- Inf
}
min.ess <- min(c(ess_kLambda,ess_kMu,ess_kPhi,ess_spec,ess_ext,ess_fos,ess_kMassExtinctions))
}
if ( SAMPLE ) {
# Progress by iteration
iteration.progress <- i / MAX_ITERATIONS
# Progress by time
time.progress <- ( (Sys.time() - startTime) / MAX_TIME )
# Progress by ESS
ess.progress <- min.ess / MIN_ESS
# Max progress
max.progress <- max(c(iteration.progress,time.progress,ess.progress),na.rm=TRUE)
if( length(max.progress) == 0 ){
max.progress <- 0
}
if ( verbose ) {
setTxtProgressBar(bar,pmin(max.progress,1))
}
if( max.progress >= 1 ) {
finished <- TRUE
}
} else {
if ( verbose ){
setTxtProgressBar(bar,pmin(1,i/burn))
}
}
i <- i + 1
if ( i == burn & !SAMPLE ) {
startTime <- Sys.time()
i <- 1
SAMPLE <- TRUE
if ( verbose ) {
cat("\n\nRunning the chain ... \n")
cat("0--------25--------50--------75--------100\n")
bar <- txtProgressBar(style=1,width=42)
}
}
}
cat("\n")
setwd(file.path(orgDir))
}
hoehna/TESS documentation built on Feb. 3, 2022, 5:59 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions tess.analysis.efbd
Documented in tess.analysis.efbd
#' tess.analysis.efbd
#'
#' @description Running an analysis on a given tree and estimating the diversification rates including rate-shifts and mass-extinction events.
#'
#' @param tree the phylogeny
#' @param initialSpeciationRate The initial value of the speciation rate when the MCMC is started.
#' @param initialExtinctionRate The initial value of the extinction rate when the MCMC is started.
#' @param initialFossilizationRate The initial value of the fossilization rate when the MCMC is started.
#' @param empiricalHyperPriors .
#' @param empiricalHyperPriorInflation .
#' @param empiricalHyperPriorForm .
#' @param speciationRatePriorMean mean parameter for the lognormal prior on lambda
#' @param speciationRatePriorStDev standard deviation for the lognormal prior on lambda
#' @param extinctionRatePriorMean mean parameter for the lognormal prior on mu
#' @param extinctionRatePriorStDev standard deviation for the lognormal prior on mu
#' @param fossilizationRatePriorMean mean parameter for the lognormal prior on phi
#' @param fossilizationRatePriorStDev standard deviation for the lognormal prior on phi
#' @param initialSpeciationRateChangeTime The initial value of the time points when speciation rate-shifts occur.
#' @param initialExtinctionRateChangeTime The initial value of the time points when extinction rate-shifts occur.
#' @param initialFossilizationRateChangeTime The initial value of the time points when fossilization rate-shifts occur.
#' @param estimateNumberSpeciationRateChanges Do we estimate the number of rate shifts?
#' @param estimateNumberExtinctionRateChanges Do we estimate the number of rate shifts?
#' @param estimateNumberFossilizationRateChanges Do we estimate the number of rate shifts?
#' @param numExpectedRateChanges Expected number of rate changes which follow a Poisson process.
#' @param samplingProbability probability of uniform sampling at present
#' @param samplingStrategy Which strategy was used to obtain the samples (taxa). Options are: uniform|diversified|age
#' @param missingSpecies The number of species missed which originated in a given time interval (empirical taxon sampling)
#' @param timesMissingSpecies The times intervals of the missing species (empirical taxon sampling)
#' @param tInitialMassExtinction The initial value of the vector of times of the mass-extinction events.
#' @param pInitialMassExtinction The initial value of the vector of survival probabilities of the mass-extinction events.
#' @param pMassExtinctionPriorShape1 The alpha (first shape) parameter of the Beta prior distribution for the survival probability of a mass-extinction event.
#' @param pMassExtinctionPriorShape2 The beta (second shape) parameter of the Beta prior distribution for the survival probability of a mass-extinction event.
#' @param estimateMassExtinctionTimes Do we estimate the times of mass-extinction events?
#' @param numExpectedMassExtinctions Expected number of mass-extinction events which follow a Poisson process.
#' @param estimateNumberMassExtinctions Do we estimate the number of mass-extinction events?
#' @param MRCA does the tree start at the mrca?
#' @param CONDITION do we condition the process on nothing|survival|taxa?
#' @param BURNIN number of starting iterations to be dropped
#' @param MAX_ITERATIONS The maximum number of iteration of the MCMC.
#' @param THINNING The frequency how often samples are recorded during the MCMC.
#' @param OPTIMIZATION_FREQUENCY The frequency how often the MCMC moves are optimized
#' @param CONVERGENCE_FREQUENCY The frequency how often we check for convergence?
#' @param MAX_TIME The maximum time the MCMC is allowed to run in seconds.
#' @param MIN_ESS The minimum number of effective samples (ESS) to assume convergence.
#' @param ADAPTIVE Do we use auto-tuning of the MCMC moves?
#' @param dir The subdirectory in which the output will be stored.
#' @param priorOnly Do we sample from the prior only?
#' @param verbose Do we want to print the progress of the MCMC?
#'
#' @return nothing. outputs are written to the file system
#' @export
#'
#' @examples
#' # we load the conifers as the test data set
#' data(conifers)
#'
#' # for the conifers we know what the total number of species is
#' total <- 630
#' # thus, we can compute what the sampling fraction is
#' rho <- (conifers$Nnode+1)/total
#'
#'
#' # next, we specify the prior mean and standard deviation
#' # for the speciation and extinction rate
#' mu_lambda = 0.15
#' std_lambda = 0.02
#' mu_mu = 0.09
#' std_mu = 0.02
#' mu_phi = 0.10
#' std_phi = 0.02
#'
#'
#'
#' # now we can run the entire analysis.
#' # note that a full analyses should be run much longer
#' tess.analysis.efbd( tree=conifers,
#' initialSpeciationRate = exp(mu_lambda),
#' initialExtinctionRate = exp(mu_mu),
#' initialFossilizationRate = exp(mu_phi),
#' empiricalHyperPriors = FALSE,
#' speciationRatePriorMean = mu_lambda,
#' speciationRatePriorStDev = std_lambda,
#' extinctionRatePriorMean = mu_mu,
#' extinctionRatePriorStDev = std_mu,
#' fossilizationRatePriorMean = mu_phi,
#' fossilizationRatePriorStDev = std_phi,
#' numExpectedRateChanges = 2,
#' samplingProbability = rho,
#' numExpectedMassExtinctions = 2,
#' BURNIN = 100,
#' MAX_ITERATIONS = 200,
#' THINNING = 10,
#' dir = "analysis_conifer")
tess.analysis.efbd <- function( tree,
initialSpeciationRate,
initialExtinctionRate,
initialFossilizationRate,
empiricalHyperPriors = TRUE,
empiricalHyperPriorInflation = 10.0,
empiricalHyperPriorForm = c("lognormal","normal","gamma"),
speciationRatePriorMean = 0.0,
speciationRatePriorStDev = 1.0,
extinctionRatePriorMean = 0.0,
extinctionRatePriorStDev = 1.0,
fossilizationRatePriorMean = 0.0,
fossilizationRatePriorStDev = 1.0,
initialSpeciationRateChangeTime = c(),
initialExtinctionRateChangeTime = c(),
initialFossilizationRateChangeTime = c(),
estimateNumberSpeciationRateChanges = TRUE,
estimateNumberExtinctionRateChanges = TRUE,
estimateNumberFossilizationRateChanges = TRUE,
numExpectedRateChanges = 2,
samplingProbability = 1,
samplingStrategy = "uniform",
missingSpecies = c(),
timesMissingSpecies = c(),
tInitialMassExtinction = c(),
pInitialMassExtinction = c(),
pMassExtinctionPriorShape1 = 5,
pMassExtinctionPriorShape2 = 95,
estimateMassExtinctionTimes = TRUE,
numExpectedMassExtinctions = 2,
estimateNumberMassExtinctions = TRUE,
MRCA = TRUE,
CONDITION = "survival",
BURNIN = 10000,
MAX_ITERATIONS = 200000,
THINNING = 100,
OPTIMIZATION_FREQUENCY = 500,
CONVERGENCE_FREQUENCY = 1000,
MAX_TIME = Inf, MIN_ESS = 500,
ADAPTIVE = TRUE,
dir = "" ,
priorOnly = FALSE,
verbose = TRUE) {
USE_EBDSTP = FALSE
orgDir <- getwd()
if ( dir != "" ) {
if (file.exists(dir)){
setwd(file.path(dir))
} else {
dir.create(file.path(dir))
setwd(file.path(dir))
}
}
if ( BURNIN > 0 & BURNIN < 1 ) {
burn <- MAX_ITERATIONS * BURNIN
} else {
burn <- BURNIN
}
MAX_TRIES <- 1000
if ( class(tree) == "phylo" ) {
use_tree_sample = FALSE
times <- tess.branching.times(tree)
} else if ( class(tree) == "multiPhylo" ) {
use_tree_sample = TRUE
times <- list()
for (i in 1:length(tree)) {
times[[i]] <- tess.branching.times(tree[[i]])
}
NUM_SAMPLED_TREES <- length(tree)
}
if ( CONDITION != "time" && CONDITION != "survival" && CONDITION != "taxa" ) {
stop("Wrong choice of argument for \"CONDITION\". Possible option are time|survival|taxa.")
}
if ( samplingStrategy != "uniform" && samplingStrategy != "diversified") {
stop("Wrong choice of argument for \"samplingStrategy\". Possible option are uniform|diversified.")
}
write.nexus(tree,file="input_tree.tre")
if ( empiricalHyperPriors == TRUE ) {
if( verbose ) {
cat("Estimating empirical hyper-parameters.\n\n")
}
empirical.prior.likelihood <- function(params) {
b <- params[1]
d <- params[2]
f <- params[3]
if ( use_tree_sample == TRUE ) {
ln_probs <- c()
max_prob <- -Inf
for (i in 1:NUM_SAMPLED_TREES) {
if ( USE_EBDSTP == FALSE ) {
ln_probs[i] <- tess.likelihood.efbd(nodes=times[[i]],
lambda=b,
mu=d,
phi=f,
samplingStrategy = samplingStrategy,
samplingProbability = samplingProbability,
MRCA = MRCA,
CONDITION = CONDITION,
log=TRUE)
} else {
ln_probs[i] <- tess.likelihood.ebdstp(nodes=times[[i]],
lambda=b,
mu=d,
phi=f,
r=0.0,
samplingStrategy = samplingStrategy,
samplingProbability = samplingProbability,
MRCA = MRCA,
CONDITION = CONDITION,
log=TRUE)
}
if ( max_prob < ln_probs[i] ) max_prob <- ln_probs[i]
}
sum <- 0.0
for (i in 1:NUM_SAMPLED_TREES) {
sum <- sum + exp( ln_probs[i] - max_prob ) / NUM_SAMPLED_TREES
}
lnl <- log( sum ) + max_prob
} else {
if ( USE_EBDSTP == FALSE ) {
lnl <- tess.likelihood.efbd(nodes=times,
lambda=b,
mu=d,
phi=f,
samplingStrategy = samplingStrategy,
samplingProbability = samplingProbability,
MRCA = MRCA,
CONDITION = CONDITION,
log=TRUE)
} else {
lnl <- tess.likelihood.ebdstp(nodes=times,
lambda=b,
mu=d,
phi=f,
r=0.0,
samplingStrategy = samplingStrategy,
samplingProbability = samplingProbability,
MRCA = MRCA,
CONDITION = CONDITION,
log=TRUE)
}
}
return(lnl)
}
if ( use_tree_sample == TRUE ) {
prior.speciation <- function(x) { dexp(x,rate=max(times[[1]]$age)/log(length(times[[1]]$age)),log=TRUE) }
} else {
prior.speciation <- function(x) { dexp(x,rate=max(times$age)/log(length(times$age)),log=TRUE) }
}
prior.speciation <- function(x) { dexp(x,rate=1.0,log=TRUE) }
prior.extinction <- function(x) { dexp(x,rate=1.0,log=TRUE) }
prior.fossilization <- function(x) { dexp(x,rate=1.0,log=TRUE) }
priors <- c("speciation"=prior.speciation,"extinction"=prior.extinction,"fossilization"=prior.fossilization)
tries <- 1
repeat {
starting.values <- runif(3,0,1)
starting.lnl <- empirical.prior.likelihood(starting.values)
if ( is.finite(starting.lnl) || tries >= MAX_TRIES ) {
break
}
tries <- tries + 1
}
if ( tries >= MAX_TRIES ) {
stop("Could not find valid starting values after ",tries," attempts.\n")
}
samples <- tess.mcmc( empirical.prior.likelihood,
priors,
starting.values,
logTransforms=c(TRUE,TRUE,TRUE),
delta=c(0.1,0.1,0.1),
iterations=20000,
burnin=2000,
verbose=verbose)
samples.lambda <- samples[,"speciation"]
m.lambda <- mean(samples.lambda)
v.lambda <- empiricalHyperPriorInflation * var(samples.lambda)
samples.mu <- samples[,"extinction"]
m.mu <- mean(samples.mu)
v.mu <- empiricalHyperPriorInflation * var(samples.mu)
samples.phi <- samples[,"fossilization"]
m.phi <- mean(samples.phi)
v.phi <- empiricalHyperPriorInflation * var(samples.phi)
lambda.hyperprior.fits <- c("lognormal"=Inf,"normal"=Inf,"gamma"=Inf)
mu.hyperprior.fits <- c("lognormal"=Inf,"normal"=Inf,"gamma"=Inf)
phi.hyperprior.fits <- c("lognormal"=Inf,"normal"=Inf,"gamma"=Inf)
if ( "lognormal" %in% empiricalHyperPriorForm ) {
lambda.hyperprior.fits["lognormal"] <- suppressWarnings(ks.test(samples.lambda,'plnorm',meanlog=log((m.lambda^2)/sqrt(v.lambda+m.lambda^2)),sdlog=sqrt( log(1+v.lambda/(m.lambda^2)))))$statistic
mu.hyperprior.fits["lognormal"] <- suppressWarnings(ks.test(samples.mu,'plnorm',meanlog=log((m.mu^2)/sqrt(v.mu+m.mu^2)),sdlog=sqrt( log(1+v.mu/(m.mu^2)))))$statistic
phi.hyperprior.fits["lognormal"] <- suppressWarnings(ks.test(samples.phi,'plnorm',meanlog=log((m.phi^2)/sqrt(v.phi+m.phi^2)),sdlog=sqrt( log(1+v.phi/(m.phi^2)))))$statistic
}
if ( "normal" %in% empiricalHyperPriorForm ) {
lambda.hyperprior.fits["normal"] <- suppressWarnings(ks.test(samples.lambda,'pnorm',mean=m.lambda,sd=sqrt(v.lambda)))$statistic
mu.hyperprior.fits["normal"] <- suppressWarnings(ks.test(samples.mu,'pnorm',mean=m.mu,sd=sqrt(v.mu)))$statistic
phi.hyperprior.fits["normal"] <- suppressWarnings(ks.test(samples.phi,'pnorm',mean=m.phi,sd=sqrt(v.phi)))$statistic
}
if ( "gamma" %in% empiricalHyperPriorForm ) {
lambda.hyperprior.fits["gamma"] <- suppressWarnings(ks.test(samples.lambda,'pgamma',shape=(m.lambda^2)/v.lambda,rate=m.lambda/v.lambda))$statistic
mu.hyperprior.fits["gamma"] <- suppressWarnings(ks.test(samples.mu,'pgamma',shape=(m.mu^2)/v.mu,rate=m.mu/v.mu))$statistic
phi.hyperprior.fits["gamma"] <- suppressWarnings(ks.test(samples.phi,'pgamma',shape=(m.phi^2)/v.phi,rate=m.phi/v.phi))$statistic
}
lambda.hyper.prior.form <- names(which.min(lambda.hyperprior.fits[empiricalHyperPriorForm]))
mu.hyper.prior.form <- names(which.min(mu.hyperprior.fits[empiricalHyperPriorForm]))
phi.hyper.prior.form <- names(which.min(phi.hyperprior.fits[empiricalHyperPriorForm]))
if ( lambda.hyper.prior.form == "lognormal" ) {
speciationRatePriorMean <- mean( log(samples.lambda) )
speciationRatePriorStDev <- sqrt( empiricalHyperPriorInflation * mean( (log(samples.lambda)-speciationRatePriorMean)^2 ) )
} else if ( lambda.hyper.prior.form == "normal" ) {
speciationRatePriorMean <- m.lambda
speciationRatePriorStDev <- sqrt(v.lambda)
} else if ( lambda.hyper.prior.form == "gamma" ) {
speciationRatePriorMean <- (m.lambda^2)/v.lambda
speciationRatePriorStDev <- m.lambda/v.lambda
}
if ( mu.hyper.prior.form == "lognormal" ) {
extinctionRatePriorMean <- mean( log(samples.mu) )
extinctionRatePriorStDev <- sqrt( empiricalHyperPriorInflation * mean( (log(samples.mu)-extinctionRatePriorMean)^2 ) )
} else if ( mu.hyper.prior.form == "normal" ) {
extinctionRatePriorMean <- m.mu
extinctionRatePriorStDev <- sqrt(v.mu)
} else if ( mu.hyper.prior.form == "gamma" ) {
extinctionRatePriorMean <- (m.mu^2)/v.mu
extinctionRatePriorStDev <- m.mu/v.mu
}
if ( phi.hyper.prior.form == "lognormal" ) {
fossilizationRatePriorMean <- mean( log(samples.phi) )
fossilizationRatePriorStDev <- sqrt( empiricalHyperPriorInflation * mean( (log(samples.phi)-fossilizationRatePriorMean)^2 ) )
} else if ( phi.hyper.prior.form == "normal" ) {
fossilizationRatePriorMean <- m.phi
fossilizationRatePriorStDev <- sqrt(v.phi)
} else if ( phi.hyper.prior.form == "gamma" ) {
fossilizationRatePriorMean <- (m.phi^2)/v.phi
fossilizationRatePriorStDev <- m.phi/v.phi
}
initialSpeciationRate <- m.lambda
initialExtinctionRate <- m.mu
initialFossilizationRate <- m.phi
empirical.hyperprior.samples <- data.frame(1:nrow(samples),samples[,"likelihood"],samples[,"prior"],samples.lambda,samples.mu,samples.phi)
write.table(empirical.hyperprior.samples[-1,],file='samplesHyperprior.txt',row.names=FALSE,col.names=c('Iteration','Likelihood','Prior','speciation','extinction','fossilization'),sep='\t')
summary.file <- "empiricalHyperpriorSummary.txt"
cat("Summary of the empirical hyperprior analysis",file=summary.file)
cat("\n============================================\n\n",file=summary.file,append=TRUE)
cat("Specation rate hyperprior\n",file=summary.file,append=TRUE)
cat("-------------------------\n",file=summary.file,append=TRUE)
cat("Form:\t",lambda.hyper.prior.form,'\n',file=summary.file,append=TRUE)
if ( lambda.hyper.prior.form == "lognormal" ) {
cat("meanlog:\t",speciationRatePriorMean,"\n",file=summary.file,append=TRUE)
cat("sdlog:\t\t",speciationRatePriorStDev,"\n",file=summary.file,append=TRUE)
} else if ( lambda.hyper.prior.form == "normal" ) {
cat("mean:\t",speciationRatePriorMean,"\n",file=summary.file,append=TRUE)
cat("sd:\t\t",speciationRatePriorStDev,"\n",file=summary.file,append=TRUE)
} else if ( lambda.hyper.prior.form == "gamma" ) {
cat("shape:\t",speciationRatePriorMean,"\n",file=summary.file,append=TRUE)
cat("rate:\t\t",speciationRatePriorStDev,"\n",file=summary.file,append=TRUE)
}
cat("\nExtinction rate hyperprior\n",file=summary.file,append=TRUE)
cat("-------------------------\n",file=summary.file,append=TRUE)
cat("Form:\t",mu.hyper.prior.form,'\n',file=summary.file,append=TRUE)
if ( mu.hyper.prior.form == "lognormal" ) {
cat("meanlog:\t",extinctionRatePriorMean,"\n",file=summary.file,append=TRUE)
cat("sdlog:\t\t",extinctionRatePriorStDev,"\n",file=summary.file,append=TRUE)
} else if ( mu.hyper.prior.form == "normal" ) {
cat("mean:\t",extinctionRatePriorMean,"\n",file=summary.file,append=TRUE)
cat("sd:\t\t",extinctionRatePriorStDev,"\n",file=summary.file,append=TRUE)
} else if ( mu.hyper.prior.form == "gamma" ) {
cat("shape:\t",extinctionRatePriorMean,"\n",file=summary.file,append=TRUE)
cat("rate:\t\t",extinctionRatePriorStDev,"\n",file=summary.file,append=TRUE)
}
cat("\nFossilization rate hyperprior\n",file=summary.file,append=TRUE)
cat("-------------------------\n",file=summary.file,append=TRUE)
cat("Form:\t",phi.hyper.prior.form,'\n',file=summary.file,append=TRUE)
if ( phi.hyper.prior.form == "lognormal" ) {
cat("meanlog:\t",fossilizationRatePriorMean,"\n",file=summary.file,append=TRUE)
cat("sdlog:\t\t",fossilizationRatePriorStDev,"\n",file=summary.file,append=TRUE)
} else if ( phi.hyper.prior.form == "normal" ) {
cat("mean:\t",fossilizationRatePriorMean,"\n",file=summary.file,append=TRUE)
cat("sd:\t\t",fossilizationRatePriorStDev,"\n",file=summary.file,append=TRUE)
} else if ( phi.hyper.prior.form == "gamma" ) {
cat("shape:\t",fossilizationRatePriorMean,"\n",file=summary.file,append=TRUE)
cat("rate:\t\t",fossilizationRatePriorStDev,"\n",file=summary.file,append=TRUE)
}
} else { # no empirical hyper priors
lambda.hyper.prior.form <- mu.hyper.prior.form <- phi.hyper.prior.form <- empiricalHyperPriorForm[1]
}
likelihood <- function(lambda,lambdaTimes,mu,muTimes,phi,phiTimes,tMassExtinction,pSurvival) {
if ( priorOnly == FALSE ) {
if( any( lambda < 0 ) | any( mu < 0 ) ) {
lnl <- -Inf
} else {
if ( use_tree_sample == TRUE ) {
ln_probs <- c()
max_prob <- -Inf
for (i in 1:NUM_SAMPLED_TREES) {
if ( USE_EBDSTP == FALSE ) {
ln_probs[i] <- tess.likelihood.efbd(node = times[[i]],
lambda = lambda,
mu = mu,
phi = phi,
rateChangeTimesLambda = lambdaTimes,
rateChangeTimesMu = muTimes,
rateChangeTimesPhi = phiTimes,
massExtinctionTimes= tMassExtinction,
massExtinctionSurvivalProbabilities = pSurvival,
samplingStrategy = samplingStrategy,
samplingProbability = samplingProbability,
MRCA = MRCA,
CONDITION = CONDITION,
log = TRUE)
} else {
ln_probs[i] <- tess.likelihood.ebdstp(node = times[[i]],
lambda = lambda,
mu = mu,
phi = phi,
r = 0,
rateChangeTimesLambda = lambdaTimes,
rateChangeTimesMu = muTimes,
rateChangeTimesPhi = phiTimes,
massDeathTimes= tMassExtinction,
massDeathProbabilities = 1.0 - pSurvival,
samplingStrategy = samplingStrategy,
samplingProbability = samplingProbability,
MRCA = MRCA,
CONDITION = CONDITION,
log = TRUE)
}
if ( max_prob < ln_probs[i] ) max_prob <- ln_probs[i]
}
sum <- 0.0
for (i in 1:NUM_SAMPLED_TREES) {
sum <- sum + exp( ln_probs[i] - max_prob ) / NUM_SAMPLED_TREES
}
lnl <- log( sum ) + max_prob
} else {
if ( USE_EBDSTP == FALSE ) {
lnl <- tess.likelihood.efbd(node = times,
lambda = lambda,
mu = mu,
phi = phi,
rateChangeTimesLambda = lambdaTimes,
rateChangeTimesMu = muTimes,
rateChangeTimesPhi = phiTimes,
massExtinctionTimes= tMassExtinction,
massExtinctionSurvivalProbabilities = pSurvival,
samplingStrategy = samplingStrategy,
samplingProbability = samplingProbability,
MRCA = MRCA,
CONDITION = CONDITION,
log = TRUE)
} else {
lnl <- tess.likelihood.ebdstp(node = times,
lambda = lambda,
mu = mu,
phi = phi,
r = 0,
rateChangeTimesLambda = lambdaTimes,
rateChangeTimesMu = muTimes,
rateChangeTimesPhi = phiTimes,
massDeathTimes= tMassExtinction,
massDeathProbabilities = 1.0 - pSurvival,
samplingStrategy = samplingStrategy,
samplingProbability = samplingProbability,
MRCA = MRCA,
CONDITION = CONDITION,
log = TRUE)
}
}
}
} else {
lnl <- 0
}
return (lnl)
}
prior <- function(timesLambda,valuesLambda,timesMu,valuesMu,timesPhi,valuesPhi,timesMassExtinction,valuesMassExtinction) {
if( any( valuesLambda < 0 ) | any( valuesMu < 0 ) | any( valuesPhi < 0 ) ) {
return (-Inf)
}
#####################
## Speciation Rate ##
#####################
# prior on number of changes for lambda
k <- length(timesLambda)
lnp <- dpois(k,numExpectedRateChanges,log=TRUE)
if ( k > 0 ) {
# prior on change times
lnp <- lnp + sum( dunif(timesLambda,0,AGE,log=TRUE) )
}
# prior on change values
if ( lambda.hyper.prior.form == "lognormal" ) {
lnp <- lnp + sum( dlnorm(valuesLambda,meanlog=speciationRatePriorMean,speciationRatePriorStDev,log=TRUE) )
} else if ( lambda.hyper.prior.form == "normal" ) {
lnp <- lnp + sum( dnorm(valuesLambda,speciationRatePriorMean,speciationRatePriorStDev,log=TRUE) )
} else if ( lambda.hyper.prior.form == "gamma" ) {
lnp <- lnp + sum( dgamma(valuesLambda,speciationRatePriorMean,speciationRatePriorStDev,log=TRUE) )
}
#####################
## Extinction Rate ##
#####################
# prior on number of changes for mu
k <- length(timesMu)
lnp <- lnp + dpois(k,numExpectedRateChanges,log=TRUE)
if ( k > 0 ) {
# prior on change times
lnp <- lnp + sum( dunif(timesMu,0,AGE,log=TRUE) )
}
# prior on change values
if ( mu.hyper.prior.form == "lognormal" ) {
lnp <- lnp + sum( dlnorm(valuesMu,meanlog=extinctionRatePriorMean,extinctionRatePriorStDev,log=TRUE) )
} else if ( mu.hyper.prior.form == "normal" ) {
lnp <- lnp + sum( dnorm(valuesMu,extinctionRatePriorMean,extinctionRatePriorStDev,log=TRUE) )
} else if ( mu.hyper.prior.form == "gamma" ) {
lnp <- lnp + sum( dgamma(valuesMu,extinctionRatePriorMean,extinctionRatePriorStDev,log=TRUE) )
}
########################
## Fossilization Rate ##
########################
# prior on number of changes for phi
k <- length(timesPhi)
lnp <- lnp + dpois(k,numExpectedRateChanges,log=TRUE)
if ( k > 0 ) {
# prior on change times
lnp <- lnp + sum( dunif(timesPhi,0,AGE,log=TRUE) )
}
# prior on change values
if ( phi.hyper.prior.form == "lognormal" ) {
lnp <- lnp + sum( dlnorm(valuesPhi,fossilizationRatePriorMean,fossilizationRatePriorStDev,log=TRUE) )
} else if ( phi.hyper.prior.form == "normal" ) {
lnp <- lnp + sum( dnorm(valuesPhi,fossilizationRatePriorMean,fossilizationRatePriorStDev,log=TRUE) )
} else if ( phi.hyper.prior.form == "gamma" ) {
lnp <- lnp + sum( dgamma(valuesPhi,fossilizationRatePriorMean,fossilizationRatePriorStDev,log=TRUE) )
}
#####################
## Mass-Extinction ##
#####################
# prior on number of changes for mass-extinction events
k <- length(timesMassExtinction)
lnp <- lnp + dpois(k,numExpectedMassExtinctions,log=TRUE)
if ( k > 0 ) {
# prior on change times
lnp <- lnp + sum( dunif(timesMassExtinction,0,AGE,log=TRUE) )
# prior on change values
lnp <- lnp + sum( dbeta(valuesMassExtinction,pMassExtinctionPriorShape1,pMassExtinctionPriorShape2,log=TRUE) )
}
return (lnp)
}
if ( use_tree_sample == TRUE ) {
AGE <- max( as.numeric( tess.branching.times( tree[[1]] )$age ) )
if ( length(tree) > 1 ) {
for ( i in 2:length(tree)) AGE <- min( c(AGE, max(as.numeric( tess.branching.times( tree[[i]] )$age) ) ) )
}
} else {
AGE <- max( as.numeric( tess.branching.times( tree )$age ) )
}
# these are the parameters we sample
posterior <- c()
kLambda <- c()
kMu <- c()
kPhi <- c()
kMassExtinction <- c()
speciationRateValues <- list()
speciationRateChangeTimes <- list()
extinctionRateValues <- list()
extinctionRateChangeTimes <- list()
fossilizationRateValues <- list()
fossilizationRateChangeTimes <- list()
survivalProbability <- list()
massExtinctionTime <- list()
cat("SpeciationRates\n",file="SpeciationRates.txt")
cat("SpeciationRateChanges\n",file="SpeciationRateChanges.txt")
cat("ExtinctionRates\n",file="ExtinctionRates.txt")
cat("ExtinctionRateChangeTimes\n",file="ExtinctionRateChanges.txt")
cat("FossilizationRates\n",file="FossilizationRates.txt")
cat("FossilizationRateChangeTimes\n",file="FossilizationRateChanges.txt")
cat("SurvivalProbabilities\n",file="SurvivalProbabilities.txt")
cat("MassExtinctionTimes\n",file="MassExtinctionTimes.txt")
cat(paste("Iteration","posterior","NumSpeciation","numExtinction","numMassExtinctions",sep="\t"),"\n",file="samples_numCategories.txt")
# the "random" starting values
lambda <- initialSpeciationRate
mu <- initialExtinctionRate
phi <- initialFossilizationRate
lambdaChangeTimes <- initialSpeciationRateChangeTime
muChangeTimes <- initialExtinctionRateChangeTime
phiChangeTimes <- initialFossilizationRateChangeTime
pMassExtinction <- pInitialMassExtinction
tMassExtinction <- tInitialMassExtinction
initialPP <- prior(lambdaChangeTimes,lambda,muChangeTimes,mu,phiChangeTimes,phi,tMassExtinction,pMassExtinction) + likelihood(lambda,lambdaChangeTimes,mu,muChangeTimes,phi,phiChangeTimes,tMassExtinction,pMassExtinction)
currentPP <- initialPP
## check if the prior and posterior gives valid probabilities
MAX_TRIES <- 1000
for ( i in 1:MAX_TRIES ) {
tmpPP <- prior(lambdaChangeTimes,lambda,muChangeTimes,mu,phiChangeTimes,phi,tMassExtinction,pMassExtinction) + likelihood(lambda,lambdaChangeTimes,mu,muChangeTimes,phi,phiChangeTimes,tMassExtinction,pMassExtinction)
if ( is.finite(tmpPP) ) break
lambdaChangeTimes <- c()
lambda <- rlnorm(1,speciationRatePriorMean,speciationRatePriorStDev)
muChangeTimes <- c()
mu <- rlnorm(1,extinctionRatePriorMean,extinctionRatePriorStDev)
phiChangeTimes <- c()
phi <- rlnorm(1,fossilizationRatePriorMean,fossilizationRatePriorStDev)
tMassExtinction <- c()
pMassExtinction <- c()
}
## MCMC proposal sizes
deltaLambdaTime <- 0.1
lambdaTimeAccepted <- 0
lambdaTimeTried <- 0
deltaLambdaValue <- 0.1
lambdaValueAccepted <- 0
lambdaValueTried <- 0
deltaMuTime <- 0.1
muTimeAccepted <- 0
muTimeTried <- 0
deltaMuValue <- 0.1
muValueAccepted <- 0
muValueTried <- 0
deltaPhiTime <- 0.1
phiTimeAccepted <- 0
phiTimeTried <- 0
deltaPhiValue <- 0.1
phiValueAccepted <- 0
phiValueTried <- 0
deltaMassExtinctionTime <- 0.1
massExtinctionTimeAccepted <- 0
massExtinctionTimeTried <- 0
deltaMassExtinctionValue <- 0.1
massExtinctionValueAccepted <- 0
massExtinctionValueTried <- 0
finished <- FALSE
SAMPLE <- FALSE
min.ess <- 0
startTime <- Sys.time()
i <- 1
if ( verbose ) {
cat("\nPerforming CoMET analysis.\n\n")
if( burn > 0 ) {
cat("Burning-in the chain ...\n")
} else {
cat("Running the chain ... \n")
}
cat("0--------25--------50--------75--------100\n")
bar <- txtProgressBar(style=1,width=42)
}
while ( !finished ) {
# change the time of an event
if ( length(lambdaChangeTimes) > 0 ) {
# compute index of value that will be changed
idx <- floor( runif(1,0,length(lambdaChangeTimes)) ) + 1
# store current value
t <- lambdaChangeTimes[idx]
# compute current posterior
oldPP <- currentPP
# propose new value
tPrime <- rnorm(1,t,deltaLambdaTime)
lambdaChangeTimes[idx] <- tPrime
# compute new posterior
if ( tPrime > 0 && tPrime < AGE ) {
newPP <- prior(lambdaChangeTimes,lambda,muChangeTimes,mu,phiChangeTimes,phi,tMassExtinction,pMassExtinction) + likelihood(lambda,lambdaChangeTimes,mu,muChangeTimes,phi,phiChangeTimes,tMassExtinction,pMassExtinction)
} else {
newPP <- -Inf
}
# accept/reject
if ( is.finite(newPP - oldPP) == FALSE || log( runif(1,0,1) ) > ( newPP - oldPP ) ) {
# reject
lambdaChangeTimes[idx] <- t
} else {
currentPP <- newPP
lambdaTimeAccepted <- lambdaTimeAccepted + 1
}
lambdaTimeTried <- lambdaTimeTried + 1
}
# change the value of an event
# compute index of value that will be changed
idx <- floor( runif(1,0,length(lambda)) ) + 1
# store current value
v <- lambda[idx]
# compute current posterior
oldPP <- currentPP
# propose new value
u <- rnorm(1,0,deltaLambdaValue)
scalingFactor <- exp( u )
vPrime <- v * scalingFactor
lambda[idx] <- vPrime
# compute the Hastings ratio
lnHastingsratio <- log( scalingFactor )
# compute new posterior
newPP <- prior(lambdaChangeTimes,lambda,muChangeTimes,mu,phiChangeTimes,phi,tMassExtinction,pMassExtinction) + likelihood(lambda,lambdaChangeTimes,mu,muChangeTimes,phi,phiChangeTimes,tMassExtinction,pMassExtinction)
# accept/reject
if ( is.finite(newPP - oldPP + lnHastingsratio) == FALSE || log( runif(1,0,1) ) > ( newPP - oldPP + lnHastingsratio ) ) {
# reject
lambda[idx] <- v
} else {
currentPP <- newPP
lambdaValueAccepted <- lambdaValueAccepted + 1
}
lambdaValueTried <- lambdaValueTried + 1
if ( estimateNumberSpeciationRateChanges == TRUE ) {
# We need to randomly pick a birth or death move
# Otherwise we might give birth and die every time
u <- runif(1,0,1)
if ( u > 0.5 ) {
# birth move
# compute current posterior
oldPP <- currentPP
# randomly pick a new time
t <- runif(1,0,AGE)
# randomly pick a new value
if ( lambda.hyper.prior.form == "lognormal" ) {
v <- rlnorm(1,speciationRatePriorMean,speciationRatePriorStDev)
} else if ( lambda.hyper.prior.form == "normal" ) {
v <- rnorm(1,speciationRatePriorMean,speciationRatePriorStDev)
} else if ( lambda.hyper.prior.form == "gamma" ) {
v <- rgamma(1,speciationRatePriorMean,speciationRatePriorStDev)
}
# construct the new parameters
lambdaChangeTimes[length(lambdaChangeTimes)+1] <- t
lambda[length(lambda)+1] <- v
# compute new posterior
newPP <- prior(lambdaChangeTimes,lambda,muChangeTimes,mu,phiChangeTimes,phi,tMassExtinction,pMassExtinction) + likelihood(lambda,lambdaChangeTimes,mu,muChangeTimes,phi,phiChangeTimes,tMassExtinction,pMassExtinction)
# compute proposal ratio
kPrime <- length(lambdaChangeTimes)
if ( lambda.hyper.prior.form == "lognormal" ) {
pr <- 1 / ( dunif(t,0,AGE) * dlnorm(v,speciationRatePriorMean,speciationRatePriorStDev) )
} else if ( lambda.hyper.prior.form == "normal" ) {
pr <- 1 / ( dunif(t,0,AGE) * dnorm(v,speciationRatePriorMean,speciationRatePriorStDev) )
} else if ( lambda.hyper.prior.form == "gamma" ) {
pr <- 1 / ( dunif(t,0,AGE) * dgamma(v,speciationRatePriorMean,speciationRatePriorStDev) )
}
# accept/reject
if ( is.finite(newPP - oldPP + log(pr)) == FALSE || log( runif(1,0,1) ) > ( newPP - oldPP + log( pr ) ) ) {
# reject
lambdaChangeTimes <- lambdaChangeTimes[-length(lambdaChangeTimes)]
lambda <- lambda[-length(lambda)]
} else {
currentPP <- newPP
}
} else {
# death move
if ( length(lambdaChangeTimes) > 0 ) {
# compute current posterior
oldPP <- currentPP
# randomly pick an index
idx <- floor( runif(1,0,length(lambdaChangeTimes)) ) + 1
# store the current value
t <- lambdaChangeTimes
v <- lambda
# construct the new parameters
lambdaChangeTimes <- lambdaChangeTimes[-idx]
lambda <- lambda[-(idx+1)]
# compute new posterior
newPP <- prior(lambdaChangeTimes,lambda,muChangeTimes,mu,phiChangeTimes,phi,tMassExtinction,pMassExtinction) + likelihood(lambda,lambdaChangeTimes,mu,muChangeTimes,phi,phiChangeTimes,tMassExtinction,pMassExtinction)
# compute proposal ratio
kPrime <- length(lambdaChangeTimes) + 1
if ( lambda.hyper.prior.form == "lognormal" ) {
pr2 <- dunif(t[idx],0,AGE) * dlnorm(v[idx+1],speciationRatePriorMean,speciationRatePriorStDev)
} else if ( lambda.hyper.prior.form == "normal" ) {
pr2 <- dunif(t[idx],0,AGE) * dnorm(v[idx+1],speciationRatePriorMean,speciationRatePriorStDev)
} else if ( lambda.hyper.prior.form == "gamma" ) {
pr2 <- dunif(t[idx],0,AGE) * dgamma(v[idx+1],speciationRatePriorMean,speciationRatePriorStDev)
}
# accept/reject
if ( is.finite(newPP - oldPP + log(pr2)) == FALSE || log( runif(1,0,1) ) > ( newPP - oldPP + log( pr2 ) ) ) {
# reject
lambdaChangeTimes <- t
lambda <- v
} else {
currentPP <- newPP
}
}
}
}
#####################
## Extinction Rate ##
#####################
# change the time of an event
if ( length(muChangeTimes) > 0 ) {
# compute index of value that will be changed
idx <- floor( runif(1,0,length(muChangeTimes)) ) + 1
# store current value
t <- muChangeTimes[idx]
# compute current posterior
oldPP <- currentPP
# propose new value
tPrime <- rnorm(1, t, deltaMuTime)
muChangeTimes[idx] <- tPrime
# compute new posterior
if ( tPrime > 0 && tPrime < AGE ) {
newPP <- prior(lambdaChangeTimes,lambda,muChangeTimes,mu,phiChangeTimes,phi,tMassExtinction,pMassExtinction) + likelihood(lambda,lambdaChangeTimes,mu,muChangeTimes,phi,phiChangeTimes,tMassExtinction,pMassExtinction)
} else {
newPP <- -Inf
}
# accept/reject
if ( is.finite(newPP - oldPP) == FALSE || log( runif(1,0,1) ) > ( newPP - oldPP ) ) {
# reject
muChangeTimes[idx] <- t
} else {
currentPP <- newPP
muTimeAccepted <- muTimeAccepted + 1
}
muTimeTried <- muTimeTried + 1
}
# change the value of an event
# compute index of value that will be changed
idx <- floor( runif(1,0,length(mu)) ) + 1
# store current value
v <- mu[idx]
# compute current posterior
oldPP <- currentPP
# propose new value
u <- rnorm(1,0,deltaMuValue)
scalingFactor <- exp( u )
vPrime <- v * scalingFactor
mu[idx] <- vPrime
# compute the Hastings ratio
lnHastingsratio <- log( scalingFactor )
# compute new posterior
newPP <- prior(lambdaChangeTimes,lambda,muChangeTimes,mu,phiChangeTimes,phi,tMassExtinction,pMassExtinction) + likelihood(lambda,lambdaChangeTimes,mu,muChangeTimes,phi,phiChangeTimes,tMassExtinction,pMassExtinction)
# accept/reject
if ( is.finite(newPP - oldPP + lnHastingsratio) == FALSE || log( runif(1,0,1) ) > ( newPP - oldPP + lnHastingsratio ) ) {
# reject
mu[idx] <- v
} else {
currentPP <- newPP
muValueAccepted <- muValueAccepted + 1
}
muValueTried <- muValueTried + 1
if ( estimateNumberExtinctionRateChanges == TRUE ) {
# We need to randomly pick a birth or death move
# Otherwise we might give birth and die every time
u <- runif(1,0,1)
if ( u > 0.5 ) {
# birth move
# compute current posterior
oldPP <- currentPP
# randomly pick a new time
t <- runif(1,0,AGE)
# randomly pick a new value
if ( mu.hyper.prior.form == "lognormal" ) {
v <- rlnorm(1,extinctionRatePriorMean,extinctionRatePriorStDev)
} else if ( mu.hyper.prior.form == "normal" ) {
v <- rnorm(1,extinctionRatePriorMean,extinctionRatePriorStDev)
} else if ( mu.hyper.prior.form == "gamma" ) {
v <- rgamma(1,extinctionRatePriorMean,extinctionRatePriorStDev)
}
# construct the new parameters
muChangeTimes[length(muChangeTimes)+1] <- t
mu[length(mu)+1] <- v
# compute new posterior
newPP <- prior(lambdaChangeTimes,lambda,muChangeTimes,mu,phiChangeTimes,phi,tMassExtinction,pMassExtinction) + likelihood(lambda,lambdaChangeTimes,mu,muChangeTimes,phi,phiChangeTimes,tMassExtinction,pMassExtinction)
# compute proposal ratio
kPrime <- length(muChangeTimes)
if ( mu.hyper.prior.form == "lognormal" ) {
pr <- 1 / ( dunif(t,0,AGE) * dlnorm(v,extinctionRatePriorMean,extinctionRatePriorStDev) )
} else if ( mu.hyper.prior.form == "normal" ) {
pr <- 1 / ( dunif(t,0,AGE) * dnorm(v,extinctionRatePriorMean,extinctionRatePriorStDev) )
} else if ( mu.hyper.prior.form == "gamma" ) {
pr <- 1 / ( dunif(t,0,AGE) * dgamma(v,extinctionRatePriorMean,extinctionRatePriorStDev) )
}
# accept/reject
if ( is.finite(newPP - oldPP + log(pr)) == FALSE || log( runif(1,0,1) ) > ( newPP - oldPP + log( pr ) ) ) {
# reject
muChangeTimes <- muChangeTimes[-length(muChangeTimes)]
mu <- mu[-length(mu)]
} else {
currentPP <- newPP
}
} else {
# death move
if ( length(muChangeTimes) > 0 ) {
# compute current posterior
oldPP <- currentPP
# randomly pick an index
idx <- floor( runif(1,0,length(muChangeTimes)) ) + 1
# store the current value
t <- muChangeTimes
v <- mu
# construct the new parameters
muChangeTimes <- muChangeTimes[-idx]
mu <- mu[-(idx+1)]
# compute new posterior
newPP <- prior(lambdaChangeTimes,lambda,muChangeTimes,mu,phiChangeTimes,phi,tMassExtinction,pMassExtinction) + likelihood(lambda,lambdaChangeTimes,mu,muChangeTimes,phi,phiChangeTimes,tMassExtinction,pMassExtinction)
# compute proposal ratio
kPrime <- length(muChangeTimes) + 1
if ( mu.hyper.prior.form == "lognormal" ) {
pr2 <- dunif(t[idx],0,AGE) * dlnorm(v[idx+1],extinctionRatePriorMean,extinctionRatePriorStDev)
} else if ( mu.hyper.prior.form == "normal" ) {
pr2 <- dunif(t[idx],0,AGE) * dnorm(v[idx+1],extinctionRatePriorMean,extinctionRatePriorStDev)
} else if ( mu.hyper.prior.form == "gamma" ) {
pr2 <- dunif(t[idx],0,AGE) * dgamma(v[idx+1],extinctionRatePriorMean,extinctionRatePriorStDev)
}
# accept/reject
if ( is.finite(newPP - oldPP + log(pr2)) == FALSE || log( runif(1,0,1) ) > ( newPP - oldPP + log( pr2 ) ) ) {
# reject
muChangeTimes <- t
mu <- v
} else {
currentPP <- newPP
}
}
}
}
########################
## Fossilization Rate ##
########################
# change the time of an event
if ( length(phiChangeTimes) > 0 ) {
# compute index of value that will be changed
idx <- floor( runif(1,0,length(phiChangeTimes)) ) + 1
# store current value
t <- phiChangeTimes[idx]
# compute current posterior
oldPP <- currentPP
# propose new value
tPrime <- rnorm(1, t, deltaPhiTime)
phiChangeTimes[idx] <- tPrime
# compute new posterior
if ( tPrime > 0 && tPrime < AGE ) {
newPP <- prior(lambdaChangeTimes,lambda,muChangeTimes,mu,phiChangeTimes,phi,tMassExtinction,pMassExtinction) + likelihood(lambda,lambdaChangeTimes,mu,muChangeTimes,phi,phiChangeTimes,tMassExtinction,pMassExtinction)
} else {
newPP <- -Inf
}
# accept/reject
if ( is.finite(newPP - oldPP) == FALSE || log( runif(1,0,1) ) > ( newPP - oldPP ) ) {
# reject
phiChangeTimes[idx] <- t
} else {
currentPP <- newPP
phiTimeAccepted <- phiTimeAccepted + 1
}
phiTimeTried <- phiTimeTried + 1
}
# change the value of an event
# compute index of value that will be changed
idx <- floor( runif(1,0,length(phi)) ) + 1
# store current value
v <- phi[idx]
# compute current posterior
oldPP <- currentPP
# propose new value
u <- rnorm(1,0,deltaPhiValue)
scalingFactor <- exp( u )
vPrime <- v * scalingFactor
phi[idx] <- vPrime
# compute the Hastings ratio
lnHastingsratio <- log( scalingFactor )
# compute new posterior
newPP <- prior(lambdaChangeTimes,lambda,muChangeTimes,mu,phiChangeTimes,phi,tMassExtinction,pMassExtinction) + likelihood(lambda,lambdaChangeTimes,mu,muChangeTimes,phi,phiChangeTimes,tMassExtinction,pMassExtinction)
# accept/reject
if ( is.finite(newPP - oldPP + lnHastingsratio) == FALSE || log( runif(1,0,1) ) > ( newPP - oldPP + lnHastingsratio ) ) {
# reject
phi[idx] <- v
} else {
currentPP <- newPP
phiValueAccepted <- phiValueAccepted + 1
}
phiValueTried <- phiValueTried + 1
if ( estimateNumberFossilizationRateChanges == TRUE ) {
# We need to randomly pick a birth or death move
# Otherwise we might give birth and die every time
u <- runif(1,0,1)
if ( u > 0.5 ) {
# birth move
# compute current posterior
oldPP <- currentPP
# randomly pick a new time
t <- runif(1,0,AGE)
# randomly pick a new value
if ( phi.hyper.prior.form == "lognormal" ) {
v <- rlnorm(1,fossilizationRatePriorMean,fossilizationRatePriorStDev)
} else if ( phi.hyper.prior.form == "normal" ) {
v <- rnorm(1,fossilizationRatePriorMean,fossilizationRatePriorStDev)
} else if ( phi.hyper.prior.form == "gamma" ) {
v <- rgamma(1,fossilizationRatePriorMean,fossilizationRatePriorStDev)
}
# construct the new parameters
phiChangeTimes[length(phiChangeTimes)+1] <- t
phi[length(phi)+1] <- v
# compute new posterior
newPP <- prior(lambdaChangeTimes,lambda,muChangeTimes,mu,phiChangeTimes,phi,tMassExtinction,pMassExtinction) + likelihood(lambda,lambdaChangeTimes,mu,muChangeTimes,phi,phiChangeTimes,tMassExtinction,pMassExtinction)
# compute proposal ratio
kPrime <- length(phiChangeTimes)
if ( phi.hyper.prior.form == "lognormal" ) {
pr <- 1 / ( dunif(t,0,AGE) * dlnorm(v,fossilizationRatePriorMean,fossilizationRatePriorStDev) )
} else if ( phi.hyper.prior.form == "normal" ) {
pr <- 1 / ( dunif(t,0,AGE) * dnorm(v,fossilizationRatePriorMean,fossilizationRatePriorStDev) )
} else if ( phi.hyper.prior.form == "gamma" ) {
pr <- 1 / ( dunif(t,0,AGE) * dgamma(v,fossilizationRatePriorMean,fossilizationRatePriorStDev) )
}
# accept/reject
if ( is.finite(newPP - oldPP + log(pr)) == FALSE || log( runif(1,0,1) ) > ( newPP - oldPP + log( pr ) ) ) {
# reject
phiChangeTimes <- phiChangeTimes[-length(phiChangeTimes)]
phi <- phi[-length(phi)]
} else {
currentPP <- newPP
}
} else {
# death move
if ( length(phiChangeTimes) > 0 ) {
# compute current posterior
oldPP <- currentPP
# randomly pick an index
idx <- floor( runif(1,0,length(phiChangeTimes)) ) + 1
# store the current value
t <- phiChangeTimes
v <- phi
# construct the new parameters
phiChangeTimes <- phiChangeTimes[-idx]
phi <- phi[-(idx+1)]
# compute new posterior
newPP <- prior(lambdaChangeTimes,lambda,muChangeTimes,mu,phiChangeTimes,phi,tMassExtinction,pMassExtinction) + likelihood(lambda,lambdaChangeTimes,mu,muChangeTimes,phi,phiChangeTimes,tMassExtinction,pMassExtinction)
# compute proposal ratio
kPrime <- length(phiChangeTimes) + 1
if ( phi.hyper.prior.form == "lognormal" ) {
pr2 <- dunif(t[idx],0,AGE) * dlnorm(v[idx+1],fossilizationRatePriorMean,fossilizationRatePriorStDev)
} else if ( phi.hyper.prior.form == "normal" ) {
pr2 <- dunif(t[idx],0,AGE) * dnorm(v[idx+1],fossilizationRatePriorMean,fossilizationRatePriorStDev)
} else if ( phi.hyper.prior.form == "gamma" ) {
pr2 <- dunif(t[idx],0,AGE) * dgamma(v[idx+1],fossilizationRatePriorMean,fossilizationRatePriorStDev)
}
# accept/reject
if ( is.finite(newPP - oldPP + log(pr2)) == FALSE || log( runif(1,0,1) ) > ( newPP - oldPP + log( pr2 ) ) ) {
# reject
phiChangeTimes <- t
phi <- v
} else {
currentPP <- newPP
}
}
}
}
#####################
## Mass-Extinction ##
#####################
# change the time of an event
if ( estimateMassExtinctionTimes == TRUE && length(tMassExtinction) > 0 ) {
# compute index of value that will be changed
idx <- floor( runif(1,0,length(tMassExtinction)) ) + 1
# store current value
t <- tMassExtinction[idx]
# compute current posterior
oldPP <- currentPP
# propose new value
tPrime <- rnorm(1, t, deltaMassExtinctionTime)
tMassExtinction[idx] <- tPrime
# compute new posterior
if ( tPrime > 0 && tPrime < AGE ) {
newPP <- prior(lambdaChangeTimes,lambda,muChangeTimes,mu,phiChangeTimes,phi,tMassExtinction,pMassExtinction) + likelihood(lambda,lambdaChangeTimes,mu,muChangeTimes,phi,phiChangeTimes,tMassExtinction,pMassExtinction)
} else {
newPP <- -Inf
}
# accept/reject
if ( is.finite(newPP - oldPP) == FALSE || log( runif(1,0,1) ) > ( newPP - oldPP ) ) {
# reject
tMassExtinction[idx] <- t
} else {
currentPP <- newPP
massExtinctionTimeAccepted <- massExtinctionTimeAccepted + 1
}
massExtinctionTimeTried <- massExtinctionTimeTried + 1
}
# change the value of an event
if ( length(pMassExtinction) > 0 ) {
# compute index of value that will be changed
idx <- floor( runif(1,0,length(pMassExtinction)) ) + 1
# store current value
v <- pMassExtinction[idx]
# compute current posterior
oldPP <- currentPP
# propose new value
vPrime <- rnorm(1, v, deltaMassExtinctionValue)
pMassExtinction[idx] <- vPrime
# compute new posterior
if ( vPrime > 0 && vPrime < 1 ) {
newPP <- prior(lambdaChangeTimes,lambda,muChangeTimes,mu,phiChangeTimes,phi,tMassExtinction,pMassExtinction) + likelihood(lambda,lambdaChangeTimes,mu,muChangeTimes,phi,phiChangeTimes,tMassExtinction,pMassExtinction)
} else {
newPP <- -Inf
}
# accept/reject
if ( is.finite(newPP - oldPP) == FALSE || log( runif(1,0,1) ) > ( newPP - oldPP ) ) {
# reject
pMassExtinction[idx] <- v
} else {
currentPP <- newPP
massExtinctionValueAccepted <- massExtinctionValueAccepted + 1
}
massExtinctionValueTried <- massExtinctionValueTried + 1
}
if ( estimateNumberMassExtinctions == TRUE ) {
# We need to randomly pick a birth or death move
# Otherwise we might give birth and die every time
u <- runif(1,0,1)
if ( u > 0.5 ) {
# birth move
# compute current posterior
oldPP <- currentPP
# randomly pick a new time
t <- runif(1,0,AGE)
# randomly pick a new value
v <- rbeta(1,pMassExtinctionPriorShape1,pMassExtinctionPriorShape2)
# construct the new parameters
tMassExtinction[length(tMassExtinction)+1] <- t
pMassExtinction[length(pMassExtinction)+1] <- v
# compute new posterior
newPP <- prior(lambdaChangeTimes,lambda,muChangeTimes,mu,phiChangeTimes,phi,tMassExtinction,pMassExtinction) + likelihood(lambda,lambdaChangeTimes,mu,muChangeTimes,phi,phiChangeTimes,tMassExtinction,pMassExtinction)
# compute proposal ratio
kPrime <- length(pMassExtinction)
pr <- 1 / ( dunif(t,0,AGE) * dbeta(v,pMassExtinctionPriorShape1,pMassExtinctionPriorShape2) )
# accept/reject
if ( is.finite(newPP - oldPP + log(pr)) == FALSE || log( runif(1,0,1) ) > ( newPP - oldPP + log( pr ) ) ) {
# reject
tMassExtinction <- tMassExtinction[-length(tMassExtinction)]
pMassExtinction <- pMassExtinction[-length(pMassExtinction)]
} else {
currentPP <- newPP
}
} else {
# death move
if ( length(tMassExtinction) > 0 ) {
# compute current posterior
oldPP <- currentPP
# randomly pick an index
idx <- floor( runif(1,0,length(tMassExtinction)) ) + 1
# store the current value
t <- tMassExtinction
v <- pMassExtinction
# construct the new parameters
tMassExtinction <- tMassExtinction[-idx]
pMassExtinction <- pMassExtinction[-idx]
# compute new posterior
newPP <- prior(lambdaChangeTimes,lambda,muChangeTimes,mu,phiChangeTimes,phi,tMassExtinction,pMassExtinction) + likelihood(lambda,lambdaChangeTimes,mu,muChangeTimes,phi,phiChangeTimes,tMassExtinction,pMassExtinction)
# compute proposal ratio
kPrime <- length(tMassExtinction) + 1
pr2 <- dunif(t[idx],0,AGE) * dbeta(v[idx],pMassExtinctionPriorShape1,pMassExtinctionPriorShape2)
# accept/reject
if ( is.finite(newPP - oldPP + log(pr2)) == FALSE || log( runif(1,0,1) ) > ( newPP - oldPP + log( pr2 ) ) ) {
# reject
tMassExtinction <- t
pMassExtinction <- v
} else {
currentPP <- newPP
}
}
}
}
#################
## Auto-Tuning ##
#################
# if we use adaptive MCMC we might have to modify our parameters
if ( ADAPTIVE == TRUE & SAMPLE == FALSE ) {
if ( i %% OPTIMIZATION_FREQUENCY == 0 ) {
# delta for time of lambda change events
if ( lambdaTimeTried > 0 ) {
rate <- lambdaTimeAccepted / lambdaTimeTried
if ( rate > 0.44 ) {
deltaLambdaTime <- deltaLambdaTime * (1.0 + ((rate-0.44)/0.56) )
} else {
deltaLambdaTime <- deltaLambdaTime / (2.0 - rate/0.44 )
}
lambdaTimeTried <- 0
lambdaTimeAccepted <- 0
}
# delta for value of lambda
rate <- lambdaValueAccepted / lambdaValueTried
if ( rate > 0.44 ) {
deltaLambdaValue <- deltaLambdaValue * (1.0 + ((rate-0.44)/0.56) )
} else {
deltaLambdaValue <- deltaLambdaValue / (2.0 - rate/0.44 )
}
lambdaValueTried <- 0
lambdaValueAccepted <- 0
# delta for time of mu change events
if ( muTimeTried > 0 ) {
rate <- muTimeAccepted / muTimeTried
if ( rate > 0.44 ) {
deltaMuTime <- deltaMuTime * (1.0 + ((rate-0.44)/0.56) )
} else {
deltaMuTime <- deltaMuTime / (2.0 - rate/0.44 )
}
muTimeTried <- 0
muTimeAccepted <- 0
}
# delta for values of mu
rate <- muValueAccepted / muValueTried
if ( rate > 0.44 ) {
deltaMuValue <- deltaMuValue * (1.0 + ((rate-0.44)/0.56) )
} else {
deltaMuValue <- deltaMuValue / (2.0 - rate/0.44 )
}
muValueTried <- 0
muValueAccepted <- 0
# delta for time of phi change events
if ( phiTimeTried > 0 ) {
rate <- phiTimeAccepted / phiTimeTried
if ( rate > 0.44 ) {
deltaPhiTime <- deltaPhiTime * (1.0 + ((rate-0.44)/0.56) )
} else {
deltaPhiTime <- deltaPhiTime / (2.0 - rate/0.44 )
}
phiTimeTried <- 0
phiTimeAccepted <- 0
}
# delta for values of phi
rate <- phiValueAccepted / phiValueTried
if ( rate > 0.44 ) {
deltaPhiValue <- deltaPhiValue * (1.0 + ((rate-0.44)/0.56) )
} else {
deltaPhiValue <- deltaPhiValue / (2.0 - rate/0.44 )
}
phiValueTried <- 0
phiValueAccepted <- 0
# delta for time of massExtinction change events
if ( massExtinctionTimeTried > 0 ) {
rate <- massExtinctionTimeAccepted / massExtinctionTimeTried
if ( rate > 0.44 ) {
deltaMassExtinctionTime <- deltaMassExtinctionTime * (1.0 + ((rate-0.44)/0.56) )
} else {
deltaMassExtinctionTime <- deltaMassExtinctionTime / (2.0 - rate/0.44 )
}
massExtinctionTimeTried <- 0
massExtinctionTimeAccepted <- 0
}
# delta for values of mu
if ( massExtinctionValueTried > 0 ) {
rate <- massExtinctionValueAccepted / massExtinctionValueTried
if ( rate > 0.44 ) {
deltaMassExtinctionValue <- deltaMassExtinctionValue * (1.0 + ((rate-0.44)/0.56) )
} else {
deltaMassExtinctionValue <- deltaMassExtinctionValue / (2.0 - rate/0.44 )
}
massExtinctionValueTried <- 0
massExtinctionValueAccepted <- 0
}
}
}
# store the values
if ( i %% THINNING == 0 & SAMPLE ) {
sampleIndex <- length(kLambda)+1
kLambda[sampleIndex] <- length(lambda)
speciationRateValues[[sampleIndex]] <- lambda
speciationRateChangeTimes[[sampleIndex]] <- lambdaChangeTimes
kMu[sampleIndex] <- length(mu)
extinctionRateValues[[sampleIndex]] <- mu
extinctionRateChangeTimes[[sampleIndex]] <- muChangeTimes
kPhi[sampleIndex] <- length(phi)
fossilizationRateValues[[sampleIndex]] <- phi
fossilizationRateChangeTimes[[sampleIndex]] <- phiChangeTimes
kMassExtinction[sampleIndex] <- length(pMassExtinction)
survivalProbability[[sampleIndex]] <- pMassExtinction
massExtinctionTime[[sampleIndex]] <- tMassExtinction
posterior[sampleIndex] <- currentPP
# write.table(speciationRateValues,"SpeciationRates.txt",sep="\t")
cat(lambda,sep="\t",file="SpeciationRates.txt",append=TRUE)
cat("\n",file="SpeciationRates.txt",append=TRUE)
# write.table(speciationRateChangeTimes,"SpeciationRateChanges.txt",sep="\t")
cat(lambdaChangeTimes,sep="\t",file="SpeciationRateChanges.txt",append=TRUE)
cat("\n",file="SpeciationRateChanges.txt",append=TRUE)
# write.table(extinctionRateValues,"ExtinctionRates.txt",sep="\t")
cat(mu,sep="\t",file="ExtinctionRates.txt",append=TRUE)
cat("\n",file="ExtinctionRates.txt",append=TRUE)
# write.table(extinctionRateChangeTimes,"ExtinctionRateChanges.txt",sep="\t")
cat(muChangeTimes,sep="\t",file="ExtinctionRateChanges.txt",append=TRUE)
cat("\n",file="ExtinctionRateChanges.txt",append=TRUE)
# write.table(fossilizationRateValues,"FossilizationRates.txt",sep="\t")
cat(phi,sep="\t",file="FossilizationRates.txt",append=TRUE)
cat("\n",file="FossilizationRates.txt",append=TRUE)
# write.table(fossilizationRateChangeTimes,"FossilizationRateChanges.txt",sep="\t")
cat(phiChangeTimes,sep="\t",file="FossilizationRateChanges.txt",append=TRUE)
cat("\n",file="FossilizationRateChanges.txt",append=TRUE)
# write.table(survivalProbability,"SurvivalProbabilities.txt",sep="\t")
cat(pMassExtinction,sep="\t",file="SurvivalProbabilities.txt",append=TRUE)
cat("\n",file="SurvivalProbabilities.txt",append=TRUE)
# write.table(massExtinctionTime,"MassExtinctionTimes.txt",sep="\t")
cat(tMassExtinction,sep="\t",file="MassExtinctionTimes.txt",append=TRUE)
cat("\n",file="MassExtinctionTimes.txt",append=TRUE)
nSamples <- sampleIndex
write.table(list(Iteration=nSamples*THINNING,posterior=posterior[nSamples],NumSpeciation=kLambda[nSamples],numExtinction=kMu[nSamples],numFossilization=kPhi[nSamples],numMassExtinctions=kMassExtinction[nSamples]),"samples_numCategories.txt", sep="\t", row.names = FALSE, append=TRUE, quote = FALSE, col.names = FALSE)
}
# Progress by convergence
if ( i %% CONVERGENCE_FREQUENCY == 0 & SAMPLE & length(posterior) > 2 ) {
samples <- length(posterior)
if ( estimateNumberSpeciationRateChanges == TRUE ) {
spectralDensityLambda <- spectrum0.ar(kLambda)$spec
if ( spectralDensityLambda > 0 ) {
ess_kLambda <- (var(kLambda) * samples/spectralDensityLambda)
} else {
ess_kLambda <- Inf
}
} else {
ess_kLambda <- Inf
}
if ( estimateNumberExtinctionRateChanges == TRUE ) {
spectralDensityMu <- spectrum0.ar(kMu)$spec
if ( spectralDensityMu > 0 ) {
ess_kMu <- (var(kMu) * samples/spectralDensityMu)
} else {
ess_kMu <- Inf
}
} else {
ess_kMu <- Inf
}
if ( estimateNumberFossilizationRateChanges == TRUE ) {
spectralDensityPhi <- spectrum0.ar(kPhi)$spec
if ( spectralDensityPhi > 0 ) {
ess_kPhi <- (var(kPhi) * samples/spectralDensityPhi)
} else {
ess_kPhi <- Inf
}
} else {
ess_kPhi <- Inf
}
spec <- c()
for ( j in 1:length(speciationRateValues) ) {
spec[j] <- speciationRateValues[[j]][1]
}
spectralDensity <- spectrum0.ar(spec)$spec
if ( spectralDensity > 0 ) {
ess_spec <- (var(spec) * samples/spectralDensity)
} else {
ess_spec <- Inf
}
ext <- c()
for ( j in 1:length(extinctionRateValues) ) {
ext[j] <- extinctionRateValues[[j]][1]
}
spectralDensity <- spectrum0.ar(ext)$spec
if ( spectralDensity > 0 ) {
ess_ext <- (var(ext) * samples/spectralDensity)
} else {
ess_ext <- Inf
}
fos <- c()
for ( j in 1:length(fossilizationRateValues) ) {
fos[j] <- fossilizationRateValues[[j]][1]
}
spectralDensity <- spectrum0.ar(fos)$spec
if ( spectralDensity > 0 ) {
ess_fos <- (var(fos) * samples/spectralDensity)
} else {
ess_fos <- Inf
}
if ( estimateNumberMassExtinctions == TRUE ) {
ess_kMassExtinctions <- (var(kMassExtinction) * samples/spectrum0.ar(kMassExtinction)$spec)
} else {
ess_kMassExtinctions <- Inf
}
min.ess <- min(c(ess_kLambda,ess_kMu,ess_kPhi,ess_spec,ess_ext,ess_fos,ess_kMassExtinctions))
}
if ( SAMPLE ) {
# Progress by iteration
iteration.progress <- i / MAX_ITERATIONS
# Progress by time
time.progress <- ( (Sys.time() - startTime) / MAX_TIME )
# Progress by ESS
ess.progress <- min.ess / MIN_ESS
# Max progress
max.progress <- max(c(iteration.progress,time.progress,ess.progress),na.rm=TRUE)
if( length(max.progress) == 0 ){
max.progress <- 0
}
if ( verbose ) {
setTxtProgressBar(bar,pmin(max.progress,1))
}
if( max.progress >= 1 ) {
finished <- TRUE
}
} else {
if ( verbose ){
setTxtProgressBar(bar,pmin(1,i/burn))
}
}
i <- i + 1
if ( i == burn & !SAMPLE ) {
startTime <- Sys.time()
i <- 1
SAMPLE <- TRUE
if ( verbose ) {
cat("\n\nRunning the chain ... \n")
cat("0--------25--------50--------75--------100\n")
bar <- txtProgressBar(style=1,width=42)
}
}
}
cat("\n")
setwd(file.path(orgDir))
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.