From bacd803833391903a1781333c89063b91bbc22bc Mon Sep 17 00:00:00 2001 From: schuemie Date: Mon, 30 Sep 2024 12:03:39 +0200 Subject: [PATCH] Starting work on simulating bias caused by event-dependent observation --- extras/TestEndOfObsDiagnostic.R | 95 +++++++++++++++++++++++++++++++++ src/SccsSimulator.cpp | 7 ++- 2 files changed, 100 insertions(+), 2 deletions(-) diff --git a/extras/TestEndOfObsDiagnostic.R b/extras/TestEndOfObsDiagnostic.R index 0be1d8c..39b01aa 100644 --- a/extras/TestEndOfObsDiagnostic.R +++ b/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 + diff --git a/src/SccsSimulator.cpp b/src/SccsSimulator.cpp index 40e049a..dabc8ac 100644 --- a/src/SccsSimulator.cpp +++ b/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); }