|
1 |
+ |
rt <- function (phy) { |
|
2 |
+ |
phy$tip.label <- phy$tip.label[sample(length(phy$tip.label))] |
|
3 |
+ |
return(phy) |
|
4 |
+ |
} |
|
5 |
+ |
|
|
6 |
+ |
#' Phylogenetic beta diversity standardized for species beta diversity |
|
7 |
+ |
#' |
|
8 |
+ |
#' This function computes the standard effect size of phylogenetic beta |
|
9 |
+ |
#' diversity by correcting for changes in species beta diversity. The |
|
10 |
+ |
#' novelty of this function is its ability to utilize sparse community |
|
11 |
+ |
#' matrix making it possible to efficiently randomize very large |
|
12 |
+ |
#' community matrices spanning thousands of taxa and sites. |
|
13 |
+ |
#' |
|
14 |
+ |
#' @param x a (sparse) community matrix, i.e., an object of class |
|
15 |
+ |
#' matrix or Matrix. |
|
16 |
+ |
#' @param phy a phylogenetic tree (object of class phylo). |
|
17 |
+ |
#' @param reps Number of replications. |
|
18 |
+ |
#' @param index.family the family of dissimilarity indices including |
|
19 |
+ |
#' \dQuote{simpson}, \dQuote{sorensen} and \dQuote{jaccard}. |
|
20 |
+ |
#' @param model The null model for separating patterns from processes and |
|
21 |
+ |
#' for contrasting against alternative hypotheses. Available null models |
|
22 |
+ |
#' include: |
|
23 |
+ |
#' \itemize{ |
|
24 |
+ |
#' \item \dQuote{tipshuffle}: shuffles phylogenetic tip labels multiple times. |
|
25 |
+ |
#' \item \dQuote{rowwise}: shuffles sites (i.e., varying richness) and |
|
26 |
+ |
#' keeping species occurrence frequency constant. |
|
27 |
+ |
#' \item \dQuote{colwise}: shuffles species occurrence frequency and |
|
28 |
+ |
#' keeping site richness constant.} |
|
29 |
+ |
#' |
|
30 |
+ |
#' @param \dots Further arguments passed to or from other methods. |
|
31 |
+ |
#' @rdname phylobeta_ses |
|
32 |
+ |
#' @importFrom stats sd |
|
33 |
+ |
#' @importFrom ape keep.tip |
|
34 |
+ |
#' @references |
|
35 |
+ |
#' Proches, S., Wilson, J.R.U. & Cowling, R.M. (2006) How much evolutionary |
|
36 |
+ |
#' history in a 10 x 10m plot? \emph{Proceedings of Royal Society B} |
|
37 |
+ |
#' \strong{273}: 1143-1148. |
|
38 |
+ |
#' @examples |
|
39 |
+ |
#' library(ape) |
|
40 |
+ |
#' library(Matrix) |
|
41 |
+ |
#' tree <- read.tree(text ="((t1:1,t2:1)N2:1,(t3:1,t4:1)N3:1)N1;") |
|
42 |
+ |
#' com <- sparseMatrix(c(1,3,4,1,4,5,1,2,3,4,5,6,3,4,6), |
|
43 |
+ |
#' c(1,1,1,2,2,2,3,3,3,3,3,3,4,4,4),x=1, |
|
44 |
+ |
#' dimnames = list(paste0("g", 1:6), tree$tip.label)) |
|
45 |
+ |
#' |
|
46 |
+ |
#' phylobeta_ses(com, tree, model="rowwise") |
|
47 |
+ |
#' |
|
48 |
+ |
#' @return A data frame of results for each community or grid cell |
|
49 |
+ |
#' \itemize{ |
|
50 |
+ |
#' \item phylobeta_obs: Observed phylobeta in community |
|
51 |
+ |
#' \item phylobeta_rand_mean: Mean phylobeta in null communities |
|
52 |
+ |
#' \item phylobeta_rand_sd: Standard deviation of phylobeta in null communities |
|
53 |
+ |
#' \item phylobeta_obs_z: Standardized effect size of phylobeta vs. |
|
54 |
+ |
#' null communities \eqn{= (phylobeta_obs - phylobeta_rand_mean) |
|
55 |
+ |
#' / phylobeta_rand_sd} |
|
56 |
+ |
#' \item reps: Number of replicates |
|
57 |
+ |
#' } |
|
58 |
+ |
#' @export |
|
59 |
+ |
phylobeta_ses <- function(x, phy, index.family="simpson", |
|
60 |
+ |
model = c("tipshuffle", "rowwise", "colwise"), |
|
61 |
+ |
reps = 1000, ...) { |
|
62 |
+ |
|
|
63 |
+ |
colnames(x) <- gsub(" ", "_", colnames(x)) |
|
64 |
+ |
p <- keep.tip(phy, intersect(phy$tip.label, colnames(x))) |
|
65 |
+ |
x <- x[, intersect(p$tip.label, colnames(x))] |
|
66 |
+ |
|
|
67 |
+ |
new_phylobeta <- function(mat, tree, index.family) { |
|
68 |
+ |
switch(index.family, |
|
69 |
+ |
simpson = phylobeta(mat, tree, index.family = "sorensen")[[1]], |
|
70 |
+ |
sorensen = phylobeta(mat, tree, index.family= "sorensen")[[3]], |
|
71 |
+ |
jaccard = phylobeta(mat, tree, index.family = "jaccard")[[3]]) |
|
72 |
+ |
} |
|
73 |
+ |
|
|
74 |
+ |
pbd_obs <- new_phylobeta(x, p, index.family) |
|
75 |
+ |
mean_obs_z <- pbd_obs |
|
76 |
+ |
|
|
77 |
+ |
|
|
78 |
+ |
dnam <- labels(pbd_obs) |
|
79 |
+ |
sort_nam <- function(x, dnam) { |
|
80 |
+ |
d <- as.matrix(x)[dnam,dnam] |
|
81 |
+ |
as.dist(d) |
|
82 |
+ |
} |
|
83 |
+ |
|
|
84 |
+ |
pbd.rand <- switch(model, |
|
85 |
+ |
tipshuffle = sapply(seq_len(reps), function(i) |
|
86 |
+ |
as.vector(new_phylobeta(x, rt(p), index.family))), |
|
87 |
+ |
rowwise = sapply(seq_len(reps), function(i) |
|
88 |
+ |
as.vector(sort_nam(new_phylobeta(x[sample(nrow(x)),], |
|
89 |
+ |
p, index.family), dnam))), |
|
90 |
+ |
colwise = sapply(seq_len(reps), function(i) |
|
91 |
+ |
as.vector(new_phylobeta(x[,sample(ncol(x))], |
|
92 |
+ |
p, index.family)))) |
|
93 |
+ |
|
|
94 |
+ |
pbd_rand_mean <- apply(X = pbd.rand, MARGIN = 1, FUN = mean, na.rm = TRUE) |
|
95 |
+ |
mean_obs_z[] <- pbd_rand_mean |
|
96 |
+ |
pbd_rand_mean <- mean_obs_z |
|
97 |
+ |
|
|
98 |
+ |
pbd_rand_sd <- apply(X = pbd.rand, MARGIN = 1, FUN = sd, na.rm = TRUE) |
|
99 |
+ |
pbd_rand_sd[pbd_rand_sd == 0] <- 1 |
|
100 |
+ |
mean_obs_z[] <- pbd_rand_sd |
|
101 |
+ |
pbd_rand_sd <- mean_obs_z |
|
102 |
+ |
|
|
103 |
+ |
pbd_obs_z <- (as.vector(pbd_obs) - pbd_rand_mean)/pbd_rand_sd |
|
104 |
+ |
mean_obs_z[] <- pbd_obs_z |
|
105 |
+ |
pbd_obs_z <- mean_obs_z |
|
106 |
+ |
|
|
107 |
+ |
return(list(phylobeta_obs = pbd_obs, |
|
108 |
+ |
phylobeta_rand_mean = pbd_rand_mean, |
|
109 |
+ |
phylobeta_rand_sd = pbd_rand_sd, |
|
110 |
+ |
phylobeta_obs_z = pbd_obs_z, reps = reps)) |
|
111 |
+ |
|
|
112 |
+ |
} |