/spacodiR/R/spacodi.calc.R

  1spacodi.calc <-function(sp.plot, phy = NULL, sp.traits = NULL, all.together=TRUE, prune=TRUE, pairwise=FALSE,...){
  2	if(!missing(phy) && !missing(sp.traits)) stop("Please supply either a phylogeny or a trait matrix, but not both.")
  3
  4	rematrix=function(arr, names){
  5		n=sqrt(length(arr))
  6		counter=1
  7		M=matrix(0,n,n)
  8		for(i in 1:n){
  9			j=i+1
 10			while(j<=n) {
 11				M[i,j]=arr[counter]
 12				counter=counter+1
 13				j=j+1
 14			}
 15		}
 16		M[lower.tri(M)]=t(M)[lower.tri(t(M))]
 17		dimnames(M)[[1]]<-dimnames(M)[[2]]<-names
 18		return(M)
 19	}
 20	
 21	# determine type of abundance
 22	stripped=unname(unlist(c(sp.plot)))
 23	if(all(!is.na(match(stripped, c(0,1))))) {
 24		abundt = 0		# presence|absence
 25	} else if(all(!is.na(match(range(stripped), c(0,1)))) && length(unique(stripped))>2) {
 26		abundt = 1		# relative abundance
 27	} else {
 28		abundt = 2		# abundance (n.individuals)
 29	}
 30	
 31	# iter is 1 unless more than a single trait is being used for separate analyses
 32	iter=1
 33	traits.tmp<-distmat<-NULL
 34	
 35	# INTERPRET whether to compute trait diversity or phylogenetic diversity & prepare data
 36	if(!missing(phy)) {
 37		# distmat is phylogenetic: Pst
 38		sp.data=match.spacodi.data(sp.plot=sp.plot, phy=phy, prune=prune, ...)
 39		sp.plot=sp.data$sp.plot
 40		phy=sp.data$sp.tree
 41		if(check.distmat(phy)) distmat <- phy else distmat <- cophen(phy)
 42		
 43	} else if(missing(sp.traits) & missing(phy)) {
 44		# distmat is null: Ist
 45		distmat <- matrix(1, ncol = nrow(sp.plot), nrow = nrow(sp.plot))
 46		
 47	} else if(!missing(sp.traits)){
 48		# distmat is trait-based: Tst
 49		if(class(sp.traits)=="dist") sp.traits=as.matrix(sp.traits)
 50		if(ncol(sp.traits)==1) all.together=TRUE
 51		if(all(is.null(names(sp.traits)))) names(sp.traits)=paste("trt",seq(1:ncol(sp.traits)),sep="")	
 52		if(all.together==TRUE | check.distmat(sp.traits)) {
 53			sp.data=match.spacodi.data(sp.plot=sp.plot, sp.traits=sp.traits, prune=prune, ...)
 54			sp.plot=sp.data$sp.plot
 55			sp.traits=sp.data$sp.traits
 56			if(check.distmat(sp.traits)) distmat=sp.traits else distmat=as.matrix(dist(sp.traits))
 57		} else {	
 58			iter=ncol(sp.traits)
 59			traits.tmp=lapply(1:ncol(sp.traits), function(x) {
 60				   trait.tt<-data.frame(sp.traits[,x])
 61				   row.names(trait.tt)<-row.names(sp.traits)
 62				   sp.data<-match.spacodi.data(sp.plot=sp.plot, sp.traits=trait.tt, prune=prune, ...)
 63				   distmat <- as.matrix(dist(sp.data$sp.traits))
 64				   return(list(distmat=distmat, sp.plot=sp.data$sp.plot))
 65				   }
 66				   )
 67		}
 68	} else if(is.null(distmat)){ 
 69		stop("Cannot decipher input object(s).")
 70	}
 71	
 72	if(is.null(names(sp.plot))) pnames=paste("plt",seq(1,ncol(sp.plot)),sep=".") else pnames=names(sp.plot)
 73
 74	
 75	# PREPARE output
 76	gen.out=list()
 77	prw.out=list()
 78	
 79	for(tt in 1:iter) {	
 80		if(is.null(traits.tmp)) {
 81			sp.plot <-	as.matrix(sp.plot)
 82		} else {
 83			sp.plot <-	as.matrix(traits.tmp[[tt]]$sp.plot)
 84			distmat <-	as.matrix(traits.tmp[[tt]]$distmat)
 85		}
 86		
 87		diag(distmat) <- 0
 88		np <- ncol(sp.plot)
 89		out <- NA
 90		out <- .C("spacodi", 
 91				  np = as.integer(np),
 92				  ns = as.integer(nrow(sp.plot)),
 93				  sp.plot = as.double(as.vector(sp.plot)),
 94				  distmat = as.double(as.vector(as.matrix(distmat))),
 95				  abundtype = as.integer(abundt),
 96				  Ndclass = as.integer(0),
 97				  dclass = as.double(c(0,0)),
 98				  Ist = 0,
 99				  Pst = 0,
100				  Bst = 0,
101				  PIst = 0,
102				  pairwiseIst = as.double(as.vector(matrix(0,np,np))),
103				  pairwisePst = as.double(as.vector(matrix(0,np,np))),
104				  pairwiseBst = as.double(as.vector(matrix(0,np,np))),
105				  pairwisePIst = as.double(as.vector(matrix(0,np,np))),
106				  PACKAGE="spacodiR"
107				  )
108		
109		# compile results for phylogenetic turnover
110		if(missing(phy) & missing(sp.traits)) {
111			r.out=as.numeric(out[8])
112			names(r.out)="Ist"
113			gen.out[[tt]]=as.data.frame(t(r.out))
114			if(pairwise) {
115				prw.out[[tt]]=list(rematrix(unlist(out[12]),pnames))
116				names(prw.out[[tt]])<-paste("pairwise","Ist",sep=".")
117			}
118		} else if(!missing(phy)){
119			r.out=as.numeric(c(out[8:11]))
120			names(r.out)<-c("Ist","Pst","Bst","PIst")
121			gen.out[[tt]]=as.data.frame(t(r.out))
122			if(pairwise) {
123				prw.out[[tt]]=lapply(out[12:15], function(x) rematrix(x,pnames))
124				names(prw.out[[tt]])<-paste("pairwise",c("Ist","Pst","Bst","PIst"),sep=".")
125			}
126		} else if(!missing(sp.traits)) {
127			r.out=as.numeric(c(out[8:11]))
128			names(r.out)=c("Ist","Tst","Ust","TAUst")
129			gen.out[[tt]]=as.data.frame(t(r.out))
130			if(pairwise) {
131				prw.out[[tt]]=lapply(out[12:15], function(x) rematrix(x,pnames))
132				names(prw.out[[tt]])<-paste("pairwise",c("Ist","Tst","Ust","TAUst"),sep=".")
133			}
134		} 
135	} 
136	if(iter>1 & !all.together) {
137		RES=lapply(1:iter, function(x) {
138				   if(pairwise) return(c(gen.out[[x]], prw.out[[x]])) else return(c(gen.out[[x]]))
139				   }
140				   )
141		names(RES)<-names(sp.traits)
142	} else {
143		prw.out=unlist(prw.out,recursive=FALSE)
144		gen.out=unlist(gen.out,recursive=FALSE)
145		if(pairwise) RES=c(gen.out,prw.out) else RES=gen.out
146	}
147	
148	return(unlist(list(RES,sp.plot=list(sp.plot),sp.tree=list(phy),sp.traits=list(sp.traits)),recursive=FALSE)) 
149
150} 
151
152