Introduction

In this vignette we focus on running several different analyses on several exposure-outcome pairs This can be useful when we want to explore the sensitivity to analyses choices, include controls, or run an experiment similar to the OMOP experiment to empirically identify the optimal analysis choices for a particular research question.

This vignette assumes you are already familiar with the SelfControlledCaseSeries package and are able to perform single studies. We will walk through all the steps needed to perform an exemplar set of analyses, and we have selected the well-studied topic of the effect of nonsteroidal anti-inflammatory drugs (NSAIDs) on gastrointestinal (GI) bleeding-related hospitalization. For simplicity, we focus on one NSAID: diclofenac. We will execute various variations of an analysis for the primary exposure pair and a large set of negative control exposures.

General approach

The general approach to running a set of analyses is that you specify all the function arguments of the functions you would normally call, and create sets of these function arguments. The final outcome models as well as intermediate data objects will all be saved to disk for later extraction.

An analysis will be executed by calling these functions in sequence:

  1. getDbSccsData()
  2. createSccsEraData()
  3. fitSccsModel()

When you provide several analyses to the SelfControlledCaseSeries package, it will determine whether any of the analyses and exposure-outcome pairs have anything in common, and will take advantage of this fact. For example, if we specify several exposure-outcome pairs with the same outcome, the data for the outcome will be extracted only once.

The function arguments you need to define have been divided into four groups:

  1. Hypothesis of interest: arguments that are specific to a hypothesis of interest, in the case of the self-controlled case series this is a combination of exposure and outcome.
  2. Analyses: arguments that are not directly specific to a hypothesis of interest, such as the washout window, whether to adjust for age and seasonality, etc.
  3. Arguments that are the output of a previous function in the SelfControlledCaseSeries package, such as the sccsData argument of the createSccsEraData function. These cannot be specified by the user.
  4. Arguments that are specific to an environment, such as the connection details for connecting to the server, and the name of the schema holding the CDM data.

Preparation for the example

We need to tell R how to connect to the server where the data are. SelfControlledCaseSeries uses the DatabaseConnector package, which provides the createConnectionDetails function. Type ?createConnectionDetails for the specific settings required for the various database management systems (DBMS). For example, one might connect to a PostgreSQL database using this code:

connectionDetails <- createConnectionDetails(dbms = "postgresql", 
                                             server = "localhost/ohdsi", 
                                             user = "joe", 
                                             password = "supersecret")

outputFolder <- "s:/temp/sccsVignette2"

cdmDatabaseSchema <- "my_cdm_data"
cohortDatabaseSchema <- "my_cohorts"
cdmVersion <- "5"

The last three lines define the cdmDatabaseSchema and cohortDatabaseSchema variables, as well as the CDM version. We’ll use these later to tell R where the data in CDM format live, where we want to store the (outcome) cohorts, and what version CDM is used. Note that for Microsoft SQL Server, databaseschemas need to specify both the database and the schema, so for example cdmDatabaseSchema <- "my_cdm_data.dbo".

We also need to prepare our exposures and outcomes of interest. The drug_era table in the OMOP Common Data Model already contains prespecified cohorts of users at the ingredient level, so we will use that for the exposures. For the outcomes, we want to restrict our analysis only to those events that are recorded in an inpatient setting, so we will need to create a custom cohort table. For this example, we are only interested in GI bleed (concept ID 192671) .

We create a text file called vignette.sql with the following content:

/***********************************
File vignette.sql 
***********************************/

IF OBJECT_ID('@cohortDatabaseSchema.@outcomeTable', 'U') IS NOT NULL
  DROP TABLE @cohortDatabaseSchema.@outcomeTable;

SELECT 1 AS cohort_definition_id,
    condition_start_date AS cohort_start_date,
    condition_end_date AS cohort_end_date,
    condition_occurrence.person_id AS subject_id
INTO @cohortDatabaseSchema.@outcomeTable
FROM @cdmDatabaseSchema.condition_occurrence
INNER JOIN @cdmDatabaseSchema.visit_occurrence
    ON condition_occurrence.visit_occurrence_id = visit_occurrence.visit_occurrence_id
