Rev Language Reference

dnMultiSpeciesCoalescentInverseGamma - Multispecies coalescent distribution with inverse gamma prior on effective population sizes

Multispecies coalescent distribution describing how gene trees can be generated from within a species tree given effective population sizes. Requires an ultrametric species tree, parameters of an inverse gamma prior on effective population sizes, and taxa with species and individual names.

Usage

dnMultiSpeciesCoalescentInverseGamma(TimeTree speciesTree, RealPos shape, RealPos scale, Natural num_genes, Taxon[][] taxa)

Arguments

speciesTree : TimeTree (pass by const reference)
The species tree in which the gene trees evolve.
shape : RealPos (pass by const reference)
The shape of the inverse gamma prior distribution on the effective population sizes.
scale : RealPos (pass by const reference)
The scale of the inverse gamma prior distribution on the effective population sizes.
num_genes : Natural (pass by const reference)
The number of genes.
taxa : Taxon[][] (pass by value)
The vector of taxa which have species and individual names.

Domain Type

TimeTree[]

Details

The species tree must be ultrametric. Parameters of an inverse gamma prior on effective population sizes must be provided. This distribution uses a conjugate prior on effective population sizes. As a consequence, effective population sizes are analytically integrated out and treated as nuisance parameters (Jones 2016). If you are interested in inferring ancestral effective population sizes, use `dnMultiSpeciesCoalescent`.

Example

# We are going to save the trees we simulate in the folder simulatedTrees:
dataFolder = "simulatedTrees/"

# Let's simulate a species tree with 10 taxa, 2 gene trees, 3 alleles per species:
n_species <- 10
n_genes <- 2
n_alleles <- 3

# We simulate an ultrametric species tree.
# Species names:
for (i in 1:n_species) {
    species[i] <- taxon(taxonName="Species_"+i, speciesName="Species_"+i)
}
spTree ~ dnBirthDeath(lambda=0.3, mu=0.2, rootAge=10, rho=1, samplingStrategy="uniform", condition="nTaxa", taxa=species)
print(spTree)

# Let's pick constant parameters for the inverse gamma distribution:
alpha <- 3
beta <- 0.003
# Let's simulate gene trees now.
# Taxon names:
for (g in 1:n_genes) {
    for (i in 1:n_species) {
        for (j in 1:n_alleles) {
            taxa[g][(i-1)*n_alleles+j] <- taxon(taxonName="Species_"+i+"_"+j, speciesName="Species_"+i)
        }
    }
    geneTrees[g] ~ dnMultiSpeciesCoalescentInverseGamma(speciesTree=spTree, shape=alpha, scale=beta, taxa=taxa[g])
    print(geneTrees[g])
}

# We can save the species tree:
write(spTree, filename=dataFolder+"speciesTree")

# Saving the gene trees:
for (i in 1:(n_genes)) {
    write(geneTrees[i], filename=dataFolder+"geneTree_"+i+".tree")
}

See Also

Edit this page