In enaY15/TabulationAutomation: Auomated analysis and tabulation of results

library(TabulationAutomation)

#=============================================================================================================#
##### Specify Regression Analyses Outcome Sets ####
#=============================================================================================================#

TabulationAutomation::addModelContent(setName = 'drinking_outcomes', 
                list = 'outcomesList', 
                setContent = c('MAXDRINKS', 'DRINKFREQ', 'USUALAMT', 'BINGE', 'INTOX'),
                setNameAttribute = c('Maximum Drinks', 'Drinking Frequency', 'Usual Amount', 'Binge Frequency', 'Intoxication Frequency'))

TabulationAutomation::addModelContent(setName = 'sf12_outcomes', 
                list = 'outcomesList', 
                setContent = c('gh', 'calm', 'depressed', 'accomplishedless', 'carefulless') ,
                setNameAttribute = c('General Health', 'Calm/Peaceful', 'Down/Depressed', 'Less Accomplished', 'Less Careful'))

TabulationAutomation::addModelContent(setName = 'dx_outcomes', 
                list = 'outcomesList', 
                setContent = c('pyptsd', 'lptsd', 'pymddisorder' , 'lmddisorder', 
                               'bpddx1', 'pygadind', 'lgadind', 'pyspeind', 'lspeind', 
                               'pypanicind', 'lpanicind', 'pyaud5', 'lifeaud5', 
                               'lmaud5', 'lnicdep5'),
                setNameAttribute = c('Past Year PTSD', 'Lifetime PTSD', 'Past Year MDD', 'Life MDD', 
                                     'BPD', 'Past Year GAD', 'Life GAD', 'Past Year Phobia', 'Life Phobia', 
                                     'Past Year Panic', 'Life Panic', 'Past Year AUD', 'Life AUD', 
                                     'Last Month AUD', 'Life Nicotine Dependence'))


#=============================================================================================================#
##### Specify Regression Analyses Predictor Sets (e.g. Sets of Factors in Different Models) ####
#=============================================================================================================#

TabulationAutomation::addModelContent(setName = 'dsm5', 
                list = 'modelContentsList',
                setContent = c('INTRUSIONS', 'AVOIDANCE', 'COGMOOD', 'HYPERAROUSAL'),
                setNameAttribute = c('Intrusions', 'Avoidance', 'Cognition/Mood', 'Hyperarousal'))

TabulationAutomation::addModelContent(setName = 'dysphoria', 
                list = 'modelContentsList',
                setContent = c('INTRUSIONS', 'AVOIDANCE', 'DYSPHORIA', 'HYPERAROUSAL'),
                setNameAttribute = c('Intrusions', 'Avoidance', 'Dysphoria', 'Hyperarousal'))

TabulationAutomation::addModelContent(setName = 'dysphoricarousal', 
                list = 'modelContentsList',
                setContent = c('INTRUSIONS', 'AVOIDANCE', 'COGMOOD', 'DYSPHAROUSAL', 'ANXAROUSAL'),
                setNameAttribute = c('Intrusions', 'Avoidance', 'Dysphoria', 'Dysphoric Arousal', 'Anxious Arousal'))

TabulationAutomation::addModelContent(setName = 'anhedonia', 
                list = 'modelContentsList',
                setContent = c('INTRUSIONS', 'AVOIDANCE', 'NEGAFF', 'ANHEDONIA', 'DYSPHAROUSAL', 'ANXAROUSAL'),
                setNameAttribute = c('Intrusions', 'Avoidance', 'Negative Affect', 'Anhedonia', 'Dysphoric Arousal', 'Anxious Arousal'))

TabulationAutomation::addModelContent(setName = 'externalizing',
                list = 'modelContentsList', 
                setContent = c('INTRUSIONS', 'AVOIDANCE', 'NUMBING', 'EXTERNALIZING', 'DYSPHAROUSAL', 'ANXAROUSAL'), 
                setNameAttribute = c('Intrusions', 'Avoidance', 'Numbing', 'Externalizing', 'Dysphoric Arousal', 'Anxious Arousal'))

