R/phylodiversity.R

In picante: Integrating Phylogenies and Ecology

`comm.phylo.cor` <-
function(samp,phylo,metric=c("cij","checkerboard","jaccard","doij"),
		null.model=c("sample.taxa.labels","pool.taxa.labels",
					"frequency","richness","independentswap","trialswap"),
					runs=999, ...)
{
	metric <- match.arg(metric)
	null.model <- match.arg(null.model)
	results <- list("obs.corr"=NA,"obs.corr.p"=NA,"obs.rank"=NA,"runs"=runs,
			"obs.rand.p"=NA,"random.corrs"=vector(length=runs))
	phylo.dist <- as.dist(cophenetic(prune.sample(samp,phylo)))
	pool.phylo.dist <- as.dist(cophenetic(phylo))
	taxa.names <- rownames(as.matrix(phylo.dist))
	samp.dist <- as.dist(as.matrix(species.dist(samp,metric))[taxa.names,taxa.names])
	results$obs.corr <- cor(phylo.dist,samp.dist,use="pairwise")
	results$obs.corr.p <- cor.test(phylo.dist,samp.dist)$p.value
	if (null.model=="sample.taxa.labels") for (run in 1:runs)
	{
		phylo.dist <- as.dist(taxaShuffle(as.matrix(phylo.dist))[taxa.names,taxa.names])
		results$random.corrs[run] <- cor(phylo.dist,samp.dist,use="pairwise")
	}
	else if (null.model=="pool.taxa.labels") for (run in 1:runs)
	{
		phylo.dist <- as.dist(taxaShuffle(as.matrix(pool.phylo.dist))[taxa.names,taxa.names])
		results$random.corrs[run] <- cor(phylo.dist,samp.dist,use="pairwise")
	}
	else for (run in 1:runs)
	{
		samp.dist <- species.dist(randomizeMatrix(samp,null.model,...),metric)
		results$random.corrs[run] <- cor(phylo.dist,samp.dist,use="pairwise")
	}
	results$obs.rank <- rank(as.vector(c(results$obs.corr,results$random.corrs)))[1]
	results$obs.rand.p <- results$obs.rank/(runs+1)
	results
}


`comm.phylo.qr` <-
function(samp,phylo,metric=c("cij","checkerboard","jaccard","doij"),
		null.model=c("sample.taxa.labels","pool.taxa.labels",
					"frequency","richness","independentswap","trialswap"),
					quant=0.75, runs=999, show.plot=FALSE, ...)
{

    if (!requireNamespace("quantreg")) {
        stop("The 'quantreg' package is required to use this function.")
    }

	metric <- match.arg(metric)
	null.model <- match.arg(null.model)
	results <- list("obs.qr.intercept"=NA,"obs.qr.slope"=NA,"obs.qr.slope.p"=NA,"obs.rank"=NA,"runs"=runs,
			"random.qr.slopes"=vector(length=runs))
	phylo.dist <- as.dist(cophenetic(prune.sample(samp,phylo)))
	pool.phylo.dist <- as.dist(cophenetic(phylo))
	taxa.names <- rownames(as.matrix(phylo.dist))
	samp.dist <- as.dist(as.matrix(species.dist(samp,metric))[taxa.names,taxa.names])
    results$quantile <- quant
    qrres <- coef(quantreg::rq(samp.dist~phylo.dist,tau=quant, na.action=na.omit))
    names(qrres) <- NULL
    results$obs.qr.intercept <- qrres[1]
	results$obs.qr.slope <- qrres[2]

	if (show.plot) {
        plot(samp.dist~phylo.dist,xlab="Phylogenetic distance",ylab="Co-occurrence")
    }

	if (null.model=="sample.taxa.labels") for (run in 1:runs)
	{
		phylo.dist <- as.dist(taxaShuffle(as.matrix(phylo.dist))[taxa.names,taxa.names])
		results$random.qr.slopes[run] <- coef(quantreg::rq(samp.dist~phylo.dist,tau=quant,
		    na.action=na.omit))[2]
	}
	else if (null.model=="pool.taxa.labels") for (run in 1:runs)
	{
		phylo.dist <- as.dist(taxaShuffle(as.matrix(pool.phylo.dist))[taxa.names,taxa.names])
		results$random.qr.slopes[run] <- coef(quantreg::rq(samp.dist~phylo.dist,tau=quant,
		    na.action=na.omit))[2]
	}
	else for (run in 1:runs)
	{
		samp.dist <- species.dist(randomizeMatrix(samp,null.model,...),metric)
		results$random.qr.slopes[run] <- coef(quantreg::rq(samp.dist~phylo.dist,tau=quant,
		    na.action=na.omit))[2]
	}
	results$obs.rank <- rank(as.vector(c(results$obs.qr.slope,results$random.qr.slopes)))[1]
	results$obs.qr.slope.p <- results$obs.rank/(runs+1)

	if (show.plot) {
        abline(results$obs.qr.intercept,results$obs.qr.slope)
        legend("topleft",paste("q",as.character(quant),sep=""))
    }

	results
}


