R/ddiwas_regression_progress.R

In patrickwu510/ddiwas: What the Package Does (One Line, Title Case)

#' DDIWAS regression with progress bar
#'
#' @param regression_table A tibble
#' @param regression_drugs A tibble
#' @param core_div A scalar 
#' @return A dataframe \code{ddiwas_results}
#' @export

ddiwas_regression_progress <- function(regression_table, regression_drugs, core_div=2) {

  library(tidyverse)
  library(logistf)
  library(parallel)
  library(MASS)
  library(pbmcapply)
  library(varhandle)
  numCores <- (detectCores()/core_div)

  regression_table <- subset(regression_table, select = -c(person_id))
  regression_drugs<- subset(regression_drugs, select = -c(person_id))

  df.cases <- filter(regression_table, groupc == 1)
  df.controls <- filter(regression_table, groupc == 0)
  nCases <- nrow(df.cases)
  nControls <- nrow(df.controls)

  ddiwas_fx <- function(drug) try({
    regression_table["drug"] <- regression_drugs[drug]
    # run model to get statistics
    firthModel <- logistf(groupc~drug+obs_length_n+age_n+is_m+is_w, data=regression_table, family=binomial)
    # calculate 2x2 contingency table numbers
    regression_table_cases <- filter(regression_table, groupc == 1)
    regression_table_controls <- filter(regression_table, groupc == 0)
    nA <- sum(regression_table_cases["drug"])
    nB <- sum(regression_table_controls["drug"])
    nC <- nCases-nA
    nD <- nControls-nB
    # store results
    results <- c(drug,
                 firthModel$coefficients[2],
                 sqrt(diag(vcov(firthModel)))[2],
                 firthModel$prob[2],
                 exp(firthModel$coefficients[2]),
                 nA,nB,nC,nD)
    })
  drugs <- names(regression_drugs)
  ddiwas_results <- pbmclapply(drugs, ddiwas_fx, mc.cores=numCores)
  ddiwas_results1 <- as.data.frame(do.call(rbind, ddiwas_results))
  names(ddiwas_results1) <- c("rxcui_ingr","coef","se","pval","or","nA","nB","nC","nD")
  ddiwas_results1 <- unfactor(as_tibble(ddiwas_results1))

  return(ddiwas_results1)
}