In enaY15/TabulationAutomation: Auomated analysis and tabulation of results

options(tinytex.verbose = TRUE)

source('combinedLaTeXtableR.R', encoding = 'UTF-8')
source('combinedLaTeXtableR.R', encoding = 'UTF-8')
source('combinedLaTeXtableR.R', encoding = 'UTF-8')

require(pscl, quietly = TRUE); require(tidyr, quietly = TRUE); require(MASS, quietly = TRUE); require(brms, quietly = TRUE); require(magrittr, quietly = TRUE); require(knitr, quietly = TRUE); require(dplyr, quietly = TRUE); require(htmlTable, quietly = TRUE); require(tibble, quietly = TRUE); require(psych, quietly = TRUE); require(stats, quietly = TRUE); require(MplusAutomation, quietly = TRUE); require(papaja, quietly = TRUE); require(Hmisc, quietly = TRUE); require(pscl, quietly = TRUE)

homepath <- paste0('C:/Users/', Sys.getenv('USERNAME'))

#### Function to read in files ####
dataRead <- function(homepath, file) eval(parse(text = paste0('read.csv("', homepath, file, ')')))

n3.a1 <- dataRead(homepath, file = '/Box Sync/Research/0 Data 0/NESARCIII/nesarcA1final.csv", na.strings = (c("9999", "NA")))')
n3.a1$CASEID <- as.character(n3.a1$CASEID)
n3 <- dataRead(homepath, file = '/Box Sync/Research/0 Data 0/NESARCIII/n3csv.csv", na.strings = (c("9999", "NA")))')

sf12 <- read.csv('sf12_n3.csv', na.strings = '9')
  colnames(sf12)[1] <- 'CASEID'

sf12$calm <- sf12$N1Q33A 
sf12$depressed <- sf12$N1Q33C 
sf12$accomplishedless <- sf12$N1Q32A
sf12$carefulless <- sf12$N1Q32B
sf12$GeneralHealth <- sf12$N1Q25
sf12$gh <- ifelse(sf12$N1Q25 == levels(sf12$N1Q25)[1], 1, ifelse(sf12$N1Q25 == levels(sf12$N1Q25)[2], 2, ifelse(sf12$N1Q25 == levels(sf12$N1Q25)[3], 3, ifelse(sf12$N1Q25 == levels(sf12$N1Q25)[4], 4, 5))))

diagnoses <- read.csv('diagnoses.csv'); diagnoses <- diagnoses[, -1]
#impairment <- read.csv('impairment.csv'); impairment <- impairment[, -1]
outcomes <- read.csv('outcomes.csv'); outcomes <- outcomes[, -1]
outcomes <- outcomes[, -which(names(outcomes) %in% names(diagnoses))]
outcomes <- data.frame(n3.a1$CASEID, outcomes); colnames(outcomes)[1] <- 'CASEID'

outcomes <- cbind(sf12[sf12$CASEID %in% outcomes$CASEID, -1])

#### Drinking Data Cleaning ####
nesarc3 <- dataRead(homepath, file = '/Box Sync/Research/0 Data 0/NESARCIII.csv")')
nesarc3.a1 <- nesarc3[nesarc3$CASEID %in% n3.a1$CASEID, ]
rm(nesarc3)

## MAXDRINKS
MAXDRINKS <- ifelse(is.na(nesarc3.a1$S2AQ4C), 0, ifelse(nesarc3.a1$S2AQ4C==99, NA, nesarc3.a1$S2AQ4C))

## FREQENCY DRINKING
DRINKFREQ <- ifelse(is.na(nesarc3.a1$S2AQ4A), 0, ifelse(nesarc3.a1$S2AQ4A==99, NA, 11-nesarc3.a1$S2AQ4A))

## USUAL AMOUNT DRINKING
USUALAMT <- ifelse(is.na(nesarc3.a1$S2AQ4B), 0, ifelse(nesarc3.a1$S2AQ4B==99, NA, nesarc3.a1$S2AQ4B))

