User-supplied bug fix for fastAnc

A phytools user, David Labonte from the University of Cambridge, recently reported the following bug with the phytools function fastAnc:

“I am using your fastAnc function to estimate ancestral states of a continuous variable, and it runs smoothly, however only if CI=FALSE. I also noted that it returns a longer variance than ancestral state estimation vector. Notably, this problem only arises for non-dichotomous trees.”

David also was kind enough to provide the following example which reproduces the error very nicely (modified slightly):

library(phytools)
set.seed(1)
tree<-pbtree(n=10,scale=1)
plotTree(tree,node.numbers=TRUE)

x<-fastBM(tree)
fastAnc(tree,x,vars=TRUE,CI=TRUE) ## runs smoothly without error

## $ace
##          11          12          13          14          15          16 
## -0.35538758 -0.39739912 -0.32825378 -0.51263348 -0.32987276 -0.13938453 
##          17          18          19 
## -0.14780235 -0.07371257 -0.09712277 
## 
## $var
##          11          12          13          14          15          16 
## 0.081863455 0.063948796 0.045742058 0.019817948 0.006889495 0.027865560 
##          17          18          19 
## 0.003468311 0.022557094 0.009799787 
## 
## $CI95
##          [,1]        [,2]
## 11 -0.9161787  0.20540351
## 12 -0.8930459  0.09824762
## 13 -0.7474467  0.09093912
## 14 -0.7885549 -0.23671206
## 15 -0.4925586 -0.16718692
## 16 -0.4665670  0.18779789
## 17 -0.2632314 -0.03237331
## 18 -0.3680854  0.22066021
## 19 -0.2911508  0.09690522

tree<-collapse.to.star(tree,fastMRCA(tree,"t6","t9"))
plotTree(tree,node.numbers=TRUE)

fastAnc(tree,x) ## works no problem

##         11         12         13         14         15         16 
## -0.3498724 -0.3911156 -0.3188186 -0.5111098 -0.3295891 -0.1216447

fastAnc(tree,x,vars=TRUE) ## works, but vars is wrong length

## $ace
##         11         12         13         14         15         16 
## -0.3498724 -0.3911156 -0.3188186 -0.5111098 -0.3295891 -0.1216447 
## 
## $var
##          11          12          13          14          15          16 
## 0.077702044 0.060153127 0.040920299 0.018985301 0.006625218 0.015899879

fastAnc(tree,x,vars=TRUE,CI=TRUE) ## doesn't work at all

## $ace
##         11         12         13         14         15         16 
## -0.3498724 -0.3911156 -0.3188186 -0.5111098 -0.3295891 -0.1216447 
## 
## $var
##          11          12          13          14          15          16 
## 0.077702044 0.060153127 0.040920299 0.018985301 0.006625218 0.015899879 
## 
## $CI95
##          [,1]        [,2]
## 11 -0.8962241  0.19647933
## 12 -0.8718278  0.08959662
## 13 -0.7153024  0.07766520
## 14 -0.7811727 -0.24104701
## 15 -0.4891242 -0.17005403
## 16 -0.3687903  0.12550100

Even better, David solved the bug by correctly identifying the error in the code. In his own words, he says:

“I looked through the code of the function, and while I am certainly not a R coding expert, I believe the problem lies in line 28:

27    if (vars || CI) {
28        v[as.character(ancNames[, 2])]
29        names(v) <- ancNames[, 1]
30    }

analogous to the previous lines, I think this should read

27    if (vars || CI) {
28        v <- v[as.character(ancNames[, 2])]
29        names(v) <- ancNames[, 1]
30    }

This is exactly correct, and if we fix this then we find that the function now works perfectly:

source("fastAnc.R")
fastAnc(tree,x,vars=TRUE,CI=TRUE)

## $ace
##         11         12         13         14         15         16 
## -0.3498724 -0.3911156 -0.3188186 -0.5111098 -0.3295891 -0.1216447 
## 
## $var
##          11          12          13          14          15          16 
## 0.077702044 0.060153127 0.040920299 0.018985301 0.006625218 0.015899879 
## 
## $CI95
##          [,1]        [,2]
## 11 -0.8962241  0.19647933
## 12 -0.8718278  0.08959662
## 13 -0.7153024  0.07766520
## 14 -0.7811727 -0.24104701
## 15 -0.4891242 -0.17005403
## 16 -0.3687903  0.12550100

Cool! If only fixing all the bugs in phytools was this easy!