R/micro_alpha_reg.R

In CharlieCarpenter/tidy.micro: A Pipeline For Microbiome Analysis And Visualization

#' @title Linear regression on alpha diversities within a micro_set
#' @name micro_alpha_reg
#' @description A simple wrapper to run standard linear regression though the \code{\link[stats]{lm}} function. Will only use alpha diversities distinct libraries (Lib) from the specified table as to not inflate the sample size
#' @param alpha_set A tidy_micro data set with alpha diversities calculated by \code{alpha_div}
#' @param table OTU table of interest
#' @param ... Covariates of interest. Can include interaction terms such as \code{Group*Age}
#' @return A data frame containing the model estimates for each alpha diversity
#' @note Be aware of your minimal sequencing depth as this will be the size of all bootstrapped resamples (rarefied).
#' @examples
#' data(phy); data(cla); data(ord); data(fam); data(met)
#'
#' otu_tabs <- list(Phylum = phy, Class = cla, Order = ord, Family = fam)
#'
#' set <- tidy_micro(otu_tabs = otu_tabs, meta = met) \%>\%
#' filter(day == 7) \%>\% ## Only including first week
#' mutate(bpd1 = factor(bpd1))
#'
#' set_fam_alpha <- set \%>\% alpha_div(table = "Family", min_depth = 5000, min_goods = 90)
#' set_fam_alpha \%>\% micro_alpha_reg(table = "Family", bpd1)
#' @export
micro_alpha_reg <- function(alpha_set, table, ..., RE = NULL){

  if(table %nin% unique(alpha_set$Table)) stop("Specified table is not in supplied micro_set")

  alpha_set %<>%
    dplyr::filter(Table == table) %>%
    dplyr::distinct(Lib, .keep_all = TRUE)

  alphas <- c("Goods", "Sobs", "Chao1", "ShannonE", "ShannonH", "SimpsonD")
  alpha_tab <- NULL
  error_tab <- NULL

  if(is.null(RE)){

    alpha_set %<>% dplyr::distinct(Lib, .keep_all = TRUE)

    for(i in 1:length(alphas)){

      f <- suppressWarnings(
        as.formula(paste0(alphas[i], "~", adonis_formula(...)))
      )

      mod <- stats::lm(f, data = alpha_set)
      CI <- stats::confint(mod) %>% round(4)

      alpha_tab <- rbind(alpha_tab,
                         cbind(Alpha_Div = alphas[i],
                               mod %>% broom::tidy(),
                               CI_95 = paste0("(", CI[,1], ", ", CI[,2], ")"))
      )
    }

    alpha_tab %<>% dplyr::rename(Coef = term, Beta = estimate, t.stat = statistic)

  } else{
    if(!is.character(RE)) stop("RE must be written as a character string")

    alpha_set %<>% filter(Taxa == unique(Taxa)[1])


    for(i in 1:length(alphas)){


      f <- suppressWarnings(
        as.formula(
          paste(
            paste0(alphas[i], "~", adonis_formula(...)),
            RE, sep = "+"
            )
          )
        )

      mod <- lme4::lmer(f, data = alpha_set)

      CI <- stats::confint(mod) %>% round(4) %>%
        as.data.frame() %>% tibble::rownames_to_column()

      fe <- mod %>% broom::tidy() %>%
        dplyr::filter(group == "fixed") %>%
        dplyr::select(-group)

      fe_CI <- CI %>% filter(rowname %in% fe$term) %>% dplyr::select(-rowname)
      er_CI <- CI %>% filter(rowname %nin% fe$term) %>% dplyr::select(-rowname)

      alpha_tab <- rbind(alpha_tab,
                         cbind(Alpha_Div = alphas[i],
                               fe,
                               CI_95 = paste0("(", fe_CI[,1], ", ", fe_CI[,2], ")"))
                         )

      error_tab <- rbind(error_tab,
                         cbind(Alpha_Div = alphas[i],
                               mod %>% broom::tidy() %>%
                                 dplyr::filter(group != "fixed") %>%
                                 dplyr::select(term, estimate, group),
                               CI_95 = paste0("(", er_CI[,1], ", ", er_CI[,2], ")")
                               )
                         )
    }

    alpha_tab %<>% dplyr::rename(Coef = term, Beta = estimate, t.stat = statistic)
    error_tab %<>% dplyr::rename(Sigma = estimate)

    alpha_tab <- list(FE = alpha_tab, RE = error_tab)
  }

  alpha_tab
}