## FREQENCY MAX
FREQMAX <- ifelse(is.na(nesarc3.a1$S2AQ4E), 0, ifelse(nesarc3.a1$S2AQ4E==99, NA, 11-nesarc3.a1$S2AQ4E))
#HOW OFTEN DRANK LARGEST NUMBER OF DRINKS IN A SINGLE DAY IN LAST 12 MONTHS

## BINGING WITH DifFERENT CUTOFFS FOR MEN AND WOMEN
BINGE <- ifelse((nesarc3.a1$SEX==1 & is.na(nesarc3.a1$S2AQ4H)), 0, ifelse((nesarc3.a1$SEX==1 & nesarc3.a1$S2AQ4H==99), NA, ifelse(nesarc3.a1$SEX==1, 11-nesarc3.a1$S2AQ4H, ifelse((nesarc3.a1$SEX==2 & is.na(nesarc3.a1$S2AQ4F)), 0, ifelse((nesarc3.a1$SEX == 2 & nesarc3.a1$S2AQ4F  == 99), NA, 11-nesarc3.a1$S2AQ4F)))))

## HOW OFTEN INTOXICATED
INTOX <- ifelse(is.na(nesarc3.a1$S2AQ9), 0, ifelse(nesarc3.a1$S2AQ9==99, NA, 11-nesarc3.a1$S2AQ9))

## FREQUENCY OF EXCEEDING DAILY DRINKING LIMITS
#if (is.na(nesarc3.a1$FEXMAX)) FREQEX <- 0 else if (nesarc3.a1$FEXMAX==99) FREQEX <- NA else FREQEX <- FEXMAX

#### CREATE DRINKING + OUTCOME DF ####
drinking <- data.frame(MAXDRINKS, DRINKFREQ, USUALAMT, FREQMAX, BINGE, INTOX
                       #, FREQEX
                       )
outcomes <- cbind(outcomes, drinking)

## Factors and Models
factorModelNameVector <- c('dsm5', 'dysphoria', 'dysphoricarousal', 'anhedonia', 'externalizing', 'hybrid')
modelFactorNamesList <- list(c('INTRUSIONS', 'AVOIDANCE', 'COGMOOD', 'HYPERAROUSAL'), 
             c('INTRUSIONS', 'AVOIDANCE', 'DYSPHORIA', 'HYPERAROUSAL'), 
             c('INTRUSIONS', 'AVOIDANCE', 'DYSPHORIA', 'DYSPHORIC_AROUSAL', 'ANXIOUS_AROUSAL'), 
             c('INTRUSIONS', 'AVOIDANCE', 'NEGAFF', 'ANHEDONIA', 'DYSPHORIC_AROUSAL', 'ANXIOUS_AROUSAL'), 
             c('INTRUSIONS', 'AVOIDANCE', 'NUMBING', 'EXTERNALIZING', 'DYSPHORIC_AROUSAL', 'ANXIOUS_AROUSAL'), 
             c('INTRUSIONS', 'AVOIDANCE', 'NEGAFF', 'ANHEDONIA', 'EXTERNALIZING', 'DYSPHORIC_AROUSAL', 'ANXIOUS_AROUSAL'))
names(modelFactorNamesList) <- factorModelNameVector

## Outcomes
sf12_outcomes <- c('gh', 'calm', 'depressed', 'accomplishedless', 'carefulless')
#dx_outcomes <- c('pyptsd', 'lptsd', 'pymdepind', 'lmdepind', 'pymddisorder' , 'lmddisorder', 'ldysind', 'bpddx1', 'pygadind', 'lgadind', 'pyspeind', 'lspeind', 'pypanicind', 'lpanicind', 'lifeaud5', 'pyaud5', 'lnicdep5', 'lmaud5')
dx_outcomes <- c('pyptsd', 'lptsd', 'pymddisorder' , 'lmddisorder', 'bpddx1', 'pygadind', 'lgadind', 'pyspeind', 'lspeind', 'pypanicind', 'lpanicind', 'pyaud5', 'lifeaud5', 'lmaud5', 'lnicdep5')
dvs <- outcomes
#drinking_outcomes <- c('numwhendrink', 'agefirstdrink', 'freqyearlydrunk')
drinking_outcomes <- c('MAXDRINKS', 'DRINKFREQ', 'USUALAMT', 
                       #'FREQMAX', 
                       'BINGE', 'INTOX'
                       #, 'FREQEX'
                       )
