R/bb_mods_helper.R

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

###############################
##
## Project: tidy.micro
##
## Purpose: Helper functions for bb_mods
##
## Author: Charlie Carpenter
## Email: charles.carpenter@cuanschutz.edu
##
## Date Created: 2019-10-10
##
## ---------------------------
## Notes:
##
##
## ---------------------------
## HEADER ----

bb_coef_est <- function(m,...){

  covs <- data.frame(Coef = names(m@coefficients),
                     Est = m@coefficients,
                     Cov_Type = m@coefficients %>%
                       names %>%
                       cov_type(!!!rlang::quos(...)))

  for(i in 1:nrow(covs)){
    ## detecting interactions
    if( stringr::str_detect(covs$Cov_Type[i], ".int") ){
      ## interaction pieces
      int. <- covs$Coef[i] %>% as.character %>% stringr::str_split(":") %>% unlist

      ## Exact tring matches of the main effects
      ## (there will only be 2 for any 2x2 interaction piece)
      me1 <- covs$Est[stringr::str_detect(covs$Coef, paste0("^", int.[1], "$"))]
      me2 <- covs$Est[stringr::str_detect(covs$Coef, paste0("^", int.[2], "$"))]

      ## interaction detected from if statement
      covs$Est[i] <- covs$Est[i] + me1 + me2
    }
  }

  covs
}


## Estimates profile likelihood CIs using confint for main effects
## and Wald CIs for interactions after adjusting estiamtes using coef_est
bb_ci_est <- function(m, ...){

  # Updating estimates in case there are interactions
  covs <- bb_coef_est(m = m, ...)

  CI <- matrix(rep(0, nrow(covs)*2), nrow = nrow(covs), ncol = 2,
               dimnames = list(rownames(covs), c("2.5 %", "97.5 %")))

  # covariance matrix
  s <- VGAM::vcov.vlm(m)

  # Confidence intervals
  for(r in 1:nrow(covs)){
    ## detecting interactions
    if( stringr::str_detect(covs$Cov_Type[r], ".int") ){

      # interaction pieces
      int. <- c(covs$Coef[r] %>% as.character %>%
                  stringr::str_split(":") %>% unlist,
                covs$Coef[r] %>% as.character)

      # pullin out rows and cols with interaction and main effect covariances
      cl <- colnames(s) %in% int. ; ro <- rownames(s) %in% int. ; s_ <- s[ro,cl]

      # summing vars and covars
      v <- 0
      for(i in 1:nrow(s_)){
        for(j in 1:ncol(s_)){

          if(j >= i){
            p <- ifelse(i == j, s_[i,j], 2*s_[i,j]) # Var or 2*Cov

            v <- v + p
          }
        }
      }

      # Wald interval for interaction
      ci <- covs$Est[r] + c(-1,1)*1.96*sqrt(v)

    } else{

      # Wald interval for main effects
      ci <- m@coefficients[r] + c(-1,1)*1.96*sqrt(diag(s))[r]
    }

    CI[r,] <- ci
  }

  ## Exponentiate and round
  CI %<>% exp %>% round(4)

  ## Make neat for the table
  paste("(", CI[,1], ", ", CI[,2], ")", sep = "")

}

## Function to run multivariable LRT and attach to output df
Anova_SS_bb <- function(conv, bb, SS_type){

  ## Coefs and LRT p-value
  aa <- VGAM::anova.vglm(bb, type = SS_type)
  a <- data.frame(coef = rownames(aa), Anova = aa$`Pr(>Chi)`)

  ## Adding LRT p-values from anova.vglm
  conv$Anova <- NA
  for(i in 1:nrow(a)){
    conv$Anova[grepl(a$coef[i], conv$Coef)] <- a$Anova[i]
  }

  conv
}

## Makes appropriate betabinomial formula from input
formula_fun_bb <- function(...){
  rlang::quos(...) %>%           ## Making "..." a quosure
    rlang::splice() %>%            ## Splicing quosure
    unlist %>%            ## Unlisting
    stringr::str_flatten(collapse = " + ") %>%            ## Combining elements of list with "+"
    stringr::str_replace_all(pattern = "~", replacement = "") %>%     ## Removing "~" (why is it there?)
    paste("cbind(cts, Total - cts) ~ ",., sep = "") %>%            ## pasting with "cts ~ " for final formula
    return()
}