/spacodiR/R/spacodi.calc.R

http://github.com/eastman/spacodiR · R · 152 lines · 128 code · 16 blank · 8 comment · 43 complexity · 4ae3a5c75c1be88ba2caf75608d90a14 MD5 · raw file

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