Finished updating single study vignette

.github/workflows/R_CMD_check_Hades.yml

@@ -20,8 +20,6 @@ jobs:
       fail-fast: false
       matrix:
         config:
-          - {os: windows-latest, r: 'release'}
-          - {os: macOS-latest, r: 'release'}
           - {os: ubuntu-20.04, r: 'release', rspm: "https://packagemanager.rstudio.com/cran/__linux__/focal/latest"}
           - {os: windows-latest, r: '4.2.3'}
           - {os: macOS-latest, r: '4.2.3'}

NAMESPACE

@@ -41,6 +41,7 @@ export(getDbCohortMethodData)
 export(getDefaultCmTable1Specifications)
 export(getFileReference)
 export(getFollowUpDistribution)
+export(getGeneralizabilityTable)
 export(getInteractionResultsSummary)
 export(getOutcomeModel)
 export(getPsModel)

R/Balance.R

@@ -188,7 +188,7 @@ computeMeansPerGroup <- function(cohorts, cohortMethodData, covariateFilter) {
            -"overallSumWSqrTarget",
            -"overallSumWSqrComparator")
 
-    return(result)
+  return(result)
 }
 
 #' Compute covariate balance before and after PS adjustment
@@ -221,7 +221,8 @@ computeMeansPerGroup <- function(cohorts, cohortMethodData, covariateFilter) {
 #'
 #' @return
 #' Returns a tibble describing the covariate balance before and after PS adjustment,
-#' with one row per covariate and the following columns:
+#' with one row per covariate, with the same data as the `covariateRef` table in the `CohortMethodData` object,
+#' and the following additional columns:
 #'
 #' - beforeMatchingMeanTarget: The (weighted) mean value in the target before PS adjustment.
 #' - beforeMatchingMeanComparator: The (weighted) mean value in the comparator before PS adjustment.
@@ -397,6 +398,10 @@ computeCovariateBalance <- function(population,
     ) %>%
     arrange(desc(abs(.data$beforeMatchingStdDiff)))
 
+  metaData <- attr(population, "metaData")
+  if (!is.null(metaData) && !is.null(metaData$targetEstimator)) {
+    attr(balance, "targetEstimator") <- metaData$targetEstimator
+  }
   delta <- Sys.time() - start
   message(paste("Computing covariate balance took", signif(delta, 3), attr(delta, "units")))
   return(balance)
@@ -761,3 +766,102 @@ plotCovariatePrevalence <- function(balance,
   }
   return(plot)
 }
+
+#' Get information on generalizability
+#'
+#' @description
+#' to assess generalizability we compare the distribution of covariates before and after
+#' any (propensity score) adjustments. We compute the standardized difference of mean as
+#' our metric of generalizability. (Lipton et al., 2017)
+#'
+#' Depending on our target estimand, we need to consider a different base population for
+#' generalizability. For example, if we aim to estimate the average treatment effect in
+#' thetreated (ATT), our base population should be the target population, meaning we
+#' should consider the covariate distribution before and after PS adjustment in the target
+#' population only. By default this function will attempt to select the right base
+#' population based on what operations have been performed on the population. For example,
+#' if PS matching has been performed we assume the target estimand is the ATT, and the
+#' target population is selected as base.
+#'
+#' Requires running [computeCovariateBalance()]` first.
+#'
+#' @param balance       A data frame created by the `computeCovariateBalance` function.
+#' @param baseSelection The selection of the population to consider for generalizability.
+#'                      Options are "auto", "target", "comparator", and "both". The "auto"
+#'                      option will attempt to use the balance meta-data to pick the most
+#'                      appropriate population based on the target estimator.
+#'
+#' @return
+#' A tibble with the following columns:
+#'
+#' - covariateId: The ID of the covariate. Can be linked to the `covariates` and `covariateRef`
+#'   tables in the `CohortMethodData` object.
+#' - covariateName: The name of the covariate.
+#' - beforeMatchingMean: The mean covariate value before any (propensity score) adjustment.
+#' - afterMatchingMean: The mean covariate value after any (propensity score) adjustment.
+#' - stdDiff: The standardized difference of means between before and after adjustment.
+#'
+#' The tibble also has a 'baseSelection' attribute, documenting the base population used
+#' to assess generalizability.
+#'
+#' @references Tipton E, Hallberg K, Hedges LV, Chan W (2017) Implications of Small Samples
+#' for Generalization: Adjustments and Rules of Thumb, Eval Rev. Oct;41(5):472-505.
+#'
+#' @export
+getGeneralizabilityTable <- function(balance, baseSelection = "auto") {
+  errorMessages <- checkmate::makeAssertCollection()
+  checkmate::assertDataFrame(balance, add = errorMessages)
+  checkmate::assertCharacter(baseSelection, len = 1, add = errorMessages)
+  checkmate::assertChoice(baseSelection, c("auto", "target", "comparator",  "both"), add = errorMessages)
+  checkmate::reportAssertions(collection = errorMessages)
+
+  if (baseSelection == "auto") {
+    targetEstimator <- attr(balance, "targetEstimator")
+    if (is.null(targetEstimator)) {
+      stop("The baseSelection is set to 'auto' but the balance object does not contain a target estimator attribute. ",
+           "Please set the baseSelection manually.")
+    }
+    if (targetEstimator == "ate" | targetEstimator == "ato") {
+      baseSelection <- "both"
+      message("Selecting both target and comparator as base for generalizability")
+    } else if (targetEstimator == "att") {
+      baseSelection <- "target"
+      message("Selecting target as base for generalizability")
+    } else if (targetEstimator == "atu") {
+      baseSelection <- "comparator"
+      message("Selecting comparator as base for generalizability")
+    } else {
+      stop("Unkown target estimator: ", targetEstimator)
+    }
+  }
+  if (baseSelection == "target") {
+    generalizability <- balance %>%
+      mutate(absGeneralizabilityStdDiff = abs(.data$targetStdDiff)) %>%
+      arrange(desc(.data$absGeneralizabilityStdDiff)) %>%
+      select("covariateId",
+             "covariateName",
+             beforeMatchingMean = "beforeMatchingMeanTarget",
+             afterMatchingMean = "afterMatchingMeanTarget",
+             stdDiff = "targetStdDiff")
+  } else if (baseSelection == "comparator") {
+    generalizability <- balance %>%
+      mutate(absGeneralizabilityStdDiff = abs(.data$comparatorStdDiff)) %>%
+      arrange(desc(.data$absGeneralizabilityStdDiff)) %>%
+      select("covariateId",
+             "covariateName",
+             beforeMatchingMean = "beforeMatchingMeanComparator",
+             afterMatchingMean = "afterMatchingMeanComparator",
+             stdDiff = "comparatorStdDiff")
+  } else {
+    generalizability <- balance %>%
+      mutate(absGeneralizabilityStdDiff = abs(.data$targetComparatorStdDiff)) %>%
+      arrange(desc(.data$absGeneralizabilityStdDiff)) %>%
+      select("covariateId",
+             "covariateName",
+             "beforeMatchingMean",
+             "afterMatchingMean",
+             stdDiff = "targetComparatorStdDiff")
+  }
+  attr(generalizability, "baseSelection") <- baseSelection
+  return(generalizability)
+}