TabulationAutomation::addModelContent(setName = 'hybrid', 
                list = 'modelContentsList',
                setContent = c('INTRUSIONS', 'AVOIDANCE', 'NEGAFF', 'ANHEDONIA', 'EXTERNALIZING', 'DYSPHAROUSAL', 'ANXAROUSAL'),
                setNameAttribute = c('Intrusions', 'Avoidance', 'Negative Affect', 'Anhedonia', 'Externalizing', 'Dysphoric Arousal', 'Anxious Arousal'))

#=============================================================================================================#
##### Make More Variables from Specifications ####
#=============================================================================================================#

##### NESARC 3: Indicators and Factors ####

indicatorModels <- matrix(c('B1', 'In', 'In', 'In', 'In', 'In', 'In',
                            'B2', 'In', 'In', 'In', 'In', 'In', 'In',
                            'B3', 'In', 'In', 'In', 'In', 'In', 'In',
                            'B4', 'In', 'In', 'In', 'In', 'In', 'In',
                            'B5', 'In', 'In', 'In', 'In', 'In', 'In',
                            'C1', 'Av', 'Av', 'Av', 'Av', 'Av', 'Av',
                            'C2', 'Av', 'Av', 'Av', 'Av', 'Av', 'Av',
                            'D1', 'CM', 'Dy', 'CM', 'Nu', 'Na', 'Na',
                            'D2', 'CM', 'Dy', 'CM', 'Nu', 'Na', 'Na',
                            'D3', 'CM', 'Dy', 'CM', 'Nu', 'Na', 'Na',
                            'D4', 'CM', 'Dy', 'CM', 'Nu', 'Na', 'Na',
                            'D5', 'CM', 'Dy', 'CM', 'Nu', 'An', 'An',
                            'D6', 'CM', 'Dy', 'CM', 'Nu', 'An', 'An',
                            'D7', 'CM', 'Dy', 'CM', 'Nu', 'An', 'An',
                            'E1', 'HA', 'Dy', 'DA', 'EB', 'DA', 'EB',
                            'E2', 'HA', 'Dy', 'DA', 'EB', 'DA', 'EB',
                            'E3', 'HA', 'HA', 'AA', 'AA', 'AA', 'AA',
                            'E4', 'HA', 'HA', 'AA', 'AA', 'AA', 'AA',
                            'E5', 'HA', 'Dy', 'DA', 'DA', 'DA', 'DA',
                            'E6', 'HA', 'Dy', 'DA', 'DA', 'DA', 'DA'), nrow = 20, byrow = TRUE)

indicatorModels <- gsub(x = indicatorModels, pattern = 'In', replacement = 'INTRUSIONS', fixed = TRUE)
indicatorModels <- gsub(x = indicatorModels, pattern = 'Av', replacement = 'AVOIDANCE', fixed = TRUE)
indicatorModels <- gsub(x = indicatorModels, pattern = 'CM', replacement = 'COGMOOD', fixed = TRUE)
indicatorModels <- gsub(x = indicatorModels, pattern = 'HA', replacement = 'HYPERAROUSAL', fixed = TRUE)
indicatorModels <- gsub(x = indicatorModels, pattern = 'Dy', replacement = 'DYSPHORIA', fixed = TRUE)
indicatorModels <- gsub(x = indicatorModels, pattern = 'Nu', replacement = 'NUMBING', fixed = TRUE)
indicatorModels <- gsub(x = indicatorModels, pattern = 'AA', replacement = 'ANXAROUSAL', fixed = TRUE)
indicatorModels <- gsub(x = indicatorModels, pattern = 'Na', replacement = 'NEGAFF', fixed = TRUE)
indicatorModels <- gsub(x = indicatorModels, pattern = 'An', replacement = 'ANHEDONIA', fixed = TRUE)
indicatorModels <- gsub(x = indicatorModels, pattern = 'EB', replacement = 'EXTERNALIZING', fixed = TRUE)
indicatorModels[which(indicatorModels == 'DA')] <- 'DYSPHAROUSAL'
indicatorModels <- data.frame(indicatorModels)
colnames(indicatorModels) <- c('mplusSymptom', 'dsm5', 'dysphoria', 'dysphoricarousal', 'externalizing', 'anhedonia', 'hybrid')

