Rev Language Reference

mvSpeciesNarrow - Narrow-exchange joint move on species tree and gene trees for multispecies coalescent models.

Makes a narrow-exchange move both in the species tree and in the gene trees that contain nodes of the relevant populations.

Usage

mvSpeciesNarrow(TimeTree speciesTree, RealPos weight, Probability tuneTarget)

Arguments

speciesTree : TimeTree (<stochastic> pass by reference)
The species tree variable on which this move operates.
weight : RealPos (pass by value)
The weight determines the relative frequency with which this move will be attempted. For details, see the description of the 'moveschedule' parameter on the documentation page for 'mcmc()'.
Default : 1
tuneTarget : Probability (pass by value)
The acceptance probability targeted by auto-tuning.
Default : 0.44

Details

The species tree must be ultrametric. All the gene trees that evolved along the species tree according to some form of multispecies coalescent must be added to the move using the addGeneTreeVariable method. This move jointly performs narrow exchange moves (Nearest-Neighbor Interchanges without branch length alterations) on the species tree and on gene trees, all of which must be ultrametric.

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 a constant effective population size of 50:
popSize <- 50

# 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] ~ dnMultiSpeciesCoalescent(speciesTree=spTree, Ne=popSize, taxa=taxa[g])
    print(geneTrees[g])
}

# Set my move index:
mi = 0
move_species_narrow_exchange = mvSpeciesNarrow( speciesTree=spTree, weight=5 )
for (i in 1:n_genes) {
    move_species_narrow_exchange.addGeneTreeVariable( geneTrees[i] )
}
moves[++mi] = move_species_narrow_exchange

# We get a handle on our model.
# We can use any node of our model as a handle; here we choose to use the topology.
mymodel = model(spTree)

# Monitors to check the progression of the program:
monitors[1] = mnScreen(printgen=10, spTree)

# Here we use a plain MCMC. You could also use mcmcmc with heated chains.
mymcmc = mcmc(mymodel, monitors, moves, nruns=4)
mymcmc.run(generations=1000)
mymcmc.operatorSummary()

Methods

  • addGeneTreeVariable(TimeTree geneTree)

See Also

Edit this page