make.prior: Make a prior function for bayou

In uyedaj/bayou: Bayesian Fitting of Ornstein-Uhlenbeck Models to Phylogenies

Make a prior function for bayou

Description

This function generates a prior function to be used for bayou according to user specifications.

Usage

make.prior(
  tree,
  dists = list(),
  param = list(),
  fixed = list(),
  plot.prior = TRUE,
  model = "OU"
)

Arguments

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 "OU" (alpha & sig2), "QG" (h2, P, w2, Ne), or "OUrepar" (halflife,Vy)

Details

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.

Value

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.

Examples

## 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

uyedaj/bayou documentation built on Jan. 28, 2024, 5:09 a.m.