/packages/archive/2010/05.2010/05.01.2010/spacodi/R/K.by.nodes.R

  1K.by.nodes<-function(sp.traits, phy, obs.only=FALSE, return.all=TRUE, n.rep=10, rand.test=TRUE, r.rep=10000) {
  2	
  3# parameter settings
  4	if(is.null(sp.traits) && is.null(phy)) {
  5		stop("Must supply both sp.traits and tree.")
  6	}
  7	if(is.null(r.rep) && rand.test)stop("Must specify 'r.rep' if using the randomization test.")
  8	
  9	if(obs.only==TRUE){
 10		rand.test=FALSE 
 11		exp.K=FALSE
 12	} else {exp.K=TRUE}
 13	
 14# internal function calling Kcalc()	
 15	K.all.nodes <-function(sp.traits, phy, return.all=TRUE){
 16		if(ncol(sp.traits)!=1)stop("Cannot handle more than a single trait at a time.")
 17		if (length(phy$tip.label)!=nrow(sp.traits)) stop("Tree cannot be matched to data.")
 18		tt=as.data.frame(as.numeric(as.vector(sp.traits[,1])))
 19		row.names(tt)=row.names(sp.traits)
 20		names(tt)="trait"
 21		sp.traits=tt
 22		sub=subtrees(phy)
 23		if(is.null(rownames(sp.traits))||!all(row.names(sp.traits)%in%phy$tip.label)) stop("Trait matrix must have row names corresponding to tip labels of the phylogeny.")
 24		n.nodes <- phy$Nnode
 25		out <- array(dim=c(length(sub), 4)) ### prep an output file
 26		for (i in 1:n.nodes) {### cycle over subtrees
 27			trt.i <- data.frame(sp.traits[match(sub[[i]]$tip.label, rownames(sp.traits)),]) #extract the relevant individual from trait matrix, put them into tiplabels order. 
 28			nn<-sub[[i]]$Nnode
 29			if(length(trt.i) > 0 & var(trt.i) != 0 & nn > 1){ # eliminate basal polytomies and cases with no trait variance
 30				K.i <- Kcalc(trt.i, sub[[i]], FALSE)
 31			} else {
 32				K.i <- NA	
 33			}
 34			out[i, 1] <- K.i
 35			out[i, 2] <- nn+1
 36			out[i, 3] <- max(branching.times(sub[[i]]))
 37			out[i, 4] <- min(sub[[i]]$node)
 38		}
 39		array.cols=c("Blomberg.K","tips","node.time","node.ID")
 40		dimnames(out)=list(NULL, array.cols)
 41		if(!return.all) {
 42			out=out[which(!is.na(out[,1])),]
 43		}
 44		return(out)
 45	}
 46	
 47# internal function preparing permuted datasets for K.all.nodes()
 48	K.exp.nodes <- function(phy, n.rep = 10, Brownian=TRUE,...) {
 49		out <- array(dim=c(phy$Nnode, 4, n.rep)) ### prep an output file
 50	
 51		if(Brownian==TRUE){
 52			trt.sims <- as.data.frame(sim.char(phy=phy, model.matrix=as.matrix(1), nsims=n.rep, model="brownian"))
 53			names(trt.sims)=paste("trt",seq(1:n.rep),sep=".")
 54		} else {
 55			trt.sims <- as.data.frame(replicate(n.rep, runif(length(phy$tip.label))))
 56			row.names(trt.sims)=phy$tip.label
 57			names(trt.sims)=paste("trt",seq(1:n.rep),sep=".")
 58		}
 59		for(k in 1:n.rep) {
 60			k.sim=as.data.frame(trt.sims[,k])
 61			row.names(k.sim)=row.names(trt.sims)
 62			out[,,k] <- foo <- K.all.nodes(phy=phy, sp.traits=k.sim, return.all=TRUE)
 63		}	
 64		dimnames(out)=list(NULL, names(as.data.frame(foo)), paste("iter",seq(1:n.rep),sep="."))
 65		return(out)
 66	}
 67	
 68# get K for nodes							
 69	o.foo=K.all.nodes(sp.traits=sp.traits, phy=phy, return.all=ifelse(obs.only==TRUE, TRUE, FALSE))
 70	o=o.foo[,"Blomberg.K"]
 71	names(o)=o.foo[,"node.ID"]
 72	o=o[order(o)]
 73	n.o=names(o)
 74	
 75	o.orig=as.data.frame(o.foo)
 76	row.names(o.orig)=o.orig$node.ID
 77	
 78	
 79	if(exp.K==TRUE) {	
 80		
 81# K permutation								
 82		e.out=array(dim=c(length(o), n.rep))		
 83		for(bb in 1:n.rep) {
 84			b.foo=as.data.frame(K.exp.nodes(sp.traits=sp.traits, phy=phy, n.rep=1)[,,1])
 85			b.match=match(as.numeric(n.o), b.foo$node.ID)
 86			e.out[,bb]=as.numeric(b.foo[b.match,"Blomberg.K"])
 87		}
 88		
 89# dataset matching
 90		e.out=as.data.frame(e.out)
 91		row.names(e.out)=as.numeric(n.o)
 92		names(e.out)=paste("iter",seq(1:n.rep),sep="")
 93		match(as.numeric(o.orig$node.ID),as.numeric(row.names(e.out)))->orig.match
 94		e.out=e.out[orig.match,]
 95		
 96# randomization test
 97		if(rand.test==TRUE) {
 98			rand.array=array(dim=c(nrow(o.orig),5))
 99			for(node in 1:nrow(o.orig)) {
100				obs=o.orig$Blomberg.K[node]
101				nn=o.orig$node.ID[node]
102				exp=e.out[node,which(!is.na(e.out[node,]))]
103				rand.array[node,1]=ifelse(length(exp)>0, randomization.test.sp(obs=obs,exp=exp, iter=r.rep, return.all=FALSE, two.tailed=TRUE), NA)
104				rand.array[node,2]=nn
105				rand.array[node,3]=obs
106				rand.array[node,4]=ifelse(length(exp)>0, mean(exp), NA)
107				rand.array[node,5]=length(exp)
108			}
109			r.test=as.data.frame(rand.array)
110			names(r.test)=c("p.value","node.ID","obs.K","m.exp.K","valid.comparisons")
111			row.names(r.test)=r.test$node.ID
112			if(return.all==TRUE) {
113				o.orig=as.data.frame(K.all.nodes(sp.traits=sp.traits, phy=phy, return.all=TRUE))
114				row.names(o.orig)=o.orig$node.ID
115			} else {o.orig=o.orig}
116			K.out=list(o.orig, e.out, r.test)
117			names(K.out)=c("observed.K","expected.K","randomization.test")
118			return(K.out)
119		} else {
120			if(return.all==TRUE) {
121				o.orig=as.data.frame(K.all.nodes(sp.traits=sp.traits, phy=phy, return.all=TRUE))
122				row.names(o.orig)=o.orig$node.ID
123			} else {o.orig=o.orig}
124			K.out=list(o.orig, e.out)
125			names(K.out)=c("observed.K","expected.K")
126			return(K.out)
127		}
128	} else {
129		K.out=list(o.orig)
130		names(K.out)="observed.K"
131		return(K.out)
132	}
133}
134
135