R/ltt.R

In liamrevell/phytools: Phylogenetic Tools for Comparative Biology (and Other Things)

## This function computes the data for a lineage-through-time plot and 
## (optionally) creates this plot the function does not require a tree 
## that is ultrametric.  Optionally, the function can remove extinct
## species from the phylogeny. If the input tree is an object of class 
## "multiPhylo" then the function will simultaneously plot all ltts.
## written by Liam J. Revell 2010-2015, 2022, 2023

## ltt method (added 2022)

ltt<-function(tree,...) UseMethod("ltt")

ltt.default<-function(tree,...){
	warning(paste(
		"ltt does not know how to handle objects of class ",
		class(tree),"."))
}

print.ltt.multiSimmap<-function(x,...){
	cat(paste(length(x),"objects of class \"ltt.simmap\" in a list\n"))
}

ltt.multiSimmap<-function(tree,gamma=TRUE,...){
	if(!inherits(tree,"multiSimmap")){
		stop("tree must be object of class \"multiSimmap\".")
	} else {
		obj<-lapply(tree,ltt,plot=FALSE,log.lineages=FALSE,gamma=gamma)
		class(obj)<-"ltt.multiSimmap"
		return(obj)
	}
}

## internally used functions
BRANCHING<-function(phy,is_ultrametric){
	x<-if(is_ultrametric) branching.times(phy)
	else { 
		sort(setNames(max(nodeHeights(phy))-sapply(1:phy$Nnode+Ntip(phy),
			nodeheight,tree=phy),1:phy$Nnode+Ntip(phy)))
	}
	x
}
TIPHEIGHTS<-function(phy,is_ultrametric){
	x<-if(is_ultrametric) {
		min(setNames(sapply(1:Ntip(phy),nodeheight,tree=phy),1:Ntip(phy)))
	} else setNames(sapply(1:Ntip(phy),nodeheight,tree=phy),1:Ntip(phy))
	x
}

ltt.simmap<-function(tree,plot=TRUE,log.lineages=FALSE,gamma=TRUE,...){
	## set tolerance
	if(hasArg(tol)) tol<-list(...)$tol else tol<-.Machine$double.eps^0.5
	## check if ultrametric
	is_ultrametric<-is.ultrametric(tree)
	if(!inherits(tree,"simmap")){
		stop("tree must be an object of class \"simmap\".")
	} else {
		levs<-sort(unique(c(getStates(tree,"tips"),
			getStates(tree,"nodes"))))
		tt<-map.to.singleton(tree)
		H<-nodeHeights(tt)
		h<-c(0,max(H)-BRANCHING(tt,is_ultrametric),TIPHEIGHTS(tt,is_ultrametric))
		ss<-setNames(as.factor(names(tt$edge.length)),
			tt$edge[,2])
		lineages<-matrix(0,length(h),length(levs),
			dimnames=list(names(h),levs))
		lineages[1,getStates(tree,"nodes")[1]]<-1
		ROOT<-Ntip(tt)+1
		for(i in 2:length(h)){
			if(h[i]<max(h)) ii<-intersect(which(h[i]>=H[,1]),which(h[i]<H[,2]))
			else ii<-intersect(which(h[i]>=H[,1]),which(h[i]<=H[,2]))
			lineages[i,]<-summary(ss[ii])
		}
		ii<-order(h)
		times<-h[ii]
		lineages<-lineages[ii,,drop=FALSE]
		lineages<-cbind(lineages,total=rowSums(lineages))
		obj<-list(times=times,ltt=lineages)	
		if(gamma==FALSE){
			obj<-list(ltt=lineages,times=times,tree=tree)
			class(obj)<-"ltt.simmap"
		} else {
			gam<-gammatest(ltt(as.phylo(tree),plot=FALSE))
			obj<-list(ltt=lineages,times=times,gamma=gam$gamma,
				p=gam$p,tree=tree)
			class(obj)<-"ltt.simmap"
		}
	}
	if(plot) plot(obj,log.lineages=log.lineages,...)
	obj
}

