1

I'm simulating discrete character data using the function rTraitDisc {ape} in R using a variety of model matrices. I've not encountered any issues with scaling when all state changes are possible. However when I supply an ordered model with 8 or more possible states, the function breaks down and returns the following error:

## library
library(ape)

## read in tree
data("bird.orders")

## build model
model.matrix <- matrix(c(0,0.1,0,0,0,0,0,0,
                     0.1,0,0.1,0,0,0,0,0,
                     0,0.1,0,0.1,0,0,0,0,
                     0,0,0.1,0,0.1,0,0,0,
                     0,0,0,0.1,0,0.1,0,0,
                     0,0,0,0,0.1,0,0.1,0,
                     0,0,0,0,0,0.1,0,0.1,
                     0,0,0,0,0,0,0.1,0), 8)

## run function
rTraitDisc(phy = bird.orders, model = model.matrix)

Error message: Error in sample.int(k, size = 1, FALSE, prob = p) : negative probability

Having dug a little deeper, it seems that when there are 8 or more states but only one possible transition (e.g. if the ancestral state is 0, only a transition to state 1 should be possible in an ordered matrix), the function matexpo produces a probability matrix with negative values for the shortest branch of the tree (0.5). As these probabilities are used by sample.int as the "prob" argument, the negative probabilities cause the function to break down.

## get number of states
k <- ncol(model.matrix)

## get equilibrium relative frequencies
freq = rep(1/k, k)

## match number of elements in model
freq <- rep(freq, each = k)

## get Q matrix
Q <- model.matrix * freq
diag(Q) <- 0
diag(Q) <- -rowSums(Q)

## get minimum edge length
min.el <- min(bird.orders$edge.length)

## run matexpo
matexpo(Q*min.el)

How do I deal with these negative values in this context? Is there a correction I can/should apply?

SmithOfToms
  • 11
  • 2
  • If the negative probabilities are something like `-1.2e-14` they can perhaps be regarded as round-off error and replaced by `0`. On the other hand, if they are genuinely negative numbers in a context where probabilities are expected then there seems to be a fundamental problem with your approach and perhaps you need to rethink your algorithm. – John Coleman Nov 23 '21 at 13:05
  • Thanks for the quick reply. The negative values are tiny and are being returned where I would expect a probability of 0 (impossible state changes), so I think it might be a rounding error. – SmithOfToms Nov 23 '21 at 13:21

0 Answers0