Skip to content

Commit

Permalink
Fix bug in ses.pd
Browse files Browse the repository at this point in the history
  • Loading branch information
@skembel
skembel committed Mar 21, 2019
1 parent 59d28e2 commit 8a69998
Show file tree
Hide file tree
Showing 2 changed files with 47 additions and 17 deletions.
60 changes: 45 additions & 15 deletions R/phylodiversity.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -614,28 +614,58 @@ pd <- function (samp, tree, include.root = TRUE) {
return(PDout)
}

`ses.pd` <-
function (samp, tree, null.model = c("taxa.labels", "richness", "frequency", "sample.pool",
"phylogeny.pool", "independentswap", "trialswap"), runs = 999, iterations = 1000, ...)
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)
{
pd.obs <- as.vector(pd(samp, tree, ...)$PD)
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), ...)$PD))),
richness = t(replicate(runs, as.vector(pd(randomizeMatrix(samp,null.model="richness"), tree, ...)$PD))),
frequency = t(replicate(runs, as.vector(pd(randomizeMatrix(samp,null.model="frequency"), tree, ...)$PD))),
sample.pool = t(replicate(runs, as.vector(pd(randomizeMatrix(samp,null.model="richness"), tree, ...)$PD))),
phylogeny.pool = t(replicate(runs, as.vector(pd(randomizeMatrix(samp,null.model="richness"),
tipShuffle(tree), ...)$PD))),
independentswap = t(replicate(runs, as.vector(pd(randomizeMatrix(samp,null.model="independentswap", iterations), tree, ...)$PD))),
trialswap = t(replicate(runs, as.vector(pd(randomizeMatrix(samp,null.model="trialswap", iterations), tree, ...)$PD)))
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, ]
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)
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))
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)))

}


}

4 changes: 2 additions & 2 deletions man/ses.pd.Rd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@
\usage{
ses.pd(samp, tree, null.model = c("taxa.labels", "richness", "frequency",
"sample.pool", "phylogeny.pool", "independentswap", "trialswap"),
runs = 999, iterations = 1000, ...)
runs = 999, iterations = 1000, include.root=TRUE)
}
\arguments{
Expand All @@ -17,7 +17,7 @@ ses.pd(samp, tree, null.model = c("taxa.labels", "richness", "frequency",
\item{ null.model }{ Null model to use (see Details section for description) }
\item{ runs }{ Number of randomizations }
\item{ iterations }{ Number of iterations to use for each randomization (for independent swap and trial null models) }
\item{ ... }{ Additional arguments to \code{\link{pd}} function }
\item{ include.root }{ Include distance to root node in calculation of PD (see documentation in \code{\link{pd}} function }
}
\value{
Expand Down

0 comments on commit 8a69998

Please sign in to comment.