WHERE condition_concept_id IN (
        SELECT descendant_concept_id
        FROM @cdmDatabaseSchema.concept_ancestor
        WHERE ancestor_concept_id = 192671 -- GI - Gastrointestinal haemorrhage
        )
    AND visit_occurrence.visit_concept_id IN (9201, 9203);

Note for CDM V4 visit_concept_id should be place_of_service_concept_id, and cohort_definition_id should be cohort_concept_id.

This is parameterized SQL which can be used by the SqlRender package. We use parameterized SQL so we do not have to pre-specify the names of the CDM and result schemas. That way, if we want to run the SQL on a different schema, we only need to change the parameter values; we do not have to change the SQL code. By also making use of translation functionality in SqlRender, we can make sure the SQL code can be run in many different environments.

library(SqlRender)
sql <- readSql("vignette.sql")
sql <- render(sql,
              cdmDatabaseSchema = cdmDatabaseSchema, 
              cohortDatabaseSchema = cohortDatabaseSchema)
sql <- translate(sql, targetDialect = connectionDetails$dbms)

connection <- connect(connectionDetails)
executeSql(connection, sql)

In this code, we first read the SQL from the file into memory. In the next line, we replace the two parameter names with the actual values. We then translate the SQL into the dialect appropriate for the DBMS we already specified in the connectionDetails. Next, we connect to the server, and submit the rendered and translated SQL.

Specifying hypotheses of interest

The first group of arguments define the exposure and outcome. Here we demonstrate how to create a list of exposure-outcome pairs:

negativeControls <- c(705178, 705944, 710650, 714785, 719174, 719311, 735340, 742185, 
                      780369, 781182, 924724, 990760, 1110942, 1111706, 1136601, 
                      1317967, 1501309, 1505346, 1551673, 1560278, 1584910, 19010309, 
                      40163731)
diclofenac <- 1124300
ppis <- c(911735, 929887, 923645, 904453, 948078, 19039926)

exposureOutcomeList <- list()
for (exposureId in c(diclofenac, negativeControls)){
  exposureOutcome <- createExposureOutcome(exposureId = exposureId,
                                           outcomeId = 1,
                                           prophylactics = ppis)
  exposureOutcomeList[[length(exposureOutcomeList) + 1]] <- exposureOutcome
}

We defined the outcome of interest to be the custom cohort with ID 1 we defined in the SQL above. The exposures include diclofenac (concept ID 1124300) and a large number of negative control exposures. We furthermore specify a custom variable prophylactics that contains the concept IDs of drugs that are used in preventing the outcome of interest, in this case proton pump inhibitors (PPIs).

A convenient way to save exposureOutcomeList to file is by using the saveExposureOutcomeList function, and we can load it again using the loadExposureOutcomeList function.

Specifying analyses

The second group of arguments are not specific to a hypothesis of interest, and comprise the majority of arguments. For each function that will be called during the execution of the analyses, a companion function is available that has (almost) the same arguments. For example, for the fitSccsModel() function there is the createFitSccsModelArgs() function. These companion functions can be used to create the arguments to be used during execution:

getDbSccsDataArgs1 <- createGetDbSccsDataArgs(useCustomCovariates = FALSE,
                                              deleteCovariatesSmallCount = 100,
                                              studyStartDate = "",
                                              studyEndDate = "",
                                              exposureIds = c())

covarEoi <- createCovariateSettings(label = "Exposure of interest",
                                              includeCovariateIds = "exposureId",
                                              start = 0,
                                              end = 0,
                                              addExposedDaysToEnd = TRUE)

createEraDataArgs1 <- createCreateSccsEraDataArgs(naivePeriod = 180,
                                                  firstOutcomeOnly = FALSE,
                                                  covariateSettings = covarEoi)

fitSccsModelArgs <- createFitSccsModelArgs()

Any argument that is not explicitly specified by the user will assume the default value specified in the function. Note that for several arguments for concept or cohort definition IDs we can use the name of a variable in the exposureOutcome objects. In this case, we defined the argument includeCovariateIds to get the value of the exposureId variable.

