MkV model
Setting variable = TRUE ensures all simulated characters
are variable across taxa, matching the MkV model of Lewis (2001).
mkv_data <- sim.morpho(time.tree = tree, k = 2, trait.num = 10,
br.rates = 0.1, variable = TRUE)
Among-character rate variation (ACRV)
Not all morphological characters evolve at the same rate. MorphSim
can model this using among-character rate variation (ACRV), where each
trait is assigned a rate drawn from a discretized distribution. Two
distributions are available.
Gamma distribution — controlled by two parameters:
alpha.gamma, the shape parameter of the gamma distribution
(smaller values produce greater rate heterogeneity, with more characters
evolving very slowly or very quickly; larger values concentrate rates
closer to the mean), and ACRV.ncats, the number of discrete
rate categories used to approximate the continuous distribution
(commonly set to 4, following standard practice in phylogenetic
inference).
gamma_data <- sim.morpho(
time.tree = tree, k = 2, trait.num = 10, br.rates = 0.1,
ACRV = "gamma", alpha.gamma = 1, ACRV.ncats = 4
)
# The discrete rate categories
gamma_data$model$RateVar
#> [1] 0.1369538 0.4767519 1.0000000 2.3862944
# Which category was assigned to each trait
gamma_data$model$RateVarTrait
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
#> [1,] 4 1 2 4 3 4 3 4 3 2
Lognormal distribution — controlled by
meanlog and sdlog, the mean and standard
deviation on the log scale, and ACRV.ncats as above. The
lognormal can produce heavier tails than the gamma, meaning a greater
spread between the slowest and fastest evolving characters.
lgn_data <- sim.morpho(
time.tree = tree, k = 2, trait.num = 10, br.rates = 0.1,
ACRV = "lgn", meanlog = 0, sdlog = 1, ACRV.ncats = 4
)
# The discrete rate categories
lgn_data$model$RateVar
#> [1] 0.2269731 0.5214124 0.9861614 2.2654531
User-defined rates — you can also provide your own
rate categories directly using ACRV = "user" and passing a
vector of rates to define.ACRV.rates.
user_data <- sim.morpho(
time.tree = tree, k = 2, trait.num = 10, br.rates = 0.1,
ACRV = "user", define.ACRV.rates = c(0.5, 1, 1.5, 2),
ACRV.ncats = 4
)
MkV + Gamma
These can be combined to simulate under an MkV + Gamma model,
consistent with models implemented in RevBayes and BEAST.
mkvg_data <- sim.morpho(
time.tree = tree, k = 3, trait.num = 10, br.rates = 0.1,
variable = TRUE, ACRV = "gamma", alpha.gamma = 1, ACRV.ncats = 4
)
Custom Q matrices
You can provide your own transition rate matrix. For example, an
ordered model where transitions can only occur between adjacent states
(0 ↔︎ 1 ↔︎ 2):
ord_Q <- matrix(c(
-0.5, 0.5, 0.0,
0.3333333, -0.6666667, 0.3333333,
0.0, 0.5, -0.5
), nrow = 3, byrow = TRUE)
ordered_data <- sim.morpho(
time.tree = tree, k = 3, trait.num = 10, br.rates = 0.1,
define.Q = ord_Q
)
Note that when using a custom Q matrix, all traits in that simulation
share the same matrix. If you want to simulate different partitions
under different Q matrices, for example, some ordered and some unordered
characters, you can simulate them separately and combine the results.
See the section on combining morpho objects below.
Specifying the root state
By default the root state is sampled uniformly. You can fix it using
root.state:
fixed_root <- sim.morpho(
time.tree = tree, k = 3, trait.num = 10,
br.rates = 0.1, root.state = 0
)
Ensure full state space
Setting full.states = TRUE ensures that all character
states specified by k are present in at least one tip for
each trait. This is useful when you want to guarantee that the full
range of states specified in the Q matrix is observed at the tips.
Without this, it is possible for a character with k = 3 to
only show states 0 and 1 at the tips if state 2 was never reached or was
lost before the present.
strict_data <- sim.morpho(
time.tree = tree, k = 3, trait.num = 10,
br.rates = 0.1, full.states = TRUE
)
Partitions
Different groups of characters can be simulated with different
numbers of states. The partition argument specifies how
many traits belong to each partition, and k gives the
number of states per partition.
part_data <- sim.morpho(
time.tree = tree,
k = c(2, 3, 4),
trait.num = 20,
partition = c(10, 5, 5),
br.rates = 0.1
)
# Model specification per partition
part_data$model$Specified
#> [1] "Mk_Part:10states:2" "Mk_Part:5states:3" "Mk_Part:5states:4"
The total number of traits defined by trait.num and partition must
match or the function will return an error.
Combining morpho objects
It is not possible to simulate multiple partitions under different
models in a single call to sim.morpho. For example, you
might want some characters to evolve under an ordered model with a
custom Q matrix and others under a standard Mk model with gamma rate
variation. To do this, simulate each partition separately and combine
them using combine.morpho.
# Partition 1: ordered characters with a custom Q matrix
ord_Q <- matrix(c(
-0.5, 0.5, 0.0,
0.3333333, -0.6666667, 0.3333333,
0.0, 0.5, -0.5
), nrow = 3, byrow = TRUE)
partition_1 <- sim.morpho(
time.tree = tree, k = 3, trait.num = 10,
br.rates = 0.1, define.Q = ord_Q
)
# Partition 2: unordered binary characters with MkV + Gamma
partition_2 <- sim.morpho(
time.tree = tree, k = 2, trait.num = 10,
br.rates = 0.1, variable = TRUE,
ACRV = "gamma", alpha.gamma = 1, ACRV.ncats = 4
)
# Combine into a single morpho object
combined <- combine.morpho(partition_1, partition_2)
# The combined object has 20 traits total
length(combined$sequences$tips[[1]])
#> [1] 20
# Model specifications from both partitions are preserved
combined$model$Specified
#> [1] "Mk_Part:10states:3" "MkgammaV_Part:10states:2"
Both morpho objects must share the same tree, branch lengths, and
fossil object (if present). The combined object merges the sequences,
transition histories, root states, and model information from both
inputs.
