Copy functions from phyloregion

to avoid loading betapart.
See futureverse/future#554
Closes #2

Author: joelnitta

DESCRIPTION

@@ -33,7 +33,9 @@ Imports:
     assertthat,
     dplyr,
     future.apply,
-    phyloregion,
+    Matrix,
+    methods,
+    phangorn,
     progressr,
     purrr,
     stats,
@@ -51,7 +53,8 @@ Suggests:
     stringr,
     magrittr,
     covr,
-    picante
+    picante,
+    phyloregion
 Config/testthat/edition: 3
 Depends: 
     R (>= 3.5.0)

R/calc_biodiv_random.R

@@ -115,12 +115,12 @@ calc_biodiv_random <- function(comm, phy, phy_alt,
 
   # Calculations ----
 
-  # Convert comm to sparse matrix format for phyloregions
-  comm_sparse <- phyloregion::dense2sparse(comm)
+  # Convert comm to sparse matrix format
+  comm_sparse <- dense2sparse(comm)
 
   # Generate random community
   random_comm <- cpr_rand_comm(comm, null_model = null_model, n_iterations = n_iterations, thin = thin, seed = seed)
-  random_comm_sparse <- phyloregion::dense2sparse(random_comm)
+  random_comm_sparse <- dense2sparse(random_comm)
 
   # Calculate statistics for random community
   # - set up null vectors first
@@ -132,20 +132,20 @@ calc_biodiv_random <- function(comm, phy, phy_alt,
   rpe <- NULL
 
   # - calculate selected metrics
-  if ("pd" %in% metrics) pd <- phyloregion::PD(random_comm_sparse, phy)
-  if ("pd_alt" %in% metrics) pd_alt <- phyloregion::PD(random_comm_sparse, phy_alt)
+  if ("pd" %in% metrics) pd <- PD(random_comm_sparse, phy)
+  if ("pd_alt" %in% metrics) pd_alt <- PD(random_comm_sparse, phy_alt)
   # pd_alt is inferred by rpd
   if ("rpd" %in% metrics) {
-    if (is.null(pd)) pd <- phyloregion::PD(random_comm_sparse, phy)
-    if (is.null(pd_alt)) pd_alt <- phyloregion::PD(random_comm_sparse, phy_alt)
+    if (is.null(pd)) pd <- PD(random_comm_sparse, phy)
+    if (is.null(pd_alt)) pd_alt <- PD(random_comm_sparse, phy_alt)
     rpd <- pd / pd_alt
   }
   # pe_alt is inferred by rpe
-  if ("pe" %in% metrics) pe <- phyloregion::phylo_endemism(random_comm_sparse, phy, weighted = TRUE)
-  if ("pe_alt" %in% metrics) pe_alt <- phyloregion::phylo_endemism(random_comm_sparse, phy_alt, weighted = TRUE)
+  if ("pe" %in% metrics) pe <- phylo_endemism(random_comm_sparse, phy, weighted = TRUE)
+  if ("pe_alt" %in% metrics) pe_alt <- phylo_endemism(random_comm_sparse, phy_alt, weighted = TRUE)
   if ("rpe" %in% metrics) {
-    if (is.null(pe)) pe <- phyloregion::phylo_endemism(random_comm_sparse, phy, weighted = TRUE)
-    if (is.null(pe_alt)) pe_alt <- phyloregion::phylo_endemism(random_comm_sparse, phy_alt, weighted = TRUE)
+    if (is.null(pe)) pe <- phylo_endemism(random_comm_sparse, phy, weighted = TRUE)
+    if (is.null(pe_alt)) pe_alt <- phylo_endemism(random_comm_sparse, phy_alt, weighted = TRUE)
     rpe <- pe / pe_alt
   }
 

R/cpr_rand_test.R

@@ -151,8 +151,8 @@ cpr_rand_test <- function(comm, phy, null_model,
     #' convert to integer before numeric comparisons
     # - null_model
     assertthat::assert_that(
-    	assertthat::is.string(null_model) | inherits(null_model, "commsim"),
-    	msg = "'null_model' must be a string (character vector of length 1) or an object of class 'commsim'"
+      assertthat::is.string(null_model) | inherits(null_model, "commsim"),
+      msg = "'null_model' must be a string (character vector of length 1) or an object of class 'commsim'"
     )
     assertthat::assert_that(assertthat::noNA(null_model))
     # - n_reps
@@ -229,7 +229,7 @@ cpr_rand_test <- function(comm, phy, null_model,
       comm <- comm_df
     }
     #' @srrstats {UL1.2} Check for default-looking rownames
-    # Default rownames not allowed because phyloregion::dense2sparse() will convert them to NULL
+    # Default rownames not allowed because dense2sparse() will convert them to NULL
     assertthat::assert_that(
       !identical(rownames(comm), as.character(seq(nrow(comm)))),
       msg = "'comm' cannot have default row names (consecutive integers from 1 to the number of rows)"
@@ -364,7 +364,7 @@ cpr_rand_test <- function(comm, phy, null_model,
     phy$edge.length <- phy$edge.length / sum(phy$edge.length)
 
     # Make sparse community df
-    comm_sparse <- phyloregion::dense2sparse(comm)
+    comm_sparse <- dense2sparse(comm)
 
     # Calculate biodiversity metrics ----
 
@@ -398,26 +398,26 @@ cpr_rand_test <- function(comm, phy, null_model,
 
     # - calculate selected metrics
     if ("pd" %in% metrics) {
-      pd_obs <- phyloregion::PD(comm_sparse, phy)
+      pd_obs <- PD(comm_sparse, phy)
       ses_pd <- get_ses(random_vals, pd_obs, "pd")
     }
 
     if ("rpd" %in% metrics) {
-      if (!exists("pd_obs")) pd_obs <- phyloregion::PD(comm_sparse, phy)
-      pd_alt_obs <- phyloregion::PD(comm_sparse, phy_alt)
+      if (!exists("pd_obs")) pd_obs <- PD(comm_sparse, phy)
+      pd_alt_obs <- PD(comm_sparse, phy_alt)
       ses_pd_alt <- get_ses(random_vals, pd_alt_obs, "pd_alt")
       rpd_obs <- pd_obs / pd_alt_obs
       ses_rpd <- get_ses(random_vals, rpd_obs, "rpd")
     }
 
     if ("pe" %in% metrics) {
-      pe_obs <- phyloregion::phylo_endemism(comm_sparse, phy, weighted = TRUE)
+      pe_obs <- phylo_endemism(comm_sparse, phy, weighted = TRUE)
       ses_pe <- get_ses(random_vals, pe_obs, "pe")
     }
 
     if ("rpe" %in% metrics) {
-      if (!exists("pe_obs")) pe_obs <- phyloregion::phylo_endemism(comm_sparse, phy, weighted = TRUE)
-      pe_alt_obs <- phyloregion::phylo_endemism(comm_sparse, phy_alt, weighted = TRUE)
+      if (!exists("pe_obs")) pe_obs <- phylo_endemism(comm_sparse, phy, weighted = TRUE)
+      pe_alt_obs <- phylo_endemism(comm_sparse, phy_alt, weighted = TRUE)
       ses_pe_alt <- get_ses(random_vals, pe_alt_obs, "pe_alt")
       rpe_obs <- pe_obs / pe_alt_obs
       ses_rpe <- get_ses(random_vals, rpe_obs, "rpe")

R/get_ses.R

@@ -28,8 +28,8 @@
 #'     1:100,
 #'     ~ calc_biodiv_random(comm, phy, phy_alt, "independentswap", 1000L, metrics = "pe")
 #'   )
-#' comm_sparse <- phyloregion::dense2sparse(comm)
-#' pe_obs <- phyloregion::phylo_endemism(comm_sparse, phy, weighted = TRUE)
+#' comm_sparse <- dense2sparse(comm)
+#' pe_obs <- phylo_endemism(comm_sparse, phy, weighted = TRUE)
 #' get_ses(random_vals, pe_obs, "pe")
 #' }
 #' @autoglobal

R/utils.R

@@ -247,3 +247,61 @@ match_phylo_comm <- function(phy, comm) {
   res$comm <- comm[, res$phy$tip.label]
   return(res)
 }
+
+# Functions copied from phyloregion v1.0.6 under AGPL-3 ----
+
+# Corresponds to phyloregion::dense2sparse()
+dense2sparse <- function(x) {
+  x <- as.matrix(x)
+  Matrix::Matrix(x, sparse = TRUE)
+}
+
+# Corresponds to phyloregion:::phylo_community()
+phylo_community <- function(x, phy) {
+  el <- numeric(max(phy$edge))
+  el[phy$edge[, 2]] <- phy$edge.length
+  x <- x[, phy$tip.label]
+  anc <- phangorn::Ancestors(phy, seq_along(phy$tip.label))
+  anc <- mapply(c, seq_along(phy$tip.label), anc, SIMPLIFY = FALSE)
+  M <- Matrix::sparseMatrix(as.integer(rep(
+    seq_along(anc),
+    lengths(anc)
+  )), as.integer(unlist(anc)), x = 1L)
+  commphylo <- x %*% M
+  commphylo@x[commphylo@x > 1e-08] <- 1
+  list(Matrix = commphylo, edge.length = el)
+}
+
+# Corresponds to phyloregion::PD()
+PD <- function(x, phy) {
+  if (!methods::is(x, "sparseMatrix")) {
+    stop("x needs to be a sparse matrix!")
+  }
+  if (length(setdiff(colnames(x), phy$tip.label)) > 0) {
+    stop("There are species labels in community matrix missing in the tree!")
+  }
+  if (length(setdiff(phy$tip.label, colnames(x))) > 0) {
+    phy <- ape::keep.tip(phy, intersect(phy$tip.label, colnames(x)))
+  }
+  x <- x[, intersect(phy$tip.label, colnames(x))]
+  z <- phylo_community(x, phy)
+  (z$Matrix %*% z$edge.length)[, 1]
+}
+
+# Corresponds to phyloregion::phylo_endemism()
+phylo_endemism <- function(x, phy, weighted = TRUE) {
+  if (length(setdiff(colnames(x), phy$tip.label)) > 0) {
+    stop("There are species labels in community matrix missing in the tree!")
+  }
+  if (length(setdiff(phy$tip.label, colnames(x))) > 0) {
+    phy <- ape::keep.tip(phy, intersect(phy$tip.label, colnames(x)))
+  }
+  comm_phylo <- phylo_community(x, phy)
+  weights <- comm_phylo$Matrix %*% Matrix::Diagonal(x = 1 / Matrix::colSums(comm_phylo$Matrix))
+  if (weighted == FALSE) {
+    weights[weights < 1] <- 0
+  }
+  pd <- (weights %*% comm_phylo$edge.length)[, 1]
+  pd <- pd[row.names(x)]
+  return(pd)
+}

README.Rmd

@@ -53,7 +53,7 @@ library(canaper)
 
 data(phylocom)
 
-# Example community matrix including 4 "clumped" communities, 
+# Example community matrix including 4 "clumped" communities,
 # one "even" community, and one "random" community
 phylocom$comm
 
@@ -71,7 +71,7 @@ rand_test_results <- cpr_rand_test(phylocom$comm, phylocom$phy, null_model = "sw
 `cpr_rand_test` produces **a lot** of columns (nine per metric), so let's just look at a subset of them:
 
 ```{r rand-test-res}
-rand_test_results[,1:9]
+rand_test_results[, 1:9]
 ```
 
 This is a summary of the columns:
@@ -133,9 +133,10 @@ You can find DOIs for older versions by viewing the "Releases" menu on the right
 
 ## Licenses
 
-- Code: [MIT](LICENSE.md)
+- Original code: [MIT](LICENSE.md)
+- Functions from [phyloregion](https://github.com/darunabas/phyloregion): [AGPL-3](https://github.com/joelnitta/canaper/blob/main/data-raw/LICENSE-agpl-3.txt)
 - Example datasets
-  - `acacia`, `biod_example`: [GNU General Public License v3.0](https://github.com/joelnitta/canaper/blob/main/data-raw/LICENSE-gpl.txt)
+  - `acacia`, `biod_example`: [GPL-3](https://github.com/joelnitta/canaper/blob/main/data-raw/LICENSE-gpl-3.txt)
   - `phylocom`: [BSD-3-Clause](https://github.com/joelnitta/canaper/blob/main/data-raw/LICENSE-bsd3.txt)
 
 ## References

README.md

@@ -49,7 +49,7 @@ library(canaper)
 
 data(phylocom)
 
-# Example community matrix including 4 "clumped" communities, 
+# Example community matrix including 4 "clumped" communities,
 # one "even" community, and one "random" community
 phylocom$comm
 #>         sp1 sp10 sp11 sp12 sp13 sp14 sp15 sp17 sp18 sp19 sp2 sp20 sp21 sp22
@@ -91,29 +91,29 @@ value to the alternative value (relative PD, relative PE).
 ``` r
 set.seed(071421)
 rand_test_results <- cpr_rand_test(phylocom$comm, phylocom$phy, null_model = "swap")
-#> [1] "Dropping tips from the tree because they are not present in the community data:"
-#> [1] "sp16" "sp23" "sp27" "sp28" "sp30" "sp31" "sp32"
+#> Warning in match_phylo_comm(phy = phy, comm = comm): Dropping tips from the tree because they are not present in the community data: 
+#>  sp16, sp23, sp27, sp28, sp30, sp31, sp32
 ```
 
 `cpr_rand_test` produces **a lot** of columns (nine per metric), so
 let’s just look at a subset of them:
 
 ``` r
-rand_test_results[,1:9]
+rand_test_results[, 1:9]
 #>            pd_obs pd_rand_mean pd_rand_sd  pd_obs_z pd_obs_c_upper
-#> clump1  0.3018868    0.4692453 0.03214267 -5.206739              0
-#> clump2a 0.3207547    0.4762264 0.03263836 -4.763465              0
-#> clump2b 0.3396226    0.4681132 0.03462444 -3.710978              0
-#> clump4  0.4150943    0.4667925 0.03180131 -1.625660              3
-#> even    0.5660377    0.4660377 0.03501739  2.855724            100
-#> random  0.5094340    0.4733962 0.03070539  1.173662             79
+#> clump1  0.3018868    0.4675472 0.03623666 -4.571624              0
+#> clump2a 0.3207547    0.4684906 0.03116570 -4.740335              0
+#> clump2b 0.3396226    0.4684906 0.03150994 -4.089754              0
+#> clump4  0.4150943    0.4664151 0.03307178 -1.551799              3
+#> even    0.5660377    0.4641509 0.03517108  2.896891            100
+#> random  0.5094340    0.4713208 0.03295196  1.156629             80
 #>         pd_obs_c_lower pd_obs_q pd_obs_p_upper pd_obs_p_lower
 #> clump1             100      100           0.00           1.00
 #> clump2a            100      100           0.00           1.00
 #> clump2b            100      100           0.00           1.00
-#> clump4              91      100           0.03           0.91
+#> clump4              90      100           0.03           0.90
 #> even                 0      100           1.00           0.00
-#> random               6      100           0.79           0.06
+#> random               7      100           0.80           0.07
 ```
 
 This is a summary of the columns:
@@ -145,7 +145,7 @@ canape_results[, "endem_type", drop = FALSE]
 #> clump2a not significant
 #> clump2b not significant
 #> clump4  not significant
-#> even              mixed
+#> even              super
 #> random            mixed
 ```
 
@@ -192,10 +192,13 @@ the right.
 
 ## Licenses
 
--   Code: [MIT](LICENSE.md)
+-   Original code: [MIT](LICENSE.md)
+-   Functions from
+    [phyloregion](https://github.com/darunabas/phyloregion):
+    [AGPL-3](https://github.com/joelnitta/canaper/blob/main/data-raw/LICENSE-agpl-3.txt)
 -   Example datasets
-    -   `acacia`, `biod_example`: [GNU General Public License
-        v3.0](https://github.com/joelnitta/canaper/blob/main/data-raw/LICENSE-gpl.txt)
+    -   `acacia`, `biod_example`:
+        [GPL-3](https://github.com/joelnitta/canaper/blob/main/data-raw/LICENSE-gpl-3.txt)
     -   `phylocom`:
         [BSD-3-Clause](https://github.com/joelnitta/canaper/blob/main/data-raw/LICENSE-bsd3.txt)
 

codemeta.json

@@ -205,6 +205,18 @@
         "url": "https://cran.r-project.org"
       },
       "sameAs": "https://CRAN.R-project.org/package=picante"
+    },
+    {
+      "@type": "SoftwareApplication",
+      "identifier": "phyloregion",
+      "name": "phyloregion",
+      "provider": {
+        "@id": "https://cran.r-project.org",
+        "@type": "Organization",
+        "name": "Comprehensive R Archive Network (CRAN)",
+        "url": "https://cran.r-project.org"
+      },
+      "sameAs": "https://CRAN.R-project.org/package=phyloregion"
     }
   ],
   "softwareRequirements": [
@@ -270,15 +282,32 @@
     },
     {
       "@type": "SoftwareApplication",
-      "identifier": "phyloregion",
-      "name": "phyloregion",
+      "identifier": "Matrix",
+      "name": "Matrix",
       "provider": {
         "@id": "https://cran.r-project.org",
         "@type": "Organization",
         "name": "Comprehensive R Archive Network (CRAN)",
         "url": "https://cran.r-project.org"
       },
-      "sameAs": "https://CRAN.R-project.org/package=phyloregion"
+      "sameAs": "https://CRAN.R-project.org/package=Matrix"
+    },
+    {
+      "@type": "SoftwareApplication",
+      "identifier": "methods",
+      "name": "methods"
+    },
+    {
+      "@type": "SoftwareApplication",
+      "identifier": "phangorn",
+      "name": "phangorn",
+      "provider": {
+        "@id": "https://cran.r-project.org",
+        "@type": "Organization",
+        "name": "Comprehensive R Archive Network (CRAN)",
+        "url": "https://cran.r-project.org"
+      },
+      "sameAs": "https://CRAN.R-project.org/package=phangorn"
     },
     {
       "@type": "SoftwareApplication",
@@ -341,7 +370,7 @@
     }
   ],
   "releaseNotes": "https://github.com/joelnitta/canaper/blob/master/NEWS.md",
-  "fileSize": "7891.09KB",
+  "fileSize": "7891.857KB",
   "contIntegration": ["https://github.com/joelnitta/canaper/actions", "https://codecov.io/gh/joelnitta/canaper?branch=main"],
   "developmentStatus": "https://www.repostatus.org/#wip",
   "keywords": [