We can now combine the arguments for the various functions into a single analysis:

sccsAnalysis1 <- createSccsAnalysis(analysisId = 1,
                                    description = "Simplest model",
                                    getDbSccsDataArgs = getDbSccsDataArgs1,
                                    createSccsEraDataArgs = createEraDataArgs1,
                                    fitSccsModelArgs = fitSccsModelArgs)

Note that we have assigned an analysis ID (1) to this set of arguments. We can use this later to link the results back to this specific set of choices. We also include a short description of the analysis.

We can easily create more analyses, for example by incuding adjustments for age and seasonality, or for including other drugs in the model:

covarProp <- createCovariateSettings(label = "Prophylactics",
                                              includeCovariateIds = "prophylactics",
                                              start = 0,
                                              end = 0,
                                              addExposedDaysToEnd = TRUE)

createSccsEraDataArgs2 <- createCreateSccsEraDataArgs(naivePeriod = 180,
                                                      firstOutcomeOnly = FALSE,
                                                      covariateSettings = list(covarEoi,
                                                                               covarProp))

sccsAnalysis2 <- createSccsAnalysis(analysisId = 2,
                                    description = "Including prophylactics",
                                    getDbSccsDataArgs = getDbSccsDataArgs1,
                                    createSccsEraDataArgs = createSccsEraDataArgs2,
                                    fitSccsModelArgs = fitSccsModelArgs)

ageSettings <- createAgeSettings(includeAge = TRUE, ageKnots = 5)

seasonSettings <- createSeasonalitySettings(includeSeasonality = TRUE, seasonKnots = 5)

covarPreExp <- createCovariateSettings(label = "Pre-exposure",
                                       includeCovariateIds = "exposureId",
                                       start = -30,
                                       end = -1)

createSccsEraDataArgs3 <- createCreateSccsEraDataArgs(naivePeriod = 180,
                                                      firstOutcomeOnly = FALSE,
                                                      covariateSettings = list(covarEoi,
                                                                               covarPreExp,
                                                                               covarProp),
                                                      ageSettings = ageSettings,
                                                      seasonalitySettings = seasonSettings,
                                                      eventDependentObservation = TRUE)

sccsAnalysis3 <- createSccsAnalysis(analysisId = 3,
                                    description = "Inc. prop., age, season, pre-exp., censor.",
                                    getDbSccsDataArgs = getDbSccsDataArgs1,
                                    createSccsEraDataArgs = createSccsEraDataArgs3,
                                    fitSccsModelArgs = fitSccsModelArgs)

covarAllDrugs <- createCovariateSettings(label = "Other exposures",
                                         excludeCovariateIds = "exposureId",
                                         stratifyById = TRUE,
                                         start = 1,
                                         end = 0,
                                         addExposedDaysToEnd = TRUE,
                                         allowRegularization = TRUE)

createSccsEraDataArgs4 <- createCreateSccsEraDataArgs(naivePeriod = 180,
                                                      firstOutcomeOnly = FALSE,
                                                      covariateSettings = list(covarEoi,
                                                                               covarPreExp,
                                                                               covarAllDrugs),
                                                      ageSettings = ageSettings,
                                                      seasonalitySettings = seasonSettings,
                                                      eventDependentObservation = TRUE)

sccsAnalysis4 <- createSccsAnalysis(analysisId = 4,
                                    description = "Including all other drugs",
                                    getDbSccsDataArgs = getDbSccsDataArgs1,
                                    createSccsEraDataArgs = createSccsEraDataArgs4,
                                    fitSccsModelArgs = fitSccsModelArgs)

These analyses can be combined in a list:

sccsAnalysisList <- list(sccsAnalysis1, sccsAnalysis2, sccsAnalysis3, sccsAnalysis4)

A convenient way to save sccsAnalysisList to file is by using the saveSccsAnalysisList function, and we can load it again using the loadSccsAnalysisList function.

Exposure and outcome selection strategies