`taxaShuffle` <-
function(x) {
    #TODO replace with vegan's permuted.index?
    if (!is.matrix(x)) x <- as.matrix(x)
	rand.names <- sample(rownames(x))
	rownames(x) <- rand.names
	colnames(x) <- rand.names
	return(x)
}

`tipShuffle` <-
function(phy) {
    phy$tip.label <- phy$tip.label[sample(length(phy$tip.label))]
    return(phy)
}


mpd <- function(samp, dis, abundance.weighted=FALSE)
{
    N <- dim(samp)[1]
    mpd <- numeric(N)
    for (i in 1:N) {
        sppInSample <- names(samp[i, samp[i, ] > 0])
        if (length(sppInSample) > 1) {
            sample.dis <- dis[sppInSample, sppInSample]
            if (abundance.weighted) {
                sample.weights <- t(as.matrix(samp[i,sppInSample,drop=FALSE])) %*% as.matrix(samp[i,sppInSample,drop=FALSE])
                mpd[i] <- weighted.mean(sample.dis,sample.weights)

            }
            else {
                mpd[i] <- mean(sample.dis[lower.tri(sample.dis)])
            }
        }
        else{
            mpd[i] <- NA
        }
    }
    mpd
}


mntd <- function(samp, dis, abundance.weighted=FALSE)
{
	N <- dim(samp)[1]
	mntd <- numeric(N)
	for (i in 1:N) {
		sppInSample <- names(samp[i,samp[i,] > 0])
		if (length(sppInSample) > 1) {
            sample.dis <- dis[sppInSample,sppInSample]
            diag(sample.dis) <- NA
            if (abundance.weighted)
            {
                mntds <- apply(sample.dis,2,min,na.rm=TRUE)
                sample.weights <- samp[i,sppInSample]
                mntd[i] <- weighted.mean(mntds, sample.weights)
            }
            else
            {
		        mntd[i] <- mean(apply(sample.dis,2,min,na.rm=TRUE))
		    }
		}
		else {
		    mntd[i] <- NA
		}
	}
	mntd
}


