Add diagnostic for reverse causality (pre-exposure increase in risk)

[schuemie]

The SCCS can detect an increase in risk after exposure initiation, but if that increase in risk was already present (and even higher) right before the exposure this means that likely the exposure did not cause the increased risk. For example, if the exposure is indicated for a disease, and the disease causes the increased risk, we may see this pattern. We therefore need a diagnostic that evaluates whether the risk of the outcome is greater just before exposure compared to just after exposure initiation.

Ideally, this diagnostic would reverse if a decrease in risk was observed during exposure instead of an increase. Perhaps this behavior can be manually switched by the user.

Below is messy code from a private study package that computes a p-value for H1: risk before > risk after against H0: risk before =< risk after.

computePreExposureGainP <- function(sccsData, studyPopulation, exposureEraId) {
  # Copied from plotExposureCentered:
  cases <- studyPopulation$cases %>%
    select(.data$caseId, caseEndDay = .data$endDay, .data$offset)

  exposures <- sccsData$eras %>%
    filter(.data$eraId == exposureEraId & .data$eraType == "rx") %>%
    group_by(.data$caseId) %>%
    inner_join(cases, by = "caseId", copy = TRUE) %>%
    mutate(
      startDay = .data$startDay - .data$offset,
      endDay = .data$endDay - .data$offset
    ) %>%
    filter(.data$startDay >= 0, .data$startDay < .data$caseEndDay) %>%
    collect()

  if (nrow(exposures) == 0) {
    warning("No exposures found with era ID ", exposureEraId)
    return(NA)
  }
  firstExposures <- exposures %>%
    group_by(.data$caseId, .data$caseEndDay) %>%
    summarise(
      startDay = min(.data$startDay, na.rm = TRUE),
      endDay = min(.data$endDay, na.rm = TRUE),
      .groups = "drop_last"
    )

  outcomes <- studyPopulation$outcomes %>%
    inner_join(firstExposures, by = "caseId") %>%
    mutate(delta = .data$outcomeDay - .data$startDay) %>%
    select(.data$caseId, .data$outcomeDay, .data$delta)

  # end copy

  # Restrict to 30 days before and after exposure start:
  outcomes <- outcomes %>%
    filter(.data$delta >= -30 & .data$delta <= 30) %>%
    mutate(
      beforeExposure = .data$delta < 0,
      y = 1
    ) %>%
    select(
      .data$caseId,
      .data$beforeExposure,
      .data$y
    )

  observed <- bind_rows(
    firstExposures %>%
      mutate(
        daysObserved = if_else(.data$startDay > 30, 30, .data$startDay),
        beforeExposure = TRUE
      ) %>%
      select(
        .data$caseId,
        .data$daysObserved,
        .data$beforeExposure
      ),
    firstExposures %>%
      mutate(daysAfterExposure = .data$caseEndDay - .data$startDay) %>%
      mutate(
        daysObserved = if_else(.data$daysAfterExposure > 30, 30, .data$daysAfterExposure),
        beforeExposure = FALSE
      ) %>%
      select(
        .data$caseId,
        .data$daysObserved,
        .data$beforeExposure
      )
  ) %>%
    filter(.data$daysObserved > 0)

  poissonData <- observed %>%
    left_join(outcomes, by = c("caseId", "beforeExposure")) %>%
    mutate(
      y = if_else(is.na(.data$y), 0, .data$y),
      logDays = log(.data$daysObserved)
    )
  # library(survival)
  cyclopsData <- Cyclops::createCyclopsData(y ~ beforeExposure + strata(caseId) + offset(logDays),
                                            modelType = "cpr",
                                            data = poissonData
  )
  fit <- Cyclops::fitCyclopsModel(cyclopsData)
  if (fit$return_flag == "ILLCONDITIONED") {
    return(NA)
  }
  # compute one-sided p-value:
  llNull <- Cyclops::getCyclopsProfileLogLikelihood(
    object = fit,
    parm = "beforeExposureTRUE",
    x = 0
  )$value
  llr <- fit$log_likelihood - llNull
  p <- EmpiricalCalibration:::computePFromLlr(llr, coef(fit))
  return(p)
}