#########################################################################
# Function part.calc provides alpha, beta, and gamma components of 
# species richness and shannon diversity by habitat types. Both additive 
# and multiplicative components of beta diversity are calculated for 
# species richness and shannon diversity
#########################################################################

part.calc<-function(specsamp.mat,totspec,totabun,nsamp,nhab) {

# sample statistics for abundance, richness, and diversity
richcount<-function(x) length(x[x>0])
sampabun<-apply(specsamp.mat,MARGIN=1,sum)
samprich<-apply(specsamp.mat,MARGIN=1,richcount)
specabun<-apply(specsamp.mat,MARGIN=2,sum)
propabun.mat<-specsamp.mat/sampabun
shannon.mat<-propabun.mat*(log(propabun.mat))
shannon.mat[shannon.mat=="NaN"]<-0
sampshan<-(-(apply(shannon.mat,MARGIN=1,sum)))
sampshan.mult<-exp(sampshan)

# habitat vectors and overall measures of shannon diversity
totprop.vec<-specabun/totabun
totshan.vec<-totprop.vec*(log(totprop.vec))
totshan.vec[totshan.vec=="NaN"]<-0
shan.tot<-(-sum(totshan.vec))
alpha.samp<-mean(samprich)
hab.vec<-as.factor(specsamp$habitat)
hab.nsamp<-table(specsamp$habitat)

# alpha diversity measures by habitat 
alpha.hab<-tapply(samprich,hab.vec,mean)
abund.hab<-tapply(sampabun,hab.vec,mean)
shan.hab<-tapply(sampshan,hab.vec,mean)
shan.mult.hab<-exp(tapply(sampshan,hab.vec,mean))

# abundance and gamma diversity by habitat
totabun.hab<-tapply(sampabun,hab.vec,sum)
hab.mat<-rowsum(specsamp.mat,hab.vec)
gamma.hab<-apply(hab.mat,MARGIN=1,richcount)
habprop.mat<-hab.mat
for (i in 1:totspec) habprop.mat[1:nhab,i]<-hab.mat[,i]/totabun.hab[1:nhab]
shan.hab.mat<-habprop.mat*(log(habprop.mat))
shan.hab.mat[shan.hab.mat=="NaN"]<-0
gamma.shan.hab<-(-(apply(shan.hab.mat,MARGIN=1,sum)))
gamma.shan.mult.hab<-exp(gamma.shan.hab)

# beta diversity by habitat
beta.hab<-gamma.hab-alpha.hab
beta.mult.hab<-gamma.hab/alpha.hab
beta.shan.hab<-gamma.shan.hab-shan.hab
beta.shan.mult.hab<-gamma.shan.mult.hab/shan.mult.hab

# set up arrays for F-statistics
beta.samp<-array(nsamp)
beta.shan.samp<-array(nsamp)
beta.mult.samp<-array(nsamp)
beta.shan.mult.samp<-array(nsamp)
for (n in 1:nsamp) {
  beta.samp[n]<-gamma.hab[hab.vec[n]]-samprich[n]
  beta.shan.samp[n]<-gamma.shan.hab[hab.vec[n]]-sampshan[n]
  beta.mult.samp[n]<-gamma.hab[hab.vec[n]]/samprich[n]
  beta.shan.mult.samp[n]<-gamma.shan.mult.hab[hab.vec[n]]/sampshan.mult[n]
}

# F-statistics for alpha and beta
alpha.Fstat<-summary(lm(samprich~hab.vec))$fstatistic[1]
beta.Fstat<-summary(lm(beta.samp~hab.vec))$fstatistic[1]
beta.mult.Fstat<-summary(lm(beta.mult.samp~hab.vec))$fstatistic[1]
shan.alpha.Fstat<-summary(lm(sampshan~hab.vec))$fstatistic[1]
shan.beta.Fstat<-summary(lm(beta.shan.samp~hab.vec))$fstatistic[1]
shan.mult.alpha.Fstat<-summary(lm(sampshan.mult~hab.vec))$fstatistic[1]
shan.mult.beta.Fstat<-summary(lm(beta.shan.mult.samp~hab.vec))$fstatistic[1]

# Pairwise differences among means
pair.alpha<-array(dim=c(nhab,nhab))
pair.beta<-array(dim=c(nhab,nhab))
pair.mult.beta<-array(dim=c(nhab,nhab))
pair.shan.alpha<-array(dim=c(nhab,nhab))
pair.shan.beta<-array(dim=c(nhab,nhab))
pair.shan.mult.alpha<-array(dim=c(nhab,nhab))
pair.shan.mult.beta<-array(dim=c(nhab,nhab))
pair.gamma<-array(dim=c(nhab,nhab))
pair.shan.gamma<-array(dim=c(nhab,nhab))
pair.shan.mult.gamma<-array(dim=c(nhab,nhab))
for (j in 1:nhab) {
  for (i in 1:nhab) {
    pair.alpha[i,j]<-alpha.hab[i]-alpha.hab[j]
    pair.beta[i,j]<-beta.hab[i]-beta.hab[j]
    pair.mult.beta[i,j]<-beta.mult.hab[i]-beta.mult.hab[j]
    pair.shan.alpha[i,j]<-shan.hab[i]-shan.hab[j]
    pair.shan.beta[i,j]<-beta.shan.hab[i]-beta.shan.hab[j]
    pair.shan.mult.alpha[i,j]<-shan.mult.hab[i]-shan.mult.hab[j]
    pair.shan.mult.beta[i,j]<-beta.shan.mult.hab[i]-beta.shan.mult.hab[j]
    pair.gamma[i,j]<-gamma.hab[i]-gamma.hab[j]
    pair.shan.gamma[i,j]<-gamma.shan.hab[i]-gamma.shan.hab[j]
    pair.shan.mult.gamma[i,j]<-gamma.shan.mult.hab[i]-gamma.shan.mult.hab[j]
  }
}

#Summarize results in data frames, combine and return as a list
results.overall<-data.frame(TotalSamples=nsamp,TotalAbundance=totabun,
    TotalSpecies=totspec)
results.part<-data.frame(Samples=as.vector(hab.nsamp),Abundance=totabun.hab,
    Alpha=alpha.hab,BetaAdd=beta.hab,BetaMult=beta.mult.hab,Gamma=gamma.hab,
    Shannon=shan.hab,BetaShanAdd=beta.shan.hab,
    AlphaShanMult=shan.mult.hab,BetaShanMult=beta.shan.mult.hab,
    GammaShanAdd=gamma.shan.hab,GammaShanMult=gamma.shan.mult.hab)
results.fstat<-data.frame(AlphaFstat=alpha.Fstat,BetaFstatAdd=beta.Fstat,
    BetaFstatMult=beta.mult.Fstat,ShanAlphaFstat=shan.alpha.Fstat,
    ShanBetaFstat=shan.beta.Fstat,
    ShanAlphaFstatMult=shan.mult.alpha.Fstat,
    ShanBetaFstatMult=shan.mult.beta.Fstat)
results.pair<-list(AlphaDif=pair.alpha,BetaAddDif=pair.beta,
    BetaMultDif=pair.mult.beta,ShanAlphaDif=pair.shan.alpha,
    ShanBetaDif=pair.shan.beta,ShanAlphaMultDif=pair.shan.mult.alpha,
    ShanBetaMultDif=pair.shan.mult.beta,GammaDif=pair.gamma,
    ShanGammaDif=pair.shan.gamma,ShanMultGammaDif=pair.shan.mult.gamma)
results.comb<-list(Overall=results.overall,Partitions=results.part,
    Fstats=results.fstat,PairDiff=results.pair)
return(results.comb)
}
# End of function

