Newer nonstatic version of phytools, with a bug fix and a new function to rescale SIMMAP style trees

I just posted a newer nonstatic version of 'phytools' (v 0.1-85) with a bug fix to the function export.as.xml (which creates an input file for the stand-alone program SIMMAP v1.5), as well as new function to rescale SIMMAP style modified "phylo" objects in R, sensibly called rescaleSimmap. A direct link to the package build is here

First, the bug fix. Yesterday, I realized that the traditional semi-colon had been cropped from the end of the Newick string in XML format. This makes sense, to some extent, because the end of the string is indicated by a </tree> tag anyway. Pulling the semicolon off the end of the tree string created by ape::write.tree is not quite as simple as it sounds. That is because in R a character string is stored as a single object, rather than a vector of individual characters. To do it, I used strsplit to split up the string, and then paste to put it back together (although in hindsight, I might have used toString instead). This was done basically as follows:

temp<-write.tree(tree) # write Newick to temp
temp<-unlist(strsplit(temp,NULL)) # split into character vector
# paste together, dropping the last character (a ";")
temp<-paste(temp[1:(length(temp)-1)],collapse="")
# write to file with tags
write(paste("\t\t",temp,"",sep=""),file,append=TRUE)

I also added a new function, rescaleSimmap (code here), which rescales a SIMMAP style "phylo" object in R. This function is analagous to rescaleTree in the 'geiger' package. This function takes the input tree and a new total depth (height) for the tree, and then rescales the branch lengths and the mapping elements ($maps and $mapped.edge) proportionally. I programmed this today because I discovered that (at least by default), for a given input tree and character value, SIMMAP samples not only mutational histories, but also total tree depth. Not sure why this is - I may investigate.

Bug in export.as.xml

In trying to use export.as.xml to export data and trees to read into SIMMAP I have discovered a small bug. Evidently (flaunting years of convention) the Newick-style string containing the tree structure in XML format does not contain a semi-colon character. Who new? I will fix this tomorrow.

New nonstatic version of phytools (v 0.1-84)

I just built (for myself) and posted (for any interested 'phytools' users) a new non-static version of phytools (0.1-84). A direct link to the build is here. To install, just download and then type:

> install.packages("phytools_0.1-84.tar.gz",type="source", repos=NULL)

The package has two new functions over the previous non-static version: export.as.xml (described here) and pgls.Ives (described here), as well as probably a few more new function relative to the latest CRAN release of phytools (v 0.1-8).

PGLS regression with sampling error

Following some discussion of the topic on the R-sig-phylo special interest group email listserve, I decided to try and program the PGLS regression method for sampling error in X & Y described originally in a paper by Tony Ives, Peter Midford, and Ted Garland in Systematic Biology (Ives et al. 2007). Even though Ives et al. worked out all the hard parts, this was not a completely trivial undertaking as some of the details were a little difficult to work through. In addition, I first tried to program this as custom correlation structure to use with nlme::gls. It's wholly possible that one might be able to do this in theory, but I was not able to (although I did learn a bit about how these correlation structures are written).

In the end, I did this the old fashioned way, and it seems to work (at least, it converges on standard PGLS when sampling error is zero; and it recovers estimates close to the generating parameter values when sampling error is included for in simulation and accounted for during estimation). Unfortunately, I have so far only programmed bivariate regression. It is theoretically straightforward to extend to multivariable regression, with the caveat that optimizing the likelihood function requires inversion of a covariance matrix whose dimensions scale with N × (p+1) for p independent variables. The optimization is also multidimensional because we simultaneously need to maximize the likelihood the for the rates of evolution in x and y (σ_x², σ_y²), the regression slope (b₁) and the states at the root for x and y (a). We don't need to separately estimate b₀ - we can just get it from the root values and the slope.

Code for the function is here. Let's try it:

> # require dependencies and load source
> require(phytools)
> source("pgls.Ives.R")
> # simulate a phylogeny
> tree<-pbtree(n=200,scale=1)
> # simulate under some arbitrary regression model
> x<-fastBM(tree)
> y<-5+3.5*x+fastBM(tree,0.2)
> # create vectors of zeroes for the "no error" case
> Vx0<-Vy0<-Cxy0<-rep(0,length(tree$tip))
> names(Vx0)<-names(Vy0)<-names(Cxy0)<-tree$tip.label
> # fit the no error model with simulated data
> res1<-pgls.Ives(tree,x,y,Vx0,Vy0,Cxy0)
> res1
$beta
[1] 4.874868 3.392148
$sig2x
[1] 0.9823115
$sig2y
[1] 1.008072
$a
[1] -0.2179666 4.1354924
$logL
[1] -224.2079
> # for reference, we can also fit this with gls()
> res2<-gls(y~x,data.frame(x,y),correlation=corBrownian(1,tree))
> res2
Generalized least squares fit by REML
  Model: y ~ x
  Data: data.frame(x, y)
  Log-restricted-likelihood: -115.1862
Coefficients:
(Intercept)           x
   4.874825    3.392136
