inst/doc/EmpiricalMaxSprtCalibrationVignette.R

In EmpiricalCalibration: Routines for Performing Empirical Calibration of Observational Study Estimates

## ---- echo = FALSE, message = FALSE, warning = FALSE--------------------------
library(EmpiricalCalibration)
allCvsAndLlrs <- readRDS("allCvsAndLlrs.rds")
set.seed(123)

doEval <- require("Cyclops", quietly = TRUE)
cv <- NA

## ----eval=doEval--------------------------------------------------------------
maxSprtSimulationData <- simulateMaxSprtData(
  n = 10000,
  pExposure = 0.5,
  backgroundHazard = 0.001,
  tar = 10,
  nullMu = 0.2,
  nullSigma = 0.2,
  maxT = 100,
  looks = 10,
  numberOfNegativeControls = 50,
  numberOfPositiveControls = 1,
  positiveControlEffectSize = 4
)
head(maxSprtSimulationData)

## ----warning=FALSE,eval=doEval------------------------------------------------
library(Cyclops)
library(survival)

dataOutcome51 <- maxSprtSimulationData[maxSprtSimulationData$outcomeId == 51, ]
dataOutcome51LookT50 <- dataOutcome51[dataOutcome51$lookTime == 50, ]

cyclopsData <- createCyclopsData(
  Surv(time, outcome) ~ exposure , 
  modelType = "cox", 
  data = dataOutcome51LookT50
)
fit <- fitCyclopsModel(cyclopsData)
coef(fit)

## ----eval=doEval--------------------------------------------------------------
# Maximum log likelihood:
fit$log_likelihood

## ----eval=doEval--------------------------------------------------------------
llNull <- getCyclopsProfileLogLikelihood(
  object = fit,
  parm = "exposureTRUE",
  x = 0
)$value
llNull

## ----eval=doEval--------------------------------------------------------------
if (fit$return_flag == "ILLCONDITIONED" || coef(fit) < 0) {
  llr <- 0
} else {
  llr <- fit$log_likelihood - llNull
}
llr

## ----eval=doEval--------------------------------------------------------------
outcomesPerLook <- aggregate(outcome ~ lookTime, dataOutcome51, sum)
# Need incremental outcomes per look:
outcomesPerLook <- outcomesPerLook$outcome[order(outcomesPerLook$lookTime)]
outcomesPerLook[2:10] <- outcomesPerLook[2:10] - outcomesPerLook[1:9]

exposedTime <- sum(dataOutcome51$time[dataOutcome51$exposure == TRUE & 
                                        dataOutcome51$lookTime == 100])
unexposedTime <- sum(dataOutcome51$time[dataOutcome51$exposure == FALSE & 
                                          dataOutcome51$lookTime == 100])
cv <- computeCvBinomial(
  groupSizes = outcomesPerLook,
  z = unexposedTime / exposedTime,
  minimumEvents = 1,
  alpha = 0.05
)
cv

## ----eval=doEval--------------------------------------------------------------
llr > cv

## ----eval=doEval--------------------------------------------------------------
llProfileOutcome51LookT50 <- getCyclopsProfileLogLikelihood(
  object = fit,
  parm = "exposureTRUE",
  bounds = log(c(0.1, 10))
)
head(llProfileOutcome51LookT50)

## ----eval=doEval--------------------------------------------------------------
library(ggplot2)
ggplot(llProfileOutcome51LookT50, aes(x = point, y = value)) +
  geom_line()

## ----eval=doEval--------------------------------------------------------------
negativeControlProfilesLookT50 <- list()
dataLookT50 <- maxSprtSimulationData[maxSprtSimulationData$lookTime == 50, ]
for (i in 1:50) {
  dataOutcomeIlookT50 <- dataLookT50[dataLookT50$outcomeId == i, ]
  cyclopsData <- createCyclopsData(
    Surv(time, outcome) ~ exposure , 
    modelType = "cox", 
    data = dataOutcomeIlookT50
  )
  fit <- fitCyclopsModel(cyclopsData)
  llProfile <- getCyclopsProfileLogLikelihood(
    object = fit,
    parm = "exposureTRUE",
    bounds = log(c(0.1, 10))
  )
  negativeControlProfilesLookT50[[i]] <- llProfile
}

## ----eval=doEval--------------------------------------------------------------
nullLookT50 <- fitNullNonNormalLl(negativeControlProfilesLookT50)
nullLookT50