`ses.mpd` <-
function (samp, dis, null.model = c("taxa.labels", "richness", "frequency", "sample.pool",
    "phylogeny.pool", "independentswap", "trialswap"),
    abundance.weighted = FALSE, runs = 999, iterations = 1000)
{
    dis <- as.matrix(dis)
    mpd.obs <- mpd(samp, dis, abundance.weighted = abundance.weighted)
    null.model <- match.arg(null.model)
    mpd.rand <- switch(null.model,
    	taxa.labels = t(replicate(runs, mpd(samp, taxaShuffle(dis), abundance.weighted=abundance.weighted))),
    	richness = t(replicate(runs, mpd(randomizeMatrix(samp,null.model="richness"), dis, abundance.weighted))),
    	frequency = t(replicate(runs, mpd(randomizeMatrix(samp,null.model="frequency"), dis, abundance.weighted))),
    	sample.pool = t(replicate(runs, mpd(randomizeMatrix(samp,null.model="richness"), dis, abundance.weighted))),
    	phylogeny.pool = t(replicate(runs, mpd(randomizeMatrix(samp,null.model="richness"),
    		taxaShuffle(dis), abundance.weighted))),
    	independentswap = t(replicate(runs, mpd(randomizeMatrix(samp,null.model="independentswap", iterations), dis, abundance.weighted))),
    	trialswap = t(replicate(runs, mpd(randomizeMatrix(samp,null.model="trialswap", iterations), dis, abundance.weighted)))
    )
    mpd.rand.mean <- apply(X = mpd.rand, MARGIN = 2, FUN = mean, na.rm=TRUE)
    mpd.rand.sd <- apply(X = mpd.rand, MARGIN = 2, FUN = sd, na.rm=TRUE)
    mpd.obs.z <- (mpd.obs - mpd.rand.mean)/mpd.rand.sd
    mpd.obs.rank <- apply(X = rbind(mpd.obs, mpd.rand), MARGIN = 2,
        FUN = rank)[1, ]
    mpd.obs.rank <- ifelse(is.na(mpd.rand.mean),NA,mpd.obs.rank)
    data.frame(ntaxa=specnumber(samp),mpd.obs, mpd.rand.mean, mpd.rand.sd, mpd.obs.rank,
        mpd.obs.z, mpd.obs.p=mpd.obs.rank/(runs+1),runs=runs, row.names = row.names(samp))
}


`ses.mntd` <-
function (samp, dis, null.model = c("taxa.labels", "richness", "frequency", "sample.pool",
    "phylogeny.pool", "independentswap", "trialswap"),
    abundance.weighted = FALSE, runs = 999, iterations = 1000)
{
    dis <- as.matrix(dis)
    mntd.obs <- mntd(samp, dis, abundance.weighted)
    null.model <- match.arg(null.model)
    mntd.rand <- switch(null.model,
    	taxa.labels = t(replicate(runs, mntd(samp, taxaShuffle(dis), abundance.weighted))),
    	richness = t(replicate(runs, mntd(randomizeMatrix(samp,null.model="richness"), dis, abundance.weighted))),
    	frequency = t(replicate(runs, mntd(randomizeMatrix(samp,null.model="frequency"), dis, abundance.weighted))),
    	sample.pool = t(replicate(runs, mntd(randomizeMatrix(samp,null.model="richness"), dis, abundance.weighted))),
    	phylogeny.pool = t(replicate(runs, mntd(randomizeMatrix(samp,null.model="richness"),
    		taxaShuffle(dis), abundance.weighted))),
    	independentswap = t(replicate(runs, mntd(randomizeMatrix(samp,null.model="independentswap", iterations), dis, abundance.weighted))),
    	trialswap = t(replicate(runs, mntd(randomizeMatrix(samp,null.model="trialswap", iterations), dis, abundance.weighted)))
    )
    mntd.rand.mean <- apply(X = mntd.rand, MARGIN = 2, FUN = mean, na.rm=TRUE)
    mntd.rand.sd <- apply(X = mntd.rand, MARGIN = 2, FUN = sd, na.rm=TRUE)
    mntd.obs.z <- (mntd.obs - mntd.rand.mean)/mntd.rand.sd
    mntd.obs.rank <- apply(X = rbind(mntd.obs, mntd.rand), MARGIN = 2,
        FUN = rank)[1, ]
    mntd.obs.rank <- ifelse(is.na(mntd.rand.mean),NA,mntd.obs.rank)
    data.frame(ntaxa=specnumber(samp),mntd.obs, mntd.rand.mean, mntd.rand.sd, mntd.obs.rank,
        mntd.obs.z, mntd.obs.p=mntd.obs.rank/(runs+1),runs=runs, row.names = row.names(samp))
}


