R/Main.R

Defines functions makeFriendlyNames execute

Documented in execute

# Copyright 2020 Observational Health Data Sciences and Informatics
#
# This file is part of dGEMcovid
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

#' Execute the Study
#'
#' @details
#' This function executes the dGEMcovid Study.
#' 
#' @param databaseDetails      The connection details and OMOP CDM details. Created using \code{PatientLevelPrediction::createDatabaseDetails}.
#' @param siteId               The name of your site (can be the university name, site with a number, etc.) tnis needs to be shared with the study administrator
#' @param outputFolder         Name of local folder to place results - make sure control is in here; make sure to use forward slashes
#'                             (/). Do not use a folder on a network drive since this greatly impacts
#'                             performance.
#' @param createCohorts        Create the cohortTable table with the target population and outcome cohorts?
#' @param createData           Create the labelled data set (this is required before runAnalysis)
#' @param createPlot           Plot the number of covid hospitalizations by month 
#' @param sampleSize           The number of patients in the target cohort to sample (if NULL uses all patients)
#' @param createControl        (for the lead site only) Run this code to create the study control
#' @param siteIds              (for the lead site only) vector with the names of all the sites contributing to the study. Required when createControl = TRUE
#' @param runInitialize        Runs the initialization step - this require downloading the json control from https://pda-ota.pdamethods.org/ 
#' @param runDerive            Runs the derive step - this require downloading the updated json control from https://pda-ota.pdamethods.org/ 
#' @param runEstimate          Runs the estimate step - this require downloading the updated json control from https://pda-ota.pdamethods.org/ 
#' @param runSynthesize        Once the site estimates are returned, it is now possible to apply each model to the data to calculate
#'                             predictions.  This step requires downloading the updated json control from https://pda-ota.pdamethods.org/                        
#' @param verbosity            Sets the level of the verbosity. If the log level is at or higher in priority than the logger threshold, a message will print. The levels are:
#'                                         \itemize{
#'                                         \item{DEBUG}{Highest verbosity showing all debug statements}
#'                                         \item{TRACE}{Showing information about start and end of steps}
#'                                         \item{INFO}{Show informative information (Default)}
#'                                         \item{WARN}{Show warning messages}
#'                                         \item{ERROR}{Show error messages}
#'                                         \item{FATAL}{Be silent except for fatal errors}
#'                                         }  
#'  @param alphaStart     The start date for the alpha wave in the data in "yyyymmdd" format    
#'  @param alphaEnd     The end date for the alpha wave in the data in "yyyymmdd" format  
#'  @param deltaStart     The start date for the delta wave in the data in "yyyymmdd" format  
#'  @param deltaEnd     The end date for the detla wave in the data in "yyyymmdd" format                                                                 
#'
#' @examples
#' \dontrun{
#' connectionDetails <- createConnectionDetails(dbms = "postgresql",
#'                                              user = "joe",
#'                                              password = "secret",
#'                                              server = "myserver")
#'
#' execute(databaseDetails,
#'         siteId = 'site_1',
#'         outputFolder = "c:/temp/study_results", 
#'         createCohorts = T,
#'         createData = T,
#'         sampleSize = 10000,
#'         createControl = F,
#'         siteIds = c('site_1', 'site_2'),
#'         leadSiteNextStep = F,
#'         runAnalysis = F,
#'         runSynthesize = F,
#'         verbosity = "INFO"
#'         )
#' }
#'
#' @export
execute <- function(
  databaseDetails,
  siteId = 'site_name',
  outputFolder,
  createCohorts = F,
  createData = F,
  createPlot = createData,
  sampleSize = NULL,
  createControl = F,
  siteIds = '', # only needed if lead site
  leadSiteNextStep = F,
  runInitialize = F,
  runDerive = F,
  runEstimate = F,
  runSynthesize = F,
  verbosity = "INFO",
  alphaStart = '20201101',
  alphaEnd = '20210201',
  deltaStart = '20210701',
  deltaEnd = '20211001'
) {
  
  
  dates <- list(
    alpha = list(startDate = alphaStart, endDate = alphaEnd),
    delta = list(startDate = deltaStart, endDate = deltaEnd)
  )
  
  waves <- c('alpha','delta')
  
  if (!file.exists(outputFolder)){
    dir.create(outputFolder, recursive = TRUE)
  }
  
  website <- 'https://pda-ota.pdamethods.org/'
  
  # load the analysis
  analysisListFile <- system.file(
    "settings",
    "analysisList.json",
    package = "dGEMcovid"
  )
  
  if(file.exists(analysisListFile)){
    analysisList <- loadJson(analysisListFile)
    analysisList <- fromJsonFormat(analysisList)
    analysisList$studySettings$covariateSettings <- addCohortSettings(
      covariateSettings = analysisList$studySettings$covariateSettings,
      cohortDatabaseSchema = databaseDetails$cohortDatabaseSchema , 
      cohortTable = databaseDetails$cohortTable
    )
  }else{
    stop('No analysisList available')
  }
  
  if(createCohorts) {
    ParallelLogger::logInfo("Creating cohorts")
    createCohorts(
      databaseDetails = databaseDetails,
      outputFolder = outputFolder
    )
  } 
  
  if(createPlot){
    ParallelLogger::logInfo('Generating date plot')
    generateDatePlot(
      databaseDetails = databaseDetails,
      outputFolder = outputFolder
    )
  }
  
  if(createData){
    
    # save the settings
    
    timeSettings <- data.frame(
      wave = names(dates),
      starts = unlist(lapply(dates, function(x) x$startDate)),
      ends = unlist(lapply(dates, function(x) x$endDate))
    )
    write.csv(timeSettings, file.path(outputFolder, 'timeSettings.csv'), row.names = F)
    
    for(wave in waves){
      
      if(!dir.exists(file.path(outputFolder, wave))){
        dir.create(file.path(outputFolder, wave), recursive = T)
      }
      # restrict to the input date:
      ParallelLogger::logInfo(paste0(wave, " date: ",dates[[wave]]$startDate, "-", dates[[wave]]$endDate ))
      analysisList$studySettings$restrictPlpDataSettings$studyStartDate <- dates[[wave]]$startDate
      analysisList$studySettings$restrictPlpDataSettings$studyEndDate <- dates[[wave]]$endDate
      
      
      ParallelLogger::logInfo("Creating labelled dataset")
      
      databaseDetails$cohortId <- analysisList$studySettings$targetId
      databaseDetails$outcomeIds <- analysisList$studySettings$outcomeId
      
      plpData <- tryCatch(
        {
        PatientLevelPrediction::getPlpData(
        databaseDetails = databaseDetails, 
        covariateSettings = analysisList$studySettings$covariateSettings, 
        restrictPlpDataSettings = analysisList$studySettings$restrictPlpDataSettings
        )}, 
        error = function(e){ParallelLogger::logInfo(e); return(NULL)}
      )
      
      if(!is.null(plpData)){
      labels <- PatientLevelPrediction::createStudyPopulation(
        plpData = plpData, 
        outcomeId = analysisList$studySettings$outcomeId, 
        populationSettings = analysisList$studySettings$populationSettings
      )
      
      # convert to matrix
      
      dataObject <- PatientLevelPrediction::toSparseM(
        plpData = plpData, 
        cohort = labels
      )
      
      #sparse matrix: dataObject$dataMatrix
      #labels: dataObject$labels
      
      columnDetails <- as.data.frame(dataObject$covariateRef)
      
      cnames <- columnDetails$covariateName[order(columnDetails$columnId)]
      
      ipMat <- as.matrix(dataObject$dataMatrix)
      ipdata <- as.data.frame(ipMat)
      colnames(ipdata) <-  makeFriendlyNames(cnames)
      ipdata$outcome <- dataObject$labels$outcomeCount
      
      # modify the covariates
      # Charlson comorbidity categories: 0-1, 2-4, and 5
      if('Charlson_index___Romano_adaptation' %in% colnames(ipdata)){
        ParallelLogger::logInfo('Processing Charlson index into categories')
        # Charlson comorbidity categories: 0-1, 2-4, and 5
        
        # make 0-1 the reference
        ##ipdata$Charlson_index___Romano_adaptation0_1 <- rep(0, nrow(ipdata))
        ##ipdata$Charlson_index___Romano_adaptation0_1[ipdata$Charlson_index___Romano_adaptation <= 1] <- 1
        
        ipdata$Charlson_index___Romano_adaptation2_4 <- rep(0, nrow(ipdata))
        ipdata$Charlson_index___Romano_adaptation2_4[ipdata$Charlson_index___Romano_adaptation <= 4 & ipdata$Charlson_index___Romano_adaptation >=2] <- 1
        
        ipdata$Charlson_index___Romano_adaptation5p <- rep(0, nrow(ipdata))
        ipdata$Charlson_index___Romano_adaptation5p[ipdata$Charlson_index___Romano_adaptation >= 5] <- 1
        
        ipdata <- ipdata[!colnames(ipdata) == 'Charlson_index___Romano_adaptation']
      }
      
      # Age categories: 18-65, 65-80, and 80
      if('age_in_years' %in% colnames(ipdata)){
        ParallelLogger::logInfo('Processing age into categories')
        
        # make 18-64 the reference
        ##ipdata$age18_64 <- rep(0, nrow(ipdata))
        ##ipdata$age18_64[ipdata$age_in_years >= 0 & ipdata$age_in_years <= 64] <- 1
        
        ipdata$age65_80 <- rep(0, nrow(ipdata))
        ipdata$age65_80[ipdata$age_in_years >= 65 & ipdata$age_in_years <= 79] <- 1
        
        ipdata$age80p <- rep(0, nrow(ipdata))
        ipdata$age80p[ipdata$age_in_years >= 80] <- 1
        
        ipdata <- ipdata[,colnames(ipdata) != 'age_in_years']
      }
      
      # make the female gender the reference
      ipdata <- ipdata[,colnames(ipdata) != 'gender___FEMALE']
      
      # save the data:
      utils::write.csv(
        x = ipdata, 
        file = file.path(outputFolder,wave,'data.csv'), 
        row.names = F
      )
      
      # add the data covariates
      
      if(!is.null(analysisList$studySettings$dataCovariateSettings)){
        dataSummary <- do.call(
          'dataSummary', 
          list(
            dataCovariateSettings = analysisList$studySettings$dataCovariateSettings,
            databaseDetails = databaseDetails
          )
        )
        
        utils::write.csv(
          x = dataSummary, 
          file = file.path(outputFolder,wave, 'dataSummary.csv'), 
          row.names = F
        )
        
      } 
        
      } # end if not null
      
    } # end wave
    
  }
  
  # Step 1: lead site creates the controls
  if(createControl){
    
    for(wave in waves){
      ParallelLogger::logInfo(paste0('Creating the control settings for wave ', wave))
      
      #check the data exist to get the names
      if(file.exists(file.path(outputFolder, wave, 'dataSummary.csv'))){
        dataSummary <- utils::read.csv(file.path(outputFolder, wave, 'dataSummary.csv'))
      } else{
        dataSummary <- NULL
      }
      
      #check the data exist to get the names
      if(file.exists(file.path(outputFolder, wave, 'data.csv'))){
        data <- utils::read.csv(file.path(outputFolder, wave, 'data.csv'))
      } else{
        stop('Please generate data before creating control')
      }
      
      control <- list(
        project_name = paste(analysisList$packageName, wave, sep = '_'),
        step = 'initialize',
        sites = siteIds,
        heterogeneity = analysisList$studySettings$control$heterogeneity,
        model = analysisList$studySettings$control$model,
        family = analysisList$studySettings$control$family,
        outcome = "outcome",
        variables = colnames(data)[colnames(data)!='outcome'],
        variables_site_level = names(dataSummary),
        optim_maxit = analysisList$studySettings$control$optim_maxit,
        lead_site = siteId,
        upload_date = as.character(Sys.time()) 
      )
      
      control <- jsonlite::toJSON(
        x = control, 
        pretty = T, 
        digits = 23, 
        auto_unbox=TRUE, 
        null = "null"
      )
      
      write(control, file.path(outputFolder, wave, 'control.json'))
      
      ParallelLogger::logInfo(
        paste0(
          'Initial control created at: ',
          file.path(outputFolder, wave, 'control.json'),
          ' please upload this to ', website
          
        )
      )
      
    }
    
  }
  
  if(runInitialize | runDerive | runEstimate | runSynthesize){
    # json control file
    jsonFileLocation <- list(
      alpha = file.path(outputFolder,'alpha', 'control.json'),
      delta = file.path(outputFolder,'delta', 'control.json')
      )
    
  }

  
  if(runInitialize){
    
    for(wave in waves){
    
    control <- tryCatch(
      {readChar(jsonFileLocation[[wave]], file.info(jsonFileLocation[[wave]])$size)},
      error= function(cond) {
        ParallelLogger::logInfo('Issue with loading json file...');
        ParallelLogger::logError(cond);
        return(NULL)
      })
    
    if(!is.null(control)){
      
      control <- jsonlite::fromJSON(control)
      
      # check control name 
      if(length(grep(wave, control$project_name))==0){
        ParallelLogger::logWarn('Control project_name suggests incorrect wave')
      }
      
      ipdata <- utils::read.csv(file.path(outputFolder, wave, 'data.csv'))
      
      # get summary variables
      if(file.exists(file.path(outputFolder, wave, 'dataSummary.csv'))){
        dataSummary <- utils::read.csv(file.path(outputFolder, wave, 'dataSummary.csv'))
      } else{
        dataSummary <- NULL
      }
      
      
      if(!is.null(control)){
        if(control$step == 'initialize'){
          
          ParallelLogger::logInfo(paste0('At step ', control$step))
          pda::pda(
            ipdata = ipdata, 
            hosdata = dataSummary,
            site_id = siteId, 
            dir = file.path(outputFolder, wave)
          )
          
          ParallelLogger::logInfo(paste0('result json ready to check in ', file.path(outputFolder, wave)))
          ParallelLogger::logInfo(paste0('If you are happy to share, please upload this to ', website))
          
        } # control exists
      } else{
        ParallelLogger::logInfo('Control is not for initialize - please check control stage')
      }
      
    } else {
      ParallelLogger::logInfo(paste0('No control for wave ', wave))
      }
    
    
  } # end wave
    
  }
  
  if(runDerive){
    
    for(wave in waves){
      
      control <- tryCatch(
        {readChar(jsonFileLocation[[wave]], file.info(jsonFileLocation[[wave]])$size)},
        error= function(cond) {
          ParallelLogger::logInfo('Issue with loading json file...');
          ParallelLogger::logError(cond);
          return(NULL)
        })
      
      
      if(!is.null(control)){
        
        control <- jsonlite::fromJSON(control)
        
        # check control name 
        if(length(grep(wave, control$project_name))==0){
          ParallelLogger::logWarn('Control project_name suggests incorrect wave')
        }
        
        ipdata <- utils::read.csv(file.path(outputFolder, wave, 'data.csv'))
        
        # get summary variables
        if(file.exists(file.path(outputFolder, wave, 'dataSummary.csv'))){
          dataSummary <- utils::read.csv(file.path(outputFolder, wave, 'dataSummary.csv'))
        } else{
          dataSummary <- NULL
        }
        
        if(!is.null(control)){
          if(control$step == 'derive'){
            
            ParallelLogger::logInfo(paste0('At step ', control$step))
            pda::pda(
              ipdata = ipdata, 
              hosdata = dataSummary,
              site_id = siteId, 
              dir = file.path(outputFolder, wave)
            )
            
            ParallelLogger::logInfo(paste0('result json ready to check in ', file.path(outputFolder, wave)))
            ParallelLogger::logInfo(paste0('If you are happy to share, please upload this to ', website))
            
          } # control exists
        } else{
          ParallelLogger::logInfo('Control is not for derive - please check control stage')
        }
        
      } else {
        ParallelLogger::logInfo(paste0('No control for wave ', wave))
      }
      
    } # end wave
    
  }
  
  if(runEstimate){
    
    for(wave in waves){
      
      control <- tryCatch(
        {readChar(jsonFileLocation[[wave]], file.info(jsonFileLocation[[wave]])$size)},
        error= function(cond) {
          ParallelLogger::logInfo('Issue with loading json file...');
          ParallelLogger::logError(cond);
          return(NULL)
        })
      
      if(!is.null(control)){
        
        control <- jsonlite::fromJSON(control)
        
        # check control name 
        if(length(grep(wave, control$project_name))==0){
          ParallelLogger::logWarn('Control project_name suggests incorrect wave')
        }
        
        ipdata <- utils::read.csv(file.path(outputFolder, wave, 'data.csv'))
        
        # get summary variables
        if(file.exists(file.path(outputFolder, wave, 'dataSummary.csv'))){
          dataSummary <- utils::read.csv(file.path(outputFolder, wave, 'dataSummary.csv'))
        } else{
          dataSummary <- NULL
        }
        
        if(!is.null(control)){
          if(control$step == 'estimate'){
            
            ParallelLogger::logInfo(paste0('At step ', control$step))
            pda::pda(
              ipdata = ipdata, 
              hosdata = dataSummary,
              site_id = siteId, 
              dir = file.path(outputFolder, wave)
            )
            
            ParallelLogger::logInfo(paste0('result json ready to check in ', file.path(outputFolder, wave)))
            ParallelLogger::logInfo(paste0('If you are happy to share, please upload this to ', website))
            
          } # control exists
        } else{
          ParallelLogger::logInfo('Control is not for estimate - please check control stage')
        }
        
      } else {
        ParallelLogger::logInfo(paste0('No control for wave ', wave))
      }
      
    } # end wave
    
  }
  
  if(runSynthesize){
    #ipdata <- utils::read.csv(file.path(outputFolder, 'data'))
    # do something with this?
    #cbind(b.odal=fit.odal$btilde,
    #  sd.odal=sqrt(diag(solve(fit.odal$Htilde)/nrow(ipdata))))
  
  }
  
  invisible(NULL)
}


makeFriendlyNames <- function(columnNames){
  
  columnNames <- gsub("[[:punct:]]", " ", columnNames)
  columnNames <- gsub(" ", "_", columnNames)
  return(columnNames)
  
}
ohdsi-studies/dGEMcovid documentation built on May 10, 2023, 3:38 a.m.