############################################################################

# read data file: column 1 is sample id, column 2 is habitat id,
# and columns 3:totspec are species id

setwd()
specsamp<-read.table("data.txt",header=TRUE)

# Determine number of species, number of samples, and overall abundance
ncolm<-length(names(specsamp))
totspec<-ncolm-2
nsamp<-length(specsamp$trap)
nhab<-length(levels(as.factor(specsamp$habitat)))
hab.names<-levels(as.factor(specsamp$habitat))
specsamp.mat<-as.matrix(specsamp[,3:ncolm])
totabun<-sum(specsamp.mat)

# Call function to calculate observed values
obs.results<-part.calc(specsamp.mat,totspec,totabun,nsamp,nhab)

# Set number of randomizations
krand<-1000

# Initialize data arrays
specsamp.mat.rand=array(0,dim=c(nsamp,totspec))
nvec<-1:nsamp
abund.null<-array(0,dim=c(nhab,krand))
alpha.null<-array(0,dim=c(nhab,krand))
beta.null<-array(0,dim=c(nhab,krand))
beta.mult.null<-array(0,dim=c(nhab,krand))
gamma.null<-array(0,dim=c(nhab,krand))
alpha.shan.null<-array(0,dim=c(nhab,krand))
beta.shan.null<-array(0,dim=c(nhab,krand))
gamma.shan.null<-array(0,dim=c(nhab,krand))
alpha.mult.shan.null<-array(0,dim=c(nhab,krand))
beta.mult.shan.null<-array(0,dim=c(nhab,krand))
gamma.mult.shan.null<-array(0,dim=c(nhab,krand))
alpha.fstat.null<-array(0,krand)
beta.fstat.null<-array(0,krand)
beta.mult.fstat.null<-array(0,krand)
alpha.shan.fstat.null<-array(0,krand)
beta.shan.fstat.null<-array(0,krand)
alpha.mult.shan.fstat.null<-array(0,krand)
beta.mult.shan.fstat.null<-array(0,krand)
pair.alpha.null<-array(0,dim=c(nhab,nhab,krand))
pair.beta.null<-array(0,dim=c(nhab,nhab,krand))
pair.beta.mult.null<-array(0,dim=c(nhab,nhab,krand))
pair.shan.alpha.null<-array(0,dim=c(nhab,nhab,krand))
pair.shan.beta.null<-array(0,dim=c(nhab,nhab,krand))
pair.shan.alpha.mult.null<-array(0,dim=c(nhab,nhab,krand))
pair.shan.beta.mult.null<-array(0,dim=c(nhab,nhab,krand))
pair.gamma.null<-array(0,dim=c(nhab,nhab,krand))
pair.gamma.shan.null<-array(0,dim=c(nhab,nhab,krand))
pair.gamma.shan.mult.null<-array(0,dim=c(nhab,nhab,krand))