R/OutcomeModels.R

@@ -647,6 +647,7 @@ print.OutcomeModel <- function(x, ...) {
   writeLines(paste("Stratified:", x$outcomeModelStratified))
   writeLines(paste("Use covariates:", x$outcomeModelUseCovariates))
   writeLines(paste("Use inverse probability of treatment weighting:", x$inversePtWeighting))
+  writeLines(paste("Target estimand:", x$targetEstimator))
   writeLines(paste("Status:", x$outcomeModelStatus))
   if (!is.null(x$outcomeModelPriorVariance) && !is.na(x$outcomeModelPriorVariance)) {
     writeLines(paste("Prior variance:", x$outcomeModelPriorVariance))

extras/SingleStudyVignetteDataFetch.R

@@ -36,14 +36,14 @@ cohortDatabaseSchema <- "scratch_mschuemi"
 cohortTable  <- "cm_vignette"
 
 
-osteoArthritisOfKneeConceptId <- 4079750
-celecoxibConceptId <- 1118084
-diclofenacConceptId <- 1124300
+
 
 # Define exposure cohorts ------------------------------------------------------
 library(Capr)
-library(CirceR)
 
+osteoArthritisOfKneeConceptId <- 4079750
+celecoxibConceptId <- 1118084
+diclofenacConceptId <- 1124300
 osteoArthritisOfKnee <- cs(
   descendants(osteoArthritisOfKneeConceptId),
   name = "Osteoarthritis of knee"
@@ -81,22 +81,21 @@ diclofenacCohort <- cohort(
                                      persistenceWindow = 30,
                                      surveillanceWindow = 0))
 )
+# Define outcome cohort --------------------------------------------------------
+library(PhenotypeLibrary)
+outcomeCohorts <- getPlCohortDefinitionSet(77) # GI bleed
+
+# Generate cohorts -------------------------------------------------------------
+library(CirceR)
 exposureCohorts <- tibble(cohortId = c(1,2),
                           cohortName = c("Celecoxib", "Diclofenac"),
                           json = c(as.json(celecoxibCohort), as.json(diclofenacCohort)))
 exposureCohorts$sql <- sapply(exposureCohorts$json,
                               buildCohortQuery,
                               options = createGenerateOptions())
-
-# Define outcome cohort --------------------------------------------------------
-library(PhenotypeLibrary)
-outcomeCohorts <- getPlCohortDefinitionSet(77) # GI bleed
-
-
-# Generate cohorts -------------------------------------------------------------
-library(CohortGenerator)
 allCohorts <- bind_rows(outcomeCohorts,
                         exposureCohorts)
+library(CohortGenerator)
 cohortTableNames <- getCohortTableNames(cohortTable = cohortTable)
 createCohortTables(connectionDetails = connectionDetails,
                    cohortDatabaseSchema = cohortDatabaseSchema,
@@ -216,6 +215,11 @@ saveRDS(balance, file = file.path(folder, "balance.rds"))
 table1 <- createCmTable1(balance)
 print(table1, row.names = FALSE, right = FALSE)
 plotCovariateBalanceScatterPlot(balance, showCovariateCountLabel = TRUE, showMaxLabel = TRUE, fileName = "extras/balanceScatterplot.png")
+getGeneralizabilityTable(balance)
+
+balanceIptw <- computeCovariateBalance(ps, cohortMethodData)
+saveRDS(balanceIptw, file = file.path(folder, "balanceIptw.rds"))
+getGeneralizabilityTable(balanceIptw)
 
 outcomeModel <- fitOutcomeModel(
   population = studyPop,
@@ -230,6 +234,15 @@ coef(outcomeModel)
 confint(outcomeModel)
 saveRDS(outcomeModel, file = file.path(folder, "OutcomeModel1.rds"))
 
+outcomeModel <- fitOutcomeModel(
+  population = ps,
+  modelType = "cox",
+  stratified = FALSE,
+  useCovariates = FALSE,
+  inversePtWeighting = TRUE
+)
+saveRDS(outcomeModel, file = file.path(folder, "OutcomeModel1b.rds"))
+
 outcomeModel <- fitOutcomeModel(
   population = matchedPop,
   modelType = "cox",