View(indicatorModels)

factorContentsList <- vector(mode = 'list', length(modelContentsList))
names(factorContentsList) <- names(modelContentsList)

for (m in 1:length(modelContentsList)) {
  factorContentsList[[names(modelContentsList)[m]]] <- vector(mode = 'list', length(modelContentsList[[names(modelContentsList)[m]]]))
  names(factorContentsList[[m]]) <- modelContentsList[[m]]
  model <- names(factorContentsList)[m] 
  for (f in 1:length(factorContentsList[[m]])) {
    factor <- names(factorContentsList[[m]][f]) 
    factorContentsList[[model]][factor] <- paste0(indicatorModels[which(indicatorModels[, model] == factor), 'mplusSymptom'], collapse = ' ')
  }
}

#=============================================================================================================#
##### Make More Variables from Specifications ####
#=============================================================================================================#

outcomeNames <- extractModelInfo(source = outcomesList, what = "vectorElements", topElementsSharedEnding = "_outcomes")

### Split up list of models into list (elements = list for each model) of lists (elements = vector for each factor) of vectors (elements = factor indicators)
  # Can remove dashes and spaces to split up strings of MPlus syntax
modelFactorIndicators_listolists <- extractModelInfo(factorContentsList, what = "mplusVariables")
print(modelFactorIndicators_listolists)
  # Example if just want to do for one model
print(extractModelInfo(factorContentsList[['hybrid']], what = "mplusVariables"))

n = 100

data <- data.frame(matrix(nrow = n, ncol = length(unlist(outcomesList))))
colnames(data) <- names(unlist(outcomesList))
data$row <- 1:dim(data)[1]


## Simulate SF-12 
    sf12covMat <- matrix(nrow = length(sf12_outcomes), ncol = length(sf12_outcomes))
    simulatedCorrs <- rnorm((length(sf12_outcomes)^2-length(sf12_outcomes))/2, 0.7, 0.02)
    sf12covMat[upper.tri(sf12covMat)] <- simulatedCorrs # simulated factor correlations store in upper triangular
    sf12covMat[lower.tri(sf12covMat)] <- sf12covMat[upper.tri(sf12covMat)] # fill in rest of symmetric matrix
    diag(sf12covMat) <- 1
    data[, sf12_outcomes] <- MASS::mvrnorm(n, rep(0, length(sf12_outcomes)), sf12covMat)


    baseRate <- rnorm(n = length(dx_outcomes), mean = 0.08, sd = 0.02) # randomly simulate base rates for each disorder
    for (y in 1:length(dx_outcomes)) {
      diagnosedPeople <- sample(x = data$row, floor(baseRate[y]*n)) # randomly pick people (such that correct proportion of sample) to have dx
      data[, dx_outcomes[y]] <- ifelse(data$row %in% diagnosedPeople, 1, 0)
    }


factorModelNameVector <- names(modelContentsList)

  for (m in 1:length(modelContentsList)) {

    eval(parse(text = paste0('data_cfa_', factorModelNameVector[m], ' <- data.frame(matrix(NA, nrow = ', n, ', ncol = ',  length(modelContentsList[[m]]), '))')))
    colnames <- modelContentsList[[m]]
    eval(parse(text = paste0('colnames(data_cfa_', factorModelNameVector[m], ') <- colnames')))
    eval(parse(text = paste0('data_cfa_', factorModelNameVector[m], ' <- data.frame(c(1:', n, '), data_cfa_', factorModelNameVector[m], ')')))
    eval(parse(text = paste0('colnames(data_cfa_', factorModelNameVector[m], ')[1] <- \'ID\'')))
    eval(parse(text = paste0('data_cfa_', factorModelNameVector[m], ' <- data.frame(data_cfa_', factorModelNameVector[m], ', data)')))

    n.factors <- length(modelContentsList[[m]])

    ## Simulate factor correlation matrix
    factorCovMat <- matrix(nrow = n.factors, ncol = n.factors)
    simulatedCorrs <- rnorm((n.factors^2-n.factors)/2, 0.8, 0.02)
    factorCovMat[upper.tri(factorCovMat)] <- simulatedCorrs # simulated factor correlations store in upper triangular
    factorCovMat[lower.tri(factorCovMat)] <- factorCovMat[upper.tri(factorCovMat)] # fill in rest of symmetric matrix
    diag(factorCovMat) <- 1

    ## Simulate factor scores
    eval(parse(text = paste0('data_cfa_', factorModelNameVector[m], '[, ', paste0('c(', paste0('\'', modelContentsList[[m]], collapse = '\','), '\')'), '] <- MASS::mvrnorm(', n, ', rep(0, n.factors), factorCovMat)')))

  }