# Randomization Loop - calls part.calc function each step of loop
for (k in 1:krand) {
  rvec<-sample(1:nsamp,nsamp,replace=FALSE)
  specsamp.mat.rand[nvec,]<-specsamp.mat[rvec,]
  rand.results<-part.calc(specsamp.mat.rand,totspec,totabun,nsamp,nhab)
  abund.null[,k]<-rand.results$Partitions$Abundance
  alpha.null[,k]<-rand.results$Partitions$Alpha
  beta.null[,k]<-rand.results$Partitions$BetaAdd
  gamma.null[,k]<-rand.results$Partitions$Gamma
  beta.mult.null[,k]<-rand.results$Partitions$BetaMult
  alpha.shan.null[,k]<-rand.results$Partitions$Shannon
  beta.shan.null[,k]<-rand.results$Partitions$BetaShanAdd
  gamma.shan.null[,k]<-rand.results$Partitions$GammaShanAdd
  alpha.mult.shan.null[,k]<-rand.results$Partitions$AlphaShanMult
  beta.mult.shan.null[,k]<-rand.results$Partitions$BetaShanMult
  gamma.mult.shan.null[,k]<-rand.results$Partitions$GammaShanMult
  alpha.fstat.null[k]<-rand.results$Fstats$AlphaFstat
  beta.fstat.null[k]<-rand.results$Fstats$BetaFstatAdd
  beta.mult.fstat.null[k]<-rand.results$Fstats$BetaFstatMult
  alpha.shan.fstat.null[k]<-rand.results$Fstats$ShanAlphaFstat
  beta.shan.fstat.null[k]<-rand.results$Fstats$ShanBetaFstat
  alpha.mult.shan.fstat.null[k]<-rand.results$Fstats$ShanAlphaFstatMult
  beta.mult.shan.fstat.null[k]<-rand.results$Fstats$ShanBetaFstatMult
  pair.alpha.null[,,k]<-rand.results$PairDiff$AlphaDif
  pair.beta.null[,,k]<-rand.results$PairDiff$BetaAddDif
  pair.beta.mult.null[,,k]<-rand.results$PairDiff$BetaMultDif
  pair.shan.alpha.null[,,k]<-rand.results$PairDiff$ShanAlphaDif
  pair.shan.beta.null[,,k]<-rand.results$PairDiff$ShanBetaDif
  pair.shan.alpha.mult.null[,,k]<-rand.results$PairDiff$ShanAlphaMultDif
  pair.shan.beta.mult.null[,,k]<-rand.results$PairDiff$ShanBetaMultDif
  pair.gamma.null[,,k]<-rand.results$PairDiff$GammaDif
  pair.gamma.shan.null[,,k]<-rand.results$PairDiff$ShanGammaDif
  pair.gamma.shan.mult.null[,,k]<-rand.results$PairDiff$ShanMultGammaDif
}
# End Randomization Loop

