Search code examples
rphylogenyape-phylo

Phylogenetic ancestor reconstruction: specify model of unidirectional character state change [ape][phytools]


Let's say I have a phylogenetic tree and some character data for my tree. I have one character that I know is unidirectional: the rate of transition from 0 to 1 is positive, but the rate of transition from 1 to 0 is zero (for instance, 0 is diploid and 1 is polyploid). Let's say I want to do an ancestor reconstruction of my tree. I know I can use Ace from package Ape, or make.simmap from package phytools to do ancestor reconstruction, but I don't know how to specify the model of unidirectional character change.

Example:

require(ape)
require(phytools)

# Generate example tree
set.seed(100)
tree<-rtree(10, rooted=T, tip.label=letters[1:10])
# Example characters
characters<-sample(0:1, 10, replace=T)
names(characters)<-letters[1:10]

# Equal-rates model
mod=matrix(c(0,1,1,0),2)
simmap<-make.simmap(tree, characters, model=mod)
plotSimmap(simmap) # Works fine


# My attempt at a unidirectional model: rate of change from 1 to 0 set to zero
mod = matrix(c(0,0,1,0),2)
simmap<-make.simmap(tree, characters, model=mod) # Gives me error; Error in eigen(mat) : infinite or missing values in 'x'

Anyone have any idea what to do?


Solution

  • Your code above worked fine for me using R version 3.1.2 (2014-10-31) -- "Pumpkin Helmet" under Windows with newly installed versions of ape and phytools. The make.simmap help documentation talks about some errors in previous versions.

    mod = matrix(c(0,0,1,0),2)
    simmap<-make.simmap(tree, characters, model=mod)
    make.simmap is sampling character histories conditioned on the transition matrix
    #Q =
    #           0         1
    #0 -0.8271365 0.8271365
    #1  0.0000000 0.0000000
    #(estimated using likelihood);
    #and (mean) root node prior probabilities
    #pi =
    #  0   1 
    #0.5 0.5 
    #Done.