| Title: | Bayesian Fitting of Ornstein-Uhlenbeck Models to Phylogenies |
|---|---|
| Description: | Tools for fitting and simulating multi-optima Ornstein-Uhlenbeck models to phylogenetic comparative data using Bayesian reversible-jump methods. |
| Authors: | Josef C. Uyeda, Jon Eastman and Luke Harmon |
| Maintainer: | Josef C. Uyeda <[email protected]> |
| License: | GPL (>= 2) |
| Version: | 2.3.0 |
| Built: | 2024-10-26 04:56:46 UTC |
| Source: | https://github.com/uyedaj/bayou |
A package for inferring adaptive evolution to phylogenetic comparative data using Bayesian reversible-jump estimation of multi-optima Ornstein-Uhlenbeck models.
bayou-package
Josef C Uyeda
Function for checking parameter lists, prior and models are consistent and error-free
bayou.checkModel( pars = NULL, tree, dat, pred = NULL, SE = 0, prior, model = "OU", autofix = TRUE )bayou.checkModel( pars = NULL, tree, dat, pred = NULL, SE = 0, prior, model = "OU", autofix = TRUE )
pars |
A list of parameters that will be specified as starting parameter |
tree |
An object of class “phylo” |
dat |
A named data vector that matches the tip lables in the provided tree |
pred |
A matrix or data frame with named columns with predictor data represented in the specified formula |
SE |
The standard error of the data. Either a single value applied to all the data, or a vector of length(dat). |
prior |
A prior function made using make.prior |
model |
Either one of c("OU", "QG" or "OUrepar") or a list specifying the model to be used. |
autofix |
A logical that indicates whether certain errors should be automatically fixed. |
A series of checks are performed, run internally within bayou.makeMCMC, but can also be run on provided inputs prior to this. Errors are reported.
If autofix == TRUE, then the following errors will be automatically corrected:
Branch lengths == 0; any branches of length 0 will be given length .Machine$double.eps is.binary(tree) == FALSE; runs multi2di pars do not match prior$fixed; parameters are resimulated from prior
A list of results of the checks and if 'autofix==TRUE', then ..$autofixed returns a list of all the input elements, with corrections.
Function for calculating likelihood of an OU model in bayou using the threepoint algorithm
bayou.lik(pars, cache, X, model = "OU")bayou.lik(pars, cache, X, model = "OU")
pars |
A list of parameters to calculate the likelihood |
cache |
A bayou cache object generated using .prepare.ou.univariate |
X |
A named vector giving the tip data |
model |
Parameterization of the OU model. Either "OU", "QG" or "OUrepar". |
This function implements the algorithm of Ho and Ane (2014) implemented
in the package phylolm for the OUfixedRoot model. It is faster
than the equivalent pruning algorithm in geiger, and can be used on non-
ultrametric trees (unlike OU.lik, which is based on the pruning algorithm in
geiger).
Runs a reversible-jump Markov chain Monte Carlo on continuous phenotypic data on a phylogeny, sampling possible shift locations and shift magnitudes, and shift numbers.
bayou.makeMCMC( tree, dat, pred = NULL, SE = 0, model = "OU", prior, samp = 10, chunk = 100, control = NULL, tuning = NULL, file.dir = tempdir(), plot.freq = 500, outname = "bayou", plot.fn = phenogram, ticker.freq = 1000, startpar = NULL, moves = NULL, control.weights = NULL, lik.fn = NULL, perform.checks = TRUE )bayou.makeMCMC( tree, dat, pred = NULL, SE = 0, model = "OU", prior, samp = 10, chunk = 100, control = NULL, tuning = NULL, file.dir = tempdir(), plot.freq = 500, outname = "bayou", plot.fn = phenogram, ticker.freq = 1000, startpar = NULL, moves = NULL, control.weights = NULL, lik.fn = NULL, perform.checks = TRUE )
tree |
a phylogenetic tree of class 'phylo' |
dat |
a named vector of continuous trait values matching the tips in tree |
pred |
A matrix or data frame with named columns with predictor data represented in the specified formula |
SE |
The standard error of the data. Either a single value applied to all the data, or a vector of length(dat). |
model |
The parameterization of the OU model used. Either "OU" for standard parameterization with alpha and sigma^2; "OUrepar" for phylogenetic half-life and stationary variance (Vy), or "QG" for the Lande model, with parameters h^2 (heritability), P (phenotypic variance), omega^2 (width of adaptive landscape), and Ne (effective population size) |
prior |
A prior function of class 'priorFn' that gives the prior distribution of all parameters |
samp |
The frequency at which Markov samples are retained |
chunk |
The number of samples retained in memory before being written to a file |
control |
A list providing a control object governing how often and which proposals are used |
tuning |
A named vector that governs how liberal or conservative proposals are that equals the number of proposal mechanisms. |
file.dir |
If a character string, then results are written to that working directory. If NULL, then results are not saved to files, but instead held in memory. Default is 'tempdir()', which writes to an R temporary directory. |
plot.freq |
How often plots should be made during the mcmc. If NULL, then plots are not produced |
outname |
The prefix given to files created by the mcmc |
plot.fn |
Function used in plotting, defaults to phytools::phenogram |
ticker.freq |
How often a summary log should be printed to the screen |
startpar |
A list with the starting parameters for the mcmc. If NULL, starting parameters are simulated from the prior distribution |
moves |
A named list providing the proposal functions to be used in the mcmc. Names correspond to the parameters to be modified in the parameter list. See 'details' for default values. |
control.weights |
A named vector providing the relative frequency each proposal mechanism is to be used during the mcmc |
lik.fn |
Likelihood function to be evaluated. Defaults to |
perform.checks |
A logical indicating whether to use bayou.checkModel to validate model inputs. |
By default, the alpha, sig2 (and various reparameterizations of these parameters) are adjusted with multiplier proposals, theta are adjusted with sliding window proposals, and the number of shifts is adjusted by splitting and merging, as well as sliding the shifts both within and between branches. Allowed shift locations are specified by the prior function (see make.prior()).
bayou2OUwie Converts a bayou formatted parameter list into OUwie formatted tree and data table that can be analyzed in OUwie
bayou2OUwie(pars, tree, dat)bayou2OUwie(pars, tree, dat)
pars |
A list with parameter values specifying |
tree |
A phylogenetic tree |
dat |
A vector of tip states |
A list with an OUwie formatted tree with mapped regimes and an OUwie formatted data table
cdpois calculates the probability density of a value k from a
Poisson distribution with a maximum kmax. rdpois draws random
numbers from a conditional Poisson distribution.
cdpois(k, lambda, kmax, log = TRUE) rdpois(n, lambda, kmax, ...)cdpois(k, lambda, kmax, log = TRUE) rdpois(n, lambda, kmax, ...)
k |
random variable value |
lambda |
rate parameter of the Poisson distribution |
kmax |
maximum value of the conditional Poisson distribution |
log |
log transformed density |
n |
number of samples to draw |
... |
additional parameters passed to |
cdpois(10,1,10) cdpois(11,1,10) #rdpois(5,10,10)cdpois(10,1,10) cdpois(11,1,10) #rdpois(5,10,10)
Combine mcmc chains
combine.chains(chain.list, thin = 1, burnin.prop = 0)combine.chains(chain.list, thin = 1, burnin.prop = 0)
chain.list |
The first chain to be combined |
thin |
A number or vector specifying the thinning interval to be used. If a single value, then the same proportion will be applied to all chains. |
burnin.prop |
A number or vector giving the proportion of burnin from each chain to be discarded. If a single value, then the same proportion will be applied to all chains. |
A combined bayouMCMC chain
This function simulates data for a given set of parameter values.
dataSim(pars, model, tree, map.type = "pars", SE = 0, phenogram = TRUE, ...)dataSim(pars, model, tree, map.type = "pars", SE = 0, phenogram = TRUE, ...)
pars |
A bayou formated parameter list |
model |
The type of model specified by the parameter list (either "OU", "OUrepar" or "QG"). |
tree |
A tree of class 'phylo' |
map.type |
Either "pars" if the regimes are taken from the parameter list, or "simmap" if taken from the stored simmap in the tree |
SE |
A single value or vector equal to the number of tips specifying the measurement error that should be simulated at the tips |
phenogram |
A logical indicating whether or not the simulated data should be plotted as a phenogram |
... |
Optional parameters passed to |
dataSim Simulates data for a given bayou model and parameter set
dhalfcauchy returns the probability density for a half-Cauchy distribution
dhalfcauchy(x, scale = 25, log = FALSE) phalfcauchy(q, scale = 25) qhalfcauchy(p, scale = 25) rhalfcauchy(n, scale = 25)dhalfcauchy(x, scale = 25, log = FALSE) phalfcauchy(q, scale = 25) qhalfcauchy(p, scale = 25) rhalfcauchy(n, scale = 25)
x |
A parameter value for which the density should be calculated |
scale |
The scale parameter of the half-Cauchy distributoin |
log |
A logical indicating whether the log density should be returned |
q |
A vector of quantiles |
p |
A vector of probabilities |
n |
The number of observations |
Since unequal probabilities are incorporated in calculating the
density via dsb, all branches are assumed to be of unit length.
Thus, the dloc function simply returns 0 if log=TRUE and 1 if log=FALSE.
dloc(loc, min = 0, max = 1, log = TRUE) rloc(k, min = 0, max = 1)dloc(loc, min = 0, max = 1, log = TRUE) rloc(k, min = 0, max = 1)
loc |
The location of the shift along the branch |
min |
The minimum position on the branch the shift can take |
max |
The maximum position on the branch the shift can take |
log |
A logical indicating whether the log density should be returned |
k |
The number of shifts to return along a branch |
dloc calculates the probability of a shift occuring at a given
location along the branch assuming a uniform distribution of unit length
rloc randomly generates the location of a shift along the branch
This function provides a means to specify the prior for the location of shifts across the phylogeny. Certain combinations are not allowed. For example, a maximum shift number of Inf on one branch cannot be combined with a maximum shift number of 1 on another. Thus, bmax must be either a vector of 0's and Inf's or a vector of 0's and 1's. Also, if bmax == 1, then all probabilities must be equal, as bayou cannot sample unequal probabilities without replacement.
dsb(sb, ntips = ntips, bmax = 1, prob = 1, log = TRUE) rsb(k, ntips = ntips, bmax = 1, prob = 1, log = TRUE)dsb(sb, ntips = ntips, bmax = 1, prob = 1, log = TRUE) rsb(k, ntips = ntips, bmax = 1, prob = 1, log = TRUE)
sb |
A vector giving the branch numbers (for a post-ordered tree) |
ntips |
The number of tips in the phylogeny |
bmax |
A single integer or a vector of integers equal to the number of branches in the phylogeny indicating the maximum number of shifts allowable in the phylogeny. Can take values 0, 1 and Inf. |
prob |
A single value or a vector of values equal to the number of branches in the phylogeny indicating the probability that a randomly selected shift will lie on this branch. Can take any positive value, values need not sum to 1 (they will be scaled to sum to 1) |
log |
A logical indicating whether the log probability should be returned. Default is 'TRUE' |
k |
The number of shifts to randomly draw from the distribution |
dsb calculates the probability of a particular arrangement of shifts
for a given set of assumptions.
The log density of the particular number and arrangement of shifts.
n=10 tree <- sim.bdtree(n=n) tree <- reorder(tree, "postorder") nbranch <- 2*n-2 sb <- c(1,2, 2, 3) # Allow any number of shifts on each branch, with probability # proportional to branch length dsb(sb, ntips=n, bmax=Inf, prob=tree$edge.length) # Disallow shifts on the first branch, returns -Inf because sb[1] = 1 dsb(sb, ntips=n, bmax=c(0, rep(1, nbranch-1)), prob=tree$edge.length) # Set maximum number of shifts to 1, returns -Inf because two shifts # are on branch 2 dsb(sb, ntips=n, bmax=1, prob=1) # Generate a random set of k branches rsb(5, ntips=n, bmax=Inf, prob=tree$edge.length)n=10 tree <- sim.bdtree(n=n) tree <- reorder(tree, "postorder") nbranch <- 2*n-2 sb <- c(1,2, 2, 3) # Allow any number of shifts on each branch, with probability # proportional to branch length dsb(sb, ntips=n, bmax=Inf, prob=tree$edge.length) # Disallow shifts on the first branch, returns -Inf because sb[1] = 1 dsb(sb, ntips=n, bmax=c(0, rep(1, nbranch-1)), prob=tree$edge.length) # Set maximum number of shifts to 1, returns -Inf because two shifts # are on branch 2 dsb(sb, ntips=n, bmax=1, prob=1) # Generate a random set of k branches rsb(5, ntips=n, bmax=Inf, prob=tree$edge.length)
Calculate Gelman's R statistic
gelman.R(parameter, chain1, chain2, freq = 20, start = 1, plot = TRUE, ...)gelman.R(parameter, chain1, chain2, freq = 20, start = 1, plot = TRUE, ...)
parameter |
The name or number of the parameter to calculate the statistic on |
chain1 |
The first bayouMCMC chain |
chain2 |
The second bayouMCMC chain |
freq |
The interval between which the diagnostic is calculated |
start |
The first sample to calculate the diagnostic at |
plot |
A logical indicating whether the results should be plotted |
... |
Optional arguments passed to |
This is a convenience function for mapping regimes interactively on the phylogeny. The method locates the nearest branch to where the cursor is clicked on the plot and records the branch number and the location selected on the branch.
identifyBranches(tree, n, fixed.loc = TRUE, plot.simmap = TRUE)identifyBranches(tree, n, fixed.loc = TRUE, plot.simmap = TRUE)
tree |
An object of class 'phylo' |
n |
The number of shifts to map interactively onto the phylogeny |
fixed.loc |
A logical indicating whether the exact location on the branch should be returned, or the shift will be free to move along the branch |
plot.simmap |
A logical indicating whether the resulting painting of regimes should be plotted following the selection shift location. |
identifyBranches opens an interactive phylogeny plot that allows the user to specify the location
of shifts in a phylogenetic tree.
Returns a list with elements "sb" which contains the branch numbers of all selected branches with length "n". If "fixed.loc=TRUE", then the list also contains a vector "loc" which contains the location of the selected shifts along the branch.
load.bayou loads a bayouFit object that was created using bayou.mcmc()
load.bayou(bayouFit, saveRDS = TRUE, file = NULL, cleanup = FALSE, ref = FALSE)load.bayou(bayouFit, saveRDS = TRUE, file = NULL, cleanup = FALSE, ref = FALSE)
bayouFit |
An object of class |
saveRDS |
A logical indicating whether the resulting chains should be saved as an *.rds file |
file |
An optional filename (possibly including path) for the saved *.rds file |
cleanup |
A logical indicating whether the files produced by |
ref |
A logical indicating whether a reference function is also in the output |
If both save.Rdata is FALSE and cleanup is TRUE, then load.bayou will trigger a
warning and ask for confirmation. In this case, if the results of load.bayou() are not stored in an object,
the results of the MCMC run will be permanently deleted.
## Not run: data(chelonia) tree <- chelonia$phy dat <- chelonia$dat prior <- make.prior(tree) fit <- bayou.mcmc(tree, dat, model="OU", prior=prior, new.dir=TRUE, ngen=5000) chain <- load.bayou(fit, save.Rdata=FALSE, cleanup=TRUE) plot(chain) ## End(Not run)## Not run: data(chelonia) tree <- chelonia$phy dat <- chelonia$dat prior <- make.prior(tree) fit <- bayou.mcmc(tree, dat, model="OU", prior=prior, new.dir=TRUE, ngen=5000) chain <- load.bayou(fit, save.Rdata=FALSE, cleanup=TRUE) plot(chain) ## End(Not run)
Return a posterior of shift locations
Lposterior(chain, tree, burnin = 0, simpar = NULL, mag = TRUE)Lposterior(chain, tree, burnin = 0, simpar = NULL, mag = TRUE)
chain |
A bayouMCMC chain |
tree |
A tree of class 'phylo' |
burnin |
A value giving the burnin proportion of the chain to be discarded |
simpar |
An optional bayou formatted parameter list giving the true values (if data were simulated) |
mag |
A logical indicating whether the average magnitude of the shifts should be returned |
A data frame with rows corresponding to postordered branches. pp indicates the
posterior probability of the branch containing a shift. magnitude of theta2 gives the average
value of the new optima after a shift. naive SE of theta2 gives the standard error of the new optima
not accounting for autocorrelation in the MCMC and rel location gives the average relative location
of the shift on the branch (between 0 and 1 for each branch).
This function generates a power posterior function for estimation of marginal likelihood using the stepping stone method
make.powerposteriorFn(Bk, priorFn, refFn, model)make.powerposteriorFn(Bk, priorFn, refFn, model)
Bk |
The sequence of steps to be taken from the reference function to the posterior |
priorFn |
The prior function to be used in marginal likelihood estimation |
refFn |
The reference function generated using |
model |
A string specifying the model type ("OU", "OUrepar", "QG") or a model parameter list |
For use in stepping stone estimation of the marginal likelihood using the method of Fan et al. (2011).
A function of class "powerposteriorFn" that returns a list of four values: result (the log density of the power posterior),
lik (the log likelihood), prior (the log prior), ref the log reference density.
This function generates a prior function to be used for bayou according to user specifications.
make.prior( tree, dists = list(), param = list(), fixed = list(), plot.prior = TRUE, model = "OU" )make.prior( tree, dists = list(), param = list(), fixed = list(), plot.prior = TRUE, model = "OU" )
tree |
A tree object of class "phylo" |
dists |
A list providing the function names of the distribution functions describing the prior distributions of parameters (see details). If no distributions are provided for a parameter, default values are given. Note that the names are provided as text strings, not the functions themselves. |
param |
A list providing the parameter values of the prior distributions (see details). |
fixed |
A list of parameters that are to be fixed at provided values. These are removed from calculation of the prior value. |
plot.prior |
A logical indicating whether the prior distributions should be plotted. |
model |
One of three specifications of the OU parameterization used.
Takes values |
Default distributions and parameter values are given as follows:
OU: list(dists=list("dalpha"="dlnorm","dsig2"="dlnorm",
"dk"="cdpois","dtheta"="dnorm","dsb"="dsb","dloc"="dunif"),
param=list("dalpha"=list(),"dsig2"=list(),"dtheta"=list(),
"dk"=list(lambda=1,kmax=2*ntips-2),"dloc"=list(min=0,max=1),"dsb"=list()))
QG: list(dists=list("dh2"="dbeta","dP"="dlnorm","dw2"="dlnorm","dNe"="dlnorm",
"dk"="cdpois","dtheta"="dnorm","dsb"="dsb","dloc"="dunif"),
param=list("dh2"=list(shape1=1,shape2=1),"dP"=list(),"dw2"=list(),"dNe"=list(),"dtheta"=list(),
"dk"=list(lambda=1,kmax=2*ntips-2),"dloc"=list(min=0,max=1),"dsb"=list()))
OUrepar: list(dists=list("dhalflife"="dlnorm","dVy"="dlnorm",
"dk"="cdpois","dtheta"="dnorm","dsb"="dsb","dloc"="dunif"),
param=list("dhalflife"=list("meanlog"=0.25,"sdlog"=1.5),"dVy"=list("meanlog"=1,"sdlog"=2),
"dk"=list(lambda=1,kmax=2*ntips-2),"dtheta"=list(),"dloc"=list(min=0,max=1)),"dsb"=list())
dalpha, dsig2, dh2, dP, dw2, dNe, dhalflife, and dVy must be positive continuous distributions and provide the parameters used to calculate alpha and sigma^2 of the OU model.
dtheta must be continuous and describes the prior distribution of the optima. dk is the prior distribution for the number of shifts. For Poisson and conditional Poisson (cdpois) are provided
the parameter lambda, which provides the total number of shifts expected on the tree (not the rate per unit branch length). Otherwise, dk can take any positive, discrete distribution.
dsb indicates the prior probability of a given set of branches having shifts, and is generally specified by the "dsb" function in the bayou package. See the documentation for dsb for specifying the number
of shifts allowed per branch, the probability of a branch having a shift, and specifying constraints on where shifts can occur."dloc" indicates the prior probability of the location of a shift within
a single branch. Currently, all locations are given uniform density. All distributions are set to return log-transformed probability densities.
returns a prior function of class "priorFn" that calculates the log prior density for a set of parameter values provided in a list with correctly named values.
## Load data data(chelonia) tree <- chelonia$phy dat <- chelonia$dat #Create a prior that allows only one shift per branch with equal probability #across branches prior <- make.prior(tree, dists=list(dalpha="dlnorm", dsig2="dlnorm", dsb="dsb", dk="cdpois", dtheta="dnorm"), param=list(dalpha=list(meanlog=-5, sdlog=2), dsig2=list(meanlog=-1, sdlog=5), dk=list(lambda=15, kmax=200), dsb=list(bmax=1,prob=1), dtheta=list(mean=mean(dat), sd=2))) #Evaluate some parameter sets pars1 <- list(alpha=0.1, sig2=0.1, k=5, ntheta=6, theta=rnorm(6, mean(dat), 2), sb=c(32, 53, 110, 350, 439), loc=rep(0.1, 5), t2=2:6) pars2 <- list(alpha=0.1, sig2=0.1, k=5, ntheta=6, theta=rnorm(6, mean(dat), 2), sb=c(43, 43, 432, 20, 448), loc=rep(0.1, 5), t2=2:6) prior(pars1) prior(pars2) #-Inf because two shifts on one branch #Create a prior that allows any number of shifts along each branch with probability proportional #to branch length prior <- make.prior(tree, dists=list(dalpha="dlnorm", dsig2="dlnorm", dsb="dsb", dk="cdpois", dtheta="dnorm"), param=list(dalpha=list(meanlog=-5, sdlog=2), dsig2=list(meanlog=-1, sdlog=5), dk=list(lambda=15, kmax=200), dsb=list(bmax=Inf,prob=tree$edge.length), dtheta=list(mean=mean(dat), sd=2))) prior(pars1) prior(pars2) #Create a prior with fixed regime placement and sigma^2 value prior <- make.prior(tree, dists=list(dalpha="dlnorm", dsig2="fixed", dsb="fixed", dk="fixed", dtheta="dnorm", dloc="dunif"), param=list(dalpha=list(meanlog=-5, sdlog=2), dtheta=list(mean=mean(dat), sd=2)), fixed=list(sig2=1, k=3, ntheta=4, sb=c(447, 396, 29))) pars3 <- list(alpha=0.01, theta=rnorm(4, mean(dat), 2), loc=rep(0.1, 4)) prior(pars3) ##Return a list of functions used to calculate prior attributes(prior)$functions ##Return parameter values used in prior distribution attributes(prior)$parameters## Load data data(chelonia) tree <- chelonia$phy dat <- chelonia$dat #Create a prior that allows only one shift per branch with equal probability #across branches prior <- make.prior(tree, dists=list(dalpha="dlnorm", dsig2="dlnorm", dsb="dsb", dk="cdpois", dtheta="dnorm"), param=list(dalpha=list(meanlog=-5, sdlog=2), dsig2=list(meanlog=-1, sdlog=5), dk=list(lambda=15, kmax=200), dsb=list(bmax=1,prob=1), dtheta=list(mean=mean(dat), sd=2))) #Evaluate some parameter sets pars1 <- list(alpha=0.1, sig2=0.1, k=5, ntheta=6, theta=rnorm(6, mean(dat), 2), sb=c(32, 53, 110, 350, 439), loc=rep(0.1, 5), t2=2:6) pars2 <- list(alpha=0.1, sig2=0.1, k=5, ntheta=6, theta=rnorm(6, mean(dat), 2), sb=c(43, 43, 432, 20, 448), loc=rep(0.1, 5), t2=2:6) prior(pars1) prior(pars2) #-Inf because two shifts on one branch #Create a prior that allows any number of shifts along each branch with probability proportional #to branch length prior <- make.prior(tree, dists=list(dalpha="dlnorm", dsig2="dlnorm", dsb="dsb", dk="cdpois", dtheta="dnorm"), param=list(dalpha=list(meanlog=-5, sdlog=2), dsig2=list(meanlog=-1, sdlog=5), dk=list(lambda=15, kmax=200), dsb=list(bmax=Inf,prob=tree$edge.length), dtheta=list(mean=mean(dat), sd=2))) prior(pars1) prior(pars2) #Create a prior with fixed regime placement and sigma^2 value prior <- make.prior(tree, dists=list(dalpha="dlnorm", dsig2="fixed", dsb="fixed", dk="fixed", dtheta="dnorm", dloc="dunif"), param=list(dalpha=list(meanlog=-5, sdlog=2), dtheta=list(mean=mean(dat), sd=2)), fixed=list(sig2=1, k=3, ntheta=4, sb=c(447, 396, 29))) pars3 <- list(alpha=0.01, theta=rnorm(4, mean(dat), 2), loc=rep(0.1, 4)) prior(pars3) ##Return a list of functions used to calculate prior attributes(prior)$functions ##Return parameter values used in prior distribution attributes(prior)$parameters
This function generates a reference function from a mcmc chain for use in marginal likelihood estimation.
make.refFn(chain, model, priorFn, burnin = 0.3, plot = TRUE)make.refFn(chain, model, priorFn, burnin = 0.3, plot = TRUE)
chain |
An mcmc chain produced by |
model |
A string specifying the model ("OU", "QG", "OUrepar") or a model parameter list |
priorFn |
The prior function used to generate the mcmc chain |
burnin |
The proportion of the mcmc chain to be discarded when generating the reference function |
plot |
Logical indicating whether or not a plot should be created |
Distributions are fit to each mcmc chain and the best-fitting distribution is chosen as
the reference distribution for that parameter using the method of Fan et al. (2011). For positive
continuous parameters alpha, sigma^2, halflife, Vy, w2, Ne, Log-normal, exponential, gamma and weibull
distributions are fit. For continuous distributions theta, Normal, Cauchy and Logistic distributions
are fit. For discrete distributions, k, negative binomial, poisson and geometric distributions are fit.
Best-fitting distributions are determined by AIC.
Returns a reference function of class "refFn" that takes a parameter list and returns the log density
given the reference distribution. If plot=TRUE, a plot is produced showing the density of variable parameters
and the fitted distribution from the reference function (in red).
This function makes a bayou model object that can be used for customized allometric regression models.
makeBayouModel( f, rjpars, tree, dat, pred, prior, SE = 0, slopechange = "immediate", impute = NULL, startpar = NULL, moves = NULL, control.weights = NULL, D = NULL, shiftpars = c("sb", "loc", "t2"), model = "OU" )makeBayouModel( f, rjpars, tree, dat, pred, prior, SE = 0, slopechange = "immediate", impute = NULL, startpar = NULL, moves = NULL, control.weights = NULL, D = NULL, shiftpars = c("sb", "loc", "t2"), model = "OU" )
f |
A formula describing the relationship between the data and one or more predictors (use 'dat' for the dependent variable) |
rjpars |
A character vector of parameters to split at the mapped shifts on the tree |
tree |
A phylogenetic tree |
dat |
A named vector of trait data (dependent variable) |
pred |
A matrix or data frame with named columns with predictor data represented in the specified formula |
prior |
A prior function made by the 'make.prior' function |
SE |
A single value or vector of measurement error estimates |
slopechange |
"immediate", "alphaWeighted" or "fullPGLS" |
impute |
The name of a single predictor for which missing values will be imputed using BM (see details). Default is NULL. |
startpar |
An optional list of starting parameters for the model. If not provided, the model will simulate starting values from the prior function. |
moves |
An optional list of moves to be passed on to bayou.makeMCMC. |
control.weights |
An optional list of control weights to be passed on to bayou.makeMCMC. |
D |
A vector of tuning parameters to be passed on to bayou.makeMCMC. |
shiftpars |
The names of the parameters defining the map of shifts (for now, always c("sb", "loc", "t2")). |
model |
The parameterization of the OU model, either "OU", "OUrepar" or "QG". |
This function generates a list with the '$model', which provides the specifications of the regression model and '$startpar', which provides starting values to input into bayou.makeMCMC. Note that this model assumes that predictors immediately affect trait values at a shift. In other words, regardless of the past history of the predictor, only the current value affects the current expected trait value. This is only reasonable for allometric models, although it may be appropriate for other models if phylogenetic inertia is very low (short half-lives).
One predictor variable may include missing data (coded as "NA"). The model will assume the maximum-likelihood best-fit BM model and simulate the missing predictor values throughout the course of the MCMC. These values will then be used to calculate the likelihood given the parameters for each MCMC step.
Make a color transparent (Taken from an answer on StackOverflow by Nick Sabbe)
makeTransparent(someColor, alpha = 100)makeTransparent(someColor, alpha = 100)
someColor |
A color, either a number, string or hexidecimal code |
alpha |
The alpha transparency. The maxColorValue is set to 255. |
Default bayou models. New models may be specified by providing a set of moves, control weights, tuning parameters, parameter names, RJ parameters and a likelihood function.
model.OUmodel.OU
An object of class list of length 9.
Function for calculating likelihood of an OU model in bayou using pruning algorithm or matrix inversion
OU.lik(pars, tree, X, SE = 0, model = "OU", invert = FALSE)OU.lik(pars, tree, X, SE = 0, model = "OU", invert = FALSE)
pars |
A list of parameters to calculate the likelihood |
tree |
A phylogenetic tree of class 'phylo' |
X |
A named vector giving the tip data |
SE |
A named vector or single number giving the standard errors of the data |
model |
Parameterization of the OU model. Either "OU", "QG" or "OUrepar". |
invert |
A logical indicating whether the likelihood should be solved by matrix inversion, rather than the pruning algorithm. This is primarily present to test that calculation of the likelihood is correct. |
This function can be used for calculating single likelihoods using previously
implemented methods. It is likely to become deprecated and replaced by bayou.lik
in the future, which is based on phylolm's threepoint algorithm, which works on
non-ultrametric trees and is substantially faster.
A list returning the log likelihood ("loglik"), the weight matrix ("W"), the optima ("theta"), the residuals ("resid") and the expected values ("Exp").
Calculates the alpha and sigma^2 from a parameter list with supplied phylogenetic half-life and stationary variance
OU.repar(pars)OU.repar(pars)
pars |
A bayou formatted parameter list with parameters halflife (phylogenetic halflife) and Vy (stationary variance) |
A list with values for alpha and sig2.
Experimental phenogram plotting function for set of model of model parameters
OUphenogram(pars, tree, dat, SE = 0, regime.col = NULL, ...)OUphenogram(pars, tree, dat, SE = 0, regime.col = NULL, ...)
pars |
A bayou formatted parameter list |
tree |
A tree of class 'phylo' |
dat |
A named vector of tip data |
SE |
Standard error of the tip states |
regime.col |
A named vector of colors equal in length to the number of regimes |
... |
Optional arguments passed to |
This is an experimental plotting utility that can plot a phenogram with a given regime painting from a parameter list. Note that it uses optimization of internal node states using matrix inversion, which is very slow for large trees. However, what is returned is the maximum likelihood estimate of the internal node states given the model, data and the parameter values.
## Not run: tree <- sim.bdtree(n=50) tree$edge.length <- tree$edge.length/max(branching.times(tree)) prior <- make.prior(tree, dists=list(dk="cdpois", dsig2="dnorm", dtheta="dnorm"), param=list(dk=list(lambda=5, kmax=10), dsig2=list(mean=1, sd=0.01), dtheta=list(mean=0, sd=3)), plot.prior=FALSE) pars <- priorSim(prior, tree, plot=FALSE, nsim=1)$pars[[1]] pars$alpha <- 4 dat <- dataSim(pars, model="OU", phenogram=FALSE, tree)$dat OUphenogram(pars, tree, dat, ftype="off") ## End(Not run)## Not run: tree <- sim.bdtree(n=50) tree$edge.length <- tree$edge.length/max(branching.times(tree)) prior <- make.prior(tree, dists=list(dk="cdpois", dsig2="dnorm", dtheta="dnorm"), param=list(dk=list(lambda=5, kmax=10), dsig2=list(mean=1, sd=0.01), dtheta=list(mean=0, sd=3)), plot.prior=FALSE) pars <- priorSim(prior, tree, plot=FALSE, nsim=1)$pars[[1]] pars$alpha <- 4 dat <- dataSim(pars, model="OU", phenogram=FALSE, tree)$dat OUphenogram(pars, tree, dat, ftype="off") ## End(Not run)
OUwie2bayou Converts OUwie formatted data into a bayou formatted parameter list
OUwie2bayou(tree, trait)OUwie2bayou(tree, trait)
tree |
A phylogenetic tree with states at internal nodes as node labels |
trait |
A data frame in OUwie format |
A bayou formatted parameter list
These functions calculate weight matrices from regimes specified by a bayou formatted parameter list
parmap.W calculates the weight matrix for a set of regimes from a phylogeny
with a stored regime history. .parmap.W calculates the same matrix, but without checks and is
generally run internally.
parmap.W(tree, pars)parmap.W(tree, pars)
tree |
either a tree of class "phylo" or a cache object produced by bayOU's internal functions. Must include list element 'maps' which is a simmap reconstruction of regime history. |
pars |
a list of the parameters used to calculate the weight matrix. Only pars$alpha is necessary to calculate the matrix, but others can be present. |
.parmap.W is more computationally efficient within a mcmc and is used internally.
This function converts a bayou formatted parameter list specifying regime locations into a simmap formatted tree that can
be plotted using plotSimmap from phytools or the plotRegimes function from bayou.
pars2simmap(pars, tree)pars2simmap(pars, tree)
pars |
A list that contains |
tree |
A tree of class 'phylo' |
pars2simmap takes a list of parameters and converts it to simmap format
A list with elements: tree A simmap formatted tree, pars bayou formatted parameter list, and cols A named vector of colors.
tree <- reorder(sim.bdtree(n=100), "postorder") pars <- list(k=5, sb=c(195, 196, 184, 138, 153), loc=rep(0, 5), t2=2:6) tr <- pars2simmap(pars, tree) plotRegimes(tr$tree, col=tr$col)tree <- reorder(sim.bdtree(n=100), "postorder") pars <- list(k=5, sb=c(195, 196, 184, 138, 153), loc=rep(0, 5), t2=2:6) tr <- pars2simmap(pars, tree) plotRegimes(tr$tree, col=tr$col)
Plots a phenogram and the posterior density for optima values
phenogram.density( tree, dat, burnin = 0, chain, colors = NULL, pp.cutoff = NULL, K = NULL, ... )phenogram.density( tree, dat, burnin = 0, chain, colors = NULL, pp.cutoff = NULL, K = NULL, ... )
tree |
A phylogeny of class 'phylo' |
dat |
A named vector of tip data |
burnin |
The initial proportion of the MCMC to be discarded |
chain |
A bayouMCMC object that contains the results of an MCMC chain |
colors |
An optional named vector of colors to assign to regimes, |
pp.cutoff |
The posterior probability cutoff value. Branches with posterior probabilities of having a shift above this value will have the average location of the regime shift painted onto the branches. |
K |
A list with the values of K to be plotted. If |
... |
Additional parameters passed to |
S3 method for plotting bayouMCMC objects
## S3 method for class 'bayouMCMC' plot(x, ...)## S3 method for class 'bayouMCMC' plot(x, ...)
x |
A mcmc chain of class 'bayouMCMC' produced by the function bayou.mcmc and loaded into the environment using load.bayou |
... |
Additional arguments passed to |
S3 method for plotting ssMCMC objects
## S3 method for class 'ssMCMC' plot(x, ...)## S3 method for class 'ssMCMC' plot(x, ...)
x |
An 'ssMCMC' object |
... |
Additional arguments passed to |
Produces 4 plots. The first 3 plot the prior, reference function and likelihood. Different colors indicate different power posteriors for each. These chains should appear to be well mixed. The final plot shows the sum of the marginal likelihood across each of the steps in the stepping stone algorithm.
Plot parameter list as a simmap tree
plotBayoupars(pars, tree, ...)plotBayoupars(pars, tree, ...)
pars |
A bayou formatted parameter list |
tree |
A tree of class 'phylo' |
... |
Additional arguments passed to plotRegimes |
A function to plot a heatmap of reconstructed parameter values on the branches of the tree
plotBranchHeatMap( tree, chain, variable, burnin = 0, nn = NULL, pal = heat.colors, legend_ticks = NULL, legend_settings = list(plot = TRUE), ... )plotBranchHeatMap( tree, chain, variable, burnin = 0, nn = NULL, pal = heat.colors, legend_ticks = NULL, legend_settings = list(plot = TRUE), ... )
tree |
A phylogenetic tree |
chain |
A bayou MCMC chain |
variable |
The parameter to reconstruct across the tree |
burnin |
The initial proportion of burnin samples to discard |
nn |
The number of discrete categories to divide the variable into |
pal |
A color palette function that produces nn colors |
legend_ticks |
The sequence of values to display a legend for |
legend_settings |
A list of legend attributes (passed to bayou:::.addColorBar) |
... |
Additional options passed to plot.phylo |
legend_settings is an optional list of any of the following:
legend - a logical indicating whether a legend should be plotted
x - the x location of the legend
y - the y location of the legend
height - the height of the legend
width - the width of the legend
n - the number of gradations in color to plot from the palette
adjx - an x adjustment for placing text next to the legend bar
cex.lab - the size of text labels next to the legend bar
text.col - The color of text labels
locator - if TRUE, then x and y coordinates are ignored and legend is placed interactively.
A function to visualize a multi-optimum OU process evolving on a phylogeny
plotOUtreesim(pars, tree, ptsperunit = 100, pal = rainbow, aph = 255, lwd = 1)plotOUtreesim(pars, tree, ptsperunit = 100, pal = rainbow, aph = 255, lwd = 1)
pars |
A bayou parameter list to simulate the OU process from |
tree |
A phylogenetic tree |
ptsperunit |
A number giving the number of points to simulate per unit time |
pal |
A color palette function |
aph |
The alpha value for transparency of the lines |
lwd |
The width of the lines |
Function to plot the regimes from a simmap tree
plotRegimes(tree, col = NULL, lwd = 1, pal = rainbow, ...)plotRegimes(tree, col = NULL, lwd = 1, pal = rainbow, ...)
tree |
A simmap tree of class phylo or simmap with a tree$maps list |
col |
A named vector of colors to assign to character states, if NULL, then colors are generated from pal |
lwd |
A numeric value indicating the width of the edges |
pal |
A color palette function to generate colors if col=NULL |
... |
Optional arguments that are passed to plot.phylo |
This function uses plot.phylo to generate coordinates and plot the tree, but plots the 'maps' element of phytools' simmap format. This provides much of the functionality of plot.phylo from the ape package. Currently, only types 'phylogram', 'unrooted', 'radial', and 'cladogram' are allowed. Phylogenies must have branch lengths.
shiftSummaries
A function to plot a list produced by shiftSummaries
plotShiftSummaries( summaries, pal = rainbow, ask = FALSE, single.plot = FALSE, label.pts = TRUE, ... )plotShiftSummaries( summaries, pal = rainbow, ask = FALSE, single.plot = FALSE, label.pts = TRUE, ... )
summaries |
A list produced by the function |
pal |
A color palette function |
ask |
Whether to wait for the user between plotting each shift summary |
single.plot |
A logical indicating whether to summarize all shifts in a single plot. |
label.pts |
A logical indicating whether to label the scatter plot. |
... |
Additional parameters passed to the function par(...) |
For each shift, this function plots the taxa on the phylogeny that are (usually) in this regime (each taxon is assigned to the specified shifts, thus some descendent taxa may not always be in indicated regime if the shift if they are sometimes in another tipward shift with low posterior probability). The function then plots the distribution of phenotypic states and the predicted regression line, as well as density plots for the intercept and any regression coefficients in the model.
Plot a phylogenetic tree with posterior probabilities from a bayouMCMC chain (function adapted from phytools' plotSimmap)
plotSimmap.mcmc( chain, burnin = NULL, lwd = 1, edge.type = c("regimes", "theta", "none", "pp"), pal = rainbow, pp.cutoff = 0.3, circles = TRUE, circle.cex.max = 3, circle.col = "red", circle.pch = 21, circle.lwd = 0.75, circle.alpha = 100, pp.labels = FALSE, pp.col = 1, pp.alpha = 255, pp.cex = 0.75, edge.color = 1, parameter.sample = 1000, ... )plotSimmap.mcmc( chain, burnin = NULL, lwd = 1, edge.type = c("regimes", "theta", "none", "pp"), pal = rainbow, pp.cutoff = 0.3, circles = TRUE, circle.cex.max = 3, circle.col = "red", circle.pch = 21, circle.lwd = 0.75, circle.alpha = 100, pp.labels = FALSE, pp.col = 1, pp.alpha = 255, pp.cex = 0.75, edge.color = 1, parameter.sample = 1000, ... )
chain |
A bayouMCMC chain |
burnin |
The proportion of runs to be discarded, if NULL, then the value stored in the bayouMCMC chain's attributes is used |
lwd |
The width of the edges |
edge.type |
Either "theta" (branches will be colored according to their median value of theta), "regimes" (clades will be assigned to distinct regimes if the posterior probability of a shift on that branch is > pp.cutoff), or "pp" (branches will be colored according to the probability of a shift on that branch). If "none" then edge.color will be assigned to all branches. |
pal |
A color palette function used to paint the branches (unless edge.type="none") |
pp.cutoff |
If edge.type=="regimes", the posterior probability above which a shift should be reconstructed on the tree. |
circles |
a logical value indicating whether or not a circle should be plotted at the base of the node with values that correspond to the posterior probability of having a shift. |
circle.cex.max |
The cex value of a circle with a posterior probability of 1 |
circle.col |
The color used to fill the circles |
circle.pch |
the type of symbol used to plot at the node to indicate posterior probability |
circle.lwd |
the line width of the points plotted at the nodes |
circle.alpha |
a value between 0 and 255 that indicates the transparency of the circles (255 is completely opaque). |
pp.labels |
a logical indicating whether the posterior probability for each branch should be printed above the branch |
pp.col |
The color used for the posterior probability labels |
pp.alpha |
a logical or numeric value indicating transparency of posterior probability labels. If TRUE, then transparency is ramped from invisible (pp=0), to black (pp=1). If numeric, all labels are given the same transparency. If NULL, then no transparency is given. |
pp.cex |
the size of the posterior probability labels |
edge.color |
The color of edges if edge.type="none" |
parameter.sample |
When edge.type=="theta", the number of samples used to estimate the median "theta" value from each branch. Since this is computationally intensive, this enables you to downsample the chain. |
... |
Additional arguments passed to ape's plot.phylo |
S3 method for printing bayouFit objects
## S3 method for class 'bayouFit' print(x, ...)## S3 method for class 'bayouFit' print(x, ...)
x |
A 'bayouFit' object produced by |
... |
Additional parameters passed to |
S3 method for printing bayouMCMC objects
## S3 method for class 'bayouMCMC' print(x, ...)## S3 method for class 'bayouMCMC' print(x, ...)
x |
A mcmc chain of class 'bayouMCMC' produced by the function bayou.mcmc and loaded into the environment using load.bayou |
... |
Additional arguments |
S3 method for printing priorFn objects
## S3 method for class 'priorFn' print(x, ...)## S3 method for class 'priorFn' print(x, ...)
x |
A function of class 'priorFn' produced by |
... |
Additional arguments passed to |
S3 method for printing refFn objects
## S3 method for class 'refFn' print(x, ...)## S3 method for class 'refFn' print(x, ...)
x |
A function of class 'refFn' produced by make.refFn |
... |
Additional arguments passed to |
S3 method for printing ssMCMC objects
## S3 method for class 'ssMCMC' print(x, ...)## S3 method for class 'ssMCMC' print(x, ...)
x |
An ssMCMC object |
... |
Optional arguments passed to print |
priorSim Simulates parameters from the prior distribution specified by make.prior
priorSim(prior, tree, plot = TRUE, nsim = 1, shiftpars = "theta", ...)priorSim(prior, tree, plot = TRUE, nsim = 1, shiftpars = "theta", ...)
prior |
A prior function created by |
tree |
A tree of class 'phylo' |
plot |
A logical indicating whether the simulated parameters should be plotted |
nsim |
The number of parameter sets to be simulated |
shiftpars |
A vector of parameters that split upon a shift, default is "theta" |
... |
Parameters passed on to |
A list of bayou parameter lists
Utility function for retrieving parameters from an MCMC chain
pull.pars(i, chain, model = "OU")pull.pars(i, chain, model = "OU")
i |
An integer giving the sample to retrieve |
chain |
A bayouMCMC chain |
model |
The parameterization used, either "OU", "QG" or "OUrepar" |
A bayou formatted parameter list
## Not run: tree <- sim.bdtree(n=30) tree$edge.length <- tree$edge.length/max(branching.times(tree)) prior <- make.prior(tree, dists=list(dk="cdpois", dsig2="dnorm", dtheta="dnorm"), param=list(dk=list(lambda=15, kmax=32), dsig2=list(mean=1, sd=0.01), dtheta=list(mean=0, sd=3)), plot.prior=FALSE) pars <- priorSim(prior, tree, plot=FALSE, nsim=1)$pars[[1]] dat <- dataSim(pars, model="OU", phenogram=FALSE, tree)$dat fit <- bayou.mcmc(tree, dat, model="OU", prior=prior, new.dir=TRUE, ngen=5000, plot.freq=NULL) chain <- load.bayou(fit, save.Rdata=TRUE, cleanup=TRUE) plotBayoupars(pull.pars(300, chain), tree) ## End(Not run)## Not run: tree <- sim.bdtree(n=30) tree$edge.length <- tree$edge.length/max(branching.times(tree)) prior <- make.prior(tree, dists=list(dk="cdpois", dsig2="dnorm", dtheta="dnorm"), param=list(dk=list(lambda=15, kmax=32), dsig2=list(mean=1, sd=0.01), dtheta=list(mean=0, sd=3)), plot.prior=FALSE) pars <- priorSim(prior, tree, plot=FALSE, nsim=1)$pars[[1]] dat <- dataSim(pars, model="OU", phenogram=FALSE, tree)$dat fit <- bayou.mcmc(tree, dat, model="OU", prior=prior, new.dir=TRUE, ngen=5000, plot.freq=NULL) chain <- load.bayou(fit, save.Rdata=TRUE, cleanup=TRUE) plotBayoupars(pull.pars(300, chain), tree) ## End(Not run)
Calculates the alpha parameter from a QG model
QG.alpha(pars)QG.alpha(pars)
pars |
A bayou formatted parameter list with parameters h2 (heritability), P (phenotypic variance) and w2 (width of adaptive landscape) |
An alpha value according to the equation alpha = h2*P/(P+w2+P).
Calculates the sigma^2 parameter from a QG model
QG.sig2(pars)QG.sig2(pars)
pars |
A bayou formatted parameter list with parameters h2 (heritability), P (phenotypic variance) and Ne (Effective population size) |
An sig2 value according to the equation alpha = h2*P/(Ne).
Adds visualization of regimes to a plot
regime.plot(pars, tree, cols, type = "rect", transparency = 100)regime.plot(pars, tree, cols, type = "rect", transparency = 100)
pars |
A bayou formatted parameter list |
tree |
A tree of class 'phylo' |
cols |
A vector of colors to give to regimes, in the same order as pars$sb |
type |
Either "rect", "density" or "lines". "rect" plots a rectangle for the 95% CI for the stationary distribution of a regime. "density" varies the transparency of the rectangles according to the probability density from the stationary distribution. "lines" plots lines for the mean and 95% CI's without filling them. |
transparency |
The alpha transparency value for the maximum density, max value is 255. |
Set the burnin proportion for bayouMCMC objects
set.burnin(chain, burnin = 0.3)set.burnin(chain, burnin = 0.3)
chain |
A bayouMCMC chain or an ssMCMC chain |
burnin |
The burnin proportion of samples to be discarded from downstream analyses. |
A bayouMCMC chain or ssMCMC chain with burnin proportion stored in the attributes.
A function for summarizing the state of a model after a shift
shiftSummaries(chain, mcmc, pp.cutoff = 0.3, branches = NULL)shiftSummaries(chain, mcmc, pp.cutoff = 0.3, branches = NULL)
chain |
A bayouMCMC chain |
mcmc |
A bayou mcmc object |
pp.cutoff |
The threshold posterior probability for shifts to summarize, if 'branches' specified than this is ignored. |
branches |
The specific branches with shifts to summarize, assuming postordered tree |
shiftSummaries summarizes the immediate parameter values after a shift on a particular branch. Parameters are summarized only for the duration that the particular shift exists. Thus, even global parameters will be different for particular shifts.
A list with elements:
pars = a bayoupars list giving the location of shifts specified;
tree = The tree;
pred = Predictor variable matrix;
dat = A vector of the data;
SE = A vector of standard errors;
PP = Posterior probabilities of the specified shifts;
model = A list specifying the model used;
variables = The variables summarized;
cladesummaries = A list providing the medians and densities of the distributions of regression
variables for each shift;
descendents = A list providing the taxa that belong to each regime
regressions = A matrix providing the regression coefficients for each regime.
These functions calculate weight matrices from regimes specified in phytools' simmap format.
simmapW calculates the weight matrix for a set of regimes from a phylogeny
with a stored regime history. .simmap.W calculates the same matrix, but without checks and is
generally run internally.
simmapW(tree, pars)simmapW(tree, pars)
tree |
either a tree of class "phylo" or a cache object produced by bayOU's internal functions. Must include list element 'maps' which is a simmap reconstruction of regime history. |
pars |
a list of the parameters used to calculate the weight matrix. Only pars$alpha is necessary to calculate the matrix, but others can be present. |
.simmap.W is more computationally efficient within a mcmc and is used internally. The value
of TotExp is supplied to speed computation and reduce redundancy, and cache objects must be supplied as
the phylogeny, and the parameter ntheta must be present in the list pars.
S3 method for summarizing bayouMCMC objects
## S3 method for class 'bayouMCMC' summary(object, ...)## S3 method for class 'bayouMCMC' summary(object, ...)
object |
A bayouMCMC object |
... |
Additional arguments passed to |
An invisible list with two elements: statistics which provides
summary statistics for a bayouMCMC chain, and branch.posteriors which summarizes
branch specific data from a bayouMCMC chain.