contMap with user-supplied node states

A phytools user contacted me earlier today about generating a contMap style plot for data in which internal node states were estimated using a different model, such as an OU model.

This opens up the general question of whether or not we can supply input states for internal nodes to the function. In fact, this is not currently permitted - however I did previously write a post describing a hack in which this was done.

OK, well now I have just posted an update to contMap that introduces a new method (method="user") in which the user can supply ancestral values for some or all of the internal nodes of the tree. The states along the edges are interpolated using the expression of Felsenstein - so these are not guaranteed to be the ancestral states along the edges that one would obtain by reconstructing under the model via which internal node states have been obtained (as they are for the BM model)!

Here is a quick demo:

library(devtools)
install_github("liamrevell/phytools")

library(phytools)

## Loading required package: ape
## Loading required package: maps
## 
##  # ATTENTION: maps v3.0 has an updated 'world' map.        #
##  # Many country borders and names have changed since 1990. #
##  # Type '?world' or 'news(package="maps")'. See README_v3. #

tree<-pbtree(n=26,scale=10,tip.label=LETTERS)
## simulate Brownian motion data with a trend
x<-fastBM(tree,sig2=0.1,internal=TRUE,mu=0.4)
## ancestral states
a<-x[1:tree$Nnode+Ntip(tree)]
a

##        27        28        29        30        31        32        33 
## 0.0000000 0.5423934 2.8023807 2.7520736 4.7089698 3.3488780 3.4507386 
##        34        35        36        37        38        39        40 
## 3.7591624 3.6837011 1.7904929 1.8667323 3.1817058 3.1827339 2.3905485 
##        41        42        43        44        45        46        47 
## 3.1213984 4.3085840 5.3381130 3.4971804 3.9120370 4.1876108 4.1280892 
##        48        49        50        51 
## 2.8458030 5.3843788 3.7128553 4.0341649

x<-x[1:Ntip(tree)] ## tip states
x

##        A        B        C        D        E        F        G        H 
## 4.561221 4.807578 4.882611 4.097495 3.585475 3.685503 4.165546 4.490254 
##        I        J        K        L        M        N        O        P 
## 4.911537 3.123553 3.299806 3.162068 5.551774 5.572570 5.760226 4.119402 
##        Q        R        S        T        U        V        W        X 
## 4.143924 4.168398 4.483631 3.657520 5.075105 5.692124 5.472482 4.478056 
##        Y        Z 
## 4.085486 4.756476

## some ancestral states
anc.states<-sample(a,8)
anc.states

##       35       41       27       30       34       45       49       36 
## 3.683701 3.121398 0.000000 2.752074 3.759162 3.912037 5.384379 1.790493

Now we're ready to use contMap. First, without any internal node states:

fit.none<-contMap(tree,x,plot=FALSE)
fit.none

## Object of class "contMap" containing:
## 
## (1) A phylogenetic tree with 26 tips and 25 internal nodes.
## 
## (2) A mapped continuous trait on the range (3.123553, 5.760226).

plot(fit.none)

Now with all states:

fit.all<-contMap(tree,x,plot=FALSE,method="user",anc.states=a)
fit.all

## Object of class "contMap" containing:
## 
## (1) A phylogenetic tree with 26 tips and 25 internal nodes.
## 
## (2) A mapped continuous trait on the range (0, 5.760226).

plot(fit.all)

Finally, with just a subset:

fit.some<-contMap(tree,x,plot=FALSE,method="user",anc.states=anc.states)
fit.some

## Object of class "contMap" containing:
## 
## (1) A phylogenetic tree with 26 tips and 25 internal nodes.
## 
## (2) A mapped continuous trait on the range (0, 5.760226).

plot(fit.some)

That's it. Just did this, so please report any issues that might be related to this update. Thanks!