From 0eb74767d65e590d1547bc28f751478c3cec5ccc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rasmus=20Skytte=20Randl=C3=B8v?= Date: Fri, 12 Jan 2024 14:50:22 +0100 Subject: [PATCH 01/13] docs(popIBM): add additional code context --- manuscripts/popIBM/popIBM.R | 47 +++++++++++++++++++++++++++++-------- 1 file changed, 37 insertions(+), 10 deletions(-) diff --git a/manuscripts/popIBM/popIBM.R b/manuscripts/popIBM/popIBM.R index 04fce6d0..ab964e5a 100644 --- a/manuscripts/popIBM/popIBM.R +++ b/manuscripts/popIBM/popIBM.R @@ -226,8 +226,14 @@ sim_list <- foreach(run_this = (first_run - 1 + 1:n_runs), .packages = "data.tab times <- seq(start_denmark, end_times, 1) xdates <- as.Date(times, origin = "2020-01-01") - # Initialise spatial heterogeneity in parishes - + # Initialise spatial heterogeneity in parishes. + # During initialisation, rel_risk_parish is set to the observed cumulative incidence from begining of epidemic until 26 April 2021. + # However, distribution of risks .... + # + # For reference, see + # Græsbøll, K., Eriksen, R. S., Kirkeby, C., & Christiansen, L. E., & (2023). + # Mass testing and local lockdowns effectively controls COVID-19. Manuscript + # accepted for publication in Communications Medicine ibm[, rel_risk_parish := rel_risk_parish^(1 / 3)] ibm[, rel_risk_parish := rel_risk_parish * .N / sum(rel_risk_parish)] @@ -314,10 +320,14 @@ sim_list <- foreach(run_this = (first_run - 1 + 1:n_runs), .packages = "data.tab # The compute the probability of test when adjusting for the number of tests p_test_corr <- p_test_inc(inc) - if (day <= day_fix_p_test) { # The number of tests is known + if (day <= day_fix_p_test) { # The number of tests is known until this date # Group population by age group and vaccinated / un-vaccinated tmp <- ibm[, .N, keyby = .(age_groups, !(vac_time < breaks_vac[[1]] | is.na(vac_time)))] + + # Then merge population counts with testing behavior in age and vaccine groups (n_test_age_groups_int_vac). + # This gives the ??? + # (n_test_age_groups_int_vac is loaded as part of init) t_pop_age_vac <- tmp[n_test_age_groups_int_vac, , on = c("age_groups", "vac_time")] t_pop_age_vac[is.na(N), N := 0] @@ -435,8 +445,10 @@ sim_list <- foreach(run_this = (first_run - 1 + 1:n_runs), .packages = "data.tab # Recover from disease I -> R ibm[disease == 2L & tt == 0, disease := 3L] - # When all age groups have same disease progression E-> I - # Also draw time to being symptomatic + # All age groups have same disease progression E-> I (disease states 1L -> 2L). + # Upon entering the I state, the time to recovery (tt) and time to symptoms (tt_symp) are drawn. + # Around 50% of infected people don't develop symptoms to the degree that they take a test. + # Therefore, we set tt_symp to NA for 50% of infected. ibm[ disease == 1L & tt == 0, `:=`( @@ -447,7 +459,8 @@ sim_list <- foreach(run_this = (first_run - 1 + 1:n_runs), .packages = "data.tab ) ] - # Do not double count people found in E states + # Persons already isolated (non_isolated == 0) are assumed to not test themselves again once they develop symptoms. + # We therefore set tt_symp to -1L to prevent such tests. ibm[disease == 2L & non_isolated == 0L & tt_symp > 0, tt_symp := -1L] @@ -457,20 +470,34 @@ sim_list <- foreach(run_this = (first_run - 1 + 1:n_runs), .packages = "data.tab # Lockdown i_lock <- sum(day > day_lockdown_change) - # Parish - n_cases <- colSums(sim_parish[(day - 7):(day - 1), ], na.rm = TRUE) + # Compute and store the incidences for the local lockdown rules. + # 1) Compute 7-day incidences per 100k leading up to current day + # 2) Store the incidence in the history tracker + # 3) Compute maximum 7-day incidence over the last 7 days + # 4) Apply the current lockdown thresholds given this max incidence + + ## Parish + n_7d_cases_parish <- colSums(sim_parish[(day - 7):(day - 1), ], na.rm = TRUE) - inc_his_parish[day, ] <- (n_cases >= 20) * n_cases / pop_parish * 1e5 + # NOTE: for parishes, local lockdown required at least 20 positive cases in the parish. + # We therefore record incidences as 0 in these cases. + inc_his_parish[day, ] <- (n_7d_cases_parish >= 20) * n_7d_cases_parish / pop_parish * 1e5 max_7d_inc <- apply(inc_his_parish[(day - 6):day, ], 2, max, na.rm = TRUE) + lockdown_parish_fac <- lockdown_parish_fun[[i_lock]](max_7d_inc) - # Municipality + + ## Municipality inc_his_municipality[day, ] <- sim_municipality[(day - 7):(day - 1), ] |> colSums(na.rm = TRUE) / pop_municipality * 1e5 + max_7d_inc <- apply(inc_his_municipality[(day - 6):day, ], 2, max, na.rm = TRUE) + lockdown_municipality_fac <- lockdown_municipality_fun[[i_lock]](max_7d_inc) + + population[ data.table(parish_id = u_parish_ids, lockdown_parish_fac), on = "parish_id", From e867fc3a9490b402b255b33c2d1500b28b598cfc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rasmus=20Skytte=20Randl=C3=B8v?= Date: Fri, 12 Jan 2024 14:51:52 +0100 Subject: [PATCH 02/13] docs(popIBM): rename day_lockdown_change to day_limit_change --- manuscripts/popIBM/popIBM.R | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/manuscripts/popIBM/popIBM.R b/manuscripts/popIBM/popIBM.R index ab964e5a..9060863d 100644 --- a/manuscripts/popIBM/popIBM.R +++ b/manuscripts/popIBM/popIBM.R @@ -126,8 +126,8 @@ day_restriction_change <- as.numeric(c(activity_scenario$list_beta_dates) - as. list_beta <- activity_scenario$list_beta # Days of changing incidence limits for imposing local lockdown (relative to start date) -day_lockdown_change <- c("2021-03-01", "2021-04-30", "2021-05-28", "2021-07-16", "2021-09-10", "2021-11-15") -day_lockdown_change <- as.numeric(as.Date(day_lockdown_change) - as.Date("2020-01-01")) - start_denmark +day_limit_change <- c("2021-03-01", "2021-04-30", "2021-05-28", "2021-07-16", "2021-09-10", "2021-11-15") +day_limit_change <- as.numeric(as.Date(day_limit_change) - as.Date("2020-01-01")) - start_denmark # Functions for lockdown (corresponding to Danish policy): lockdown_parish_fun <- list( @@ -465,10 +465,10 @@ sim_list <- foreach(run_this = (first_run - 1 + 1:n_runs), .packages = "data.tab # Implement the effects of local lockdown - if (activate_lockdown && day > day_lockdown_change[1]) { + if (activate_lockdown && day > day_limit_change[1]) { # Lockdown - i_lock <- sum(day > day_lockdown_change) + i_lock <- sum(day > day_limit_change) # Compute and store the incidences for the local lockdown rules. # 1) Compute 7-day incidences per 100k leading up to current day From b5aecfb20ccb1280ecc70b8fb14a3d492d8b6bfd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rasmus=20Skytte=20Randl=C3=B8v?= Date: Mon, 15 Jan 2024 10:01:18 +0100 Subject: [PATCH 03/13] docs(popIBM): revert comment --- manuscripts/popIBM/popIBM.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/manuscripts/popIBM/popIBM.R b/manuscripts/popIBM/popIBM.R index 9060863d..cef394db 100644 --- a/manuscripts/popIBM/popIBM.R +++ b/manuscripts/popIBM/popIBM.R @@ -320,7 +320,7 @@ sim_list <- foreach(run_this = (first_run - 1 + 1:n_runs), .packages = "data.tab # The compute the probability of test when adjusting for the number of tests p_test_corr <- p_test_inc(inc) - if (day <= day_fix_p_test) { # The number of tests is known until this date + if (day <= day_fix_p_test) { # The number of tests is known # Group population by age group and vaccinated / un-vaccinated tmp <- ibm[, .N, keyby = .(age_groups, !(vac_time < breaks_vac[[1]] | is.na(vac_time)))] From a0a67bce6cd0b781386dbc96feaa8fe78750bfdc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rasmus=20Skytte=20Randl=C3=B8v?= Date: Mon, 15 Jan 2024 10:10:24 +0100 Subject: [PATCH 04/13] docs(popIBM): use long variable names in for loops --- manuscripts/popIBM/popIBM.R | 50 +++++++++++++++++++------------------ 1 file changed, 26 insertions(+), 24 deletions(-) diff --git a/manuscripts/popIBM/popIBM.R b/manuscripts/popIBM/popIBM.R index cef394db..815b7816 100644 --- a/manuscripts/popIBM/popIBM.R +++ b/manuscripts/popIBM/popIBM.R @@ -208,8 +208,8 @@ tic <- Sys.time() sim_list <- foreach(run_this = (first_run - 1 + 1:n_runs), .packages = "data.table", .verbose = TRUE) %dopar% { scenario <- unlist(sce_combi[run_this, ]) - for (i in seq_along(scenario)) { - assign(x = names(scenario)[i], value = scenario[i]) + for (scenario_id in seq_along(scenario)) { + assign(x = names(scenario)[scenario_id], value = scenario[scenario_id]) } cat("\t run: ", run_this, " ") @@ -386,21 +386,22 @@ sim_list <- foreach(run_this = (first_run - 1 + 1:n_runs), .packages = "data.tab # Collect data on the number of test positives each day - by variant, age and vaccination status - for (k in 1:n_variants) { + for (variant_id in 1:n_variants) { - sim_test_positive[day, , k] <- ibm[ - tt_symp == 0L & variant == k, .N, by = .(age_groups) + sim_test_positive[day, , variant_id] <- ibm[ + tt_symp == 0L & variant == variant_id, .N, by = .(age_groups) ][.(age_groups = 1:9), on = "age_groups"]$N - sim_test_positive_vac[day, , k, 1] <- ibm[ - tt_symp == 0L & variant == k & (vac_time < breaks_vac[1] | is.na(vac_time)), + sim_test_positive_vac[day, , variant_id, 1] <- ibm[ + tt_symp == 0L & variant == variant_id & (vac_time < breaks_vac[1] | is.na(vac_time)), .N, by = .(age_groups) ][.(age_groups = 1:9), on = "age_groups"]$N - for (kk in 2:n_vac_groups_out) { - sim_test_positive_vac[day, , k, kk] <- ibm[ - tt_symp == 0L & variant == k & vac_time >= breaks_vac[kk - 1] & vac_time < breaks_vac[kk], .N, + for (vac_id in 2:n_vac_groups_out) { + sim_test_positive_vac[day, , variant_id, vac_id] <- ibm[ + tt_symp == 0L & variant == variant_id & vac_time >= breaks_vac[vac_id - 1] & vac_time < breaks_vac[vac_id], + .N, by = .(age_groups) ][.(age_groups = 1:9), on = "age_groups"]$N } @@ -420,19 +421,20 @@ sim_list <- foreach(run_this = (first_run - 1 + 1:n_runs), .packages = "data.tab # Collect probability of hospitalisation each day - by variant, age and vaccination status - for (k in 1:n_variants) { - sim_hospital[day, , k] <- ibm[disease == 2L & tt == 0 & variant == k, sum(prob_hospital), + for (variant_id in 1:n_variants) { + sim_hospital[day, , variant_id] <- ibm[disease == 2L & tt == 0 & variant == variant_id, sum(prob_hospital), by = .(age_groups)][.(age_groups = 1:9), on = "age_groups"]$V1 - sim_hospital_vac[day, , k, 1] <- ibm[ - disease == 2L & tt == 0 & variant == k & (vac_time < breaks_vac[1] | is.na(vac_time)), + sim_hospital_vac[day, , variant_id, 1] <- ibm[ + disease == 2L & tt == 0 & variant == variant_id & (vac_time < breaks_vac[1] | is.na(vac_time)), sum(prob_hospital), by = .(age_groups) ][.(age_groups = 1:9), on = "age_groups"]$V1 - for (kk in 2:n_vac_groups_out) { - sim_hospital_vac[day, , k, kk] <- ibm[ - disease == 2L & tt == 0 & variant == k & vac_time >= breaks_vac[kk - 1] & vac_time < breaks_vac[kk], + for (vac_id in 2:n_vac_groups_out) { + sim_hospital_vac[day, , variant_id, vac_id] <- ibm[ + disease == 2L & tt == 0 & variant == variant_id & + vac_time >= breaks_vac[vac_id - 1] & vac_time < breaks_vac[vac_id], sum(prob_hospital), by = .(age_groups) ][.(age_groups = 1:9), on = "age_groups"]$V1 @@ -521,17 +523,17 @@ sim_list <- foreach(run_this = (first_run - 1 + 1:n_runs), .packages = "data.tab } # Infected individuals with different strains - for (k in 1:n_variants) { + for (variant_id in 1:n_variants) { # Calculate the infection pressure inf_persons_municipality <- ibm[ - disease == 2L & variant == k, + disease == 2L & variant == variant_id, .(inf_persons = sum(lockdown_fac * non_isolated * vac_fac_trans)), by = .(municipality_id, age_groups) ][all_pop_combi, on = c("municipality_id", "age_groups")] inf_persons_municipality[is.na(inf_persons), inf_persons := 0] - inf_persons_municipality[, inf_persons := inf_persons * v_rel_beta[k]] + inf_persons_municipality[, inf_persons := inf_persons * v_rel_beta[variant_id]] inf_pressure <- inf_persons_municipality[ , current_beta %*% inf_persons, by = .(municipality_id) @@ -571,7 +573,7 @@ sim_list <- foreach(run_this = (first_run - 1 + 1:n_runs), .packages = "data.tab `:=`(disease = 1L, tt = pmax(1L, round(rgamma(n = .N, shape = v_shape_e[age_groups], scale = v_scale_e[age_groups]))), - variant = k)] + variant = variant_id)] } # Count down to change in disease state or symptoms @@ -583,10 +585,10 @@ sim_list <- foreach(run_this = (first_run - 1 + 1:n_runs), .packages = "data.tab # Implement the effect of vaccination - for (k in 1:n__vac){ + for (vac_id in 1:n_vac){ - ibm[vac_type == k & vac_time == v_vac_tt_effect[k], - `:=`(vacF_ec = delta_vac_effect[k], + ibm[vac_type == vac_id & vac_time == v_vac_tt_effect[vac_id], + `:=`(vacF_ec = delta_vac_effect[vac_id], prob_hospital = delta_red_hospital_risk * prob_hospital, vac_fac_trans = red_transmission_vac)] From 0ef368d6166550140cfaef2f8628eb40a4887ab6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rasmus=20Skytte=20Randl=C3=B8v?= Date: Mon, 15 Jan 2024 10:12:54 +0100 Subject: [PATCH 05/13] chore(popIBM): align line breaks between loops --- manuscripts/popIBM/popIBM.R | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/manuscripts/popIBM/popIBM.R b/manuscripts/popIBM/popIBM.R index 815b7816..db40cdc1 100644 --- a/manuscripts/popIBM/popIBM.R +++ b/manuscripts/popIBM/popIBM.R @@ -422,8 +422,11 @@ sim_list <- foreach(run_this = (first_run - 1 + 1:n_runs), .packages = "data.tab # Collect probability of hospitalisation each day - by variant, age and vaccination status for (variant_id in 1:n_variants) { - sim_hospital[day, , variant_id] <- ibm[disease == 2L & tt == 0 & variant == variant_id, sum(prob_hospital), - by = .(age_groups)][.(age_groups = 1:9), on = "age_groups"]$V1 + sim_hospital[day, , variant_id] <- ibm[ + disease == 2L & tt == 0 & variant == variant_id, + sum(prob_hospital), + by = .(age_groups) + ][.(age_groups = 1:9), on = "age_groups"]$V1 sim_hospital_vac[day, , variant_id, 1] <- ibm[ disease == 2L & tt == 0 & variant == variant_id & (vac_time < breaks_vac[1] | is.na(vac_time)), From 8fdbb9eac61ab4dfb0313e9449f4f290744bcd14 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rasmus=20Skytte=20Randl=C3=B8v?= Date: Mon, 15 Jan 2024 10:16:30 +0100 Subject: [PATCH 06/13] chore(popIBM): white space fix --- manuscripts/popIBM/popIBM.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/manuscripts/popIBM/popIBM.R b/manuscripts/popIBM/popIBM.R index db40cdc1..0110a55d 100644 --- a/manuscripts/popIBM/popIBM.R +++ b/manuscripts/popIBM/popIBM.R @@ -450,7 +450,7 @@ sim_list <- foreach(run_this = (first_run - 1 + 1:n_runs), .packages = "data.tab # Recover from disease I -> R ibm[disease == 2L & tt == 0, disease := 3L] - # All age groups have same disease progression E-> I (disease states 1L -> 2L). + # All age groups have same disease progression E-> I (disease states 1L -> 2L). # Upon entering the I state, the time to recovery (tt) and time to symptoms (tt_symp) are drawn. # Around 50% of infected people don't develop symptoms to the degree that they take a test. # Therefore, we set tt_symp to NA for 50% of infected. @@ -475,7 +475,7 @@ sim_list <- foreach(run_this = (first_run - 1 + 1:n_runs), .packages = "data.tab # Lockdown i_lock <- sum(day > day_limit_change) - # Compute and store the incidences for the local lockdown rules. + # Compute and store the incidences for the local lockdown rules. # 1) Compute 7-day incidences per 100k leading up to current day # 2) Store the incidence in the history tracker # 3) Compute maximum 7-day incidence over the last 7 days From bbe2efa858c692fda91271815ca6f7cc16b561f9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rasmus=20Skytte=20Randl=C3=B8v?= Date: Mon, 15 Jan 2024 10:16:58 +0100 Subject: [PATCH 07/13] docs(popIBM): use longer variable names - part 6 --- manuscripts/popIBM/popIBM.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/manuscripts/popIBM/popIBM.R b/manuscripts/popIBM/popIBM.R index 0110a55d..e1d412d9 100644 --- a/manuscripts/popIBM/popIBM.R +++ b/manuscripts/popIBM/popIBM.R @@ -512,10 +512,10 @@ sim_list <- foreach(run_this = (first_run - 1 + 1:n_runs), .packages = "data.tab population[ data.table(municipality_id = u_municipality_ids, lockdown_municipality_fac), on = "municipality_id", - kom_fac := lockdown_municipality_fac + municipality_fac := lockdown_municipality_fac ] - population[, lockdown_max := pmax(parish_fac, kom_fac)] + population[, lockdown_max := pmax(parish_fac, municipality_fac)] # Weighted sum as lockdown factor population[, lockdown_fac := 1 * (1 - lockdown_max) + lockdown_factor * lockdown_max] From 321f8c0b2d65911bbc89d69bb9c343c33366511e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rasmus=20Skytte=20Randl=C3=B8v?= Date: Mon, 15 Jan 2024 10:17:19 +0100 Subject: [PATCH 08/13] docs(popIBM): add context comments --- manuscripts/popIBM/popIBM.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/manuscripts/popIBM/popIBM.R b/manuscripts/popIBM/popIBM.R index e1d412d9..ebb1f907 100644 --- a/manuscripts/popIBM/popIBM.R +++ b/manuscripts/popIBM/popIBM.R @@ -447,7 +447,7 @@ sim_list <- foreach(run_this = (first_run - 1 + 1:n_runs), .packages = "data.tab - # Recover from disease I -> R + # Recover from disease I -> R (disease states 2L -> 3L) ibm[disease == 2L & tt == 0, disease := 3L] # All age groups have same disease progression E-> I (disease states 1L -> 2L). @@ -502,7 +502,7 @@ sim_list <- foreach(run_this = (first_run - 1 + 1:n_runs), .packages = "data.tab lockdown_municipality_fac <- lockdown_municipality_fun[[i_lock]](max_7d_inc) - + # Transfer degree of local lockdown (lockdown_fac) to population tracker population[ data.table(parish_id = u_parish_ids, lockdown_parish_fac), on = "parish_id", From b39fef61f5b96c0a586bc9eb8704d19cf7c5f007 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rasmus=20Skytte=20Randl=C3=B8v?= Date: Mon, 15 Jan 2024 10:22:13 +0100 Subject: [PATCH 09/13] docs(popIBM): add context comments - part 2 --- manuscripts/popIBM/popIBM.R | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/manuscripts/popIBM/popIBM.R b/manuscripts/popIBM/popIBM.R index ebb1f907..68873f88 100644 --- a/manuscripts/popIBM/popIBM.R +++ b/manuscripts/popIBM/popIBM.R @@ -277,11 +277,17 @@ sim_list <- foreach(run_this = (first_run - 1 + 1:n_runs), .packages = "data.tab # Set "beta" based on restriction levels and seasonal change beta_season <- (1 - season_fac * (1 - seasonal_rel_beta(as.Date(start_denmark, origin = "2020-01-01"), day))) - if (day > day_restriction_change[1]) { # Activity different from initial activity (restriction change) + # When activity is different from initial activity (i.e, due to an restriction change) + # We store some helper variables + if (day > day_restriction_change[1]) 0{ + + # Index for current beta matrix i_beta <- max(which(day_restriction_change <= day)) + # Value of current beta current_beta <- beta_season * r_ref * 0.35 * list_beta[[i_beta]] + # Factor needed to reduce activity to the initial level (i.e. full lockdown) lockdown_factor <- sqrt(eigen(list_beta[[1]])$values[1] / eigen(list_beta[[i_beta]])$values[1]) } else { # Initial activity level From 2fcaab3ef094962638ba3fbd32093c1696685cbb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rasmus=20Skytte=20Randl=C3=B8v?= Date: Mon, 15 Jan 2024 10:32:41 +0100 Subject: [PATCH 10/13] docs(popIBM): add trace back info --- manuscripts/popIBM/popIBM.R | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/manuscripts/popIBM/popIBM.R b/manuscripts/popIBM/popIBM.R index 68873f88..4c16c093 100644 --- a/manuscripts/popIBM/popIBM.R +++ b/manuscripts/popIBM/popIBM.R @@ -45,7 +45,7 @@ load_info <- FALSE # Simulation dates new_end_date <- as.Date("2021-06-30") end_times <- as.numeric(new_end_date) - as.numeric(as.Date("2020-01-01")) -times <- seq(start_denmark, end_times, 1) +times <- seq(start_denmark, end_times, 1) # start_denmark is simulation start date (loaded in init file) # Proportion of transmission within municipality w_municipality <- 0.9 @@ -161,6 +161,7 @@ variant_id_delta <- 3 # Variant id for delta # Determine the last day where the number of tests is known # after this day, we fix p_test to the last known value +# ntal is loaded in init file day_fix_p_test <- as.numeric(ntal[, as.Date(max(pr_date))] - as.Date("2020-01-01")) - start_denmark sce_test_red <- 1 # Factor for probability of taking a test @@ -324,6 +325,7 @@ sim_list <- foreach(run_this = (first_run - 1 + 1:n_runs), .packages = "data.tab inc <- colSums(sim_municipality[(day - 7):(day - 1), ], na.rm = TRUE) / pop_municipality * 1e5 # The compute the probability of test when adjusting for the number of tests + # p_test_inc loaded from init p_test_corr <- p_test_inc(inc) if (day <= day_fix_p_test) { # The number of tests is known From ebe4a4486753710147444740415db8ad11673df2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rasmus=20Skytte=20Randl=C3=B8v?= Date: Mon, 15 Jan 2024 10:32:50 +0100 Subject: [PATCH 11/13] chore(popIBM): fix formatting --- manuscripts/popIBM/popIBM.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/manuscripts/popIBM/popIBM.R b/manuscripts/popIBM/popIBM.R index 4c16c093..54756dab 100644 --- a/manuscripts/popIBM/popIBM.R +++ b/manuscripts/popIBM/popIBM.R @@ -280,7 +280,7 @@ sim_list <- foreach(run_this = (first_run - 1 + 1:n_runs), .packages = "data.tab # When activity is different from initial activity (i.e, due to an restriction change) # We store some helper variables - if (day > day_restriction_change[1]) 0{ + if (day > day_restriction_change[1]) { # Index for current beta matrix i_beta <- max(which(day_restriction_change <= day)) @@ -596,7 +596,7 @@ sim_list <- foreach(run_this = (first_run - 1 + 1:n_runs), .packages = "data.tab # Implement the effect of vaccination - for (vac_id in 1:n_vac){ + for (vac_id in 1:n_vac) { ibm[vac_type == vac_id & vac_time == v_vac_tt_effect[vac_id], `:=`(vacF_ec = delta_vac_effect[vac_id], From 5a51f6fe3a778713247cdaecfdd5078a5431e87a Mon Sep 17 00:00:00 2001 From: kaare-gr Date: Mon, 15 Jan 2024 15:44:06 +0100 Subject: [PATCH 12/13] docs(popIBM): add context comment --- manuscripts/popIBM/popIBM.R | 25 +++++++++++++++++-------- 1 file changed, 17 insertions(+), 8 deletions(-) diff --git a/manuscripts/popIBM/popIBM.R b/manuscripts/popIBM/popIBM.R index 54756dab..5a188c0f 100644 --- a/manuscripts/popIBM/popIBM.R +++ b/manuscripts/popIBM/popIBM.R @@ -334,7 +334,7 @@ sim_list <- foreach(run_this = (first_run - 1 + 1:n_runs), .packages = "data.tab tmp <- ibm[, .N, keyby = .(age_groups, !(vac_time < breaks_vac[[1]] | is.na(vac_time)))] # Then merge population counts with testing behavior in age and vaccine groups (n_test_age_groups_int_vac). - # This gives the ??? + # This merge is to be able to make groups with zeroes (.N ignores them) # (n_test_age_groups_int_vac is loaded as part of init) t_pop_age_vac <- tmp[n_test_age_groups_int_vac, , on = c("age_groups", "vac_time")] t_pop_age_vac[is.na(N), N := 0] @@ -523,6 +523,7 @@ sim_list <- foreach(run_this = (first_run - 1 + 1:n_runs), .packages = "data.tab municipality_fac := lockdown_municipality_fac ] + # an individual is affected by the biggest lockdown either in parish or municipality population[, lockdown_max := pmax(parish_fac, municipality_fac)] # Weighted sum as lockdown factor @@ -533,10 +534,10 @@ sim_list <- foreach(run_this = (first_run - 1 + 1:n_runs), .packages = "data.tab } - # Infected individuals with different strains + # Infected individuals with different variants for (variant_id in 1:n_variants) { - # Calculate the infection pressure + # Calculate the number of infected persons per municipality + age_group inf_persons_municipality <- ibm[ disease == 2L & variant == variant_id, .(inf_persons = sum(lockdown_fac * non_isolated * vac_fac_trans)), @@ -546,6 +547,7 @@ sim_list <- foreach(run_this = (first_run - 1 + 1:n_runs), .packages = "data.tab inf_persons_municipality[is.na(inf_persons), inf_persons := 0] inf_persons_municipality[, inf_persons := inf_persons * v_rel_beta[variant_id]] + # Calculate the infection pressure by multiplying betas onto inf_persons inf_pressure <- inf_persons_municipality[ , current_beta %*% inf_persons, by = .(municipality_id) ][, age_groups := rep(1:9, n_municipality)] @@ -555,8 +557,10 @@ sim_list <- foreach(run_this = (first_run - 1 + 1:n_runs), .packages = "data.tab inf_pressure <- dt_pop_municipality[inf_pressure, , on = "municipality_id"] inf_pressure[, r_inf_municipality := r_inf_municipality / pop] + # calculate infected persons nationwide inf_persons <- inf_persons_municipality[, .(inf_persons = sum(inf_persons)), by = .(age_groups)] + # calculate nationwide infection pressure tmp <- inf_persons$inf_persons inf_pressure_dk <- inf_persons[ , @@ -564,19 +568,22 @@ sim_list <- foreach(run_this = (first_run - 1 + 1:n_runs), .packages = "data.tab by = .(age_groups) ] + # weigth nationwide and municipal infection pressure together + # each age_groups and municipality now have a risk of infection inf_pressure_total <- inf_pressure_dk[inf_pressure, , on = "age_groups"] inf_pressure_total[ , prob_inf := (1 - exp(-(1 - w_municipality) * r_inf_denmark - w_municipality * r_inf_municipality)) ] + # TODO replace indices c(1, 3, 6) with names tmp2 <- inf_pressure_total[, c(1, 3, 6)] names(tmp2)[3] <- "prob_inf_new" - # Evaluate probability of infection + # merge probability of infection to the individual + # NB for future use - should maybe include a check on reasonable values [0;1] ibm[tmp2, on = c("municipality_id", "age_groups"), prob_inf := prob_inf_new] - # NB for future use - should maybe include a check on reasonable values [0;1] ibm[, prob_inf := prob_inf * vac_fac * rel_risk_parish * lockdown_fac] # Randomly infect some individuals based on probability @@ -598,10 +605,12 @@ sim_list <- foreach(run_this = (first_run - 1 + 1:n_runs), .packages = "data.tab # Implement the effect of vaccination for (vac_id in 1:n_vac) { + # this is binary vaccination effect, works for alpha variant, from delta should be different ibm[vac_type == vac_id & vac_time == v_vac_tt_effect[vac_id], - `:=`(vacF_ec = delta_vac_effect[vac_id], - prob_hospital = delta_red_hospital_risk * prob_hospital, - vac_fac_trans = red_transmission_vac)] + ':='(vac_fac = sample( c(1,0) , .N , prob = + c(1-v_vac_effect[vac_id],v_vac_effect[vac_id]),replace=TRUE), + prob_hospital = red_hospital_risk * prob_hospital, + vac_fac_trans = red_transmission_vac)] } From 5582eb0b96b88c15b13fba839d29dbc86d03f117 Mon Sep 17 00:00:00 2001 From: kaare-gr Date: Fri, 19 Jan 2024 14:54:18 +0100 Subject: [PATCH 13/13] docs(popIBM): add context comment --- manuscripts/popIBM/popIBM.R | 25 +++++++++++++++++-------- 1 file changed, 17 insertions(+), 8 deletions(-) diff --git a/manuscripts/popIBM/popIBM.R b/manuscripts/popIBM/popIBM.R index 5a188c0f..51f31142 100644 --- a/manuscripts/popIBM/popIBM.R +++ b/manuscripts/popIBM/popIBM.R @@ -342,35 +342,40 @@ sim_list <- foreach(run_this = (first_run - 1 + 1:n_runs), .packages = "data.tab # Count population by municipality, age group and vaccinated / un-vaccinated tmp <- ibm[, .(pop = .N), keyby = .(municipality_id, age_groups, !(vac_time < breaks_vac[[1]] | is.na(vac_time)))] + # ensure that there are no missing groups with zero + naming tmp <- tmp[t_pop_age_vac, , on = c("age_groups", "vac_time")] names(tmp)[c(3, 5)] <- c("vac_status", "t_pop") names(n_test_age_vac)[1] <- "age_groups" + + # round numbers of test n_test_age_vac$age_groups <- as.integer(n_test_age_vac$age_groups) + # merge number of tests onto subpopulations tmp_test <- tmp[n_test_age_vac, , on = c("age_groups", "vac_status")] + # calculate test probability tmp_test[, p_test_corr := w_test / t_pop] } - - - + # ensure that there are no missing groups with zero tmp2 <- data.table(municipality_id = u_municipality_ids, p_test_fac = p_test_corr) - tmp <- tmp_test[tmp2, , on = "municipality_id"] + # adjust test probability and ensure number of test match expected tmp[, p_test_corr := p_test_corr * p_test_fac] tmp[, p_test_corr := p_test_corr * sum(n_test_age_vac$w_test) / sum(p_test_corr * pop)] tmp[, vac_eff_dose := as.integer(vac_status)] + # tests must be distributed to both 1 and 2 dose vaccinated tmp2 <- copy(tmp[vac_eff_dose == 1L, ]) tmp2[, vac_eff_dose := 2L] tmp <- rbindlist(list(tmp, tmp2)) + # merge test probability onto individuals ibm[tmp, on = c("municipality_id", "age_groups", "vac_eff_dose"), p_test := p_test_corr] } @@ -526,8 +531,8 @@ sim_list <- foreach(run_this = (first_run - 1 + 1:n_runs), .packages = "data.tab # an individual is affected by the biggest lockdown either in parish or municipality population[, lockdown_max := pmax(parish_fac, municipality_fac)] - # Weighted sum as lockdown factor - population[, lockdown_fac := 1 * (1 - lockdown_max) + lockdown_factor * lockdown_max] + # local lockdown reduces activity, but proportionally more when national restriction are low + population[, lockdown_fac := 1 - lockdown_max * (1 - lockdown_factor) ] # Merging on ibm: ibm[population, on = c("parish_id", "municipality_id"), lockdown_fac := lockdown_fac] @@ -537,7 +542,7 @@ sim_list <- foreach(run_this = (first_run - 1 + 1:n_runs), .packages = "data.tab # Infected individuals with different variants for (variant_id in 1:n_variants) { - # Calculate the number of infected persons per municipality + age_group + # Calculate the number of effectively infectious persons per municipality + age_group inf_persons_municipality <- ibm[ disease == 2L & variant == variant_id, .(inf_persons = sum(lockdown_fac * non_isolated * vac_fac_trans)), @@ -545,6 +550,8 @@ sim_list <- foreach(run_this = (first_run - 1 + 1:n_runs), .packages = "data.tab ][all_pop_combi, on = c("municipality_id", "age_groups")] inf_persons_municipality[is.na(inf_persons), inf_persons := 0] + + # accounting for variant transmissibility inf_persons_municipality[, inf_persons := inf_persons * v_rel_beta[variant_id]] # Calculate the infection pressure by multiplying betas onto inf_persons @@ -555,9 +562,11 @@ sim_list <- foreach(run_this = (first_run - 1 + 1:n_runs), .packages = "data.tab names(inf_pressure)[2] <- "r_inf_municipality" inf_pressure <- dt_pop_municipality[inf_pressure, , on = "municipality_id"] + + # convert to true rate by adjusting for population size inf_pressure[, r_inf_municipality := r_inf_municipality / pop] - # calculate infected persons nationwide + # calculate effectively infectious persons nationwide inf_persons <- inf_persons_municipality[, .(inf_persons = sum(inf_persons)), by = .(age_groups)] # calculate nationwide infection pressure