## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----eval = FALSE-------------------------------------------------------------
# # if (!require("remotes", quietly = TRUE))
# #     install.packages("remotes")
# #
# # remotes::install_github("statdivlab/radEmu")

## ----setup, message = FALSE---------------------------------------------------
library(radEmu)
library(dplyr)

## -----------------------------------------------------------------------------
data("wirbel_sample")
# encode control as the baseline level
wirbel_sample$Group <- factor(wirbel_sample$Group, levels = c("CTR","CRC"))
dim(wirbel_sample)
data("wirbel_otu")
dim(wirbel_otu)

## -----------------------------------------------------------------------------
mOTU_names <- colnames(wirbel_otu)
chosen_genera <- c("Eubacterium", "Faecalibacterium", "Fusobacterium", "Porphyromonas")
mOTU_name_df <- data.frame(name = mOTU_names) %>% 
  mutate(base_name = stringr::str_remove(mOTU_names, "unknown ") %>%
                      stringr::str_remove("uncultured ")) %>%
  mutate(genus_name = stringr::word(base_name, 1))
restricted_mOTU_names <- mOTU_name_df %>%
  filter(genus_name %in% chosen_genera) %>%
  pull(name)

## -----------------------------------------------------------------------------
ch_study_obs <- which(wirbel_sample$Country %in% c("CHI"))

## -----------------------------------------------------------------------------
small_Y <- wirbel_otu[ch_study_obs, restricted_mOTU_names]
sum(rowSums(small_Y) == 0) # no samples have a count sum of 0 
sum(colSums(small_Y) == 0) # one category has a count sum of 0

category_to_rm <- which(colSums(small_Y) == 0)
small_Y <- small_Y[, -category_to_rm]
sum(colSums(small_Y) == 0)

## -----------------------------------------------------------------------------
mod <- emuFit(data = wirbel_sample[ch_study_obs, ],
              Y = small_Y,
              formula = ~ Group + Gender, 
              test_kj = data.frame(k = 2, j = 1))

## -----------------------------------------------------------------------------
head(mod$coef)