optimalCollapse: Determine which clades can be collapsed

In eliotmiller/shiftPlot: Phylogenetic shift plotting

Determine which clades can be collapsed

Description

Optimal collapse

Usage

optimalCollapse(orig.tree, states.df, flip.tips)

Arguments

orig.tree	Phylogeny in ape format, corresponding to states.df.
states.df	Data frame in the specified shiftPlot format. Should contain one column named "state", and number for every node in the phylogeny, with the tips above the internal nodes, and no row names. See details and examples.
flip.tips	Whether or not to flip tips to the state of their parent nodes.

Details

states.df should have one column titled "present". This column should be coded either as a 0 or a 1, indicating the presence of the trait. states.df should have as many rows as there are nodes in phylogeny, and the tip nodes should come first in the data frame. For example, you might rbind the $tip.states and $states objects from a corHMM output together to create states.df The main reason one would use this function is to preclude shifts on terminal branches. You might want to do that if, for example, your phylogeny was so large that corHMM was having trouble calculating the marginal probabilities of a tip having a hidden trait, and the result would otherwise be that it appeared there were many shifts in the hidden trait towards the tips.

Value

An updated states.df with the tip values flipped as needed.

Author(s)

Eliot Miller

Examples

#start with a corHMM output and build up a states.df
#load data. these are the results of following the corHMM precursor model vignette
data(Precur_res.corHMM)
data(phy)
nodeStates <- data.frame(state=Precur_res.corHMM$states[,3]+Precur_res.corHMM$states[,4])
tipStates <- data.frame(state=Precur_res.corHMM$tip.states[,3]+Precur_res.corHMM$tip.states[,4])

#note that tip states comes first here!
states <- rbind(tipStates, nodeStates)

#binarize this. choosing to call 0.5 chance of having trait present
states$state[states$state >= 0.5] <- 1
states$state[states$state < 0.5] <- 0

#flip node 103 and all nodes towards tips from there to trait = absent
induced <- states
induced[geiger:::.get.descendants.of.node(103, phy),] <- 0
induced[103,] <- 0

#get rid of row names
row.names(induced) <- NULL

#run the function and don't flip those tips
ocResult <- optimalCollapse(orig.tree=phy, states.df=induced, flip.tips=FALSE)

eliotmiller/shiftPlot documentation built on March 30, 2023, 4:26 a.m.