# Initialize arrays to calculate p-values and quantiles by habitats
pval.abund=array(nhab)
pval.alpha=array(nhab)
pval.beta=array(nhab)
pval.beta.mult=array(nhab)
pval.gamma=array(nhab)
pval.alpha.shan=array(nhab)
pval.beta.shan=array(nhab)
pval.alpha.shan.mult=array(nhab)
pval.beta.shan.mult=array(nhab)
pval.gamma.shan=array(nhab)
quant.abund=array(0,dim=c(nhab,5))
quant.alpha=array(0,dim=c(nhab,5))
quant.beta=array(0,dim=c(nhab,5))
quant.beta.mult=array(0,dim=c(nhab,5))
quant.gamma=array(0,dim=c(nhab,5))
quant.alpha.shan=array(0,dim=c(nhab,5))
quant.beta.shan=array(0,dim=c(nhab,5))
quant.alpha.shan.mult=array(0,dim=c(nhab,5))
quant.beta.shan.mult=array(0,dim=c(nhab,5))
quant.gamma.shan=array(0,dim=c(nhab,5))

# Calculate P-values for within-habitat abundance and diversity and the 
# 95% quantiles of the null distributions
for (h in 1:nhab) {
  pval.abund[h]<-length(abund.null[h,abund.null[h,]>
    obs.results$Partitions$Abundance[h]])/krand
  pval.alpha[h]<-length(alpha.null[h,alpha.null[h,]>
    obs.results$Partitions$Alpha[h]])/krand
  pval.beta[h]<-length(beta.null[h,beta.null[h,]>
    obs.results$Partitions$BetaAdd[h]])/krand
  pval.beta.mult[h]<-length(beta.mult.null[h,beta.mult.null[h,]>
    obs.results$Partitions$BetaMult[h]])/krand
  pval.gamma[h]<-length(gamma.null[h,gamma.null[h,]>
    obs.results$Partitions$Gamma[h]])/krand
  pval.alpha.shan[h]<-length(alpha.shan.null[h,alpha.shan.null[h,]>
    obs.results$Partitions$Shannon[h]])/krand
  pval.beta.shan[h]<-length(beta.shan.null[h,beta.shan.null[h,]>
    obs.results$Partitions$BetaShanAdd[h]])/krand
  pval.gamma.shan[h]<-length(gamma.shan.null[h,gamma.shan.null[h,]>
    obs.results$Partitions$GammaShanAdd[h]])/krand
  pval.alpha.shan.mult[h]<-
    length(alpha.mult.shan.null[h,alpha.mult.shan.null[h,]>
    obs.results$Partitions$AlphaShanMult[h]])/krand
  pval.beta.shan.mult[h]<-
    length(beta.mult.shan.null[h,beta.mult.shan.null[h,]>
    obs.results$Partitions$BetaShanMult[h]])/krand
  quant.abund[h,]<-quantile(abund.null[h,],probs=c(0,0.025,0.5,0.975,1))
  quant.alpha[h,]<-quantile(alpha.null[h,],probs=c(0,0.025,0.5,0.975,1))
  quant.beta[h,]<-quantile(beta.null[h,],probs=c(0,0.025,0.5,0.975,1))
  quant.gamma[h,]<-quantile(gamma.null[h,],probs=c(0,0.025,0.5,0.975,1))
  quant.alpha.shan[h,]<-
    quantile(alpha.shan.null[h,],probs=c(0,0.025,0.5,0.975,1))
  quant.beta.shan[h,]<-
    quantile(beta.shan.null[h,],probs=c(0,0.025,0.5,0.975,1))
  quant.gamma.shan[h,]<-
    quantile(gamma.shan.null[h,],probs=c(0,0.025,0.5,0.975,1))
}