Often we would like to evaluate different definitions of the exposure and/or outcome. We could include these by created extra exposure-outcome pairs, but that would mean that all defined analyses would be executed against these variations of the definitions, and this may not be what we want. Perhaps we would like to define just a single sensitivity analyses with a different outcome definition, in which case we could argue that the strategy of selecting the outcome becomes part of the analysis.

In such a case, we can define the multiple strategies using a list:

outcomeIds = list(narrowDefinition = 1,
                  broadDefinition = 2)

exposureOutcome <- createExposureOutcome(exposureId = 1124300,
                                         outcomeId = outcomeIds)

When we specify an analysis, we can then refer to one definition or the other:

sccsAnalysisA <- createSccsAnalysis(analysisId = 1,
                                    description = "Simplest model, using narrow def.",
                                    outcomeType = "narrowDefinition",
                                    getDbSccsDataArgs = getDbSccsDataArgs1,
                                    createSccsEraDataArgs = createEraDataArgs1,
                                    fitSccsModelArgs = fitSccsModelArgs)

sccsAnalysisB <- createSccsAnalysis(analysisId = 2,
                                    description = "Simplest model, using broad def.",
                                    outcomeType = "broadDefinition",
                                    getDbSccsDataArgs = getDbSccsDataArgs1,
                                    createSccsEraDataArgs = createEraDataArgs1,
                                    fitSccsModelArgs = fitSccsModelArgs)

sccsAnalysisList2 <- list(sccsAnalysisA, sccsAnalysisB)

In this example, the first analysis (analysisID = 1) will use cohort definition 1 as outcome, whilst the second analysis analysis (analysisID = 2) will use cohort definition 2 as outcome.

The same mechanism can be used to specifiy types for the exposureId.

Executing multiple analyses

We can now run the analyses against the hypotheses of interest using the runScsAnalyses()function. This function will run all specified analyses against all hypotheses of interest, meaning that the total number of outcome models is length(sccsAnalysisList) * length(exposureOutcomeList).

result <- runSccsAnalyses(connectionDetails = connectionDetails,
                          cdmDatabaseSchema = cdmDatabaseSchema,
                          oracleTempSchema = cdmDatabaseSchema,
                          exposureDatabaseSchema = cdmDatabaseSchema,
                          exposureTable = "drug_era",
                          outcomeDatabaseSchema = cohortDatabaseSchema,
                          outcomeTable = outcomeTable,
                          cdmVersion = cdmVersion,
                          outputFolder = outputFolder,
                          combineDataFetchAcrossOutcomes = TRUE,
                          exposureOutcomeList = exposureOutcomeList,
                          sccsAnalysisList = sccsAnalysisList,
                          getDbSccsDataThreads = 1,
                          createSccsEraDataThreads = 5,
                          fitSccsModelThreads = 3,
                          cvThreads = 10)

In the code above, we provide the arguments for connecting to the database, which schemas and tables to use, as well as the analyses and hypotheses of interest. The outputFolder specifies where the outcome models and intermediate files will be written. We also instruct SelfControlledCaseSeries to use multiple threads for various stages in the analyses, meaning these will be executed in parallel on multiple CPUs in the computer. Multithreading can significantly reduce execution time, but will require more system resources such as memory and temporary disk space.

Restarting

If for some reason the execution was interrupted, you can restart by re-issuing the runSccsAnalyses() command. Any intermediate and final products that have already been completed and written to disk will be skipped.

Retrieving the results

The result of the runSccsAnalyses() is a data frame with one row per exposure-outcome-analysis combination. It provides the file names of the intermediate and end-result files that were constructed. For example, we can retrieve the fitted model for the combination of our drug of interest, outcome, and first analysis:

sccsModelFile <- result$sccsModelFile[result$exposureId == 1124300 & 
                                      result$outcomeId == 1 &
                                      result$analysisId == 1]
sccsModel <- readRDS(file.path(outputFolder, sccsModelFile))
summary(sccsModel)
#> sccsModel object summary
#> 
#> Outcome ID: 1
#> 
#> Outcome count:
#>   Event count Case count
#> 1       60118      26323
#> 
#> Estimates:
#>                               Name    ID  Estimate  lower .95  upper .95
#>   Exposure of interest: Diclofenac  1000     1.277      1.234      1.321
#>    logRr  seLogRr
#>   0.2447   0.0174