psv<-function(samp,tree,compute.var=TRUE,scale.vcv=TRUE){
  # Make samp matrix a pa matrix
  samp[samp>0]<-1

  flag=0
  if (is.null(dim(samp))) #if the samp matrix only has one site
  {
    samp<-rbind(samp,samp)
    flag=2
  }

  if(is(tree)[1]=="phylo")
  {
    if(is.null(tree$edge.length)){tree<-compute.brlen(tree, 1)}  #If phylo has no given branch lengths
    tree<-prune.sample(samp,tree)
    # Make sure that the species line up
    samp<-samp[,tree$tip.label]
    # Make a correlation matrix of the species pool phylogeny
    Cmatrix<-vcv.phylo(tree,corr=scale.vcv)
  } else {
    Cmatrix<-tree
    species<-colnames(samp)
    preval<-colSums(samp)/sum(samp)
    species<-species[preval>0]
    Cmatrix<-Cmatrix[species,species]
    samp<-samp[,colnames(Cmatrix)]
  }

    # numbers of locations and species
    SR<-rowSums(samp)
    nlocations<-dim(samp)[1]
    nspecies<-dim(samp)[2]

    ##################################
    # calculate observed PSVs
    #
    PSVs<-NULL

    for(i in 1:nlocations)
    {
      index<-seq(1:nrow(Cmatrix))[samp[i,]==1]	#species present
      n<-length(index)			#number of species present
      if(n>1)
      {
        C<-Cmatrix[index,index]	#C for individual locations
        PSV<-(n*sum(diag(as.matrix(C)))-sum(C))/(n*(n-1))
      } else {
        PSV<-NA
      }
        PSVs<-c(PSVs,PSV)
    }
    PSVout<-cbind(PSVs,SR)

  if(flag==2)
  {
    PSVout<-PSVout[-2,]
    return(PSVout)
  } else {
    if(compute.var==FALSE)
    {
      return(data.frame(PSVout))
    } else {

      PSVvar<-NULL
      X<-Cmatrix-(sum(sum(Cmatrix-diag(nspecies))))/(nspecies*(nspecies-1));
      X<-X-diag(diag(X))

      SS1<-sum(X*X)/2

      SS2<-0;
      for(i in 1:(nspecies-1)){
        for(j in (i+1):nspecies){
          SS2<-SS2+X[i,j]*(sum(X[i,])-X[i,j])
        }
      }
      SS3<--SS1-SS2

      S1<-SS1*2/(nspecies*(nspecies-1))
      S2<-SS2*2/(nspecies*(nspecies-1)*(nspecies-2))

      if(nspecies==3){
        S3<-0
      }else{
      S3<-SS3*2/(nspecies*(nspecies-1)*(nspecies-2)*(nspecies-3))
      }

      for(n in 2:nspecies){
        approxi<-2/(n*(n-1))*(S1 + (n-2)*S2 + (n-2)*(n-3)*S3)
        PSVvar<-rbind(PSVvar, c(n, approxi))
      }

      vars<-rep(0,nlocations)
      PSVout<-cbind(PSVout,vars)

      for (g in 1:nlocations)
      {
        if(PSVout[g,2]>1)
        {
          PSVout[g,3]<-PSVvar[PSVout[g,2]-1,2]
        } else {
          PSVout[g,3]<-NA
        }
      }
      return(data.frame(PSVout))
    }
  }
}

