Discrete morphology - Multistate characters

Introduction

So far in this tutorial series we have only worked with binary characters. Real morphological matrices very often mix characters with different numbers of observed states in the same matrix: some characters are simple presence/absence codings, while others (like “shape of element X: round / oval / triangular”) have three, four, or more states. This tutorial shows how to combine the Mk model and ACRV (Discrete morphology - Among-character rate variation) with a dataset like this.

We will use a dataset of 27 extinct echinoderm (cinctan) taxa and 60 binary and multistate characters, from Zamora et al. (2013).

The Key Idea: Partitioning by Number of States

A single Q matrix has a fixed number of rows and columns, matching the number of states of the characters it describes. This means we cannot use one Q matrix for a matrix that mixes 2-state and 4-state characters. Instead, RevBayes lets us split (“partition”) the data matrix into blocks that each contain only characters with the same number of observed states, and give each block its own, appropriately-sized Q matrix (Khakurel et al. 2024). All of the blocks still share the same tree and, if we want, the same ACRV hyperparameters — only the dimension of the Q matrix changes from block to block.

Example: Inferring a Phylogeny of Extinct Cinctans

Data and Files

Download the data file called Cinctans.nex and save it to your data directory.

Creating Rev Files

Create the Rev-script files model_multistate_partition.Rev and mcmc_multistate.Rev in your scripts directory.

We continue to reuse model_tree.Rev, unchanged, from Discrete morphology - The Mk model.

Load Data Matrices

Loading the data matrix works exactly as before.

morpho <- readDiscreteCharacterData("data/Cinctans.nex")

taxa <- morpho.names()
num_taxa <- taxa.size()
num_branches <- 2 * num_taxa - 2

moves    = VectorMoves()
monitors = VectorMonitors()

source("scripts/model_tree.Rev")

Partitioning by Number of States

In model_multistate_partition.Rev, we again start with Gamma-distributed ACRV, exactly as in Discrete morphology - Among-character rate variation.

alpha_morpho ~ dnUniform( 0.0, 1E6 )
rates_morpho := fnDiscretizeGamma( alpha_morpho, alpha_morpho, 4 )
moves.append( mvScale(alpha_morpho, lambda=1, weight=2.0) )

We do not know in advance how many states the most variable character in our matrix has, so we pick a generous upper bound (n_max_states) and loop over every possible number of states from 2 up to that bound. For each candidate state count i, we make a local copy of the data, restrict it to only the characters that actually have i states using setNumStatesPartition(), and skip the block entirely if it turns out to be empty.

n_max_states <- 7
idx = 1
morpho_bystate[1] <- morpho
for (i in 2:n_max_states) {

morpho_bystate[i] <- morpho
    morpho_bystate[i].setNumStatesPartition(i)

nc = morpho_bystate[i].nchar()

For every non-empty block, we build an i-by-i Jukes-Cantor matrix with fnJC(i), and give it its own dnPhyloCTMC, sharing the same tree and the same discretized rates across all blocks:

if (nc > 0) {

        # make an i-by-i rate matrix
        q[idx] <- fnJC(i)

        # create the model of evolution for the character block
        m_morph[idx] ~ dnPhyloCTMC( tree=phylogeny,
                                    Q=q[idx],
                                    nSites=nc,
                                    siteRates=rates_morpho,
                                    type="Standard")

        # attach the data
        m_morph[idx].clamp(morpho_bystate[i])

        # increment counter
        idx = idx + 1
        idx
    }
}

It would not make sense to model a 5-state character using a 2-state Q matrix (or vice versa), which is exactly what this partitioning scheme avoids.

Back in mcmc_multistate.Rev, source this file, build the model, and run the MCMC exactly as before.

source("scripts/model_multistate_partition.Rev")

mymodel = model(phylogeny)

monitors.append( mnModel(filename="output/multistate.log", printgen=10) )
monitors.append( mnFile(filename="output/multistate.trees", printgen=10, phylogeny) )
monitors.append( mnScreen(printgen=100) )

mymcmc = mcmc(mymodel, monitors, moves, nruns=2, combine="mixed")
mymcmc.run(generations=50000, tuningInterval=200)

Once the analysis is complete, we can compute several different types of summary trees:

trace = readTreeTrace("output/multistate.trees", treetype="non-clock")
trace.setBurnin(0.25)

# maximum clade credibility tree
mccTree(trace, file="output/multistate.mcc.tre")
# maximum a posteriori tree
mapTree(trace, file="output/multistate.map.tre")
# majority rule consensus tree
consensusTree(trace, file="output/multistate.majrule.tre")

Choose Your Own Adventure

For this next portion of the tutorial, we will edit some pieces of the model. Feel free to add your own innovations. Make copies of the relevant scripts and modify them; make sure to change the output filenames accordingly before running the MCMC!

Ascertainment Bias

As discussed in Discrete morphology - The Mk model, we should correct for ascertainment bias whenever our matrix has no invariant characters — as is typical of morphological data assembled under parsimony. Because our substitution model is now split across several dnPhyloCTMC calls (one per state-count partition), the fix is to add coding="variable" to each of them. Try adding it to model_multistate_partition.Rev yourself, and compare the resulting tree-length estimates to your original, uncorrected run — you should see the same inflation-without-correction pattern we saw for the binary bears dataset in Discrete morphology - The Mk model.

A Lognormal Alternative

As in Discrete morphology - Among-character rate variation, you can also try substituting a discretized lognormal distribution for the Gamma distribution in model_multistate_partition.Rev:

rates_morpho := fnDiscretizeDistribution( dnLognormal(ln(alpha_morpho), 0.01), 4 )

Have each person at your table try a different one of these variants, load all of the resulting log files into Tracer, and see how they differ.

When you are ready, move on to Discrete morphology - Ordered characters, where we use the same partitioning scheme for a matrix in which some characters are ordered, so that they can only change between adjacent states.