diff --git a/analysis/0-lib/campaign_info.json b/analysis/0-lib/campaign_info.json deleted file mode 100644 index fbf95dd..0000000 --- a/analysis/0-lib/campaign_info.json +++ /dev/null @@ -1,158 +0,0 @@ -{ - "1900-01-01": { - "campaign_label": "Pre-2020-07-01", - "campaign_start_date": "1900-01-01", - "primary_milestone_date": "1900-01-01", - "age_date": "1900-01-01", - "age_threshold": 16, - "clinical_priority": "primis_atrisk", - "early_milestone_date": "1900-02-25", - "final_milestone_date": "2020-06-30", - "early_milestone_days": 56, - "primary_milestone_days": 1, - "final_milestone_days": 44011 - }, - "2020-07-01": { - "campaign_label": "Pre-roll-out", - "campaign_start_date": "2020-07-01", - "primary_milestone_date": "2020-07-01", - "age_date": "2020-07-01", - "age_threshold": 16, - "clinical_priority": "primis_atrisk", - "early_milestone_date": "2020-08-25", - "final_milestone_date": "2020-12-06", - "early_milestone_days": 56, - "primary_milestone_days": 1, - "final_milestone_days": 159 - }, - "2020-12-07": { - "campaign_label": "Primary series", - "campaign_start_date": "2020-12-07", - "primary_milestone_date": "2021-06-30", - "age_date": "2021-03-31", - "age_threshold": 16, - "clinical_priority": "primis_atrisk", - "early_milestone_date": "2021-01-31", - "final_milestone_date": "2021-09-05", - "early_milestone_days": 56, - "primary_milestone_days": 206, - "final_milestone_days": 273 - }, - "2021-09-06": { - "campaign_label": "Autumn 2021", - "campaign_start_date": "2021-09-06", - "primary_milestone_date": "2022-02-28", - "age_date": "2021-08-31", - "age_threshold": 16, - "clinical_priority": "primis_atrisk", - "early_milestone_date": "2021-10-31", - "final_milestone_date": "2022-03-20", - "early_milestone_days": 56, - "primary_milestone_days": 176, - "final_milestone_days": 196 - }, - "2022-03-21": { - "campaign_label": "Spring 2022", - "campaign_start_date": "2022-03-21", - "primary_milestone_date": "2022-06-30", - "age_date": "2022-06-30", - "age_threshold": 75, - "clinical_priority": "immunosuppressed", - "early_milestone_date": "2022-05-15", - "final_milestone_date": "2022-08-28", - "early_milestone_days": 56, - "primary_milestone_days": 102, - "final_milestone_days": 161 - }, - "2022-08-29": { - "campaign_label": "Autumn 2022", - "campaign_start_date": "2022-08-29", - "primary_milestone_date": "2023-02-28", - "age_date": "2023-03-31", - "age_threshold": 50, - "clinical_priority": "primis_atrisk", - "early_milestone_date": "2022-10-23", - "final_milestone_date": "2023-04-02", - "early_milestone_days": 56, - "primary_milestone_days": 184, - "final_milestone_days": 217 - }, - "2023-04-03": { - "campaign_label": "Spring 2023", - "campaign_start_date": "2023-04-03", - "primary_milestone_date": "2023-06-30", - "age_date": "2023-06-30", - "age_threshold": 75, - "clinical_priority": "immunosuppressed", - "early_milestone_date": "2023-05-28", - "final_milestone_date": "2023-08-27", - "early_milestone_days": 56, - "primary_milestone_days": 89, - "final_milestone_days": 147 - }, - "2023-08-28": { - "campaign_label": "Autumn 2023", - "campaign_start_date": "2023-08-28", - "primary_milestone_date": "2024-02-28", - "age_date": "2024-03-31", - "age_threshold": 65, - "clinical_priority": "primis_atrisk", - "early_milestone_date": "2023-10-22", - "final_milestone_date": "2024-04-14", - "early_milestone_days": 56, - "primary_milestone_days": 185, - "final_milestone_days": 231 - }, - "2024-04-15": { - "campaign_label": "Spring 2024", - "campaign_start_date": "2024-04-15", - "primary_milestone_date": "2024-06-30", - "age_date": "2024-06-30", - "age_threshold": 75, - "clinical_priority": "immunosuppressed", - "early_milestone_date": "2024-06-09", - "final_milestone_date": "2024-09-29", - "early_milestone_days": 56, - "primary_milestone_days": 77, - "final_milestone_days": 168 - }, - "2024-09-30": { - "campaign_label": "Autumn 2024", - "campaign_start_date": "2024-09-30", - "primary_milestone_date": "2025-02-28", - "age_date": "2025-03-31", - "age_threshold": 65, - "clinical_priority": "primis_atrisk", - "early_milestone_date": "2024-11-24", - "final_milestone_date": "2025-03-30", - "early_milestone_days": 56, - "primary_milestone_days": 152, - "final_milestone_days": 182 - }, - "2025-03-31": { - "campaign_label": "Spring 2025", - "campaign_start_date": "2025-03-31", - "primary_milestone_date": "2025-06-30", - "age_date": "2025-06-30", - "age_threshold": 75, - "clinical_priority": "immunosuppressed", - "early_milestone_date": "2025-05-25", - "final_milestone_date": "2025-09-28", - "early_milestone_days": 56, - "primary_milestone_days": 92, - "final_milestone_days": 182 - }, - "2025-09-29": { - "campaign_label": "Autumn 2025", - "campaign_start_date": "2025-09-29", - "primary_milestone_date": "2026-02-28", - "age_date": "2026-03-31", - "age_threshold": 75, - "clinical_priority": "primis_atrisk", - "early_milestone_date": "2025-11-23", - "final_milestone_date": "2026-01-31", - "early_milestone_days": 56, - "primary_milestone_days": 153, - "final_milestone_days": 125 - } -} diff --git a/analysis/0-lib/cohort_info.json b/analysis/0-lib/cohort_info.json new file mode 100644 index 0000000..f4f613e --- /dev/null +++ b/analysis/0-lib/cohort_info.json @@ -0,0 +1,102 @@ +{ + "covid_2025_26": { + "cohort_id": "covid_2025_26", + "target": "COVID-19", + "cohort": "2025/26", + "cohort_start_date": "2025-09-01", + "cohort_end_date": "2026-03-31", + "age_threshold": 75, + "clinical_priority": "primis_atrisk", + "cohort_end_days": 212 + }, + "covid_ba1_ba2": { + "cohort_id": "covid_ba1_ba2", + "target": "COVID-19", + "cohort": "BA.1/BA.2 wave", + "cohort_start_date": "2021-12-15", + "cohort_end_date": "2022-04-29", + "age_threshold": "NA", + "clinical_priority": "primis_atrisk", + "cohort_end_days": 136 + }, + "covid_jn1": { + "cohort_id": "covid_jn1", + "target": "COVID-19", + "cohort": "JN.1 wave", + "cohort_start_date": "2023-12-04", + "cohort_end_date": "2024-03-31", + "age_threshold": "NA", + "clinical_priority": "primis_atrisk", + "cohort_end_days": 119 + }, + "covid_wave_1": { + "cohort_id": "covid_wave_1", + "target": "COVID-19", + "cohort": "Wave 1", + "cohort_start_date": "2020-03-23", + "cohort_end_date": "2020-05-30", + "age_threshold": "NA", + "clinical_priority": "primis_atrisk", + "cohort_end_days": 69 + }, + "flu_2023_24": { + "cohort_id": "flu_2023_24", + "target": "Influenza", + "cohort": "2023/24", + "cohort_start_date": "2023-09-01", + "cohort_end_date": "2024-03-31", + "age_threshold": 65, + "clinical_priority": "primis_atrisk", + "cohort_end_days": 213 + }, + "flu_2024_25": { + "cohort_id": "flu_2024_25", + "target": "Influenza", + "cohort": "2024/25", + "cohort_start_date": "2024-09-01", + "cohort_end_date": "2025-03-31", + "age_threshold": 65, + "clinical_priority": "primis_atrisk", + "cohort_end_days": 212 + }, + "flu_2025_26": { + "cohort_id": "flu_2025_26", + "target": "Influenza", + "cohort": "2025/26", + "cohort_start_date": "2025-09-01", + "cohort_end_date": "2026-03-31", + "age_threshold": 65, + "clinical_priority": "primis_atrisk", + "cohort_end_days": 212 + }, + "rsv_2023_24": { + "cohort_id": "rsv_2023_24", + "target": "RSV", + "cohort": "2023/24", + "cohort_start_date": "2023-09-01", + "cohort_end_date": "2024-03-31", + "age_threshold": "NA", + "clinical_priority": "primis_atrisk", + "cohort_end_days": 213 + }, + "rsv_2024_25": { + "cohort_id": "rsv_2024_25", + "target": "RSV", + "cohort": "2024/25", + "cohort_start_date": "2024-09-01", + "cohort_end_date": "2025-03-31", + "age_threshold": "NA", + "clinical_priority": "primis_atrisk", + "cohort_end_days": 212 + }, + "rsv_2025_26": { + "cohort_id": "rsv_2025_26", + "target": "RSV", + "cohort": "2025/26", + "cohort_start_date": "2025-09-01", + "cohort_end_date": "2026-03-31", + "age_threshold": 75, + "clinical_priority": "primis_atrisk", + "cohort_end_days": 212 + } +} diff --git a/analysis/0-lib/design.R b/analysis/0-lib/design.R index 973cb49..211e2fa 100644 --- a/analysis/0-lib/design.R +++ b/analysis/0-lib/design.R @@ -84,13 +84,15 @@ fct_case_when <- function(...) { # - start_date is when we start the observational period proper, at the start of the mass vax programme # - end_date is when we stop the observation period. This must be extended as the study progresses -study_dates <- - list( - firstpossiblevax_date = "2020-07-01", - start_date = "2020-12-07", - end_date = "2026-03-31" - ) |> - lapply(as.Date) +# study_dates <- +# list( +# # firstpossiblevax_date = "2020-07-01", +# start_date1 = "2010-09-01", +# end_date1 = "2013-08-31", +# start_date2 = "2020-03-23", +# end_date2 = "2026-08-31" +# ) |> +# lapply(as.Date) # make these available in the global environment # so we don't have to use `study_dates$start_date` or `start_date <- study_dates$start_date` in each script @@ -99,32 +101,39 @@ study_dates <- # statistical disclosure control rounding precision sdc_threshold <- 10 -# covid-19 vaccine campaign dates -campaign_info <- +# cohort dates +cohort_info <- tribble( - ~campaign_label, ~campaign_start_date, ~primary_milestone_date, ~age_date, ~age_threshold, ~clinical_priority, - "Pre-2020-07-01", "1900-01-01", "1900-01-01", "1900-01-01", 16, "primis_atrisk", - "Pre-roll-out", as.character(study_dates$firstpossiblevax_date), as.character(study_dates$firstpossiblevax_date), as.character(study_dates$firstpossiblevax_date), 16, "primis_atrisk", - "Primary series", "2020-12-07", "2021-06-30", "2021-03-31", 16, "primis_atrisk", - "Autumn 2021", "2021-09-06", "2022-02-28", "2021-08-31", 16, "primis_atrisk", - "Spring 2022", "2022-03-21", "2022-06-30", "2022-06-30", 75, "immunosuppressed", - "Autumn 2022", "2022-08-29", "2023-02-28", "2023-03-31", 50, "primis_atrisk", - "Spring 2023", "2023-04-03", "2023-06-30", "2023-06-30", 75, "immunosuppressed", - "Autumn 2023", "2023-08-28", "2024-02-28", "2024-03-31", 65, "primis_atrisk", - "Spring 2024", "2024-04-15", "2024-06-30", "2024-06-30", 75, "immunosuppressed", - "Autumn 2024", "2024-09-30", "2025-02-28", "2025-03-31", 65, "primis_atrisk", - "Spring 2025", "2025-03-31", "2025-06-30", "2025-06-30", 75, "immunosuppressed", - "Autumn 2025", "2025-09-29", "2026-02-28", "2026-03-31", 75, "primis_atrisk", - ) |> + ~cohort_id, ~target, ~cohort, ~cohort_start_date, ~cohort_end_date, ~age_threshold, ~clinical_priority, + + "flu_2023_24", "Influenza", "2023/24", "2023-09-01", "2024-03-31", 65, "primis_atrisk", + "flu_2024_25", "Influenza", "2024/25", "2024-09-01", "2025-03-31", 65, "primis_atrisk", + "flu_2025_26", "Influenza", "2025/26", "2025-09-01", "2026-03-31", 65, "primis_atrisk", + + "covid_wave_1", "COVID-19", "Wave 1", "2020-03-23", "2020-05-30", NA, "primis_atrisk", + "covid_ba1_ba2", "COVID-19", "BA.1/BA.2 wave", "2021-12-15", "2022-04-29", NA, "primis_atrisk", + "covid_jn1", "COVID-19", "JN.1 wave", "2023-12-04", "2024-03-31", NA, "primis_atrisk", + "covid_2025_26", "COVID-19", "2025/26", "2025-09-01", "2026-03-31", 75, "primis_atrisk", + + "rsv_2023_24", "RSV", "2023/24", "2023-09-01", "2024-03-31", NA, "primis_atrisk", + "rsv_2024_25", "RSV", "2024/25", "2024-09-01", "2025-03-31", NA, "primis_atrisk", + "rsv_2025_26", "RSV", "2025/26", "2025-09-01", "2026-03-31", 75, "primis_atrisk" #, + + # "zoster_2010_11", "Zoster", "2010/11", "2010-09-01", "2011-08-31", 70, "primis_atrisk", + # "zoster_2011_12", "Zoster", "2011/12", "2011-09-01", "2012-08-31", 70, "primis_atrisk", + # "zoster_2012_13", "Zoster", "2012/13", "2012-09-01", "2013-08-31", 70, "primis_atrisk", + # "zoster_2023_24", "Zoster", "2023/24", "2023-09-01", "2024-08-31", 65, "primis_atrisk", + # "zoster_2024_25", "Zoster", "2024/25", "2024-09-01", "2025-08-31", 65, "primis_atrisk", + # "zoster_2025_26", "Zoster", "2025/26", "2025-09-01", "2026-08-31", 65, "primis_atrisk" + )|> mutate( - across(c(campaign_start_date, primary_milestone_date, age_date), as.Date), - early_milestone_date = campaign_start_date + (7 * 8) - 1, # end of eighth week after campaign_start_date - final_milestone_date = lead(campaign_start_date, 1, as.Date("2026-02-01")) - 1 # day before next campaign date (or some arbitrary future date if last campaign) - ) |> + across( + c(cohort_start_date, cohort_end_date), + as.Date + ) + ) |> mutate( - early_milestone_days = as.integer(early_milestone_date - campaign_start_date) + 1L, - primary_milestone_days = as.integer(primary_milestone_date - campaign_start_date) + 1L, - final_milestone_days = as.integer(final_milestone_date - campaign_start_date) + 1L + cohort_end_days = as.integer(cohort_end_date - cohort_start_date) + 1L ) # output from https://jobs.opensafely.org/opensafely-internal/tpp-vaccination-names/ workspace @@ -353,34 +362,6 @@ standardise_demographic_characteristics <- ) ## CKD-RRT classification -# adapted from https://github.com/opensafely/covid_mortality_over_time/blob/main/analysis/utils/kidney_functions.R - -# egfr <- function(rrt_cat, sex, creatinine_umol, creatinine_age) { -# # creatinine mg/dL from µmol/L -# SCR_adj <- creatinine_umol / 88.4 - -# # min/max creatinine -# min_creat <- if_else( -# sex == "male", -# pmin(SCR_adj / 0.9, 1)^-0.411, -# pmin(SCR_adj / 0.7, 1)^-0.329 -# ) -# max_creat <- if_else( -# sex == "male", -# pmax(SCR_adj / 0.9, 1)^-1.209, -# pmax(SCR_adj / 0.7, 1)^-1.209 -# ) -# # eGFR (male/female) -# egfr_male <- case_when( -# is.na(creatinine_umol) ~ NA_real_, -# is.na(creatinine_age) ~ NA_real_, -# TRUE ~ (min_creat * max_creat * 141) * (0.993^creatinine_age) -# ) -# egfr <- if_else(sex == "female", 1.018 * egfr_male, egfr_male) - -# return(egfr) -# } - standardise_primis_and_extended_characteristics <- rlang::quos( @@ -425,30 +406,30 @@ level1_group <- c( "all", # level 1B (age4) - "ageband4", - - # level 1C (age13) - "ageband13", - - # Level 1D (eligibility) - "any_eligibility", - "age_above_eligiblity_threshold", - "clinical_priority", - "clinical_priority_only", - "carehome_status", - - # Level 1E (clinical risk) - "primis_atrisk", - "crd", - "chd", - "ckd", - "cld", - "cns_learndis", - "diabetes", - "immunosuppressed", - "asplenia", - "severe_obesity", - "smi" + "ageband4" #, + + # # level 1C (age13) + # "ageband13", + + # # Level 1D (eligibility) + # "any_eligibility", + # "age_above_eligiblity_threshold", + # "clinical_priority", + # "clinical_priority_only", + # "carehome_status", + + # # Level 1E (clinical risk) + # "primis_atrisk", + # "crd", + # "chd", + # "ckd", + # "cld", + # "cns_learndis", + # "diabetes", + # "immunosuppressed", + # "asplenia", + # "severe_obesity", + # "smi" ) level2_group <- c( @@ -459,7 +440,7 @@ level2_group <- c( "ageband13", "sex", "ethnicity5", - # "ethnicity16", + "ethnicity16", "region", "imd_quintile", "carehome_status", @@ -508,15 +489,16 @@ localrun <- Sys.getenv("OPENSAFELY_BACKEND") %in% c("", "expectations") if (localrun) { - jsonlite::write_json( - study_dates, - path = here::here("analysis", "0-lib", "study_dates.json"), - pretty = TRUE, auto_unbox = TRUE - ) + # jsonlite::write_json( + # study_dates, + # path = here::here("analysis", "0-lib", "study_dates.json"), + # pretty = TRUE, auto_unbox = TRUE + # ) jsonlite::write_json( - split(campaign_info, f = campaign_info$campaign_start_date) |> lapply(as.list), - path = here::here("analysis", "0-lib", "campaign_info.json"), - pretty = TRUE, auto_unbox = TRUE, + split(cohort_info, f = cohort_info$cohort_id) |> lapply(as.list), + path = here::here("analysis", "0-lib", "cohort_info.json"), + pretty = TRUE, + auto_unbox = TRUE ) } diff --git a/analysis/0-lib/study_dates.json b/analysis/0-lib/study_dates.json deleted file mode 100644 index 82a13bc..0000000 --- a/analysis/0-lib/study_dates.json +++ /dev/null @@ -1,5 +0,0 @@ -{ - "firstpossiblevax_date": "2020-07-01", - "start_date": "2020-12-07", - "end_date": "2026-03-31" -} diff --git a/analysis/1-extract/codelists.py b/analysis/1-extract/codelists.py index 62e65d6..57c8935 100644 --- a/analysis/1-extract/codelists.py +++ b/analysis/1-extract/codelists.py @@ -264,8 +264,25 @@ column="code", ) +# Outcomes ----------------------------- +#TODO mild outcomes # COVID-19 -covid_icd10 = ["U071", "U072", "U109"] +covid_icd10 = ["U071", "U072", "U109", "U075"] + +# Flu +flu_icd10 = codelist_from_csv( + "codelists/opensafely-influenza-identification-secondary-care.csv", + column="code", +) + +# RSV +rsv_icd10 = codelist_from_csv( + "codelists/opensafely-rsv-identification-secondary-care.csv", + column="code", +) + +#TODO: zoster + # Extended subgroups diff --git a/analysis/1-extract/dataset_definition.py b/analysis/1-extract/dataset_definition.py new file mode 100644 index 0000000..fa39eef --- /dev/null +++ b/analysis/1-extract/dataset_definition.py @@ -0,0 +1,124 @@ +########################## +# purpose: +# extract for a given cohort +# - patient-level baseline demographic + clinical variables +# - disease-specific vaccination variables, and +# - disease-specific outcomes +# - death date and deregistration +# for a given cohort +########################## + +from json import loads +from pathlib import Path +from datetime import date + +from ehrql import ( + get_parameter, + create_dataset, + days, + weeks, +) + +from ehrql.tables.tpp import ( + patients, + practice_registrations, + ons_deaths, + # apcs, + # vaccinations, +) + +import codelists +from variables_function import * + + +# ------------------------------------------------------------------------------ +# Cohort parameters +# ------------------------------------------------------------------------------ + +cohort_id = get_parameter(name="cohort_id") + +all_cohort_info = loads( + Path("analysis/0-lib/cohort_info.json").read_text() +) + +cohort_info = all_cohort_info[cohort_id] + +target = cohort_info["target"] +cohort_start_date = date.fromisoformat(cohort_info["cohort_start_date"]) +cohort_end_date = date.fromisoformat(cohort_info["cohort_end_date"]) + +# Age calculated the day before cohort start +age_calculation_date = cohort_start_date - days(1) + + +# ------------------------------------------------------------------------------ +# Initialise dataset +# ------------------------------------------------------------------------------ + +dataset = create_dataset() +dataset.configure_dummy_data(population_size=1000) + + +# ------------------------------------------------------------------------------ +# Population definition +# ------------------------------------------------------------------------------ + +registered_patients = practice_registrations.for_patient_on(cohort_start_date) + +registered = registered_patients.exists_for_patient() +registered_start_date = registered_patients.start_date + +alive_ONS = (ons_deaths.date > cohort_start_date) | ons_deaths.date.is_null() +alive_GP = (patients.date_of_death > cohort_start_date) | patients.date_of_death.is_null() +alive = alive_ONS & alive_GP + +eligibility_age = patients.age_on(age_calculation_date) + +dataset.define_population( + registered + & (registered_start_date <= (cohort_start_date - weeks(12))) + & alive + & (eligibility_age >= 12) + & (eligibility_age <= 104) + & patients.sex.is_in(["male", "female"]) +) + + +# ------------------------------------------------------------------------------ +# Baseline variables +# ------------------------------------------------------------------------------ + +demographic_variables(dataset=dataset, index_date=cohort_start_date) + +primis_variables(dataset=dataset, index_date=cohort_start_date) + +extended_subgroups(dataset=dataset, index_date=cohort_start_date) + +# Registration end date for the registration active at cohort start +dataset.deregistered_date = registered_patients.end_date + +# All-cause death date +# TODO: include ehr date? +dataset.death_date = ons_deaths.date + + +# ------------------------------------------------------------------------------ +# Disease-specific variables: vaccination + outcomes +# ------------------------------------------------------------------------------ +# TODO: add mild outcomes +if target == "COVID-19": + add_target_vaccines(dataset, ["SARS-2 CORONAVIRUS"], cohort_start_date, cohort_end_date) + add_target_sev_outcomes(dataset, codelists.covid_icd10, cohort_start_date, cohort_end_date) + +elif target == "Influenza": + add_target_vaccines(dataset, ["INFLUENZA"], cohort_start_date, cohort_end_date) + add_target_sev_outcomes(dataset, codelists.flu_icd10, cohort_start_date, cohort_end_date) + +elif target == "RSV": + add_target_vaccines(dataset, ["HUMAN RESPIRATORY SYNCYTIAL VIRUS"], cohort_start_date, cohort_end_date) + add_target_sev_outcomes(dataset, codelists.rsv_icd10, cohort_start_date, cohort_end_date) + +#TODO: Zoster + +else: + raise ValueError(f"Unknown target: {target}") \ No newline at end of file diff --git a/analysis/1-extract/dataset_definition_fixed.py b/analysis/1-extract/dataset_definition_fixed.py deleted file mode 100644 index d181cbf..0000000 --- a/analysis/1-extract/dataset_definition_fixed.py +++ /dev/null @@ -1,101 +0,0 @@ - -########################## -# purpose: -# extract time-invariant patient information as at a fixed time at the end of observation period -# time-invariant info: sex, ethnicity (assumed invariant), death date -# alsoextract count of previous vaccines as at that date -########################## - -# import libraries and functions - -from json import loads -from pathlib import Path - -from ehrql import ( - case, - create_dataset, - days, - when, -# minimum_of, -# maximum_of -) -from ehrql.tables.tpp import ( - patients, - practice_registrations, - vaccinations, - clinical_events, - ons_deaths -) -# import codelists -import codelists - -study_dates = loads( - Path("analysis/0-lib/study_dates.json").read_text(), -) - -# Change these in ./0-lib/deisgn.R if necessary -start_date = study_dates["start_date"] -end_date = study_dates["end_date"] - -# all covid-19 vaccination events -covid_vaccinations = ( - vaccinations - .where(vaccinations.target_disease.is_in(["SARS-2 CORONAVIRUS"])) - .where(vaccinations.date.is_on_or_before(end_date)) - .sort_by(vaccinations.date) -) - -registered_at_any_time_in_observation_period = practice_registrations.where( - # starting during period - practice_registrations.start_date.is_on_or_between(start_date, end_date) | - - # ending during period - practice_registrations.end_date.is_on_or_between(start_date, end_date) | - - # starting before and ending after - ( - practice_registrations.start_date.is_on_or_before(start_date) & - (practice_registrations.end_date.is_on_or_after(end_date + days(1)) | practice_registrations.end_date.is_null()) - ) -) - -# initialise dataset -dataset = create_dataset() -dataset.configure_dummy_data(population_size=1000) - -# define dataset poppulation -dataset.define_population( - registered_at_any_time_in_observation_period.exists_for_patient() - & (patients.age_on(end_date) >= 12) # only include people who are aged 12 or over during at least one season -) - -# patient sex -dataset.sex = patients.sex - -# patient ethnicity -# note that enthicity is documented using codelists, not as a categorical variable -# if ethnicity was transferred from another practice, the date may not have been captured, and will default to 1900-01-01 -# we choose here to look at the last known ethnicity recorded _across the entire record_ -# rather than as known/recorded on the vaccination date, even though this looks into the future. - -ethnicity = (clinical_events - .where(clinical_events.snomedct_code.is_in(codelists.ethnicity16)) - .sort_by(clinical_events.date) - .where(clinical_events.date.is_on_or_before(end_date)) - .last_for_patient() -) - -# # ethnicity using 5 groups + unknown -dataset.ethnicity5 = ethnicity.snomedct_code.to_category(codelists.ethnicity5) -# -# # ethnicity using 16 groups + unknown -dataset.ethnicity16 = ethnicity.snomedct_code.to_category(codelists.ethnicity16) - -# patient death date -dataset.death_date = ons_deaths.date -dataset.covid_death_date = case( - when(ons_deaths.cause_of_death_is_in(codelists.covid_icd10)).then(ons_deaths.date), - otherwise=None -) -# number of covid vaccines as at end of observation period -dataset.covid_vax_count = covid_vaccinations.count_for_patient() diff --git a/analysis/1-extract/dataset_definition_snapshot.py b/analysis/1-extract/dataset_definition_snapshot.py deleted file mode 100644 index 26096b6..0000000 --- a/analysis/1-extract/dataset_definition_snapshot.py +++ /dev/null @@ -1,188 +0,0 @@ - -########################## -# extract patient information as at a fixed date -# vaccination dates and info will be joined onto this later using the "varying" dataset -########################## - -# import libraries and functions - -from json import loads -from pathlib import Path -from datetime import datetime - -from ehrql import ( -# show, - get_parameter, -# case, - create_dataset, - days, - weeks, -# when, -# minimum_of, -# maximum_of -) -from ehrql.tables.tpp import ( - patients, - practice_registrations, - apcs, - ons_deaths, -) - -# import codelists -from codelists import * - -from variables_function import * - -# get arguments supplied in the yaml file, stored in argv[1], argv[2], etc -snapshot_date_string = get_parameter(name="snapshot_date") -snapshot_date = datetime.strptime(get_parameter(name="snapshot_date"), '%Y%m%d').strftime('%Y-%m-%d') - -# extract campaign-specific dates -all_campaign_info = loads( - Path("analysis/0-lib/campaign_info.json").read_text(), -) - -campaign_info = all_campaign_info[snapshot_date] - -# get date on which age should be calculated -# this is different to snapshot date because we allow people to be eligible if they reach the age during the campaign -age_calculation_date = campaign_info["age_date"] - -# end of follow-up for campaign (= day before next campaign)) -end_of_campaign_date = campaign_info["final_milestone_date"] - -# initialise dataset - -dataset = create_dataset() - -dataset.configure_dummy_data(population_size=1000) - -registered_patients = practice_registrations.for_patient_on(snapshot_date) - -registered = registered_patients.exists_for_patient() -registered_start_date = registered_patients.start_date -alive_ONS = (ons_deaths.date>snapshot_date) | ons_deaths.date.is_null() -alive_GP = (patients.date_of_death>snapshot_date) | patients.date_of_death.is_null() -alive = (alive_ONS & alive_GP) -eligibility_age = patients.age_on(age_calculation_date) - -# define dataset poppulation -dataset.define_population( - registered - & (registered_start_date <= (snapshot_date - weeks(12))) - & alive - & (eligibility_age >= 12) & (eligibility_age <= 104) - & (patients.sex.is_in(["male", "female"])) -) - -## demographic variables ---- - -dataset.eligibility_age = eligibility_age - -demographic_variables(dataset = dataset, index_date = snapshot_date) - -## primis variables ---- - -primis_variables(dataset = dataset, index_date = snapshot_date) - -## extended subgroups ---- - -extended_subgroups(dataset = dataset, index_date = snapshot_date) - -# other_cx_variables(dataset = dataset, index_date = snapshot_date) - -# COVID-19 vaccination history ---- - -covid_vaccinations = ( - vaccinations - .where(vaccinations.target_disease.is_in(["SARS-2 CORONAVIRUS"])) - .sort_by(vaccinations.date) -) - -# retrieve first 3 vaccines before and first 2 vaccines on or after snapshot date - - -# add next 2 recorded covid vaccination -add_n_vaccines( - dataset = dataset, - index_date = snapshot_date, - target_disease = "SARS-2 CORONAVIRUS", - name = "covid_vax", - direction = "on_or_after", - number_of_vaccines = 2, - minimum_gap = 13 -) - -# add most recent 3 vaccines prior to snapshot date -add_n_vaccines( - dataset = dataset, - index_date = snapshot_date, - target_disease = "SARS-2 CORONAVIRUS", - name = "covid_vax_prior", - direction = "before", - number_of_vaccines = 3, - minimum_gap = 13 -) - -# count total number of vaccines prior to vaccine date -dataset.covid_vax_prior_count = ( - covid_vaccinations - .where(vaccinations.date.is_before(snapshot_date)) - .count_for_patient() -) - -# Deregistration dates after the snapshot date -dataset.deregistered_date = registered_patients.end_date - -# Covid-19 outcomes - -all_covid_admissions = ( - apcs - .where(apcs.all_diagnoses.contains_any_of(codelists.covid_icd10)) - .where(apcs.admission_method.is_in(["21", "22", "23", "24", "25", "2A", "2B", "2C", "2D", "28"])) - .where(apcs.patient_classification == "1") # Ordinary admissions only - .where(apcs.admission_date.is_on_or_between(snapshot_date,end_of_campaign_date)) - .sort_by(apcs.admission_date) -) - -# covid-related admission -dataset.covid_admitted_date = ( - all_covid_admissions - .first_for_patient() - .admission_date -) - -# covid-related admission, primary diagnosis only -dataset.covid_admitted_primary_date = ( - all_covid_admissions - .where(apcs.primary_diagnosis.is_in(codelists.covid_icd10)) - .first_for_patient() - .admission_date -) - -# first covid-related admission, length of stay -#dataset.covid_admitted_first_los = (all_covid_admissions.first_for_patient().discharge_date) - dataset.covid_admitted_date - - -# length of stay across all covid admissions with admission within the campaign date -dataset.covid_admitted_los = ( - all_covid_admissions.discharge_date - all_covid_admissions.admission_date -).days.sum_for_patient() - -# covid-related critical care admission -dataset.covid_critcare_date = ( - all_covid_admissions - .where(apcs.days_in_critical_care>0) - .first_for_patient() - .admission_date -) - - -# covid-related death -#dataset.covid_death_date = ( -# ons_deaths -# .cause_of_death_is_in(codelists.covid_icd10) -# .date -#) - - diff --git a/analysis/1-extract/dataset_definition_varying.py b/analysis/1-extract/dataset_definition_varying.py deleted file mode 100644 index 92b5e2f..0000000 --- a/analysis/1-extract/dataset_definition_varying.py +++ /dev/null @@ -1,120 +0,0 @@ -# _________________________________________________ -# Purpose: -# extract patient information as at each vaccination date -# column set 1 = details of first vaccination event -# column set 2 = details of second vaccination event -# _________________________________________________ - -# import libraries and functions - -from json import loads -from pathlib import Path - -from ehrql import ( - # case, - create_dataset, - # days, - # when, - # minimum_of, - # maximum_of, - claim_permissions -) -from ehrql.tables.tpp import ( - patients, - practice_registrations, - vaccinations, -# clinical_events, -# ons_deaths, -# addresses, -) -# import codelists -from codelists import * - - -study_dates = loads( - Path("analysis/0-lib/study_dates.json").read_text(), -) - -# Change these in ./0-lib/design.R if necessary -start_date = study_dates["start_date"] -end_date = study_dates["end_date"] - -#import function for clinical variables -from variables_function import * - -# all covid-19 vaccination events -covid_vaccinations = ( - vaccinations - .where(vaccinations.target_disease.is_in(["SARS-2 CORONAVIRUS"])) - .sort_by(vaccinations.date) -) - -# initialise dataset -dataset = create_dataset() -dataset.configure_dummy_data(population_size=1000) - -# define dataset poppulation -dataset.define_population( - covid_vaccinations.exists_for_patient() - & (patients.age_on(end_date) >= 12) # only include people who are aged 2 or over during at least one season -) - - -# Arbitrary date guaranteed to be before any vaccination events of interest -previous_vax_date = "1899-01-01" - -# loop over first, second, ..., nth vaccination event for each person -# extract info on vaccination date and product, and basic demographics -for i in range(1, 16+1): - - suffix = f"_{i}" - - # vaccine variables - current_vax = covid_vaccinations.where(covid_vaccinations.date>previous_vax_date).first_for_patient() - dataset.add_column(f"covid_vax_{i}_date", current_vax.date) - dataset.add_column(f"covid_vax_product_{i}", current_vax.product_name) - - # registration variables - registration = practice_registrations.for_patient_on(current_vax.date) - dataset.add_column(f"registered_{i}", registration.exists_for_patient()) - dataset.add_column(f"deregistered_{i}_date", registration.end_date) - - # deomgraphic variables - demographic_variables(dataset = dataset, index_date = current_vax.date, var_name_suffix = suffix) - - # primis variables -# primis_variables(dataset = dataset, index_date = current_vax.date, var_name_suffix = suffix) -# dataset.add_column(f"primis_atrisk_{i}", primis_atrisk(current_vax.date)) # at risk - - # Extended subgroups -# extended_subgroups(dataset = dataset, index_date = current_vax.date, var_name_suffix = suffix) - -# other_cx_variables(dataset = dataset, index_date = current_vax.date, var_name_suffix = suffix) - - previous_vax_date = current_vax.date - - - -### Add event-level data extract to test equivalence of ELD extract versus patient-level extract - -# all covid-19 vaccination events -# i don't think it's possible to restrict to first 16 events (as above for patient level data) -# so I'll do this post-extract - -# covid_vaccinations_ELD = covid_vaccinations.where(covid_vaccinations.date>"1899-01-01") - -# Some of the tables or features you are using require special permission to use with real -# patient data. The permissions needed are:`dataset.add_event_table()` method -# You can continue to work on your code using dummy data by “claiming” the required permisions: -# Note that you will only be able to run your code against real data if you actually have these -# permissions assigned by the OpenSAFELY team. For more information see: https://docs.opensafely.org/ehrql/reference/language/#permissions - -claim_permissions("event_level_data") - -covid_vaccinations_ELD = covid_vaccinations -dataset.add_event_table( - "vaccinations", - vax_date = covid_vaccinations_ELD.date, - vax_product = covid_vaccinations_ELD.product_name, - age = patients.age_on(covid_vaccinations_ELD.date), -) diff --git a/analysis/1-extract/dummy-data/dummy_covid_2025_26.arrow b/analysis/1-extract/dummy-data/dummy_covid_2025_26.arrow new file mode 100644 index 0000000..8e4c4f8 Binary files /dev/null and b/analysis/1-extract/dummy-data/dummy_covid_2025_26.arrow differ diff --git a/analysis/1-extract/dummy-data/dummy_covid_ba1_ba2.arrow b/analysis/1-extract/dummy-data/dummy_covid_ba1_ba2.arrow new file mode 100644 index 0000000..8f5a74d Binary files /dev/null and b/analysis/1-extract/dummy-data/dummy_covid_ba1_ba2.arrow differ diff --git a/analysis/1-extract/dummy-data/dummy_covid_jn1.arrow b/analysis/1-extract/dummy-data/dummy_covid_jn1.arrow new file mode 100644 index 0000000..af867a6 Binary files /dev/null and b/analysis/1-extract/dummy-data/dummy_covid_jn1.arrow differ diff --git a/analysis/1-extract/dummy-data/dummy_covid_wave_1.arrow b/analysis/1-extract/dummy-data/dummy_covid_wave_1.arrow new file mode 100644 index 0000000..7634741 Binary files /dev/null and b/analysis/1-extract/dummy-data/dummy_covid_wave_1.arrow differ diff --git a/analysis/1-extract/dummy-data/dummy_fixed.arrow b/analysis/1-extract/dummy-data/dummy_fixed.arrow deleted file mode 100644 index 8c2314f..0000000 Binary files a/analysis/1-extract/dummy-data/dummy_fixed.arrow and /dev/null differ diff --git a/analysis/1-extract/dummy-data/dummy_flu_2023_24.arrow b/analysis/1-extract/dummy-data/dummy_flu_2023_24.arrow new file mode 100644 index 0000000..dd2834d Binary files /dev/null and b/analysis/1-extract/dummy-data/dummy_flu_2023_24.arrow differ diff --git a/analysis/1-extract/dummy-data/dummy_flu_2024_25.arrow b/analysis/1-extract/dummy-data/dummy_flu_2024_25.arrow new file mode 100644 index 0000000..0070851 Binary files /dev/null and b/analysis/1-extract/dummy-data/dummy_flu_2024_25.arrow differ diff --git a/analysis/1-extract/dummy-data/dummy_flu_2025_26.arrow b/analysis/1-extract/dummy-data/dummy_flu_2025_26.arrow new file mode 100644 index 0000000..8e4c4f8 Binary files /dev/null and b/analysis/1-extract/dummy-data/dummy_flu_2025_26.arrow differ diff --git a/analysis/1-extract/dummy-data/dummy_rsv_2023_24.arrow b/analysis/1-extract/dummy-data/dummy_rsv_2023_24.arrow new file mode 100644 index 0000000..dd2834d Binary files /dev/null and b/analysis/1-extract/dummy-data/dummy_rsv_2023_24.arrow differ diff --git a/analysis/1-extract/dummy-data/dummy_rsv_2024_25.arrow b/analysis/1-extract/dummy-data/dummy_rsv_2024_25.arrow new file mode 100644 index 0000000..0070851 Binary files /dev/null and b/analysis/1-extract/dummy-data/dummy_rsv_2024_25.arrow differ diff --git a/analysis/1-extract/dummy-data/dummy_rsv_2025_26.arrow b/analysis/1-extract/dummy-data/dummy_rsv_2025_26.arrow new file mode 100644 index 0000000..8e4c4f8 Binary files /dev/null and b/analysis/1-extract/dummy-data/dummy_rsv_2025_26.arrow differ diff --git a/analysis/1-extract/dummy-data/dummy_snapshot_20201207.arrow b/analysis/1-extract/dummy-data/dummy_snapshot_20201207.arrow deleted file mode 100644 index 3fe4729..0000000 Binary files a/analysis/1-extract/dummy-data/dummy_snapshot_20201207.arrow and /dev/null differ diff --git a/analysis/1-extract/dummy-data/dummy_snapshot_20201208.arrow b/analysis/1-extract/dummy-data/dummy_snapshot_20201208.arrow deleted file mode 100644 index da14540..0000000 Binary files a/analysis/1-extract/dummy-data/dummy_snapshot_20201208.arrow and /dev/null differ diff --git a/analysis/1-extract/dummy-data/dummy_snapshot_20210901.arrow b/analysis/1-extract/dummy-data/dummy_snapshot_20210901.arrow deleted file mode 100644 index 91c8c62..0000000 Binary files a/analysis/1-extract/dummy-data/dummy_snapshot_20210901.arrow and /dev/null differ diff --git a/analysis/1-extract/dummy-data/dummy_snapshot_20210906.arrow b/analysis/1-extract/dummy-data/dummy_snapshot_20210906.arrow deleted file mode 100644 index 5dafc26..0000000 Binary files a/analysis/1-extract/dummy-data/dummy_snapshot_20210906.arrow and /dev/null differ diff --git a/analysis/1-extract/dummy-data/dummy_snapshot_20220321.arrow b/analysis/1-extract/dummy-data/dummy_snapshot_20220321.arrow deleted file mode 100644 index 8247a3d..0000000 Binary files a/analysis/1-extract/dummy-data/dummy_snapshot_20220321.arrow and /dev/null differ diff --git a/analysis/1-extract/dummy-data/dummy_snapshot_20220829.arrow b/analysis/1-extract/dummy-data/dummy_snapshot_20220829.arrow deleted file mode 100644 index d8119e3..0000000 Binary files a/analysis/1-extract/dummy-data/dummy_snapshot_20220829.arrow and /dev/null differ diff --git a/analysis/1-extract/dummy-data/dummy_snapshot_20230403.arrow b/analysis/1-extract/dummy-data/dummy_snapshot_20230403.arrow deleted file mode 100644 index 691a033..0000000 Binary files a/analysis/1-extract/dummy-data/dummy_snapshot_20230403.arrow and /dev/null differ diff --git a/analysis/1-extract/dummy-data/dummy_snapshot_20230828.arrow b/analysis/1-extract/dummy-data/dummy_snapshot_20230828.arrow deleted file mode 100644 index 59ed617..0000000 Binary files a/analysis/1-extract/dummy-data/dummy_snapshot_20230828.arrow and /dev/null differ diff --git a/analysis/1-extract/dummy-data/dummy_snapshot_20240415.arrow b/analysis/1-extract/dummy-data/dummy_snapshot_20240415.arrow deleted file mode 100644 index a0faea6..0000000 Binary files a/analysis/1-extract/dummy-data/dummy_snapshot_20240415.arrow and /dev/null differ diff --git a/analysis/1-extract/dummy-data/dummy_snapshot_20240930.arrow b/analysis/1-extract/dummy-data/dummy_snapshot_20240930.arrow deleted file mode 100644 index d2a1166..0000000 Binary files a/analysis/1-extract/dummy-data/dummy_snapshot_20240930.arrow and /dev/null differ diff --git a/analysis/1-extract/dummy-data/dummy_snapshot_20250331.arrow b/analysis/1-extract/dummy-data/dummy_snapshot_20250331.arrow deleted file mode 100644 index 7a1e6b3..0000000 Binary files a/analysis/1-extract/dummy-data/dummy_snapshot_20250331.arrow and /dev/null differ diff --git a/analysis/1-extract/dummy-data/dummy_snapshot_20250929.arrow b/analysis/1-extract/dummy-data/dummy_snapshot_20250929.arrow deleted file mode 100644 index a9f6800..0000000 Binary files a/analysis/1-extract/dummy-data/dummy_snapshot_20250929.arrow and /dev/null differ diff --git a/analysis/1-extract/dummy-data/dummy_varying.arrow b/analysis/1-extract/dummy-data/dummy_varying.arrow deleted file mode 100644 index ecffe6e..0000000 Binary files a/analysis/1-extract/dummy-data/dummy_varying.arrow and /dev/null differ diff --git a/analysis/1-extract/dummy_dataset_definition_snapshot.R b/analysis/1-extract/dummy_dataset_definition.R similarity index 54% rename from analysis/1-extract/dummy_dataset_definition_snapshot.R rename to analysis/1-extract/dummy_dataset_definition.R index 3df7c12..64bbea8 100644 --- a/analysis/1-extract/dummy_dataset_definition_snapshot.R +++ b/analysis/1-extract/dummy_dataset_definition.R @@ -1,6 +1,6 @@ # _________________________________________________ # Purpose: -# Create a dummy dataset for`dataset_definition_fixed.py` script +# Create a dummy dataset for the disease-specific dataset definition # _________________________________________________ @@ -43,29 +43,28 @@ population_size <- 10000 # we do this because some dplyr operations to not preserve date attributes, so dates will be converted to numerics -all_snapshot_dates <- campaign_info |> - filter( - campaign_start_date >= as.Date("2020-12-07"), - campaign_start_date <= as.Date("2026-03-31") - ) |> - pull(campaign_start_date) +cohort_ids <- cohort_info$cohort_id -for (snapshot_date in all_snapshot_dates) { - # use these if selecting a specific date to run - # snapshot_date <- as.Date("2020-12-07") - # snapshot_date <- campaign_info |> filter(campaign_label == "Spring 2025") |> pull(campaign_start_date) +for (cohort_id in cohort_ids) { - snapshot_date <- as.Date(snapshot_date) + cohort_row <- cohort_info |> + filter(cohort_id == .env$cohort_id) - print(snapshot_date) + target <- cohort_row$target + cohort <- cohort_row$cohort + snapshot_date <- as.Date(cohort_row$cohort_start_date) + cohort_end_date <- as.Date(cohort_row$cohort_end_date) + + print(glue("{cohort_id}: {target} {cohort}")) snapshot_day <- 0L + cohort_end_day <- as.integer(cohort_end_date - snapshot_date) # set the variables and functions that are known a-priori to the simulation engine # ie, defined and accessible outside of the scope of the dataset known_variables <- c( - "index_date", - "snapshot_day" + "snapshot_day", + "cohort_end_day" ) @@ -75,39 +74,69 @@ for (snapshot_date in all_snapshot_dates) { # ..n as a place holder for the length of the variable sim_list <- lst( - deregistered_day = bn_node( - ~ as.integer(runif(n = ..n, snapshot_day, snapshot_day + 1200)), - missing_rate = ~0.99 - ), - age = bn_node( - ~ as.integer(rnormt(n = ..n, mean = 60, sd = 14, range = c(12, 104))) - ), - eligibility_age = bn_node( - ~ age + rbernoulli(n = ..n, p = 0.25) - ), - imd = bn_node( - ~ as.integer(plyr::round_any(runif(n = ..n, 100, 32000), 100)), - missing = ~0.05 - ), - stp = bn_node( - ~ factor(as.integer(runif(n = ..n, 1, 36)), levels = 1:36) - ), - region = bn_node( - variable_formula = ~ rfactor(n = ..n, levels = c( - "North East", - "North West", - "Yorkshire and The Humber", - "East Midlands", - "West Midlands", - "East", - "London", - "South East", - "South West" - ), p = c(0.2, 0.2, 0.3, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05)) - ), - carehome_status = bn_node( - ~ rbernoulli(n = ..n, p = if_else(age < 65, 0.01, 0.2)), - ), + # demographics + sex = bn_node( + ~ rfactor(n = ..n, levels = c("female", "male", "intersex", "unknown"), p = c(0.51, 0.49, 0, 0)) + ), + + age = bn_node( + ~ as.integer(rnormt(n = ..n, mean = 60, sd = 14, range = c(12, 104))) + ), + + ethnicity5 = bn_node(variable_formula = ~ ethnicity_16_to_5(ethnicity16), needs = "ethnicity16"), + ethnicity16 = bn_node( + variable_formula = ~ rfactor( + n = ..n, + levels = c( + "White - British", + "White - Irish", + "White - Any other White background", + "Mixed - White and Black Caribbean", + "Mixed - White and Black African", + "Mixed - White and Asian", + "Mixed - Any other mixed background", + "Asian or Asian British - Indian", + "Asian or Asian British - Pakistani", + "Asian or Asian British - Bangladeshi", + "Asian or Asian British - Any other Asian background", + "Black or Black British - Caribbean", + "Black or Black British - African", + "Black or Black British - Any other Black background", + "Other Ethnic Groups - Chinese", + "Other Ethnic Groups - Any other ethnic group" + ), + p = c( + 0.5, 0.05, 0.05, # White + 0.025, 0.025, 0.025, 0.025, # Mixed + 0.025, 0.025, 0.025, 0.025, # Asian + 0.033, 0.033, 0.034, # Black + 0.05, 0.05 # Other + ) + ), + missing_rate = ~0.1, + ), + + region = bn_node( + variable_formula = ~ rfactor(n = ..n, levels = c( + "North East", "North West", "Yorkshire and The Humber", + "East Midlands", "West Midlands", "East", + "London", "South East", "South West" + ), p = c(0.2, 0.2, 0.3, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05)) + ), + + stp = bn_node( + ~ factor(as.integer(runif(n = ..n, 1, 36)), levels = 1:36) + ), + + imd = bn_node( + ~ as.integer(plyr::round_any(runif(n = ..n, 100, 32000), 100)), + missing_rate = ~0.05 + ), + + carehome_status = bn_node( + ~ rbernoulli(n = ..n, p = if_else(age < 65, 0.01, 0.2)), + needs = "age" + ), # PRIMIS crd = bn_node( ~ rbernoulli(n = ..n, p = 0.02), @@ -197,61 +226,54 @@ for (snapshot_date in all_snapshot_dates) { homeless = bn_node( ~ rbernoulli(n = ..n, p = 0.02), ), - # covid vaccines - covid_vax_1_day = bn_node( - ~ runif(n = ..n, snapshot_day, snapshot_day + 200), - missing_rate = ~0.2, - ), - covid_vax_2_day = bn_node( - ~ runif(n = ..n, covid_vax_1_day + 14, covid_vax_1_day + 400), - missing_rate = ~0.95, - needs = "covid_vax_1_day" - ), - covid_vax_1_product = bn_node(~ rcat(n = ..n, c("pfizer", "az", "moderna"), c(0.5, 0.3, 0.2)), needs = "covid_vax_1_day"), - covid_vax_2_product = bn_node(~ if_else(runif(..n) < 0.98, covid_vax_1_product, "pfizer"), needs = "covid_vax_2_day"), - - covid_vax_prior_1_day = bn_node( - ~ runif(n = ..n, snapshot_day - 400, snapshot_day - 200), - missing_rate = ~0.1, - needs = "covid_vax_1_day" - ), - covid_vax_prior_2_day = bn_node( - ~ runif(n = ..n, covid_vax_prior_1_day - 400, covid_vax_prior_1_day - 200), - missing_rate = ~0.1, - needs = "covid_vax_prior_1_day" - ), - covid_vax_prior_3_day = bn_node( - ~ runif(n = ..n, covid_vax_prior_2_day - 400, covid_vax_prior_2_day - 200), - missing_rate = ~0.1, - needs = "covid_vax_prior_2_day" - ), - - covid_vax_prior_1_product = bn_node(~ rcat(n = ..n, c("pfizer", "az", "moderna"), c(0.5, 0.3, 0.2)), needs = "covid_vax_prior_1_day"), - covid_vax_prior_2_product = bn_node(~ rcat(n = ..n, c("pfizer", "az", "moderna"), c(0.5, 0.3, 0.2)), needs = "covid_vax_prior_2_day"), - covid_vax_prior_3_product = bn_node(~ rcat(n = ..n, c("pfizer", "az", "moderna"), c(0.5, 0.3, 0.2)), needs = "covid_vax_prior_3_day"), - - covid_vax_prior_count = bn_node(~ as.integer(runif(n = ..n, 0, 15))), # in dummy data, this will not match total vax count in "time-varying" dataset, but that's ok - - covid_admitted_day = bn_node( - ~ as.integer(runif(n = ..n, snapshot_day, snapshot_day + 300)), - missing_rate = ~0.7 - ), - - covid_admitted_los = bn_node( - ~ as.integer(runif(n = ..n, 1, 10)), - needs = "covid_admitted_day" - ), - - covid_admitted_primary_day = bn_node( - ~ if_else(rbernoulli(n = ..n, p = 0.5) == 1, covid_admitted_day, NA_integer_) - ), - - covid_critcare_day = bn_node( - ~covid_admitted_day, - needs = "covid_admitted_day", - missing_rate = ~0.7 - ), + # registration + deregistered_day = bn_node( + ~ as.integer(runif(n = ..n, snapshot_day, snapshot_day + 1200)), + missing_rate = ~0.99 + ), + + # all-cause death + death_day = bn_node( + ~ as.integer(runif(n = ..n, snapshot_day, snapshot_day + 3000)), + missing_rate = ~0.98 + ), + + # disease-specific vaccination variables + last_vax_date_before_start_day = bn_node( + ~ runif(n = ..n, snapshot_day - 400, snapshot_day - 1), + missing_rate = ~0.5 + ), + + last_vax_product_before_start = bn_node( + ~ rcat(n = ..n, c("prodA", "prodB"), c(0.5, 0.5)), + needs = "last_vax_date_before_start_day" + ), + + first_vax_date_after_start_day = bn_node( + ~ runif(n = ..n, snapshot_day, snapshot_day + 200), + missing_rate = ~0.5 + ), + + first_vax_product_after_start = bn_node( + ~ rcat(n = ..n, c("prodA", "prodB"), c(0.5, 0.5)), + needs = "first_vax_date_after_start_day" + ), + # disease-specific severe outcomes + admitted_day = bn_node( + ~ as.integer(runif(n = ..n, snapshot_day, snapshot_day + 300)), + missing_rate = ~0.7 + ), + + admitted_primary_day = bn_node( + ~ if_else(rbernoulli(n = ..n, p = 0.5) == 1, admitted_day, NA_integer_), + needs = "admitted_day" + ), + + disease_death_day = bn_node( + ~ as.integer(runif(n = ..n, snapshot_day, snapshot_day + 300)), + missing_rate = ~0.95 ) +) # check and create the simulation object, including all dependencies, topological orders, etc bn <- bn_create(sim_list, known_variables = known_variables) @@ -270,15 +292,17 @@ for (snapshot_date in all_snapshot_dates) { # do some post simulation processing for features that are not easily handled by the simulation configuration dummydata_processed <- dummydata %>% - # convert integer days to dates since index date and rename vars + mutate( + cohort_start_date = snapshot_date, + cohort_end_date = cohort_end_date + ) %>% mutate(across(ends_with("_day"), ~ as.Date(as.character(snapshot_date + .)))) %>% rename_with(~ str_replace(., "_day", "_date"), ends_with("_day")) # save the dataset in arrow format - snapshot_date_compact <- format(snapshot_date, "%Y%m%d") - write_feather(dummydata_processed, sink = here("analysis", "1-extract", "dummy-data", glue("dummy_snapshot_{snapshot_date_compact}.arrow"))) + write_feather(dummydata_processed, sink = here("analysis", "1-extract", "dummy-data", glue("dummy_{cohort_id}.arrow"))) } diff --git a/analysis/1-extract/dummy_dataset_definition_fixed.R b/analysis/1-extract/dummy_dataset_definition_fixed.R deleted file mode 100644 index 3d5d888..0000000 --- a/analysis/1-extract/dummy_dataset_definition_fixed.R +++ /dev/null @@ -1,121 +0,0 @@ -# _________________________________________________ -# Purpose: -# Create a dummy dataset for`dataset_definition_fixed.py` script -# _________________________________________________ - - -# Import libraries and functions ---- - -library("tidyverse") -library("arrow") -library("here") -library("glue") - -# remotes::install_github("https://github.com/wjchulme/dd4d") #package for more convenient data simulation -library("dd4d") - -# Import custom functions -source(here("analysis", "0-lib", "design.R")) - -# Define and simulate the dataset ---- - -# Set the size of the dataset -population_size <- 10000 - -# set the index date for date variables -# all variables will be defined as the number of days before or after this day -# and then at the end of the script they are transformed into dates -# we do this because some dplyr operations to not preserve date attributes, so dates will be converted to numerics -index_date <- study_dates$end_date - -index_day <- 0L - -# set the variables and functions that are known a-priori to the simulation engine -# ie, defined and accessible outside of the scope of the dataset -known_variables <- c( - "index_date", - "index_day" -) - - -# define the simulation configuration -# ie, a list of variables to simulate -# use the form _variable_name_ = bn_node(~ _formula_for_simulating_variable_, ) see help("bn_node") -# ..n as a place holder for the length of the variable - - - -sim_list <- lst( - sex = bn_node( - ~ rfactor(n = ..n, levels = c("female", "male", "intersex", "unknown"), p = c(0.51, 0.49, 0, 0)) - # missing_rate = ~0.001 # this is shorthand for ~(rbernoulli(n=..n, p = 0.2)) - ), - ethnicity5 = bn_node(variable_formula = ~ ethnicity_16_to_5(ethnicity16), needs = "ethnicity16"), - ethnicity16 = bn_node( - variable_formula = ~ rfactor( - n = ..n, - levels = c( - "White - British", - "White - Irish", - "White - Any other White background", - "Mixed - White and Black Caribbean", - "Mixed - White and Black African", - "Mixed - White and Asian", - "Mixed - Any other mixed background", - "Asian or Asian British - Indian", - "Asian or Asian British - Pakistani", - "Asian or Asian British - Bangladeshi", - "Asian or Asian British - Any other Asian background", - "Black or Black British - Caribbean", - "Black or Black British - African", - "Black or Black British - Any other Black background", - "Other Ethnic Groups - Chinese", - "Other Ethnic Groups - Any other ethnic group" - ), - p = c( - 0.5, 0.05, 0.05, # White - 0.025, 0.025, 0.025, 0.025, # Mixed - 0.025, 0.025, 0.025, 0.025, # Asian - 0.033, 0.033, 0.034, # Black - 0.05, 0.05 # Other - ) - ), - missing_rate = ~0.1, - ), - death_day = bn_node( - ~ as.integer(runif(n = ..n, index_day, index_day + 3000)), - missing_rate = ~0.98 - ), - - covid_death_day = bn_node( - ~death_day, - missing_rate = ~0.7, - needs = "death_day" - ), - covid_vax_count = bn_node(~ as.integer(runif(n = ..n, 0, 15))) # in dummy data, this will not match total vax count in "time-varying" dataset, but that's ok -) - -# check and create the simulation object, including all dependencies, topological orders, etc -bn <- bn_create(sim_list, known_variables = known_variables) - -# plot the network -bn_plot(bn) - -# plot the network (connected nodes only) -bn_plot(bn, connected_only = TRUE) - -# set the seed for the simulation -set.seed(10) - -# simulate the dataset -dummydata <- bn_simulate(bn, pop_size = population_size, keep_all = FALSE, .id = "patient_id") - -# do some post simulation processing for features that are not easily handled by the simulation configuration -dummydata_processed <- dummydata %>% - # convert integer days to dates since index date and rename vars - mutate(across(ends_with("_day"), ~ as.Date(as.character(index_date + .)))) %>% - rename_with(~ str_replace(., "_day", "_date"), ends_with("_day")) - - -# save the datasetin arrow format -write_feather(dummydata_processed, sink = here("analysis", "1-extract", "dummy-data", "dummy_fixed.arrow")) diff --git a/analysis/1-extract/dummy_dataset_definition_varying.R b/analysis/1-extract/dummy_dataset_definition_varying.R deleted file mode 100644 index b1863d1..0000000 --- a/analysis/1-extract/dummy_dataset_definition_varying.R +++ /dev/null @@ -1,354 +0,0 @@ -# _________________________________________________ -# Purpose: -# Create a dummy dataset for `dataset_definition_varying.py` script -# column set 1 = details of first vaccination event -# column set 2 = details of second vaccination event -# and so on -# _________________________________________________ - - -# Import libraries and functions ---- - -library("tidyverse") -library("arrow") -library("here") -library("glue") - -# remotes::install_github("https://github.com/wjchulme/dd4d") #package for more convenient data simulation -library("dd4d") - -# Import custom functions -source(here("analysis", "0-lib", "design.R")) - -# Define and simulate the dataset ---- - -# Set the size of the dataset -population_size <- 1000 - -# set the index date for date variables -# all variables will be defined as the number of days before or after this day -# and then at the end of the script they are transformed into dates -# we do this because some dplyr operations to not preserve date attributes, so dates will be converted to numerics - -# choose the june 2020 to pick up possible trial participants, but definitely exclude errors or unknown dates (1900-01-01, 2018-03-24, etc) -index_date <- study_dates$firstpossiblevax_date - -index_day <- 0L - - -start_day <- as.integer(study_dates$start_date - study_dates$index_date) - -# set the variables that are known a-priori to the simulation engine -# ie, defined ad accessible outside of the scope of the dataset -known_variables <- c( - "index_date", "start_date", - "index_day", "start_day" -) - -# define the simulation configuration -# ie, a list of variables to simulate -# use the form _variable_name_ = bn_node(~ _formula_for_simulating_variable_, ) see help("bn_node") -# ..n as a place holder for the length of the variable - - - -# define the list of vaccination-related variables to extract for each vaccination event -# this is essentially a date-product pair for up to the Nth vaccination -# it uses some very roughly realistic date ranges and product types -sim_list_vax_info <- lst( - - ## vaccination variables - - # covid vax any - covid_vax_1_day = bn_node(~ runif(n = ..n, start_day, start_day + 365), - missing_rate = ~0.05, ), - covid_vax_2_day = bn_node( - ~ runif(n = ..n, covid_vax_1_day + 14, covid_vax_1_day + 100), - missing_rate = ~0.001, - needs = "covid_vax_1_day" - ), - covid_vax_3_day = bn_node( - ~ runif(n = ..n, covid_vax_2_day + 250, covid_vax_2_day + 350), - missing_rate = ~0.3, - needs = "covid_vax_2_day" - ), - covid_vax_4_day = bn_node( - ~ runif(n = ..n, covid_vax_3_day + 300, covid_vax_3_day + 400), - missing_rate = ~0.7, - needs = "covid_vax_3_day" - ), - covid_vax_5_day = bn_node( - ~ runif(n = ..n, covid_vax_4_day + 300, covid_vax_4_day + 400), - missing_rate = ~0.9, - needs = "covid_vax_4_day" - ), - covid_vax_6_day = bn_node( - ~ runif(n = ..n, covid_vax_5_day + 300, covid_vax_5_day + 400), - missing_rate = ~0.99, - needs = "covid_vax_5_day" - ), - covid_vax_7_day = bn_node( - ~ runif(n = ..n, covid_vax_6_day + 300, covid_vax_6_day + 400), - missing_rate = ~0.99, - needs = "covid_vax_6_day" - ), - covid_vax_8_day = bn_node( - ~ runif(n = ..n, covid_vax_7_day + 300, covid_vax_7_day + 400), - missing_rate = ~0.99, - needs = "covid_vax_7_day" - ), - covid_vax_9_day = bn_node( - ~ runif(n = ..n, covid_vax_8_day + 300, covid_vax_8_day + 400), - missing_rate = ~0.99, - needs = "covid_vax_8_day" - ), - covid_vax_10_day = bn_node( - ~ runif(n = ..n, covid_vax_9_day + 300, covid_vax_9_day + 400), - missing_rate = ~0.99, - needs = "covid_vax_9_day" - ), - covid_vax_11_day = bn_node( - ~ runif(n = ..n, covid_vax_10_day + 300, covid_vax_10_day + 400), - missing_rate = ~0.99, - needs = "covid_vax_10_day" - ), - covid_vax_12_day = bn_node( - ~ runif(n = ..n, covid_vax_11_day + 300, covid_vax_11_day + 400), - missing_rate = ~0.99, - needs = "covid_vax_11_day" - ), - covid_vax_13_day = bn_node( - ~ runif(n = ..n, covid_vax_12_day + 300, covid_vax_12_day + 400), - missing_rate = ~0.99, - needs = "covid_vax_12_day" - ), - covid_vax_14_day = bn_node( - ~ runif(n = ..n, covid_vax_13_day + 300, covid_vax_13_day + 400), - missing_rate = ~0.99, - needs = "covid_vax_13_day" - ), - covid_vax_15_day = bn_node( - ~ runif(n = ..n, covid_vax_14_day + 300, covid_vax_14_day + 400), - missing_rate = ~0.99, - needs = "covid_vax_14_day" - ), - covid_vax_16_day = bn_node( - ~ runif(n = ..n, covid_vax_15_day + 300, covid_vax_15_day + 400), - missing_rate = ~0.99, - needs = "covid_vax_15_day" - ), - covid_vax_product_1 = bn_node(~ rcat(n = ..n, c("pfizer", "az"), c(0.5, 0.5)), needs = "covid_vax_1_day"), - covid_vax_product_2 = bn_node(~ if_else(runif(..n) < 0.98, covid_vax_product_1, "az"), needs = "covid_vax_2_day"), - covid_vax_product_3 = bn_node(~ rcat(n = ..n, c("pfizer", "moderna"), c(0.5, 0.5)), needs = "covid_vax_3_day"), - covid_vax_product_4 = bn_node(~ rcat(n = ..n, c("pfizer", "moderna"), c(0.5, 0.5)), needs = "covid_vax_4_day"), - covid_vax_product_5 = bn_node(~ rcat(n = ..n, c("pfizer", "moderna"), c(0.5, 0.5)), needs = "covid_vax_5_day"), - covid_vax_product_6 = bn_node(~ rcat(n = ..n, c("pfizer", "moderna"), c(0.5, 0.5)), needs = "covid_vax_6_day"), - covid_vax_product_7 = bn_node(~ rcat(n = ..n, c("pfizer", "moderna"), c(0.5, 0.5)), needs = "covid_vax_7_day"), - covid_vax_product_8 = bn_node(~ rcat(n = ..n, c("pfizer", "moderna"), c(0.5, 0.5)), needs = "covid_vax_8_day"), - covid_vax_product_9 = bn_node(~ rcat(n = ..n, c("pfizer", "moderna"), c(0.5, 0.5)), needs = "covid_vax_9_day"), - covid_vax_product_10 = bn_node(~ rcat(n = ..n, c("pfizer", "moderna"), c(0.5, 0.5)), needs = "covid_vax_10_day"), - covid_vax_product_11 = bn_node(~ rcat(n = ..n, c("pfizer", "moderna"), c(0.5, 0.5)), needs = "covid_vax_11_day"), - covid_vax_product_12 = bn_node(~ rcat(n = ..n, c("pfizer", "moderna"), c(0.5, 0.5)), needs = "covid_vax_12_day"), - covid_vax_product_13 = bn_node(~ rcat(n = ..n, c("pfizer", "moderna"), c(0.5, 0.5)), needs = "covid_vax_13_day"), - covid_vax_product_14 = bn_node(~ rcat(n = ..n, c("pfizer", "moderna"), c(0.5, 0.5)), needs = "covid_vax_14_day"), - covid_vax_product_15 = bn_node(~ rcat(n = ..n, c("pfizer", "moderna"), c(0.5, 0.5)), needs = "covid_vax_15_day"), - covid_vax_product_16 = bn_node(~ rcat(n = ..n, c("pfizer", "moderna"), c(0.5, 0.5)), needs = "covid_vax_16_day"), -) - - -# define a function to create a sublist containing all demographic / clinical variables to extract for the nth vaccination event -# function will be iterated over first second third ... Nth vaccines -sim_list_varying_i <- function(i) { - vax_variable <- glue("covid_vax_{i}_day") - - lst( - "deregistered_{i}_day" := bn_node( - ~ as.integer(runif(n = ..n, index_day, index_day + 1200)), - missing_rate = ~0.99, - needs = vax_variable - ), - "age_{i}" := bn_node(~ as.integer(rnorm( - n = ..n, mean = 60, sd = 14 - )), - needs = vax_variable), - "imd_{i}" := bn_node( - ~ as.integer(plyr::round_any(runif(n = ..n, 100, 32000), 100)), - missing = ~0.05, - needs = vax_variable - ), - "registered_{i}" := bn_node(~ rbernoulli(n = ..n, p = 0.99), - needs = vax_variable), - "stp_{i}" := bn_node(~ factor(as.integer( - runif(n = ..n, 1, 36) - ), levels = 1:36), - needs = vax_variable), - "region_{i}" := bn_node( - variable_formula = ~ rfactor( - n = ..n, - levels = c( - "North East", - "North West", - "Yorkshire and The Humber", - "East Midlands", - "West Midlands", - "East", - "London", - "South East", - "South West" - ), - p = c(0.2, 0.2, 0.3, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05) - ), - needs = vax_variable - ), - # PRIMIS - "crd_{i}" := bn_node(~ rbernoulli(n = ..n, p = 0.02), - needs = vax_variable), - # chronic respiratory disease - "chd_{i}" := bn_node(~ rbernoulli(n = ..n, p = 0.02), - needs = vax_variable), - # chronic heart disease - "ckd_{i}" := bn_node(~ rbernoulli(n = ..n, p = 0.02), - needs = vax_variable), - # chronic kidney disease - "cld_{i}" := bn_node(~ rbernoulli(n = ..n, p = 0.02), - needs = vax_variable), - # chronic liver disease - "cns_{i}" := bn_node(~ rbernoulli(n = ..n, p = 0.02), - needs = vax_variable), - # chronic neurological disease - "learndis_{i}" := bn_node(~ rbernoulli(n = ..n, p = 0.02), - needs = vax_variable), - # learning Disability - "diabetes_{i}" := bn_node(~ rbernoulli(n = ..n, p = 0.02), - needs = vax_variable), - # diabetes - "immunosuppressed_{i}" := bn_node(~ rbernoulli(n = ..n, p = 0.02), - needs = vax_variable), - # immunosuppress grouped - "asplenia_{i}" := bn_node(~ rbernoulli(n = ..n, p = 0.02), - needs = vax_variable), - # asplenia or Dysfunction of the Spleen - "severe_obesity_{i}" := bn_node(~ rbernoulli(n = ..n, p = 0.02), - needs = vax_variable), - # immunosuppress grouped - "smi_{i}" := bn_node(~ rbernoulli(n = ..n, p = 0.02), - needs = vax_variable), - # severe mental illness - "primis_atrisk_{i}" := bn_node(~ rbernoulli(n = ..n, p = 0.02), - needs = vax_variable), - # clinically vulnerable - - # extended subgroups - "rrt_{i}" := bn_node( - variable_formula = ~ rfactor( - n = ..n, - levels = c("0", - "1", - "2"), - p = c(0.98, 0.01, 0.01) - ), - needs = vax_variable - ), - "ckd_stage_3to5_{i}" := bn_node( - variable_formula = ~ rfactor( - n = ..n, - levels = c("no ckd", "3", "4", "5", "ckd, without ckd3-5 code"), - p = c(0.90, 0.06, 0.02, 0.01, 0.01) - ), - needs = vax_variable), - - # "creatinine_umol_{i}" := bn_node( - # ~ as.numeric(runif(n = ..n, 20.0, 3000.0)), - # missing_rate = ~ 0.60, - # needs = vax_variable - # ), - # "creatinine_age_{i}" := bn_node( - # ~ as.integer(rnorm(n = ..n, mean = 60, sd = 14)), - # needs = vax_variable - # ), - "copd_{i}" := bn_node(~ rbernoulli(n = ..n, p = 0.02), - needs = vax_variable), - "learndis_cat_{i}" := bn_node( - variable_formula = ~ rfactor( - n = ..n, - levels = c("No learning disability", "Down’s syndrome", "Other learning disability", "Learning disability register"), - p = c(0.80, 0.05, 0.10, 0.05)), - needs = vax_variable), - "sickle_cell_{i}" := bn_node(~ rbernoulli(n = ..n, p = 0.02), - needs = vax_variable), - "cochlear_implant_{i}" := bn_node(~ rbernoulli(n = ..n, p = 0.02), - needs = vax_variable), - - "cystic_fibrosis_{i}" := bn_node(~ rbernoulli(n = ..n, p = 0.02), - needs = vax_variable), - - "csfl_{i}" := bn_node(~ rbernoulli(n = ..n, p = 0.02), - needs = vax_variable), - ) -} - - -# compose simulation list for all variables of interest for each vaccination event -sim_list_varying <- splice( - sim_list_varying_i(1), - sim_list_varying_i(2), - sim_list_varying_i(3), - sim_list_varying_i(4), - sim_list_varying_i(5), - sim_list_varying_i(6), - sim_list_varying_i(7), - sim_list_varying_i(8), - sim_list_varying_i(9), - sim_list_varying_i(10), - sim_list_varying_i(11), - sim_list_varying_i(12), - sim_list_varying_i(13), - sim_list_varying_i(14), - sim_list_varying_i(15), - sim_list_varying_i(16) -) - - -# put vax info and other info together in the same list -sim_list <- splice(sim_list_vax_info, - sim_list_varying) - -# check and create the simulation object, including all dependencies, topological orders, etc -bn <- bn_create(sim_list, known_variables = known_variables) - -# plot the network -bn_plot(bn) - -# plot the network (connected nodes only) -bn_plot(bn, connected_only = TRUE) - -# set the seed for the simulation -set.seed(10) - -# simulate the dataset -dummydata <- - bn_simulate(bn, - pop_size = population_size, - keep_all = FALSE, - .id = "patient_id") - -# do some post simulation processing for features that are not easily handled by the simulation configuration -dummydata_processed <- dummydata %>% - # convert integer days to dates since index date and rename vars - mutate(across(ends_with("_day"), ~ as.Date(as.character(index_date + .)))) %>% - rename_with(~ str_replace(., "_day", "_date"), ends_with("_day")) %>% - # convert vaccine product short names to full product names - mutate(across( - starts_with("covid_vax_product_"), - ~ factor( - ., - levels = names(vax_product_lookup), - labels = unname(vax_product_lookup) - ) - )) - -# save the datasetin arrow format -write_feather( - dummydata_processed, - sink = here("analysis", "1-extract", "dummy-data", "dummy_varying.arrow") -) diff --git a/analysis/1-extract/variables_function.py b/analysis/1-extract/variables_function.py index a473278..7054dac 100644 --- a/analysis/1-extract/variables_function.py +++ b/analysis/1-extract/variables_function.py @@ -16,11 +16,11 @@ from ehrql.tables.tpp import ( addresses, -# opa_cost, clinical_events, practice_registrations, -# appointments, -# vaccinations + vaccinations, + ons_deaths, + apcs ) @@ -80,9 +80,6 @@ def last_prior_meds(codelist, index_date, where=True): # PRIMIS ####################################################### - - - # Asthma def has_asthma(index_date): # Asthma diagnosis @@ -481,8 +478,6 @@ def extended_subgroups(dataset, index_date, var_name_suffix=""): ## extended subgroups dataset.add_column(f"rrt_cat{var_name_suffix}", rrt_cat(index_date)) # rrt dataset.add_column(f"ckd_stage_3to5{var_name_suffix}", ckd_stage_3to5(index_date)) # ckd 3-5 - # dataset.add_column(f"creatinine_umol{var_name_suffix}", last_creatinine_event(index_date).numeric_value) - # dataset.add_column(f"creatinine_age{var_name_suffix}", patients.age_on(last_creatinine_event(index_date).date)) dataset.add_column(f"copd{var_name_suffix}", has_prior_event(codelists.copd, index_date)) # Chronic obstructive pulmonary disease dataset.add_column(f"learndis_cat{var_name_suffix}", learndis_cat(index_date)) # Learning disabilities categories dataset.add_column(f"sickle_cell{var_name_suffix}", has_prior_event(codelists.sickle_cell, index_date)) # Sickle cell anaemia @@ -491,74 +486,165 @@ def extended_subgroups(dataset, index_date, var_name_suffix=""): dataset.add_column(f"cystic_fibrosis{var_name_suffix}", has_prior_event(codelists.cystic_fibrosis, index_date)) # cystic fibrosis dataset.add_column(f"csfl{var_name_suffix}", has_prior_event(codelists.csfl, index_date)) # Cerebrospinal fluid leak dataset.add_column(f"homeless{var_name_suffix}", homeless(index_date)) # Homeless - - -# def other_cx_variables(dataset, index_date, var_name_suffix=""): - ## others of interest -# dataset.add_column(f"sol_org_trans{var_name_suffix}", has_prior_event(solid_organ_transplant, index_date)) # Organs transplant -# dataset.add_column(f"hiv{var_name_suffix}", has_prior_event(hiv_aids, index_date)) #HIV/AIDS -# dataset.add_column(f"cancer{var_name_suffix}", -# has_prior_event(cancer_nonhaem_snomed, index_date, where=clinical_events.date.is_after(index_date - days(int(3 * 365.25))))| -# has_prior_event(cancer_haem_snomed, index_date, where=clinical_events.date.is_after(index_date - days(int(3 * 365.25)))) -# ) #cancer + ############################################################ # demographic variables ############################################################ +# Ethnicity +def ethnicity5(index_date): + ethnicity5 = last_prior_event(codelists.ethnicity16, index_date).snomedct_code.to_category(codelists.ethnicity5) + return ethnicity5 + +def ethnicity16(index_date): + ethnicity16 = last_prior_event(codelists.ethnicity16, index_date).snomedct_code.to_category(codelists.ethnicity16) + return ethnicity16 + def demographic_variables(dataset, index_date, var_name_suffix=""): registration = practice_registrations.for_patient_on(index_date) - dataset.add_column(f"age{var_name_suffix}", patients.age_on(index_date)) + dataset.add_column(f"sex{var_name_suffix}", patients.sex) + dataset.add_column(f"age{var_name_suffix}", patients.age_on(index_date - days(1))) + dataset.add_column(f"ethnicity5{var_name_suffix}", ethnicity5(index_date)) + dataset.add_column(f"ethnicity16{var_name_suffix}", ethnicity16(index_date)) dataset.add_column(f"region{var_name_suffix}", registration.practice_nuts1_region_name) dataset.add_column(f"stp{var_name_suffix}", registration.practice_stp) dataset.add_column(f"imd{var_name_suffix}", addresses.for_patient_on(index_date).imd_rounded) dataset.add_column(f"carehome_status{var_name_suffix}", carehome_status(index_date)) + + +# Target disease specific functions ------------------------------- + +## Vaccination +# TODO: should we add vaccine product name? +# Identify target disease vaccinations +def target_vaccinations(target_disease_values): + + target_vaccinations = ( + vaccinations + .where(vaccinations.target_disease.is_in(target_disease_values)) + .sort_by(vaccinations.date) + ) + + return target_vaccinations + + +# Most recent vaccine before campaign +def last_vax_before_start_date(target_disease_values, cohort_start_date): + + last_vax_before_start_date = ( + target_vaccinations(target_disease_values) + .where(vaccinations.date.is_before(cohort_start_date)) + .last_for_patient() + ) + + return last_vax_before_start_date + + +# First vaccine during campaign +def first_vax_after_start_date(target_disease_values, cohort_start_date, cohort_end_date): + + first_vax_after_start_date = ( + target_vaccinations(target_disease_values) + .where(vaccinations.date.is_on_or_between(cohort_start_date, cohort_end_date)) + .first_for_patient() + ) + + return first_vax_after_start_date + + +def add_target_vaccines(dataset, target_disease_values, cohort_start_date, cohort_end_date, var_name_suffix=""): + + last_vax = last_vax_before_start_date( + target_disease_values, + cohort_start_date, + ) + + first_vax = first_vax_after_start_date( + target_disease_values, + cohort_start_date, + cohort_end_date, + ) + + dataset.add_column( + f"last_vax_date_before_start_date{var_name_suffix}", + last_vax.date, + ) + + dataset.add_column( + f"last_vax_product_before_start{var_name_suffix}", + last_vax.product_name, + ) + + dataset.add_column( + f"first_vax_date_after_start_date{var_name_suffix}", + first_vax.date, + ) + + dataset.add_column( + f"first_vax_product_after_start{var_name_suffix}", + first_vax.product_name, + ) + +# Target disease +## Mild outcome +# TODO: add +# Severe outcomes +def target_admissions(disease_icd10_codelist, cohort_start_date, cohort_end_date): + target_admissions = ( + apcs + .where(apcs.all_diagnoses.contains_any_of(disease_icd10_codelist)) + .where(apcs.admission_method.is_in([ + "21", "22", "23", "24", "25", + "2A", "2B", "2C", "2D", "28" + ])) + .where(apcs.patient_classification == "1") + .where(apcs.admission_date.is_on_or_between(cohort_start_date, cohort_end_date)) + .sort_by(apcs.admission_date) + ) + + return target_admissions -# See https://github.com/opensafely/reusable-variables/blob/main/analysis/vaccine-history/vaccine_variables.py -# this is an adpated version that only selects vaccines _near_ the index date -# without extracting the entire vacciantion history -from ehrql.tables.tpp import ( - vaccinations -) +def admitted_date(disease_icd10_codelist, cohort_start_date, cohort_end_date): + + admitted_date = ( + target_admissions(disease_icd10_codelist, cohort_start_date, cohort_end_date) + .first_for_patient() + .admission_date + ) + + return admitted_date -def add_n_vaccines(dataset, index_date, target_disease, name, direction = None, number_of_vaccines = 3, minimum_gap = 1): - assert direction in ["after", "on_or_after", "before", "on_or_before"], "direction value must be after, on_or_after, before, on_or_before" +def admitted_primary_date(disease_icd10_codelist, cohort_start_date, cohort_end_date): - assert minimum_gap > 0, "minimum_gap must be at least 1 to ensure that same-day vaccinations are not stuck in a loop" - - if direction == "after": - current_date = index_date - days(minimum_gap-1) - elif direction == "on_or_after": - current_date = index_date - days(1) - days(minimum_gap-1) - elif direction == "before": - current_date = index_date + days(minimum_gap-1) - elif direction == "on_or_before": - current_date = index_date + days(1) + days(minimum_gap-1) - else: - raise ValueError("direction must be 'before' or 'after'") + admitted_primary_date = ( + target_admissions(disease_icd10_codelist, cohort_start_date, cohort_end_date) + .where(apcs.primary_diagnosis.is_in(disease_icd10_codelist)) + .first_for_patient() + .admission_date + ) - # select all vaccination events that target {target_disease} - covid_vaccinations = ( - vaccinations - .where(vaccinations.target_disease == target_disease) - .sort_by(vaccinations.date) + return admitted_primary_date + + +def disease_death_date(disease_icd10_codelist, cohort_start_date, cohort_end_date): + + disease_death_date = case( + when( + ons_deaths.cause_of_death_is_in(disease_icd10_codelist) + & ons_deaths.date.is_on_or_between(cohort_start_date, cohort_end_date) + ).then(ons_deaths.date), + otherwise=None, ) - # loop over first, second, ..., nth vaccination event ON OR AFTER or ON OR BEFORE index date for each person - # extract info on vaccination date and product - for i in range(1, number_of_vaccines + 1): - - # vaccine variables - if direction in ["after", "on_or_after"]: - current_vax = covid_vaccinations.where(covid_vaccinations.date > (current_date + days(minimum_gap-1))).first_for_patient() - if direction in ["before", "on_or_before"]: - current_vax = covid_vaccinations.where(covid_vaccinations.date < (current_date - days(minimum_gap-1))).last_for_patient() - - dataset.add_column(f"{name}_{i}_date", current_vax.date) - dataset.add_column(f"{name}_{i}_product", current_vax.product_name) - - current_date = current_vax.date - + return disease_death_date + + +def add_target_sev_outcomes(dataset, disease_icd10_codelist, cohort_start_date, cohort_end_date, var_name_suffix=""): + + dataset.add_column(f"admitted_date{var_name_suffix}", admitted_date(disease_icd10_codelist, cohort_start_date, cohort_end_date)) + dataset.add_column(f"admitted_primary_date{var_name_suffix}", admitted_primary_date(disease_icd10_codelist, cohort_start_date, cohort_end_date)) + dataset.add_column(f"disease_death_date{var_name_suffix}", disease_death_date(disease_icd10_codelist, cohort_start_date, cohort_end_date)) diff --git a/analysis/2-prepare/prepare.R b/analysis/2-prepare/prepare.R index af6cb9f..69acef7 100644 --- a/analysis/2-prepare/prepare.R +++ b/analysis/2-prepare/prepare.R @@ -1,7 +1,6 @@ # _________________________________________________ # Purpose: -# import vaccination date data extracted by ehrql -# organise vaccination date data to "vax X product", "vax X date" (rather than "pfizer X date", "az X date", ...) +# Prepare one cohort extract for vaccination coverage and burden analyses # _________________________________________________ # Preliminaries ---- @@ -17,267 +16,338 @@ library("glue") # Import custom functions source(here("analysis", "0-lib", "design.R")) -# create output directory -output_dir <- here("output", "2-prepare", "prepare") +args <- commandArgs(trailingOnly = TRUE) + +if (length(args) == 0) { + cohort_id <- "flu_2023_24" +} else { + cohort_id <- args[[1]] +} + +cohort_info <- cohort_info |> + filter(cohort_id == .env$cohort_id) + +cohort_id_value <- cohort_info$cohort_id[[1]] +target_value <- cohort_info$target[[1]] +cohort_value <- cohort_info$cohort[[1]] +cohort_start_date_value <- cohort_info$cohort_start_date[[1]] +cohort_end_date_value <- cohort_info$cohort_end_date[[1]] +cohort_end_days_value <- cohort_info$cohort_end_days[[1]] +age_threshold_value <- cohort_info$age_threshold[[1]] +clinical_priority_value <- cohort_info$clinical_priority[[1]] + +output_dir <- here("output", "2-prepare", glue("prepare_{cohort_id}")) fs::dir_create(output_dir) -options(width = 200) # set output width for capture.output -# Import and process fixed dataset ---- +options(width = 200) -# Import fixed dataset -data_extract_fixed <- read_feather(here("output", "1-extract", "extract_fixed.arrow")) +# Import extract ---- -stopifnot( - "inconsistency between ethnicity5 and ethnicity 16" = identical(data_extract_fixed$ethnicity5, ethnicity_16_to_5(data_extract_fixed$ethnicity16)) +data_extract <- read_feather( + here("output", "1-extract", glue("extract_{cohort_id}.arrow")) ) -# print details about dataset capture.output( - skimr::skim_without_charts(data_extract_fixed), - file = fs::path(output_dir, "data_extract_fixed_skim.txt"), + skimr::skim_without_charts(data_extract), + file = fs::path(output_dir, "data_extract_skim.txt"), split = FALSE ) -# Process snapshot dataset -data_processed_fixed <- - data_extract_fixed |> +# Prepare dataset ---- + +data_prepared <- + data_extract |> lazy_dt() |> mutate( - sex = fct_case_when( - sex == "female" ~ "Female", - sex == "male" ~ "Male", - # sex == "intersex" ~ "Inter-sex", - # sex == "unknown" ~ "Unknown", - TRUE ~ NA_character_ - ), - ethnicity5 = factor(ethnicity5, levels = factor_levels$ethnicity5, ordered = FALSE), - ethnicity16 = factor(ethnicity16, levels = factor_levels$ethnicity16, ordered = FALSE) |> - fct_relabel(~ str_extract(.x, "(?<= - )(.*)")), # pick up everything after " - " - ) |> - as_tibble() - -# print details about dataset -capture.output( - skimr::skim_without_charts(data_processed_fixed), - file = fs::path(output_dir, "data_processed_fixed_skim.txt"), - split = FALSE -) - -# save processed fixed dataset -data_processed_fixed |> - select( - patient_id, - sex, - ethnicity5, - ethnicity16, - death_date, - covid_death_date - ) |> - write_feather(fs::path(output_dir, "data_fixed.arrow")) + cohort_id = cohort_id_value, + target = target_value, + cohort = cohort_value, + cohort_start_date = cohort_start_date_value, + cohort_end_date = cohort_end_date_value, + cohort_end_days = cohort_end_days_value, + age_threshold = age_threshold_value, + all = "All", + # demographics + !!!standardise_demographic_characteristics, + !!!standardise_primis_and_extended_characteristics) |> + mutate( -## delete in-memory objects to save space -rm(data_processed_fixed) -rm(data_extract_fixed) + # eligibility + age_above_eligibility_threshold = age >= age_threshold, + # used to choose if the at risk group is all clinical risk variables + # or just immunosuppressed people + clinical_priority = .data[[clinical_priority_value]], -# Import and process fixed dataset ---- + clinical_priority_only = clinical_priority & !age_above_eligibility_threshold, -# import -# data_extract_varying <- -# import_extract( -# here("lib", "dummydata", "dummyinput_varying.arrow"), -# here("output", "extracts", "extract_varying.arrow") -# ) + any_eligibility = age_above_eligibility_threshold | clinical_priority | carehome_status, -data_extract_varying <- read_feather(here("output", "1-extract", "extract_varying", "dataset.arrow")) + # baseline vaccination history + baseline_vax_status = case_when( + target == "Influenza" & + last_vax_date_before_start_date >= cohort_start_date - years(1) ~ "Vaccinated", -# Reshape vaccination data -data_vax <- - data_extract_varying |> - # lazy_dt() |> - select( - patient_id, - matches("covid_vax\\_\\d+\\_date"), - matches("covid_vax_product_\\d+"), - matches("registered_\\d+"), - matches("deregistered_\\d+"), - matches("age_\\d+"), - matches("region_\\d+"), - matches("stp_\\d+"), - matches("imd_\\d+"), - matches("imd_quintile_\\d+"), - ) |> - pivot_longer( - cols = -patient_id, - names_to = c(".value", "vax_index"), - names_pattern = "^(.*)_(\\d+)", - # values_drop_na = TRUE, # this causes an error in dtplyr - replace with filter(!is_na(covid_vax)) - # names_transform = list(vax_index = as.integer) # not supported by dtplyr - use vax_index = as.integer(vax_index) in a mutate step - ) |> - filter(!is.na(covid_vax)) |> - mutate( - vax_index = as.integer(vax_index) - ) |> - rename( - vax_date = covid_vax, - vax_product = covid_vax_product, - ) |> - # as_tibble() |> # insert this here to revert to standard dplyr as `cut` function doesn't work with dtplyr - mutate( - !!!standardise_demographic_characteristics, - vax_campaign = cut( - vax_date, - breaks = c(campaign_info$campaign_start_date, study_dates$end_date), - labels = campaign_info$campaign_label, - include.lowest = TRUE, right = FALSE - ) - ) |> - arrange(patient_id, vax_date) |> - mutate( - vax_product_raw = vax_product, - vax_product = fct_recode(factor(vax_product, vax_product_lookup), !!!vax_product_lookup) |> fct_na_value_to_level("UNMAPPED") - ) |> - group_by(patient_id) |> - mutate( - vax_interval = as.integer(vax_date - lag(vax_date, 1)) - ) |> - ungroup() + target == "Influenza" ~ "Unvaccinated", -capture.output( - skimr::skim_without_charts(data_vax), - file = fs::path(output_dir, "data_vax_skim.txt"), - split = FALSE -) + target == "COVID-19" & + last_vax_date_before_start_date >= cohort_start_date - months(6) ~ "<6 months", -# save dataset with all vaccines -write_feather(data_vax, fs::path(output_dir, "data_vax.arrow")) - -# remove vaccinations occurring within 14 days of a previous vaccination -data_vax_clean <- - # remove vaccine events occurring within 14 days of a previous vaccine event - data_vax |> - filter( - !is.na(vax_date), - is.na(vax_interval) | vax_interval >= 14, - vax_date >= study_dates$start_date, - vax_date <= study_dates$end_date - ) |> - group_by(patient_id) |> - mutate( - vax_index = row_number() - ) |> - ungroup() + target == "COVID-19" & + last_vax_date_before_start_date >= cohort_start_date - months(12) ~ "6-11 months", -capture.output( - skimr::skim_without_charts(data_vax_clean), - file = fs::path(output_dir, "data_vax_clean_skim.txt"), - split = FALSE -) + target == "COVID-19" & + !is.na(last_vax_date_before_start_date) ~ ">=12 months", -# save dataset with <14-day vaccines removed -write_feather(data_vax_clean, fs::path(output_dir, "data_vax_clean.arrow")) + target == "COVID-19" ~ "Unvaccinated", + target == "RSV" & + last_vax_date_before_start_date >= as.Date("2024-09-01") ~ "Vaccinated", -# extract event level data for vaccines ---- + target == "RSV" ~ "Unvaccinated", -data_vax_ELD0 <- read_feather(here("output", "1-extract", "extract_varying", "vaccinations.arrow")) + TRUE ~ NA_character_ + ) |> + factor(), + + # censoring + censor_date = pmin( + deregistered_date, + death_date, + cohort_end_date, + na.rm = TRUE + ), -# - remove rows where vaccination date is missing -# - attach info about the campaign during which the vaccination was given -# - collapse exact duplicates (where patient id, date, and product all match) -data_vax_ELD <- - data_vax_ELD0 |> - lazy_dt() |> - arrange(patient_id, vax_date) |> - filter(!is.na(vax_date)) |> - # distinct(.keep_all = TRUE) |> # remove exact duplicates # or use - count(patient_id, vax_date, vax_product, age) |> # or alternatively, capture how many duplicate vaccines there are. This creates a new variable `n` counting the duplicates + # follow-up vaccination + vax_time = as.integer( + pmin(first_vax_date_after_start_date, death_date, censor_date, na.rm = TRUE) - + cohort_start_date + ) + 1L, + + vax_indicator = + first_vax_date_after_start_date <= pmin(death_date, censor_date, na.rm = TRUE) & + !is.na(first_vax_date_after_start_date), + + # severe outcome: admitted + admitted_time = as.integer( + pmin(admitted_date, death_date, censor_date, na.rm = TRUE) - + cohort_start_date + ) + 1L, + + admitted_indicator = + admitted_date <= pmin(death_date, censor_date, na.rm = TRUE) & + !is.na(admitted_date), + + # severe outcome: primary diagnosis admitted + admitted_primary_time = as.integer( + pmin(admitted_primary_date, death_date, censor_date, na.rm = TRUE) - + cohort_start_date + ) + 1L, + + admitted_primary_indicator = + admitted_primary_date <= pmin(death_date, censor_date, na.rm = TRUE) & + !is.na(admitted_primary_date), + + # disease-specific death + disease_death_time = as.integer( + pmin(disease_death_date, death_date, censor_date, na.rm = TRUE) - + cohort_start_date + ) + 1L, + + disease_death_indicator = + disease_death_date <= pmin(death_date, censor_date, na.rm = TRUE) & + !is.na(disease_death_date), + + # all-cause death + death_time = as.integer( + pmin(death_date, censor_date, na.rm = TRUE) - + cohort_start_date + ) + 1L, + + death_indicator = + death_date <= censor_date & + !is.na(death_date), + + # deregistration + deregistration_time = as.integer( + pmin(deregistered_date, censor_date, na.rm = TRUE) - + cohort_start_date + ) + 1L, + + deregistration_indicator = + deregistered_date <= censor_date & + !is.na(deregistered_date), + + alive_and_registered = !death_indicator & !deregistration_indicator + ) |> as_tibble() |> mutate( - vax_product_raw = vax_product, - vax_product = fct_recode(factor(vax_product, vax_product_lookup), !!!vax_product_lookup) |> fct_na_value_to_level("UNMAPPED"), - campaign = cut( - vax_date, - breaks = c(campaign_info$campaign_start_date, as.Date(Inf)), - labels = campaign_info$campaign_label - ), - campaign_start = cut( - vax_date, - breaks = c(campaign_info$campaign_start_date, as.Date(Inf)), - labels = campaign_info$campaign_start_date - ), - ) |> - lazy_dt() + vax_status = case_when( + vax_indicator ~ "vaccinated", + death_date <= censor_date & !is.na(death_date) ~ "died", + TRUE ~ "censored" + ) |> + factor(levels = c("censored", "vaccinated", "died")), + + across( + where(is.factor) | where(is.character), + ~ fct_drop(fct_na_value_to_level(.x, level = "(Missing)")) + ) + ) +# Checks ---- -# Test equivalence of ELD extract ---- +time_check <- data_prepared |> + summarise( + min_vax_time = min(vax_time, na.rm = TRUE), + min_admitted_time = min(admitted_time, na.rm = TRUE), + min_admitted_primary_time = min(admitted_primary_time, na.rm = TRUE), + min_disease_death_time = min(disease_death_time, na.rm = TRUE), + min_death_time = min(death_time, na.rm = TRUE), + min_deregistration_time = min(deregistration_time, na.rm = TRUE) + ) -data_vax_ELD_filtered <- - data_vax_ELD |> - filter(vax_date > as.Date("1899-01-01")) |> - group_by(patient_id) |> - filter((vax_date != lag(vax_date)) | row_number() == 1) |> - mutate(vax_index = row_number()) |> - filter(vax_index <= 16) |> - ungroup() |> - as_tibble() +print(time_check) capture.output( - skimr::skim_without_charts(data_vax_ELD_filtered), - file = fs::path(output_dir, "data_vax_ELD_skim.txt"), + skimr::skim_without_charts(data_prepared), + file = fs::path(output_dir, "data_prepared_skim.txt"), split = FALSE ) -write_feather(data_vax_ELD_filtered, fs::path(output_dir, "data_vax_ELD.arrow")) - -data_vax_PLD <- - data_vax |> - select(patient_id, vax_date, vax_product = vax_product_raw, age, vax_index) -capture.output( - skimr::skim_without_charts(data_vax_PLD), - file = fs::path(output_dir, "data_vax_PLD_skim.txt"), - split = FALSE +write_feather( + data_prepared, + fs::path(output_dir, glue("prepare_{cohort_id}.arrow")) ) -write_feather(data_vax_PLD, fs::path(output_dir, "data_vax_PLD.arrow")) -# check equality of datasets -cat( - "\n", - "are datasets from ELD versus PLD identical after some standardisation? \n" -) +vax_prod_table <- bind_rows( + data_prepared |> + transmute( + indicador = "last_vax_before_start", + date = last_vax_date_before_start_date, + product = last_vax_product_before_start + ), -all.equal(data_vax_ELD_filtered, data_vax_PLD) - -# report multiple vaccinations on the same day -cat( - "\n", - "number of occassions where a person is vaccinated more than once in a day:\n", - data_vax_ELD |> - group_by(patient_id, vax_date) |> - summarise(n = n()) |> - filter(n > 1) |> - nrow() -) + data_prepared |> + transmute( + indicador = "first_vax_after_start", + date = first_vax_date_after_start_date, + product = first_vax_product_after_start + ) +) |> + filter(!is.na(date)) |> + mutate( + year = year(date), + period = case_when( + month(date) %in% 9:12 ~ paste0("Sep-Feb/", year(date)), + month(date) %in% 1:2 ~ paste0("Sep-Feb/", year(date) - 1), + month(date) %in% 3:8 ~ paste0("Feb-Aug/", year(date)), + TRUE ~ NA_character_ + ) + ) |> + count(indicador, period, product, name = "n_round10") |> + mutate( + n_round10 = round_any(n_round10, sdc_threshold + )) |> + arrange(indicador, period, product) -# report no vax date -cat( - "\n", - "number of occassions where a person is vaccinated with a null date:\n", - data_vax_ELD |> - as_tibble() |> - nrow() +write_csv( + vax_prod_table, + fs::path(output_dir, glue("vax_prod_table_{cohort_id}.csv")) ) - -# report no vax date -cat( - "\n", - "number of occassions where a person is vaccinated on or before 1899:\n", - data_vax_ELD |> - filter(vax_date <= as.Date("1899-01-01")) |> - as_tibble() |> - nrow() -) +# Table 1 ---- +# This information is inside each adjusted_estimates() +# table1_summary <- function(...) { +# group_names <- c(...) + +# summary_table <- +# data_prepared |> +# group_by(across(all_of(group_names))) |> +# lazy_dt() |> +# summarise( +# target = first(target), +# cohort_id = first(cohort_id), +# cohort = first(cohort), + +# number_total_subgroup = n(), + +# number_baseline_vax_status_vaccinated = +# sum(baseline_vax_status == "Vaccinated", na.rm = TRUE), + +# number_baseline_vax_status_unvaccinated = +# sum(baseline_vax_status == "Unvaccinated", na.rm = TRUE), + +# number_baseline_vax_status_less_than_6_months = +# sum(baseline_vax_status == "<6 months", na.rm = TRUE), + +# number_baseline_vax_status_6_11_months = +# sum(baseline_vax_status == "6-11 months", na.rm = TRUE), + +# number_baseline_vax_status_12_months_or_more = +# sum(baseline_vax_status == ">=12 months", na.rm = TRUE), + +# number_vax_post_index_date = +# sum(vax_indicator, na.rm = TRUE), + +# number_admitted = +# sum(admitted_indicator, na.rm = TRUE), + +# number_admitted_primary = +# sum(admitted_primary_indicator, na.rm = TRUE), + +# number_disease_death = +# sum(disease_death_indicator, na.rm = TRUE), + +# number_all_cause_death = +# sum(death_indicator, na.rm = TRUE), + +# number_deregistered = +# sum(deregistration_indicator, na.rm = TRUE), + +# .groups = "drop" +# ) |> +# as_tibble() + +# return(summary_table) +# } + +# table1 <- +# level_combos |> +# mutate( +# table1 = map2( +# group1, group2, +# .f = function(x, y) { +# if (is.na(y)) y <- NULL + +# lookup <- c( +# group1_value = x, +# group2_value = y +# ) + +# table1_summary(x, y) |> +# mutate(across(c(all_of(c(x, y))), as.character)) |> +# rename(any_of(lookup)) +# } +# ) +# ) |> +# unnest(table1) |> +# select( +# target, +# cohort_id, +# cohort, +# group1, +# group1_value, +# group2, +# group2_value, +# everything() +# ) + +# Save ---- +# write_csv( +# table1, +# fs::path(output_dir, glue("table1_{cohort_id}.csv")) +# ) diff --git a/analysis/2-prepare/vax_data_quality.R b/analysis/2-prepare/vax_data_quality.R deleted file mode 100644 index 688c13d..0000000 --- a/analysis/2-prepare/vax_data_quality.R +++ /dev/null @@ -1,186 +0,0 @@ -# _________________________________________________ -# Purpose: -# import vaccination date data extracted by ehrql -# report data quality for event-level vaccination data -# _________________________________________________ - -# Preliminaries ---- - -# Import libraries -library("tidyverse") -library("dtplyr") -library("lubridate") -library("arrow") -library("here") -library("glue") - -# Import custom functions -source(here("analysis", "0-lib", "design.R")) - -# create output directory -output_dir <- here("output", "2-prepare", "vax_data_quality") -fs::dir_create(output_dir) -options(width = 200) # set output width for capture.output - - -# extract event level data for vaccines ---- - -data_vax_ELD0 <- read_feather(here("output", "1-extract", "extract_varying", "vaccinations.arrow")) - - -# when running locally, create some dummy data to make outputs more sensible -# this should NOT EVER be run on the real data, as it will overwrite real data -if (localrun) { - data_vax_ELD0 <- - data_vax_ELD0 |> - mutate( - vax_date = as.Date(runif(n(), as.Date("2019-01-12"), max(campaign_info$early_milestone_date))), - vax_product = sample(vax_product_lookup, n(), replace = TRUE) - ) -} - -# - remove rows where vaccination date is missing -# - attach info about the campaign during which the vaccination was given -# - collapse exact duplicates (where patient id, date, and product all match) -data_vax_ELD <- - data_vax_ELD0 |> - lazy_dt() |> - arrange(patient_id, vax_date) |> - filter(!is.na(vax_date)) |> - # distinct(.keep_all = TRUE) |> # remove exact duplicates # or use - count(patient_id, vax_date, vax_product, age) |> # or alternatively, capture how many duplicate vaccines there are. This creates a new variable `n` counting the duplicates - as_tibble() |> - mutate( - vax_product_raw = vax_product, - vax_product = fct_recode(factor(vax_product, vax_product_lookup), !!!vax_product_lookup) |> fct_na_value_to_level("UNMAPPED"), - campaign = cut( - vax_date, - breaks = c(campaign_info$campaign_start_date, as.Date(Inf)), - labels = campaign_info$campaign_label - ), - campaign_start = cut( - vax_date, - breaks = c(campaign_info$campaign_start_date, as.Date(Inf)), - labels = campaign_info$campaign_start_date - ), - ) |> - lazy_dt() - -# report any unmapped product names -# and stop if there are any -unmapped_products <- data_vax_ELD |> filter(vax_product %in% "UNMAPPED") |> pull(vax_product_raw) |> unique() -cat("Unmapped product names: \n") -cat(paste0(unmapped_products, collapse = "\n")) -stopifnot("There are unmapped product names" = length(unmapped_products) == 0) - -## count products ---- - -# count the number of vaccines of each type over all time, by age (over / under 12 years) - -# overall -count_product <- - data_vax_ELD |> - mutate( - under12 = age < 12, - vax_date_onorafter20201201 = if_else(vax_date >= as.Date("2020-12-01"), vax_date, as.Date(NA)) - ) |> - group_by(under12, vax_product) |> - summarise( - count_total = roundmid_any(n(), sdc_threshold), - count_before20200101 = roundmid_any(sum(vax_date < as.Date("2020-01-01")), sdc_threshold), - count_onorafter20200101 = roundmid_any(sum(vax_date >= as.Date("2020-01-01")), sdc_threshold), - count_onorafter20201201 = roundmid_any(sum(vax_date >= as.Date("2020-12-01")), sdc_threshold), - earliest_date_onorafter20201201 = min(vax_date_onorafter20201201, na.rm = TRUE), - count_on_earlier_date = sum(vax_date_onorafter20201201 %in% min(vax_date_onorafter20201201, na.rm = TRUE)) - ) |> - as_tibble() - -write_csv(count_product, fs::path(output_dir, "count_product.csv")) - -# count the number of vaccines of each type, by age (over or under 12 years) and by campaign - -count_product_campaign <- - data_vax_ELD |> - mutate( - under12 = age < 12 - ) |> - group_by(under12, campaign, vax_product) |> - summarise( - count_total = roundmid_any(n(), sdc_threshold), - earliest_date = min(vax_date), - count_on_earliest_date = sum(vax_date %in% min(vax_date)) - ) |> - as_tibble() - -write_csv(count_product_campaign, fs::path(output_dir, "count_product_campaign.csv")) - - - -## count product same-day co-occurrence ---- - -# summary function to turn eg c("A", "A", "B", "C", "C", "C") into "2x A -- 1x B -- 3x C" -flat_table_chr <- function(x, collapse = NULL) { - rle_obj <- rle(x) - paste0(rle_obj$length, "x ", rle_obj$values, collapse = collapse) -} - -products_cooccurrence_flat <- - data_vax_ELD |> - filter(age >= 12) |> - arrange(vax_product) |> - group_by(patient_id, vax_date, campaign) |> - summarise( - # vax_product = flat_table_chr(as.character(vax_product), collapse = " --AND-- "), # this is only useful if we don't remove duplicate rows (same patient, same date, same product) - vax_product = paste0(paste0(n, "x ", as.character(vax_product)), collapse = " --AND-- "), - ) - - -# CONSIDER USING THIS FORMAT INSTEAD, IF EASIER! - -# products_cooccurrence_wide <- -# data_vax_ELD |> -# filter(age >= 12) |> -# pivot_wider( -# id_cols = c(patient_id, vax_date, campaign), -# names_from = vax_product, -# values_from = n, -# values_fill = 0L -# ) |> as_tibble() - - -# count overall - -count_products_cooccurrence <- - products_cooccurrence_flat |> - mutate( - vax_date_onorafter20201201 = if_else(vax_date >= as.Date("2020-12-01"), vax_date, as.Date(NA)) - ) |> - group_by(vax_product) |> - summarise( - count_total = roundmid_any(n(), sdc_threshold), - count_before20200101 = roundmid_any(sum(vax_date < as.Date("2020-01-01")), sdc_threshold), - count_onorafter20200101 = roundmid_any(sum(vax_date >= as.Date("2020-01-01")), sdc_threshold), - count_onorafter20201201 = roundmid_any(sum(vax_date >= as.Date("2020-12-01")), sdc_threshold), - earliest_date_onorafter20201201 = min(vax_date_onorafter20201201, na.rm = TRUE), - count_on_earliest_date = sum(vax_date_onorafter20201201 %in% min(vax_date_onorafter20201201, na.rm = TRUE)) - ) |> - as_tibble() - -write_csv(count_products_cooccurrence, fs::path(output_dir, "count_product_cooccurrence.csv")) - -# count by campaign - -count_products_cooccurrence_campaign <- - products_cooccurrence_flat |> - group_by(vax_product, campaign) |> - summarise( - count_total = roundmid_any(n(), sdc_threshold), - count_before20200101 = roundmid_any(sum(vax_date < as.Date("2020-01-01")), sdc_threshold), - count_onorafter20200101 = roundmid_any(sum(vax_date >= as.Date("2020-01-01")), sdc_threshold), - count_onorafter20201201 = roundmid_any(sum(vax_date >= as.Date("2020-12-01")), sdc_threshold), - earliest_date = min(vax_date, na.rm = TRUE), - count_on_earliest_date = sum(vax_date %in% min(vax_date, na.rm = TRUE)) - ) |> - as_tibble() - -write_csv(count_products_cooccurrence_campaign, fs::path(output_dir, "count_product_cooccurrence_campaign.csv")) diff --git a/analysis/3-history/report_history.R b/analysis/3-history/report_history.R deleted file mode 100644 index db9e370..0000000 --- a/analysis/3-history/report_history.R +++ /dev/null @@ -1,398 +0,0 @@ -# _______________________________________________________________________________________ -# Purpose: -# Report the distribution of vaccines schedules in different population subgroups -# _______________________________________________________________________________________ - -# Preliminaries ---- - -# Import libraries -library("tidyverse") -library("dtplyr") -library("lubridate") -library("glue") -library("here") -library("arrow") - -# Import custom functions -source(here("analysis", "0-lib", "design.R")) - -## Create output directory -output_dir <- here("output", "3-history", "report_history") -fs::dir_create(output_dir) - -# Import prepared data ---- -data_fixed <- read_feather(here("output", "2-prepare", "prepare", "data_fixed.arrow")) -data_varying <- read_feather(here("output", "2-prepare", "prepare", "data_vax.arrow")) -data_varying_clean <- read_feather(here("output", "2-prepare", "prepare", "data_vax_clean.arrow")) - -# _______________________________________________________________________________________ -# Process datasets ---- -# _______________________________________________________________________________________ - -## add time-invariant info from the fixed dataset to vaccine data and format variables for printing / plots ---- - -data_vax <- - left_join( - lazy_dt(data_varying), - lazy_dt(data_fixed) |> select(patient_id, sex, ethnicity5, ethnicity16, death_date), - by = "patient_id" - ) |> - mutate( - vax_dosenumber = factor(vax_index, levels = sort(unique(vax_index)), labels = paste("Dose", sort(unique(vax_index)))), - vax_week = floor_date(vax_date, unit = "week", week_start = 1), - vax_product_core = fct_other(vax_product, vax_product_core_levels, other_level = "other"), - all = "" - ) |> - as_tibble() - -data_vax_clean <- - left_join( - lazy_dt(data_varying_clean), - lazy_dt(data_fixed) |> select(patient_id, sex, ethnicity5, ethnicity16, death_date), - by = "patient_id" - ) |> - mutate( - vax_dosenumber = factor(vax_index, levels = sort(unique(vax_index)), labels = paste("Dose", sort(unique(vax_index)))), - vax_week = floor_date(vax_date, unit = "week", week_start = 1), - vax_product_core = fct_other(vax_product, vax_product_core_levels, other_level = "other"), - all = "", - all2 = "" - ) |> - as_tibble() |> - mutate( - across(where(is.factor) | where(is.character), ~ fct_explicit_na(.x, na_level = "Unknown")) - ) - - -# _______________________________________________________________________________________ -# data validation checks ---- -# _______________________________________________________________________________________ - -## output vax date validation info ---- - -summary_validation <- - data_vax |> - # group_by(vax_campaign) |> - summarise( - n = roundmid_any(n(), sdc_threshold), - n_missing_date = roundmid_any(sum(is.na(vax_date)), sdc_threshold), - pct_missing_date = n_missing_date / n, - n_earlier_than_start_date = roundmid_any(sum(vax_date < study_dates$start_date, na.rm = TRUE), sdc_threshold), - pct_earlier_than_start_date = n_earlier_than_start_date / n, - n_earlier_than_firstpossiblevax_date = roundmid_any(sum(vax_date < study_dates$firstpossiblevax_date, na.rm = TRUE), sdc_threshold), - pct_earlier_than_firstpossiblevax_date = n_earlier_than_firstpossiblevax_date / n, - n_interval_within_14days = roundmid_any(sum(vax_interval < 14, na.rm = TRUE), sdc_threshold), - pct_interval_within_14days = n_interval_within_14days / n, - ) |> - ungroup() - -write_csv(summary_validation, fs::path(output_dir, "validation.csv")) - - -## output vax date validation info, stratified by dose number and product ---- - -summary_validation_stratified <- - data_vax |> - group_by( - vax_dosenumber, vax_product_core - ) |> - summarise( - n = roundmid_any(n(), sdc_threshold), - n_missing_date = roundmid_any(sum(is.na(vax_date)), sdc_threshold), - pct_missing_date = n_missing_date / n, - n_earlier_than_start_date = roundmid_any(sum(vax_date < study_dates$start_date, na.rm = TRUE), sdc_threshold), - pct_earlier_than_start_date = n_earlier_than_start_date / n, - n_earlier_than_firstpossiblevax_date = roundmid_any(sum(vax_date < study_dates$firstpossiblevax_date, na.rm = TRUE), sdc_threshold), - pct_earlier_than_firstpossiblevax_date = n_earlier_than_firstpossiblevax_date / n, - n_interval_within_14days = roundmid_any(sum(vax_interval < 14, na.rm = TRUE), sdc_threshold), - pct_interval_within_14days = n_interval_within_14days / n, - ) |> - ungroup() - -write_csv(summary_validation_stratified, fs::path(output_dir, "validation_stratified.csv")) - - -## output frequency of total number of doses (vax_count) - -summary_vax_count <- - data_vax |> - group_by(patient_id) |> - summarise(vax_count = n()) |> - group_by(vax_count) |> - summarise(frequency = roundmid_any(n(), sdc_threshold)) - -write_csv(summary_vax_count, fs::path(output_dir, "validation_vax_count.csv")) - - -## output frequency of vaccination product by campaign - -summary_vax_product_campaign <- - data_vax |> - group_by( - vax_campaign, vax_product - ) |> - summarise( - n = roundmid_any(n(), sdc_threshold) - ) |> - ungroup() - -write_csv(summary_vax_product_campaign, fs::path(output_dir, "vax_counts_product_campaign.csv")) - - -## output frequency of vaccination product by dose number - -summary_vax_product_dosenumber <- - data_vax |> - group_by( - vax_dosenumber, vax_product - ) |> - summarise( - n = roundmid_any(n(), sdc_threshold) - ) |> - ungroup() - -write_csv(summary_vax_product_dosenumber, fs::path(output_dir, "vax_counts_product_dosenumber.csv")) - - - -# _______________________________________________________________________________________ -# Report info using characteristics recorded on each vaccination date ---- -# _______________________________________________________________________________________ - -## note that all patient characteristics are determined as at the date of vaccination. -## for example, a person who moves from london to manchester between their first and second dose will be classed as in "london" for their first dose and "north west" for their second dose. -## except for fixed vairables extracted from the dataset_definition_fixed script - -## output fully stratified vaccine counts ---- -## this is useful for anyone wanting to externally re-construct different cuts of data for plotting etc -## However, the row count current exceeds the limit for viewing outputs in Airlock, so we also break these down in uni / bivariate tables below - -summary_stratified <- - data_vax |> - group_by( - vax_dosenumber, vax_product_core, vax_campaign, - sex, ageband4, ethnicity5, region, imd_quintile, - # PRIMIS - # crd, chd, ckd, cld, cns, learndis, diabetes, immunosuppressed, asplenia, severe_obesity, smi, - # primis_atrisk - # Extended subgroups - # ckd_rrt, copd, learndis_cat, sickle_cell, cirrhosis, cochlear_implant, cystic_fibrosis, csfl - - ) |> - summarise( - n = roundmid_any(n(), sdc_threshold) - ) |> - ungroup() - -write_csv(summary_stratified, fs::path(output_dir, "vax_counts_stratified.csv")) - - - -## output plots of vaccine counts per week --- -## stratified by product, dose number, and other characteristics - -plot_vax_dates <- function(rows, cols) { - summary_by <- data_vax_clean |> - group_by(vax_product_core, vax_week) |> - group_by({{ rows }}, {{ cols }}, .add = TRUE) |> - summarise( - n = roundmid_any(n(), sdc_threshold) - ) - - temp_plot <- - ggplot(summary_by) + - geom_col( - aes(x = vax_week, y = n, fill = vax_product_core, group = vax_product_core), - alpha = 0.5, - position = position_stack(reverse = TRUE), - # position=position_identity(), - width = 7 - ) + - facet_grid( - rows = vars({{ rows }}), - cols = vars({{ cols }}), - switch = "y", - space = "free_x", - scales = "free_x" - ) + - labs( - x = "Date", - y = NULL, - fill = NULL - ) + - scale_fill_brewer( - palette = "Set2", - na.value = "grey50", - labels = function(breaks) { - breaks[is.na(breaks)] <- "Other" - breaks - } - ) + - # scale_fill_manual(values=c(RColorBrewer::brewer.pal(8, "Set2"), "grey50"))+ - scale_x_date( - breaks = as.Date(c("2021-01-01", "2022-01-01", "2023-01-01", "2024-01-01")), - date_minor_breaks = "month", - date_labels = "%Y", # labels = scales::label_date("%b"), - # sec.axis = sec_axis( - # breaks=as.Date(c("2021-01-01","2022-01-01","2023-01-01","2024-01-01")), - # transform = ~as.Date(.), - # labels = scales::label_date("%Y") - # ) - ) + - theme_minimal() + - theme( - axis.text.x.top = element_text(hjust = 0), - axis.text.x.bottom = element_text(hjust = 0), - strip.text.y.left = element_text(angle = 0, hjust = 1), - strip.placement = "outside", - axis.ticks.x = element_line(), - legend.position = "bottom" - ) - - # print(temp_plot) - - row_name <- deparse(substitute(rows)) - col_name <- deparse(substitute(cols)) - - ggsave(fs::path(output_dir, glue("vax_dates_{row_name}_{col_name}.png")), plot = temp_plot) -} - - - -plot_vax_dates(all, all2) - -plot_vax_dates(sex, all) -plot_vax_dates(ageband4, all) -plot_vax_dates(ethnicity5, all) -plot_vax_dates(region, all) -plot_vax_dates(imd_quintile, all) -plot_vax_dates(vax_campaign, all) -plot_vax_dates(vax_dosenumber, all) - -# PRIMIS -# plot_vax_dates(crd, all) # chronic respiratory disease -# plot_vax_dates(chd, all) # chronic heart disease -# plot_vax_dates(ckd, all) # chronic kidney disease -# plot_vax_dates(cld, all) # chronic liver disease -# plot_vax_dates(cns, all) #chronic neurological -# plot_vax_dates(learndis, all) # learning disability -# plot_vax_dates(diabetes, all) # diabetes -# plot_vax_dates(immunosuppressed, all) # immunosuppressed -# plot_vax_dates(asplenia, all) # asplenia or dysfunction of the spleen -# plot_vax_dates(severe_obesity, all) # obesity -# plot_vax_dates(smi, all) # severe mental illness -# plot_vax_dates(primis_atrisk, all) # clinically vulnerable - - -plot_vax_dates(sex, vax_dosenumber) -plot_vax_dates(ageband4, vax_dosenumber) -plot_vax_dates(ethnicity5, vax_dosenumber) -plot_vax_dates(region, vax_dosenumber) -plot_vax_dates(imd_quintile, vax_dosenumber) -plot_vax_dates(vax_campaign, vax_dosenumber) -# PRIMIS -# plot_vax_dates(crd, vax_dosenumber) # chronic respiratory disease -# plot_vax_dates(chd, vax_dosenumber) # chronic heart disease -# plot_vax_dates(ckd, vax_dosenumber) # chronic kidney disease -# plot_vax_dates(cld, vax_dosenumber) # chronic liver disease -# plot_vax_dates(cns, vax_dosenumber) # chronic neurological disease -# plot_vax_dates(learndis, vax_dosenumber) # learning Disability -# plot_vax_dates(diabetes, vax_dosenumber) # diabetes -# plot_vax_dates(immunosuppressed, vax_dosenumber) # immunosuppressed -# plot_vax_dates(asplenia, vax_dosenumber) # asplenia or dysfunction of the spleen -# plot_vax_dates(severe_obesity, vax_dosenumber) # obesity -# plot_vax_dates(smi, vax_dosenumber) # severe mental illness -# plot_vax_dates(primis_atrisk, vax_dosenumber) # Clinically vulnerable - -## output plots of time since previous vaccination by product, dose number, and other characteristics ---- - -plot_vax_intervals <- function(rows, cols) { - summary_by <- data_vax_clean |> - filter(vax_index != 1) |> - mutate( - vax_interval = roundmid_any(vax_interval + 1, 7), # to split into 0-6, 7-13, 14-20, 21-28, .... - vax_dosenumber = factor(vax_index, levels = sort(unique(vax_index)), labels = paste("Dose ", sort(unique(vax_index)) - 1, "-", sort(unique(vax_index)))), - ) |> - group_by(vax_dosenumber, vax_product_core, vax_interval) |> - group_by({{ rows }}, {{ cols }}, .add = TRUE) |> - summarise( - n = roundmid_any(n(), sdc_threshold), - ) - - temp_plot <- - ggplot(summary_by) + - geom_col( - aes(x = vax_interval, y = n, fill = vax_product_core, group = vax_product_core), - alpha = 0.5, - position = position_stack(reverse = TRUE), - # position=position_identity(), - width = 7 - ) + - facet_grid( - rows = vars({{ rows }}), - cols = vars({{ cols }}), - switch = "y", - space = "free_x", - scales = "free_x" - ) + - labs( - x = "Interval (days)", - y = NULL, - fill = NULL - ) + - scale_fill_brewer( - palette = "Set2", - na.value = "grey50", - labels = function(breaks) { - breaks[is.na(breaks)] <- "Other" - breaks - } - ) + - scale_x_continuous( - breaks = (0:100) * 4 * 7, - # limits = c(0, NA), - sec.axis = sec_axis( - transform = ~ . / 7 - ) - ) + - # scale_y_continuous(limits=c(0,100))+ - theme_minimal() + - theme( - strip.text.y.left = element_text(angle = 0, hjust = 1), - strip.placement = "outside", - axis.ticks.x = element_line(), - # axis.text.y = element_blank(), - legend.position = "bottom" - ) - - print(temp_plot) - - row_name <- deparse(substitute(rows)) - col_name <- deparse(substitute(cols)) - - ggsave(fs::path(output_dir, glue("vax_intervals_{row_name}_{col_name}.png")), plot = temp_plot) -} - -plot_vax_intervals(sex, vax_dosenumber) -plot_vax_intervals(ageband4, vax_dosenumber) -plot_vax_intervals(ethnicity5, vax_dosenumber) -plot_vax_intervals(region, vax_dosenumber) -plot_vax_intervals(imd_quintile, vax_dosenumber) -plot_vax_intervals(vax_campaign, vax_dosenumber) -plot_vax_intervals(vax_campaign, all) -plot_vax_intervals(vax_dosenumber, all) - -# PRIMIS - -# plot_vax_intervals(crd, vax_dosenumber) # chronic respiratory disease -# plot_vax_intervals(chd, vax_dosenumber) # chronic heart disease -# plot_vax_intervals(ckd, vax_dosenumber) # chronic kidney disease -# plot_vax_intervals(cld, vax_dosenumber) # chronic liver disease -# plot_vax_intervals(cns, vax_dosenumber) # chronic neurological -# plot_vax_intervals(learndis, vax_dosenumber) # learning disability -# plot_vax_intervals(diabetes, vax_dosenumber) # diabetes -# plot_vax_intervals(immunosuppressed, vax_dosenumber) # immunosuppressed -# plot_vax_intervals(asplenia, vax_dosenumber) # asplenia or dysfunction of the spleen -# plot_vax_intervals(severe_obesity, vax_dosenumber) # obesity -# plot_vax_intervals(smi, vax_dosenumber) # severe mental illness -# plot_vax_intervals(primis_atrisk, vax_dosenumber) # Clinically vulnerable diff --git a/analysis/3-report/report_cohort.R b/analysis/3-report/report_cohort.R new file mode 100644 index 0000000..2425696 --- /dev/null +++ b/analysis/3-report/report_cohort.R @@ -0,0 +1,752 @@ +# _______________________________________________________________________________________ +# Purpose: +# Report vaccination coverage and disease burden in different population subgroups for each campaign +# _______________________________________________________________________________________ + + +# Preliminaries ---- + +# Import libraries +library("tidyverse") +library("dtplyr") +library("lubridate") +library("glue") +library("here") +library("arrow") +library("survival") +library("splines") +library("parglm") + +# Import custom functions +source(here("analysis", "0-lib", "design.R")) + +args <- commandArgs(trailingOnly = TRUE) + + +if (length(args) == 0) { + cohort_id <- "flu_2023_24" +} else { + cohort_id <- args[[1]] +} + +cohort_info <- cohort_info |> + filter(cohort_id == .env$cohort_id) +# list2env(campaign_info, globalenv()) + +# how wide are the temporal bins for frequencies over time for Kaplan-Meier plots? in days +temporal_resolution_km <- 7L + +# overwrite cohort end to match rounded Kaplan-Meier curve +cohort_info$cohort_end_days <- ceiling_any(cohort_info$cohort_end_days, temporal_resolution_km) + +# Create output directory ---- + +output_dir <- here("output", "3-report", glue("report_{cohort_id}")) +fs::dir_create(output_dir) +options(width = 200) + +# Import prepared data ---- + +data_prepared <- read_feather( + here("output", "2-prepare", glue("prepare_{cohort_id}"), glue("prepare_{cohort_id}.arrow")) +) + +capture.output( + skimr::skim_without_charts(data_prepared), + file = fs::path(output_dir, "data_prepared_skim.txt"), + split = FALSE +) + +################################################################################################################# +# 1. Baseline characteristics +################################################################################################################# + +# ________________________________ +# 1.A Baseline vaccination history +# ________________________________ + +baseline_vax_summary <- function(data, ...) { + group_names <- c(...) + + data |> + count( + across(all_of(group_names)), + cohort_id, + baseline_vax_status, + name = "n_round10" + ) |> + group_by( + across(all_of(group_names)), + cohort_id + ) |> + mutate( + n_total_round10 = sum(n_round10), + n_total_round10 = roundmid_any(n_total_round10, sdc_threshold), + n_round10 = roundmid_any(n_round10, sdc_threshold), + pct = round(100 * n_round10 / n_total_round10, 1) + ) |> + ungroup() +} + +baseline_vax_table <- + level_combos |> + mutate( + summary = map2( + group1, + group2, + \(x, y) { + + if (is.na(y)) y <- NULL + + lookup <- c( + group1_value = x, + group2_value = y + ) + + baseline_vax_summary(data_prepared, x, y) |> + mutate(across(c(all_of(c(x, y))), as.character)) |> + rename(any_of(lookup)) + } + ) + ) |> + unnest(summary) |> + select( + cohort_id, + group1, + group1_value, + group2, + group2_value, + everything() + ) + +write_csv( + baseline_vax_table, + fs::path(output_dir, "baseline_vax_history.csv") +) + +################################################################################################################# +# 2. Follow-up descriptive analyses +################################################################################################################# + +# ________________________________ +# 2.A Vaccination uptake +# Kaplan-Meier estimates of cumulative vaccination coverage +# ________________________________ + +# This code borrows heavily from the KM reusable action https://github.com/opensafely-actions/kaplan-meier-function/blob/main/analysis/km.R +# The reason not to use the KM reusable action directly is that we need to reuse it multiple times across many different stratification variables +# This would create a large number of project.yaml actions, so it's easier to do the repeats within a single script +# The reusable action could be modified to enable multiple stratified analyses to be run, but it's a bit faffy and a bit scope-creepy + + +## tests ---- + +times_count <- table(cut(data_prepared$vax_time, c(-Inf, 0, 1, Inf), right = FALSE, labels = c("<0", "0", ">0")), useNA = "ifany") + +if (!identical(as.integer(times_count), c(0L, 0L, nrow(data_prepared)))) { + print(times_count) + stop("all event times must be strictly positive") +} + +## Function to calculate KM estimates for a given stratification variable ---- + +# group variables are provided as characters via dots (...) +# resolution argument is the precision used for the time dimension. If zero, then original resolution is used. +km_estimates <- function(data, group_name1, group_name2, event_name, event_time, event_indicator, resolution = 0) { + + group_names <- c(group_name1, group_name2) + + # if (is.na(group_name2)) { + # group_name2 <- "all" + # } + + data_outcome <- + data |> + select( + patient_id, + all_of(group_names), + event_time = any_of(event_time), + event_indicator = any_of(event_indicator) + ) |> + mutate(event_time = ceiling_any(event_time, resolution)) + + data_km <- + data_outcome |> + group_by(across(all_of(group_names))) |> + nest() |> + mutate( + surv_obj_tidy = map( + data, ~ { + survfit( + Surv(event_time, event_indicator) ~ 1, + data = .x, + conf.type = "log-log" + ) |> + broom::tidy() |> + add_row( + time = 0, # assumes time origin is zero + n.risk = 0, + n.event = 0, + n.censor = 0, + estimate = 1, + conf.low = 1, + conf.high = 1, + .before = 1L + ) |> + complete( + time = seq(0L, cohort_info$cohort_end_days, resolution), # fill in 1 row for each period (defined by resolution) of follow up + fill = list(n.event = 0L, n.censor = 0L) # fill in zero events on those days + ) |> + fill( + n.risk, + .direction = "up" + ) |> + fill( + estimate, conf.low, conf.high, + .direction = "down" + ) + } + ), + ) |> + select(-data) |> + unnest(surv_obj_tidy) |> + mutate( + + # disclosure control + n.risk = ceiling_any(n.risk, sdc_threshold), + estimate = plyr::round_any(estimate, sdc_threshold / nth(n.risk, 2)), # use 2nd value as this skips the t=0 row where n.risk=0 + conf.low = plyr::round_any(conf.low, sdc_threshold / nth(n.risk, 2)), + conf.high = plyr::round_any(conf.high, sdc_threshold / nth(n.risk, 2)), + + # cumulative incidence + cmlinc = 1 - estimate, + cmlinc.low = 1 - conf.high, + cmlinc.high = 1 - conf.low, + ) + + rm(data_outcome) + + # print(data_km) + + # write tables that capture underlying plotting data + data_km_nozero <- + data_km |> + arrange(across(all_of(group_names))) |> + ungroup() |> + filter(time != 0) |> + select( + all_of(group_names), + time, + cmlinc, + cmlinc.low, + cmlinc.high, + ) + + # commented out because outputting the entire data creates a dataset that is too large to output-check on it's own + # write_csv(data_km_nozero, fs::path(output_dir, glue("km_{event_name}_{paste0(group_names, collapse='_')}.csv"))) + + return(data_km_nozero) + +} + +# for testing function interactively +# km_estimates_vax <- partial( +# km_estimates, data = data_prepared, event_name = "vax", event_time = "vax_time", event_indicator = "vax_indicator", resolution = temporal_resolution_km +# ) +# km_estimates_vax("all", "all") +# km_estimates_vax("all", "ageband4") +# km_estimates_vax("ageband4", "all") +# km_estimates_vax("ageband4", "sex") + +## Function to get km estimates over all subgroup combinations ---- + +get_all_km_estimates <- function(data, event_name, event_time, event_indicator, resolution) { + # loop over all group1 and group2 variable combinations and combine into one big dataset + km_estimates_table <- + level_combos |> + mutate( + km_summary = map2( + group1, group2, + .f = \(x, y) { + + if (is.na(y)) y <- NULL + lookup <- c(group1_value = x, group2_value = y) + + km_estimates(data = data, group_name1 = x, group_name2 = y, event_name = event_name, event_time = event_time, event_indicator = event_indicator, resolution = temporal_resolution_km) |> + mutate(across(c(all_of(c(x, y))), as.character)) |> + rename(all_of(c(group1_value = x, group2_value = y))) + } + ) + ) |> + unnest(km_summary) |> + select(group1, group1_value, group2, group2_value, everything()) |> + mutate( + cohort_end_milestone = (time == cohort_info$cohort_end_days) * 1L, + ) + + + # Write table to CSV files + # split up by level1 grouping variables, so as not to exceed 5,000 row limit + iwalk( + split(km_estimates_table, km_estimates_table$group1), + ~ write_csv(.x, fs::path(output_dir, glue("km_estimates_{event_name}_table_{.y}.csv"))) + ) + + # Write table to a CSV file containing _only_ reporting milestones + km_estimates_milestones <- + km_estimates_table |> + filter(cohort_end_milestone == 1L) |> + mutate( + milestone_date = cohort_info$cohort_end_date, + milestone = "Cohort end" + ) |> + rename(days_since_cohort_start = time) |> + select(-cohort_end_milestone) + + # # extremely irritatingly we have to split this output up because it's too big to be shown (5000 row limit) + # iwalk( + # split(km_estimates_milestones, km_estimates_milestones$milestone), + # ~ write_csv(.x, fs::path(output_dir, glue("km_estimates_{event_name}_milestones_{.y}.csv"))) + # ) + write_csv( + km_estimates_milestones, + fs::path(output_dir, glue("km_estimates_{event_name}_milestones_cohort_end.csv")) + ) + return(km_estimates_table) +} + +## _______________________________________________________________________________________ +cat("## Report KM cumulative incidence of vaccination in a standardised table", "\n") +## _______________________________________________________________________________________ + + +km_estimates_vax_table <- get_all_km_estimates(data = data_prepared, event_name = "vax", event_time = "vax_time", event_indicator = "vax_indicator") +km_estimates_vax_alive_table <- get_all_km_estimates(filter(data_prepared, alive_and_registered), "vax_alive", event_time = "vax_time", event_indicator = "vax_indicator") + +# consider raw KM plots for disease burden too +# km_estimates_covid_admitted_table <- get_all_km_estimates(data = data_prepared, event_name = "covid_admitted", event_time = "covid_admitted_time", event_indicator = "covid_admitted_indicator") + + +# km_estimates_vax_table |> +# group_by(group1) |> +# summarise(n = n()) + + +# ________________________________ +# 2.A.i Kaplan-Meier plots (for server inspection only) +# ________________________________ +## function to make KM plots for the data ---- + +km_plot <- function(km_data, event_name, group1, group2) { + + group_names <- c(group1, group2) + + # plot km curves locally for checking (but probs not for release as these can be reconstructed from released data) + coverage_plot <- + bind_rows( + # this bit adds an extra row of data so that the firsrt line of the KM plot shows + km_data |> summarise( + time = 0, + cmlinc = 0, cmlinc.low = 0, cmlinc.high = 0, lagtime = 0, + .by = c("group1", "group1_value", "group2", "group2_value") + ), + km_data |> mutate( + lagtime = lag(time, 1L, 0L), # assumes the time-origin is zero + + .by = c("group1", "group1_value", "group2", "group2_value") + ) + ) |> + ggplot() + + geom_step(aes(x = time, y = cmlinc, group = group2_value, colour = group2_value), direction = "vh") + + geom_rect(aes(xmin = lagtime, xmax = time, ymin = cmlinc.low, ymax = cmlinc.high, group = group2_value, fill = group2_value), alpha = 0.1, colour = "transparent") + + facet_grid(rows = "group1_value") + + scale_color_brewer(type = "qual", palette = "Set1", na.value = "grey") + + scale_fill_brewer(type = "qual", palette = "Set1", guide = "none", na.value = "grey") + + scale_y_continuous(expand = expansion(mult = c(0, 0.01))) + + coord_cartesian(xlim = c(0, NA)) + + labs( + x = "Days since cohort start", + y = "Cumulative vaccination coverage", + colour = NULL, + title = NULL + ) + + theme_minimal() + + theme( + axis.line.x = element_line(colour = "black"), + panel.grid.minor.x = element_blank(), + legend.position = "inside", + legend.position.inside = c(.05, .95), + legend.justification = c(0, 1), + ) + + ggsave(fs::path(output_dir, glue("km_{event_name}_{paste0(group_names, collapse='_')}.png")), plot = coverage_plot) + # print(coverage_plot) +} + +# km_plot(km_estimates_all |> filter(group1 == "all", group2 == "crd"), "vax", "all", "crd") + +## _______________________________________________________________________________________ +cat("## Plot KM cumulative incidences of for inspection the server", "\n") +## _______________________________________________________________________________________ + + +# plot km estimates for inspection the server +walk2( + level_combos$group1, level_combos$group2, + .f = \(x, y) { + km_plot( + km_data = km_estimates_vax_table |> filter(group1 == x, group2 == y), + event_name = "vax", + group1 = x, + group2 = y + ) + } +) + + +################################################################################################################# +# 2.B Disease burden +# Incidence of disease-specific outcomes +################################################################################################################# + + +# list of reference categories, which overwrites default reference category if needed +contrasts_reference_levels <- list( + ageband4 = "65-74", + ageband13 = "65-69" +) + +## Function to output HRs and IRRs for disease burden comparing different subgroups ---- +# note that parglm is faster, but produces an annoying warning that "'mustart' will not be used" +# don't know how to get rid of it! + +adjusted_estimates <- function(data, subgroup, event_time, event_indicator, model = 1) { + # Model 1 is unadjusted + if (model == 1) { + poisson_formula <- as.formula(glue("event_indicator ~ {subgroup}")) + } + + # Model 2 adjusts for age using a restricted cubic spline with three knots and sex + if (model == 2) { + poisson_formula <- as.formula(glue("event_indicator ~ {subgroup} + sex + ns(age, 3)")) + + # use age-splines unless age is the subgroup of interest + if (subgroup %in% c("ageband4", "ageband13", "age_group")) { + poisson_formula <- as.formula(glue("event_indicator ~ {subgroup} + sex")) + } + + # do not adjust for sex if sex is the subgroup of interest + if (subgroup == "sex") { + poisson_formula <- as.formula(glue("event_indicator ~ {subgroup} + ns(age, 3)")) + } + } + # prepare dataset + data_outcome <- + data |> + mutate( + event_time = .data[[event_time]], + event_indicator = .data[[event_indicator]] + ) |> + select( + all_of(subgroup), + sex, age, + event_time, + event_indicator + ) |> + # to save memory, reduce dataset size by counting all distinct rows + # then use these counts as weights when modelling / summarising later + summarise( + count = n(), + .by = all_of(c(subgroup, "sex", "age", "event_time", "event_indicator")) + ) + + # how many possible values of the group are there + n_values <- n_distinct(data_outcome[[subgroup]]) + + # all levels of variable that exist in the data, sorted by factor (or alphanumeric if not) + # the first of these is the default reference level + all_levels <- unique(as.character(sort(data_outcome[[subgroup]]))) + + # explicitly include reference_level, and define contrasts if unusual + # for example, if we wanted to use not the default reference level for a factor + # this object is used inside a glm() call, to be passed to the `contrasts` argument + if ( + (subgroup %in% names(contrasts_reference_levels)) & + (n_values > 0) & + ifelse(is.null(contrasts_reference_levels[[subgroup]]), FALSE, contrasts_reference_levels[[subgroup]] %in% all_levels) + ) { + # create the object passed to `contrasts` argument in glm call to use a different reference category + subgroup_contrasts <- list(contr.treatment(all_levels, which(all_levels == contrasts_reference_levels[[subgroup]]))) + names(subgroup_contrasts) <- subgroup + reference_level <- contrasts_reference_levels[[subgroup]] + } else { + subgroup_contrasts <- NULL + reference_level <- all_levels[1] + } + + # cat(subgroup, "-", reference_level, " n=", n_values, " \n") + + # ------------------------------- + # Descriptive summary + # summarise total people, events, and person-time + # and add contrast label for merging with model output later + #-------------------------------- + data_summary <- + data_outcome |> + mutate(label = .data[[subgroup]]) |> + arrange(label) |> + summarise( + variable = subgroup, + n_obs = roundmid_any(sum(count), sdc_threshold), + n_event = roundmid_any(sum(event_indicator * count), sdc_threshold), + exposure = roundmid_any(sum(event_time * count), sdc_threshold), + + .by = label + ) |> + mutate( + label = as.character(label), + reference_row = label == reference_level, + contrast = glue("{subgroup}{label}") + ) + + # ------------------------------- + # Poisson regression + # Incidence rate ratios (IRR) + # ------------------------------- + + if (n_values > 1) { + + parglm_control <- parglm.control(maxit = 40, nthreads = 4) + + # fit the model + # if there is an error, just return an empty dataset rather than fail + data_poisson0 <- + tryCatch( + expr = { + data_outcome |> + parglm( + data = _, + formula = poisson_formula, + family = poisson, + offset = log(event_time), + control = parglm_control, + weights = count, + contrasts = subgroup_contrasts + ) |> + broom.helpers::tidy_and_attach(tidy_fun = broom.helpers::tidy_parameters, ci_method = "wald") |> + # broom.helpers::tidy_add_reference_rows() |> + broom.helpers::tidy_add_term_labels() |> + filter(variable == subgroup) |> + select(variable, label, estimate, std.error, conf.low, conf.high) + # note: the filter above is the same as doing marginaleffects::avg_comparisons(model, type = "link", variables = subgroup, comparison = "difference"), + # as long as there are no interaction terms between subgroup and anything else + # we use broom.helpers functions because it gives us the really nice variable and label info formatting for the outputted tidy dataset + # if we want to use avg_comparisons in future, then attach the nicely formatted meta info onto a broom::tidy(avg_comparisons) object + # or see "get_estimates_using_marginaleffects.R" script for a clue + }, + error = function(e) { + cat("error for subgroup", subgroup, ":", conditionMessage(e), "\n") + data_summary |> + select(variable, label) |> + mutate(estimate = NA_real_, std.error = NA_real_, conf.low = NA_real_, conf.high = NA_real_) + } + ) + + # combine summary and model outputs + data_poisson <- + full_join( + data_summary, + data_poisson0, + by = c("variable", "label"), + ) |> + transmute( + variable, label, reference_row, + n_obs, n_event, exposure, + ir = n_event / exposure, + irr_unadjusted = ir / ir[reference_row], + irr = exp(estimate), + irr.low = exp(conf.low), + irr.high = exp(conf.high), + irr.ln.std.error = std.error, + ) + + } else { + data_poisson <- data_summary |> select(-contrast) + } + + return(data_poisson) + + rm(data_outcome) + gc() +} + +# adjusted_estimates(data_prepared, "ageband4", "covid_admitted_time", "covid_admitted_indicator") + +################################################################################################################# +# 3. Comparative analyses +################################################################################################################# + +# 3.A Poisson regression +# Incidence rate ratios (IRRs) + +## function to loop over all group combinations and report IRR of level0 versus levelX ---- +# for a given outcome, loop over all groups combinations, obtaining contrasts for each using adjusted_estimates function, and combining into one file +# specifically, compare level2 groups amongst each other, for all people meeting level1 group criteria +get_all_estimates <- function(data, event_name, event_time, event_indicator) { + + l1ticker <<- "" + + estimates_list <- + level_combos |> + mutate( + estimates = map2( + group1, group2, + function(group1, group2) { + # just to keep track of how far along we are + if (l1ticker != group1) { + print(group1) + l1ticker <<- group1 + } + + summary_data <- + data |> + mutate( + label1 = data[[group1]], + ) |> + nest(.by = c(label1), .key = "group1_subset") |> + mutate( + estimates = map(group1_subset, \(group1_subset) { + bind_rows( + adjusted_estimates(group1_subset, group2, event_time, event_indicator, model = 1), + adjusted_estimates(group1_subset, group2, event_time, event_indicator, model = 2) + ) + }) + ) |> + select(-group1_subset) |> + unnest(estimates) |> + select(-variable) |> + rename(label2 = label) |> + mutate( + across(c(label1, label2), as.character) # to ensure the unnest() works later + ) + + return(summary_data) + rm(summary_data) + gc() + } + ) + ) |> + unnest(estimates) |> + select(group1, label1, group2, label2, everything()) # reorder columns + + write_csv(estimates_list, fs::path(output_dir, glue("contrasts_{event_name}.csv"))) + + # return(estimates_list) + +} + +## _______________________________________________________________________________________ +cat("## Get all IRR contrasts for vaccination and burden", "\n") +## _______________________________________________________________________________________ + +# IRR for vaccination, in the usual way +get_all_estimates(data_prepared, "vax", "vax_time", "vax_indicator") + +# IRR for vaccination, only looking at those who survived or stayed registered to the end of the season (collider bias! but matches UKHSA reports) +get_all_estimates( + filter(data_prepared, alive_and_registered), + "vax_alive", "vax_time", "vax_indicator" +) + +# IRR for burden, in the usual way +#TODO: add here mild +get_all_estimates(data_prepared, "admitted", "admitted_time", "admitted_indicator") +get_all_estimates(data_prepared, "admitted_primary", "admitted_primary_time", "admitted_primary_indicator") +get_all_estimates(data_prepared, "disease_death", "disease_death_time", "disease_death_indicator") + + +# ## Function to output length of stay quantiles for different subgroups ---- +# los_estimates <- function(data, subgroup, event_los) { +# +# subgroup_name <- deparse(substitute(subgroup)) +# +# # prepare dataset +# data_outcome <- +# data |> +# mutate( +# event_los = .data[[event_los]] +# ) |> +# select( +# all_of(subgroup), +# event_los +# ) +# +# # LoS summary stats +# +# data_los <- +# data_outcome |> +# mutate( +# variable = subgroup, +# label = .data[[subgroup]], +# ) |> +# summarise( +# n = roundmid_any(n(), sdc_threshold), +# n_at_least_1_event = roundmid_any(sum(!is.na(event_los)), sdc_threshold), +# median_los = quantile(event_los, 0.5, na.rm = TRUE), +# p10 = quantile(event_los, 0.1, na.rm = TRUE), +# p25 = quantile(event_los, 0.25, na.rm = TRUE), +# p75 = quantile(event_los, 0.75, na.rm = TRUE), +# p90 = quantile(event_los, 0.9, na.rm = TRUE), +# +# .by = c(variable, label) +# ) +# +# return(data_los) +# } +# +# +# los_estimates(data_prepared, "sex", "covid_admitted_los") +# +# ## function to get LoS across all group combinations ---- +# # for a given los outcome, loop over all groups combinations, obtaining los summaries for each using los_estimates function, and combining into one file +# get_all_los_estimates <- function(data, event_name, event_los) { +# +# estimates_list <- +# level_combos |> +# mutate( +# estimates = map2( +# group1, group2, +# \(group1, group2) { +# +# data |> +# mutate( +# label1 = data[[group1]], +# ) |> +# nest(.by = c(label1), .key = "group1_subset") |> +# mutate( +# estimates = map(group1_subset, \(group1_subset) { +# los_estimates(group1_subset, group2, event_los) +# }) +# ) |> +# select(-group1_subset) |> +# unnest(estimates) |> +# select(-variable) |> +# rename(label2 = label) |> +# mutate( +# across(c(label1, label2), as.character) +# ) +# } +# ) +# ) |> +# unnest(estimates) |> +# select(group1, label1, group2, label2, everything()) # reorder columns +# +# write_csv(estimates_list, fs::path(output_dir, glue("los_{event_name}.csv"))) +# +# } +# +# ## _______________________________________________________________________________________ +# cat("## Get all LOS values for burden", "\n") +# ## _______________________________________________________________________________________ +# +# get_all_los_estimates(data_prepared, "covid_admitted", "covid_admitted_los") diff --git a/analysis/4-snapshot/report_snapshot.R b/analysis/4-snapshot/report_snapshot.R deleted file mode 100644 index d397218..0000000 --- a/analysis/4-snapshot/report_snapshot.R +++ /dev/null @@ -1,1069 +0,0 @@ -# _______________________________________________________________________________________ -# Purpose: -# Report the distribution of vaccines schedules in different population subgroups as at a given date -# _______________________________________________________________________________________ - - -# Preliminaries ---- - -# Import libraries -library("tidyverse") -library("dtplyr") -library("lubridate") -library("glue") -library("here") -library("arrow") -library("survival") -library("splines") -library("parglm") - -# Import custom functions -source(here("analysis", "0-lib", "design.R")) - -args <- commandArgs(trailingOnly = TRUE) - -if (length(args) == 0) { - # use for interactive testing - # removeobjects <- FALSE - snapshot_date <- as.Date("20220321", format = "%Y%m%d") -} else { - # removeobjects <- TRUE - snapshot_date <- as.Date(args[[1]], format = "%Y%m%d") -} - -# how wide are the temporal bins for frequencies over time? in days -# use fix width, rather than months or bi-months or quarters, so that temporal denominator is even and the ratio of weekends to weekdays is fixed -temporal_resolution_history <- 28L - -# how wide are the temporal bins for frequencies over time for Kaplan-Meier plots? in days -temporal_resolution_km <- 7L - -campaign_info <- campaign_info |> filter(campaign_start_date == snapshot_date) - -# list2env(campaign_info, globalenv()) - -# overwrite primary milestone to match rounded kaplan meier curve -campaign_info$primary_milestone_days <- ceiling_any(campaign_info$primary_milestone_days, temporal_resolution_km) - - -# dates to round down to -# use this with `findInterval` until lubridate package is updated in the opensafely R image -# (then use `floor_date(date, unit=floor_dates`) -floor_dates <- seq( - as.Date("2020-06-01"), # monday - as.Date("2029-12-31"), # to monday! - by = temporal_resolution_history -) - -# create string representation of date in compact format YYYMMDD -snapshot_date_compact <- format(snapshot_date, "%Y%m%d") - -## Create output directory -output_dir <- here("output", "4-snapshot", glue("report_snapshot_{snapshot_date_compact}")) -fs::dir_create(output_dir) -options(width = 200) # set output width for capture.output - - -# Import processed data ---- -data_snapshot <- read_feather(here("output", "1-extract", glue("extract_snapshot_{snapshot_date_compact}.arrow"))) -data_fixed <- read_feather(here("output", "2-prepare", "prepare", "data_fixed.arrow")) - -capture.output( - skimr::skim_without_charts(data_snapshot), - file = fs::path(output_dir, "data_snapshot_skim.txt"), - split = FALSE -) - -# merge fixed data and vaccine data onto snapshot data -# note that in dummy data this doesn't work very well because patient IDs might not be matched across all datasets -data_combined <- - data_snapshot |> - lazy_dt() |> - left_join( - lazy_dt(data_fixed) |> select(patient_id, sex, ethnicity5, ethnicity16, death_date, covid_death_date), - by = "patient_id" - ) |> - # remove currently unused variables - select( - -covid_vax_prior_2_date, - -covid_vax_prior_2_product, - -covid_vax_prior_3_date, - -covid_vax_prior_3_product, - -stp, - -age # remove snapshot age and replace (below) with age_eligibility - ) |> - rename( - - # overwrite age to be that of a date later on in the campaign (rather than as at the start date) to allow "operational flexibility" - age = eligibility_age, - - # previous vaccine summary - # add more variables here based on covid_vax_prior_1_date, covid_vax_prior_2_date,... etc if needed - vax_count = covid_vax_prior_count, - last_vax_date = covid_vax_prior_1_date, - last_vax_product = covid_vax_prior_1_product, - - # info on first vaccine(s) received after snapshot date - next_vax_date = covid_vax_1_date, - next_vax_product = covid_vax_1_product, - # next2_vax_date = covid_vax_2_date, - ) |> - mutate( - all = "All", - - !!!standardise_demographic_characteristics, - !!!standardise_primis_and_extended_characteristics, - - age_above_eligiblity_threshold = (age >= campaign_info$age_threshold), - - # used to chose if the at risk group is all clinical risk variables or just immunosuppressed people - clinical_priority = .data[[campaign_info$clinical_priority]], - - clinical_priority_only = clinical_priority & !age_above_eligiblity_threshold, - - any_eligibility = age_above_eligiblity_threshold | clinical_priority | carehome_status, - - last_vax_product = fct_na_value_to_level(last_vax_product, "Unvaccinated"), - last_vax_date = if_else(vax_count == 0, study_dates$firstpossiblevax_date + as.integer(runif(n(), 0, 10)), last_vax_date), - # last_vax_week = floor_date(last_vax_date, unit = "week", week_start = 1), # starting on a monday - last_vax_period = floor_date(last_vax_date, unit = floor_dates), # use floor_dates[findInterval(last_vax_date, floor_dates)] if lubridate isn't working - - censor_date = pmin( - deregistered_date, - campaign_info$final_milestone_date, - study_dates$end_date, - na.rm = TRUE - ), - - # time from snapshot date until next vaccination - vax_time = as.integer(pmin(next_vax_date, death_date, censor_date, na.rm = TRUE) - snapshot_date) + 1L, # +1 because vaccination on snapshot date is allowed, but events at time zero are not - vax_indicator = (next_vax_date <= pmin(censor_date, death_date, na.rm = TRUE)) & !is.na(next_vax_date), - - # time from snapshot date until covid hospital admission - covid_admitted_time = as.integer(pmin(covid_admitted_date, death_date, censor_date, na.rm = TRUE) - snapshot_date) + 1L, - covid_admitted_indicator = (covid_admitted_date <= pmin(censor_date, death_date, na.rm = TRUE)) & !is.na(covid_admitted_date), - - # time from snapshot date until covid (primary position only) hospital admission - covid_admitted_primary_time = as.integer(pmin(covid_admitted_primary_date, death_date, censor_date, na.rm = TRUE) - snapshot_date) + 1L, - covid_admitted_primary_indicator = (covid_admitted_primary_date <= pmin(censor_date, death_date, na.rm = TRUE)) & !is.na(covid_admitted_primary_date), - - # time from snapshot date until covid critical care admission - covid_critcare_time = as.integer(pmin(covid_critcare_date, death_date, censor_date, na.rm = TRUE) - snapshot_date) + 1L, - covid_critcare_indicator = (covid_critcare_date <= pmin(censor_date, death_date, na.rm = TRUE)) & !is.na(covid_critcare_date), - - # time from snapshot date until covid death - covid_death_time = as.integer(pmin(covid_death_date, death_date, censor_date, na.rm = TRUE) - snapshot_date) + 1L, - covid_death_indicator = (covid_death_date <= pmin(censor_date, death_date, na.rm = TRUE)) & !is.na(covid_death_date), - - # time from snapshot date until death - death_time = as.integer(pmin(death_date, censor_date, na.rm = TRUE) - snapshot_date) + 1L, - death_indicator = (death_date <= pmin(censor_date, death_date, na.rm = TRUE)) & !is.na(death_date), - - # time from snapshot date until deregistration - deregistration_time = as.integer(pmin(deregistered_date, censor_date, na.rm = TRUE) - snapshot_date) + 1L, - deregistration_indicator = (deregistered_date <= pmin(censor_date, na.rm = TRUE)) & !is.na(deregistered_date), - - # indicator for if patient is alive and registered at the end of the campaign (for comparison with UKHSA reporting) - alive_and_registered = (!death_indicator) & (!deregistration_indicator) - ) |> - as_tibble() |> - mutate( - vax_status = case_when( - vax_indicator ~ "vaccinated", - (death_date <= censor_date) & !is.na(death_date) ~ "died", - (censor_date < death_date) | is.na(death_date) ~ "censored", - .default = NA_character_ - ) |> factor(levels = c("censored", "vaccinated", "death")), - vax_status_product = if_else(vax_indicator, next_vax_product, vax_status), - next_vax_product_uncensored = if_else(vax_indicator, next_vax_product, NA_character_), - - # setting missing values to explicit missing - across(where(is.factor) | where(is.character), ~ fct_drop(fct_na_value_to_level(.x, level = "(Missing)"))) - ) - -# remove snapshot and fixed datasets to save memory, as these is all captured in the data_combiend dataset. -rm(data_snapshot, data_fixed) -# recover memory -gc() - - -capture.output( - skimr::skim_without_charts(data_combined), - file = fs::path(output_dir, "data_combined_skim.txt"), - split = FALSE -) - -# ________________________________________________________________________________________ -# Pre-snapshot date vaccine history, stratified by characteristic recorded on the snapshot_date ---- -# ________________________________________________________________________________________ - - -## _______________________________________________________________________________________ -cat("## Report info on date of last vaccination", "\n") -## _______________________________________________________________________________________ - -# This produces data and a plot showing the distribution of most recent prior vaccination date (to the nearest 4 weeks) -# if no documented prior vaccination, then values are stacked on the LHS of the chart - -plot_date_of_last_dose <- function(subgroup) { - - over2years_dummy_date <- (snapshot_date - ceiling_any(365 * 2, 7)) - - summary_by <- - data_combined |> - mutate( - subgroup = .data[[subgroup]] - ) |> - lazy_dt() |> - group_by(subgroup, last_vax_product, last_vax_period) |> - summarise( - n = roundmid_any(n(), sdc_threshold) - ) |> - ungroup() |> - mutate( - # if last vaccination date was over 2 years ago, replace with dummy date - last_vax_period = if_else( - (last_vax_period < over2years_dummy_date) | is.na(last_vax_period), - over2years_dummy_date - 42, - last_vax_period - ) - ) |> - as_tibble() |> - complete( - subgroup, last_vax_product, last_vax_period, - fill = list(n = 0) - ) - - breaks <- seq(as.Date("2021-01-01"), as.Date("2028-06-01"), by = "6 month") - - temp_plot <- - summary_by |> - ggplot() + - geom_col( - aes(x = last_vax_period, y = n, fill = last_vax_product, group = last_vax_product), - alpha = 0.5, - position = position_stack(reverse = TRUE), - # position=position_identity(), - width = temporal_resolution_history - ) + - facet_grid( - rows = vars(subgroup), - switch = "y", - space = "free_x", - scales = "free_x" - ) + - labs( - x = "Date of most recent COVID-19 vaccine", - y = NULL, - fill = NULL - ) + - scale_fill_brewer( - palette = "Set2", - na.value = "grey50", - labels = function(breaks) { - breaks[is.na(breaks)] <- "Other" - breaks - } - ) + - scale_x_date( - breaks = c(over2years_dummy_date - 42, breaks), - date_minor_breaks = "month", - # labels = ~{c("Unvaccinated", scales::label_date("%Y")(.x[-1]))}, - labels = c("+2 years prior)", scales::label_date("%Y-%b")(breaks)), - ) + - theme_minimal() + - theme( - axis.text.x.top = element_text(hjust = 0), - axis.text.x.bottom = element_text(hjust = 0), - strip.text.y.left = element_text(angle = 0, hjust = 1), - axis.ticks.x = element_line(), - strip.placement = "outside", - axis.text.y = element_blank(), - legend.position = "bottom" - ) - - # print(temp_plot) - - ggsave(fs::path(output_dir, glue("last_vax_date_{subgroup}.png")), plot = temp_plot) - - # write tables that capture underlying plotting data - # write_csv(summary_by, fs::path(output_dir, glue("last_vax_date_{subgroup}.csv"))) -} - -for (group in level1_group) { - plot_date_of_last_dose(group) -} - -## _______________________________________________________________________________________ -cat("## Report info on prior dose count and product type", "\n") -## _______________________________________________________________________________________ - -plot_vax_count <- function(subgroup) { - - summary_by <- - data_combined |> - mutate( - subgroup = .data[[subgroup]], - vax_count_group = cut(vax_count, c(-Inf, 0, 2, 4, Inf), labels = c("0", "1-2", "3-4", "5+")), - ) |> - lazy_dt() |> - group_by(subgroup, vax_count_group) |> - summarise( - n = roundmid_any(n(), sdc_threshold), - ) |> - ungroup() |> - as_tibble() |> - complete( - subgroup, vax_count_group, - fill = list(n = 0) - ) |> - group_by(subgroup) |> - mutate( - row_total = sum(n), - prop = ifelse(row_total > 0, n / row_total * 100, 0), - ) |> - ungroup() - - temp_plot <- - ggplot(summary_by) + - geom_bar( - aes(x = prop, y = subgroup, width = row_total, fill = vax_count_group), - stat = "identity", # position = "fill", - position = position_stack(reverse = TRUE), - ) + - facet_grid( - rows = vars(subgroup), - scales = "free_y", - space = "free_y") + - labs( - x = "%", - y = NULL, - fill = "Vaccine count" - ) + - scale_fill_brewer( - palette = "Set2", - na.value = "grey50", - labels = function(breaks) { - breaks[is.na(breaks)] <- "(Missing)" - breaks - } - ) + - theme_minimal() + - theme( - axis.text.x.top = element_text(hjust = 0), - axis.text.x.bottom = element_text(hjust = 0), - strip.text.y.left = element_text(angle = 0, hjust = 1), - axis.ticks.x = element_line(), - strip.text = element_blank(), - legend.position = "bottom" - ) + - NULL - - # print(temp_plot) - - ggsave(fs::path(output_dir, glue("vax_count_{subgroup}.png")), plot = temp_plot) - - # write tables that capture underlying plotting data - write_csv(summary_by, fs::path(output_dir, glue("vax_count_{subgroup}.csv"))) -} - - -for (group in level1_group) { - plot_vax_count(group) -} - - -## _______________________________________________________________________________________ -cat("## Report info in a standardised table", "\n") -## _______________________________________________________________________________________ - -# function to print table for an abritrary number of grouping variables - -# groups are passed as strings via dots (...) -table_prior_vax_summary <- function(...) { - - group_names <- c(...) - - summary_table <- - data_combined |> - mutate(days_since_vax = snapshot_date - last_vax_date) |> - group_by(across(all_of(group_names))) |> - lazy_dt() |> - summarise( - # Dose counts - total = roundmid_any(n(), sdc_threshold), - count_n0 = roundmid_any(sum(vax_count == 0, na.rm = TRUE), sdc_threshold), - count_n1 = roundmid_any(sum(vax_count == 1, na.rm = TRUE), sdc_threshold), - count_n2 = roundmid_any(sum(vax_count == 2, na.rm = TRUE), sdc_threshold), - count_n3 = roundmid_any(sum(vax_count == 3, na.rm = TRUE), sdc_threshold), - count_n4 = roundmid_any(sum(vax_count == 4, na.rm = TRUE), sdc_threshold), - count_n5plus = roundmid_any(sum(vax_count >= 5, na.rm = TRUE), sdc_threshold), - # Dose summary - count_median = quantile(vax_count, probs = 0.5, na.rm = TRUE), - count_p10 = quantile(vax_count, probs = 0.10, na.rm = TRUE), - count_p25 = quantile(vax_count, probs = 0.25, na.rm = TRUE), - count_p75 = quantile(vax_count, probs = 0.75, na.rm = TRUE), - count_p90 = quantile(vax_count, probs = 0.90, na.rm = TRUE), - # Vaccination in past 12 and 24 months - days_since_n12m = roundmid_any(sum(days_since_vax <= 365, na.rm = TRUE), sdc_threshold), - days_since_n24m = roundmid_any(sum(days_since_vax <= 365 * 2, na.rm = TRUE), sdc_threshold), - # Time since last dose - days_since_median = quantile(days_since_vax, probs = 0.5, na.rm = TRUE), - days_since_p10 = quantile(days_since_vax, probs = 0.10, na.rm = TRUE), - days_since_p25 = quantile(days_since_vax, probs = 0.25, na.rm = TRUE), - days_since_p75 = quantile(days_since_vax, probs = 0.75, na.rm = TRUE), - days_since_p90 = quantile(days_since_vax, probs = 0.90, na.rm = TRUE), - - .groups = "drop" - ) |> - mutate( - # Dose percentages - put this here and not in earlier summarise step so that it works with dtplyr - count_pct0 = round(count_n0 * 100 / total, 1), - count_pct1 = round(count_n1 * 100 / total, 1), - count_pct2 = round(count_n2 * 100 / total, 1), - count_pct3 = round(count_n3 * 100 / total, 1), - count_pct4 = round(count_n4 * 100 / total, 1), - count_pct5plus = round(count_n5plus * 100 / total, 1), - # Vaccination % in past 12 and 24 months - days_since_pct12m = round(days_since_n12m * 100 / total, 1), - days_since_pct24m = round(days_since_n24m * 100 / total, 1), - ) |> - as_tibble() - - # subgroup_name <- map_chr(rlang::quos(...), rlang::as_name) |> paste0(collapse = "_") - - # Write table to a CSV file - # write_csv(summary_table, fs::path(output_dir, glue("prior_vax_table_{subgroup_name}.csv"))) - - return(summary_table) -} - - -# for testing function interactively -# table_prior_vax_summary("all", "all") -# table_prior_vax_summary("all", "ageband4") -# table_prior_vax_summary("ageband4", "all") -# table_prior_vax_summary("ageband4", "sex") - -# loop over all group1 and group2 variable combinations and combine into one big dataset -prior_vax_summary_table_all <- - level_combos |> - mutate( - prior_vax_summary = map2( - group1, group2, - .f = function(x, y) { - - if (is.na(y)) y <- NULL - lookup <- c(group1_value = x, group2_value = y) - - table_prior_vax_summary(x, y) |> - mutate(across(c(all_of(c(x, y))), as.character)) |> - rename(any_of(lookup)) - } - ) - ) |> - unnest(prior_vax_summary) |> - select(group1, group1_value, group2, group2_value, everything()) - -# Write table to a CSV file -write_csv(prior_vax_summary_table_all, fs::path(output_dir, glue("prior_vax_table.csv"))) - - -# ________________________________________________________________________________________ -# Post-snapshot vaccine coverage, stratified by characteristic recorded on the snapshot_date ---- -# ________________________________________________________________________________________ - -# This code borrows heavily from the KM reusable action https://github.com/opensafely-actions/kaplan-meier-function/blob/main/analysis/km.R -# The reason not to use the KM reusable action directly is that we need to reuse it multiple times across many different stratification variables -# This would create a large number of project.yaml actions, so it's easier to do the repeats within a single script -# The reusable action could be modified to enable multiple stratified analyses to be run, but it's a bit faffy and a bit scope-creepy - - -## tests ---- - -times_count <- table(cut(data_combined$vax_time, c(-Inf, 0, 1, Inf), right = FALSE, labels = c("<0", "0", ">0")), useNA = "ifany") - -if (!identical(as.integer(times_count), c(0L, 0L, nrow(data_combined)))) { - print(times_count) - stop("all event times must be strictly positive") -} - -## Function to calculate KM estimates for a given stratification variable ---- - -# group variables are provided as characters via dots (...) -# resolution argument is the precision used for the time dimension. If zero, then original resolution is used. -km_estimates <- function(data, group_name1, group_name2, event_name, event_time, event_indicator, resolution = 0) { - - group_names <- c(group_name1, group_name2) - - # if (is.na(group_name2)) { - # group_name2 <- "all" - # } - - data_outcome <- - data |> - select( - patient_id, - all_of(group_names), - event_time = any_of(event_time), - event_indicator = any_of(event_indicator) - ) |> - mutate(event_time = ceiling_any(event_time, resolution)) - - data_km <- - data_outcome |> - group_by(across(all_of(group_names))) |> - nest() |> - mutate( - surv_obj_tidy = map( - data, ~ { - survfit( - Surv(event_time, event_indicator) ~ 1, - data = .x, - conf.type = "log-log" - ) |> - broom::tidy() |> - add_row( - time = 0, # assumes time origin is zero - n.risk = 0, - n.event = 0, - n.censor = 0, - estimate = 1, - conf.low = 1, - conf.high = 1, - .before = 1L - ) |> - complete( - time = seq(0L, campaign_info$final_milestone_days, resolution), # fill in 1 row for each period (defined by resolution) of follow up - fill = list(n.event = 0L, n.censor = 0L) # fill in zero events on those days - ) |> - fill( - n.risk, - .direction = "up" - ) |> - fill( - estimate, conf.low, conf.high, - .direction = "down" - ) - } - ), - ) |> - select(-data) |> - unnest(surv_obj_tidy) |> - mutate( - - # disclosure control - n.risk = ceiling_any(n.risk, sdc_threshold), - estimate = plyr::round_any(estimate, sdc_threshold / nth(n.risk, 2)), # use 2nd value as this skips the t=0 row where n.risk=0 - conf.low = plyr::round_any(conf.low, sdc_threshold / nth(n.risk, 2)), - conf.high = plyr::round_any(conf.high, sdc_threshold / nth(n.risk, 2)), - - # cumulative incidence - cmlinc = 1 - estimate, - cmlinc.low = 1 - conf.high, - cmlinc.high = 1 - conf.low, - ) - - rm(data_outcome) - - # print(data_km) - - # write tables that capture underlying plotting data - data_km_nozero <- - data_km |> - arrange(across(all_of(group_names))) |> - ungroup() |> - filter(time != 0) |> - select( - all_of(group_names), - time, - cmlinc, - cmlinc.low, - cmlinc.high, - ) - - # commented out because outputting the entire data creates a dataset that is too large to output-check on it's own - # write_csv(data_km_nozero, fs::path(output_dir, glue("km_{event_name}_{paste0(group_names, collapse='_')}.csv"))) - - return(data_km_nozero) - -} - -# for testing function interactively -# km_estimates_vax <- partial( -# km_estimates, data = data_combined, event_name = "vax", event_time = "vax_time", event_indicator = "vax_indicator", resolution = temporal_resolution_km -# ) -# km_estimates_vax("all", "all") -# km_estimates_vax("all", "ageband4") -# km_estimates_vax("ageband4", "all") -# km_estimates_vax("ageband4", "sex") - -## Function to get km estimates over all subgroup combinations ---- - -get_all_km_estimates <- function(data, event_name, event_time, event_indicator, resolution) { - # loop over all group1 and group2 variable combinations and combine into one big dataset - km_estimates_table <- - level_combos |> - mutate( - km_summary = map2( - group1, group2, - .f = \(x, y) { - - if (is.na(y)) y <- NULL - lookup <- c(group1_value = x, group2_value = y) - - km_estimates(data = data, group_name1 = x, group_name2 = y, event_name = event_name, event_time = event_time, event_indicator = event_indicator, resolution = temporal_resolution_km) |> - mutate(across(c(all_of(c(x, y))), as.character)) |> - rename(all_of(c(group1_value = x, group2_value = y))) - } - ) - ) |> - unnest(km_summary) |> - select(group1, group1_value, group2, group2_value, everything()) |> - mutate( - early_milestone = (time == campaign_info$early_milestone_days) * 1L, - primary_milestone = (time == campaign_info$primary_milestone_days) * 1L, - final_milestone = (time == campaign_info$final_milestone_days) * 1L, - ) - - - # Write table to CSV files - # split up by level1 grouping variables, so as not to exceed 5,000 row limit - iwalk( - split(km_estimates_table, km_estimates_table$group1), - ~ write_csv(.x, fs::path(output_dir, glue("km_estimates_{event_name}_table_{.y}.csv"))) - ) - - # Write table to a CSV file containing _only_ reporting milestones - km_estimates_milestones <- - km_estimates_table |> - filter( - time %in% (c(campaign_info$early_milestone_days, campaign_info$primary_milestone_days, campaign_info$final_milestone_days) * 1L) - ) |> - mutate( - milestone_date = case_when( - early_milestone == 1 ~ campaign_info$early_milestone_date, - primary_milestone == 1 ~ campaign_info$primary_milestone_date, - final_milestone == 1 ~ campaign_info$final_milestone_date, - ), - milestone = case_when( - early_milestone == 1 ~ "Early", - primary_milestone == 1 ~ "Primary", - final_milestone == 1 ~ "Final", - ) - ) |> - rename(days_since_campaign_start = time) |> - select(-early_milestone, -primary_milestone, -final_milestone) - - # extremely irritatingly we have to split this output up because it's too big to be shown (5000 row limit) - iwalk( - split(km_estimates_milestones, km_estimates_milestones$milestone), - ~ write_csv(.x, fs::path(output_dir, glue("km_estimates_{event_name}_milestones_{.y}.csv"))) - ) - - return(km_estimates_table) -} - -## _______________________________________________________________________________________ -cat("## Report KM cumulative incidence of vaccination in a standardised table", "\n") -## _______________________________________________________________________________________ - - -km_estimates_vax_table <- get_all_km_estimates(data = data_combined, event_name = "vax", event_time = "vax_time", event_indicator = "vax_indicator") -km_estimates_vax_alive_table <- get_all_km_estimates(filter(data_combined, alive_and_registered), "vax_alive", event_time = "vax_time", event_indicator = "vax_indicator") - -# consider raw KM plots for disease burden too -# km_estimates_covid_admitted_table <- get_all_km_estimates(data = data_combined, event_name = "covid_admitted", event_time = "covid_admitted_time", event_indicator = "covid_admitted_indicator") - - -# km_estimates_vax_table |> -# group_by(group1) |> -# summarise(n = n()) - - -## function to make KM plots for the data ---- - -km_plot <- function(km_data, event_name, group1, group2) { - - group_names <- c(group1, group2) - - # plot km curves locally for checking (but probs not for release as these can be reconstructed from released data) - coverage_plot <- - bind_rows( - # this bit adds an extra row of data so that the firsrt line of the KM plot shows - km_data |> summarise( - time = 0, - cmlinc = 0, cmlinc.low = 0, cmlinc.high = 0, lagtime = 0, - .by = c("group1", "group1_value", "group2", "group2_value") - ), - km_data |> mutate( - lagtime = lag(time, 1L, 0L), # assumes the time-origin is zero - - .by = c("group1", "group1_value", "group2", "group2_value") - ) - ) |> - ggplot() + - geom_step(aes(x = time, y = cmlinc, group = group2_value, colour = group2_value), direction = "vh") + - geom_rect(aes(xmin = lagtime, xmax = time, ymin = cmlinc.low, ymax = cmlinc.high, group = group2_value, fill = group2_value), alpha = 0.1, colour = "transparent") + - facet_grid(rows = "group1_value") + - scale_color_brewer(type = "qual", palette = "Set1", na.value = "grey") + - scale_fill_brewer(type = "qual", palette = "Set1", guide = "none", na.value = "grey") + - scale_y_continuous(expand = expansion(mult = c(0, 0.01))) + - coord_cartesian(xlim = c(0, NA)) + - labs( - x = "Days since start of campaign", - y = "Cumulative Incidence", - colour = NULL, - title = NULL - ) + - theme_minimal() + - theme( - axis.line.x = element_line(colour = "black"), - panel.grid.minor.x = element_blank(), - legend.position = "inside", - legend.position.inside = c(.05, .95), - legend.justification = c(0, 1), - ) - - ggsave(fs::path(output_dir, glue("km_{event_name}_{paste0(group_names, collapse='_')}.png")), plot = coverage_plot) - # print(coverage_plot) -} - -# km_plot(km_estimates_all |> filter(group1 == "all", group2 == "crd"), "vax", "all", "crd") - -## _______________________________________________________________________________________ -cat("## Plot KM cumulative incidences of for inspection the server", "\n") -## _______________________________________________________________________________________ - - -# plot km estimates for inspection the server -walk2( - level_combos$group1, level_combos$group2, - .f = \(x, y) { - km_plot( - km_data = km_estimates_vax_table |> filter(group1 == x, group2 == y), - event_name = "vax", - group1 = x, - group2 = y - ) - } -) - - -# ________________________________________________________________________________________ -# Post-snapshot Covid-19 vaccination and disease burden up to final milestone, stratified by characteristics recorded on the snapshot_date ---- -# ________________________________________________________________________________________ - - -# list of reference categories, which overwrites default reference category if needed -contrasts_reference_levels <- list( - ageband4 = "65-74", - ageband13 = "65-69" -) - -## Function to output HRs and IRRs for disease burden comparing different subgroups ---- -# note that parglm is faster, but produces an annoying warning that "'mustart' will not be used" -# don't know how to get rid of it! - -adjusted_estimates <- function(data, subgroup, event_time, event_indicator) { - # use age-splines unless age is the subgroup of interest - - poisson_formula <- as.formula(glue("event_indicator ~ {subgroup} + sex + ns(age, 3)")) - if (subgroup == "ageband4") poisson_formula <- as.formula(glue("event_indicator ~ ageband4 + sex")) - if (subgroup == "ageband13") poisson_formula <- as.formula(glue("event_indicator ~ ageband13 + sex")) - - # prepare dataset - data_outcome <- - data |> - mutate( - event_time = .data[[event_time]], - event_indicator = .data[[event_indicator]] - ) |> - select( - all_of(subgroup), - sex, age, - event_time, - event_indicator - ) |> - # to save memory, reduce dataset size by counting all distinct rows - # then use these counts as weights when modelling / summarising later - summarise( - count = n(), - .by = all_of(c(subgroup, "sex", "age", "event_time", "event_indicator")) - ) - - # how many possible values of the group are there - n_values <- n_distinct(data_outcome[[subgroup]]) - - # all levels of variable that exist in the data, sorted by factor (or alphanumeric if not) - # the first of these is the default reference level - all_levels <- unique(as.character(sort(data_outcome[[subgroup]]))) - - # explicitly include reference_level, and define contrasts if unusual - # for example, if we wanted to use not the default reference level for a factor - # this object is used inside a glm() call, to be passed to the `contrasts` argument - if ( - (subgroup %in% names(contrasts_reference_levels)) & - (n_values > 0) & - ifelse(is.null(contrasts_reference_levels[[subgroup]]), FALSE, contrasts_reference_levels[[subgroup]] %in% all_levels) - ) { - # create the object passed to `contrasts` argument in glm call to use a different reference category - subgroup_contrasts <- list(contr.treatment(all_levels, which(all_levels == contrasts_reference_levels[[subgroup]]))) - names(subgroup_contrasts) <- subgroup - reference_level <- contrasts_reference_levels[[subgroup]] - } else { - subgroup_contrasts <- NULL - reference_level <- all_levels[1] - } - - # cat(subgroup, "-", reference_level, " n=", n_values, " \n") - - # summarise total people, events, and person-time - # and add contrast label for merging with model output later - data_summary <- - data_outcome |> - mutate(label = .data[[subgroup]]) |> - arrange(label) |> - summarise( - variable = subgroup, - n_obs = roundmid_any(sum(count), sdc_threshold), - n_event = roundmid_any(sum(event_indicator * count), sdc_threshold), - exposure = roundmid_any(sum(event_time * count), sdc_threshold), - - .by = label - ) |> - mutate( - label = as.character(label), - reference_row = label == reference_level, - contrast = glue("{subgroup}{label}") - ) - - # IRR model - - if (n_values > 1) { - - parglm_control <- parglm.control(maxit = 40, nthreads = 4) - - # fit the model - # if there is an error, just return an empty dataset rather than fail - data_poisson0 <- - tryCatch( - expr = { - data_outcome |> - parglm( - data = _, - formula = poisson_formula, - family = poisson, - offset = log(event_time), - control = parglm_control, - weights = count, - contrasts = subgroup_contrasts - ) |> - broom.helpers::tidy_and_attach(tidy_fun = broom.helpers::tidy_parameters, ci_method = "wald") |> - # broom.helpers::tidy_add_reference_rows() |> - broom.helpers::tidy_add_term_labels() |> - filter(variable == subgroup) |> - select(variable, label, estimate, std.error, conf.low, conf.high) - # note: the filter above is the same as doing marginaleffects::avg_comparisons(model, type = "link", variables = subgroup, comparison = "difference"), - # as long as there are no interaction terms between subgroup and anything else - # we use broom.helpers functions because it gives us the really nice variable and label info formatting for the outputted tidy dataset - # if we want to use avg_comparisons in future, then attach the nicely formatted meta info onto a broom::tidy(avg_comparisons) object - # or see "get_estimates_using_marginaleffects.R" script for a clue - }, - error = function(e) { - cat("error for subgroup", subgroup, ":", conditionMessage(e), "\n") - data_summary |> - select(variable, label) |> - mutate(estimate = NA_real_, std.error = NA_real_, conf.low = NA_real_, conf.high = NA_real_) - } - ) - - # combine summary and model outputs - data_poisson <- - full_join( - data_summary, - data_poisson0, - by = c("variable", "label"), - ) |> - transmute( - variable, label, reference_row, - n_obs, n_event, exposure, - ir = n_event / exposure, - irr_unadjusted = ir / ir[reference_row], - irr = exp(estimate), - irr.low = exp(conf.low), - irr.high = exp(conf.high), - irr.ln.std.error = std.error, - ) - - } else { - data_poisson <- data_summary |> select(-contrast) - } - - return(data_poisson) - - rm(data_outcome) - gc() -} - -# adjusted_estimates(data_combined, "ageband4", "covid_admitted_time", "covid_admitted_indicator") - -## function to loop over all group combinations and report IRR of level0 versus levelX ---- -# for a given outcome, loop over all groups combinations, obtaining contrasts for each using adjusted_estimates function, and combining into one file -# specifically, compare level2 groups amongst each other, for all people meeting level1 group criteria -get_all_estimates <- function(data, event_name, event_time, event_indicator) { - - l1ticker <<- "" - - estimates_list <- - level_combos |> - mutate( - estimates = map2( - group1, group2, - function(group1, group2) { - # just to keep track of how far along we are - if (l1ticker != group1) { - print(group1) - l1ticker <<- group1 - } - - summary_data <- - data |> - mutate( - label1 = data[[group1]], - ) |> - nest(.by = c(label1), .key = "group1_subset") |> - mutate( - estimates = map(group1_subset, \(group1_subset) { - adjusted_estimates(group1_subset, group2, event_time, event_indicator) - }) - ) |> - select(-group1_subset) |> - unnest(estimates) |> - select(-variable) |> - rename(label2 = label) |> - mutate( - across(c(label1, label2), as.character) # to ensure the unnest() works later - ) - - return(summary_data) - rm(summary_data) - gc() - } - ) - ) |> - unnest(estimates) |> - select(group1, label1, group2, label2, everything()) # reorder columns - - write_csv(estimates_list, fs::path(output_dir, glue("contrasts_{event_name}.csv"))) - - # return(estimates_list) - -} - - -## _______________________________________________________________________________________ -cat("## Get all IRR contrasts for covid vaccination and burden", "\n") -## _______________________________________________________________________________________ - -# IRR for vaccination, in the usual way -get_all_estimates(data_combined, "vax", "vax_time", "vax_indicator") - -# IRR for vaccination, only looking at those who survived or stayed registered to the end of the season (collider bias! but matches UKHSA reports) -get_all_estimates( - filter(data_combined, alive_and_registered), - "vax_alive", "vax_time", "vax_indicator" -) - -# IRR for COVID burden, in the usual way -get_all_estimates(data_combined, "covid_admitted", "covid_admitted_time", "covid_admitted_indicator") -get_all_estimates(data_combined, "covid_admitted_primary", "covid_admitted_primary_time", "covid_admitted_primary_indicator") -get_all_estimates(data_combined, "covid_critcare", "covid_critcare_time", "covid_critcare_indicator") -get_all_estimates(data_combined, "covid_death", "covid_death_time", "covid_death_indicator") - - -## Function to output length of stay quantiles for different subgroups ---- -los_estimates <- function(data, subgroup, event_los) { - - subgroup_name <- deparse(substitute(subgroup)) - - # prepare dataset - data_outcome <- - data |> - mutate( - event_los = .data[[event_los]] - ) |> - select( - all_of(subgroup), - event_los - ) - - # LoS summary stats - - data_los <- - data_outcome |> - mutate( - variable = subgroup, - label = .data[[subgroup]], - ) |> - summarise( - n = roundmid_any(n(), sdc_threshold), - n_at_least_1_event = roundmid_any(sum(!is.na(event_los)), sdc_threshold), - median_los = quantile(event_los, 0.5, na.rm = TRUE), - p10 = quantile(event_los, 0.1, na.rm = TRUE), - p25 = quantile(event_los, 0.25, na.rm = TRUE), - p75 = quantile(event_los, 0.75, na.rm = TRUE), - p90 = quantile(event_los, 0.9, na.rm = TRUE), - - .by = c(variable, label) - ) - - return(data_los) -} - - -los_estimates(data_combined, "sex", "covid_admitted_los") - -## function to get LoS across all group combinations ---- -# for a given los outcome, loop over all groups combinations, obtaining los summaries for each using los_estimates function, and combining into one file -get_all_los_estimates <- function(data, event_name, event_los) { - - estimates_list <- - level_combos |> - mutate( - estimates = map2( - group1, group2, - \(group1, group2) { - - data |> - mutate( - label1 = data[[group1]], - ) |> - nest(.by = c(label1), .key = "group1_subset") |> - mutate( - estimates = map(group1_subset, \(group1_subset) { - los_estimates(group1_subset, group2, event_los) - }) - ) |> - select(-group1_subset) |> - unnest(estimates) |> - select(-variable) |> - rename(label2 = label) |> - mutate( - across(c(label1, label2), as.character) - ) - } - ) - ) |> - unnest(estimates) |> - select(group1, label1, group2, label2, everything()) # reorder columns - - write_csv(estimates_list, fs::path(output_dir, glue("los_{event_name}.csv"))) - -} - -## _______________________________________________________________________________________ -cat("## Get all LOS values for burden", "\n") -## _______________________________________________________________________________________ - -get_all_los_estimates(data_combined, "covid_admitted", "covid_admitted_los") diff --git a/codelists/codelists.json b/codelists/codelists.json index ea0170d..47b8a51 100644 --- a/codelists/codelists.json +++ b/codelists/codelists.json @@ -317,6 +317,18 @@ "url": "https://www.opencodelists.org/codelist/primis-covid19-vacc-uptake/reside_cod/v2.5/", "downloaded_at": "2026-01-07 16:52:31.359052Z", "sha": "906a866473bd4c99a48ddb1d86c08986d228ba38" + }, + "opensafely-influenza-identification-secondary-care.csv": { + "id": "opensafely/influenza-identification-secondary-care/7d33be15", + "url": "https://www.opencodelists.org/codelist/opensafely/influenza-identification-secondary-care/7d33be15/", + "downloaded_at": "2026-06-19 17:39:13.237202Z", + "sha": "29f102e6fd667e9ed466cbfc06c310d7d781991e" + }, + "opensafely-rsv-identification-secondary-care.csv": { + "id": "opensafely/rsv-identification-secondary-care/5ebecaf4", + "url": "https://www.opencodelists.org/codelist/opensafely/rsv-identification-secondary-care/5ebecaf4/", + "downloaded_at": "2026-06-19 17:39:13.327786Z", + "sha": "ab9f9712939394b9b8f6a261ab3afb2d9f94daa4" } } } \ No newline at end of file diff --git a/codelists/codelists.txt b/codelists/codelists.txt index d265cfd..5488fe6 100644 --- a/codelists/codelists.txt +++ b/codelists/codelists.txt @@ -116,3 +116,12 @@ nhsd-primary-care-domain-refsets/csflatrisk1_cod/20250912 # Homeless primis-covid19-vacc-uptake/homeless_cod/v2.5 primis-covid19-vacc-uptake/reside_cod/v2.5 + + +# Outcomes + +# flu +opensafely/influenza-identification-secondary-care/7d33be15 + +# rsv +opensafely/rsv-identification-secondary-care/5ebecaf4 diff --git a/codelists/opensafely-influenza-identification-secondary-care.csv b/codelists/opensafely-influenza-identification-secondary-care.csv new file mode 100644 index 0000000..a02b8b9 --- /dev/null +++ b/codelists/opensafely-influenza-identification-secondary-care.csv @@ -0,0 +1,10 @@ +code,term +J09,Influenza due to identified zoonotic or pandemic influenza virus +J10,Influenza due to identified seasonal influenza virus +J100,"Influenza with pneumonia, seasonal influenza virus identified" +J101,"Influenza with other respiratory manifestations, seasonal influenza virus identified" +J108,"Influenza with other manifestations, seasonal influenza virus identified" +J11,"Influenza, virus not identified" +J110,"Influenza with pneumonia, virus not identified" +J111,"Influenza with other respiratory manifestations, virus not identified" +J118,"Influenza with other manifestations, virus not identified" diff --git a/codelists/opensafely-rsv-identification-secondary-care.csv b/codelists/opensafely-rsv-identification-secondary-care.csv new file mode 100644 index 0000000..accbec3 --- /dev/null +++ b/codelists/opensafely-rsv-identification-secondary-care.csv @@ -0,0 +1,7 @@ +code,term +B974,Respiratory syncytial virus as the cause of diseases classified to other chapters +J121,Respiratory syncytial virus pneumonia +J205,Acute bronchitis due to respiratory syncytial virus +J21,Acute bronchiolitis +J210,Acute bronchiolitis due to respiratory syncytial virus +J219,"Acute bronchiolitis, unspecified" diff --git a/project.yaml b/project.yaml index ef5f76e..fc99e29 100644 --- a/project.yaml +++ b/project.yaml @@ -1,249 +1,302 @@ -version: '4.0' +version: '5.0' actions: - extract_fixed: - run: ehrql:v1 generate-dataset analysis/1-extract/dataset_definition_fixed.py - --output output/1-extract/extract_fixed.arrow - --dummy-data-file analysis/1-extract/dummy-data/dummy_fixed.arrow +# covid: wave_1 + extract_covid_wave_1: + run: ehrql:v1 generate-dataset analysis/1-extract/dataset_definition.py + --output output/1-extract/extract_covid_wave_1.arrow + --dummy-data-file analysis/1-extract/dummy-data/dummy_covid_wave_1.arrow + -- + --cohort_id covid_wave_1 outputs: highly_sensitive: - cohort: output/1-extract/extract_fixed.arrow + cohort: output/1-extract/extract_covid_wave_1.arrow - extract_varying: - run: ehrql:v1 generate-dataset analysis/1-extract/dataset_definition_varying.py - --output output/1-extract/extract_varying:arrow - #--dummy-data-file analysis/1-extract/dummy-data/dummy_varying.arrow + prepare_covid_wave_1: + run: r:v2 analysis/2-prepare/prepare.R covid_wave_1 + needs: [extract_covid_wave_1] outputs: highly_sensitive: - cohort: output/1-extract/extract_varying/*.arrow + cohort: output/2-prepare/prepare_covid_wave_1/prepare_covid_wave_1.arrow + moderately_sensitive: + txt: output/2-prepare/prepare_covid_wave_1/*.txt - vax_data_quality: - run: r:v2 analysis/2-prepare/vax_data_quality.R - needs: [extract_varying] + report_covid_wave_1: + run: r:v2 analysis/3-report/report_cohort.R covid_wave_1 + needs: [prepare_covid_wave_1] outputs: moderately_sensitive: - csv: output/2-prepare/vax_data_quality/*.csv - #txt: output/2-prepare/vax_data_quality/*.txt + csv: output/3-report/report_covid_wave_1/*.csv + png: output/3-report/report_covid_wave_1/*.png + txt: output/3-report/report_covid_wave_1/*.txt + + +# covid: ba1_ba2 + extract_covid_ba1_ba2: + run: ehrql:v1 generate-dataset analysis/1-extract/dataset_definition.py + --output output/1-extract/extract_covid_ba1_ba2.arrow + --dummy-data-file analysis/1-extract/dummy-data/dummy_covid_ba1_ba2.arrow + -- + --cohort_id covid_ba1_ba2 + outputs: + highly_sensitive: + cohort: output/1-extract/extract_covid_ba1_ba2.arrow - prepare: - run: r:v2 analysis/2-prepare/prepare.R - needs: [extract_fixed, extract_varying] + prepare_covid_ba1_ba2: + run: r:v2 analysis/2-prepare/prepare.R covid_ba1_ba2 + needs: [extract_covid_ba1_ba2] outputs: highly_sensitive: - arrow: output/2-prepare/prepare/*.arrow + cohort: output/2-prepare/prepare_covid_ba1_ba2/prepare_covid_ba1_ba2.arrow moderately_sensitive: - #csv: output/2-prepare/prepare/*.csv - txt: output/2-prepare/prepare/*.txt + txt: output/2-prepare/prepare_covid_ba1_ba2/*.txt -# report vaccine history over entire observational period --- - report_history: - run: r:v2 analysis/3-history/report_history.R - needs: [prepare] + report_covid_ba1_ba2: + run: r:v2 analysis/3-report/report_cohort.R covid_ba1_ba2 + needs: [prepare_covid_ba1_ba2] outputs: moderately_sensitive: - csv: output/3-history/report_history/*.csv - png: output/3-history/report_history/*.png + csv: output/3-report/report_covid_ba1_ba2/*.csv + png: output/3-report/report_covid_ba1_ba2/*.png + txt: output/3-report/report_covid_ba1_ba2/*.txt - # report detailed vaccine history at specific snapshots in time ---- - # 2020-12-07 - extract_snapshot_20201207: - run: ehrql:v1 generate-dataset analysis/1-extract/dataset_definition_snapshot.py - --output output/1-extract/extract_snapshot_20201207.arrow - --dummy-data-file analysis/1-extract/dummy-data/dummy_snapshot_20201207.arrow +# covid: jn1 + extract_covid_jn1: + run: ehrql:v1 generate-dataset analysis/1-extract/dataset_definition.py + --output output/1-extract/extract_covid_jn1.arrow + --dummy-data-file analysis/1-extract/dummy-data/dummy_covid_jn1.arrow -- - --snapshot_date 20201207 + --cohort_id covid_jn1 outputs: highly_sensitive: - cohort: output/1-extract/extract_snapshot_20201207.arrow - - report_snapshot_20201207: - run: r:v2 analysis/4-snapshot/report_snapshot.R 20201207 - needs: [extract_snapshot_20201207, prepare] - outputs: - moderately_sensitive: - csv: output/4-snapshot/report_snapshot_20201207/*.csv - png: output/4-snapshot/report_snapshot_20201207/*.png - txt: output/4-snapshot/report_snapshot_20201207/*.txt + cohort: output/1-extract/extract_covid_jn1.arrow - # 2021-09-06 - extract_snapshot_20210906: - run: ehrql:v1 generate-dataset analysis/1-extract/dataset_definition_snapshot.py - --output output/1-extract/extract_snapshot_20210906.arrow - --dummy-data-file analysis/1-extract/dummy-data/dummy_snapshot_20210906.arrow - -- - --snapshot_date 20210906 + prepare_covid_jn1: + run: r:v2 analysis/2-prepare/prepare.R covid_jn1 + needs: [extract_covid_jn1] outputs: highly_sensitive: - cohort: output/1-extract/extract_snapshot_20210906.arrow + cohort: output/2-prepare/prepare_covid_jn1/prepare_covid_jn1.arrow + moderately_sensitive: + txt: output/2-prepare/prepare_covid_jn1/*.txt - report_snapshot_20210906: - run: r:v2 analysis/4-snapshot/report_snapshot.R 20210906 - needs: [extract_snapshot_20210906, prepare] + report_covid_jn1: + run: r:v2 analysis/3-report/report_cohort.R covid_jn1 + needs: [prepare_covid_jn1] outputs: moderately_sensitive: - csv: output/4-snapshot/report_snapshot_20210906/*.csv - png: output/4-snapshot/report_snapshot_20210906/*.png - txt: output/4-snapshot/report_snapshot_20210906/*.txt + csv: output/3-report/report_covid_jn1/*.csv + png: output/3-report/report_covid_jn1/*.png + txt: output/3-report/report_covid_jn1/*.txt - # 2022-03-21 - extract_snapshot_20220321: - run: ehrql:v1 generate-dataset analysis/1-extract/dataset_definition_snapshot.py - --output output/1-extract/extract_snapshot_20220321.arrow - --dummy-data-file analysis/1-extract/dummy-data/dummy_snapshot_20220321.arrow + +# covid: 2025_26 + extract_covid_2025_26: + run: ehrql:v1 generate-dataset analysis/1-extract/dataset_definition.py + --output output/1-extract/extract_covid_2025_26.arrow + --dummy-data-file analysis/1-extract/dummy-data/dummy_covid_2025_26.arrow -- - --snapshot_date 20220321 + --cohort_id covid_2025_26 outputs: highly_sensitive: - cohort: output/1-extract/extract_snapshot_20220321.arrow - - report_snapshot_20220321: - run: r:v2 analysis/4-snapshot/report_snapshot.R 20220321 - needs: [extract_snapshot_20220321, prepare] - outputs: - moderately_sensitive: - csv: output/4-snapshot/report_snapshot_20220321/*.csv - png: output/4-snapshot/report_snapshot_20220321/*.png - txt: output/4-snapshot/report_snapshot_20220321/*.txt + cohort: output/1-extract/extract_covid_2025_26.arrow - # 2022-08-29 - extract_snapshot_20220829: - run: ehrql:v1 generate-dataset analysis/1-extract/dataset_definition_snapshot.py - --output output/1-extract/extract_snapshot_20220829.arrow - --dummy-data-file analysis/1-extract/dummy-data/dummy_snapshot_20220829.arrow - -- - --snapshot_date 20220829 + prepare_covid_2025_26: + run: r:v2 analysis/2-prepare/prepare.R covid_2025_26 + needs: [extract_covid_2025_26] outputs: highly_sensitive: - cohort: output/1-extract/extract_snapshot_20220829.arrow + cohort: output/2-prepare/prepare_covid_2025_26/prepare_covid_2025_26.arrow + moderately_sensitive: + txt: output/2-prepare/prepare_covid_2025_26/*.txt - report_snapshot_20220829: - run: r:v2 analysis/4-snapshot/report_snapshot.R 20220829 - needs: [extract_snapshot_20220829, prepare] + report_covid_2025_26: + run: r:v2 analysis/3-report/report_cohort.R covid_2025_26 + needs: [prepare_covid_2025_26] outputs: moderately_sensitive: - csv: output/4-snapshot/report_snapshot_20220829/*.csv - png: output/4-snapshot/report_snapshot_20220829/*.png - txt: output/4-snapshot/report_snapshot_20220829/*.txt + csv: output/3-report/report_covid_2025_26/*.csv + png: output/3-report/report_covid_2025_26/*.png + txt: output/3-report/report_covid_2025_26/*.txt - # 2023-04-03 - extract_snapshot_20230403: - run: ehrql:v1 generate-dataset analysis/1-extract/dataset_definition_snapshot.py - --output output/1-extract/extract_snapshot_20230403.arrow - --dummy-data-file analysis/1-extract/dummy-data/dummy_snapshot_20230403.arrow + +# flu: 2023_24 + extract_flu_2023_24: + run: ehrql:v1 generate-dataset analysis/1-extract/dataset_definition.py + --output output/1-extract/extract_flu_2023_24.arrow + --dummy-data-file analysis/1-extract/dummy-data/dummy_flu_2023_24.arrow -- - --snapshot_date 20230403 + --cohort_id flu_2023_24 + outputs: + highly_sensitive: + cohort: output/1-extract/extract_flu_2023_24.arrow + + prepare_flu_2023_24: + run: r:v2 analysis/2-prepare/prepare.R flu_2023_24 + needs: [extract_flu_2023_24] outputs: highly_sensitive: - cohort: output/1-extract/extract_snapshot_20230403.arrow + cohort: output/2-prepare/prepare_flu_2023_24/prepare_flu_2023_24.arrow + moderately_sensitive: + txt: output/2-prepare/prepare_flu_2023_24/*.txt - report_snapshot_20230403: - run: r:v2 analysis/4-snapshot/report_snapshot.R 20230403 - needs: [extract_snapshot_20230403, prepare] + report_flu_2023_24: + run: r:v2 analysis/3-report/report_cohort.R flu_2023_24 + needs: [prepare_flu_2023_24] outputs: moderately_sensitive: - csv: output/4-snapshot/report_snapshot_20230403/*.csv - png: output/4-snapshot/report_snapshot_20230403/*.png - txt: output/4-snapshot/report_snapshot_20230403/*.txt + csv: output/3-report/report_flu_2023_24/*.csv + png: output/3-report/report_flu_2023_24/*.png + txt: output/3-report/report_flu_2023_24/*.txt + - # 2023-08-28 - extract_snapshot_20230828: - run: ehrql:v1 generate-dataset analysis/1-extract/dataset_definition_snapshot.py - --output output/1-extract/extract_snapshot_20230828.arrow - --dummy-data-file analysis/1-extract/dummy-data/dummy_snapshot_20230828.arrow +# flu: 2024_25 + extract_flu_2024_25: + run: ehrql:v1 generate-dataset analysis/1-extract/dataset_definition.py + --output output/1-extract/extract_flu_2024_25.arrow + --dummy-data-file analysis/1-extract/dummy-data/dummy_flu_2024_25.arrow -- - --snapshot_date 20230828 + --cohort_id flu_2024_25 outputs: highly_sensitive: - cohort: output/1-extract/extract_snapshot_20230828.arrow + cohort: output/1-extract/extract_flu_2024_25.arrow - report_snapshot_20230828: - run: r:v2 analysis/4-snapshot/report_snapshot.R 20230828 - needs: [extract_snapshot_20230828, prepare] + prepare_flu_2024_25: + run: r:v2 analysis/2-prepare/prepare.R flu_2024_25 + needs: [extract_flu_2024_25] + outputs: + highly_sensitive: + cohort: output/2-prepare/prepare_flu_2024_25/prepare_flu_2024_25.arrow + moderately_sensitive: + txt: output/2-prepare/prepare_flu_2024_25/*.txt + + report_flu_2024_25: + run: r:v2 analysis/3-report/report_cohort.R flu_2024_25 + needs: [prepare_flu_2024_25] outputs: moderately_sensitive: - csv: output/4-snapshot/report_snapshot_20230828/*.csv - png: output/4-snapshot/report_snapshot_20230828/*.png - txt: output/4-snapshot/report_snapshot_20230828/*.txt + csv: output/3-report/report_flu_2024_25/*.csv + png: output/3-report/report_flu_2024_25/*.png + txt: output/3-report/report_flu_2024_25/*.txt + - # 2024-04-15 - extract_snapshot_20240415: - run: ehrql:v1 generate-dataset analysis/1-extract/dataset_definition_snapshot.py - --output output/1-extract/extract_snapshot_20240415.arrow - --dummy-data-file analysis/1-extract/dummy-data/dummy_snapshot_20240415.arrow +# flu: 2025_26 + extract_flu_2025_26: + run: ehrql:v1 generate-dataset analysis/1-extract/dataset_definition.py + --output output/1-extract/extract_flu_2025_26.arrow + --dummy-data-file analysis/1-extract/dummy-data/dummy_flu_2025_26.arrow -- - --snapshot_date 20240415 + --cohort_id flu_2025_26 + outputs: + highly_sensitive: + cohort: output/1-extract/extract_flu_2025_26.arrow + + prepare_flu_2025_26: + run: r:v2 analysis/2-prepare/prepare.R flu_2025_26 + needs: [extract_flu_2025_26] outputs: highly_sensitive: - cohort: output/1-extract/extract_snapshot_20240415.arrow + cohort: output/2-prepare/prepare_flu_2025_26/prepare_flu_2025_26.arrow + moderately_sensitive: + txt: output/2-prepare/prepare_flu_2025_26/*.txt - report_snapshot_20240415: - run: r:v2 analysis/4-snapshot/report_snapshot.R 20240415 - needs: [extract_snapshot_20240415, prepare] + report_flu_2025_26: + run: r:v2 analysis/3-report/report_cohort.R flu_2025_26 + needs: [prepare_flu_2025_26] outputs: moderately_sensitive: - csv: output/4-snapshot/report_snapshot_20240415/*.csv - png: output/4-snapshot/report_snapshot_20240415/*.png - txt: output/4-snapshot/report_snapshot_20240415/*.txt + csv: output/3-report/report_flu_2025_26/*.csv + png: output/3-report/report_flu_2025_26/*.png + txt: output/3-report/report_flu_2025_26/*.txt - # 2024-09-30 - extract_snapshot_20240930: - run: ehrql:v1 generate-dataset analysis/1-extract/dataset_definition_snapshot.py - --output output/1-extract/extract_snapshot_20240930.arrow - --dummy-data-file analysis/1-extract/dummy-data/dummy_snapshot_20240930.arrow + +# rsv: 2023_24 + extract_rsv_2023_24: + run: ehrql:v1 generate-dataset analysis/1-extract/dataset_definition.py + --output output/1-extract/extract_rsv_2023_24.arrow + --dummy-data-file analysis/1-extract/dummy-data/dummy_rsv_2023_24.arrow -- - --snapshot_date 20240930 + --cohort_id rsv_2023_24 outputs: highly_sensitive: - cohort: output/1-extract/extract_snapshot_20240930.arrow + cohort: output/1-extract/extract_rsv_2023_24.arrow - report_snapshot_20240930: - run: r:v2 analysis/4-snapshot/report_snapshot.R 20240930 - needs: [extract_snapshot_20240930, prepare] + prepare_rsv_2023_24: + run: r:v2 analysis/2-prepare/prepare.R rsv_2023_24 + needs: [extract_rsv_2023_24] outputs: + highly_sensitive: + cohort: output/2-prepare/prepare_rsv_2023_24/prepare_rsv_2023_24.arrow moderately_sensitive: - csv: output/4-snapshot/report_snapshot_20240930/*.csv - png: output/4-snapshot/report_snapshot_20240930/*.png - txt: output/4-snapshot/report_snapshot_20240930/*.txt + txt: output/2-prepare/prepare_rsv_2023_24/*.txt - # 2025-03-31 - extract_snapshot_20250331: - run: ehrql:v1 generate-dataset analysis/1-extract/dataset_definition_snapshot.py - --output output/1-extract/extract_snapshot_20250331.arrow - --dummy-data-file analysis/1-extract/dummy-data/dummy_snapshot_20250331.arrow + report_rsv_2023_24: + run: r:v2 analysis/3-report/report_cohort.R rsv_2023_24 + needs: [prepare_rsv_2023_24] + outputs: + moderately_sensitive: + csv: output/3-report/report_rsv_2023_24/*.csv + png: output/3-report/report_rsv_2023_24/*.png + txt: output/3-report/report_rsv_2023_24/*.txt + + +# rsv: 2024_25 + extract_rsv_2024_25: + run: ehrql:v1 generate-dataset analysis/1-extract/dataset_definition.py + --output output/1-extract/extract_rsv_2024_25.arrow + --dummy-data-file analysis/1-extract/dummy-data/dummy_rsv_2024_25.arrow -- - --snapshot_date 20250331 + --cohort_id rsv_2024_25 outputs: highly_sensitive: - cohort: output/1-extract/extract_snapshot_20250331.arrow + cohort: output/1-extract/extract_rsv_2024_25.arrow - report_snapshot_20250331: - run: r:v2 analysis/4-snapshot/report_snapshot.R 20250331 - needs: [extract_snapshot_20250331, prepare] + prepare_rsv_2024_25: + run: r:v2 analysis/2-prepare/prepare.R rsv_2024_25 + needs: [extract_rsv_2024_25] outputs: + highly_sensitive: + cohort: output/2-prepare/prepare_rsv_2024_25/prepare_rsv_2024_25.arrow moderately_sensitive: - csv: output/4-snapshot/report_snapshot_20250331/*.csv - png: output/4-snapshot/report_snapshot_20250331/*.png - txt: output/4-snapshot/report_snapshot_20250331/*.txt + txt: output/2-prepare/prepare_rsv_2024_25/*.txt -# 2025-09-29 + report_rsv_2024_25: + run: r:v2 analysis/3-report/report_cohort.R rsv_2024_25 + needs: [prepare_rsv_2024_25] + outputs: + moderately_sensitive: + csv: output/3-report/report_rsv_2024_25/*.csv + png: output/3-report/report_rsv_2024_25/*.png + txt: output/3-report/report_rsv_2024_25/*.txt - extract_snapshot_20250929: - run: ehrql:v1 generate-dataset analysis/1-extract/dataset_definition_snapshot.py - --output output/1-extract/extract_snapshot_20250929.arrow - --dummy-data-file analysis/1-extract/dummy-data/dummy_snapshot_20250929.arrow + +# rsv: 2025_26 + extract_rsv_2025_26: + run: ehrql:v1 generate-dataset analysis/1-extract/dataset_definition.py + --output output/1-extract/extract_rsv_2025_26.arrow + --dummy-data-file analysis/1-extract/dummy-data/dummy_rsv_2025_26.arrow -- - --snapshot_date 20250929 + --cohort_id rsv_2025_26 outputs: highly_sensitive: - cohort: output/1-extract/extract_snapshot_20250929.arrow + cohort: output/1-extract/extract_rsv_2025_26.arrow + + prepare_rsv_2025_26: + run: r:v2 analysis/2-prepare/prepare.R rsv_2025_26 + needs: [extract_rsv_2025_26] + outputs: + highly_sensitive: + cohort: output/2-prepare/prepare_rsv_2025_26/prepare_rsv_2025_26.arrow + moderately_sensitive: + txt: output/2-prepare/prepare_rsv_2025_26/*.txt - report_snapshot_20250929: - run: r:v2 analysis/4-snapshot/report_snapshot.R 20250929 - needs: [extract_snapshot_20250929, prepare] + report_rsv_2025_26: + run: r:v2 analysis/3-report/report_cohort.R rsv_2025_26 + needs: [prepare_rsv_2025_26] outputs: moderately_sensitive: - csv: output/4-snapshot/report_snapshot_20250929/*.csv - png: output/4-snapshot/report_snapshot_20250929/*.png - txt: output/4-snapshot/report_snapshot_20250929/*.txt + csv: output/3-report/report_rsv_2025_26/*.csv + png: output/3-report/report_rsv_2025_26/*.png + txt: output/3-report/report_rsv_2025_26/*.txt \ No newline at end of file