MSCsimtester (version 1.1)

pairwiseDist: Compute and plot sample and theoretical pairwise distance densities.

Description

Computes theoretical pairwise distance densities under the MSC on a species tree and empirical pairwise distances from gene trees in a sample. A histogram of empirical values is plotted over the theoretical pdf.

Usage

pairwiseDist(
  stree,
  popSizes,
  gtSample,
  taxon1,
  taxon2,
  numSteps = 1000,
  tailProb = 0.01,
  numBreaks = 40
)

Value

A list of items needed for Anderson-Darling test(s), for use by ADtest, returned invisibly. See function code for more details.

Arguments

stree

An object of class phylo containing a rooted metric species tree. Edge lengths are in number of generations.

popSizes

A vector containing constant population sizes, one entry for each edge/population in the species tree, for a haploid population. Sizes should be doubled for diploids. If stree has k edges, then popSizes must have k+1 elements, with final entry the size of the population ancestral to the root.

gtSample

An object of class multiPhylo holding a sample of gene trees from a simulation. Taxon labels on gene trees must be identical to those on stree.

taxon1

A string specifying one taxon on stree.

taxon2

A string specifying a second taxon on stree, distinct from taxon1.

numSteps

A positive integer number of values to be computed for graphing the theoretical pairwise distance density. Default is numSteps = 1000. A larger value produces a smoother plot.

tailProb

A cutoff value, between 0 and 1, for the theoretical density, with a default of 0.01.

numBreaks

Number of breaks in histogram. Default is numBreaks = 40. The theoretical pairwise distance is plotted from (0, xMax), where xMax is the larger of the maximum pairwise distance in the gene tree sample and the value cutting off a tail of area tailProb under the pdf. A message returns the proportion of sample distances in this tail.

See Also

plotEdgeOrder, plotPops, ADtest

Examples

stree=read.tree(text="((((a:10000,b:10000):10000,c:20000):10000,d:30000):10000,e:40000);")
pops=c(15000,25000,10000,1,1,1,1,1,12000)
gts=read.tree(file=system.file("extdata","genetreeSample",package="MSCsimtester"))
pairwiseDist(stree,pops,gts,"a","b")