plot.ltt.multiSimmap<-function(x,...){
	if(hasArg(add)) add<-list(...)$add
	else add<-FALSE
	hh<-max(sapply(x,function(x) max(nodeHeights(x$tree))))
	if(hasArg(alpha)) alpha<-list(...)$alpha else alpha<-0.05
	if(hasArg(res)) res<-list(...)$res else res<-1000
	if(hasArg(log.lineages)) log.lineages<-list(...)$log.lineages
	else log.lineages<-FALSE
	levs<-sort(unique(unlist(lapply(x,function(x) 
		c(getStates(x$tree,"tips"),getStates(x$tree,"nodes"))))))
	if(hasArg(colors)) colors<-list(...)$colors	
	else colors<-setNames(c(palette()[1:length(levs)+1],par()$fg),
		c(levs,"total"))
	if(hasArg(legend)) legend<-list(...)$legend else 
		legend<-"topleft"
	plot.leg<-TRUE
	if(is.logical(legend)) if(legend) plot.leg<-TRUE else 
		plot.leg<-FALSE
	if(hasArg(lwd)) lwd<-list(...)$lwd else lwd<-5
	if(hasArg(show.total)) show.total<-list(...)$show.total else
		show.total<-TRUE
	nn<-max(sapply(x,function(x,tot) if(!tot) 
		max(x$ltt[,-which(colnames(x$ltt)=="total")]) else
		Ntip(x$tree),tot=show.total))
	xlim<-if(hasArg(xlim)) list(...)$xlim else c(0,hh)
	ylim<-if(hasArg(ylim)) list(...)$ylim else 
		if(log.lineages) log(c(1,1.05*nn)) else 
		c(0,1.05*nn)
	xlab<-if(hasArg(xlab)) list(...)$xlab else "time"
	ylab<-if(hasArg(ylab)) list(...)$ylab else if(log.lineages) 
		"log(lineages)" else "lineages"
	TIMES<-seq(0,hh,length.out=res)
	LINEAGES<-matrix(0,length(TIMES),length(levs)+1)
	colnames(LINEAGES)<-c(levs,"total")
	for(i in 1:length(TIMES)){
		for(j in 1:length(x)){
			ii<-which(x[[j]]$times<=TIMES[i])
			ADD<-if(length(ii)==0) rep(0,length(levs)) else 
				x[[j]]$ltt[max(ii),]/length(x)
			LINEAGES[i,]<-LINEAGES[i,]+ADD
		}
	}
	args<-list(...)
	args$res<-NULL
	args$alpha<-NULL
	args$log.lineages<-NULL
	args$colors<-NULL
	args$legend<-NULL
	args$show.total<-NULL
	args$add<-NULL
	args$xlim<-xlim
	args$ylim<-ylim
	args$xlab<-xlab
	args$ylab<-ylab
	args$x<-NA
	if(!add) do.call(plot,args)
	if(!show.total) dd<-1 else dd<-0
	for(i in 1:length(x)){
		for(j in 1:(ncol(LINEAGES)-dd)){
			nm<-colnames(x[[i]]$ltt)[j]
			ltt<-if(log.lineages) log(x[[i]]$ltt[,j]) else x[[i]]$ltt[,j]
			lines(x[[i]]$times,ltt,type="s",lwd=1,
				col=make.transparent(colors[nm],
				alpha))
		}
	}
	for(i in 1:(ncol(LINEAGES)-dd)){
		nm<-colnames(LINEAGES)[i]
		ltt<-if(log.lineages) log(LINEAGES[,i]) else LINEAGES[,i]
		lines(TIMES,ltt,lwd=lwd,col=colors[nm])
	}
	if(plot.leg){
		nm<-c(levs,"total")
		if(!show.total) nm<-setdiff(nm,"total")
		legend(legend,legend=nm,pch=22,pt.bg=colors[nm],pt.cex=1.2,
			cex=0.8,bty="n")
	}
	invisible(list(times=TIMES,ltt=LINEAGES))
}