psr <- function(samp,tree,compute.var=TRUE,scale.vcv=TRUE){
  PSVout<-psv(samp,tree,compute.var,scale.vcv=scale.vcv)
  if(is.null(dim(PSVout))==TRUE)
  {
    PSRout<-data.frame(cbind(PSVout[1]*PSVout[2],PSVout[2]))
    names(PSRout)<-c("PSR","SR")
    rownames(PSRout)<-""
    return(PSRout)
  } else {
    PSRout<-cbind(PSVout[,1]*PSVout[,2],PSVout[,2])
    if(compute.var!=TRUE) {
      colnames(PSRout)<-c("PSR","SR")
      rownames(PSRout)<-rownames(PSVout)
      return(data.frame(PSRout))
    } else {
      PSRout<-cbind(PSRout,PSVout[,3]*(PSVout[,2])^2)
      colnames(PSRout)<-c("PSR","SR","vars")
      rownames(PSRout)<-rownames(PSVout)
      return(data.frame(PSRout))
    }
  }
}

pse<-function(samp,tree,scale.vcv=TRUE){
  flag=0
  if (is.null(dim(samp))) #if the samp matrix only has one site
  {
    samp<-rbind(samp,samp)
    flag=2
  }

  samp<-as.matrix(samp)

  if(is(tree)[1]=="phylo")
  {
    if(is.null(tree$edge.length)){tree<-compute.brlen(tree, 1)}  #If phylo has no given branch lengths
    tree<-prune.sample(samp,tree)
    # Make sure that the species line up
    samp<-samp[,tree$tip.label, drop=FALSE]
    # Make a correlation matrix of the species pool phylogeny
    Cmatrix<-vcv.phylo(tree,corr=scale.vcv)
  } else {
    Cmatrix<-tree
    species<-colnames(samp)
    preval<-colSums(samp)/sum(samp)
    species<-species[preval>0]
    Cmatrix<-Cmatrix[species,species]
    samp<-samp[,colnames(Cmatrix),drop=FALSE]
  }

  # numbers of locations and species
  SR<-apply(samp>0,1,sum)
  nlocations<-dim(samp)[1]
  nspecies<-dim(samp)[2]

  #################################
  # calculate observed phylogenetic species evenness
  PSEs<-NULL
  for(i in 1:nlocations){
    index<-seq(1,ncol(Cmatrix))[samp[i,]>0]	  # species present
    n<-length(index)                          #location species richness
    if(n>1)
    {
    C<-Cmatrix[index,index]		                #C for individual locations
    N<-sum(samp[i,])                          #location total abundance
    M<-samp[i,samp[i,]>0]                  #species abundance column
    mbar<-mean(M)                             #mean species abundance
    PSE<-(N*t(diag(as.matrix(C)))%*%M-t(M)%*%as.matrix(C)%*%M)/(N^2-N*mbar)   #phylogenetic species evenness
    } else {PSE<-NA}
    PSEs<-c(PSEs,PSE)
  }
  PSEout=cbind(PSEs,SR)
  if (flag==2)
  {
    PSEout<-PSEout[-2,]
    return(PSEout)
  } else {
    return(data.frame(PSEout))
  }
}



psc<-function(samp,tree,scale.vcv=TRUE){
  # Make samp matrix a pa matrix
  samp[samp>0]<-1
  flag=0
  if (is.null(dim(samp))) #if the samp matrix only has one site
  {
    samp<-rbind(samp,samp)
    flag=2
  }

  if(is(tree)[1]=="phylo")
  {
    if(is.null(tree$edge.length)){tree<-compute.brlen(tree, 1)}  #If phylo has no given branch lengths
    tree<-prune.sample(samp,tree)
    # Make sure that the species line up
    samp<-samp[,tree$tip.label]
    # Make a correlation matrix of the species pool phylogeny
    Cmatrix<-vcv.phylo(tree,corr=scale.vcv)
  } else {
    Cmatrix<-tree
    species<-colnames(samp)
    preval<-colSums(samp)/sum(samp)
    species<-species[preval>0]
    Cmatrix<-Cmatrix[species,species]
    samp<-samp[,colnames(Cmatrix)]
  }

  # numbers of locations and species
  SR<-rowSums(samp)
  nlocations<-dim(samp)[1]
  nspecies<-dim(samp)[2]

  ##################################
  # calculate observed PSCs
  #
  PSCs<-NULL

  for(i in 1:nlocations)
  {
    index<-seq(1:nrow(Cmatrix))[samp[i,]==1]	#species present
    n<-length(index)			#number of species present
    if(n>1)
    {
      C<-Cmatrix[index,index]	#C for individual locations
      diag(C)<--1
      PSC<- 1-(sum(apply(C,1,max))/n)
    } else {PSC<-NA}
      PSCs<-c(PSCs,PSC)
  }
    PSCout<-cbind(PSCs,SR)
 if (flag==2)
  {
    PSCout<-PSCout[-2,]
    return(PSCout)
  } else {
    return(data.frame(PSCout))
  }
}

