Completely different approach to event-dep-exposure diagnostic

schuemie · schuemie · commit 8a6068d97b04 · 2024-10-08T13:47:46.000+02:00
diff --git a/R/Diagnostics.R b/R/Diagnostics.R
@@ -362,99 +362,6 @@ computeEventDependentObservationP <- function(sccsModel) {
 }
 
 
-computeExposureDaysToEvent <- function(studyPopulation, sccsData, exposureEraId) {
-  # This number of days before exposure start are assumed to be dealt with and are removed from
-  # both numerator (exposure days) and denominator (observation days):
-  preExposureDays <- 60 + 1
-
-  cases <- studyPopulation$cases |>
-    select("caseId", "startDay", "endDay")
-
-  # Keep only exposures that overlap with the observation periods of the study population (also
-  # truncate exposures to the observation period):
-  exposures <- sccsData$eras |>
-    filter(.data$eraId == exposureEraId & .data$eraType == "rx") |>
-    inner_join(cases,
-               by = join_by("caseId", "eraEndDay" >= "startDay", "eraStartDay" < "endDay"),
-               copy = TRUE) |>
-    collect() |>
-    mutate(eraStartDay = pmax(eraStartDay, startDay),
-           eraEndDay = pmin(eraEndDay, endDay))
-
-  if (nrow(exposures) == 0) {
-    warning("No exposures found with era ID ", exposureEraId)
-    return(NULL)
-  }
-  firstOutcomes <- studyPopulation$outcomes |>
-    group_by(.data$caseId) |>
-    filter(row_number(.data$outcomeDay) == 1)
-
-  # Merge overlapping exposures if needed:
-  # exposures <- exposures |>
-  #   arrange(caseId, eraStartDay) |>
-  #   group_by(caseId) |>
-  #   mutate(newGroup = cumsum(lag(eraEndDay, default = first(eraEndDay)) < eraStartDay)) |>
-  #   group_by(caseId, newGroup) |>
-  #   summarise(
-  #     eraStartDay = min(eraStartDay),
-  #     eraEndDay = max(eraEndDay),
-  #     .groups = 'drop'
-  #   ) |>
-  #   select(caseId, eraStartDay, eraEndDay)
-
-  # Ensure at least <preExposureDays> before each exposure start, by moving end day back:
-  truncatedExposures <- exposures |>
-    arrange(caseId, eraStartDay) |>
-    group_by(caseId) |>
-    mutate(
-      eraEndDay = ifelse(lead(eraStartDay, default = Inf) - eraEndDay < preExposureDays, lead(eraStartDay) - preExposureDays, eraEndDay)
-    ) |>
-    filter(eraEndDay > eraStartDay) |>
-    select(caseId, eraStartDay, eraEndDay)
-
-  exposureDeltas <- truncatedExposures |>
-    inner_join(firstOutcomes, by = join_by("caseId")) |>
-    mutate(deltaExposureStart = .data$eraStartDay - .data$outcomeDay,
-           deltaExposureEnd = .data$eraEndDay - .data$outcomeDay) |>
-    select("caseId", "deltaExposureStart", "deltaExposureEnd")
-
-  # Remove pre-exposure time from observation periods:
-  joined <- studyPopulation$cases |>
-    select(caseId, startDay, endDay) |>
-    # left_join(truncatedExposures, by = "caseId") |>
-    left_join(exposures |>
-                select("caseId", "eraStartDay", "eraEndDay"),
-              by = "caseId") |>
-    arrange(caseId, eraStartDay)
-  truncatedObservationPeriods <- joined |>
-    group_by(caseId) |>
-    mutate(
-      periodStart = lag(eraStartDay, default = first(startDay)),
-      periodEnd = pmin(eraStartDay - preExposureDays, endDay),
-      lastPeriodStart = eraStartDay,
-      lastPeriodEnd = endDay
-    ) |>
-    filter(periodStart <= periodEnd) |>
-    select(caseId, start = periodStart, end = periodEnd) |>
-    bind_rows(
-      joined |>
-        group_by(caseId) |>
-        slice(n()) |>
-        transmute(caseId, start = if_else(is.na(eraStartDay), startDay, eraStartDay), end = endDay) |>
-        filter(start <= end)
-    ) |>
-    arrange(caseId, start) |>
-    ungroup()
-
-  observationPeriodDeltas <- truncatedObservationPeriods |>
-    inner_join(firstOutcomes, by = join_by("caseId")) |>
-    mutate(deltaStart = .data$start - .data$outcomeDay,
-           deltaEnd = .data$end - .data$outcomeDay) |>
-    select("caseId", "deltaStart", "deltaEnd")
-
-  return(list(exposureDeltas = exposureDeltas, observationPeriodDeltas = observationPeriodDeltas))
-}
-
 #' Compute p for whether exposure probability changed following the outcome
 #'
 #' @param exposureEraId       The exposure to create the era data for. If not specified it is
