Problem with simplifyBeta()

[msmcfarlin]

Hello,

I am trying to use simplifyBeta() and running into the error...

Error: Can't subset columns that don't exist. x Columns 'A017', 'A032', 'A034', 'A038', 'A045', etc. don't exist.

This error doesn't seem to be related to the sample names like in issue #73 as my sample names do not have any special characters. The error seems to be related to vctrs::stop_subscript(). Do yo you have any thoughts of what the issue is and how to fix this?

Here are the outputs from rlang::last_error() and rlang::last_trace() in case they are helpful. Let me know if any other information would be helpful!

Thank you!

rlang::last_error()

  1. `%>%`(...)
 19. tidyr:::gather.data.frame(., "Sample2", "var", names(vars))
 21. tidyselect::vars_select(tbl_vars(data), !!!quos)
 22. tidyselect:::eval_select_impl(...)
 30. tidyselect:::vars_select_eval(...)
 31. tidyselect:::walk_data_tree(expr, data_mask, context_mask)
 32. tidyselect:::eval_c(expr, data_mask, context_mask)
 33. tidyselect:::reduce_sels(node, data_mask, context_mask, init = init)
 34. tidyselect:::walk_data_tree(new, data_mask, context_mask)
 35. tidyselect:::as_indices_sel_impl(...)
 36. tidyselect:::as_indices_impl(x, vars, strict = strict)
 37. tidyselect:::chr_as_locations(x, vars)
 38. vctrs::vec_as_location(x, n = length(vars), names = vars)
 40. vctrs:::stop_subscript_oob(...)
 41. vctrs:::stop_subscript(...)

rlang::last_trace()

1. +-`%>%`(...)
  2. +-dplyr::pull(., beta_est)
  3. +-DivNet::simplifyBeta(...)
  4. | \-`%>%`(...)
  5. +-dplyr::mutate(...)
  6. +-dplyr::distinct(., beta_est, .keep_all = TRUE)
  7. +-base::unique(.)
  8. +-dplyr::filter(., beta_est > 1e-16)
  9. +-dplyr::select(., Covar1, Covar2, beta_est, beta_var)
 10. +-dplyr::mutate(., Covar1 = vars[Sample1], Covar2 = vars[Sample2])
 11. +-tibble::add_column(...)
 12. | \-tibble::tibble(..., .name_repair = .name_repair)
 13. |   \-tibble:::tibble_quos(xs, .rows, .name_repair)
 14. |     \-rlang::eval_tidy(xs[[j]], mask)
 15. +-`%>%`(...)
 16. +-base::unlist(.)
 17. +-dplyr::select(., "var")
 18. +-tidyr::gather(., "Sample2", "var", names(vars))
 19. \-tidyr:::gather.data.frame(., "Sample2", "var", names(vars))
 20.   +-base::unname(tidyselect::vars_select(tbl_vars(data), !!!quos))
 21.   \-tidyselect::vars_select(tbl_vars(data), !!!quos)
 22.     \-tidyselect:::eval_select_impl(...)
 23.       +-tidyselect:::with_subscript_errors(...)
 24.       | +-base::tryCatch(...)
 25.       | | \-base:::tryCatchList(expr, classes, parentenv, handlers)
 26.       | |   \-base:::tryCatchOne(expr, names, parentenv, handlers[[1L]])
 27.       | |     \-base:::doTryCatch(return(expr), name, parentenv, handler)
 28.       | \-tidyselect:::instrument_base_errors(expr)
 29.       |   \-base::withCallingHandlers(...)
 30.       \-tidyselect:::vars_select_eval(...)
 31.         \-tidyselect:::walk_data_tree(expr, data_mask, context_mask)
 32.           \-tidyselect:::eval_c(expr, data_mask, context_mask)
 33.             \-tidyselect:::reduce_sels(node, data_mask, context_mask, init = init)
 34.               \-tidyselect:::walk_data_tree(new, data_mask, context_mask)
 35.                 \-tidyselect:::as_indices_sel_impl(...)
 36.                   \-tidyselect:::as_indices_impl(x, vars, strict = strict)
 37.                     \-tidyselect:::chr_as_locations(x, vars)
 38.                       \-vctrs::vec_as_location(x, n = length(vars), names = vars)
 39.                         \-(function () ...
 40.                           \-vctrs:::stop_subscript_oob(...)
 41.                             \-vctrs:::stop_subscript(...)

[ailurophilia]

Hi @msmcfarlin,

Thanks for letting us know about this. Are you able to put together a minimal reproducible example you can share?

Thanks!

Best, David

[msmcfarlin]

Hi @ailurophilia,

Yes, here is a reproducible sample. Let me know if you have any issues with it!

#divnet phyloseq data
py <- B2_phyloseq

#Tax glom to phylum for time and update sample names to S1, S2...
py.phy <- py %>% tax_glom(taxrank="phylum")
sample_names(py.phy) <- paste("S",c(1:31),sep = "")

#Set diagonal design matrix for beta analysis
X <- diag(nrow(sample_data(py.phy)))

colnames(X) <- rownames(X) <-  as.character(1:ncol(X))

#Run divnet with diagonal
divnet_phylum <-  divnet(W = py.phy,
                         X = X,
                         variance = "none")

#Change metadata variable V1 to character
py.phy@sam_data$V1 <- as.character(py.phy@sam_data$V1)

#Pull estimates
estimates <- simplifyBeta(divnet_phylum,
                          py.phy,
                          "euclidean", "V1") %>% pull(beta_est)

[ailurophilia]

Hi @msmcfarlin,

Sorry for the delay getting back to you. It looks like the issue was a result of simplifyBeta expecting non-NULL variance estimates in the DivNet object passed to it. I've patched simplifyBeta so it no longer requires this, so your code should work now with the latest version of DivNet.

Best, David

[msmcfarlin]

Hi @ailurophilia,

Thank you for the update!

I updated DivNet and simplifyBeta is working when "euclidean" or "bray-curtis" are the measurements that are specified. Though I run into the same issue when "aitchison" is used.

[msmcfarlin]

Hi @ailurophilia,

Just to check, you said the issue was...

a result of simplifyBeta expecting non-NULL variance estimates in the DivNet object passed to it.

I had a NULL variance estimate because I was following the recommendation of the beta diversity vignette. Is it still recommended to set variance equal to "none" in the divnet function then bootstrapping in a subsequent step for beta diversity analysis?

Thank you for the assistance!

[ailurophilia]

Yes – this is still recommended. The behavior of simplifyBeta here was a bug.

David