new function mccr & methods

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: phytools
-Version: 0.6-65
-Date: 2018-12-17
+Version: 0.6-66
+Date: 2019-1-19
 Title: Phylogenetic Tools for Comparative Biology (and Other Things)
 Author: Liam J. Revell
 Maintainer: Liam J. Revell <liam.revell@umb.edu>
@@ -13,6 +13,6 @@ ZipData: no
 Description: A wide range of functions for phylogenetic analysis. Functionality is concentrated in phylogenetic comparative biology, but also includes a diverse array of methods for visualizing, manipulating, reading or writing, and even inferring phylogenetic trees and data. Included among the functions in phylogenetic comparative biology are various for ancestral state reconstruction, model-fitting, simulation of phylogenies and data, and multivariate analysis. There are a broad range of plotting methods for phylogenies and comparative data which include, but are not restricted to, methods for mapping trait evolution on trees, for projecting trees into phenotypic space or a geographic map, and for visualizing correlated speciation between trees. Finally, there are a number of functions for reading, writing, analyzing, inferring, simulating, and manipulating phylogenetic trees and comparative data not covered by other packages. For instance, there are functions for randomly or non-randomly attaching species or clades to a phylogeny, for estimating supertrees or consensus phylogenies from a set, for simulating trees and phylogenetic data under a range of models, and for a wide variety of other manipulations and analyses that phylogenetic biologists might find useful in their research.
 License: GPL (>= 2)
 URL: http://github.com/liamrevell/phytools
-Packaged: 2018-12-17 12:00:00 EST
+Packaged: 2019-1-19 12:00:00 EST
 Repository: 
-Date/Publication: 2018-12-17 12:00:00 EST
+Date/Publication: 2019-1-19 12:00:00 EST

NAMESPACE

@@ -14,7 +14,7 @@ export(getParent, getSisters, getStates, gtt)
 export(labelnodes, ladderize.simmap, lambda.transform, lik.bd, likMlambda, likSurface.rateshift, linklabels, locate.fossil, locate.yeti)
 export(ls.consensus, ls.tree, ltt, ltt95)
 export(make.era.map, make.simmap, make.transparent, map.overlap, Map.Overlap, map.to.singleton, mapped.states, markChanges)
-export(matchLabels, matchNodes, mcmcMk, mergeMappedStates, midpoint.root, minRotate, minSplit, minTreeDist, modified.Grafen)
+export(matchLabels, matchNodes, mccr, mcmcMk, mergeMappedStates, midpoint.root, minRotate, minSplit, minTreeDist, modified.Grafen)
 export(mrp.supertree, multi.mantel, multiC, multiOU, multiRF)
 export(nodeheight, nodeHeights, nodelabels.cophylo, node.paths)
 export(optim.phylo.ls, orderMappedEdge)
@@ -146,6 +146,8 @@ S3method(plot, gtt)
 S3method(print, gtt)
 S3method(print, mcmcMk)
 S3method(print, Dtest)
+S3method(print, mccr)
+S3method(plot, mccr)
 
 importFrom(animation, ani.options, ani.record, ani.replay, saveVideo)
 importFrom(ape, .PlotPhyloEnv, .uncompressTipLabel, ace, all.equal.phylo, as.DNAbin, as.phylo, bind.tree, branching.times, collapse.singles)

R/ltt.R