plot.ltt.simmap<-function(x,...){
	if(hasArg(add)) add<-list(...)$add else add<-FALSE
	if(hasArg(log.lineages)) log.lineages<-list(...)$log.lineages
	else log.lineages<-FALSE
	if(hasArg(colors)) colors<-list(...)$colors	
	else colors<-setNames(c(palette()[2:ncol(x$ltt)],par()$fg),
		colnames(x$ltt))
	if(hasArg(legend)) legend<-list(...)$legend else 
		legend<-"topleft"
	plot.leg<-TRUE
	if(is.logical(legend)) if(legend) plot.leg<-TRUE else 
		plot.leg<-FALSE
	if(hasArg(show.tree)) show.tree<-list(...)$show.tree else
		show.tree<-FALSE
	if(hasArg(lwd)) lwd<-list(...)$lwd else lwd<-3
	if(hasArg(outline)) outline<-list(...)$outline else 
		outline<-show.tree
	if(hasArg(show.total)) show.total<-list(...)$show.total else
		show.total<-TRUE
	xlim<-if(hasArg(xlim)) list(...)$xlim else range(x$times)
	ylim<-if(hasArg(ylim)) list(...)$ylim else 
		if(log.lineages) log(c(1,1.05*max(x$ltt))) else 
		c(0,1.1*max(x$ltt))
	xlab<-if(hasArg(xlab)) list(...)$xlab else "time"
	ylab<-if(hasArg(ylab)) list(...)$ylab else if(log.lineages) 
		"log(lineages)" else "lineages"
	args<-list(...)
	args$log.lineages<-NULL
	args$colors<-NULL
	args$legend<-NULL
	args$show.tree<-NULL
	args$show.total<-NULL
	args$add<-NULL
	args$xlim<-xlim
	args$ylim<-ylim
	args$xlab<-xlab
	args$ylab<-ylab
	args$x<-NA
	if(!add) do.call(plot,args)
	tips<-if(log.lineages) seq(0,log(Ntip(x$tree)),
		length.out=Ntip(x$tree)) else 1:Ntip(x$tree)
	if(show.tree){
		mar<-par()$mar
		plot(x$tree,
			make.transparent(colors[1:(ncol(x$ltt)-1)],0.5),
			tips=tips,xlim=xlim,ylim=ylim,
			ftype="off",add=TRUE,lwd=1,mar=mar)
		plot.window(xlim=xlim,ylim=ylim)
	}
	if(!show.total) dd<-1 else dd<-0
	for(i in 1:(ncol(x$ltt)-dd)){
		LTT<-if(log.lineages) log(x$ltt) else x$ltt
		if(outline) lines(x$times,LTT[,i],type="s",lwd=lwd+2,
				col=if(par()$bg=="transparent") "white" else 
				par()$bg)
		lines(x$times,LTT[,i],type="s",lwd=lwd,col=colors[i])
	}
	if(plot.leg){
		nn<-colnames(x$ltt)
		if(!show.total) nn<-setdiff(nn,"total")
		legend(legend,nn,pch=22,pt.bg=colors[nn],pt.cex=1.2,
		cex=0.8,bty="n")
	}
}

print.ltt.simmap<-function(x,digits=4,...){
	ss<-sort(unique(c(getStates(x$tree,"tips"),
		getStates(x$tree,"nodes"))))
	cat("Object of class \"ltt.simmap\" containing:\n\n")
	cat(paste("(1) A phylogenetic tree with ",Ntip(x$tree), 
		" tips, ",x$tree$Nnode," internal\n",sep= ""))
	cat(paste("    nodes, and a mapped state with ",length(ss),
		" states.\n\n",sep=""))
	cat(paste("(2) A matrix containing the number of lineages in each\n",
		"    state (ltt) and event timings (times) on the tree.\n\n",sep=""))
	if(!is.null(x$gamma))
		cat(paste("(3) A value for Pybus & Harvey's \"gamma\"",
			" statistic of\n    gamma = ",round(x$gamma,digits),
			", p-value = ",
			round(x$p,digits),".\n\n",sep=""))
}

ltt.multiPhylo<-function(tree,drop.extinct=FALSE,gamma=TRUE,...){
	if(!inherits(tree,"multiPhylo")){
		stop("tree must be object of class \"multiPhylo\".")
	} else {
		obj<-lapply(tree,ltt,plot=FALSE,drop.extinct=drop.extinct,
			log.lineages=FALSE,gamma=gamma)
		class(obj)<-"multiLtt"
		return(obj)
	}
}