psv.spp<-function(samp,tree){
  # Make samp matrix a pa matrix
  samp[samp>0]<-1
  if (is.null(dim(samp))) #if the samp matrix only has one site
  {
    samp<-rbind(samp,samp)
  }
  if(is(tree)[1]=="phylo")
  {
    if(is.null(tree$edge.length)){tree<-compute.brlen(tree, 1)}  #If phylo has no given branch lengths
    tree<-prune.sample(samp,tree)
    # Make sure that the species line up
    samp<-samp[,tree$tip.label]
    # Make a correlation matrix of the species pool phylogeny
    Cmatrix<-vcv.phylo(tree,corr=TRUE)
  } else {
    Cmatrix<-tree
    species<-colnames(samp)
    preval<-colSums(samp)/sum(samp)
    species<-species[preval>0]
    Cmatrix<-Cmatrix[species,species]
    samp<-samp[,colnames(Cmatrix)]
  }
  # reduce given Cmatrix to the species observed in samp
  samp<-samp[rowSums(samp)>1,] # prune out locations with <2 species

  #cull the species that are not found in the samp set after all communities with 1 species are removed
  preval<-colSums(samp)/sum(samp)
  indexcov<-preval>0
  Cmatrix<-Cmatrix[indexcov,indexcov]
  samp<-samp[,indexcov]

  obs.PSV<-mean(psv(samp,Cmatrix,compute.var=FALSE)[,1],na.rm=TRUE)

  # numbers of locations and species
  nlocations<-dim(samp)[1]
  nspecies<-dim(samp)[2]

  spp.PSVs<-NULL
  for (j in 1:nspecies)
  {
    spp.samp<-samp[,-j]
    spp.Cmatrix<-Cmatrix[-j,-j]
    spp.PSV<-mean(psv(spp.samp,spp.Cmatrix,compute.var=FALSE)[,1],na.rm=TRUE)
    spp.PSVs<-c(spp.PSVs,spp.PSV)
  }
  spp.PSVout<-(spp.PSVs-obs.PSV)/sum(abs(spp.PSVs-obs.PSV))
  names(spp.PSVout)<-colnames(samp)
  return(spp.PSVout)
}

psd<-function(samp,tree,compute.var=TRUE,scale.vcv=TRUE){
  if (is.null(dim(samp))) #if the samp matrix only has one site
  {
    PSDout<-data.frame(c(psv(samp,tree,compute.var,scale.vcv)[1],psc(samp,tree,scale.vcv)[1],psr(samp,tree,compute.var,scale.vcv)[1],pse(samp,tree,scale.vcv)))
    names(PSDout)<-c("PSV","PSC","PSR","PSE","SR")
    return(PSDout)
  } else {
    if (compute.var==TRUE)
    {
      PSDout<-cbind(psv(samp,tree,compute.var,scale.vcv)[,c(1,3)],psc(samp,tree,scale.vcv)[,1],psr(samp,tree,compute.var,scale.vcv)[,c(1,3)],pse(samp,tree,scale.vcv))
      colnames(PSDout)<-c("PSV","var.PSV","PSC","PSR","var.PSR","PSE","SR")
    } else {
      PSDout<-cbind(psv(samp,tree,compute.var,scale.vcv)[,1],psc(samp,tree,scale.vcv)[,1],psr(samp,tree,compute.var,scale.vcv)[,1],pse(samp,tree,scale.vcv))
      colnames(PSDout)<-c("PSV","PSC","PSR","PSE","SR")
    }
    return(data.frame(PSDout))
  }
}

