# simulation based test for a correlation between the state of x & the rate of y
# written by Liam J. Revell 2013

ratebystate<-function(tree,x,y,nsim=100,corr=c("pearson","spearman"),...){
	corr<-corr[1]
	if(hasArg(sim.method)) sim.method<-list(...)$sim.method
	else sim.method<-"sim.corrs"
	if(hasArg(method)) method<-list(...)$method
	else method<-"by.node"
	if(hasArg(message)) message<-list(...)$message
	else message<-TRUE
	if(!is.binary.tree(tree)) tree<-multi2di(tree)
	V<-phyl.vcv(cbind(x,y),vcv(tree),lambda=1)$R
	if(method=="by.branch"){
		aa<-c(x[tree$tip.label],fastAnc(tree,x))
		names(aa)[1:length(tree$tip)]<-1:length(tree$tip)
		aa<-rowMeans(matrix(aa[tree$edge],nrow(tree$edge),2))
		a<-vector()
		for(i in 1:tree$Nnode+length(tree$tip)){
			j<-which(tree$edge[,1]==i)
			a[i-length(tree$tip)]<-sum(aa[j]*tree$edge.length[j])/sum(tree$edge.length[j])
		}
		names(a)<-1:tree$Nnode+length(tree$tip)
	}	
	else a<-fastAnc(tree,x)
	b<-pic(y,tree)[names(a)]^2
	r<-cor(a,b,method=corr)
	beta<-setNames(lm(b~a)$coefficients[2],NULL)
	foo<-function(tree,V){
		if(sim.method=="fastBM") XY<-fastBM(tree,nsim=2)%*%sqrt(diag(diag(V)))
		else if(sim.method=="sim.corrs") XY<-sim.corrs(tree,V)
		a<-fastAnc(tree,XY[,1])
		b<-pic(XY[,2],tree)[names(a)]^2
		r<-cor(a,b,method=corr)
		return(r)
	}
	for(i in 1:99) foo(tree,V)
	r.null<-c(r,replicate(nsim-1,foo(tree,V)))
	P<-mean(abs(r.null)>=abs(r))
	if(message) return(list(beta=beta,r=r,P=P,corr=corr,method=method))
	else  return(list(beta=beta,r=r,P=P))
}