#Kcalc randomization test, imputs: phylogeny object, com (phylocom format:community name, abundance, species name), number of randomizations (n.rand), whether a subtree needs to be created (default sub=FALSE), and to also output an histogram (default fig=TRUE). It can be used in a loop or tapply with a data file with mulitple communities.
#Marc Cadotte
#Nov 2008

Krand<- function (phy,com,n.rand=100,sub=FALSE,fig=TRUE) {
	
	if (sub==TRUE) {
		source("/Volumes/cadotte/Documents/R/phylogeny code/subTree.R")
		phy<-subTree(phy,com[,3])
		}
		
	s.order <- match(phy$tip.label,com[,3])
	trait <- com[,2][s.order]
	
	#variable has to be a matrix
	trait <- matrix(log(trait),ncol=1)
		
	K.obs<-Kcalc(trait,phy)
	K.obs<-data.frame(K.obs)	
	
	kvalues<-numeric(length(n.rand))
	
	for (i in 1:n.rand) {
 		kvalues[i]<-Kcalc(sample(trait),phy)
		}

	#null stats
	K.null<-kvalues
	K.null.mean<-mean(K.null)
	K.null.95<-quantile(K.null,probs=c(0.025,0.975))
	names(K.null.95)<-NULL
	
	results<-cbind(K.obs,K.null.mean,
		K.null.lwr=K.null.95[1],
		K.null.upr=K.null.95[2],
		n.rand)
					
	if (fig==TRUE) {
		par(cex.lab=1.5, cex.axis=1.2, cex.main=1.5)
		hist(K.null, xlab="K values",
		main="Null K values; dahsed line = K obs; dotted = 95% CI")
		abline(v=K.obs, lty="dashed")
		abline(v=K.null.95[1],lty="dotted")
		abline(v=K.null.95[2],lty="dotted")
		}	
	
	return(results)	
	}