Ancestral state estimation for continuous characters with some nodes fixed

I just added the method to fix some ancestral states during ML ancestral state estimation for continuous characters that I used in contMap and previously described in the pages of this blog to the phytools function fastAnc. Interested readers can look at the update here, and it can also be obtained by updating phytools from GitHub.

Here's a quick example for simulated data in which the internal states are known:

library(phytools)
## simulate a tree
tree<-pbtree(n=26,tip.label=LETTERS)
plotTree(tree)

## simulate data for tips & nodes
xa<-fastBM(tree,internal=TRUE)
a<-xa[as.character(1:tree$Nnode+Ntip(tree))]
a ## simulated ancestral states

##          27          28          29          30          31          32 
##  0.00000000 -0.56601239 -0.14105673 -1.19661359 -1.66581973 -0.93686310 
##          33          34          35          36          37          38 
## -0.97710948 -2.14034585 -0.53578250  0.09301371 -0.46628240 -0.94499655 
##          39          40          41          42          43          44 
## -1.56438497  0.24762819 -2.29351327 -1.04066230 -1.40266486 -2.00346322 
##          45          46          47          48          49          50 
##  0.08526085  0.28998116  1.42932620  1.77921782  1.58183599  1.31419370 
##          51 
## -1.16905491

x<-xa[tree$tip.label] ## simulated tip data
x

##          A          B          C          D          E          F 
##  0.5683050 -1.5702611 -0.6006650 -0.8855694 -2.3290004 -2.2543046 
##          G          H          I          J          K          L 
##  0.2865251 -0.7192771 -1.2468277 -1.5512149  0.1103867  0.5530111 
##          M          N          O          P          Q          R 
## -2.4353686 -2.5456411 -1.3246120 -0.7733470 -1.4208267 -0.1892435 
##          S          T          U          V          W          X 
##  0.6573859  1.3248485  1.7718648  1.9577381  1.3256986  0.9330483 
##          Y          Z 
## -1.5227223 -0.7031367

## visualize
obj<-contMap(tree,x,anc.states=a,plot=FALSE)
plot(obj)

Now let's proceed to estimate ancestral states - first with no known states, then with some, then (in the trivial case) with all:

## first with no internal state data
anc1<-fastAnc(tree,x,CI=TRUE)
anc1

## Ancestral character estimates using fastAnc:
##        27        28        29        30        31        32        33 
## -0.458331 -0.555336 -0.531588 -0.613868 -0.944556 -1.006013 -1.299123 
##        34        35        36        37        38        39        40 
## -2.246533 -0.306729 -0.248593 -0.687350 -0.724508 -1.116524 -1.087959 
##        41        42        43        44        45        46        47 
## -2.407758 -0.965980 -1.127005 -1.085506  0.054002  0.370668  1.315666 
##        48        49        50        51 
##  1.616995  1.585510  1.231706 -0.886161 
## 
## Lower & upper 95% CIs:
##        lower     upper
## 27 -1.732581  0.815918
## 28 -1.618844  0.508172
## 29 -1.568306  0.505131
## 30 -1.680216  0.452481
## 31 -1.940270  0.051157
## 32 -1.984038 -0.027989
## 33 -2.253199 -0.345047
## 34 -2.587840 -1.905225
## 35 -1.420263  0.806805
## 36 -1.348681  0.851494
## 37 -1.331019 -0.043681
## 38 -1.280478 -0.168539
## 39 -1.587909 -0.645139
## 40 -2.170496 -0.005422
## 41 -2.909779 -1.905736
## 42 -1.978196  0.046235
## 43 -1.928715 -0.325295
## 44 -1.771882 -0.399130
## 45 -0.966783  1.074788
## 46 -0.593459  1.334794
## 47  0.864155  1.767176
## 48  1.310525  1.923465
## 49  1.317635  1.853386
## 50  0.829624  1.633787
## 51 -1.607091 -0.165230

## now with 8 (of 25) internal nodes known
anc.states<-sample(a,8)
anc.states

##         37         44         49         28         30         39 
## -0.4662824 -2.0034632  1.5818360 -0.5660124 -1.1966136 -1.5643850 
##         47         34 
##  1.4293262 -2.1403459

anc2<-fastAnc(tree,x,CI=TRUE,anc.states=anc.states)
anc2

