R: Compute and plot sample and theoretical pairwise distance...

pairwiseDist {MSCsimtester}

R Documentation

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
)

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

Value

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

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")