pd <- function (samp, tree, include.root = TRUE) {
		if (is.null(tree$edge.length)) {
      stop("Tree has no branch lengths, cannot compute pd")
    }
		if (include.root) {
			# Make sure tree is rooted if needed
			if (!is.rooted(tree)) {
				stop("Rooted tree required to calculate PD with include.root=TRUE argument")
				}
    	tree <- node.age(tree)
			}
    species <- colnames(samp)
    SR <- rowSums(ifelse(samp > 0, 1, 0))
    nlocations <- dim(samp)[1]
    nspecies <- dim(samp)[2]
    PDs <- rep(NA, nlocations)

    for (i in 1:nlocations) {
        present <- species[samp[i, ] > 0]
        treeabsent <- tree$tip.label[which(!(tree$tip.label %in% present))]
        # Check that sample has species; If no species are present, PD = 0
        if (length(present) == 0) {
            PDs[i] <- 0
        }
        # Check if there is only ONE species in sample
        else if (length(present) == 1) {
            # If tree is not rooted, parse error message
            if (!is.rooted(tree) || !include.root) {
                warning("Rooted tree and include.root=TRUE argument required to calculate PD of single-species communities. Single species community assigned PD value of NA.")
                PDs[i] <- NA
            }
            # Else the PD is the node age of that single species present
            else {
                PDs[i] <- tree$ages[which(tree$edge[, 2] ==
                  which(tree$tip.label == present))]
            }
        }
        # If there are no absent species (all tips are in sample) then PD
        # is the sum of all edges in the tree
        else if (length(treeabsent) == 0) {
            PDs[i] <- sum(tree$edge.length)
        }
        # Otherwise, we will need to remove absent species
        else {
            # Make a SubTree with only present species
            sub.tree <- drop.tip(tree, treeabsent)
            # If the tree is rooted, you need to account for different in
            # maximum branch length between subtree and original
            if (include.root) {
                # Make sure tree is rooted if needed
                if (!is.rooted(tree)) {
                  stop("Rooted tree required to calculate PD with include.root=TRUE argument")
                }
                # Calculate diff btw max depth original and max depth subtree
                sub.tree.depth <- max(node.age(sub.tree)$ages)
                orig.tree.depth <- max(tree$ages[which(tree$edge[,2] %in% which(tree$tip.label %in% present))])
                # PD is the sum of edges in subtree, plus any difference in the
                # max distance (of only present tips) between original and sub
                PDs[i] <- sum(sub.tree$edge.length) + (orig.tree.depth -
                  sub.tree.depth)
            }
            # If root is not included, just add all edges
            else {
                PDs[i] <- sum(sub.tree$edge.length)
            }
        }
    }
    PDout <- data.frame(PD = PDs, SR = SR)
    rownames(PDout) <- rownames(samp)
    return(PDout)
}

