From 080d82bbcb78a0caad790cb507aee55665313ed0 Mon Sep 17 00:00:00 2001 From: benjaminhlina Date: Fri, 27 Sep 2024 07:58:35 -0300 Subject: [PATCH 01/35] Moved isotope_names default and error checking to be in SIBER call not nicherover --- R/extract_mu.R | 37 ++++++++++++++++--------------------- 1 file changed, 16 insertions(+), 21 deletions(-) diff --git a/R/extract_mu.R b/R/extract_mu.R index 63c1555..234a862 100644 --- a/R/extract_mu.R +++ b/R/extract_mu.R @@ -98,28 +98,9 @@ extract_mu <- function(data, 'nicheROVER' or 'SIBER'.") } - # defaults of isotpoe a and b - if (is.null(isotope_names)) { - isotope_names <- c("d13c", "d15n") - } - - # Check if isotope_names is a character vector - if (!is.vector(isotope_names) || !is.character(isotope_names)) { - cli::cli_abort("The supplied argument for 'isotope_names' must be a vector of characters.") - } - - # Check if isotope_names has exactly 2 elements - if (length(isotope_names) != 2) { - cli::cli_abort("The 'isotope_names' vector must have exactly 2 elements, representing isotope_a and isotope_b.") - } - - # # Check if isotope_b is character - # if (!is.character(isotope_b)) { - # cli::cli_abort("The supplied argument for 'isotope_b' must be a character.") - # } - # sett data formatt + # sett data format if (is.null(data_format)) { data_format <- "long" } @@ -191,6 +172,20 @@ extract_mu <- function(data, } + # defaults of isotpoe a and b + if (is.null(isotope_names)) { + isotope_names <- c("d13c", "d15n") + } + + # Check if isotope_names is a character vector + if (!is.vector(isotope_names) || !is.character(isotope_names)) { + cli::cli_abort("The supplied argument for 'isotope_names' must be a vector of characters.") + } + + # Check if isotope_names has exactly 2 elements + if (length(isotope_names) != 2) { + cli::cli_abort("'isotope_names' must have exactly 2 elements, representing isotope_a and isotope_b") + } id_isotope <- isotope_names @@ -199,7 +194,7 @@ extract_mu <- function(data, df <- .x[, 5:6] |> t() |> as.numeric() |> - matrix(ncol = 1, byrow = T) |> + matrix(ncol = 1, byrow = TRUE) |> as.data.frame() |> tibble::as_tibble() From ab64428c3e392db96416713b179862ca7fc621b3 Mon Sep 17 00:00:00 2001 From: benjaminhlina Date: Fri, 27 Sep 2024 07:59:29 -0300 Subject: [PATCH 02/35] Moved tests and updated them to check extract_mu when pkg is set to SIBER specificallly checking if isotope_names works --- tests/testthat/test-extract_mu.R | 55 ++++++++++++++++++-------------- 1 file changed, 31 insertions(+), 24 deletions(-) diff --git a/tests/testthat/test-extract_mu.R b/tests/testthat/test-extract_mu.R index 42d29ff..6058a88 100644 --- a/tests/testthat/test-extract_mu.R +++ b/tests/testthat/test-extract_mu.R @@ -91,30 +91,6 @@ test_that("Check if column names extracted are correct", { }) -test_that("test that isotope a and b error when not characters ", { - expect_error(extract_mu( - data = niw_fish_post, - data_format = "wide", - isotope_names = 1 - ), "The supplied argument for 'isotope_names' must be a vector of characters.") - -}) -test_that("test that isotope a and b error when not characters ", { - expect_error(extract_mu( - data = niw_fish_post, - data_format = "wide", - isotope_names = c(1, 2) - ), "The supplied argument for 'isotope_names' must be a vector of characters.") - -}) -test_that("test that isotope a and b error when not characters ", { - expect_error(extract_mu( - data = niw_fish_post, - data_format = "wide", - isotope_names = c("d13c") - ), "The 'isotope_names' vector must have exactly 2 elements, representing isotope_a and isotope_b.") - -}) # test_that("test that isotope a and b error when not characters ", { # expect_error(extract_mu( @@ -306,3 +282,34 @@ test_that("extract_mu throws an error if community_df is missing 'community' or community_df = invalid_community_df_2), regexp = "The data frame does not contain a column named 'community' and 'group'.") }) + +test_that("test that isotope a and b error when not characters ", { + expect_error(extract_mu( + data = post_sam_siber, + pkg = "SIBER", + community_df = cg_names, + data_format = "wide", + isotope_names = 1 + ), "The supplied argument for 'isotope_names' must be a vector of characters.") + +}) +test_that("test that isotope a and b error when not characters ", { + expect_error(extract_mu( + data = post_sam_siber, + pkg = "SIBER", + community_df = cg_names, + data_format = "wide", + isotope_names = c(1, 2) + ), "The supplied argument for 'isotope_names' must be a vector of characters.") + +}) +test_that("test that isotope a and b error when not characters ", { + expect_error(extract_mu( + data = post_sam_siber, + pkg = "SIBER", + community_df = cg_names, + data_format = "wide", + isotope_names = c("d13c") + ), "'isotope_names' must have exactly 2 elements, representing isotope_a and isotope_b") + +}) From 856a2c069ae63c5d4b605425e1885c9a9b49aef4 Mon Sep 17 00:00:00 2001 From: benjaminhlina Date: Fri, 27 Sep 2024 07:59:58 -0300 Subject: [PATCH 03/35] Updated the manual regarding isotope_names --- R/extract_mu.R | 4 +++- man/extract_mu.Rd | 6 ++++-- 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/R/extract_mu.R b/R/extract_mu.R index 234a862..9703583 100644 --- a/R/extract_mu.R +++ b/R/extract_mu.R @@ -12,7 +12,8 @@ #' you're using. Defaults to `"nicheROVER"`. #' Alternatively the user can supply the argument with `"SIBER"`. #' @param isotope_names is a vector of `character` string used change the column name -#' of isotopes used in the analysis. Defaults to `c("d13c", "d15n")`. +#' of isotopes used in the analysis. Defaults to `c("d13c", "d15n")`. To be used when +#' `pkg` is set to `"SIBER"`. #' @param data_format a `character` string that decides whether the returned object is #' in long or wide format. Default is `"long"`, with the alternative supplied #' being `"wide"`. @@ -27,6 +28,7 @@ #' from [{SIBER}](https://CRAN.R-project.org/package=SIBER). #' The third and fourth columns contains the actual names of the communities #' and groups the user is working with (e.g., `"region"`, `"common_name"`). +#' To be used when `pkg` is set to `"SIBER"`. #' #' @return Returns a `tibble` of extracted estimates of \eqn{\mu} created by the #' function `niw.post()` or `siberMVN()` in the packages diff --git a/man/extract_mu.Rd b/man/extract_mu.Rd index 7c1b669..5222b3a 100644 --- a/man/extract_mu.Rd +++ b/man/extract_mu.Rd @@ -23,7 +23,8 @@ you're using. Defaults to \code{"nicheROVER"}. Alternatively the user can supply the argument with \code{"SIBER"}.} \item{isotope_names}{is a vector of \code{character} string used change the column name -of isotopes used in the analysis. Defaults to \code{c("d13c", "d15n")}.} +of isotopes used in the analysis. Defaults to \code{c("d13c", "d15n")}. To be used when +\code{pkg} is set to \code{"SIBER"}.} \item{data_format}{a \code{character} string that decides whether the returned object is in long or wide format. Default is \code{"long"}, with the alternative supplied @@ -39,7 +40,8 @@ The second column will be the names of the groups that are needed to supply required by the function, \code{createSiberObject()} from \href{https://CRAN.R-project.org/package=SIBER}{{SIBER}}. The third and fourth columns contains the actual names of the communities -and groups the user is working with (e.g., \code{"region"}, \code{"common_name"}).} +and groups the user is working with (e.g., \code{"region"}, \code{"common_name"}). +To be used when \code{pkg} is set to \code{"SIBER"}.} } \value{ Returns a \code{tibble} of extracted estimates of \eqn{\mu} created by the From b45f7064947dc8d820eb6b262144692d619d68a4 Mon Sep 17 00:00:00 2001 From: benjaminhlina Date: Mon, 30 Sep 2024 12:23:00 -0400 Subject: [PATCH 04/35] Trying to make niche_ellipse work with 3 --- R/niche_ellipse.R | 39 +++++++++++++++++++-------------------- 1 file changed, 19 insertions(+), 20 deletions(-) diff --git a/R/niche_ellipse.R b/R/niche_ellipse.R index d6ac53b..dc2fb4e 100644 --- a/R/niche_ellipse.R +++ b/R/niche_ellipse.R @@ -14,10 +14,10 @@ #' @param p_ell is the confidence interval of each ellipse estimate. #' Default is 0.95 (i.e., 95% confidence interval). #' This value is bound by 0 and 1 and has to be a `numeric`. -#' @param isotope_a character string that is the column name of the first -#' isotope used in `dat_sigma`. Defaults to `"d13c"`. -#' @param isotope_b character string that is the column name of the second -#' isotope used in `dat_sigma`. Defaults to `"d15n"`. +#' @param isotope_n a `numeric` either `2` or `3` that is the number of isotopes +#' used in the anlsysis. Will default to `2`. +#' @param isotope_names is a vector of `character` string used change the column name +#' of isotopes used in the analysis. Defaults to `c("d13c", "d15n")`. #' @param random logical value indicating whether or not to randomly sample #' posterior distributions for \eqn{\mu} and \eqn{\Sigma} to create a sub-sample #' of ellipse. Default is `TRUE`. @@ -55,8 +55,8 @@ niche_ellipse <- function( dat_mu, dat_sigma, - isotope_a = NULL, - isotope_b = NULL, + isotope_n = NULL, + isotope_names = NULL, p_ell = NULL, random = NULL, set_seed = NULL, @@ -75,26 +75,22 @@ niche_ellipse <- function( if (!inherits(dat_sigma, what = c("tbl_df", "tbl", "data.frame"))) { cli::cli_abort("Input 'dat_sigma' must be class data.frame.") } - - - if (is.null(isotope_a)) { - isotope_a <- "d13c" - } - if (is.null(isotope_b)) { - isotope_b <- "d15n" + # create name vector that will be used to id isotopes. + if (is.null(isotope_n)) { + isotope_n <- 2 } - # Validate if isotope_a is a character - if (!is.character(isotope_a)) { - cli::cli_abort("Argument 'isotope_a' must be a character.") + if (!is.numeric(isotope_n) || !(isotope_n %in% c(2, 3))) { + cli::cli_abort("Argument 'isotope_n' must be a numeric value and either 2 or 3.") } - # Validate if isotope_b is a character - if (!is.character(isotope_b)) { - cli::cli_abort("Argument 'isotope_b' must be a character.") + if (isotope_n == 2) { + if (is.null(isotope_names)) { + isotope_names <- c("d13c", "d15n") + } } - if(is.null(p_ell)) { + if (is.null(p_ell)) { p_ell <- 0.95 } # Additional parameter validation, if needed @@ -136,6 +132,9 @@ niche_ellipse <- function( set.seed(set_seed) sample_numbers <- sample(dat_mu$sample_number, n) + + isotope_a <- isotope_name[1] + isotope_b <- isotope_name[2] # prepare mu for ellipse mu <- dat_mu |> dplyr::select(sample_name, sample_number, mu_est) |> From b0b813fb91476f697c7797146478c3d2bd706a75 Mon Sep 17 00:00:00 2001 From: benjaminhlina Date: Mon, 30 Sep 2024 12:23:08 -0400 Subject: [PATCH 05/35] Update manual --- man/niche_ellipse.Rd | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/man/niche_ellipse.Rd b/man/niche_ellipse.Rd index 5ebf041..3426a67 100644 --- a/man/niche_ellipse.Rd +++ b/man/niche_ellipse.Rd @@ -7,8 +7,8 @@ niche_ellipse( dat_mu, dat_sigma, - isotope_a = NULL, - isotope_b = NULL, + isotope_n = NULL, + isotope_names = NULL, p_ell = NULL, random = NULL, set_seed = NULL, @@ -27,11 +27,11 @@ This \code{data.frame} needs be in wide format, that is \eqn{\Sigma} (covariance matrices stacked on top of each other. See example of how to convert to wide format. This can be produced using \code{extract_sigma()}.} -\item{isotope_a}{character string that is the column name of the first -isotope used in \code{dat_sigma}. Defaults to \code{"d13c"}.} +\item{isotope_n}{a \code{numeric} either \code{2} or \code{3} that is the number of isotopes +used in the anlsysis. Will default to \code{2}.} -\item{isotope_b}{character string that is the column name of the second -isotope used in \code{dat_sigma}. Defaults to \code{"d15n"}.} +\item{isotope_names}{is a vector of \code{character} string used change the column name +of isotopes used in the analysis. Defaults to \code{c("d13c", "d15n")}.} \item{p_ell}{is the confidence interval of each ellipse estimate. Default is 0.95 (i.e., 95\% confidence interval). From 27344d83e89f39af2e08accb9957accfb84531f0 Mon Sep 17 00:00:00 2001 From: benjaminhlina Date: Mon, 7 Oct 2024 11:12:41 -0400 Subject: [PATCH 06/35] Updated version number --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index ba7eae9..1701622 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: nichetools Type: Package Title: Complementary Package to 'nicheROVER' and 'SIBER' -Version: 0.3.1 +Version: 0.3.3 Authors@R: c(person(given = "Benjamin L.", family = "Hlina", From e5ed6d692abc46d0b44d91e89b9ae9b7c95f5a8f Mon Sep 17 00:00:00 2001 From: benjaminhlina Date: Mon, 7 Oct 2024 11:12:56 -0400 Subject: [PATCH 07/35] Updated version number and version name --- R/zzz.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/zzz.R b/R/zzz.R index bacd416..a9f960f 100644 --- a/R/zzz.R +++ b/R/zzz.R @@ -1,5 +1,5 @@ .onAttach <- function(libname, pkgname) { - packageStartupMessage("version 0.3.1 ('visiting-friends').\nHave you loaded {nicheROVER} or {SIBER}? If not, please do so.") + packageStartupMessage("version 0.3.3 ('let-the-snow-fly').\nHave you loaded {nicheROVER} or {SIBER}? If not, please do so.") } utils::globalVariables(c("sample_name", "sample_number", "metric", "x", "y", "V1", "V2", From 3ceeb4fed1d5beb8498bbabc8f83a3ff1426e39e Mon Sep 17 00:00:00 2001 From: benjaminhlina Date: Mon, 7 Oct 2024 11:24:28 -0400 Subject: [PATCH 08/35] Fixed format spacing issue for extract_layman function --- R/extract_layman.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/extract_layman.R b/R/extract_layman.R index 19b80d0..18f737b 100644 --- a/R/extract_layman.R +++ b/R/extract_layman.R @@ -183,7 +183,7 @@ extract_layman <- function(data, element_y <- "N" } - if(type %in% "bay") { + if (type %in% "bay") { df_layman <- data |> purrr::map(~ as_tibble(.x)) |> From 020844e74f5e32dc1d8308564b75a5775ec68d25 Mon Sep 17 00:00:00 2001 From: benjaminhlina Date: Mon, 7 Oct 2024 11:25:24 -0400 Subject: [PATCH 09/35] Minor format changes --- R/extract_group_metrics.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/extract_group_metrics.R b/R/extract_group_metrics.R index 8d2d544..0ab83ef 100644 --- a/R/extract_group_metrics.R +++ b/R/extract_group_metrics.R @@ -90,7 +90,7 @@ extract_group_metrics <- function(data = NULL, # set data formatt - if(is.null(data_format)) { + if (is.null(data_format)) { data_format <- "long" } From faa4477a0d35dca07e8b53d45d5ae6440f9ef192 Mon Sep 17 00:00:00 2001 From: benjaminhlina Date: Mon, 7 Oct 2024 11:25:41 -0400 Subject: [PATCH 10/35] replaced magrittr pipes with base pipe --- R/extract_mu.R | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/R/extract_mu.R b/R/extract_mu.R index 9703583..b0ee74d 100644 --- a/R/extract_mu.R +++ b/R/extract_mu.R @@ -216,12 +216,12 @@ extract_mu <- function(data, dplyr::rename( mu_est = V1 ) |> - dplyr::select(metric, sample_name,sample_number, isotope, mu_est) %>% + dplyr::select(metric, sample_name,sample_number, isotope, mu_est) |> separate_wider_delim(sample_name, cols_remove = FALSE, delim = ".", names = c("community", - "group")) %>% + "group")) |> left_join(community_df, by = c("community", - "group")) %>% + "group")) |> dplyr::select(metric:sample_number, community_name, group_name, isotope, mu_est) From 2c03f8bb04c245f90257a98a5ba66d0c6d780e0b Mon Sep 17 00:00:00 2001 From: benjaminhlina Date: Mon, 7 Oct 2024 11:25:59 -0400 Subject: [PATCH 11/35] Minor formating changes to the code --- R/extract_overlap.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/extract_overlap.R b/R/extract_overlap.R index 4e3d952..f2cd002 100644 --- a/R/extract_overlap.R +++ b/R/extract_overlap.R @@ -32,11 +32,11 @@ extract_overlap <- function(data, } # set name_a null - if (is.null(name_a)){ + if (is.null(name_a)) { name_a <- "sample_name_a" } # set name_b null - if (is.null(name_b)){ + if (is.null(name_b)) { name_b <- "sample_name_b" } From 77efd631257f4afc5f8a6eb4aa502690b1858af2 Mon Sep 17 00:00:00 2001 From: benjaminhlina Date: Mon, 7 Oct 2024 11:26:22 -0400 Subject: [PATCH 12/35] replaced separate with separte_wider_delim --- R/extract_overlap.R | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/R/extract_overlap.R b/R/extract_overlap.R index f2cd002..a623433 100644 --- a/R/extract_overlap.R +++ b/R/extract_overlap.R @@ -60,12 +60,12 @@ extract_overlap <- function(data, tidyr::pivot_longer(cols = -c(id, species_a), names_to = "species_b", values_to = "niche_overlap") |> - tidyr::separate(species_b, into = c("species_b", "sample_number"), - sep = "\\.") |> + tidyr::separate_wider_delim(species_b, names = c("species_b", "sample_number"), + delim = ".") |> dplyr::rename( {{name_a}} := species_a, {{name_b}} := species_b - ) %>% + ) |> dplyr::mutate( niche_overlap_perc = niche_overlap * 100 ) From e07381ea5970914a29bed001d29af335eb8f7b05 Mon Sep 17 00:00:00 2001 From: benjaminhlina Date: Mon, 7 Oct 2024 11:26:52 -0400 Subject: [PATCH 13/35] replaced magrittr pipes with base pipes --- R/extract_sigma.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/extract_sigma.R b/R/extract_sigma.R index 09b399d..624a24e 100644 --- a/R/extract_sigma.R +++ b/R/extract_sigma.R @@ -109,7 +109,7 @@ extract_sigma <- function(data, dplyr::bind_rows(.id = "isotope_num") |> dplyr::mutate( id = 1 - ) %>% + ) |> tidyr::pivot_longer(-id) isotope_length <- unique(isotope_number$value) @@ -145,7 +145,7 @@ extract_sigma <- function(data, dplyr::bind_rows(.id = "isotope_num") |> dplyr::mutate( id = 1 - ) %>% + ) |> tidyr::pivot_longer(-id) isotope_length <- unique(isotope_number$value) From b9860818451aad997d0831c4b3252a1cc4208691 Mon Sep 17 00:00:00 2001 From: benjaminhlina Date: Mon, 7 Oct 2024 11:27:19 -0400 Subject: [PATCH 14/35] replaced separate with separate_wider_delim --- R/extract_sigma.R | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/R/extract_sigma.R b/R/extract_sigma.R index 624a24e..92c9894 100644 --- a/R/extract_sigma.R +++ b/R/extract_sigma.R @@ -171,8 +171,9 @@ extract_sigma <- function(data, names_to = "isotope", values_to = "post_sample" ) |> - tidyr::separate(isotope, into = c("isotope", "sample_number"), - sep = "\\.") + tidyr::separate_wider_delim(isotope, names = c("isotope", "sample_number"), + + delim = ".") From 0e9e393f3fa54253842b7b19d0d45e9ffb67c005 Mon Sep 17 00:00:00 2001 From: benjaminhlina Date: Mon, 7 Oct 2024 11:27:27 -0400 Subject: [PATCH 15/35] Replaced T with TRUE --- R/extract_sigma.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/extract_sigma.R b/R/extract_sigma.R index 92c9894..786a539 100644 --- a/R/extract_sigma.R +++ b/R/extract_sigma.R @@ -221,7 +221,7 @@ extract_sigma <- function(data, df <- .x[, 1:4] |> t() |> as.numeric() |> - matrix(ncol = 2, byrow = T) |> + matrix(ncol = 2, byrow = TRUE) |> as.data.frame() |> tibble::as_tibble() From 64c0fe3a4e571be32266282b91decdf904af5fa6 Mon Sep 17 00:00:00 2001 From: benjaminhlina Date: Mon, 7 Oct 2024 11:28:41 -0400 Subject: [PATCH 16/35] Update version number --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 5f41d5c..0c14f13 100644 --- a/README.md +++ b/README.md @@ -73,4 +73,4 @@ To cite this package please cite the following publications - Layman, C.A., Arrington, D.A., Montaña, C.G., and Post, D.M. 2007. Can stable isotope ratios provide for community-wide measures of trophic structure? Ecology 88(1): 42–48. [link]( https://doi.org/10.1890/0012-9658(2007)88[42:CSIRPF]2.0.CO;2) -- Hlina B.L. 2024. nichetools: Complementary package to nicheROVER and SIBER. R package version 0.3.1. https://benjaminhlina.github.io/nichetools/ +- Hlina B.L. 2024. nichetools: Complementary package to nicheROVER and SIBER. R package version 0.3.3. https://benjaminhlina.github.io/nichetools/ From 98424b20eae7ea0709ee3e39ee7f3ec0daf166be Mon Sep 17 00:00:00 2001 From: benjaminhlina Date: Mon, 7 Oct 2024 11:28:47 -0400 Subject: [PATCH 17/35] update version number --- inst/CITATION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/inst/CITATION b/inst/CITATION index 365f361..f877870 100644 --- a/inst/CITATION +++ b/inst/CITATION @@ -44,6 +44,6 @@ bibentry( title = "nichetools: Complementary Package to 'nicheROVER' and 'SIBER'", author = c("Benjamin L. Hlina"), year = "2024", - note = "R package version 0.3.1", + note = "R package version 0.3.3", url = "https://benjaminhlina.github.io/nichetools/" ) From 23638957535346a410886c7addb6d3b418830062 Mon Sep 17 00:00:00 2001 From: benjaminhlina Date: Mon, 7 Oct 2024 11:44:04 -0400 Subject: [PATCH 18/35] Minor formatting changes --- R/extract_layman.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/extract_layman.R b/R/extract_layman.R index 18f737b..dc61339 100644 --- a/R/extract_layman.R +++ b/R/extract_layman.R @@ -224,7 +224,7 @@ extract_layman <- function(data, return(df_layman) } - if (data_format %in% "wide"){ + if (data_format %in% "wide") { return(df_layman) } } From f4fc53be419b6dfedeeb3092baa71c427d19f818 Mon Sep 17 00:00:00 2001 From: benjaminhlina Date: Tue, 12 Nov 2024 14:29:11 -0500 Subject: [PATCH 19/35] Change from isotope_a and b to isotope_names --- R/niche_ellipse.R | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/R/niche_ellipse.R b/R/niche_ellipse.R index dc2fb4e..ff0d569 100644 --- a/R/niche_ellipse.R +++ b/R/niche_ellipse.R @@ -66,7 +66,6 @@ niche_ellipse <- function( # options(error = recover) start_time <- Sys.time() - # Check if dat_mu is a data.frame if (!inherits(x = dat_mu, what = c("tbl_df", "tbl", "data.frame"))) { cli::cli_abort("Input 'dat_mu' must be class data.frame.") @@ -87,6 +86,14 @@ niche_ellipse <- function( if (is.null(isotope_names)) { isotope_names <- c("d13c", "d15n") } + if (!is.character(isotope_names)) { + cli::cli_abort("Argument 'isotope_names' must be a character.") + } + } + if (isotope_n == 3) { + if (is.null(isotope_names)) { + isotope_names <- c("d13c", "d15n", "d34s") + } } From b3b09ae89a85529d2c6d79b26a865321cb859f7f Mon Sep 17 00:00:00 2001 From: benjaminhlina Date: Tue, 12 Nov 2024 14:29:58 -0500 Subject: [PATCH 20/35] Added new method to create ellipse for 3 isotopes --- R/niche_ellipse.R | 537 ++++++++++++++++++++++++++++++++++++---------- 1 file changed, 425 insertions(+), 112 deletions(-) diff --git a/R/niche_ellipse.R b/R/niche_ellipse.R index ff0d569..f32274c 100644 --- a/R/niche_ellipse.R +++ b/R/niche_ellipse.R @@ -138,138 +138,444 @@ niche_ellipse <- function( set.seed(set_seed) - sample_numbers <- sample(dat_mu$sample_number, n) + if (isotope_n == 2) { + sample_numbers <- sample(dat_mu$sample_number, n) + + isotope_a <- isotope_names[1] + isotope_b <- isotope_names[2] + # prepare mu for ellipse + mu <- dat_mu |> + dplyr::select(sample_name, sample_number, mu_est) |> + filter(sample_number %in% sample_numbers) |> + dplyr::group_split(sample_name, sample_number) |> + purrr::map(~ .x$mu_est, + .progress = "Prepare mu for ellipse") + + # preppare sigama for epplipse + isotope_a_sym <- rlang::sym(isotope_a) + isotope_b_sym <- rlang::sym(isotope_b) + # + # + # Erroring at isotope names fix will need make flexible qith {{{}}} + sigma <- dat_sigma |> + filter(sample_number %in% sample_numbers) |> + dplyr::select(sample_name, sample_number, !!isotope_a_sym, + !!isotope_b_sym) |> + dplyr::group_split(sample_name, sample_number) |> + purrr::map(~ cbind(.x[[rlang::as_name(isotope_a_sym)]], + .x[[rlang::as_name(isotope_b_sym)]]) |> + as.matrix(2, 2), .progress = "Prepare sigma for ellipse" + ) + + + # create ellipses + ellipse_dat <- purrr::pmap(list(sigma, mu), function(first, second) + ellipse::ellipse(x = first, + centre = second, + which = c(1, 2), + level = p_ell), + .progress = "Create ellipses") + + # Converting ellipse estimates into tibble + ellipse_dat <- ellipse_dat |> + purrr::map(~ .x |> + tibble::as_tibble(), + .progress = "Converting ellipse estimates into tibble" + ) + # create names to join name each ellimate of the list + list_names <- dat_sigma |> + filter(sample_number %in% sample_numbers) |> + dplyr::group_by(sample_name, sample_number) |> + dplyr::group_keys() |> + dplyr::mutate(sample_name_num = paste(sample_name, sample_number, + sep = ":")) + + names(ellipse_dat) <- list_names$sample_name_num + + # bind and rename columns + all_ellipses <- dplyr::bind_rows(ellipse_dat, .id = "ellipse_name") |> + dplyr::rename( + {{isotope_a}} := x, + {{isotope_b}} := y + ) |> + tidyr::separate(ellipse_name, + into = c("sample_name", "sample_number"), sep = ":") |> + dplyr::mutate( + sample_number = as.numeric(sample_number) + ) + + + end_time <- Sys.time() + + time_spent <- round((end_time - start_time), digits = 2) + if (message) { + cli::cli_alert(paste("Total time processing was", time_spent, units(time_spent), + sep = " ")) + } + + return(all_ellipses) + } + if (isotope_n == 3) { + + # create sample numbers + sample_numbers <- sample(dat_mu$sample_number, n) + + # ---- create every combo in table ---- + iso_combo <- tidyr::expand_grid(iso_a = isotope_names, + iso_b = isotope_names) |> + dplyr::filter(iso_a < iso_b) |> + dplyr::mutate( + iso_ab = paste(iso_a, iso_b, sep = "_") + ) |> + dplyr::distinct(iso_ab) |> + tidyr::separate_wider_delim(iso_ab, names = c("iso_a", "iso_b"), + delim = "_") + + # created id column + iso_combo$id <- 1:nrow(iso_combo) + + # pivot longer to create vector that can be used to filter combinations + split_iso <- split(iso_combo, iso_combo$id) |> + purrr::map(~ tidyr::pivot_longer(.x, cols = -id, + values_to = "iso_name") |> + getElement("iso_name") + ) + # ----- mu ---- + + # select columns we need + mu_select <- df_mu |> + dplyr::select(sample_name, sample_number, isotope, mu_est) |> + filter(sample_number %in% sample_numbers) + + # split into three different dataframes based on the iso combos + mu <- split_iso |> + purrr::map(~ mu_select |> + dplyr::filter(isotope %in% .x) |> + dplyr::group_split(sample_name, sample_number) + ) + + # extract the two est mu for each sample for each combo + mu_2 <- mu |> + purrr::map(~ purrr::map(.x, ~ .x$mu_est, + .progress = "Prepare mu for ellipse") + ) + + # ---- sigma ----- + + # transform to long as it is easier to keep the format + # required for the function + # to be wide + sigma_long <- df_sigma |> + tidyr::pivot_longer(cols = -c(metric, sample_name, sample_number, + isotope), + names_to = "id", + values_to = "est") |> + dplyr::filter(sample_number %in% sample_numbers) + + # filter based on split and then split into each group + sigma <- split_iso |> + purrr::map(~ sigma_long |> + dplyr::filter(id %in% .x & + isotope %in% .x) |> + tidyr::pivot_wider(names_from = "id", + values_from = "est") |> + dplyr::group_split(sample_name, sample_number) + ) + + # grab sample name and number + group_names <- split_iso |> + purrr::map(~ sigma_long |> + dplyr::filter(id %in% .x & + isotope %in% .x) |> + pivot_wider(names_from = "id", + values_from = "est") |> + dplyr::group_by(sample_name, sample_number) |> + dplyr::group_keys() |> + dplyr::mutate( + id = dplyr::row_number() |> + as.character() + ) + ) + + # empty lists to dump mattrix + sigma_list <- list() + + sm_list <- list() + + # for loop that converts every variance and covariance combo into + # matricies for ellipse + for (i in 1:length(sigma)) { + + st <- sigma[[i]] + + for (k in 1:length(st)) { + + sm <- st[[k]] |> + dplyr:: select(tidyselect::any_of(split_iso[[i]])) |> + as.matrix() + + sm_list[[k]] <- sm + } + sigma_list[[i]] <- sm_list + } + + # ----- make ellipses ----- + + ell_map <- map2( + .x = sigma_list, + .y = mu_2, + .f = map2, + function(x, y) ellipse::ellipse(x = x, + centre = y, + which = c(1, 2), + level = p_ell) |> + tibble::as_tibble(), + .progress = "Create ellipses" + ) - isotope_a <- isotope_name[1] - isotope_b <- isotope_name[2] - # prepare mu for ellipse - mu <- dat_mu |> - dplyr::select(sample_name, sample_number, mu_est) |> - filter(sample_number %in% sample_numbers) |> - dplyr::group_split(sample_name, sample_number) |> - purrr::map(~ .x$mu_est, - .progress = "Prepare mu for ellipse") - - # preppare sigama for epplipse - isotope_a_sym <- rlang::sym(isotope_a) - isotope_b_sym <- rlang::sym(isotope_b) - # - # - # Erroring at isotope names fix will need make flexible qith {{{}}} - sigma <- dat_sigma |> - filter(sample_number %in% sample_numbers) |> - dplyr::select(sample_name, sample_number, !!isotope_a_sym, - !!isotope_b_sym) |> - dplyr::group_split(sample_name, sample_number) |> - purrr::map(~ cbind(.x[[rlang::as_name(isotope_a_sym)]], - .x[[rlang::as_name(isotope_b_sym)]]) |> - as.matrix(2, 2), .progress = "Prepare sigma for ellipse" + # ---- add in sample names ---- + ell_sam <- map2( + .x = ell_map, + .y = group_names, + .f = function(x, y) x |> + dplyr::bind_rows(.id = "id") |> + dplyr::left_join(y, by = "id") |> + dplyr::select(-id) ) + # convert id of iso_combo to character for joining + iso_combo$id <- as.character(iso_combo$id) - # create ellipses - ellipse_dat <- purrr::pmap(list(sigma, mu), function(first, second) - ellipse::ellipse(x = first, - centre = second, - which = c(1, 2), - level = p_ell), - .progress = "Create ellipses") - # Converting ellipse estimates into tibble - ellipse_dat <- ellipse_dat |> - purrr::map(~ .x |> - tibble::as_tibble(), - .progress = "Converting ellipse estimates into tibble" - ) - # create names to join name each ellimate of the list - list_names <- dat_sigma |> - filter(sample_number %in% sample_numbers) |> - dplyr::group_by(sample_name, sample_number) |> - dplyr::group_keys() |> - dplyr::mutate(sample_name_num = paste(sample_name, sample_number, - sep = ":")) - - names(ellipse_dat) <- list_names$sample_name_num - - # bind and rename columns - all_ellipses <- dplyr::bind_rows(ellipse_dat, .id = "ellipse_name") |> - dplyr::rename( - {{isotope_a}} := x, - {{isotope_b}} := y - ) |> - tidyr::separate(ellipse_name, - into = c("sample_name", "sample_number"), sep = ":") |> - dplyr::mutate( - sample_number = as.numeric(sample_number) - ) + ell_final <- ell_sam |> + dplyr::bind_rows(.id = "id") |> + dplyr::left_join(iso_combo, by = "id") |> + dplyr::select(-id) |> + dplyr::mutate( + iso_combos = paste(iso_a, iso_b, sep = " - ") + ) |> + dplyr::select(sample_name:iso_combos, x, y) + end_time <- Sys.time() - end_time <- Sys.time() + time_spent <- round((end_time - start_time), digits = 2) + if (message) { + cli::cli_alert(paste("Total time processing was", time_spent, units(time_spent), + sep = " ")) + } + + return(ell_final) - time_spent <- round((end_time - start_time), digits = 2) - if (message) { - cli::cli_alert(paste("Total time processing was", time_spent, units(time_spent), - sep = " ")) - } - return(all_ellipses) + } } + # ---- put in random sample for 10 random samples if (random %in% FALSE) { # prepare mu for ellipse - mu <- dat_mu |> - dplyr::select(sample_name, sample_number, mu_est) |> - dplyr::group_split(sample_name, sample_number) |> - purrr::map(~ .x$mu_est, - .progress = "Prepare mu for ellipse") - - # preppare sigama for epplipse - # - # Erroring at isotope names fix will need make flexible qith {{{}}} - sigma <- dat_sigma |> - dplyr::select(sample_name, sample_number, d13c, d15n) |> - dplyr::group_split(sample_name, sample_number) |> - purrr::map(~ cbind(.x$d13c, .x$d15n) |> - as.matrix(2, 2), .progress = "Prepare sigma for ellipse" + if (isotope_n == 2) { + + isotope_a <- isotope_names[1] + isotope_b <- isotope_names[2] + + mu <- dat_mu |> + dplyr::select(sample_name, sample_number, mu_est) |> + dplyr::group_split(sample_name, sample_number) |> + purrr::map(~ .x$mu_est, + .progress = "Prepare mu for ellipse") + + # preppare sigama for epplipse + # + # Erroring at isotope names fix will need make flexible qith {{{}}} + sigma <- dat_sigma |> + dplyr::select(sample_name, sample_number, d13c, d15n) |> + dplyr::group_split(sample_name, sample_number) |> + purrr::map(~ cbind(.x$d13c, .x$d15n) |> + as.matrix(2, 2), .progress = "Prepare sigma for ellipse" + ) + + + # create ellipses + ellipse_dat <- purrr::pmap(list(sigma, mu), function(first, second) + ellipse::ellipse(x = first, + centre = second, + which = c(1, 2), + level = p_ell), + .progress = "Create ellipses") + + # Converting ellipse estimates into tibble + ellipse_dat <- ellipse_dat |> + purrr::map(~ .x |> + tibble::as_tibble(), + .progress = "Converting ellipse estimates into tibble" + ) + # create names to join name each ellimate of the list + list_names <- dat_sigma |> + dplyr::group_by(sample_name, sample_number) |> + dplyr::group_keys() |> + dplyr::mutate(sample_name_num = paste(sample_name, sample_number, + sep = ":")) + + names(ellipse_dat) <- list_names$sample_name_num + + # bind and rename columns + all_ellipses <- dplyr::bind_rows(ellipse_dat, .id = "ellipse_name") |> + dplyr::rename( + {{isotope_a}} := x, + {{isotope_b}} := y + ) |> + tidyr::separate_wider_delim(ellipse_name, + names = c("sample_name", "sample_number"), delim = ":") |> + dplyr::mutate( + sample_number = as.numeric(sample_number) + ) + end_time <- Sys.time() + + time_spent <- round((end_time - start_time), digits = 2) + if (message) { + cli::cli_alert(paste("Total time processing was", time_spent, units(time_spent), + sep = " ")) + } + + return(all_ellipses) + + } + if (isotope_n == 3) { + + # ---- create every combo in table ---- + iso_combo <- tidyr::expand_grid(iso_a = isotope_names, + iso_b = isotope_names) |> + dplyr::filter(iso_a < iso_b) |> + dplyr::mutate( + iso_ab = paste(iso_a, iso_b, sep = "_") + ) |> + dplyr::distinct(iso_ab) |> + tidyr::separate_wider_delim(iso_ab, names = c("iso_a", "iso_b"), + delim = "_") + + # created id column + iso_combo$id <- 1:nrow(iso_combo) + + # pivot longer to create vector that can be used to filter combinations + split_iso <- split(iso_combo, iso_combo$id) |> + purrr::map(~ tidyr::pivot_longer(.x, cols = -id, + values_to = "iso_name") |> + getElement("iso_name") ) + # ----- mu ---- + # select columns we need + mu_select <- df_mu |> + dplyr::select(sample_name, sample_number, isotope, mu_est) - # create ellipses - ellipse_dat <- purrr::pmap(list(sigma, mu), function(first, second) - ellipse::ellipse(x = first, - centre = second, - which = c(1, 2), - level = p_ell), - .progress = "Create ellipses") + # split into three different dataframes based on the iso combos + mu <- split_iso |> + purrr::map(~ mu_select |> + dplyr::filter(isotope %in% .x) |> + dplyr::group_split(sample_name, sample_number) + ) - # Converting ellipse estimates into tibble - ellipse_dat <- ellipse_dat |> - purrr::map(~ .x |> - tibble::as_tibble(), - .progress = "Converting ellipse estimates into tibble" + # extract the two est mu for each sample for each combo + mu_2 <- mu |> + purrr::map(~ purrr::map(.x, ~ .x$mu_est, + .progress = "Prepare mu for ellipse") ) - # create names to join name each ellimate of the list - list_names <- dat_sigma |> - dplyr::group_by(sample_name, sample_number) |> - dplyr::group_keys() |> - dplyr::mutate(sample_name_num = paste(sample_name, sample_number, - sep = ":")) - - names(ellipse_dat) <- list_names$sample_name_num - - # bind and rename columns - all_ellipses <- dplyr::bind_rows(ellipse_dat, .id = "ellipse_name") |> - dplyr::rename( - {{isotope_a}} := x, - {{isotope_b}} := y - ) |> - tidyr::separate(ellipse_name, - into = c("sample_name", "sample_number"), sep = ":") |> - dplyr::mutate( - sample_number = as.numeric(sample_number) + + # ---- sigma ----- + + # transform to long as it is easier to keep the format + # required for the function + # to be wide + sigma_long <- df_sigma |> + tidyr::pivot_longer(cols = -c(metric, sample_name, sample_number, + isotope), + names_to = "id", + values_to = "est") + + # filter based on split and then split into each group + sigma <- split_iso |> + purrr::map(~ sigma_long |> + dplyr::filter(id %in% .x & + isotope %in% .x) |> + tidyr::pivot_wider(names_from = "id", + values_from = "est") |> + dplyr::group_split(sample_name, sample_number) ) + # grab sample name and number + group_names <- split_iso |> + purrr::map(~ sigma_long |> + dplyr::filter(id %in% .x & + isotope %in% .x) |> + pivot_wider(names_from = "id", + values_from = "est") |> + dplyr::group_by(sample_name, sample_number) |> + dplyr::group_keys() |> + dplyr::mutate( + id = dplyr::row_number() |> + as.character() + ) + ) + + # empty lists to dump mattrix + sigma_list <- list() + + sm_list <- list() + + # for loop that converts every variance and covariance combo into + # matricies for ellipse + for (i in 1:length(sigma)) { + + st <- sigma[[i]] + + for (k in 1:length(st)) { + + sm <- st[[k]] |> + dplyr:: select(tidyselect::any_of(split_iso[[i]])) |> + as.matrix() + + sm_list[[k]] <- sm + } + sigma_list[[i]] <- sm_list + } + + # ----- make ellipses ----- + + ell_map <- map2( + .x = sigma_list, + .y = mu_2, + .f = map2, + function(x, y) ellipse::ellipse(x = x, + centre = y, + which = c(1, 2), + level = p_ell) |> + tibble::as_tibble(), + .progress = "Create ellipses" + ) + + # ---- add in sample names ---- + ell_sam <- map2( + .x = ell_map, + .y = group_names, + .f = function(x, y) x |> + dplyr::bind_rows(.id = "id") |> + dplyr::left_join(y, by = "id") |> + dplyr::select(-id) + ) + + # convert id of iso_combo to character for joining + iso_combo$id <- as.character(iso_combo$id) + + + ell_final <- ell_sam |> + dplyr::bind_rows(.id = "id") |> + dplyr::left_join(iso_combo, by = "id") |> + dplyr::select(-id) |> + dplyr::mutate( + iso_combos = paste(iso_a, iso_b, sep = " - ") + ) |> + dplyr::select(sample_name:iso_combos, x, y) end_time <- Sys.time() @@ -279,7 +585,14 @@ niche_ellipse <- function( sep = " ")) } - return(all_ellipses) + return(ell_final) + + + # return(all_ellipses) + } + } + + } From 72acee9b234b5a0c40551290e088347facf9c2ec Mon Sep 17 00:00:00 2001 From: benjaminhlina Date: Tue, 12 Nov 2024 14:30:20 -0500 Subject: [PATCH 21/35] Updated tests to include isotope_names not isotope_a and b --- tests/testthat/test-niche_ellipse.R | 14 ++------------ 1 file changed, 2 insertions(+), 12 deletions(-) diff --git a/tests/testthat/test-niche_ellipse.R b/tests/testthat/test-niche_ellipse.R index 47888f2..5370d7a 100644 --- a/tests/testthat/test-niche_ellipse.R +++ b/tests/testthat/test-niche_ellipse.R @@ -89,8 +89,7 @@ test_that("Check if naming works", { dat_sigma = sigma_est_wide, set_seed = 4, message = FALSE, - isotope_a = "cal_d13c", - isotope_b = "cal_d15n", + isotope_names = c("cal_d13c", "cal_d15n"), ) ) @@ -110,19 +109,10 @@ test_that("Test if isotope_a errors if not a charcters", { dat_mu = mu_est_long, dat_sigma = sigma_est_wide, message = FALSE, - isotope_a = 6), regexp = "Argument 'isotope_a' must be a character." + isotope_names = 6), regexp = "Argument 'isotope_names' must be a character." ) }) -test_that("Test if isotope_b errors if not a charcters", { - expect_error(niche_ellipse( - dat_mu = mu_est_long, - dat_sigma = sigma_est_wide, - set_seed = 4, - message = FALSE, - isotope_b = 4), regexp = "Argument 'isotope_b' must be a character." - ) -}) test_that("Parameter 'p_ell' is can take other values than 0.95", { From 7c8563bae537c790d67647d6ac672aa40d88d74d Mon Sep 17 00:00:00 2001 From: benjaminhlina Date: Tue, 12 Nov 2024 17:37:39 -0500 Subject: [PATCH 22/35] Added tidyselect --- DESCRIPTION | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 6802009..0b18e7c 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -27,7 +27,8 @@ Imports: rlang, SIBER, tibble, - tidyr + tidyr, + tidyselect License: CC0 Encoding: UTF-8 LazyData: true From 51f622de80ff7a873f9f3d51a5ebcab7e96d8748 Mon Sep 17 00:00:00 2001 From: benjaminhlina Date: Tue, 12 Nov 2024 17:37:49 -0500 Subject: [PATCH 23/35] Added tidyselect --- NAMESPACE | 1 + 1 file changed, 1 insertion(+) diff --git a/NAMESPACE b/NAMESPACE index bfe1bb9..7347d72 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -15,6 +15,7 @@ import(ellipse) import(purrr) import(tibble) import(tidyr) +import(tidyselect) importFrom(lifecycle,deprecated) importFrom(nicheROVER,niche.size) importFrom(rlang,":=") From 907f17aced2059b266a23ef5d5aebd97d9e080de Mon Sep 17 00:00:00 2001 From: benjaminhlina Date: Tue, 12 Nov 2024 17:38:29 -0500 Subject: [PATCH 24/35] Fixed typo in documents --- R/niche_ellipse.R | 2 +- man/niche_ellipse.Rd | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/R/niche_ellipse.R b/R/niche_ellipse.R index f32274c..0ddfabd 100644 --- a/R/niche_ellipse.R +++ b/R/niche_ellipse.R @@ -15,7 +15,7 @@ #' Default is 0.95 (i.e., 95% confidence interval). #' This value is bound by 0 and 1 and has to be a `numeric`. #' @param isotope_n a `numeric` either `2` or `3` that is the number of isotopes -#' used in the anlsysis. Will default to `2`. +#' used in the analysis. Will default to `2`. #' @param isotope_names is a vector of `character` string used change the column name #' of isotopes used in the analysis. Defaults to `c("d13c", "d15n")`. #' @param random logical value indicating whether or not to randomly sample diff --git a/man/niche_ellipse.Rd b/man/niche_ellipse.Rd index 3426a67..22698c3 100644 --- a/man/niche_ellipse.Rd +++ b/man/niche_ellipse.Rd @@ -28,7 +28,7 @@ matrices stacked on top of each other. See example of how to convert to wide format. This can be produced using \code{extract_sigma()}.} \item{isotope_n}{a \code{numeric} either \code{2} or \code{3} that is the number of isotopes -used in the anlsysis. Will default to \code{2}.} +used in the analysis. Will default to \code{2}.} \item{isotope_names}{is a vector of \code{character} string used change the column name of isotopes used in the analysis. Defaults to \code{c("d13c", "d15n")}.} From 0532baaded0347570f06d3f799cb7163f902a335 Mon Sep 17 00:00:00 2001 From: benjaminhlina Date: Tue, 12 Nov 2024 17:38:47 -0500 Subject: [PATCH 25/35] added tidyselect --- R/niche_ellipse.R | 1 + 1 file changed, 1 insertion(+) diff --git a/R/niche_ellipse.R b/R/niche_ellipse.R index 0ddfabd..3ce57df 100644 --- a/R/niche_ellipse.R +++ b/R/niche_ellipse.R @@ -48,6 +48,7 @@ #' @importFrom stats runif #' @import tibble #' @import tidyr +#' @import tidyselect #' @export From 6c614ab72c7a12e1b637b883423cd4df5559fc72 Mon Sep 17 00:00:00 2001 From: benjaminhlina Date: Tue, 12 Nov 2024 17:39:06 -0500 Subject: [PATCH 26/35] Fixed df_ to dat_ --- R/niche_ellipse.R | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/R/niche_ellipse.R b/R/niche_ellipse.R index 3ce57df..4070c04 100644 --- a/R/niche_ellipse.R +++ b/R/niche_ellipse.R @@ -244,7 +244,7 @@ niche_ellipse <- function( # ----- mu ---- # select columns we need - mu_select <- df_mu |> + mu_select <- dat_mu |> dplyr::select(sample_name, sample_number, isotope, mu_est) |> filter(sample_number %in% sample_numbers) @@ -266,7 +266,7 @@ niche_ellipse <- function( # transform to long as it is easier to keep the format # required for the function # to be wide - sigma_long <- df_sigma |> + sigma_long <- dat_sigma |> tidyr::pivot_longer(cols = -c(metric, sample_name, sample_number, isotope), names_to = "id", @@ -467,7 +467,7 @@ niche_ellipse <- function( # ----- mu ---- # select columns we need - mu_select <- df_mu |> + mu_select <- dat_mu |> dplyr::select(sample_name, sample_number, isotope, mu_est) # split into three different dataframes based on the iso combos @@ -488,7 +488,7 @@ niche_ellipse <- function( # transform to long as it is easier to keep the format # required for the function # to be wide - sigma_long <- df_sigma |> + sigma_long <- dat_sigma |> tidyr::pivot_longer(cols = -c(metric, sample_name, sample_number, isotope), names_to = "id", From 478767274505276d477c18e517533847466035e3 Mon Sep 17 00:00:00 2001 From: benjaminhlina Date: Tue, 12 Nov 2024 17:39:39 -0500 Subject: [PATCH 27/35] Added new variables --- R/zzz.R | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/R/zzz.R b/R/zzz.R index 659d506..d3f0e4a 100644 --- a/R/zzz.R +++ b/R/zzz.R @@ -12,4 +12,6 @@ utils::globalVariables(c("sample_name", "sample_number", "metric", "x", "y", "overlap", "area_1", "area_2", "niche_overlap_perc", "niche_overlap", "species_a", "species_b", "isotope", - "post_sample", "mu_est", "ellipse_name")) + "post_sample", "mu_est", "ellipse_name", + "iso_a", "iso_ab", "iso_b", "iso_combos" + )) From 87afefa1f177c02d84f1b8ab1fd9c479a3325537 Mon Sep 17 00:00:00 2001 From: benjaminhlina Date: Tue, 12 Nov 2024 17:48:34 -0500 Subject: [PATCH 28/35] Update version number to 0.3.3 --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 0b18e7c..8c482d0 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: nichetools Type: Package Title: Complementary Package to 'nicheROVER' and 'SIBER' -Version: 0.3.2 +Version: 0.3.3 Authors@R: c(person(given = "Benjamin L.", family = "Hlina", From 62c96c5dea65203b52652fe3d0a01e4ae008b4c5 Mon Sep 17 00:00:00 2001 From: benjaminhlina Date: Tue, 12 Nov 2024 17:48:47 -0500 Subject: [PATCH 29/35] Update package version --- R/zzz.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/zzz.R b/R/zzz.R index d3f0e4a..d5f5a89 100644 --- a/R/zzz.R +++ b/R/zzz.R @@ -1,5 +1,5 @@ .onAttach <- function(libname, pkgname) { - packageStartupMessage("version 0.3.2 ('fall-fieldwork').\nHave you loaded {nicheROVER} or {SIBER}? If not, please do so.") + packageStartupMessage("version 0.3.3 ('let-the-snow-fly').\nHave you loaded {nicheROVER} or {SIBER}? If not, please do so.") } utils::globalVariables(c("sample_name", "sample_number", "metric", "x", "y", "V1", "V2", From c7aba11b4776307f3d5d208e205c4af04d77a1e2 Mon Sep 17 00:00:00 2001 From: benjaminhlina Date: Wed, 18 Dec 2024 11:28:52 -0500 Subject: [PATCH 30/35] Export isotope_pairs --- NAMESPACE | 1 + 1 file changed, 1 insertion(+) diff --git a/NAMESPACE b/NAMESPACE index 7347d72..98d868d 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,6 +1,7 @@ # Generated by roxygen2: do not edit by hand export(create_comparisons) +export(create_isotope_pairs) export(extract_group_metrics) export(extract_layman) export(extract_mu) From 70a8693c111507be5af9d48eacc9753c6a937bf9 Mon Sep 17 00:00:00 2001 From: benjaminhlina Date: Wed, 18 Dec 2024 11:29:29 -0500 Subject: [PATCH 31/35] create isotope_pairs function --- R/create_isotope_pairs.R | 55 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 55 insertions(+) create mode 100644 R/create_isotope_pairs.R diff --git a/R/create_isotope_pairs.R b/R/create_isotope_pairs.R new file mode 100644 index 0000000..4f05bf6 --- /dev/null +++ b/R/create_isotope_pairs.R @@ -0,0 +1,55 @@ +#' Create isotope pairs +#' +#' @param isotope_n number of isotopes +#' @param isotope_names names of isotopes +#' +#' @export + +create_isotope_pairs <- function(isotope_n = NULL, + isotope_names = NULL) { + + # create name vector that will be used to id isotopes. + if (is.null(isotope_n)) { + isotope_n <- 3 + } + if (!is.numeric(isotope_n) || !(isotope_n %in% c(3))) { + cli::cli_abort("Argument 'isotope_n' must 3.") + } + + if (isotope_n == 3) { + if (is.null(isotope_names)) { + isotope_names <- c("d13c", "d15n", "d34s") + } + } + if (!is.character(isotope_names)) { + cli::cli_abort("Argument 'isotope_names' must be a character.") + } + + # ---- create every combo in table ---- + iso_combo <- tidyr::expand_grid(iso_a = isotope_names, + iso_b = isotope_names) |> + dplyr::filter(iso_a < iso_b) |> + dplyr::mutate( + iso_ab = paste(iso_a, iso_b, sep = "_") + ) |> + dplyr::distinct(iso_ab) |> + tidyr::separate_wider_delim(iso_ab, names = c("iso_a", "iso_b"), + delim = "_") + + # created id column + iso_combo$id <- 1:nrow(iso_combo) + + # pivot longer to create vector that can be used to filter combinations + split_iso <- split(iso_combo, iso_combo$id) |> + purrr::map(~ tidyr::pivot_longer(.x, cols = -id, + values_to = "iso_name") |> + dplyr::mutate( + iso_name = factor(iso_name, level = c("d34s", + "d13c", + "d15n")) + ) |> + dplyr::arrange(iso_name) |> + getElement("iso_name") + ) + return(split_iso) +} From f73c4dc3a80238689723f2dc5862b3ae59eb88b9 Mon Sep 17 00:00:00 2001 From: benjaminhlina Date: Wed, 18 Dec 2024 11:29:47 -0500 Subject: [PATCH 32/35] Updated documentatyion --- man/create_isotope_pairs.Rd | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) create mode 100644 man/create_isotope_pairs.Rd diff --git a/man/create_isotope_pairs.Rd b/man/create_isotope_pairs.Rd new file mode 100644 index 0000000..37e195f --- /dev/null +++ b/man/create_isotope_pairs.Rd @@ -0,0 +1,16 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/create_isotope_pairs.R +\name{create_isotope_pairs} +\alias{create_isotope_pairs} +\title{Create isotope pairs} +\usage{ +create_isotope_pairs(isotope_n = NULL, isotope_names = NULL) +} +\arguments{ +\item{isotope_n}{number of isotopes} + +\item{isotope_names}{names of isotopes} +} +\description{ +Create isotope pairs +} From 5df73d1895483d8d6679b4c0b7938c5e0ff4718b Mon Sep 17 00:00:00 2001 From: YourName Date: Wed, 18 Dec 2024 11:34:51 -0500 Subject: [PATCH 33/35] Create vignettes on three isotopes --- ...ols-and-nicheROVER-with-three-isotopes.Rmd | 532 ++++++++++++++++++ 1 file changed, 532 insertions(+) create mode 100644 vignettes/using-nichetools-and-nicheROVER-with-three-isotopes.Rmd diff --git a/vignettes/using-nichetools-and-nicheROVER-with-three-isotopes.Rmd b/vignettes/using-nichetools-and-nicheROVER-with-three-isotopes.Rmd new file mode 100644 index 0000000..3f1571f --- /dev/null +++ b/vignettes/using-nichetools-and-nicheROVER-with-three-isotopes.Rmd @@ -0,0 +1,532 @@ +--- +title: "Using {nichetools} and {nicheROVER} with three isotopes" +author: Benjamin L. Hlina +date: "2024-11-13" +output: rmarkdown::html_vignette + +vignette: > + %\VignetteIndexEntry{Using {nichetools} and {nicheROVER} with three isotopes} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r, include = FALSE} + knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>" + ) +``` + +### Our Objectives + +The purpose of this vignette is to use [{nicheROVER}](https://cran.r-project.org/package=nicheROVER) and [{nichetools}](https://benjaminhlina.github.io/nichetools/) to extract and then visualize estimates of isotopic niche size and similarities using three isotopes (i.e., $\delta$13C, $\delta$15N, and $\delta$34S) for multiple freshwater fish using [{ggplot2}](https://ggplot2.tidyverse.org/). + +This vignette can be used for additional purposes including estimating niche size and similarities among different groups of aquatic and/or terrestrial species. Furthermore, niche size and similarities for different behaviours exhibited within a population can be made using behavioural data generated from acoustic telemetry (e.g., differences in habitat occupancy). + + +### Bring in trophic niche data +First we will load the necessary packages to preform the analysis and visualization. We will use [{nicheROVER}](https://cran.r-project.org/package=nicheROVER) and [{nichetools}](https://benjaminhlina.github.io/nichetools/) to preform the analysis. We will use [{dplyr}](https://dplyr.tidyverse.org/), [{tidyr}](https://tidyr.tidyverse.org/), and [{purrr}](https://purrr.tidyverse.org/) to manipulate data and iterate processes. Lastly, we will use [{ggplot2}](https://ggplot2.tidyverse.org/), [{ggtext}](https://wilkelab.org/ggtext/), and [{patchwork}](https://patchwork.data-imaginist.com/) to plot, add labels, and arrange plots. + + I will add that many of the `{dplyr}` and `{tidyr}` functions and processes can be replaced using [{data.table}](https://CRAN.R-project.org/package=data.table) which is great when working with large data sets. + +```{r, message = FALSE} + { + library(dplyr) + library(ggplot2) + library(ggtext) + library(ggh4x) + library(nicheROVER) + library(nichetools) + library(patchwork) + library(purrr) + library(stringr) + library(tidyr) + } +``` + + For the purpose of the vignette we will be using the `fish` data frame that is available within `{nicheROVER}`. + + We will first use the function `janitor::clean_names()` to clean up column names. For your purposes you will need to replace fish with your data frame either by loading a csv, rds, or qs, with your data. You can do this multiple ways, I prefer using `readr::read_csv()` but base R's `read.csv()` works perfectly fine. + +```{r, message = FALSE} + df <- fish %>% + janitor::clean_names() +``` + +If there are any isotopic values that did not run and are `NA`, they will need to be removed because functions in `{nicheROVER}` will not accommodate values of `NA`. + + ### Estimate posterior distribution with Normal-Inverse-Wishart (NIW) priors. + + We will take 1,000 posterior samples for each group. You can change this but suggest nothing less than 1,000. + +```{r, message = FALSE} + nsample <- 1000 +``` + +We will then split the data frame into a list with each species as a data frame object within the list, We will then iterate over the list, using `map()` from [{purrr}](https://purrr.tidyverse.org/), to estimate posterior distribution using Normal-Inverse-Wishart (NIW) priors. + +You will notice that $\delta$34S has been selected for this analysis which differs compared to the other vignette focused on using `{nichetools}` with `{nicheROVER}` for two isotopes. + +```{r, message = FALSE} +fish_par <- df %>% + split(.$species) %>% + map(~ select(., d13c, d15n, d34s)) %>% + map(~ niw.post(nsample = nsample, X = .)) +``` + +### Extract μ values + +We will use `extract_mu()`to extract posteriors for $\mu$ estimates. The default output of `extract_mu()` is long format +which works for plotting with {ggplot2} and other functions in {nichetools}. If we want wide format we can specify the argument `format` with `"wide"`, however, it is unlikely you will need this data in wide format. + +```{r, message = FALSE} +df_mu <- extract_mu(fish_par) +``` + +The default output will be lacking some info for plotting. We will need to add in a column that is the element abbreviation and neutron number to be used in axis labeling. + +```{r} +df_mu <- df_mu %>% + mutate( + element = case_when( + isotope == "d15n" ~ "N", + isotope == "d13c" ~ "C", + isotope == "d34s" ~ "S", + ), + neutron = case_when( + isotope == "d15n" ~ 15, + isotope == "d13c" ~ 13, + isotope == "d34s" ~ 34, + ) + ) +``` + +### Extract Σ values +We will use `extract_sigma()` to extract posterior estimates for $\Sigma$. The default output of `extract_sigma()` is wide format which doesn't work for plotting with {ggplot2} but does work other functions in {nichetools}. If we want long for plotting we can specify the argument `format` with `"long"`. + +Remember to change the argument `isotope_n` from `2` to `3`, considering we are working with three isotopes in this vignette. + +```{r, message = FALSE} +df_sigma <- extract_sigma(fish_par, isotope_n = 3) +``` + +For plotting we will need the extracted $\Sigma$ values to be in long format. We also need to remove $\Sigma$ values for when the both isotope columns are the same isotope. + +```{r, message = FALSE} +df_sigma_cn <- extract_sigma(fish_par, + data_format = "long", + isotope_n = 3 +) %>% + filter(id != isotope) +``` + + +### Plot posterior distrubtion of μ and Σ + +For most plotting within this vignette, I will `split()` the data frame by isotope, creating a list that I will then use `imap()` to iterate over the list to create plots. We will use `geom_density()` to represent densities for both $\mu$ and $\Sigma$. Plot objects will then be stored in a list. + +First we will plot $\mu$ for each isotope. We will use [{patchwork}](https://patchwork.data-imaginist.com/) to configure plots for multi-panel figures. This package is phenomenal and uses math operators to configure and manipulate the plots to create multi-panel figures. + +For labeling we are also going to use `element_markdown()` from [{ggtext}](https://wilkelab.org/ggtext/) to work with the labels that are needed to correctly display the isotopic signature. If you are working other data please replace. + +```{r, warning = FALSE} +#| out-width: 100% +posterior_plots <- df_mu %>% + split(.$isotope) %>% + imap( + ~ ggplot(data = ., aes(x = mu_est)) + + geom_density(aes(fill = sample_name), alpha = 0.5) + + scale_fill_viridis_d(begin = 0.25, end = 0.75, + option = "D", name = "Species") + + theme_bw() + + theme(panel.grid = element_blank(), + axis.title.x = element_markdown(), + axis.title.y = element_markdown(), + legend.position = "none", + legend.background = element_blank() + ) + + labs( + x = paste("\u00b5\U03B4", "", + unique(.$neutron), "", + "",unique(.$element), "", sep = ""), + y = paste0("p(\u00b5 \U03B4","", + unique(.$neutron), "", + "",unique(.$element),"", + " | X)"), sep = "") + ) + +posterior_plots$d15n + + theme(legend.position = c(0.18, 0.82)) + + posterior_plots$d13c + + posterior_plots$d34s +``` + + +For labeling purposes we need to add columns that are the element abbreviation and neutron number. I do this by using `case_when()` which are vectorized if else statements. + +```{r, message = FALSE} +df_sigma_cn <- df_sigma_cn %>% + mutate( + element_id = case_when( + id == "d15n" ~ "N", + id == "d13c" ~ "C", + id == "d34s" ~ "S", + ), + neutron_id = case_when( + id == "d15n" ~ 15, + id == "d13c" ~ 13, + id == "d34s" ~ 34, + ), + element_iso = case_when( + isotope == "d15n" ~ "N", + isotope == "d13c" ~ "C", + isotope == "d34s" ~ "S", + ), + neutron_iso = case_when( + isotope == "d15n" ~ 15, + isotope == "d13c" ~ 13, + isotope == "d34s" ~ 34, + ) + ) +``` + +Next we will plot the posteriors for $\Sigma$. +```{r, message = FALSE} +#| out-width: 100% +sigma_plots <- df_sigma_cn %>% + group_split(id, isotope) %>% + imap( + ~ ggplot(data = ., aes(x = post_sample)) + + geom_density(aes(fill = sample_name), alpha = 0.5) + + scale_fill_viridis_d(begin = 0.25, end = 0.75, + option = "D", name = "Species") + + theme_bw() + + theme(panel.grid = element_blank(), + axis.title.x = element_markdown(), + axis.title.y = element_markdown(), + legend.position = "none" + ) + + labs( + x = paste("\U03A3","\U03B4", + "", unique(.$neutron_id), "", + "",unique(.$element_id),""," ", + "\U03B4", + "", unique(.$neutron_iso), "", + "",unique(.$element_iso),"", sep = ""), + y = paste("p(", "\U03A3","\U03B4", + "", unique(.$neutron_id), "", + "",unique(.$element_id),""," ", + "\U03B4", + "", unique(.$neutron_iso), "", + "",unique(.$element_iso),"", " | X)", sep = ""), + ) + ) + +sigma_plots[[1]] + + theme(legend.position = c(0.1, 0.82)) + + sigma_plots[[2]] + + sigma_plots[[4]] + +``` + + +### Estimate niche ellipse + +We then will use `niche_ellipse()` to easily extract ellipse for each $\Sigma$ estimate (i.e., 1000). The function will tell you how long it took to process as with large sets of isotope data it is nice to know the time it takes for the function to work. The function also has the argument `random` which by default is set to `TRUE`. This argument randomly subsamples and returns 10 ellipse estimates out of the total number of samples taken in this case it is 1,000. The `set_seed` argument allows you to change the `set.seed` value by giving a numerical value to make the results of the function reproducible, default is a random value. I highly suggest using set_seed otherwise it will subsample different values, there is no default value because CRAN will not allow for a default value. You can change the number of subsamples but 10 seems pretty standard. If you'd like all 1,000 ellipses you can set `random` to `FALSE`. + +Remember to change the argument `isotope_n` from `2` to `3`, considering we are working with three isotopes in this vignette. + + +```{r} +ellipse_df <- niche_ellipse(dat_mu = df_mu, + dat_sigma = df_sigma, + isotope_n = 3, + set_seed = 4) +``` + +```{r, message = FALSE} +ellipse_df <- ellipse_df %>% + mutate( + element_id = case_when( + iso_a == "d15n" ~ "N", + iso_a == "d13c" ~ "C", + iso_a == "d34s" ~ "S", + ), + neutron_id = case_when( + iso_a == "d15n" ~ 15, + iso_a == "d13c" ~ 13, + iso_a == "d34s" ~ 34, + ), + element_iso = case_when( + iso_b == "d15n" ~ "N", + iso_b == "d13c" ~ "C", + iso_b == "d34s" ~ "S", + ), + neutron_iso = case_when( + iso_b == "d15n" ~ 15, + iso_b == "d13c" ~ 13, + iso_b == "d34s" ~ 34, + ) + ) +``` + + +### Plot ellipses, densities of each istope, and isotope biplot + +We will first plot the ellipse for each sample_name +```{r} +ellipse_plots <- + ellipse_df %>% + split(.$iso_combos) %>% + map(~ .x %>% + ggplot() + + geom_polygon(data = .x, + mapping = aes(x = x, y = y, + group = interaction(sample_number, sample_name), + color = sample_name), + fill = NA, + linewidth = 0.5) + + + scale_colour_viridis_d(begin = 0.25, end = 0.75, + option = "D", name = "species", + ) + + theme_bw(base_size = 10) + + theme(axis.text = element_text(colour = "black"), + axis.title.x = element_markdown(), + axis.title.y = element_markdown(), + panel.grid = element_blank(), + legend.position = "none", + legend.title = element_text(hjust = 0.5), + legend.background = element_blank() + ) + + labs( + x = paste("\U03B4","", unique(.$neutron_id), "", + unique(.$element_id), sep = ""), + y = paste("\U03B4","", unique(.$neutron_iso), "", + unique(.$element_iso), sep = ""), + ) + ) +``` +We need to turn `df` into long format to iterate over using `imap()` to easily create density plots. You will notice that I again use `case_when()` to make columns of element abbreviations and neutron numbers that will be used in plot labeling. + +```{r} +iso_long <- df %>% + pivot_longer(cols = -species, + names_to = "isotope", + values_to = "value") %>% + mutate( + element = case_when( + isotope == "d15n" ~ "N", + isotope == "d13c" ~ "C", + isotope == "d34s" ~ "S", + ), + neutron = case_when( + isotope == "d15n" ~ 15, + isotope == "d13c" ~ 13, + isotope == "d34s" ~ 34, + ) + ) +``` + +We will then make density plots for each isotope using `geom_density()` +```{r, warning=FALSE} +iso_density <- iso_long %>% + group_split(isotope) %>% + imap( + ~ ggplot(data = .) + + geom_density(aes(x = value, + fill = species), + alpha = 0.35, + linewidth = 0.8) + + scale_fill_viridis_d(begin = 0.25, end = 0.75, + option = "D", name = "Species") + + theme_bw(base_size = 10) + + theme(axis.text = element_text(colour = "black"), + panel.grid = element_blank(), + legend.position = c(0.15, 0.55), + legend.background = element_blank(), + axis.title.x = element_markdown(family = "sans")) + + labs(x = paste("\U03B4", + "", unique(.$neutron), "",unique(.$element), + sep = ""), + y = "Density") + ) + +d13c_density <- iso_density[[1]] + + scale_x_continuous(breaks = rev(seq(-20, -34, -2)), + limits = rev(c(-20, -34))) + +d15n_density <- iso_density[[2]] + + scale_x_continuous(breaks = seq(5, 15, 2.5), + limits = c(5, 15)) + + theme( + legend.position = "none" + ) + +d34s_density <- iso_density[[3]] + + theme( + legend.position = "none" + ) + +``` + + +```{r, warning=FALSE} +split_iso <- create_isotope_pairs() + +iso_split <- split_iso |> + purrr::map(~ iso_long |> + mutate( + isotope = factor(isotope, level = c("d34s", + "d13c", + "d15n")) + ) |> + arrange(isotope) |> + filter(isotope %in% .x) |> + pivot_wider(id_cols = species, + names_from = "isotope", + values_from = c("value", "element", "neutron")) |> + unnest() + ) + +``` + +Lastly we will use `geom_point()` to make isotopic biplot. + +```{r} +iso_biplot <- iso_split %>% + purrr::map(~ ggplot(data = .x, aes(x = .x[[2]], y = .x[[3]])) + + geom_point(aes(fill = .x$species), size = 3, + shape = 21, colour = "black", + stroke = 0.8, alpha = 0.70) + + scale_fill_viridis_d(begin = 0.25, end = 0.75, + option = "D", name = "species") + + theme_bw(base_size = 10) + + theme(axis.text = element_text(colour = "black"), + axis.title.x = element_markdown(), + axis.title.y = element_markdown(), + panel.grid = element_blank(), + legend.position = "none", + legend.title = element_text(hjust = 0.5), + legend.background = element_blank() + ) + + labs( + x = paste("\U03B4","", unique(.x[[6]]), "", + unique(.x[[4]]), sep = ""), + y = paste("\U03B4","", unique(.x[[7]]), "", + unique(.x[[5]]), sep = ""), + ) + ) +``` + +### Use {patchwork} to make ellipse, density, and biplots into a paneled figure. + +We can also use the function `plot_annotation()` to add lettering to the figure that can be used in the figure description. To maneuver where `plot_annotation()` places the lettering, we need to add `plot.tag.position = c(x, y)` to the `theme()` call in every plot. + +```{r} +#| out-width: 100% +d13c_density + ellipse_plots[[1]] + ellipse_plots[[2]] + + iso_biplot[[1]] + d15n_density + ellipse_plots[[3]] + + iso_biplot[[2]] + iso_biplot[[3]] + d34s_density + + plot_annotation(tag_levels = "a", + tag_suffix = ")") +``` + +### Determine the 95% niche similarties for each species + +We will use the `overlap()` function from [{nicheROVER}](https://cran.r-project.org/package=nicheROVER)) to estimate the percentage of similarity among species. We will set overlap to assess based on 95% similarities. + +```{r} +over_stat <- overlap(fish_par, nreps = nsample, + nprob = 1000, + alpha = 0.95) +``` + +We then are going transform this output to a data frame using `extract_overlap()` plotting so we can assess overall similarities among species. + +```{r} +over_stat_df <- extract_overlap(data = over_stat) + +``` + +We then are going to take our newly made data frame and extract out the median percentage of similarities and the 2.5% and 97.5% quantiles. + +```{r, message = FALSE} +over_sum <- over_stat_df %>% + group_by(sample_name_a, sample_name_b) %>% + summarise( + median_niche_overlap = round(median(niche_overlap_perc), digits = 2), + qual_2.5 = round(quantile(niche_overlap_perc, + probs = 0.025, na.rm = TRUE), digits = 2), + qual_97.5 = round(quantile(niche_overlap_perc, + probs = 0.975, na.rm = TRUE), digits = 2) + ) %>% + ungroup() %>% + pivot_longer(cols = -c(sample_name_a, sample_name_b, median_niche_overlap), + names_to = "percentage", + values_to = "niche_overlap_qual") %>% + mutate( + percentage = as.numeric(str_remove(percentage, "qual_")) + ) +``` + + + +We are now going to use `ggplot()`, `geom_violin()`, and `stat_summary()` to represent the posterior distributions and the median of the posterior distributions. Representations of posterior distributions should either use the mode or median, not the mean. See the following publication as a guide to representation of posterior distributions in vizialations and in text. You can extract credible intervals or the equal-tailed intervals using [{bayestestR}](https://easystats.github.io/bayestestR/). + +```{r, warning = FALSE} +#| out-width: 100% +ggplot(data = over_stat_df, aes(x = sample_name_a, + y = niche_overlap_perc, + fill = sample_name_b)) + + geom_violin() + + stat_summary(fun.y = median, geom = "point", + size = 3, + position = position_dodge(width = 0.9)) + + geom_vline(xintercept = seq(1.5, 3.5, 1), + linetype = 2) + + scale_fill_viridis_d(begin = 0.25, end = 0.75, + option = "D", name = "Species", + alpha = 0.35) + + theme_bw() + + theme( + panel.grid = element_blank(), + axis.text = element_text(colour = "black"), + legend.background = element_blank(), + strip.background = element_blank() + ) + + labs(x = paste("Overlap Probability (%)", "\u2013", + "Niche Region Size: 95%"), + y = "p(Percent Overlap | X)") +``` + +### Estimate overall niche size + +We are now going to estimate the overall size of the niche for each posterior sample by using the function `extract_niche_size()` which is a wrapper around `niche.size()` and some data manipulation functions. + +```{r} +niche_size <- extract_niche_size(fish_par) +``` + + +### Plot niche size + +We will now use `geom_violin()`, `geom_point()`, and `geom_errorbar()` to plot the distribution for niche size for each species. + +```{r} +ggplot(data = niche_size, + aes(x = sample_name, y = niche_size)) + + geom_violin( + + width = 0.2) + + stat_summary(fun.y = median, geom = "point", + size = 3, + position = position_dodge(width = 0.9)) + + theme_bw(base_size = 15) + + theme(panel.grid = element_blank(), + axis.text = element_text(colour = "black")) + + labs(x = "Species", + y = "Niche Size") +``` + +Now that we have our niche sizes and similarities determined we can make inferences about the species, trophic similarities, and the ecosystem. From 9ff4155fe76f68830bbf5555e2392265f29e8721 Mon Sep 17 00:00:00 2001 From: YourName Date: Wed, 18 Dec 2024 11:35:04 -0500 Subject: [PATCH 34/35] renamed vigignettes --- vignettes/using-nichetools-with-the-package-SIBER.Rmd | 4 ++-- .../using-nichetools-with-the-package-nicheROVER.Rmd | 10 +++++----- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/vignettes/using-nichetools-with-the-package-SIBER.Rmd b/vignettes/using-nichetools-with-the-package-SIBER.Rmd index 3c61d4b..50a033d 100644 --- a/vignettes/using-nichetools-with-the-package-SIBER.Rmd +++ b/vignettes/using-nichetools-with-the-package-SIBER.Rmd @@ -393,7 +393,7 @@ viridis_colors_s <- viridis(4, begin = 0.25, end = 0.85, Now we can use point interval plots to visually represent posterior distributions. -```{r} +```{r, eval = FALSE} #| out-width: 100% ggplot() + @@ -483,7 +483,7 @@ viridis_colors <- viridis(2, begin = 0.25, end = 0.85, Lastly, we can visualize the distributions of the posterior estimates using `stat_pointinterval()` from {ggdist}. -```{r} +```{r, eval = FALSE} #| out-width: 100% ggplot() + diff --git a/vignettes/using-nichetools-with-the-package-nicheROVER.Rmd b/vignettes/using-nichetools-with-the-package-nicheROVER.Rmd index 9ea9442..a11dec1 100644 --- a/vignettes/using-nichetools-with-the-package-nicheROVER.Rmd +++ b/vignettes/using-nichetools-with-the-package-nicheROVER.Rmd @@ -1,11 +1,11 @@ --- -title: "Using {nichetools} with nicheROVER" +title: "Using {nichetools} and {nicheROVER} with two isotopes" author: Benjamin L. Hlina date: "2024-03-15" output: rmarkdown::html_vignette vignette: > - %\VignetteIndexEntry{Using {nichetools} with nicheROVER} + %\VignetteIndexEntry{Using {nichetools} and {nicheROVER} with two isotopes} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- @@ -19,7 +19,7 @@ knitr::opts_chunk$set( ### Our Objectives -The purpose of this vignette is to use [{nicheROVER}](https://cran.r-project.org/package=nicheROVER) and [{nichetools}](https://benjaminhlina.github.io/nichetools/) to extract and then visualize estimates of trophic niche size and similarities for multiple freshwater fish using [{ggplot2}](https://ggplot2.tidyverse.org/). +The purpose of this vignette is to use [{nicheROVER}](https://cran.r-project.org/package=nicheROVER) and [{nichetools}](https://benjaminhlina.github.io/nichetools/) to extract and then visualize estimates of isotopic niche size and similarities using two isotopes for multiple freshwater fish using [{ggplot2}](https://ggplot2.tidyverse.org/). This vignette can be used for additional purposes including estimating niche size and similarities among different groups of aquatic and/or terrestrial species. Furthermore, niche size and similarities for different behaviours exhibited within a population can be made using behavioural data generated from acoustic telemetry (e.g., differences in habitat occupancy). @@ -44,7 +44,7 @@ I will add that many of the `{dplyr}` and `{tidyr}` functions and processes can } ``` -For the purpose of the vignette we will be using the `fish` data frame that is available within `{nicheROVER}`. We will remove $\delta$34S for simplicity of the vignette. If more than two isotopes or metrics are being used to compare niche sizes and similarities, please use the functions for each pairing. Right now some functions (i.e., `niche_ellipse()`) in `{nichetools}` doesn't have the ability to work with more than two isotopes. This will become a feature at some point but for now. Please be patient and use the functions for each pairing you have. +For the purpose of the vignette we will be using the `fish` data frame that is available within `{nicheROVER}`. We will remove $\delta$34S for simplicity of the vignette. If more than two isotopes or metrics are being used to compare niche sizes and similarities, please remember to use the `isotope_n` argument in the following functions `extract_sigma()` and `niche_ellipse()`. We will first use the function `janitor::clean_names()` to clean up column names. For your purposes you will need to replace fish with your data frame either by loading a csv, rds, or qs, with your data. You can do this multiple ways, I prefer using `readr::read_csv()` but base R's `read.csv()` works perfectly fine. @@ -217,7 +217,7 @@ sigma_plots[[1]] + ### Estimate niche ellipse -We then will use `niche_ellipse()` to easily extract ellipse for each $\Sigma$ estimate (i.e., 1000). If you are to have additional isotopes or metrics, you will need to create mu and sigma objects for each pairing, as currently this function only handles two isotopes. In the future, there likely will be the ability to specify the number of isotopes you have with the default being two. The reason for the lack of functionality is `ellipse::ellipse()` can only work within two-dimensions, not three, so you will have to create multiple `ellipse()` calls for each combination of isotopes or metrics and I haven't had the time to implement this. The function will also tell you how long it took to process as with large sets of isotope data it is nice to know the time it takes for the function to work. The function also has the argument `random` which by default is set to `TRUE`. This argument randomly subsamples and returns 10 ellipse estimates out of the total number of samples taken in this case it is 1,000. The `set_seed` argument allows you to change the `set.seed` value by giving a numerical value to make the results of the function reproducable, default is a random value. I highly suggest using set_seed otherwise it will subsample different values, there is no default value because CRAN will not allow for a default value. You can change the number of subsamples but 10 seems pretty standard. If you'd like all 1,000 ellipses you can set `random` to `FALSE`. +We then will use `niche_ellipse()` to easily extract ellipse for each $\Sigma$ estimate (i.e., 1000). The function will tell you how long it took to process as with large sets of isotope data it is nice to know the time it takes for the function to work. The function also has the argument `random` which by default is set to `TRUE`. This argument randomly subsamples and returns `10` ellipse estimates out of the total number of samples taken in this case it is 1,000. The `set_seed` argument allows you to change the `set.seed` value by giving a numerical value to make the results of the function reproducible, default is a random value. I highly suggest using set_seed otherwise it will subsample different values, there is no default value because CRAN will not allow for a default value. You can change the number of subsamples but 10 seems pretty standard. If you'd like all 1,000 ellipses you can set `random` to `FALSE`. ```{r} ellipse_df <- niche_ellipse(dat_mu = df_mu, dat_sigma = df_sigma, From 5728eaa579fdd2a63277f75517eae8ccec10a422 Mon Sep 17 00:00:00 2001 From: YourName Date: Wed, 18 Dec 2024 11:50:19 -0500 Subject: [PATCH 35/35] replace with create_iso_pairs --- R/niche_ellipse.R | 22 +++------------------- 1 file changed, 3 insertions(+), 19 deletions(-) diff --git a/R/niche_ellipse.R b/R/niche_ellipse.R index 4070c04..bce4029 100644 --- a/R/niche_ellipse.R +++ b/R/niche_ellipse.R @@ -222,25 +222,9 @@ niche_ellipse <- function( sample_numbers <- sample(dat_mu$sample_number, n) # ---- create every combo in table ---- - iso_combo <- tidyr::expand_grid(iso_a = isotope_names, - iso_b = isotope_names) |> - dplyr::filter(iso_a < iso_b) |> - dplyr::mutate( - iso_ab = paste(iso_a, iso_b, sep = "_") - ) |> - dplyr::distinct(iso_ab) |> - tidyr::separate_wider_delim(iso_ab, names = c("iso_a", "iso_b"), - delim = "_") - - # created id column - iso_combo$id <- 1:nrow(iso_combo) - - # pivot longer to create vector that can be used to filter combinations - split_iso <- split(iso_combo, iso_combo$id) |> - purrr::map(~ tidyr::pivot_longer(.x, cols = -id, - values_to = "iso_name") |> - getElement("iso_name") - ) + + split_iso <- create_isotope_pairs() + # ----- mu ---- # select columns we need