Starting work on simulating bias caused by event-dependent observation

extras/TestEndOfObsDiagnostic.R

@@ -318,3 +318,98 @@ ggplot(data, aes(x = startDate + 15, y = ir)) +
   geom_line() +
   geom_point()
 
+# Simulations --------------------------------------------------------------------------------------
+library(SelfControlledCaseSeries)
+
+rw <- createSimulationRiskWindow(start = 0,
+                                 end = 0,
+                                 endAnchor = "era end",
+                                 relativeRisks = 2)
+
+settings <- createSccsSimulationSettings(minBaselineRate = 0.001,
+                                         maxBaselineRate = 0.01,
+                                         eraIds = 1,
+                                         simulationRiskWindows = list(rw),
+                                         includeAgeEffect = FALSE,
+                                         includeSeasonality = FALSE,
+                                         includeCalendarTimeEffect = FALSE)
+
+simulateOne <- function(seed, settings, censorFunction) {
+  set.seed(seed)
+  sccsData <- simulateSccsData(1000, settings)
+  sccsData$cases |>
+    count()
+  firstOutcomeEra <- sccsData$eras |>
+    filter(eraId == 10) |>
+    group_by(caseId) |>
+    filter(row_number(eraStartDay) == 1) |>
+    ungroup() |>
+    select(caseId, outcomeDay = eraStartDay)
+
+  # Censor observation at outcome:
+  sccsData$cases <- sccsData$cases |>
+    inner_join(firstOutcomeEra, by = join_by(caseId)) |>
+    collect() |>
+    mutate(endDay = censorFunction(endDay, outcomeDay)) |>
+    select(-outcomeDay)
+
+
+  covarSettings <- createEraCovariateSettings(label = "Exposures of interest",
+                                              includeEraIds = 1,
+                                              stratifyById = FALSE,
+                                              start = 0,
+                                              end = 0,
+                                              endAnchor = "era end")
+
+  studyPop <- createStudyPopulation(sccsData = sccsData,
+                                    outcomeId = settings$outcomeId,
+                                    firstOutcomeOnly = TRUE,
+                                    naivePeriod = 0)
+
+  sccsIntervalData <- createSccsIntervalData(studyPopulation = studyPop,
+                                             sccsData = sccsData,
+                                             eraCovariateSettings = covarSettings)
+
+  model <- fitSccsModel(sccsIntervalData)
+  estimates <- model$estimates
+  return(tibble(ci95Lb = exp(estimates$logLb95[2]),
+                ci95Ub = exp(estimates$logUb95[2]),
+                diagnosticEstimate = exp(estimates$logRr[1])))
+}
+
+# One change at dying in next week:
+censorFunction <- function(endDay, outcomeDay) {
+  if_else(runif(length(endDay)) < 0.9, round(pmin(endDay, outcomeDay + runif(length(endDay), 0, 7))), endDay)
+}
+# Added hazard of dying for rest of time:
+censorFunction <- function(endDay, outcomeDay) {
+  pmin(endDay, outcomeDay + rexp(length(endDay), 0.05))
+}
+
+cluster <- ParallelLogger::makeCluster(10)
+ParallelLogger::clusterRequire(cluster, "SelfControlledCaseSeries")
+
+
+
+results <- ParallelLogger::clusterApply(cluster, 1:100, simulateOne, settings = settings, censorFunction = censorFunction)
+results <- bind_rows(results)
+coverage <- results |>
+  mutate(coverage = ci95Lb < rw$relativeRisks & ci95Ub > rw$relativeRisks) |>
+  summarise(mean(coverage)) |>
+  pull()
+diagnosticEstimate <- results |>
+  mutate(diagnosticEstimate = log(diagnosticEstimate)) |>
+  summarise(exp(mean(diagnosticEstimate))) |>
+  pull()
+writeLines(sprintf("Coverage: %0.3f, mean diagnostic estimate: %0.2f", coverage, diagnosticEstimate))
+
+
+
+
+ParallelLogger::stopCluster(cluster)
+
+print(sprintf("True IRR = %0.2f", rw$relativeRisks))
+
+h <- 0.0003
+1-(1-h)^365
+

src/SccsSimulator.cpp

@@ -110,11 +110,14 @@ List SccsSimulator::simulateOutcomes(){
   int eraEndIndex = 0;
   for (int caseIndex = 0; caseIndex < casesCaseId.size(); caseIndex++){
     int caseId = casesCaseId[caseIndex];
-    while (erasCaseId[eraEndIndex] == caseId){
+    while (eraEndIndex < erasCaseId.size() && erasCaseId[eraEndIndex] < caseId){
       eraEndIndex++;
     }
-    processPerson(caseIndex, eraStartIndex, eraEndIndex);
     eraStartIndex = eraEndIndex;
+    while (eraEndIndex < erasCaseId.size() && erasCaseId[eraEndIndex] == caseId){
+      eraEndIndex++;
+    }
+    processPerson(caseIndex, eraStartIndex, eraEndIndex);
   }
   return List::create(Named("caseId") = outcomeCaseIds, Named("startDay") = outcomeStartDays);
 }