diff --git a/DESCRIPTION b/DESCRIPTION index 6802009..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", @@ -27,7 +27,8 @@ Imports: rlang, SIBER, tibble, - tidyr + tidyr, + tidyselect License: CC0 Encoding: UTF-8 LazyData: true diff --git a/NAMESPACE b/NAMESPACE index bfe1bb9..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) @@ -15,6 +16,7 @@ import(ellipse) import(purrr) import(tibble) import(tidyr) +import(tidyselect) importFrom(lifecycle,deprecated) importFrom(nicheROVER,niche.size) importFrom(rlang,":=") 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) +} 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" } diff --git a/R/extract_layman.R b/R/extract_layman.R index 19b80d0..dc61339 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)) |> @@ -224,7 +224,7 @@ extract_layman <- function(data, return(df_layman) } - if (data_format %in% "wide"){ + if (data_format %in% "wide") { return(df_layman) } } diff --git a/R/extract_mu.R b/R/extract_mu.R index 63c1555..b0ee74d 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 @@ -98,28 +100,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 +174,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 +196,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() @@ -219,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) diff --git a/R/extract_overlap.R b/R/extract_overlap.R index 4e3d952..a623433 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" } @@ -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 ) diff --git a/R/extract_sigma.R b/R/extract_sigma.R index 09b399d..786a539 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) @@ -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 = ".") @@ -220,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() diff --git a/R/niche_ellipse.R b/R/niche_ellipse.R index d6ac53b..bce4029 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 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 #' posterior distributions for \eqn{\mu} and \eqn{\Sigma} to create a sub-sample #' of ellipse. Default is `TRUE`. @@ -48,6 +48,7 @@ #' @importFrom stats runif #' @import tibble #' @import tidyr +#' @import tidyselect #' @export @@ -55,8 +56,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, @@ -66,7 +67,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.") @@ -75,26 +75,30 @@ 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.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") + } } - if(is.null(p_ell)) { + if (is.null(p_ell)) { p_ell <- 0.95 } # Additional parameter validation, if needed @@ -135,135 +139,428 @@ niche_ellipse <- function( set.seed(set_seed) - sample_numbers <- sample(dat_mu$sample_number, n) - # 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" + 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 ---- + + split_iso <- create_isotope_pairs() + + # ----- mu ---- + + # select columns we need + mu_select <- dat_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 <- dat_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" ) + # ---- 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) + ) - # 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") + # convert id of iso_combo to character for joining + iso_combo$id <- as.character(iso_combo$id) - # 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 <- dat_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 <- dat_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() @@ -273,7 +570,14 @@ niche_ellipse <- function( sep = " ")) } - return(all_ellipses) + return(ell_final) + + + # return(all_ellipses) + } + } + + } diff --git a/R/zzz.R b/R/zzz.R index 659d506..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", @@ -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" + )) 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/ 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/" ) 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 +} 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 diff --git a/man/niche_ellipse.Rd b/man/niche_ellipse.Rd index 5ebf041..22698c3 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 analysis. 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). 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") + +}) 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", { 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. 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,