Note that some of the file names will appear several times in the table. For example, all analysis share the same sccsData object.

We can create a summary of the results using summarizeSccsAnalyses():

analysisSum <- summarizeSccsAnalyses(result, outputFolder)
head(analysisSum)
#>   analysisId exposureId outcomeId caseCount eventCount rr(Exposure of interest)
#> 1          1    1124300         1     26323      60118                1.2772953
#> 2          1     705178         1      1116       2743                1.0715129
#> 3          1     705944         1      5420       9244                0.9511726
#> 4          1     710650         1       681       1454                0.9329491
#> 5          1     714785         1      6509      11199                1.0374991
#> 6          1     719174         1      8855      22823                1.0544470
#>   ci95lb(Exposure of interest) ci95ub(Exposure of interest)
#> 1                    1.2343060                     1.321447
#> 2                    0.9495034                     1.207337
#> 3                    0.8907972                     1.015086
#> 4                    0.7873739                     1.100697
#> 5                    0.9839866                     1.093620
#> 6                    1.0116936                     1.098813
#>   logRr(Exposure of interest) seLogRr(Exposure of interest) rr(Pre-exposure)
#> 1                  0.24474478                    0.01740293               NA
#> 2                  0.06907160                    0.06128518               NA
#> 3                 -0.05005978                    0.03331990               NA
#> 4                 -0.06940460                    0.08545961               NA
#> 5                  0.03681306                    0.02694863               NA
#> 6                  0.05301646                    0.02107298               NA
#>   ci95lb(Pre-exposure) ci95ub(Pre-exposure) logRr(Pre-exposure)
#> 1                   NA                   NA                  NA
#> 2                   NA                   NA                  NA
#> 3                   NA                   NA                  NA
#> 4                   NA                   NA                  NA
#> 5                   NA                   NA                  NA
#> 6                   NA                   NA                  NA
#>   seLogRr(Pre-exposure)
#> 1                    NA
#> 2                    NA
#> 3                    NA
#> 4                    NA
#> 5                    NA
#> 6                    NA

This tells us, per exposure-outcome-analysis combination, the estimated relative risk and 95% confidence interval, as well as the number of subjects (cases) and the number of events observed for those subjects.

Empirical calibration

Now that we have produced estimates for all outcomes including our negative controls, we can perform empirical calibration to estimate the bias of the various analyses included in our study. We will create the calibration effect plots for every analysis ID. In each plot, the blue dots represent our negative control exposures, and the yellow diamond represents our exposure of interest: diclofenac. An unbiased, well-calibrated analysis should have 95% of the negative controls between the dashed lines (ie. 95% should have p > .05).

install.packages("EmpiricalCalibration")
library(EmpiricalCalibration)

# Analysis 1: Simplest model
negCons <- analysisSum[analysisSum$analysisId == 1 & analysisSum$exposureId != 1124300, ]
ei <-  analysisSum[analysisSum$analysisId == 1 & analysisSum$exposureId == 1124300, ]
null <- fitNull(negCons$`logRr(Exposure of interest)`, 
                negCons$`seLogRr(Exposure of interest)`)
plotCalibrationEffect(logRrNegatives = negCons$`logRr(Exposure of interest)`, 
                      seLogRrNegatives = negCons$`seLogRr(Exposure of interest)`, 
                      logRrPositives = ei$`logRr(Exposure of interest)`, 
                      seLogRrPositives = ei$`seLogRr(Exposure of interest)`, 
                      null)

# Analysis 2: Including prophylactics
negCons <- analysisSum[analysisSum$analysisId == 2 & analysisSum$exposureId != 1124300, ]
ei <-  analysisSum[analysisSum$analysisId == 2 & analysisSum$exposureId == 1124300, ]
null <- fitNull(negCons$`logRr(Exposure of interest)`, 
                negCons$`seLogRr(Exposure of interest)`)