#outcomeNameVector <- c('SF-12 General Health', 'SF-12 Calm or Peaceful', 'SF-12 Down or Depressed', 'SF-12 Less Accomplished', 'SF-12 Less Careful', 'Past Year PTSD', 'Lifetime PTSD', 'Past Year MDD', 'Life MDD', 'BPD', 'Past Year GAD', 'Life GAD', 'pyspeind', 'lspeind', 'Past Year Panic', 'Life Panic', 'Past Year AUD', 'Life AUD', 'lmaud5', 'Life Nicotine Dependence', 'Number Typical Drinks', 'Age 1st Drink', 'Intoxication Frequency')
outcomeNameVector <- c('SF-12 General Health', 'SF-12 Calm or Peaceful', 'SF-12 Down or Depressed', 'SF-12 Less Accomplished', 'SF-12 Less Careful', 'Past Year PTSD', 'Lifetime PTSD', 'Past Year MDD', 'Life MDD', 'BPD', 'Past Year GAD', 'Life GAD', 'pyspeind', 'lspeind', 'Past Year Panic', 'Life Panic', 'Past Year AUD', 'Life AUD', 'lmaud5', 'Life Nicotine Dependence', 'Maximum Drinks', 'Drinking Frequency', 'Usual Drink Number', 'Frequency Drink Maximum', 'Binge Frequency', 'Intoxication Frequency', 'Frequency Exceed Daily Limits')
#names(outcomeNameVector) <- c('gh', 'calm', 'depressed', 'accomplishedless', 'carefulless', 'lptsd', 'pyptsd', 'laud5', 'pyaud5', 'numwhendrink', 'agefirstdrink', 'freqyearlydrunk.1')
names(outcomeNameVector) <- c(sf12_outcomes, dx_outcomes, drinking_outcomes)

lca_cprob <- read.csv('lca_cprob_nesarc.csv')
lca <- data.frame(lca_cprob, outcomes, diagnoses) 

#### BCFA FScores ####
#fscoresMakeR <- function() {
  for (m in 1:length(factorModelNameVector)) {
    model <- factorModelNameVector[m]
    factors <- modelFactorNamesList[[model]] # factors in structural model    
    filePath <- paste0(getwd(), '/FScores_N3_BCFA_', length(factors), '_', model, '.csv')
    fscores <- read.delim(file = filePath, header = FALSE, sep = '', na.strings = '*')
    #fscores <- eval(parse(text = paste0('fscores_cfa_', model))) # Mplus fscores output
    fNames <- seq(from = -1*(dim(fscores)[2] - 4), by = 5, length.out = length(factors))
        # negative numbers so can count up starting with 5th from end of df (last factor score) is reference point, 
    fNames <- fNames[(order(fNames*-1))] # then reorder by *-1
    fNames <- gsub(fNames, pattern = '-', replacement = 'V') # capitalize on negative signs left over to easily replace with "V" to make variable names 
    data <- data.frame(fscores[, fNames])
    colnames(data) <- factors
    eval(envir = globalenv(), parse(text = paste0('data_bcfa_', model, ' <- data')))
  }
#}