ltt.phylo<-function(tree,plot=TRUE,drop.extinct=FALSE,log.lineages=TRUE,
	gamma=TRUE,...){
	# set tolerance
	tol<-1e-6
	# check "phylo" object
	if(!inherits(tree,"phylo")){
		stop("tree must be object of class \"phylo\".")
	} else {
		# reorder the tree
		tree<-reorder.phylo(tree,order="cladewise")
		if(!is.null(tree$node.label)){ 
			node.names<-setNames(tree$node.label,1:tree$Nnode+Ntip(tree))
			tree$node.label<-NULL
		} else node.names<-NULL
		## check if tree is ultrametric & if yes, then make *precisely* ultrametric
		if(is.ultrametric(tree)){
			h<-max(nodeHeights(tree))
			time<-c(0,h-sort(branching.times(tree),decreasing=TRUE),h)
			nodes<-as.numeric(names(time)[2:(length(time)-1)])
			ltt<-c(cumsum(c(1,sapply(nodes,function(x,y) sum(y==x)-1,y=tree$edge[,1]))),
				length(tree$tip.label))
			names(ltt)<-names(time)
		} else {
			# internal functions
			# drop extinct tips
			drop.extinct.tips<-function(phy){
				temp<-diag(vcv(phy))
				if(length(temp[temp<(max(temp)-tol)])>0)
					pruned.phy<-drop.tip(phy,names(temp[temp<(max(temp)-tol)]))
				else
					pruned.phy<-phy
				return(pruned.phy)
			}
			# first, if drop.extinct==TRUE
			if(drop.extinct==TRUE)
				tree<-drop.extinct.tips(tree)
			# compute node heights
			root<-length(tree$tip)+1
			node.height<-matrix(NA,nrow(tree$edge),2)
			for(i in 1:nrow(tree$edge)){
				if(tree$edge[i,1]==root){
					node.height[i,1]<-0.0
					node.height[i,2]<-tree$edge.length[i]
				} else {
					node.height[i,1]<-node.height[match(tree$edge[i,1],tree$edge[,2]),2]
					node.height[i,2]<-node.height[i,1]+tree$edge.length[i]
				}
			}
			ltt<-vector()
			tree.length<-max(node.height) # tree length
			n.extinct<-sum(node.height[tree$edge[,2]<=length(tree$tip),2]<(tree.length-tol))
			# fudge things a little bit
			node.height[tree$edge[,2]<=length(tree$tip),2]<-
				node.height[tree$edge[,2]<=length(tree$tip),2]+1.1*tol
			time<-c(0,node.height[,2]); names(time)<-as.character(c(root,tree$edge[,2]))
			temp<-vector()
			time<-time[order(time)]
			time<-time[1:(tree$Nnode+n.extinct+1)] # times
			# get numbers of lineages
			for(i in 1:(length(time)-1)){
				ltt[i]<-0
				for(j in 1:nrow(node.height))
					ltt[i]<-ltt[i]+(time[i]>=(node.height[j,1]-
						tol)&&time[i]<=(node.height[j,2]-tol))
			}
			ltt[i+1]<-0
			for(j in 1:nrow(node.height))
				ltt[i+1]<-ltt[i+1]+(time[i+1]<=(node.height[j,2]+tol))
			# set names (these are the node indices)
			names(ltt)<-names(time)
			# append 0,1 for time 0
			ltt<-c(1,ltt)
			time<-c(0,time)
			# subtract fudge factor
			time[length(time)]<-time[length(time)]-1.1*tol
		}
		if(!is.null(node.names)){ nn<-sapply(names(time),function(x,y) 
			if(any(names(y)==x)) y[which(names(y)==x)] else "",y=node.names)
			names(ltt)<-names(time)<-nn
		}
		if(gamma==FALSE){
			obj<-list(ltt=ltt,times=time,tree=tree)
			class(obj)<-"ltt"
		} else {
			gam<-gammatest(list(ltt=ltt,times=time))
			obj<-list(ltt=ltt,times=time,gamma=gam$gamma,p=gam$p,tree=tree)
			class(obj)<-"ltt"
		}
	}
	if(plot) plot(obj,log.lineages=log.lineages,...)
	obj
}

## function computes the gamma-statistic & a two-tailed P-value
## written by Liam J. Revell 2011, 2019

gammatest<-function(x){
	n<-max(x$ltt)
	g<-vector()
	for(i in 2:(length(x$times))) g[i-1]<-x$times[i]-x$times[i-1]
	T<-sum((2:n)*g[2:n])
	doublesum<-0
	for(i in 1:(n-1)) for(k in 1:i) doublesum<-doublesum+k*g[k]
	gamma<-(1/(n-2)*doublesum-T/2)/(T*sqrt(1/(12*(n-2))))
	p<-2*pnorm(abs(gamma),lower.tail=F)
	object<-list(gamma=gamma,p=p)
	class(object)<-"gammatest"
	object
}