@@ -463,12 +370,19 @@ computeExposureDaysToEvent <- function(studyPopulation, sccsData, exposureEraId)
 #' @template StudyPopulation
 #' @template SccsData
 #' @param bounds              Bounds for the null of no change in the exposure rate.
+#' @param ignoreExposureStarts Ignore exposure starts when computing the diagnostic. This makes the
+#'                             diagnostic robust against the outcome temporarily preventing exposure
+#'                             starting, which should be dealt with by the pre-exposure window.
 #'
 #' @return
 #' The p-value
 #'
 #' @export
-computeExposureChangeP <- function(sccsData, studyPopulation, exposureEraId = NULL, bounds = log(c(0.5, 2))) {
+computeExposureChangeP <- function(sccsData,
+                                   studyPopulation,
+                                   exposureEraId = NULL,
+                                   bounds = log(c(0.5, 2)),
+                                   ignoreExposureStarts = TRUE) {
   errorMessages <- checkmate::makeAssertCollection()
   checkmate::assertClass(sccsData, "SccsData", add = errorMessages)
   checkmate::assertList(studyPopulation, min.len = 1, add = errorMessages)
@@ -481,60 +395,131 @@ computeExposureChangeP <- function(sccsData, studyPopulation, exposureEraId = NU
       stop("No exposure ID specified, but multiple exposures found")
     }
   }
-  data <- computeExposureDaysToEvent(studyPopulation = studyPopulation,
-                                     sccsData = sccsData,
-                                     exposureEraId = exposureEraId)
-  if (is.null(data)) {
-    return(NA)
-  }
-  periods <- dplyr::tibble(afterOutcome = c(0,1),
-                           start = c(-30, 0),
-                           end = c(-1, 30))
-
-  exposure <- periods |>
-    cross_join(data$exposureDeltas) |>
-    mutate(daysExposure = pmax(0, pmin(end, deltaExposureEnd) - pmax(start, deltaExposureStart) + 1)) |>
-    group_by(caseId, afterOutcome) |>
-    summarise(daysExposure = sum(daysExposure), .groups = "drop") |>
-    select(caseId, afterOutcome, daysExposure)
-
-  observation <- periods |>
-    cross_join(data$observationPeriodDeltas) |>
-    mutate(daysObserved = pmax(0, pmin(end, deltaEnd) - pmax(start, deltaStart) + 1)) |>
-    group_by(caseId, afterOutcome) |>
-    summarise(daysObserved = sum(daysObserved), .groups = "drop") |>
-    select(caseId, afterOutcome, daysObserved)
-
-  casesWithExposure <- exposure |>
-    distinct(caseId) |>
-    pull()
-
-  poissonData <- observation |>
-    filter(caseId %in% casesWithExposure & daysObserved > 0) |>
-    left_join(exposure, by = join_by(caseId, afterOutcome)) |>
+  cases <- studyPopulation$cases |>
+    select("caseId", "startDay", "endDay")
+
+  # Keep only exposures that overlap with the observation periods of the study population (also
+  # truncate exposures to the observation period):
+  exposures <- sccsData$eras |>
+    filter(.data$eraId == exposureEraId & .data$eraType == "rx") |>
+    inner_join(cases,
+               by = join_by("caseId", "eraEndDay" >= "startDay", "eraStartDay" < "endDay"),
+               copy = TRUE) |>
+    collect() |>
+    mutate(eraStartDay = pmax(eraStartDay, startDay),
+           eraEndDay = pmin(eraEndDay, endDay))
+
+  exposures <- exposures |>
+    arrange(caseId, eraStartDay) |>
+    group_by(caseId) |>
+    mutate(newGroup = cumsum(lag(eraEndDay, default = first(eraEndDay)) < eraStartDay)) |>
+    group_by(caseId, newGroup) |>
+    summarise(
+      eraStartDay = min(eraStartDay),
+      eraEndDay = max(eraEndDay),
+      .groups = 'drop'
+    ) |>
+    select(caseId, eraStartDay, eraEndDay)
+
+  firstOutcomes <- studyPopulation$outcomes |>
+    group_by(.data$caseId) |>
+    filter(row_number(.data$outcomeDay) == 1)
+
+  # Compute exposure days before and after outcome, after removing exposures starting in the
+  # respective windows.
+  joined <- exposures |>
+    inner_join(firstOutcomes, by = join_by("caseId")) |>
+    mutate(deltaExposureStart = .data$eraStartDay - .data$outcomeDay,
+           deltaExposureEnd = .data$eraEndDay - .data$outcomeDay)
+
+  exposureBefore <- joined |>
+    filter(deltaExposureEnd >= -30 & deltaExposureStart <= -1) |>
+    filter(!ignoreExposureStarts | deltaExposureStart < -30 | deltaExposureStart > -1) |>
+    mutate(deltaExposureStart = pmax(deltaExposureStart, -30),
+           deltaExposureEnd = pmin(deltaExposureEnd, -1)) |>
+    group_by(caseId) |>
+    summarise(daysExposed = sum(deltaExposureEnd - deltaExposureStart + 1)) |>
+    select(caseId, daysExposed)
+
+  exposureAfter <- joined |>
+    filter(deltaExposureEnd >= 0 & deltaExposureStart <= 29) |>
+    filter(deltaExposureStart < 0 | deltaExposureStart > 29) |>
+    mutate(deltaExposureStart = pmax(deltaExposureStart, 0),
+           deltaExposureEnd = pmin(deltaExposureEnd, 29)) |>
+    group_by(caseId) |>
+    summarise(daysExposed = sum(deltaExposureEnd - deltaExposureStart + 1)) |>
+    select(caseId, daysExposed)
+
+  # Compute days observed
+  joined <- firstOutcomes |>
+    inner_join(studyPopulation$cases, by = join_by("caseId")) |>
+    mutate(deltaObservationStart = .data$startDay - .data$outcomeDay,
+           deltaObservationEnd = .data$endDay - .data$outcomeDay)
+
+  observationBefore <- joined |>
+    filter(deltaObservationEnd >= -30 & deltaObservationStart <= -1) |>
+    filter(deltaObservationStart < -30 | deltaObservationStart > -1) |>
+    mutate(deltaObservationStart = pmax(deltaObservationStart, -30),
+           deltaObservationEnd = pmin(deltaObservationEnd, -1)) |>
+    group_by(caseId) |>
+    summarise(daysObserved = sum(deltaObservationEnd - deltaObservationStart + 1)) |>
+    select(caseId, daysObserved)
+
+  observationAfter <- joined |>
+    filter(deltaObservationEnd >= 0 & deltaObservationStart <= 29) |>
+    filter(deltaObservationStart < 0 | deltaObservationStart > 29) |>
+    mutate(deltaObservationStart = pmax(deltaObservationStart, 0),
+           deltaObservationEnd = pmin(deltaObservationEnd, 29)) |>
+    group_by(caseId) |>
+    summarise(daysObserved = sum(deltaObservationEnd - deltaObservationStart + 1)) |>
+    select(caseId, daysObserved)
+
+  poissonData <- left_join(
+    bind_rows(
+      observationBefore |>
+        mutate(afterOutcome = 0),
+      observationAfter |>
+        mutate(afterOutcome = 1)
+    ),
+    bind_rows(
+      exposureBefore |>
+                mutate(afterOutcome = 0),
+              exposureAfter |>
+                mutate(afterOutcome = 1)
+    ),
+    by = join_by("caseId", "afterOutcome")
+  ) |>
+    filter(daysObserved > 0) |>
     mutate(
       rowId = row_number(),
-      covariateId = 1
+      covariateId = 1,
+      daysExposed = if_else(is.na(daysExposed), 0, daysExposed)
     ) |>
     select(
       "rowId",
       stratumId = "caseId",
       "covariateId",
       covariateValue = "afterOutcome",
       time = "daysObserved",
-      y = "daysExposure"
+      y = "daysExposed"
     )
 
+  casesWithExposure <- poissonData |>
+    filter(y > 0) |>
+    pull(stratumId)
+
   poissonData <- poissonData |>
-    filter((covariateValue == 0 & time == 30) | (covariateValue == 1 & time == 31))
+    filter(stratumId %in% casesWithExposure)
 
   cyclopsData <- Cyclops::convertToCyclopsData(outcomes = poissonData,
                                                covariates = poissonData,
                                                addIntercept = FALSE,
                                                modelType = "cpr",
                                                quiet = TRUE)
   fit <- Cyclops::fitCyclopsModel(cyclopsData)
-  fit$log_likelihood
+  if (fit$return_flag != "SUCCESS") {
+    return(NA)
+  }
   logRr <- coef(fit)
   if (logRr >= bounds[1] && logRr <= bounds[2]) {
     llNull <- fit$log_likelihood
diff --git a/extras/EndOfExposureSimulations.R b/extras/EndOfExposureSimulations.R
@@ -55,7 +55,7 @@ writeLines(sprintf("Number of simulation scenarios: %d", length(scenarios)))
 # Run simulations ----------------------------------------------------------------------------------
 folder <- "e:/SccsEdeSimulations100"
 
-scenario = scenarios[[30]]
+scenario = scenarios[[34]]
 scenario$censorType
 
 simulateOne <- function(seed, scenario) {
@@ -180,14 +180,16 @@ simulateOne <- function(seed, scenario) {
   idx2 <- which(estimates$covariateId == 1001)
   p <- computeExposureChangeP(sccsData, studyPop, 1)
   p
+  p2 <- computeExposureChangeP(sccsData, studyPop, 1, ignoreExposureStarts = FALSE)
   # plotExposureCentered(studyPop, sccsData, 1)
   # plotOutcomeCentered(studyPop, sccsData, 1)
 
   row <- tibble(logRr = estimates$logRr[idx1],
                 ci95Lb = exp(estimates$logLb95[idx1]),
                 ci95Ub = exp(estimates$logUb95[idx1]),
                 diagnosticEstimate = exp(estimates$logRr[idx2]),
-                diagnosticP = p)
+                diagnosticP = p,
+                diagnosticP2 = p2)
   return(row)
 }
 
@@ -213,11 +215,13 @@ for (i in seq_along(scenarios)) {
   metrics <- results |>
     mutate(coverage = ci95Lb < scenario$trueRr & ci95Ub > scenario$trueRr,
            diagnosticEstimate = log(diagnosticEstimate),
-           failDiagnostic = diagnosticP < 0.05) |>
+           failDiagnostic = diagnosticP < 0.05,
+           failDiagnostic2 = diagnosticP2 < 0.05) |>
     summarise(coverage = mean(coverage, na.rm = TRUE),
               bias = mean(logRr - log(scenario$trueRr), na.rm = TRUE),
               meanDiagnosticEstimate = exp(mean(diagnosticEstimate, na.rm = TRUE)),
-              fractionFailingDiagnostic = mean(failDiagnostic, na.rm = TRUE))
+              fractionFailingDiagnostic = mean(failDiagnostic, na.rm = TRUE),
+              fractionFailingDiagnostic2 = mean(failDiagnostic2, na.rm = TRUE))
   metrics
   row <- as_tibble(scenarioKey) |>
     bind_cols(metrics)
diff --git a/man/computeExposureChangeP.Rd b/man/computeExposureChangeP.Rd