#### Bifactor FScores ####
#fscoresMakeR_bf <- function() {
  for (m in 1:length(factorModelNameVector)) {
    #options(browser())
    model <- factorModelNameVector[m]
    factors <- modelFactorNamesList[[model]] # factors in structural model    
    filePath <- paste0(getwd(), '/FScores_N3_BF_', length(factors), '_', model, '.csv')
    factors_bf <- c('GENERAL', modelFactorNamesList[[model]]) # if bifactor
    fscores_bf <- read.delim(file = filePath, header = FALSE, sep = '', na.strings = '*')
    #fscores_bf <- eval(parse(text = paste0('fscores_cfa_', model))) # Mplus fscores output
    fNames_bf <- seq(from = -1*(dim(fscores_bf)[2] - 4), by = 5, length.out = length(factors_bf))
        # negative numbers so can count up starting with 5th from end of df (last factor score) is reference point, 
    fNames_bf <- fNames_bf[(order(fNames_bf*-1))] # then reorder by *-1
    fNames_bf <- gsub(fNames_bf, pattern = '-', replacement = 'V') # capitalize on negative signs left over to easily replace with "V" to make variable names 
    data <- data.frame(fscores_bf[, fNames_bf])
    colnames(data) <- factors_bf
    eval(envir = globalenv(), parse(text = paste0('data_bf_', model, ' <- data')))
  }
#}

# fscoresMakeR()
# fscoresMakeR_bf()

  for (m in 1:length(factorModelNameVector)) {
    eval(envir = globalenv(), parse(text = paste0('data_bcfa_', factorModelNameVector[m], ' <- data.frame(dvs, data_bcfa_', factorModelNameVector[m], ')')))
  }
rm(fscores)

  for (m in 1:length(factorModelNameVector)) {
    eval(envir = globalenv(), parse(text = paste0('data_bf_', factorModelNameVector[m], ' <- data.frame(dvs, data_bf_', factorModelNameVector[m], ')')))
  }
rm(fscores_bf) 

#### Frequentist CFA FScores ####
filepathcfa <- gsub(x = paste0(getwd(), '/FScores_N3_'), pattern = '/', replacement = '//')
for (m in 1:length(factorModelNameVector)) {
  model <- factorModelNameVector[m]
  factors <- modelFactorNamesList[[model]] # factors in structural model    
  factors_se <- rep(paste0(modelFactorNamesList[[model]], '_se')) 
  for(f in 1:length(factors)) c(factors[f], factors_se[f])
  filePath <- paste0(model, '.csv')
  eval(parse(text = paste0('fscores_', model, ' <- read.delim(file = \"', filepathcfa, filePath, '\", header = FALSE, sep = \'\', na.strings = \'*\')')))
  fscores <- eval(parse(text = paste0('fscores_', model))) # Mplus fscores output
  fscores <- fscores[, -c(1:20)] # take out symptom columns
  fscores <- fscores[, c(seq(from = 1, to = dim(fscores)[2], by = 2))] # drop factor SE cols; leave just factors scores
  colnames(fscores) <- factors
  eval(envir = globalenv(), parse(text = paste0('data_cfa_', factorModelNameVector[m], ' <- data.frame(dvs , fscores)')))
}

