Description Usage Arguments Value Author(s) References Examples
View source: R/tess.PosteriorPrediction.R
tess.PosteriorPrediction calls the simulation function exactly once for each sampled parameter combination. In that way, posterior predictive simulations can be obtained which then in turn can be used to compute summary statistics based on these posterior predictive simulations. Fore more information see the vignette.
1 | tess.PosteriorPrediction(simulationFunction,parameters,burnin)
|
simulationFunction |
The simulation function which will be called internally by simulationFunction(parameters). |
parameters |
A matrix of parameters where the rows represent samples of parameters and the column the different parameters. |
burnin |
The fraction of samples to be discarded as burnin. This is 0.25 by default |
Returns samples simulated from the posterior predictive distribution.
Sebastian Hoehna
S. Hoehna: Fast simulation of reconstructed phylogenies under global, time-dependent birth-death processes. 2013, Bioinformatics, 29:1367-1374
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 | # We first run an MCMC to obtain samples from the posterior distribution
# and then simulate the posterior predictive distribution.
# The bird phylogeny as the test data set
data(cettiidae)
times <- as.numeric( branching.times(cettiidae) )
# The log-likelihood function
likelihood <- function(params) {
# We use the parameters as diversification rate and turnover rate.
# Thus we need to transform first
b <- params[1] + params[2]
d <- params[2]
lnl <- tess.likelihood(times,b,d,samplingProbability=1.0,log=TRUE)
return (lnl)
}
prior_diversification <- function(x) { dexp(x,rate=0.1,log=TRUE) }
prior_turnover <- function(x) { dexp(x,rate=0.1,log=TRUE) }
priors <- c(prior_diversification,prior_turnover)
# Note, the number of iterations and the burnin is too small here
# and should be adapted for real analyses
samples <- tess.mcmc(likelihood,priors,c(1,0.1),c(TRUE,TRUE),c(0.1,0.1),10,10)
tmrca <- max(branching.times(cettiidae))
# The simulation function
sim <- function(params) {
# We use the parameters as diversification rate and turnover rate.
# Thus we need to transform first
b <- params[1] + params[2]
d <- params[2]
tree <- tess.sim.age(n=1,age=tmrca,b,d,samplingProbability=1.0)[[1]]
return (tree)
}
trees <- tess.PosteriorPrediction(sim,samples)
# compute the posterior predictive test statistic
ppt <- tess.PosteriorPredictiveTest(trees,cettiidae,gammaStat)
# get the p-value of the observed test-statistic
ppt[[2]]
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.