## Ancestral character estimates using fastAnc:
##        27        28        29        30        31        32        33 
## -0.462569 -0.566012 -0.626395 -1.196614 -1.243811 -1.227162 -1.383056 
##        34        35        36        37        38        39        40 
## -2.140346 -0.324441 -0.199393 -0.466282 -0.779036 -1.564385 -1.262843 
##        41        42        43        44        45        46        47 
## -2.418075 -1.243738 -1.710369 -2.003463  0.083770  0.411412  1.429326 
##        48        49        50        51 
##  1.631835  1.581836  1.294140 -0.880376 
## 
## Lower & upper 95% CIs:
##        lower     upper
## 27 -1.419954  0.494815
## 28 -0.566012 -0.566012
## 29 -1.133676 -0.119115
## 30 -1.196614 -1.196614
## 31 -2.136108 -0.351515
## 32 -2.181969 -0.272355
## 33 -2.388912 -0.377201
## 34 -2.140346 -2.140346
## 35 -1.385475  0.736594
## 36 -1.326371  0.927585
## 37 -0.466282 -0.466282
## 38 -1.195829 -0.362243
## 39 -1.564385 -1.564385
## 40 -2.354270 -0.171416
## 41 -2.965873 -1.870277
## 42 -2.295524 -0.191951
## 43 -2.441551 -0.979187
## 44 -2.003463 -2.003463
## 45 -1.002757  1.170297
## 46 -0.615960  1.438785
## 47  1.429326  1.429326
## 48  1.379314  1.884356
## 49  1.581836  1.581836
## 50  0.948393  1.639887
## 51 -1.666475 -0.094277

It's probably worth noting that for any node with a fixed state the 95% CI has a width of zero - just as we'd expect.

Finally:

anc3<-fastAnc(tree,x,CI=TRUE,anc.states=a)
anc3

## Ancestral character estimates using fastAnc:
##        27        28        29        30        31        32        33 
##  0.000000 -0.566012 -0.141057 -1.196614 -1.665820 -0.936863 -0.977109 
##        34        35        36        37        38        39        40 
## -2.140346 -0.535783  0.093014 -0.466282 -0.944997 -1.564385  0.247628 
##        41        42        43        44        45        46        47 
## -2.293513 -1.040662 -1.402665 -2.003463  0.085261  0.289981  1.429326 
##        48        49        50        51 
##  1.779218  1.581836  1.314194 -1.169055 
## 
## Lower & upper 95% CIs:
##        lower     upper
## 27  0.000000  0.000000
## 28 -0.566012 -0.566012
## 29 -0.141057 -0.141057
## 30 -1.196614 -1.196614
## 31 -1.665820 -1.665820
## 32 -0.936863 -0.936863
## 33 -0.977109 -0.977109
## 34 -2.140346 -2.140346
## 35 -0.535783 -0.535783
## 36  0.093014  0.093014
## 37 -0.466282 -0.466282
## 38 -0.944997 -0.944997
## 39 -1.564385 -1.564385
## 40  0.247628  0.247628
## 41 -2.293513 -2.293513
## 42 -1.040662 -1.040662
## 43 -1.402665 -1.402665
## 44 -2.003463 -2.003463
## 45  0.085261  0.085261
## 46  0.289981  0.289981
## 47  1.429326  1.429326
## 48  1.779218  1.779218
## 49  1.581836  1.581836
## 50  1.314194  1.314194
## 51 -1.169055 -1.169055

## now
plot(a,anc1$ace,xlab="true states",ylab="estimated states",
    main="No internal nodes fixed") ## no internal node information
lines(c(min(par()$usr[c(1,3)]),max(par()$usr[c(2,4)])),
    c(min(par()$usr[c(1,3)]),max(par()$usr[c(2,4)])),
    lty="dashed") ## 1:1 line

plot(a,anc2$ace,xlab="true states",ylab="estimated states",
    main="Some internal nodes fixed") ## some internal node information
lines(c(min(par()$usr[c(1,3)]),max(par()$usr[c(2,4)])),
    c(min(par()$usr[c(1,3)]),max(par()$usr[c(2,4)])),
    lty="dashed") ## 1:1 line

plot(a,anc3$ace,xlab="true states",ylab="estimated states",
    main="All internal nodes fixed") ## all internal node information
lines(c(min(par()$usr[c(1,3)]),max(par()$usr[c(2,4)])),
    c(min(par()$usr[c(1,3)]),max(par()$usr[c(2,4)])),
    lty="dashed") ## 1:1 line

That's it.

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!