#### Write FScores to csv ####
  for (m in 1:length(factorModelNameVector)) {
    eval(parse(text = paste0('write.csv(data_bcfa_', factorModelNameVector[m], ', \'data_bcfa_', factorModelNameVector[m], '.csv\')')))
  }
  for (m in 1:length(factorModelNameVector)) {
    eval(parse(text = paste0('write.csv(data_bf_', factorModelNameVector[m], ', \'data_bf_', factorModelNameVector[m], '.csv\')')))
  }
  for (m in 1:length(factorModelNameVector)) {
    dims <- (20 + 2*length(modelFactorNamesList[[m]]))
    cols <- c(1:dims)
    eval(parse(text = paste0('write.csv(data_cfa_', factorModelNameVector[m], ', \'data_cfa_', factorModelNameVector[m], '.csv\', col.names = cols)')))
  }

  for (m in 1:length(factorModelNameVector)) {
    eval(parse(text = paste0('data_cfa_', factorModelNameVector[m], ' <- read.csv(\'data_cfa_', factorModelNameVector[m], '.csv\')')))
  }

  for (m in 1:length(factorModelNameVector)) {
    eval(parse(text = paste0('data_bcfa_', factorModelNameVector[m], ' <- read.csv(\'data_bcfa_', factorModelNameVector[m], '.csv\')')))
  }

  for (m in 1:length(factorModelNameVector)) {
    eval(parse(text = paste0('data_bf_', factorModelNameVector[m], ' <- read.csv(\'data_bf_', factorModelNameVector[m], '.csv\')')))
  }

  for (y in 1:length(drinking_outcomes)) {
    dv <- drinking_outcomes[y]
    for (m in 1:length(factorModelNameVector)) {
      eval(parse(text = paste0('data <- data_cfa_', factorModelNameVector[m])))
      dv <- paste0('data_cfa_', factorModelNameVector[m], '$', dv)
      eval(envir = globalenv(), parse(text = paste0(drinking_outcomes[y], '_cfa_', factorModelNameVector[m], '_results <- printResults(analysis = \'cfa\', outcomeVar = drinking_outcomes[', y, '], data = data, modelNum = ', m, ')')))
    }
  }

## CFA
cat('Frequentist Confirmatory Factor Analysis SF-12 Summary Table')
  for (y in 1:length(sf12_outcomes)) {
    for (m in 1:length(factorModelNameVector)) {
      eval(parse(text = paste0('data <- data_cfa_', factorModelNameVector[m])))
      eval(envir = globalenv(), parse(text = paste0(sf12_outcomes[y], '_cfa_', factorModelNameVector[m], ' <- summaryTableFx(family = \'gaussian\', analysis = \'cfa\', outcomeVar = sf12_outcomes[', y, '], data = data, modelNum = ', m, ')')))
    }
  }

## BCFA
cat('Bayesian Confirmatory Factor Analysis Short Form Health Survey (SF-12) Summary Tables')
  for (y in 1:length(sf12_outcomes)) {
    for (m in 1:length(factorModelNameVector)) {
      eval(parse(text = paste0('data <- data_bcfa_', factorModelNameVector[m])))
      eval(envir = globalenv(), parse(text = paste0(sf12_outcomes[y], '_bcfa_', factorModelNameVector[m], ' <- summaryTableFx(family = \'gaussian\', analysis = \'bcfa\', outcomeVar = sf12_outcomes[', y, '], data = data, modelNum = ', m, ')')))
    }
  }

cat('Bifactor Analysis Short Form Health Survey (SF-12) Summary Tables')
## Bi-Factor
  for (y in 1:length(sf12_outcomes)) {
    for (m in 1:length(factorModelNameVector)) {
      eval(parse(text = paste0('data <- data_bf_', factorModelNameVector[m])))
      eval(envir = globalenv(), parse(text = paste0(sf12_outcomes[y], '_bf_', factorModelNameVector[m], ' <- summaryTableFx(family = \'gaussian\', analysis = \'bf\', outcomeVar = sf12_outcomes[', y, '], data = data, modelNum = ', m, ')')))
    }
  }

cat('Frequentist Confirmatory Factor Analysis and Psychiatric Diagnoses Summary Tables')
  for (y in 1:length(dx_outcomes)) {
    for (m in 1:length(factorModelNameVector)) {
      eval(parse(text = paste0('data <- data_cfa_', factorModelNameVector[m])))
      data <- data.frame(data, diagnoses)
      eval(envir = globalenv(), parse(text = paste0(dx_outcomes[y], '_cfa_', factorModelNameVector[m], ' <- summaryTableFx(family = \'binomial\', analysis = \'cfa\', outcomeVar = dx_outcomes[', y, '], data = data, modelNum = ', m, ')')))
    }
  }