dataCols <- paste0(paste0("b", 1:5, collapse = ' '), ' ', paste0("c", 1:2, collapse = ' '), ' ', paste0("d", 1:7, collapse = ' '), ' ', paste0("e", 1:6, collapse = ' '), sep = ' ')

factorAnalyzeR(mplusDataFileColumnNames = dataCols,
               factorIndicatorsList = factorContentsList, modelFactorsList = modelContentsList,
               templateName = 'testTemplate', mplusOutputDirectory = getwd(),
               factorVarSpec = 'free', mplusDataFilePath = "C:/Users/ENA/Documents/GitHub/PTSDn2/n2.csv",
               categoricalInds = "all", na.string = "-9999",
               run = FALSE)

factorAnalyzeR(mplusDataFileColumnNames = dataCols,
               factorIndicatorsList = factorContentsList, modelFactorsList = modelContentsList,
               templateName = 'testTemplate1_fixed', mplusOutputDirectory = getwd(),
               factorVarSpec = 'standardized', mplusDataFilePath = "C:/Users/ENA/Documents/GitHub/PTSDn2/n2.csv",
               categoricalInds = "all", na.string = "-9999",
               run = FALSE)

factorAnalyzeR(mplusDataFileColumnNames = dataCols,
               factorIndicatorsList = factorContentsList, modelFactorsList = modelContentsList,
               templateName = 'testTemplate_noCat_oneFixedAllMods', mplusOutputDirectory = getwd(),
               factorVarSpec = c('@5', rep('free', length(factorContentsList)-1)),
               mplusDataFilePath = "C:/Users/ENA/Documents/GitHub/PTSDn2/n2.csv",
               categoricalInds = NULL, na.string = "-9999",
               run = FALSE)

factorAnalyzeR(mplusDataFileColumnNames = dataCols,
               factorIndicatorsList = factorContentsList, modelFactorsList = modelContentsList,
               templateName = 'testTemplate_someCat', mplusOutputDirectory = getwd(),
               factorVarSpec = 'free',
               mplusDataFilePath = "C:/Users/ENA/Documents/GitHub/PTSDn2/n2.csv",
               categoricalInds = 'a1 a2 b5', na.string = "-9999",
               run = FALSE)

fVarSpecMixed <- vector(mode = 'list', length(modelContentsList))
names(fVarSpecMixed) <- names(modelContentsList)
for (m in 1:(length(modelContentsList))) fVarSpecMixed[[m]] <- rep('free', length(modelContentsList[[m]]))
TabulationAutomation::addModelContent(setName = 'hybrid', 
                list = 'fVarSpecMixed',
                setContent = c("free", "free", "@1", "@1", "free", "free", "free"),
                setNameAttribute = modelContentsList[['hybrid']])

factorAnalyzeR(mplusDataFileColumnNames = dataCols,
               factorIndicatorsList = factorContentsList, modelFactorsList = modelContentsList,
               templateName = 'testTemplate_oneFixedSomeMods', mplusOutputDirectory = getwd(),
               factorVarSpec = fVarSpecMixed,
               mplusDataFilePath = "C:/Users/ENA/Documents/GitHub/PTSDn2/n2.csv",
               categoricalInds = "all", na.string = "-9999",
               run = FALSE)

outputList <- TabulationAutomation::autoAnalyzeR(corrReturn = TRUE, modelPredictorNamesList = modelContentsList, outcomeVector = sf12_outcomes, outcomesSetName = 'sf12', analysis = 'cfa', minDecimals = 2, n.roundDigits = 2)