ses.pd <- function(samp, tree, null.model = c("taxa.labels", "richness", "frequency",
    "sample.pool", "phylogeny.pool", "independentswap", "trialswap"),
    runs = 999, iterations = 1000, include.root=TRUE)
{
  if(include.root == TRUE) {
    pd.obs <- as.vector(pd(samp, tree, include.root=TRUE)$PD)
    null.model <- match.arg(null.model)
    pd.rand <- switch(null.model,
                      taxa.labels = t(replicate(runs, as.vector(pd(samp, tipShuffle(tree), include.root=TRUE)$PD))),
                      richness = t(replicate(runs, as.vector(pd(randomizeMatrix(samp,null.model="richness"), tree, include.root=TRUE)$PD))),
                      frequency = t(replicate(runs, as.vector(pd(randomizeMatrix(samp,null.model="frequency"), tree, include.root=TRUE)$PD))),
                      sample.pool = t(replicate(runs, as.vector(pd(randomizeMatrix(samp,null.model="richness"), tree, include.root=TRUE)$PD))),
                      phylogeny.pool = t(replicate(runs, as.vector(pd(randomizeMatrix(samp,null.model="richness"),
                                                                      tipShuffle(tree), include.root=TRUE)$PD))),
                      independentswap = t(replicate(runs, as.vector(pd(randomizeMatrix(samp,null.model="independentswap", iterations), tree, include.root=TRUE)$PD))),
                      trialswap = t(replicate(runs, as.vector(pd(randomizeMatrix(samp,null.model="trialswap", iterations), tree, include.root=TRUE)$PD)))
    )
    pd.rand.mean <- apply(X = pd.rand, MARGIN = 2, FUN = mean, na.rm=TRUE)
    pd.rand.sd <- apply(X = pd.rand, MARGIN = 2, FUN = sd, na.rm=TRUE)
    pd.obs.z <- (pd.obs - pd.rand.mean)/pd.rand.sd
    pd.obs.rank <- apply(X = rbind(pd.obs, pd.rand), MARGIN = 2,
                         FUN = rank)[1, ]
    pd.obs.rank <- ifelse(is.na(pd.rand.mean),NA,pd.obs.rank)
    return(data.frame(ntaxa=specnumber(samp),pd.obs, pd.rand.mean, pd.rand.sd, pd.obs.rank,
               pd.obs.z, pd.obs.p=pd.obs.rank/(runs+1),runs=runs, row.names = row.names(samp)))
    
  }
  
  if(include.root == FALSE) {
    pd.obs <- as.vector(pd(samp, tree, include.root=FALSE)$PD)
    null.model <- match.arg(null.model)
    pd.rand <- switch(null.model,
                      taxa.labels = t(replicate(runs, as.vector(pd(samp, tipShuffle(tree), include.root=FALSE)$PD))),
                      richness = t(replicate(runs, as.vector(pd(randomizeMatrix(samp,null.model="richness"), tree, include.root=FALSE)$PD))),
                      frequency = t(replicate(runs, as.vector(pd(randomizeMatrix(samp,null.model="frequency"), tree, include.root=FALSE)$PD))),
                      sample.pool = t(replicate(runs, as.vector(pd(randomizeMatrix(samp,null.model="richness"), tree, include.root=FALSE)$PD))),
                      phylogeny.pool = t(replicate(runs, as.vector(pd(randomizeMatrix(samp,null.model="richness"),
                                                                      tipShuffle(tree), include.root=FALSE)$PD))),
                      independentswap = t(replicate(runs, as.vector(pd(randomizeMatrix(samp,null.model="independentswap", iterations), tree, include.root=FALSE)$PD))),
                      trialswap = t(replicate(runs, as.vector(pd(randomizeMatrix(samp,null.model="trialswap", iterations), tree, include.root=FALSE)$PD)))
    )
    pd.rand.mean <- apply(X = pd.rand, MARGIN = 2, FUN = mean, na.rm=TRUE)
    pd.rand.sd <- apply(X = pd.rand, MARGIN = 2, FUN = sd, na.rm=TRUE)
    pd.obs.z <- (pd.obs - pd.rand.mean)/pd.rand.sd
    pd.obs.rank <- apply(X = rbind(pd.obs, pd.rand), MARGIN = 2,
                         FUN = rank)[1, ]
    pd.obs.rank <- ifelse(is.na(pd.rand.mean),NA,pd.obs.rank)
    return(data.frame(ntaxa=specnumber(samp),pd.obs, pd.rand.mean, pd.rand.sd, pd.obs.rank,
                      pd.obs.z, pd.obs.p=pd.obs.rank/(runs+1),runs=runs, row.names = row.names(samp)))
    
  }
  
  
}