@@ -246,3 +246,43 @@ plot.gtt<-function(x,...){
 print.gtt<-function(x,...)
 	cat("Object of class \"gtt\".\n\n")
 
+## perform the MCCR test of Pybus & Harvey (2000)
+## written by Liam J. Revell 2018
+
+mccr<-function(obj,rho=1,nsim=100,...){
+	N<-round(Ntip(obj$tree)/rho)
+	tt<-pbtree(n=N,nsim=nsim)
+	foo<-function(x,m) drop.tip(x,sample(x$tip.label,m))
+	tt<-lapply(tt,foo,m=N-Ntip(obj$tree))
+	g<-sapply(tt,function(x) ltt(x,plot=FALSE)$gamma)
+	P<-if(obj$gamma>median(g)) 2*mean(g>=obj$gamma) else 2*mean(g<=obj$gamma)
+	result<-list(gamma=obj$gamma,"P(two-tailed)"=P,null.gamma=g)
+	class(result)<-"mccr"
+	result
+}
+
+## print method for "mccr" object class
+
+print.mccr<-function(x,digits=4,...){
+	cat("Object of class \"mccr\" consisting of:\n\n")
+	cat(paste("(1) A value for Pybus & Harvey's \"gamma\"",
+		" statistic of ",round(x$gamma,digits),".\n\n",sep=""))
+	cat(paste("(2) A two-tailed p-value from the MCCR test of ",
+		round(x$'P(two-tailed)',digits),".\n\n", sep = ""))
+	cat(paste("(3) A simulated null-distribution of gamma from ",
+		length(x$null.gamma)," simulations.\n\n",sep=""))
+}
+
+## plot method for "mccr" object class
+
+plot.mccr<-function(x,...){
+	hist(x$null.gamma,breaks=min(c(max(12,round(length(x$null.gamma)/10)),20)),
+		xlim=range(c(x$gamma,x$null.gamma)),
+		main=expression(paste("null distribution of ",
+		gamma)),xlab=expression(gamma),col="lightgrey")
+	arrows(x0=x$gamma,y0=par()$usr[4],y1=0,length=0.12,
+		col=make.transparent("blue",0.5),lwd=2)
+	text(x$gamma,0.98*par()$usr[4],
+		expression(paste("observed value of ",gamma)),
+		pos=if(x$gamma>mean(x$null.gamma)) 2 else 4)
+}

man/ltt.Rd

@@ -1,10 +1,12 @@
 \name{ltt}
 \alias{ltt}
 \alias{gtt}
+\alias{mccr}
 \title{Creates lineage-through-time plot (including extinct lineages)}
 \usage{
 ltt(tree, plot=TRUE, drop.extinct=FALSE, log.lineages=TRUE, gamma=TRUE, ...)
 gtt(tree, n=100, ...)
+mccr(obj, rho=1, nsim=100, ...)
 }
 \arguments{
 	\item{tree}{is a phylogenetic tree in \code{"phylo"} format, or an object of class \code{"multiPhylo"} containing a list of phylogenetic trees.}
@@ -13,12 +15,17 @@ gtt(tree, n=100, ...)
 	\item{log.lineages}{logical value indicating whether LTT plot should be on log-linear (default) or linear-linear scale.}
 	\item{gamma}{logical value indicating whether or not to compute eqn{gamma} from Pybus & Harvey (2000; \emph{Proc. Roy. Soc. B}).}
 	\item{n}{for \code{gtt} the number of time intervals to use to track \eqn{\gamma} through time.}
+	\item{obj}{for \code{mccr} an object of class \code{"ltt"}.}
+	\item{rho}{for \code{mccr} sampling fraction.}
+	\item{nsim}{for \code{mccr} number of simulations to use for the MCCR test.}
 	\item{...}{other arguments to be passed to plotting methods. See \code{\link{plot.default}}.}
 }
 \description{
 	The function \code{ltt} computes LTT plot with extant and extinct lineages, and optionally conducts \eqn{\gamma}-test of Pybus & Harvey (2000). The object returned by \code{ltt} can be plotted or re-plotted using \code{\link{plot}}.
 
 	The function \code{gtt} computes the value of Pybus & Harvey's \eqn{\gamma} statistic through time by slice the tree at various points - by default in even intervals from the time above the root at which \emph{N} = 3 to the present day.
+	
+	The function \code{mccr} performs the MCCR test of Pybus & Harvey (2000) which takes into account incomplete taxon sampling in computing a P-value of the \eqn{\gamma} statistic.
 }
 \details{
 	Although it is calculated here, it's unclear how to interpret the \eqn{\gamma}-statistic if not all the tips in the tree are contemporaneous.
@@ -32,6 +39,8 @@ gtt(tree, n=100, ...)
 	If \code{tree} is an object of class \code{"multiPhylo"}, then an object of class \code{"multiLtt"} is returned consisting of a list of object of class \code{"ltt"}.
 
 	\code{gtt} returns an object of class \code{"gtt"}.
+
+	\code{mccr} returns of object of class \code{"mccr"}.
 }
 \references{
 	Pybus, O. G., and P. H. Harvey (2000) Testing macro-evolutionary models using incomplete molecular phylogenies. \emph{Proc. R. Soc. Lond. B}, \bold{267}, 2267-2272.