Estimate Gamma Distribution for MCMCtree analysis

Estimate the shape and rate paramaters of Gamma distribution and output trees for MCMCtree input

Usage

estimateGamma(
  minAge,
  maxAge,
  phy,
  monoGroups,
  alpha = 188,
  beta = 2690,
  offset = 0.1,
  estimateAlpha = TRUE,
  estimateBeta = FALSE,
  plot = FALSE,
  pdfOutput = "gammaPlot.pdf",
  writeMCMCtree = FALSE,
  MCMCtreeName = "gammaInput.tre"
)

Arguments

minAge: vector of minimum age bounds for nodes matching order in monoGroups
maxAge: vector of maximum age bounds for nodes matching order in monoGroups
phy: fully resolved phylogeny in ape format
monoGroups: list with each element containing species that define a node of interest
alpha: alpha value for gamma distribution (default = 188)
beta: beta value for gamma distribution (default = 2690)
offset: distance of mean value from minimum bound
estimateAlpha: logical specifying whether to estimate alpha with a given beta value (default = TRUE)
estimateBeta: logical specifying whether to estimate beta with a given alpha value (default = FALSE)
plot: logical specifying whether to plot to PDF
pdfOutput: pdf output file name
writeMCMCtree: logical whether to write tree in format that is compatible with MCMCTree to file
MCMCtreeName: MCMCtree.output file name

Value

list containing node estimates for each distribution

"parameters" estimated parameters for each node
"apePhy" phylogeny in APE format with node labels showing node distributions
"MCMCtree" phylogeny in MCMCtreeR format
"nodeLabels" node labels in MCMCtreeR format

If plot=TRUE plot of distributions in file 'pdfOutput' written to current working directory

If writeMCMCtree=TRUE tree in MCMCtree format in file "MCMCtreeName" written to current working directory

Author

Mark Puttick

Examples

data(apeData)
attach(apeData)
#> The following objects are masked from apeData (pos = 3):
#> 
#>     apeTree, maximumTimes, minimumTimes, monophyleticGroups
#> The following objects are masked from apeData (pos = 4):
#> 
#>     apeTree, maximumTimes, minimumTimes, monophyleticGroups
#> The following objects are masked from apeData (pos = 5):
#> 
#>     apeTree, maximumTimes, minimumTimes, monophyleticGroups
## extract taxon descending from calibrated nodes 8, 10, 11, 13
## these nodes can be visualised using plot.phylo
## and nodelabels from ape
monophyleticGroups <- tipDes(apeData$apeTree, c(8,10,11,13))
minimumTimes <- c("8"=15, "10"=6,
"11"=8, "13"=13) / 10
maximumTimes <- c("8" = 30, "10" = 12,
"11"=12, "13" = 20) / 10
gamma.nodes <- estimateGamma(minAge=minimumTimes, maxAge=maximumTimes, 
monoGroups=monophyleticGroups, alpha=188, beta=2690, 
offset=0.1, phy=apeTree, plot=FALSE)
#> [1] "warning - alpha parameter value recycled"
#> [1] "warning - beta parameter value recycled"
#> [1] "warning - offset parameter value recycled"
#> [1] "warning - estimateAlpha argument recycled"
#> [1] "warning - estimateBeta argument recycled"
#> Warning: Length of node.label does not match number of nodes.
gamma.nodes
#> $parameters
#>        alpha beta
#> node_1  4304 2690
#> node_2  1883 2690
#> node_3  2421 2690
#> node_4  3766 2690
#> 
#> $apePhy
#> 
#> Phylogenetic tree with 7 tips and 6 internal nodes.
#> 
#> Tip labels:
#>   human, chimpanzee, bonobo, gorilla, orangutan, sumatran, ...
#> Node labels:
#>   'G[4304~2690]', NA, 'G[1883~2690]', 'G[2421~2690]', NA, 'G[3766~2690]'
#> 
#> Rooted; no branch lengths.
#> 
#> $MCMCtree
#>      
#> 1 7 1
#>                                                                                                                                 
#> 1 ((((human,(chimpanzee,bonobo))'G(2421,2690)',gorilla)'G(1883,2690)',(orangutan,sumatran)'G(3766,2690)'),gibbon)'G(4304,2690)';
#>                
#> 1 //end of file
#> 
#> $nodeLabels
#> [1] "'G[4304~2690]'" "'G[1883~2690]'" "'G[2421~2690]'" "'G[3766~2690]'"
#>