## ----eval=doEval--------------------------------------------------------------
calibratedCv <- computeCvBinomial(
  groupSizes = outcomesPerLook,
  z = unexposedTime / exposedTime,
  minimumEvents = 1,
  alpha = 0.05,
  nullMean = nullLookT50[1],
  nullSd = nullLookT50[2]
)
calibratedCv

## ----eval=doEval--------------------------------------------------------------
llr > calibratedCv

## ----eval=FALSE---------------------------------------------------------------
#  allCvsAndLlrs <- data.frame()
#  for (t in unique(maxSprtSimulationData$lookTime)) {
#  
#    # Compute likelihood profiles and LLR for all negative controls:
#    negativeControlProfilesLookTt <- list()
#    llrsLookTt <- c()
#    dataLookTt <- maxSprtSimulationData[maxSprtSimulationData$lookTime == t, ]
#    for (i in 1:50) {
#      dataOutcomeIlookTt <- dataLookTt[dataLookTt$outcomeId == i, ]
#      cyclopsData <- createCyclopsData(Surv(time, outcome) ~ exposure,
#        modelType = "cox",
#        data = dataOutcomeIlookTt
#      )
#      fit <- fitCyclopsModel(cyclopsData)
#  
#      # likelihood profile:
#      llProfile <- getCyclopsProfileLogLikelihood(
#        object = fit,
#        parm = "exposureTRUE",
#        bounds = log(c(0.1, 10))
#      )
#      negativeControlProfilesLookTt[[i]] <- llProfile
#  
#      # LLR:
#      llNull <- getCyclopsProfileLogLikelihood(
#        object = fit,
#        parm = "exposureTRUE",
#        x = 0
#      )$value
#      if (fit$return_flag == "ILLCONDITIONED" || coef(fit) < 0) {
#        llr <- 0
#      } else {
#        llr <- fit$log_likelihood - llNull
#      }
#      llrsLookTt[i] <- llr
#    }
#  
#    # Fit null distribution:
#    nullLookTt <- fitNullNonNormalLl(negativeControlProfilesLookTt)
#  
#    # Compute calibrated and uncalibrated CV for all negative controls:
#    cvs <- c()
#    calibratedCvsLookT <- c()
#    for (i in 1:50) {
#      dataOutcomeI <- maxSprtSimulationData[maxSprtSimulationData$outcomeId == i, ]
#      outcomesPerLook <- aggregate(outcome ~ lookTime, dataOutcomeI, sum)
#      # Need incremental outcomes per look:
#      outcomesPerLook <- outcomesPerLook$outcome[order(outcomesPerLook$lookTime)]
#      outcomesPerLook[2:10] <- outcomesPerLook[2:10] - outcomesPerLook[1:9]
#  
#      exposedTime <- sum(dataOutcomeI$time[dataOutcomeI$exposure == TRUE &
#                                             dataOutcomeI$lookTime == 100])
#      unexposedTime <- sum(dataOutcomeI$time[dataOutcomeI$exposure == FALSE &
#                                               dataOutcomeI$lookTime == 100])
#  
#      # Note: uncalibrated CV will be same for every t, but computing in loop
#      # over t for clarity of code:
#      cv <- computeCvBinomial(
#        groupSizes = outcomesPerLook,
#        z = unexposedTime / exposedTime,
#        minimumEvents = 1,
#        alpha = 0.05
#      )
#      cvs[i] <- cv
#  
#      calibratedCv <- computeCvBinomial(
#        groupSizes = outcomesPerLook,
#        z = unexposedTime / exposedTime,
#        minimumEvents = 1,
#        alpha = 0.05,
#        nullMean = nullLookTt[1],
#        nullSd = nullLookTt[2]
#      )
#      calibratedCvsLookT[i] <- calibratedCv
#    }
#  
#    # Store in a data frame:
#    allCvsAndLlrs <- rbind(
#      allCvsAndLlrs,
#      data.frame(
#        t = t,
#        outcomeId = 1:50,
#        llr = llrsLookTt,
#        cv = cvs,
#        calibrateCv = calibratedCvsLookT
#      )
#    )
#  }

## ----eval=doEval--------------------------------------------------------------
signals <- c()
calibratedSignals <- c()
for (i in 1:50) {
  idx <- allCvsAndLlrs$outcomeId == i
  signals[i] <- any(allCvsAndLlrs$llr[idx] > allCvsAndLlrs$cv[idx])
  calibratedSignals[i] <-  any(allCvsAndLlrs$llr[idx] > allCvsAndLlrs$calibrateCv[idx])
}
# Type 1 error when not calibrated (should be 0.05):
mean(signals)

# Type 2 error when calibrated (should be 0.05):
mean(calibratedSignals)