cat('Bayesian Confirmatory Factor Analysis and Psychiatric Diagnoses Summary Tables')
  for (y in 1:length(dx_outcomes)) {
    for (m in 1:length(factorModelNameVector)) {
      eval(parse(text = paste0('data <- data_bcfa_', factorModelNameVector[m])))
      data <- data.frame(data, diagnoses)
      eval(envir = globalenv(), parse(text = paste0(dx_outcomes[y], '_bcfa_', factorModelNameVector[m], ' <- summaryTableFx(family = \'binomial\', analysis = \'bcfa\', outcomeVar = dx_outcomes[', y, '], data = data, modelNum = ', m, ')')))
    }
  }

cat('Biactor Analysis and Psychiatric Diagnoses Summary Tables')
  for (y in 1:length(dx_outcomes)) {
    for (m in 1:length(factorModelNameVector)) {
      eval(parse(text = paste0('data <- data_bf_', factorModelNameVector[m])))
      data <- data.frame(data, diagnoses)
      eval(envir = globalenv(), parse(text = paste0(dx_outcomes[y], '_bf_', factorModelNameVector[m], ' <- summaryTableFx(family = \'binomial\', analysis = \'bf\', outcomeVar = dx_outcomes[', y, '], data = data, modelNum = ', m, ')')))
    }
  }

cat('Frequentist Confirmatory Factor Analysis and Alcohol Validators Summary Tables')
  for (y in 1:length(drinking_outcomes)) {
    dv <- drinking_outcomes[y]
    for (m in 1:length(factorModelNameVector)) {
      eval(parse(text = paste0('data <- data_cfa_', factorModelNameVector[m])))
      dv <- paste0('data_cfa_', factorModelNameVector[m], '$', dv)
      eval(envir = globalenv(), parse(text = paste0(drinking_outcomes[y], '_cfa_', factorModelNameVector[m], ' <- summaryTableFx(family = \'gaussian\', analysis = \'cfa\', outcomeVar = drinking_outcomes[', y, '], data = data, modelNum = ', m, ')')))
    }
  }

cat('Bayesian Confirmatory Factor Analysis and Alcohol Validators Summary Tables')
  for (y in 1:length(drinking_outcomes)) {
    dv <- drinking_outcomes[y]
    for (m in 1:length(factorModelNameVector)) {
      eval(parse(text = paste0('data <- data_bcfa_', factorModelNameVector[m])))
      dv <- paste0('data_bcfa_', factorModelNameVector[m], '$', dv)
      eval(envir = globalenv(), parse(text = paste0(drinking_outcomes[y], '_bcfa_', factorModelNameVector[m], ' <- summaryTableFx(family = \'gaussian\', analysis = \'bcfa\', outcomeVar = drinking_outcomes[', y, '], data = data, modelNum = ', m, ')')))
    }
  }

cat('Bifactor Analysis and Alcohol Validators Summary Tables')
for (y in 1:length(drinking_outcomes)) {
    for (m in 1:length(factorModelNameVector)) {
      eval(parse(text = paste0('data <- data_bf_', factorModelNameVector[m])))
      eval(envir = globalenv(), parse(text = paste0(drinking_outcomes[y], '_bf_', factorModelNameVector[m], ' <- summaryTableFx(family = \'gaussian\', analysis = \'bf\', outcomeVar = drinking_outcomes[', y, '], data = data, modelNum = ', m, ')')))
    }
}

#cat('Frequentist Confirmatory Factor Analysis and Alcohol Validators Correlation Tables')
for (m in 1:length(factorModelNameVector)) {
  eval(parse(text = paste0('data <- data_cfa_', factorModelNameVector[m])))
  eval(envir = globalenv(), parse(text = paste0('drinking_corr_cfa_', factorModelNameVector[m], ' <- corrTableFx(analysis = \'cfa\', outcomes = drinking_outcomes, data = data, modelNum = ', m, ')')))
}