## print method for object class "gammatest"
## written by Liam J. Revell 2019

print.gammatest<-function(x,...){
	if(hasArg(digits)) digits<-list(...)$digits
	else digits<-4
	cat("\nAn object of class \"gammatest\" with:\n")
	cat(paste("(1) Pybus & Harvey's gamma = ",
		round(x$gamma,digits),sep=""))
	cat(paste("\n(2) p-value = ",round(x$p,digits),
		"\n\n",sep=""))
}

## S3 print method for object of class "ltt"
## written by Liam J. Revell 2015

print.ltt<-function(x,digits=4,...){
	cat("Object of class \"ltt\" containing:\n\n")
	cat(paste("(1) A phylogenetic tree with ",Ntip(x$tree), 
		" tips and ",x$tree$Nnode," internal\n",sep= ""))
	cat("    nodes.\n\n")
	cat(paste("(2) Vectors containing the number of lineages (ltt) and\n",
		"    branching times (times) on the tree.\n\n",sep=""))
	if(!is.null(x$gamma))
		cat(paste("(3) A value for Pybus & Harvey's \"gamma\"",
			" statistic of\n    gamma = ",round(x$gamma,digits),
			", p-value = ",
			round(x$p,digits),".\n\n",sep=""))
}

## S3 print method for object of class "multiLtt"
## written by Liam J. Revell 2015

print.multiLtt<-function(x,...){
	cat(paste(length(x),"objects of class \"ltt\" in a list\n"))
}

## S3 plot method for object of class "ltt"
## written by Liam J. Revell 2015, 2023

plot.ltt<-function(x,...){
	args<-list(...)
	if(hasArg(show.tree_mode)) show.tree_mode<-list(...)$show.tree_mode
	else show.tree_mode<-"classic"
	args$show.tree_mode<-NULL
	args$x<-x$time
	if(!is.null(args$log.lineages)){ 
		logl<-args$log.lineages
		args$log.lineages<-NULL
	} else logl<-TRUE
	if(!is.null(args$show.tree)){
		show.tree<-args$show.tree
		args$show.tree<-NULL
	} else show.tree<-FALSE
	if(!is.null(args$transparency)){
		transparency<-args$transparency
		args$transparency<-NULL
	} else transparency<-0.3
	args$y<-if(logl) log(x$ltt) else x$ltt
	if(is.null(args$xlab)) args$xlab<-"time"
	if(is.null(args$ylab)) 
		args$ylab<-if(logl) "log(lineages)" else "lineages"
	if(is.null(args$type)) args$type<-"s"
	if(hasArg(add)){ 
		add<-list(...)$add
		args$add<-NULL
	} else add<-FALSE
	if(!add) do.call(plot,args)
	else do.call(lines,args)
	if(show.tree){
		tips<-if(par()$ylog) exp(setNames(seq(log(min(args$y)),log(max(args$y)),
			length.out=Ntip(x$tree)),x$tree$tip.label)) else setNames(1:Ntip(x$tree),
			x$tree$tip.label)
		if(show.tree_mode=="classic"){
			plotTree(x$tree,color=rgb(0,0,1,transparency),
				ftype="off",add=TRUE,mar=par()$mar,tips=tips)
		} else if(show.tree_mode=="next generation"){
			## this doesn't really work yet
			plotTree(x$tree,color=rgb(0,0,1,transparency),
				ftype="off",add=TRUE,mar=par()$mar,tips=tips,
				plot=FALSE)
			pp<-get("last_plot.phylo",envir=.PlotPhyloEnv)
			max_yy<-max(pp$yy)
			max_y<-max(args$y)
			hh<-c(tips,(pp$yy[1:x$tree$Nnode+Ntip(x$tree)]/max_yy)*
				args$y[as.character(1:x$tree$Nnode+Ntip(x$tree))]+1)
			for(i in 1:nrow(x$tree$edge)) lines(pp$xx[x$tree$edge[i,]],
				hh[x$tree$edge[i,]],col=rgb(0,0,1,transparency),lwd=2)
		}
	}
}

## S3 plot method for object of class "multiLtt"
## written by Liam J. Revell 2015, 2020

