Nothing
## ---- include=FALSE-----------------------------------------------------------
# install.packages("remotes")
# remotes::install_github("adw96/breakaway")
if (requireNamespace('dplyr', quietly = TRUE)) {
library(dplyr)
} else {
message("'dplyr' is not installed!")
}
library(breakaway)
library(magrittr)
library(phyloseq)
## -----------------------------------------------------------------------------
data("soil_phylo")
soil_phylo %>% sample_data %>% head
## -----------------------------------------------------------------------------
subset_soil <- soil_phylo %>%
subset_samples(Amdmt == 1) %>% # only biochar
subset_samples(Day %in% c(0, 2)) # only Days 0 and 82
## -----------------------------------------------------------------------------
richness_soil <- subset_soil %>% breakaway
plot(richness_soil, physeq=subset_soil, color="Day", shape = "ID")
## -----------------------------------------------------------------------------
summary(richness_soil) %>% tibble::as_tibble()
## -----------------------------------------------------------------------------
meta <- subset_soil %>%
sample_data %>%
tibble::as_tibble() %>%
dplyr::mutate("sample_names" = subset_soil %>% sample_names )
## -----------------------------------------------------------------------------
combined_richness <- meta %>%
dplyr::left_join(summary(richness_soil),
by = "sample_names")
# Old way (still works)
bt_day_fixed <- betta(chats = combined_richness$estimate,
ses = combined_richness$error,
X = model.matrix(~Day, data = combined_richness))
# Fancy new way -- thanks to Sarah Teichman for implementing!
bt_day_fixed <- betta(formula = estimate ~ Day,
ses = error, data = combined_richness)
bt_day_fixed$table
## -----------------------------------------------------------------------------
# Old way (still works)
bt_day_fixed_id_random <- betta_random(chats = combined_richness$estimate,
ses = combined_richness$error,
X = model.matrix(~Day, data = combined_richness),
groups=combined_richness$ID)
# Fancy new way
bt_day_fixed_id_random <-
betta_random(formula = estimate ~ Day | ID,
ses = error, data = combined_richness)
bt_day_fixed_id_random$table
## -----------------------------------------------------------------------------
betta_lincom(fitted_betta = bt_day_fixed_id_random,
linear_com = c(1,1),
signif_cutoff = 0.05)
## -----------------------------------------------------------------------------
subset_soil_days_1_2 <- soil_phylo %>%
subset_samples(Amdmt == 1) %>% # only biochar
subset_samples(Day %in% c(0, 1, 2)) # Days 0, 12, and 82
## -----------------------------------------------------------------------------
meta_days_1_2 <- subset_soil_days_1_2 %>%
sample_data %>%
tibble::as_tibble() %>%
dplyr::mutate("sample_names" = subset_soil_days_1_2 %>% sample_names )
## -----------------------------------------------------------------------------
richness_days_1_2 <- subset_soil_days_1_2 %>%
breakaway
combined_richness_days_1_2 <- meta_days_1_2 %>%
dplyr::left_join(summary(richness_days_1_2),
by = "sample_names")
combined_richness_days_1_2
## -----------------------------------------------------------------------------
bt_day_1_2_fixed_id_random <- betta_random(formula = estimate ~ Day | ID,
ses = error, data = combined_richness_days_1_2)
bt_day_1_2_fixed_id_random$table
## -----------------------------------------------------------------------------
set.seed(345)
submodel_test <- test_submodel(bt_day_1_2_fixed_id_random,
submodel_formula = estimate~1,
method = "bootstrap",
nboot = 100)
submodel_test$pval
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.