# Calculate the P-values for the F-statistics
pval.alpha.fstat<-length(alpha.fstat.null[alpha.fstat.null[]>
  obs.results$Fstats$AlphaFstat])/krand
pval.beta.fstat<-length(beta.fstat.null[beta.fstat.null[]>
  obs.results$Fstats$BetaFstatAdd])/krand
pval.beta.mult.fstat<-length(beta.mult.fstat.null[beta.mult.fstat.null[]>
  obs.results$Fstats$BetaFstatMult])/krand
pval.alpha.shan.fstat<-length(alpha.shan.fstat.null[alpha.shan.fstat.null[]>
  obs.results$Fstats$ShanAlphaFstat])/krand
pval.beta.shan.fstat<-length(beta.shan.fstat.null[beta.shan.fstat.null[]>
  obs.results$Fstats$ShanBetaFstat])/krand
pval.alpha.shan.mult.fstat<-
  length(alpha.mult.shan.fstat.null[alpha.mult.shan.fstat.null[]>
  obs.results$Fstats$ShanAlphaFstatMult])/krand
pval.beta.shan.mult.fstat<-
  length(beta.mult.shan.fstat.null[beta.mult.shan.fstat.null[]>
  obs.results$Fstats$ShanBetaFstatMult])/krand

# Set up arrays to get p-values for pairwise comparisons
pval.pair.alpha.null<-array(dim=c(nhab,nhab))
pval.pair.beta.null<-array(dim=c(nhab,nhab))
pval.pair.beta.mult.null<-array(dim=c(nhab,nhab))
pval.pair.shan.alpha.null<-array(dim=c(nhab,nhab))
pval.pair.shan.beta.null<-array(dim=c(nhab,nhab))
pval.pair.shan.alpha.mult.null<-array(dim=c(nhab,nhab))
pval.pair.shan.beta.mult.null<-array(dim=c(nhab,nhab))
pval.pair.gamma.null<-array(dim=c(nhab,nhab))
pval.pair.gamma.shan.null<-array(dim=c(nhab,nhab))
pval.pair.gamma.shan.mult.null<-array(dim=c(nhab,nhab))