Correlation Structure: corBrownian
Formula: ~1
Parameter estimate(s):
numeric(0)
Degrees of freedom: 200 total; 198 residual
Residual standard error: 1.00909
> # simulate sampling variances and covariances
> Vx<-rexp(length(tree$tip))/2
> Vy<-rexp(length(tree$tip))/2
> Cxy=sqrt(Vx*Vy)*runif(length(tree$tip),min=-1,max=1)
> names(Vx)<-names(Vy)<-names(Cxy)<-names(x)
> # simulate data with sampling error
> xe<-ye<-vector()
> for(i in 1:length(tree$tip)){
  temp<-mvrnorm(mu=c(x[i],y[i]),
Sigma=matrix(c(Vx[i],Cxy[i],Cxy[i],Vy[i]),2,2))
  xe[i]<-temp[1]; ye[i]<-temp[2]
  }
> names(xe)<-names(x)
> names(ye)<-names(y)
> # fit no error model to data with error
> res3<-pgls.Ives(tree,xe,ye,Vx0,Vy0,Cxy0)
> res3
$beta
[1] 4.1621231 0.2000734
$sig2x
[1] 27.1866
$sig2y
[1] 24.35597
$a
[1] -0.2924528  4.1036111
$logL
[1] -874.7336
> # fit error model
> res4<-pgls.Ives(tree,xe,ye,Vx,Vy,Cxy)
> res4
$beta
[1] 5.133185 3.485432
$sig2x
[1] 0.9533415
$sig2y
[1] 0.830365
$a
[1] -0.2760699  4.1709624
$logL
[1] -568.1378

Obviously, in this case our estimate - particularly of the regression slope - is quite bad when sampling error is ignored; but we do much better with the method of Ives et al. Cool!

R trick 1: get the frequencies of factors in a vector

Here's a quick R hint. (I had briefly forgotten how to do this, and the solution wasn't totally obvious online.) Say I have a vector of factors in memory in R and I want to get the frequency or relative frequency of the different levels of the factor, I can do this using the base generic function summary. Just to see how this works, consider a vector containing the best-fitting quantitative trait evolution model for a set of 100 trees:

> best.fit
 [1] BM     BM     OU     lambda BM     BM     OU     BM
 [9] BM     BM     BM     lambda lambda BM     OU     BM
[17] BM     BM     OU     BM     BM     BM     lambda BM
[25] BM     BM     BM     OU     lambda BM     BM     BM
[33] BM     OU     BM     BM     lambda lambda lambda BM
[41] BM     BM     BM     OU     BM     BM     BM     BM
[49] BM     OU     BM     BM     BM     BM     BM     BM
[57] lambda lambda OU     OU     BM     BM     lambda BM
[65] BM     BM     BM     BM     lambda BM     BM     BM
[73] BM     OU     OU     BM     lambda BM     lambda BM
[81] BM     lambda BM     BM     BM     OU     BM     BM
[89] BM     BM     OU     OU     lambda BM     BM     BM
[97] BM     BM     OU     BM
Levels: BM lambda OU

We can count up the number or relative frequency of trees with each best fit model as follows:

> summary(best.fit)
   BM lambda     OU
   68     16     16
> summary(best.fit)/sum(summary(best.fit))
   BM lambda     OU
 0.68   0.16   0.16

That's it.

In the field

I'm presently in the field in Puerto Rico through June 26th, so I may not be blogging much on phytools; however I did post a short field story to Anole Annals (see link here) and I may do that a few more times before I'm done. BTW, the picture below is of Anolis krugi (the "mountain" grass-bush anole in Puerto Rico, although we are only at about 200 feet above sea level), one of the species we will be studying on the island.

New function to export data & trees in XML format for SIMMAP

Phytools has a function, make.simmap, that can be used to generate stochastic character mapped trees for discretely valued character data and phylogenies. This function, though, is really a much simplified version of Jonathan Bollback's stand alone program SIMMAP. The biggest difference between stochastic mapping in his program and with phytools::make.simmap is that my function first fits a model of character evolution using ML, and then samples from the conditional posterior distribution of character histories (conditioning on the fitted model), rather than the joint posterior distribution of substitution rates and character histories.

Some phytools users might therefore be interested in exporting their character data & trees for analysis in SIMMAP. In fact, I can attest to this because I want to do it - mainly to cross-validate a result that I obtained using make.simmap in phytools. Phytools can already be used to read in stochastic histories generated using SIMMAP (now both v1.0 and v1.5).

SIMMAP (v1.5, I'll get to v1.0 later) reads trees and data in XML format. XML is a markup language, like HTML, using tags for different components of the data structure. SIMMAP can read both DNA and discrete, numerically coded morphological traits. The former is called datatype="nucleotide" by SIMMAP and the latter datatype="standard" by SIMMAP.

I have now posted a function online, export.as.xml, that can create an XML datafile using data stored in various formats and a tree or multiple trees stored as a "phylo" or "multiPhylo" object (respectively).

A direct link to the code is here. I will also most likely include it in future releases of phytools. (I may also include the capacity to export for SIMMAP v1.0. SIMMAP v1.0 does not use XML, but a modified NEXUS format, and is still in fairly wide use.)