Naiive Diversity Profiles + Rarefaction

NaiiveD

require(picante)
require(gtools)
require(reshape)




###########################################################
### For calculation of Naiive (Z = I)  diversity profiles ###
###########################################################

NaiiveD = function(comm, phy){
  #Begin
  samp1 <- t(comm)
  samp1 <- samp1[mixedsort(rownames(samp1)),]
  taxa = nrow(samp1)
  rownames(samp1) <- rep(1:taxa)
  sums <- subset(rowSums(samp1), rowSums(samp1) > 0)
  samp1 <- samp1[names(sums),]
  samples = ncol(samp1)
  taxa = nrow(samp1)

  p=matrix(0,taxa,samples)
  
  for (k in 1:samples){
    p[,k]<-samp1[,k]/sum(samp1[,k])
  }
  
  Z <- diag(taxa)
    
  lenq = 50
  qq <- seq(length=lenq, from=0, by=.11)
  
  Zp=matrix(0,taxa,samples)
  
  for (k in 1:samples) {
    for (i in 1:taxa){
      for (j in 1:taxa){
        Zp[i,k]<-Zp[i,k]+Z[i,j]*p[j,k]
      }
    }}
  
  Dqz = matrix(0, lenq ,samples)  
  
  for (k in 1:samples) {
    
    for (iq in 1:lenq)  {
      q<-qq[iq];
      
      
      for (zpi in 1:length(Zp[,k])){
        if (Zp[zpi,k]>0)(
          Dqz[iq,k]<-Dqz[iq,k]+ p[zpi,k]*(Zp[zpi,k])^(q-1))
      }
      
      Dqz[iq,k] <- Dqz[iq,k]^(1/(1-q));
    }}
  results = Dqz
  return(results)
}

###########################################################
### ChaoRarefy- diversity profiles rarefied ###
###########################################################

#This function takes your sample and phylogeny and rarefies to account for differences in sample sizes
# Y is usually the smallest sample in the community = min(rowSums(samp))
# X is the number of sampled communities to average the results over (should be >100)

NaiiveRarefy=function (samp,Y,K){
  Res<-matrix(NA,ncol=3,nrow=dim(samp)[1])
  rowsums<-rowSums(samp)
  df.list<-vector("list",K)
  
  for(a in 1:K){
    print(a)
    samplematrix1<-samp
    rare_sample<-rrarefy(samplematrix1,Y)
    rare_sample = rare_sample[,colSums(rare_sample)>0]
    dqz<-NaiiveD(rare_sample)
    df.list[[a]]<-dqz
    
  } # close for loop
  
  
  
  Res<-Reduce("+",df.list)/length(df.list)
  
  
  return(Res)
  
}


###########################################################
### An example ###
###########################################################
data(phylocom)
samp = phylocom$sample
phy = phylocom$phylo

NaiiveD(samp)
NaiiveRarefy(samp = samp, phy = phy, Y = min(rowSums(samp)), K = 100)