# Calculate p-values for pairwise comparisons
for (i in 1:nhab) {
  for (j in 1:nhab) {
    pval.pair.alpha.null[i,j]<-
      length(pair.alpha.null[i,j,pair.alpha.null[i,j,]>
      obs.results$PairDiff$AlphaDif[i,j]])/krand
    pval.pair.beta.null[i,j]<-
      length(pair.beta.null[i,j,pair.beta.null[i,j,]>
      obs.results$PairDiff$BetaAddDif[i,j]])/krand
    pval.pair.beta.mult.null[i,j]<-
      length(pair.beta.mult.null[i,j,pair.beta.mult.null[i,j,]>
      obs.results$PairDiff$BetaMultDif[i,j]])/krand
    pval.pair.shan.alpha.null[i,j]<-
      length(pair.shan.alpha.null[i,j,pair.shan.alpha.null[i,j,]>
      obs.results$PairDiff$ShanAlphaDif[i,j]])/krand
    pval.pair.shan.beta.null[i,j]<-
      length(pair.shan.beta.null[i,j,pair.shan.beta.null[i,j,]>
      obs.results$PairDiff$ShanBetaDif[i,j]])/krand
    pval.pair.shan.alpha.mult.null[i,j]<-
      length(pair.shan.alpha.mult.null[i,j,pair.shan.alpha.mult.null[i,j,]>
      obs.results$PairDiff$ShanAlphaMultDif[i,j]])/krand
    pval.pair.shan.beta.mult.null[i,j]<-
      length(pair.shan.beta.mult.null[i,j,pair.shan.beta.mult.null[i,j,]>
      obs.results$PairDiff$ShanBetaMultDif[i,j]])/krand
    pval.pair.gamma.null[i,j]<-
      length(pair.gamma.null[i,j,pair.gamma.null[i,j,]>
      obs.results$PairDiff$GammaDif[i,j]])/krand
    pval.pair.gamma.shan.null[i,j]<-
      length(pair.gamma.shan.null[i,j,pair.gamma.shan.null[i,j,]>
      obs.results$PairDiff$ShanGammaDif[i,j]])/krand
    pval.pair.gamma.shan.mult.null[i,j]<-
      length(pair.gamma.shan.mult.null[i,j,pair.gamma.shan.mult.null[i,j,]>
      obs.results$PairDiff$ShanMultGammaDif[i,j]])/krand
  }
} 

#Summarize results in data frame, and combine into a list
pval.habitat<-data.frame(Abundance=pval.abund,Alpha=pval.alpha,
   BetaAdd=pval.beta,Gamma=pval.gamma,BetaMult=pval.beta.mult,
   Shannon=pval.alpha.shan,BetaShan=pval.beta.shan,
   GammaShan=pval.gamma.shan,AlphaShanMult=pval.alpha.shan.mult,
   BetaShanMult=pval.beta.shan.mult,row.names=hab.names)
quant.habitat<-data.frame(AbundLo=quant.abund[,2],AbundHi=quant.abund[,4],
   AlphaLo=quant.alpha[,2],AlphaHi=quant.alpha[,4],BetaLo=quant.beta[,2],
   BetaHi=quant.beta[,4],GammaLo=quant.gamma[,2],GammaHi=quant.gamma[,4],
   ShanLo=quant.alpha.shan[,2],ShanHi=quant.alpha.shan[,4],
   BetaShanLo=quant.beta.shan[,2],BetaShanHi=quant.beta.shan[,4],
   GammaShanLo=quant.gamma.shan[,2],GammaShanHi=quant.gamma.shan[,4],
   row.names=hab.names)
pval.fstats<-data.frame(ProbAlphaFstat=pval.alpha.fstat,
   ProbBetaFstatAdd=pval.beta.fstat,
   ProbBetaFstatMult=pval.beta.mult.fstat,
   ProbAlphaShanFstat=pval.alpha.shan.fstat,
   ProbBetaShanFstat=pval.beta.shan.fstat,
   ProbAlphaShanMult=pval.alpha.shan.mult.fstat,
   ProbBetaShanMult=pval.beta.shan.mult.fstat)
pval.pair<-list(ProbAlphaDif=pval.pair.alpha.null,
   ProbBetaAddDif=pval.pair.beta.null,
   ProbBetaMultDif=pval.pair.beta.mult.null,
   ProbShanAlphaDif=pval.pair.shan.alpha.null,
   ProbShanBetaDif=pval.pair.shan.beta.null,
   ProbShanAlphaMultDif=pval.pair.shan.alpha.mult.null,
   ProbShanBetaMultDiff=pval.pair.shan.beta.mult.null,
   ProbGammaDif=pval.pair.gamma.null,
   ProbGammaShanDiff=pval.pair.gamma.shan.null,
   ProbGammaShanMultDif=pval.pair.gamma.shan.mult.null)

null.results<-list(PvaluesHabitats=pval.habitat,Lo025.Hi975=quant.habitat,
   PvaluesFstats=pval.fstats,PvaluesPairedDiff=pval.pair)

obs.results
null.results