plotCalibrationEffect(logRrNegatives = negCons$`logRr(Exposure of interest)`, 
                      seLogRrNegatives = negCons$`seLogRr(Exposure of interest)`, 
                      logRrPositives = ei$`logRr(Exposure of interest)`, 
                      seLogRrPositives = ei$`seLogRr(Exposure of interest)`, 
                      null)

# Analysis 3: Including prophylactics, age, season, pre-exposure, and censoring
negCons <- analysisSum[analysisSum$analysisId == 3 & analysisSum$exposureId != 1124300, ]
ei <-  analysisSum[analysisSum$analysisId == 3 & analysisSum$exposureId == 1124300, ]
null <- fitNull(negCons$`logRr(Exposure of interest)`, 
                negCons$`seLogRr(Exposure of interest)`)
plotCalibrationEffect(logRrNegatives = negCons$`logRr(Exposure of interest)`, 
                      seLogRrNegatives = negCons$`seLogRr(Exposure of interest)`, 
                      logRrPositives = ei$`logRr(Exposure of interest)`, 
                      seLogRrPositives = ei$`seLogRr(Exposure of interest)`, 
                      null)

# Analysis 4: Including all other drugs (as well as prophylactics, age, season, pre-
# exposure, and censoring)
negCons <- analysisSum[analysisSum$analysisId == 4 & analysisSum$exposureId != 1124300, ]
ei <-  analysisSum[analysisSum$analysisId == 4 & analysisSum$exposureId == 1124300, ]
null <- fitNull(negCons$`logRr(Exposure of interest)`, 
                negCons$`seLogRr(Exposure of interest)`)
plotCalibrationEffect(logRrNegatives = negCons$`logRr(Exposure of interest)`, 
                      seLogRrNegatives = negCons$`seLogRr(Exposure of interest)`, 
                      logRrPositives = ei$`logRr(Exposure of interest)`, 
                      seLogRrPositives = ei$`seLogRr(Exposure of interest)`, 
                      null)

Acknowledgments

Considerable work has been dedicated to provide the SelfControlledCaseSeries package.

citation("SelfControlledCaseSeries")
#> 
#> To cite package 'SelfControlledCaseSeries' in publications use:
#> 
#>   Martijn Schuemie, Patrick Ryan, Trevor Shaddox and Marc Suchard
#>   (2020). SelfControlledCaseSeries: Self-Controlled Case Series. R
#>   package version 1.4.1.
#>   https://github.com/OHDSI/SelfControlledCaseSeries
#> 
#> A BibTeX entry for LaTeX users is
#> 
#>   @Manual{,
#>     title = {SelfControlledCaseSeries: Self-Controlled Case Series},
#>     author = {Martijn Schuemie and Patrick Ryan and Trevor Shaddox and Marc Suchard},
#>     year = {2020},
#>     note = {R package version 1.4.1},
#>     url = {https://github.com/OHDSI/SelfControlledCaseSeries},
#>   }

Further, SelfControlledCaseSeries makes extensive use of the Cyclops package.

citation("Cyclops")
#> 
#> To cite Cyclops in publications use:
#> 
#> Suchard MA, Simpson SE, Zorych I, Ryan P, Madigan D (2013). "Massive
#> parallelization of serial inference algorithms for complex generalized
#> linear models." _ACM Transactions on Modeling and Computer Simulation_,
#> *23*, 10. <URL: http://dl.acm.org/citation.cfm?id=2414791>.
#> 
#> A BibTeX entry for LaTeX users is
#> 
#>   @Article{,
#>     author = {M. A. Suchard and S. E. Simpson and I. Zorych and P. Ryan and D. Madigan},
#>     title = {Massive parallelization of serial inference algorithms for complex generalized linear models},
#>     journal = {ACM Transactions on Modeling and Computer Simulation},
#>     volume = {23},
#>     pages = {10},
#>     year = {2013},
#>     url = {http://dl.acm.org/citation.cfm?id=2414791},
#>   }

This work is supported in part through the National Science Foundation grant IIS 1251151.