/spacodiR/R/spacodi.calc.R
R | 152 lines | 128 code | 16 blank | 8 comment | 43 complexity | 4ae3a5c75c1be88ba2caf75608d90a14 MD5 | raw file
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