plot.multiLtt<-function(x,...){
	max.lineages<-max(sapply(x,function(x) max(x$ltt)))
	max.time<-max(sapply(x,function(x) max(x$time)))
	args<-list(...)
	if(!is.null(args$log.lineages)) logl<-list(...)$log.lineages
	else logl<-TRUE	
	if(is.null(args$xlim)) args$xlim<-c(0,max.time)
	if(is.null(args$ylim)) 
		args$ylim<-if(logl) c(0,log(max.lineages)) else c(1,max.lineages)
	args$x<-x[[1]]
	do.call(plot,args)
	args$add<-TRUE
	if(!is.null(args$log)) args$log<-NULL
	if(length(x)>2){
		for(i in 2:length(x)){
			args$x<-x[[i]]
			do.call(plot,args)
		}
	} else {
		args$x<-x[[2]]
		do.call(plot,args)
	}
}

## compute the gamma statistic through time by slicing the tree n times
## written by Liam J. Revell 2017

gtt<-function(tree,n=100,...){
	if(!inherits(tree,"phylo")) 
		stop("tree must be object of class \"phylo\".")
	if(inherits(tree,"simmap")) tree<-as.phylo(tree)
	if(hasArg(plot)) plot<-list(...)$plot
	else plot<-FALSE
	obj<-ltt(tree,plot=FALSE)
	t<-obj$times[which(obj$ltt==3)[1]]
	h<-max(nodeHeights(tree))
	x<-seq(t,h,by=(h-t)/(n-1))
	trees<-lapply(x,treeSlice,tree=tree,orientation="rootwards")
	gamma<-sapply(trees,function(x,plot){ 
		obj<-unlist(gammatest(ltt<-ltt(x,plot=FALSE)));
		if(plot) plot(ltt,xlim=c(0,h),ylim=c(1,Ntip(tree)),
			log.lineages=FALSE,log="y");
		Sys.sleep(0.01);
		obj},plot=plot)
	object<-list(t=x,gamma=gamma[1,],p=gamma[2,],tree=tree)
	class(object)<-"gtt"
	object
}

## plot method for "gtt" object class
plot.gtt<-function(x,...){
	args<-list(...)
	args$x<-x$t
	args$y<-x$gamma
	if(!is.null(args$show.tree)){ 
		show.tree<-args$show.tree
		args$show.tree<-NULL
	} else show.tree<-TRUE
	if(is.null(args$xlim)) args$xlim<-c(0,max(x$t))
	if(is.null(args$xlab)) args$xlab<-"time"
	if(is.null(args$ylab)) args$ylab<-expression(gamma)
	if(is.null(args$lwd)) args$lwd<-3
	if(is.null(args$type)) args$type<-"s"
	if(is.null(args$bty)) args$bty<-"l"
	if(is.null(args$main)) args$main<-expression(paste(gamma," through time plot"))
	do.call(plot,args)
	if(show.tree) plotTree(x$tree,add=TRUE,ftype="off",mar=par()$mar,
		xlim=args$xlim,color=make.transparent("blue",0.1))
}

## print method for "gtt" object class
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 \n    gamma = ",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),"\n    simulations.\n\n",sep=""))
}

## plot method for "mccr" object class

plot.mccr<-function(x,...){
	if(hasArg(main)) main<-list(...)$main
	else main=expression(paste("null distribution of ",
		gamma))
	if(hasArg(las)) las<-list(...)$las
	else las<-par()$las
	if(hasArg(cex.axis)) cex.axis<-list(...)$cex.axis
	else cex.axis<-par()$cex.axis
	if(hasArg(cex.lab)) cex.lab<-list(...)$cex.lab
	else cex.lab<-par()$cex.lab
	hh<-hist(x$null.gamma,breaks=min(c(max(12,
		round(length(x$null.gamma)/10)),20)),
		plot=FALSE)
	if(hasArg(ylim)) ylim<-list(...)$ylim
	else ylim<-c(0,1.15*max(hh$counts))
	plot(hh,xlim=range(c(x$gamma,x$null.gamma)),
		main=main,xlab=expression(gamma),col="lightgrey",
		ylim=ylim,las=las,cex.axis=cex.axis,cex.lab=cex.lab)
	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.96*par()$usr[4],
		expression(paste("observed value of ",gamma)),
		pos=if(x$gamma>mean(x$null.gamma)) 2 else 4,
		cex=0.9)
}