#cat('Frequentist Confirmatory Factor Analysis and Alcohol Validators Correlation Tables')
for (m in 1:length(factorModelNameVector)) {
  eval(parse(text = paste0('data <- data_bcfa_', factorModelNameVector[m])))
  eval(envir = globalenv(), parse(text = paste0('drinking_corr_bcfa_', factorModelNameVector[m], ' <- corrTableFx(analysis = \'bcfa\', outcomes = drinking_outcomes, data = data, modelNum = ', m, ')')))
}

#cat('Frequentist Confirmatory Factor Analysis and Alcohol Validators Correlation Tables')
for (m in 1:length(factorModelNameVector)) {
  eval(parse(text = paste0('data <- data_bf_', factorModelNameVector[m])))
  eval(envir = globalenv(), parse(text = paste0('drinking_corr_bf_', factorModelNameVector[m], ' <- corrTableFx(analysis = \'bf\', outcomes = drinking_outcomes, data = data, modelNum = ', m, ')')))
}

combinedLaTeXtableR(analysis = 'cfa',
                    stats = c('Estimate', 'Pr(>|t|)'), statsNames = c('Estimate', 'P'),
                    outcomes = sf12_outcomes, outcomeNameVector = outcomeNameVector, 
                    basic = TRUE, stars = TRUE)
combinedLaTeXtableR(analysis = 'bcfa',
                    stats = c('Estimate', 'Pr(>|t|)'), statsNames = c('Estimate', 'P'),
                    outcomes = sf12_outcomes, outcomeNameVector = outcomeNameVector, 
                    basic = TRUE, stars = TRUE)  
combinedLaTeXtableR(analysis = 'bf',
                    stats = c('Estimate', 'Pr(>|t|)'), statsNames = c('Estimate', 'P'),
                    outcomes = sf12_outcomes, outcomeNameVector = outcomeNameVector, 
                    basic = TRUE, stars = TRUE)

bodyBuildeR(mat = c('gh_cfa_dsm5', 'depressed_cfa_dsm5'), matModelName = 'dsm5', multiOutcome = TRUE, outcomeNames = c('General Health', 'Depressed'))

combinedLaTeXtableR(analysis = 'cfa', 
                    outcomes = dx_outcomes, outcomeNameVector = outcomeNameVector, 
                    stats = c('Estimate', 'Pr(>|z|)'), statsNames = c('Estimate', 'P'), 
                    stars = TRUE, basic = TRUE)
combinedLaTeXtableR(analysis = 'bcfa', 
                    outcomes = dx_outcomes, outcomeNameVector = outcomeNameVector, 
                    stats = c('Estimate', 'Pr(>|z|)'), statsNames = c('Estimate', 'P'), 
                    stars = TRUE, basic = TRUE)
combinedLaTeXtableR(analysis = 'bf', outcomes = dx_outcomes, 
                    outcomeNameVector = outcomeNameVector, 
                    stats = c('Estimate', 'Pr(>|z|)'), statsNames = c('Estimate', 'P'), 
                    stars = TRUE, basic = TRUE)

combinedLaTeXtableR(analysis = 'cfa', 
                    outcomes = drinking_outcomes, outcomeNameVector = outcomeNameVector, 
                    stats = c('Estimate', 'Pr(>|t|)'), statsNames = c('Estimate', 'P'), 
                    stars = TRUE, basic = TRUE)  
combinedLaTeXtableR(analysis = 'bcfa', 
                    outcomes = drinking_outcomes, outcomeNameVector = outcomeNameVector, 
                    stats = c('Estimate', 'Pr(>|t|)'), statsNames = c('Estimate', 'P'), 
                    stars = TRUE, basic = TRUE)  
combinedLaTeXtableR(analysis = 'bf', 
                    outcomes = drinking_outcomes, outcomeNameVector = outcomeNameVector, 
                    stats = c('Estimate', 'Pr(>|t|)'), statsNames = c('Estimate', 'P'), 
                    stars = TRUE, basic = TRUE)