diff --git a/.gitignore b/.gitignore index 59546af7..b6c14be4 100644 --- a/.gitignore +++ b/.gitignore @@ -2,6 +2,7 @@ private* venv* *.exe +*.log **/.idea/ **/.DS_Store /wd/* diff --git a/src/analysis/multi-source_statistical/1_mcmc.R b/src/analysis/multi-source_statistical/1_mcmc.R index 51e26a8a..906d9d9f 100644 --- a/src/analysis/multi-source_statistical/1_mcmc.R +++ b/src/analysis/multi-source_statistical/1_mcmc.R @@ -28,8 +28,11 @@ dir.create(out_dir, showWarnings = F, recursive = T) #---- configure model data ----# +# country code (iso-2) +country <- "UA" + # define model name -model_name <- "4_props_covs_model" +model_name <- "5_agesex_model" # source model-specific config functions source(file.path(src_dir, "models", paste0(model_name, "_config.R"))) @@ -38,40 +41,52 @@ source(file.path(src_dir, "models", paste0(model_name, "_config.R"))) dir.create(file.path(out_dir, "modelling", model_name, "mcmc"), recursive = T, showWarnings = F) # load data -if (file.exists(file.path(in_dir, "cod-ps_2023", "DO_NOT_SHARE_UKR_ADM2_POP_2023.csv"))) { - codps_N1 <- read.csv(file.path(in_dir, "cod-ps_2023", "DO_NOT_SHARE_UKR_ADM2_POP_2023.csv")) +if (file.exists(file.path(in_dir, "COD-PS", "2023", "DO_NOT_SHARE_UKR_ADM2_POP_2023.csv"))) { + codps_N1 <- read.csv(file.path(in_dir, "COD-PS", "2023", "DO_NOT_SHARE_UKR_ADM2_POP_2023.csv")) +} else { + codps_N1 <- read.csv(file.path(in_dir, "COD-PS", "2023", "UKR_ADM2_POP_2023_sim.csv")) +} + +if (file.exists(file.path(in_dir, "COD-PS", "2024", "DO_NOT_SHARE_UKR_ADM2_POP_2024_Sept_27.csv"))) { + codps_N2 <- read.csv(file.path(in_dir, "COD-PS", "2024", "DO_NOT_SHARE_UKR_ADM2_POP_2024_Sept_27.csv")) } else { - codps_N1 <- read.csv(file.path(in_dir, "cod-ps_2023", "UKR_ADM2_POP_2023_sim.csv")) + codps_N2 <- read.csv(file.path(in_dir, "COD-PS", "2024", "UKR_ADM2_POP_2024_sim.csv")) } # create model data md <- model_data( - idx = read.csv(file.path(out_dir, "ua_master_index.csv")), - idx_F = read.csv(file.path(out_dir, "population_proxy", "social_media_audience", "ua_facebook_audience.csv")), - idx_G = read.csv(file.path(out_dir, "population_proxy", "social_media_audience", "ua_instagram_audience.csv")), + idx = read.csv(file.path(out_dir, paste0(tolower(country), "_master_index", ".csv"))), + idx_F = read.csv(file.path(out_dir, "population_proxy", "social_media_audience", paste0(tolower(country), "_facebook_audience", ".csv"))), + idx_G = read.csv(file.path(out_dir, "population_proxy", "social_media_audience", paste0(tolower(country), "_instagram_audience", ".csv"))), covs = read.csv(file.path(out_dir, "covariates", "final", "ua_covariates_oblast.csv")), + outside_border = read.csv(file.path(out_dir, "population_proxy", "crossing_borders", "dat_refugees.csv")), codps = read.csv(file.path(data_dir, "cod-ps", "population_baseline.csv")), codps_N1 = codps_N1, date_N1 = "2023-07-01", confidence_N1 = 0.1, # 95% chance true pop is within confidence_N1*100 percent of codps_N1 estimate - outside_border = read.csv(file.path(out_dir, "population_proxy", "crossing_borders", "dat_refugees.csv")), - last_date = "2023-08-31", # max(idx$t_name) + codps_N2 = codps_N2, + date_N2 = "2024-05-14", + confidence_N2 = 0.1, # 95% chance true pop is within confidence_N1*100 percent of codps_N1 estimate + last_date = "2024-05-14", # max(idx$t_name) process_drop_locations = c(), # 3782=Donetska, 3791=Luhanksa, 3788=Crimea, 3797=Sevastopol observation_drop_locations = c(3782, 3788, 3791, 3797), process_cov_select = read.csv(file.path(in_dir, "ua_r_covariates_oblast_selection.csv")), observation_cov_select = read.csv(file.path(in_dir, "ua_p_covariates_oblast_selection.csv")) ) -# idx = read.csv(file.path(out_dir, "ua_master_index.csv")) -# idx_F = read.csv(file.path(out_dir, "population_proxy", "social_media_audience", "ua_facebook_audience.csv")) -# idx_G = read.csv(file.path(out_dir, "population_proxy", "social_media_audience", "ua_instagram_audience.csv")) +# idx = read.csv(file.path(out_dir, paste0(tolower(country), "_master_index", ".csv"))) +# idx_F = read.csv(file.path(out_dir, "population_proxy", "social_media_audience", paste0(tolower(country), "_facebook_audience", ".csv"))) +# idx_G = read.csv(file.path(out_dir, "population_proxy", "social_media_audience", paste0(tolower(country), "_instagram_audience", ".csv"))) # covs = read.csv(file.path(out_dir, "covariates", "final", "ua_covariates_oblast.csv")) +# outside_border = read.csv(file.path(out_dir, "population_proxy", "crossing_borders", "dat_refugees.csv")) # codps = read.csv(file.path(data_dir, "cod-ps", "population_baseline.csv")) # codps_N1 = codps_N1 # date_N1 = "2023-07-01" -# confidence_N1 = 0.1 # 95% chance true pop is within confidence_N1*100 percent of codps_N1 estimate -# outside_border = read.csv(file.path(out_dir, "population_proxy", "crossing_borders", "dat_refugees.csv")) -# last_date = "2023-08-31" # max(idx$t_name) +# confidence_N1 = 0.1 # 95% chance true pop is within confidence_N1*100 percent of codps_N1 estimate +# codps_N2 = codps_N2 +# date_N2 = "2024-05-14" +# confidence_N2 = 0.1 # 95% chance true pop is within confidence_N1*100 percent of codps_N1 estimate +# last_date = "2024-05-14" # max(idx$t_name) # process_drop_locations = c() # 3782=Donetska, 3791=Luhanksa, 3788=Crimea, 3797=Sevastopol # observation_drop_locations = c(3782, 3788, 3791, 3797) # process_cov_select = read.csv(file.path(in_dir, "ua_r_covariates_oblast_selection.csv")) diff --git a/src/analysis/multi-source_statistical/1_mcmc.sh b/src/analysis/multi-source_statistical/1_mcmc.sh new file mode 100755 index 00000000..47f7d168 --- /dev/null +++ b/src/analysis/multi-source_statistical/1_mcmc.sh @@ -0,0 +1,13 @@ +#!/bin/bash + +# Define the path to your R script +R_SCRIPT="1_mcmc.R" + +# execute in the background with logging +nohup Rscript "$R_SCRIPT" > mcmc.log 2>&1 & + +# Get the Process ID (PID) of the process just started +PID=$! + +echo "R script '$R_SCRIPT' is running in the background with PID: $PID" +echo "You can check progress by running: tail -f mcmc.log" \ No newline at end of file diff --git a/src/analysis/multi-source_statistical/models/3_covs_model_config.R b/src/analysis/multi-source_statistical/models/3_covs_model_config.R index 71b13aca..4ca6e2bb 100644 --- a/src/analysis/multi-source_statistical/models/3_covs_model_config.R +++ b/src/analysis/multi-source_statistical/models/3_covs_model_config.R @@ -14,6 +14,7 @@ model_data <- function( # set seed for random number generators seed <- round(runif(1, 1, 1e6)) set.seed(seed) + md$seed <- seed #---- location and date filtering ----# diff --git a/src/analysis/multi-source_statistical/models/5_agesex_model b/src/analysis/multi-source_statistical/models/5_agesex_model new file mode 100755 index 00000000..7538a77b Binary files /dev/null and b/src/analysis/multi-source_statistical/models/5_agesex_model differ diff --git a/src/analysis/multi-source_statistical/models/5_agesex_model.stan b/src/analysis/multi-source_statistical/models/5_agesex_model.stan new file mode 100644 index 00000000..f88faeb8 --- /dev/null +++ b/src/analysis/multi-source_statistical/models/5_agesex_model.stan @@ -0,0 +1,265 @@ +data { + int T; + int I; + int A; + int S; + int C_full; + int C_adult; + int K_r; + int K_p; + + vector[T] y_N_tot; + matrix[T * I * (A-1) * S, K_r] X_r; + matrix[T * I * (A-1) * S, K_p] X_p; + + int n_F; + int n_G; + array[n_F] real y_F; + array[n_G] real y_G; + array[n_F] int tias_F; + array[n_G] int tias_G; + + array[T * I * (A-1) * S] int tt; + array[T * I * (A-1) * S] int ii; + array[T * I * (A-1) * S] int aa; + array[T * I * (A-1) * S] int ss; + + vector[T * I * S] ratio_children_female; + + vector[C_adult] A0; + vector[C_adult] A1; + vector[C_adult] A2; + vector[I * S] Y0; + vector[I * S] Y1; + vector[I * S] Y2; + + array[T * (I * S)] int young_idx; + array[T * C_adult] int adult_idx; + + array[T, C_full] int slice_full; + array[T, C_adult] int slice_adult; + array[T, C_adult] int slice_adult_lag; + array[T, I * S] int slice_young; + + int N_fem_a; + array[T * I, 3] int idx_female_triple; +} + +parameters { + vector[T * C_adult] log_r_raw; + real sigma_r; + real alpha_r; + vector[K_r] beta_r; + + + real sigma_delta_r; + real sigma_gamma_r_i; + real sigma_gamma_r_a; + real sigma_gamma_r_s; + + vector[T-1] delta_r_innov_raw; + vector[I] gamma_r_i_raw; + vector[A-2] gamma_r_a_raw; + vector[S] gamma_r_s_raw; + + + real alpha_p; + vector[K_p] beta_p; + real phi_p; + + real sigma_delta_p; + real sigma_gamma_p_i; + real sigma_gamma_p_a; + real sigma_gamma_p_s; + + vector[T-1] delta_p_innov_raw; + vector[I] gamma_p_i_raw; + vector[A-2] gamma_p_a_raw; + vector[S] gamma_p_s_raw; + + real sigma_F; + real sigma_G; + + real sigma_phi_p_a; + + vector[A-2] phi_p_a_raw; +} + +transformed parameters { + vector[I] gamma_p_i = sigma_gamma_p_i * gamma_p_i_raw; + vector[S] gamma_p_s = sigma_gamma_p_s * gamma_p_s_raw; + + vector[A-2] gamma_p_a_free = sigma_gamma_p_a * gamma_p_a_raw; + vector[A-1] gamma_p_a; + + gamma_p_a[1] = gamma_p_a_free[1]; + gamma_p_a[6] = gamma_p_a_free[1]; + + gamma_p_a[2] = gamma_p_a_free[2]; + gamma_p_a[3] = gamma_p_a_free[3]; + gamma_p_a[4] = gamma_p_a_free[4]; + gamma_p_a[5] = gamma_p_a_free[5]; + + gamma_p_a -= mean(gamma_p_a); + gamma_p_i -= mean(gamma_p_i); + gamma_p_s -= mean(gamma_p_s); + + + vector[T] delta_p; + delta_p[1] = 0; + for (t in 2:T) + delta_p[t] = delta_p[t-1] + sigma_delta_p * delta_p_innov_raw[t-1]; + + + vector[T * I * (A-1) * S] nu_p = alpha_p + + delta_p[tt] + + gamma_p_i[ii] + + gamma_p_a[aa] + + gamma_p_s[ss] + + X_p * beta_p; + + vector[A-2] phi_p_a_free = sigma_phi_p_a * phi_p_a_raw; + vector[A-1] phi_p_a; + + phi_p_a[1] = phi_p_a_free[1]; + phi_p_a[6] = phi_p_a_free[1]; + + phi_p_a[2] = phi_p_a_free[2]; + phi_p_a[3] = phi_p_a_free[3]; + phi_p_a[4] = phi_p_a_free[4]; + phi_p_a[5] = phi_p_a_free[5]; + + phi_p_a -= mean(phi_p_a); + + vector[T * I * (A-1) * S] p_F = inv_logit(nu_p); + vector[T * I * (A-1) * S] p_G = inv_logit(nu_p + phi_p + phi_p_a[aa]); + + +vector[T] delta_r; +delta_r[1] = 0; +for (t in 2:T) + delta_r[t] = delta_r[t-1] + sigma_delta_r * delta_r_innov_raw[t-1]; + + + vector[I] gamma_r_i = sigma_gamma_r_i * gamma_r_i_raw; + vector[S] gamma_r_s = sigma_gamma_r_s * gamma_r_s_raw; + + vector[A-2] gamma_r_a_free = sigma_gamma_r_a * gamma_r_a_raw; + vector[A-1] gamma_r_a; + + gamma_r_a[1] = gamma_r_a_free[1]; + gamma_r_a[6] = gamma_r_a_free[1]; + + gamma_r_a[2] = gamma_r_a_free[2]; + gamma_r_a[3] = gamma_r_a_free[3]; + gamma_r_a[4] = gamma_r_a_free[4]; + gamma_r_a[5] = gamma_r_a_free[5]; + + gamma_r_a -= mean(gamma_r_a); + gamma_r_i -= mean(gamma_r_i); + gamma_r_s -= mean(gamma_r_s); + + + vector[T * I * (A-1) * S] nu_r_log = alpha_r + + delta_r[tt] + + gamma_r_i[ii] + + gamma_r_a[aa] + + gamma_r_s[ss] + + X_r * beta_r; + + vector[T * I * (A-1) * S] N_adult; + vector[T * I * S] N_young; + + N_adult[slice_adult[1]] = A0; + + vector[T * C_adult] log_r = nu_r_log + sigma_r * log_r_raw; + vector[T * C_adult] r = exp(log_r); + + for (t in 2:T) { + N_adult[slice_adult[t]] = N_adult[slice_adult_lag[t]] .* r[slice_adult[t]]; + } + + vector[T * I] N_female_repof = rep_vector(0, T * I); + for (a in 1:3) + N_female_repof += N_adult[idx_female_triple[, a]]; + + + vector[T * I * S] N_female_repof_sex = + to_vector( rep_matrix(N_female_repof', S) ); + + N_young = N_female_repof_sex .* ratio_children_female; + N_young[slice_young[1]] = Y0; + + vector[T * C_full] N_full; + N_full[young_idx] = N_young; + N_full[adult_idx] = N_adult; + + vector[T] N_tot; + for (t in 1:T) + N_tot[t] = sum(N_full[slice_full[t]]); +} + +model { + alpha_r ~ normal(0, 0.2); + beta_r ~ normal(0, 0.5); + + alpha_p ~ normal(0, 2); + beta_p ~ normal(0, 0.5); + phi_p ~ normal(0, 0.3); + + + sigma_delta_r ~ normal(0, 0.01); + sigma_gamma_r_i ~ normal(0, 0.05); + sigma_gamma_r_a ~ normal(0, 0.05); + sigma_gamma_r_s ~ normal(0, 0.05); + + delta_r_innov_raw ~ std_normal(); + gamma_r_i_raw ~ std_normal(); + gamma_r_a_raw ~ std_normal(); + gamma_r_s_raw ~ std_normal(); + + sigma_delta_p ~ normal(0, 0.01); + sigma_gamma_p_i ~ normal(0, 0.05); + sigma_gamma_p_a ~ normal(0, 0.05); + sigma_gamma_p_s ~ normal(0, 0.05); + + delta_p_innov_raw ~ std_normal(); + gamma_p_i_raw ~ std_normal(); + gamma_p_a_raw ~ std_normal(); + gamma_p_s_raw ~ std_normal(); + + sigma_phi_p_a ~ normal(0, 0.05); + phi_p_a_raw ~ std_normal(); + + sigma_F ~ exponential(5); + sigma_G ~ exponential(5); + + log_r_raw ~ std_normal(); + sigma_r ~ normal(0, 0.01); + + A1 ~ lognormal(log(N_adult[slice_adult[71]]), 0.02); + A2 ~ lognormal(log(N_adult[slice_adult[117]]), 0.02); + + Y1 ~ lognormal(log(fmax(N_young[slice_young[71]], 1e-6)), 0.02); + Y2 ~ lognormal(log(fmax(N_young[slice_young[117]], 1e-6)), 0.02); + + + for (t in 2:T) + if (t != 71 && t != 117) + N_tot[t] ~ lognormal(log(y_N_tot[t]), 0.05); + + + y_F ~ lognormal(log(N_adult[tias_F] .* p_F[tias_F]), sigma_F); + y_G ~ lognormal(log(N_adult[tias_G] .* p_G[tias_G]), sigma_G); +} + +generated quantities { + array[n_F] real F_hat; + array[n_G] real G_hat; + + for (i in 1:n_F) + F_hat[i] = lognormal_rng(log(N_adult[tias_F[i]] .* p_F[tias_F[i]]), sigma_F); + + for (j in 1:n_G) + G_hat[j] = lognormal_rng(log(N_adult[tias_G[j]] .* p_G[tias_G[j]]), sigma_G); +} diff --git a/src/analysis/multi-source_statistical/models/5_agesex_model_config.R b/src/analysis/multi-source_statistical/models/5_agesex_model_config.R new file mode 100644 index 00000000..a0a7337f --- /dev/null +++ b/src/analysis/multi-source_statistical/models/5_agesex_model_config.R @@ -0,0 +1,860 @@ +library(dplyr) +library(lubridate) +library(tidyr) + +#---- function to generate model data ----# +model_data <- function( + idx, idx_F, idx_G, + covs, codps, outside_border, + codps_N1, confidence_N1, date_N1, + codps_N2, confidence_N2, date_N2, + last_date, + process_drop_locations, observation_drop_locations, + process_cov_select, observation_cov_select) { + # model data + md <- list() + + # set seed for random number generators + seed <- round(runif(1, 1, 1e6)) + set.seed(seed) + md$seed <- seed + + # drop age groups + drop_ages <- c(7) + + # drop locations + combined_drop_locations <- unique(c(process_drop_locations, observation_drop_locations)) + + # revise master index + selected_locations <- idx %>% filter(!i_key %in% combined_drop_locations) %>% distinct(i) %>% pull() + selected_locations_key <- idx %>% filter(!i_key %in% combined_drop_locations) %>% distinct(i_key) %>% pull() + + i_idx <- idx %>% + select(i_key) %>% + distinct() %>% + filter(i_key %in% selected_locations_key) %>% + arrange(i_key) %>% + mutate(i = row_number()) + + md$idx <- idx %>% + filter(i %in% selected_locations, !a %in% drop_ages) %>% + select( + t, t_key, t_name, i_key, i_name, + a, a_key, a_name, s, s_key, s_name + ) %>% + left_join(i_idx, by = "i_key") %>% + arrange(t, i, a, s) %>% + mutate(tias = row_number()) %>% + group_by(a, s) %>% + mutate(ti = row_number()) %>% + ungroup() %>% + group_by(t, i) %>% + mutate(as = row_number()) %>% + ungroup() %>% + filter(t_name <= floor_date(as.Date(last_date), "week", week_start = 1)) %>% + select( + ti, as, tias, t, i, a, s, + t_key, t_name, i_key, i_name, + a_key, a_name, s_key, s_name + ) + + md$A <- length(unique(md$idx$a)) + length(drop_ages) # 7 + md$S <- length(unique(md$idx$s)) # 2 + md$I <- length(unique(md$idx$i)) # 23 + md$T <- length(unique(md$idx$t)) # 117 + md$C_full <- md$I * md$A * md$S + md$C_adult <- md$I * (md$A - 1) * md$S + md$C_young <- md$I * md$S + + adult_map <- expand.grid( + t = seq_len(md$T), + i = seq_len(md$I), + a = 1:(md$A - 1), + s = seq_len(md$S) + ) %>% + arrange(t, i, a, s) %>% + mutate(tias_adult = row_number()) %>% + select(t, i, a, s, tias_adult) + + idx_keys <- idx %>% + select(i, i_key, a, a_key, s, s_key) %>% + distinct() + + md$idx_F <- idx_F %>% + filter(i %in% selected_locations, a < 6) %>% + select(t, i, a, s, m, value) %>% + left_join(idx_keys, by = c("i", "a", "s")) %>% + select(-i) %>% + left_join(i_idx, by = "i_key") %>% + right_join(md$idx %>% select(t, i, a, s, ti, as, tias), by = c("t", "i", "a", "s")) %>% + left_join(adult_map, by = c("t", "i", "a", "s")) %>% + filter(!i_key %in% combined_drop_locations, value > 0, is.finite(value)) %>% + group_by(ti, as, tias, tias_adult, t, i, a, s) %>% + summarise(value = mean(value, na.rm = TRUE), .groups = "drop") %>% + select(ti, as, tias_full = tias, tias_adult, t, i, a, s, value) %>% + arrange(t, i, a, s) + + md$idx_G <- idx_G %>% + filter(i %in% selected_locations, a < 6) %>% + select(t, i, a, s, m, value) %>% + left_join(idx_keys, by = c("i", "a", "s")) %>% + select(-i) %>% + left_join(i_idx, by = "i_key") %>% + right_join(md$idx %>% select(t, i, a, s, ti, as, tias), by = c("t", "i", "a", "s")) %>% + left_join(adult_map, by = c("t", "i", "a", "s")) %>% + filter(!i_key %in% combined_drop_locations, value > 0, is.finite(value)) %>% + group_by(ti, as, tias, tias_adult, t, i, a, s) %>% + summarise(value = mean(value, na.rm = TRUE), .groups = "drop") %>% + select(ti, as, tias_full = tias, tias_adult, t, i, a, s, value) %>% + arrange(t, i, a, s) + + md$y_F <- md$idx_F$value + md$n_F <- length(md$y_F) + md$tias_F <- md$idx_F$tias_adult + + md$y_G <- md$idx_G$value + md$n_G <- length(md$y_G) + md$tias_G <- md$idx_G$tias_adult + + cols_to_pivot <- c( + "F_00_04", "F_05_09", "F_10_14", "F_15_19", "F_20_24", "F_25_29", "F_30_34", "F_35_39", + "F_40_44", "F_45_49", "F_50_54", "F_55_59", "F_60_64", "F_65_69", "F_70_74", "F_75_79", + "F_80Plus", + "M_00_04", "M_05_09", "M_10_14", "M_15_19", "M_20_24", "M_25_29", "M_30_34", "M_35_39", + "M_40_44", "M_45_49", "M_50_54", "M_55_59", "M_60_64", "M_65_69", "M_70_74", "M_75_79", + "M_80Plus" + ) + + date_start <- "2022-02-25" + + time_index_expanded <- tibble( + collection_date = seq(as.Date(date_start), as.Date(last_date), by = 1) + ) |> + mutate( + t_name = floor_date(as.Date(collection_date), "week", week_start = 1), + t_key = stringr::str_replace_all(as.character(t_name), "-", "") |> as.integer(), + t = t_name |> as.character() |> haven::as_factor() |> as.integer() + ) + + time_index <- time_index_expanded |> + distinct(t_name, t_key, t) |> + arrange(t) + + ############### + # COD-PS 2022 + ############### + md$N010 <- codps %>% + select(ADM1_PCODE, fb_key, all_of(cols_to_pivot)) %>% + pivot_longer( + cols = all_of(cols_to_pivot), + names_to = "demog", + values_to = "pop0" + ) |> + separate(demog, into = c("sex", "age_min", "age_max"), sep = "_") |> + filter(sex != "T", age_min != "TL") |> + rename(pcode = ADM1_PCODE, i_key = fb_key) |> + filter(i_key %in% selected_locations_key) %>% + mutate( + age_group10y = case_when( + age_min == "00" & age_max == "04" ~ "00-14", + age_min == "05" & age_max == "09" ~ "00-14", + age_min == "10" & age_max == "14" ~ "00-14", + age_min == "15" & age_max == "19" ~ "15-19", + age_min == "20" & age_max == "24" ~ "20-29", + age_min == "25" & age_max == "29" ~ "20-29", + age_min == "30" & age_max == "34" ~ "30-39", + age_min == "35" & age_max == "39" ~ "30-39", + age_min == "40" & age_max == "44" ~ "40-49", + age_min == "45" & age_max == "49" ~ "40-49", + age_min == "50" & age_max == "54" ~ "50-59", + age_min == "55" & age_max == "59" ~ "50-59", + age_min == "60" & age_max == "64" ~ "60-999", + age_min == "65" & age_max == "69" ~ "60-999", + age_min == "70" & age_max == "74" ~ "60-999", + age_min == "75" & age_max == "79" ~ "60-999", + age_min == "80Plus" ~ "60-999" + ) + ) |> + group_by(i_key, pcode, age_group10y, sex) |> + summarise(value = sum(pop0), .groups = "drop") |> + mutate( + a = case_when( + age_group10y == "00-14" ~ 7L, + age_group10y == "15-19" ~ 6L, + age_group10y == "20-29" ~ 1L, + age_group10y == "30-39" ~ 2L, + age_group10y == "40-49" ~ 3L, + age_group10y == "50-59" ~ 4L, + age_group10y == "60-999" ~ 5L + ), + s = if_else(sex == "F", 2L, 1L) + ) |> + left_join(i_idx, by = "i_key") |> + arrange(i, a, s) |> + mutate(ias = row_number()) %>% + group_by(i, s) %>% + mutate(is = row_number()) %>% + ungroup() %>% + filter(!is.na(a)) %>% + select(ias, is, i, a, s, value) + + md$N0 <- md$N010$value + md$ias_N0 <- md$N010$ias + + md$Y0 <- md$N010 |> + arrange(i, a, s) |> + filter(a == 7) %>% + pull(value) + + md$A0 <- md$N010 |> + arrange(i, a, s) |> + filter(a != 7) %>% + pull(value) + + md$N01 <- codps |> + filter(fb_key %in% selected_locations_key) %>% + arrange(match(fb_key, i_idx$i_key)) |> + pull(T_TL) + + full_pop22 <- codps %>% + select(fb_key, all_of(cols_to_pivot)) %>% + pivot_longer( + cols = all_of(cols_to_pivot), + names_to = c("sex", "age_min", "age_max"), + names_sep = "_", + values_to = "pop0" + ) %>% + filter(sex != "T", age_min != "TL") + + total_pop_all22 <- sum(full_pop22$pop0, na.rm = TRUE) + total_pop_sel22 <- full_pop22 %>% + filter(fb_key %in% selected_locations_key) %>% + summarise(sum_sel = sum(pop0, na.rm = TRUE)) %>% + pull(sum_sel) + + prop_sel22 <- total_pop_sel22 / total_pop_all22 + prop_sel22 + + ############## + # COD-PS 2023 + ############## + md$N110 <- codps_N1 %>% + select(ADM1_PCODE, ADM2_PCODE, all_of(cols_to_pivot)) %>% + pivot_longer( + cols = all_of(cols_to_pivot), + names_to = "demog", + values_to = "pop1" + ) |> + group_by(ADM1_PCODE, demog) |> + summarise(pop1 = sum(pop1, na.rm = TRUE), .groups = "drop") |> + left_join(codps[, c("ADM1_PCODE", "fb_key")], by = "ADM1_PCODE") %>% + separate(demog, into = c("sex", "age_min", "age_max"), sep = "_") |> + filter(sex != "T", age_min != "TL") |> + rename(pcode = ADM1_PCODE, i_key = fb_key) |> + filter(i_key %in% selected_locations_key) %>% + mutate( + age_group10y = case_when( + age_min == "00" & age_max == "04" ~ "00-14", + age_min == "05" & age_max == "09" ~ "00-14", + age_min == "10" & age_max == "14" ~ "00-14", + age_min == "15" & age_max == "19" ~ "15-19", + age_min == "20" & age_max == "24" ~ "20-29", + age_min == "25" & age_max == "29" ~ "20-29", + age_min == "30" & age_max == "34" ~ "30-39", + age_min == "35" & age_max == "39" ~ "30-39", + age_min == "40" & age_max == "44" ~ "40-49", + age_min == "45" & age_max == "49" ~ "40-49", + age_min == "50" & age_max == "54" ~ "50-59", + age_min == "55" & age_max == "59" ~ "50-59", + age_min == "60" & age_max == "64" ~ "60-999", + age_min == "65" & age_max == "69" ~ "60-999", + age_min == "70" & age_max == "74" ~ "60-999", + age_min == "75" & age_max == "79" ~ "60-999", + age_min == "80Plus" ~ "60-999" + ) + ) %>% + group_by(i_key, pcode, age_group10y, sex) |> + summarise(value = sum(pop1), .groups = "drop") |> + mutate( + a = case_when( + age_group10y == "00-14" ~ 7L, + age_group10y == "15-19" ~ 6L, + age_group10y == "20-29" ~ 1L, + age_group10y == "30-39" ~ 2L, + age_group10y == "40-49" ~ 3L, + age_group10y == "50-59" ~ 4L, + age_group10y == "60-999" ~ 5L + ), + s = if_else(sex == "F", 2L, 1L) + ) |> + left_join(i_idx, by = "i_key") |> + arrange(i, a, s) |> + mutate(ias = row_number()) %>% + group_by(i, s) %>% + mutate(is = row_number()) %>% + ungroup() %>% + filter(!is.na(a)) %>% + select(ias, is, i, a, s, value) + + md$N1 <- md$N110$value + md$ias_N1 <- md$N110$ias + + md$Y1 <- md$N110 |> + arrange(i, a, s) |> + filter(a == 7) %>% + pull(value) + + md$A1 <- md$N110 |> + arrange(i, a, s) |> + filter(a != 7) %>% + pull(value) + + md$N11 <- codps_N1 %>% + select(ADM1_PCODE, ADM2_PCODE, all_of(cols_to_pivot)) %>% + pivot_longer( + cols = all_of(cols_to_pivot), + names_to = "demog", + values_to = "pop1" + ) |> + group_by(ADM1_PCODE) |> + summarise(pop1 = sum(pop1, na.rm = TRUE), .groups = "drop") |> + left_join(codps[, c("ADM1_PCODE", "fb_key")], by = "ADM1_PCODE") %>% + filter(fb_key %in% selected_locations_key) %>% + arrange(match(fb_key, i_idx$i_key)) %>% + pull(pop1) + + pop1_week <- floor_date(as.Date(date_N1), "week", week_start = 1) + t_pop1 <- time_index %>% + filter(t_name == pop1_week) %>% + pull(t) + t_pop1 + + full_pop23 <- codps_N1 %>% + select(ADM1_PCODE, ADM2_PCODE, all_of(cols_to_pivot)) %>% + pivot_longer( + cols = all_of(cols_to_pivot), + names_to = "demog", + values_to = "pop1" + ) %>% + group_by(ADM1_PCODE, demog) %>% + summarise(pop1 = sum(pop1, na.rm = TRUE), .groups = "drop") %>% + separate(demog, into = c("sex", "age_min", "age_max"), sep = "_") %>% + filter(sex != "T", age_min != "TL") + + total_pop_all23 <- sum(full_pop23$pop1, na.rm = TRUE) + + full_pop23_sel <- full_pop23 %>% + left_join(codps %>% select(ADM1_PCODE, fb_key), by = "ADM1_PCODE") %>% + filter(fb_key %in% selected_locations_key) + + total_pop_sel23 <- sum(full_pop23_sel$pop1, na.rm = TRUE) + prop_sel23 <- total_pop_sel23 / total_pop_all23 + prop_sel23 + + ################ + # COD-PS 2024 + ################ + md$N210 <- codps_N2 %>% + select(ADM1_PCODE, ADM2_PCODE, all_of(cols_to_pivot)) %>% + pivot_longer( + cols = all_of(cols_to_pivot), + names_to = "demog", + values_to = "pop2" + ) %>% + group_by(ADM1_PCODE, demog) %>% + summarise(pop2 = sum(pop2, na.rm = TRUE), .groups = "drop") %>% + left_join(codps[, c("ADM1_PCODE", "fb_key")], by = "ADM1_PCODE") %>% + separate(demog, into = c("sex", "age_min", "age_max"), sep = "_") %>% + filter(sex != "T", age_min != "TL") %>% + rename(pcode = ADM1_PCODE, i_key = fb_key) %>% + filter(i_key %in% selected_locations_key) %>% + mutate( + age_group10y = case_when( + age_min == "00" & age_max == "04" ~ "00-14", + age_min == "05" & age_max == "09" ~ "00-14", + age_min == "10" & age_max == "14" ~ "00-14", + age_min == "15" & age_max == "19" ~ "15-19", + age_min == "20" & age_max == "24" ~ "20-29", + age_min == "25" & age_max == "29" ~ "20-29", + age_min == "30" & age_max == "34" ~ "30-39", + age_min == "35" & age_max == "39" ~ "30-39", + age_min == "40" & age_max == "44" ~ "40-49", + age_min == "45" & age_max == "49" ~ "40-49", + age_min == "50" & age_max == "54" ~ "50-59", + age_min == "55" & age_max == "59" ~ "50-59", + age_min == "60" & age_max == "64" ~ "60-999", + age_min == "65" & age_max == "69" ~ "60-999", + age_min == "70" & age_max == "74" ~ "60-999", + age_min == "75" & age_max == "79" ~ "60-999", + age_min == "80Plus" ~ "60-999" + ) + ) %>% + group_by(i_key, pcode, age_group10y, sex) %>% + summarise(value = sum(pop2), .groups = "drop") %>% + mutate( + a = case_when( + age_group10y == "00-14" ~ 7L, + age_group10y == "15-19" ~ 6L, + age_group10y == "20-29" ~ 1L, + age_group10y == "30-39" ~ 2L, + age_group10y == "40-49" ~ 3L, + age_group10y == "50-59" ~ 4L, + age_group10y == "60-999" ~ 5L + ), + s = if_else(sex == "F", 2L, 1L) + ) %>% + left_join(i_idx, by = "i_key") %>% + arrange(i, a, s) %>% + mutate(ias = row_number()) %>% + group_by(i, s) %>% + mutate(is = row_number()) %>% + ungroup() %>% + filter(!is.na(a)) %>% + select(ias, is, i, a, s, value) + + md$N2 <- md$N210$value + md$ias_N2 <- md$N210$ias + + md$Y2 <- md$N210 |> + arrange(i, a, s) |> + filter(a == 7) %>% + pull(value) + + md$A2 <- md$N210 |> + arrange(i, a, s) |> + filter(a != 7) %>% + pull(value) + + md$N21 <- codps_N2 %>% + select(ADM1_PCODE, ADM2_PCODE, all_of(cols_to_pivot)) %>% + pivot_longer( + cols = all_of(cols_to_pivot), + names_to = "demog", + values_to = "pop2" + ) |> + group_by(ADM1_PCODE) |> + summarise(pop2 = sum(pop2, na.rm = TRUE), .groups = "drop") |> + left_join(codps[, c("ADM1_PCODE", "fb_key")], by = "ADM1_PCODE") %>% + arrange(match(fb_key, i_idx$i_key)) %>% + pull(pop2) + + + pop2_week <- floor_date(as.Date(date_N2), "week", week_start = 1) + t_pop2 <- time_index %>% + filter(t_name == pop2_week) %>% + pull(t) + t_pop2 + + full_pop24 <- codps_N2 %>% + select(ADM1_PCODE, ADM2_PCODE, all_of(cols_to_pivot)) %>% + pivot_longer( + cols = all_of(cols_to_pivot), + names_to = "demog", + values_to = "pop2" + ) %>% + group_by(ADM1_PCODE, demog) %>% + summarise(pop2 = sum(pop2, na.rm = TRUE), .groups = "drop") %>% + separate(demog, into = c("sex", "age_min", "age_max"), sep = "_") %>% + filter(sex != "T", age_min != "TL") + + total_pop_all24 <- sum(full_pop24$pop2, na.rm = TRUE) + + full_pop24_sel <- full_pop24 %>% + left_join(codps %>% select(ADM1_PCODE, fb_key), by = "ADM1_PCODE") %>% + filter(fb_key %in% selected_locations_key) + + total_pop_sel24 <- sum(full_pop24_sel$pop2, na.rm = TRUE) + prop_sel24 <- total_pop_sel24 / total_pop_all24 + prop_sel24 + + ###################################### + # Total population at every time step + ###################################### + weekly_avg <- outside_border |> + mutate(week = floor_date(as.Date(date), "week", week_start = 1)) |> + group_by(week) |> + summarise(avg_value = mean(individuals, na.rm = TRUE), .groups = "drop") |> + filter( + week >= min(md$idx$t_name), + week <= max(md$idx$t_name) + ) + + border_df <- time_index %>% + select(t, t_name) %>% + left_join(weekly_avg %>% rename(t_name = week), by = "t_name") %>% + mutate(avg_value = if_else(t == 1, 0, avg_value)) + + md$y_N_tot <- border_df %>% + mutate( + baseline = case_when( + t < t_pop1 ~ sum(md$N01), + t == t_pop1 ~ sum(md$N11), + t < t_pop2 ~ sum(md$N11), + t == t_pop2 ~ sum(md$N21), + TRUE ~ sum(md$N21) + ), + prop = case_when( + t <= t_pop1 ~ prop_sel22, + t_pop1 < t & t <= t_pop2 ~ prop_sel23, + TRUE ~ prop_sel24 + ), + ytmp = baseline - (avg_value * prop), + y_N_tot = as.integer(ytmp) + ) %>% + pull(y_N_tot) + + rm(weekly_avg) + + children <- md$N010 %>% + filter(a == 7) %>% + mutate(t = 1) %>% + rename(children = value) %>% + rbind( + md$N110 %>% filter(a == 7) %>% mutate(t = t_pop1) %>% rename(children = value) + ) %>% + rbind( + md$N210 %>% filter(a == 7) %>% mutate(t = t_pop2) %>% rename(children = value) + ) + + female_repo <- md$N010 %>% + filter(a >= 1 & a <= 3, s == 2) %>% + group_by(i) %>% + summarise(female_repo = sum(value, na.rm = TRUE), .groups = "drop") %>% + mutate(t = 1) %>% + rbind( + md$N110 %>% + filter(a >= 1 & a <= 3, s == 2) %>% + group_by(i) %>% + summarise(female_repo = sum(value, na.rm = TRUE), .groups = "drop") %>% + mutate(t = t_pop1) + ) %>% + rbind( + md$N210 %>% + filter(a >= 1 & a <= 3, s == 2) %>% + group_by(i) %>% + summarise(female_repo = sum(value, na.rm = TRUE), .groups = "drop") %>% + mutate(t = t_pop2) + ) + + md$ratio_children_female <- expand.grid( + t = seq_len(md$T), + i = seq_len(md$I), + s = seq_len(md$S) + ) %>% + left_join(children) %>% + left_join(female_repo) %>% + mutate(ratio_children_female = children / female_repo) %>% + group_by(i, s) %>% + arrange(t, .by_group = TRUE) %>% + mutate(ratio_children_female_im = zoo::na.approx(ratio_children_female, x = t, na.rm = FALSE)) %>% + ungroup() %>% + arrange(t, i, s) %>% + pull(ratio_children_female_im) + + table(!is.na(md$ratio_children_female)) + + t <- unique(md$idx$t) + i <- unique(md$idx$i) + a <- seq_len(md$A - 1) + s <- unique(md$idx$s) + + I <- md$I + A <- md$A + S <- md$S + T <- md$T + + C_full <- I * A * S + C_adult <- I * (A - 1) * S + C_young <- I * S + IS <- I * S + + adult_grid <- expand.grid(t = 1:T, i = 1:I, a = 1:(A - 1), s = 1:S) |> + dplyr::arrange(t, i, a, s) + + md$tt <- adult_grid$t + md$ii <- adult_grid$i + md$aa <- adult_grid$a + md$ss <- adult_grid$s + + md$adult_grid <- adult_grid[order(adult_grid$t, adult_grid$i, adult_grid$a, adult_grid$s), ] + md$female_core_mask <- as.integer(md$adult_grid$a %in% 1:7 & md$adult_grid$s == 2) + md$adult_grid$idx_tis <- (md$adult_grid$t - 1) * I * S + (md$adult_grid$i - 1) * S + md$adult_grid$s + md$tis_index <- md$adult_grid$idx_tis + + full_grid <- expand.grid(t = 1:T, i = 1:I, a = 1:A, s = 1:S) %>% + arrange(t, i, a, s) %>% + mutate(idx_full = row_number()) + + md$young_pos_full <- full_grid %>% + filter(t == 1, a == 7) %>% + arrange(i, a, s) %>% + pull(idx_full) %>% + as.integer() + + md$young_idx <- full_grid %>% + filter(a == 7) %>% + arrange(t, i, a, s) %>% + pull(idx_full) %>% + as.integer() + + md$adult_pos_full <- full_grid %>% + filter(t == 1, a != 7) %>% + arrange(i, a, s) %>% + pull(idx_full) %>% + as.integer() + + md$adult_idx <- full_grid %>% + filter(a != 7) %>% + arrange(t, i, a, s) %>% + pull(idx_full) %>% + as.integer() + + md$N_y <- T * length(md$young_pos_full) + md$N_a <- T * length(md$adult_pos_full) + + md$slice_full <- matrix(seq_len(T * C_full), nrow = T, byrow = TRUE) + md$slice_adult <- matrix(seq_len(T * C_adult), nrow = T, byrow = TRUE) + md$slice_adult_lag <- rbind(rep(1L, C_adult), md$slice_adult[-T, ]) + + md$slice_young <- matrix(seq_len(T * C_young), nrow = T, byrow = TRUE) + md$slice_young_lag <- rbind(rep(1L, C_young), md$slice_young[-T, ]) + + f_core <- (md$ss == 2L) & (md$aa %in% 1:3) + md$idx_female_a234 <- which(f_core) + md$idx_female_triple <- matrix(md$idx_female_a234, ncol = 3, byrow = TRUE) + md$N_fem_a <- length(md$idx_female_a234) + + idx_a1519 <- adult_map %>% + filter(a == 6) %>% + arrange(t, i, s) %>% + pull(tias_adult) + idx_a2029 <- adult_map %>% + filter(a == 1) %>% + arrange(t, i, s) %>% + pull(tias_adult) + + md$idx_a1519 <- as.integer(idx_a1519) + md$idx_a2029 <- as.integer(idx_a2029) + md$N_tie <- length(idx_a1519) + length(idx_a2029) + + standardise_cols <- function(df, vars) { + stats <- vector("list", length(vars)) + names(stats) <- vars + + for (v in vars) { + x <- df[[v]] + m <- mean(x, na.rm = TRUE) + s <- sd(x, na.rm = TRUE) + + if (!is.finite(s) || s <= 0) { + df[[v]] <- 0 + stats[[v]] <- c(mean = m, sd = 1) + } else { + df[[v]] <- (x - m) / s + stats[[v]] <- c(mean = m, sd = s) + } + } + + list(df = df, stats = stats) + } + + covs$time_std[is.na(covs$time_std)] <- "none" + covs$space_std[is.na(covs$space_std)] <- "none" + + rcov_select <- process_cov_select %>% + filter(select == 1) %>% + mutate(across(c(time_std, space_std), ~ replace(., is.na(.), "none"))) + + cov_names_r <- rcov_select %>% + transmute(nice_name = paste(covariate, sum_stat, time_std, space_std, sep = "_")) %>% + pull() + + idx <- expand.grid( + t = 1:md$T, + i = 1:md$I, + a = min(md$aa):max(md$aa), + s = 1:md$S + ) %>% + arrange(t, i, a, s) + + Xr_df <- idx + for (k in seq_along(cov_names_r)) { + sel <- rcov_select[k, ] + col_real <- cov_names_r[k] + + cov_k <- covs %>% + filter( + covariate == sel$covariate, + sum_stat == sel$sum_stat, + time_std == sel$time_std, + space_std == sel$space_std + ) %>% + select(t, i, value_std) %>% + rename(!!col_real := value_std) + + Xr_df <- left_join(Xr_df, cov_k, by = c("t", "i")) + } + + Xr_df[is.na(Xr_df)] <- 0 + + std_r <- standardise_cols(Xr_df, cov_names_r) + Xr_df <- std_r$df + md$X_r_means <- sapply(std_r$stats, `[[`, "mean") + md$X_r_sds <- sapply(std_r$stats, `[[`, "sd") + + md$X_r <- as.matrix(Xr_df["acled_withfatalities_raw_none_none"]) + md$K_r <- ncol(md$X_r) + + pcov_select <- observation_cov_select %>% + filter(select == 1) %>% + mutate(across(c(time_std, space_std), ~ replace(., is.na(.), "none"))) + + cov_names_p <- pcov_select %>% + transmute(nice_name = paste(covariate, sum_stat, time_std, space_std, sep = "_")) %>% + pull() + + Xp_df <- idx + for (k in seq_along(cov_names_p)) { + sel <- pcov_select[k, ] + col_real <- cov_names_p[k] + + cov_k <- covs %>% + filter( + covariate == sel$covariate, + sum_stat == sel$sum_stat, + time_std == sel$time_std, + space_std == sel$space_std + ) %>% + select(t, i, value_std) %>% + rename(!!col_real := value_std) + + Xp_df <- left_join(Xp_df, cov_k, by = c("t", "i")) + } + + Xp_df[is.na(Xp_df)] <- 0 + + std_p <- standardise_cols(Xp_df, cov_names_p) + Xp_df <- std_p$df + md$X_p_means <- sapply(std_p$stats, `[[`, "mean") + md$X_p_sds <- sapply(std_p$stats, `[[`, "sd") + + md$X_p <- as.matrix(Xp_df["pwtt_sum_24week_none_none"]) + md$K_p <- ncol(md$X_p) + + return(md) +} + +#---- function to generate initial values ----# +init_generator <- function(md, chain_id = 1) { + set.seed(2025 + chain_id) + + T <- md$T + I <- md$I + A <- md$A + S <- md$S + K_r <- md$K_r + K_p <- md$K_p + + if (is.null(md$C_adult)) { + stop("md$C_adult is NULL; please pass C_adult in md (= I*(A-1)*S).") + } + C_adult <- md$C_adult + N_cells <- T * C_adult + + A_free <- max(A - 2, 0) + Tm1 <- max(T - 1, 0) + + if (!is.null(md$tt) && length(md$tt) != N_cells) { + stop(sprintf("length(md$tt)=%d but T*C_adult=%d (T=%d, C_adult=%d).", + length(md$tt), N_cells, T, C_adult)) + } + if (!is.null(md$X_r) && nrow(md$X_r) != N_cells) { + stop(sprintf("nrow(md$X_r)=%d but T*C_adult=%d (T=%d, C_adult=%d).", + nrow(md$X_r), N_cells, T, C_adult)) + } + if (!is.null(md$X_p) && nrow(md$X_p) != N_cells) { + stop(sprintf("nrow(md$X_p)=%d but T*C_adult=%d (T=%d, C_adult=%d).", + nrow(md$X_p), N_cells, T, C_adult)) + } + + draw_sigma <- function(lo, hi) runif(1, lo, hi) + + beta_r <- if (K_r > 0) rep(0, K_r) else numeric(0) + beta_p <- if (K_p > 0) rep(0, K_p) else numeric(0) + + alpha_r <- 0 + alpha_p <- 0 + phi_p <- 0 + + sigma_F <- runif(1, 0.05, 0.20) + sigma_G <- runif(1, 0.05, 0.20) + + sigma_r <- runif(1, 0.005, 0.03) + + sigma_delta_r <- draw_sigma(0.005, 0.05) + sigma_gamma_r_i <- draw_sigma(0.005, 0.05) + sigma_gamma_r_a <- draw_sigma(0.005, 0.05) + sigma_gamma_r_s <- draw_sigma(0.005, 0.05) + + sigma_delta_p <- draw_sigma(0.005, 0.05) + sigma_gamma_p_i <- draw_sigma(0.005, 0.05) + sigma_gamma_p_a <- draw_sigma(0.005, 0.05) + sigma_gamma_p_s <- draw_sigma(0.005, 0.05) + + sigma_phi_p_a <- draw_sigma(0.005, 0.05) + + delta_r_innov_raw <- if (Tm1 > 0) rnorm(Tm1) else numeric(0) + delta_p_innov_raw <- if (Tm1 > 0) rnorm(Tm1) else numeric(0) + + gamma_r_i_raw <- if (I > 0) rnorm(I) else numeric(0) + gamma_p_i_raw <- if (I > 0) rnorm(I) else numeric(0) + + gamma_r_s_raw <- if (S > 0) rnorm(S) else numeric(0) + gamma_p_s_raw <- if (S > 0) rnorm(S) else numeric(0) + + gamma_r_a_raw <- if (A_free > 0) rnorm(A_free) else numeric(0) + gamma_p_a_raw <- if (A_free > 0) rnorm(A_free) else numeric(0) + + phi_p_a_raw <- if (A_free > 0) rnorm(A_free) else numeric(0) + + log_r_raw <- rnorm(N_cells, mean = 0, sd = 1) + log_r_raw <- pmin(pmax(log_r_raw, -3), 3) + + list( + log_r_raw = log_r_raw, + sigma_r = sigma_r, + alpha_r = alpha_r, + beta_r = beta_r, + + sigma_delta_r = sigma_delta_r, + sigma_gamma_r_i = sigma_gamma_r_i, + sigma_gamma_r_a = sigma_gamma_r_a, + sigma_gamma_r_s = sigma_gamma_r_s, + + delta_r_innov_raw = delta_r_innov_raw, + gamma_r_i_raw = gamma_r_i_raw, + gamma_r_a_raw = gamma_r_a_raw, + gamma_r_s_raw = gamma_r_s_raw, + + alpha_p = alpha_p, + beta_p = beta_p, + phi_p = phi_p, + + sigma_delta_p = sigma_delta_p, + sigma_gamma_p_i = sigma_gamma_p_i, + sigma_gamma_p_a = sigma_gamma_p_a, + sigma_gamma_p_s = sigma_gamma_p_s, + + delta_p_innov_raw = delta_p_innov_raw, + gamma_p_i_raw = gamma_p_i_raw, + gamma_p_a_raw = gamma_p_a_raw, + gamma_p_s_raw = gamma_p_s_raw, + + sigma_F = sigma_F, + sigma_G = sigma_G, + + sigma_phi_p_a = sigma_phi_p_a, + phi_p_a_raw = phi_p_a_raw + ) +} diff --git a/src/data_prep/pop_data/10_create_master_index.R b/src/data_prep/pop_data/10_create_master_index.R index a5ee94e2..8375ce5e 100644 --- a/src/data_prep/pop_data/10_create_master_index.R +++ b/src/data_prep/pop_data/10_create_master_index.R @@ -1,40 +1,16 @@ -# This script creates a master index for all unique combinations of time steps, regions, age, -# and sex groups. This index will be used in downstream modelling scripts. - -# cleanup rm(list = ls()) gc() -# Load required helpers source(file.path(here::here(), "src", "helpers", "R_helpers", "generic.R")) source(file.path(here::here(), "src", "helpers", "R_helpers", "data_querying.R")) -# Input parameters ------------------------------------------------------- - - -# agesex demographic groups - agesex <- c( - "F_20_29", "F_30_39", "F_40_49", "F_50_59", "F_60Plus", - "M_20_29", "M_30_39", "M_40_49", "M_50_59", "M_60Plus" + "F_20_29", "F_30_39", "F_40_49", "F_50_59", "F_60Plus", "F_15_19", + "M_20_29", "M_30_39", "M_40_49", "M_50_59", "M_60Plus", "M_15_19" ) -# agesex <- c("F_13Plus", "F_18Plus", "F_20Plus", "F_13_19", "F_15_49", "F_15_64", "F_18_34", "F_20_29", "F_30_39", "F_40_49", "F_50_59", "F_60Plus", "F_65Plus", -# "M_13Plus", "M_18Plus", "M_20Plus", "M_13_19", "M_15_49", "M_15_64", "M_18_34", "M_20_29", "M_30_39", "M_40_49", "M_50_59", "M_60Plus", "M_65Plus", -# "T_13Plus", "T_18Plus", "T_20Plus", "T_13_19", "T_15_49", "T_15_64", "T_18_34", "T_20_29", "T_30_39", "T_40_49", "T_50_59", "T_60Plus", "T_65Plus") - -# Function to create regular age groupings -# create_regular_group <- function(age_gap, gender = c("F", "M"), age_max=60) { -# age_group <- c(paste(seq(from = 20, to = age_max-age_gap, by = age_gap), -# seq(from = 20+age_gap-1, to = age_max, by = age_gap), sep = "_"), paste0(age_max, 'Plus')) -# lapply(gender, function(g) paste(g, age_group, sep = "_")) |> unlist() -# } -# agesex <- create_regular_group(10) - -# date date_start <- "2022-02-25" -date_end <- "2025-08-01" - +date_end <- "2024-05-14" country <- "UA" meta_keys <- read_csv(file.path(env$repo_dir, "data", "meta", paste0(tolower(country), "_meta_keys.csv"))) @@ -42,44 +18,61 @@ pcodes <- read_csv(file = file.path(env$repo_dir, "data", "cod-ps", "population_ output_label <- "" -# create master_index ---------------------------------------------------- - geo_index <- meta_keys |> rename(i_key = geo_key) |> distinct(i_key) |> left_join( pcodes |> - mutate(macroregion = case_when( - name == "Autonomous Republic of Crimea" ~ "Autonomous", - name == "Cherkasy Oblast" ~ "Center", - name == "Chernihiv Oblast" ~ "North", - name == "Chernivtsi Oblast" ~ "West", - name == "Dnipropetrovsk Oblast" ~ "East", - name == "Donetsk Oblast" ~ "East", - name == "Ivano-Frankivsk Oblast" ~ "West", - name == "Kharkiv Oblast" ~ "East", - name == "Kherson Oblast" ~ "South", - name == "Khmelnytskyi Oblast" ~ "West", - name == "Kyiv Oblast" ~ "North", - name == "Kirovohrad Oblast" ~ "Center", - name == "Kyiv" ~ "City", - name == "Luhansk Oblast" ~ "East", - name == "Lviv Oblast" ~ "West", - name == "Mykolaiv Oblast" ~ "South", - name == "Odessa Oblast" ~ "South", - name == "Poltava Oblast" ~ "Center", - name == "Rivne Oblast" ~ "West", - name == "Sevastopol" ~ "Autonomous", - name == "Sumy Oblast" ~ "North", - name == "Ternopil Oblast" ~ "West", - name == "Vinnytsia Oblast" ~ "Center", - name == "Volyn Oblast" ~ "West", - name == "Zakarpattia Oblast" ~ "West", - name == "Zaporizhia Oblast" ~ "East", - name == "Zhytomyr Oblast" ~ "North" - )) |> - select(fb_key, ADM1_PCODE, ADM1_EN, macroregion) |> - rename(i_key = fb_key, i_name = ADM1_EN), + mutate( + macroregion = case_when( + name == "Autonomous Republic of Crimea" ~ "Autonomous", + name == "Cherkasy Oblast" ~ "Center", + name == "Chernihiv Oblast" ~ "North", + name == "Chernivtsi Oblast" ~ "West", + name == "Dnipropetrovsk Oblast" ~ "East", + name == "Donetsk Oblast" ~ "East", + name == "Ivano-Frankivsk Oblast" ~ "West", + name == "Kharkiv Oblast" ~ "East", + name == "Kherson Oblast" ~ "South", + name == "Khmelnytskyi Oblast" ~ "West", + name == "Kyiv Oblast" ~ "North", + name == "Kirovohrad Oblast" ~ "Center", + name == "Kyiv" ~ "City", + name == "Luhansk Oblast" ~ "East", + name == "Lviv Oblast" ~ "West", + name == "Mykolaiv Oblast" ~ "South", + name == "Odessa Oblast" ~ "South", + name == "Poltava Oblast" ~ "Center", + name == "Rivne Oblast" ~ "West", + name == "Sevastopol" ~ "Autonomous", + name == "Sumy Oblast" ~ "North", + name == "Ternopil Oblast" ~ "West", + name == "Vinnytsia Oblast" ~ "Center", + name == "Volyn Oblast" ~ "West", + name == "Zakarpattia Oblast" ~ "West", + name == "Zaporizhia Oblast" ~ "East", + name == "Zhytomyr Oblast" ~ "North" + ), + e_name = case_when( + macroregion == "Autonomous" ~ "Autonomous & East", + macroregion == "East" ~ "Autonomous & East", + macroregion == "Center" ~ "Center", + macroregion == "City" ~ "City", + macroregion == "North" ~ "North", + macroregion == "South" ~ "South", + macroregion == "West" ~ "West" + ), + e = case_when( + e_name == "Autonomous & East" ~ 6, + e_name == "Center" ~ 2, + e_name == "City" ~ 1, + e_name == "North" ~ 5, + e_name == "South" ~ 4, + e_name == "West" ~ 3 + ) %>% as.integer() + ) |> + select(fb_key, ADM1_PCODE, ADM1_EN, macroregion, e_name, e) |> + rename(i_key = fb_key, i_name = ADM1_EN) ) |> arrange(i_key) |> mutate( @@ -92,8 +85,8 @@ time_index_expanded <- tibble( ) |> mutate( t_name = floor_date(as.Date(collection_date), "week", week_start = 1), - t_key = str_replace_all(as.character(t_name), "-", "") |> as.integer(), - t = t_name |> as.character() |> as_factor() |> as.integer() + t_key = str_replace_all(as.character(t_name), "-", "") |> as.integer(), + t = t_name |> as.character() |> as_factor() |> as.integer() ) time_index <- time_index_expanded |> @@ -108,13 +101,14 @@ agesex_index <- lapply( } ) |> bind_rows() |> - arrange(age_min, desc(agesex)) |> mutate( a_name = str_sub(agesex, 3), - a = as.integer(factor(paste0(age_min, age_max))), + a = as.integer(factor(a_name, + levels = c("20_29", "30_39", "40_49", "50_59", "60Plus", "15_19") + )), a_key = paste0(age_min, str_pad(age_max, 3, pad = "0")) |> as.integer(), s_name = str_sub(agesex, 1, 1), - s = ifelse(s_name == "F", 2, 1), # Assign 2 for F and 1 for M + s = ifelse(s_name == "F", 2, 1), s_key = s ) |> select(-age_min, -age_max, -gender) @@ -130,10 +124,9 @@ master_index <- expand_grid( left_join(geo_index, by = "i") |> left_join(agesex_index, by = c("a", "s")) |> mutate(parameter = row_number()) |> - arrange(t, i, a, s) - + arrange(t, i, a, s) |> + mutate(tias = row_number()) -# Write output ----------------------------------------------------------- write_csv(master_index, file.path(out_dir, paste0(tolower(country), "_master_index", output_label, ".csv"))) write_csv(time_index_expanded, file.path(out_dir, paste0(tolower(country), "_time_index", output_label, ".csv"))) diff --git a/src/data_prep/pop_data/20_prepare_sma_data.R b/src/data_prep/pop_data/20_prepare_sma_data.R index 57678c3c..606d983d 100644 --- a/src/data_prep/pop_data/20_prepare_sma_data.R +++ b/src/data_prep/pop_data/20_prepare_sma_data.R @@ -25,7 +25,6 @@ master_index <- read_csv(file.path(out_dir, paste0(tolower(country), "_master_in ) time_index <- read_csv(file.path(out_dir, paste0(tolower(country), "_time_index", output_label, ".csv"))) -# Retrieve parameters date_start <- min(master_index$t_name) date_end <- max(master_index$t_name) + 6 agesex <- unique(paste(master_index$s_name, master_index$a_name, sep = "_")) @@ -51,15 +50,14 @@ sma_facebook <- retrieve_sma_data("facebook") sma_instagram <- retrieve_sma_data("instagram") # Process social media audience data - master_index_without_t <- master_index |> rowwise() |> mutate( agesex = paste(s_name, a_name, sep = "_") ) -process_sma_data <- function(sma_data, metric = audience_metric) { - sma_data_ <- sma_data |> +process_sma_data <- function(sma_data, metric) { + sma_data_ <- sma_facebook |> mutate(t_name = floor_date(collection_date, "week", week_start = 1)) |> rename( value = all_of(metric), @@ -76,9 +74,12 @@ process_sma_data <- function(sma_data, metric = audience_metric) { return(sma_data_) } -sma_facebook_processed <- process_sma_data(sma_facebook) -sma_instagram_processed <- process_sma_data(sma_instagram) +sma_facebook_processed <- process_sma_data(sma_facebook, metric = audience_metric) +sma_instagram_processed <- process_sma_data(sma_instagram, metric = audience_metric) # Write processed data to CSV files +saveRDS(sma_facebook, file.path(out_dir_pop, paste0(tolower(country), "_facebook_audience_raw", output_label, ".rds"))) +saveRDS(sma_instagram, file.path(out_dir_pop, paste0(tolower(country), "_instagram_audience_raw", output_label, ".rds"))) + write_csv(sma_facebook_processed, file.path(out_dir_pop, paste0(tolower(country), "_facebook_audience", output_label, ".csv"))) write_csv(sma_instagram_processed, file.path(out_dir_pop, paste0(tolower(country), "_instagram_audience", output_label, ".csv"))) diff --git a/src/helpers/R_helpers/data_querying.R b/src/helpers/R_helpers/data_querying.R index cd44527d..a97c5cc5 100644 --- a/src/helpers/R_helpers/data_querying.R +++ b/src/helpers/R_helpers/data_querying.R @@ -39,7 +39,7 @@ query_sql <- function(con, sql, geom_col = NULL) { # Define query_api function -query_api <- function(endpoint, args, max_attempts = 3, url = "http://18.135.72.18/api/v1/", token = env$sma_API_token) { +query_api <- function(endpoint, args, max_attempts = 3, url = paste0("http://", Sys.getenv("smaDB_host"), "/api/v1/"), token = Sys.getenv("sma_API_token")) { args$token <- token # Submit query as GET request diff --git a/src/sandbox/multi-source_statistical/models_andrea/397_covs_model/1_mcmc_397.R b/src/sandbox/multi-source_statistical/models_andrea/397_covs_model/1_mcmc_397.R new file mode 100644 index 00000000..992ea863 --- /dev/null +++ b/src/sandbox/multi-source_statistical/models_andrea/397_covs_model/1_mcmc_397.R @@ -0,0 +1,46 @@ +library(cmdstanr) +library(posterior) +library(bayesplot) +library(here) +library(parallel) +library(dplyr) +library(lubridate) +library(reshape2) + +env <- new.env() +source(here::here(".env"), local = env) + +dir.create(file.path(here::here(), "UkraineNowPop/wd"), showWarnings = FALSE, recursive = TRUE) +setwd(file.path(here::here(), "UkraineNowPop/wd")) + + +src_dir <- file.path(env$repo_dir, "src") +out_dir <- file.path(env$out_dir) +dir.create(out_dir, showWarnings = FALSE, recursive = TRUE) + +country <- "UA" +model_name <- "397_covs_model" + +model_out_dir <- file.path(out_dir, "modelling", model_name, "mcmc") +dir.create(model_out_dir, recursive = TRUE, showWarnings = FALSE) + +chains <- 4 +warmup <- 500 +samples <- 1500 + +md <- readRDS(file.path(model_out_dir, paste0("md_", model_name, ".rds"))) + +mod <- cmdstan_model(file.path(src_dir, "analysis", "multi-source_statistical", "models", model_name, paste0(model_name, ".stan")), cpp_options = list(stan_threads = TRUE)) + +fit <- mod$sample( + data = md, + parallel_chains = chains, + iter_sampling = samples, + iter_warmup = warmup, + save_warmup = TRUE, + seed = md$seed, + threads_per_chain = 4 + ) + +fit$save_object(file = file.path(model_out_dir, paste0("fit_", model_name, ".rds"))) + diff --git a/src/sandbox/multi-source_statistical/models_andrea/397_covs_model/2_eval_397.R b/src/sandbox/multi-source_statistical/models_andrea/397_covs_model/2_eval_397.R new file mode 100644 index 00000000..adbc2686 --- /dev/null +++ b/src/sandbox/multi-source_statistical/models_andrea/397_covs_model/2_eval_397.R @@ -0,0 +1,876 @@ +library(bayesplot) +library(ggplot2) +library(cmdstanr) +library(posterior) +library(dplyr) +library(here) +library(patchwork) +library(purrr) + +env <- new.env() +source(here::here(".env"), local = env) + +ncores <- 10 + +dir.create(file.path(here::here(), "wd"), showWarnings = FALSE, recursive = TRUE) +setwd(file.path(here::here(), "wd")) + +# Directories +repo_dir <- env$repo_dir +src_dir <- file.path(repo_dir, "src", "data_prep", "pop_data") +in_dir <- env$in_dir +out_dir <- file.path(env$out_dir, "modelling") + +# Source functions +source(file.path(src_dir, "2_eval_fun.R")) + +country <- "UA" +model_name <- "397_covs_model" + + +fit <- readRDS(file.path(out_dir, model_name, "mcmc", paste0("fit_", model_name,".rds"))) +md <- readRDS(file.path(out_dir, model_name, "mcmc", paste0("md_", model_name, ".rds"))) + +dir.create(file.path(out_dir, model_name, "eval"), showWarnings = FALSE, recursive = TRUE) +dir.create(file.path(out_dir, model_name, "eval", "time_series_plots"), showWarnings = FALSE, recursive = TRUE) +dir.create(file.path(out_dir, model_name, "eval", "trace_plots"), showWarnings = FALSE, recursive = TRUE) +dir.create(file.path(out_dir, model_name, "eval", "population_pyramids"), showWarnings = FALSE, recursive = TRUE) + +out_dir_idx <- file.path("/data/home/andrea/git/OxfordDemSci/UkraineNowPop/wd/out") +idx = read.csv(file.path(out_dir_idx, paste0(tolower(country), "_master_index", ".csv"))) + + +A <- 7 +S <- 2 +I <- 23 +T <- 117 +C_full <- I*A*S +C_adult <- I*(A-1)*S +C_young <- I*S + + +#selected_locations <- seq(1, 27, by = 1)[-c(11, 5, 13, 19)] %>% as.integer() + +#selected_locations_key <- +# c(3800, 3801, 3781, 3804, 3802, 3803, 3783, 3790, 3787, +# 3792, 3793, 3794, 3795, 3796, 3798, 3799, 3784, 3785, +# 3786, 3778, 3780, 3779, 4290) %>% +# as.integer() + +#selected_locations_pcode <- c(#"UA01", "UA14", "UA44", "UA85", +#"UA71", "UA74", "UA73", "UA12", "UA26", "UA63", "UA65", "UA68", "UA35", +#"UA32", "UA46", "UA48", "UA51", "UA53", "UA56", "UA59", "UA61", "UA05", +#"UA07", "UA21", "UA23", "UA18", "UA80") +############################## +# Checking total population +############################## +names_N <- paste0("N_full[", 1:(md$T * md$C_full), "]") + +draws_df <- fit$draws(variables = "N_full", format = "df") +head(names(draws_df)) + +N_draws <- as.matrix(draws_df[, grep("^N_full\\[", names(draws_df))]) +N_means <- colMeans(N_draws) + + +N_t_mean <- sapply(seq_len(md$T), function(t) { + idxs <- ((t - 1) * md$C_full + 1):(t * md$C_full) + sum(N_means[idxs]) +}) + +t_names <- md$idx %>% + dplyr::distinct(t, t_name) %>% + dplyr::arrange(t) %>% + dplyr::pull(t_name) + +df_tot <- data.frame(t_name = t_names, + observed = md$y_N_tot, + predicted = N_t_mean) + +p <- ggplot(df_tot, aes(x = observed, y = predicted)) + + geom_point() + + geom_abline(slope = 1, intercept = 0, color = "red") + + labs(title = "Total Population Check", + x = "Observed y_N_tot", + y = "Posterior mean total N" + ) + theme_minimal() +p +ggsave(file.path(out_dir, model_name, "eval", "total_population_check.jpg"), + plot = p, width = 10, height = 8) +##################################################### +# Posterior predictive check +##################################################### +plot_postpred_fit( + dat = "y_F", + hat = "F_hat", + outfile = file.path(out_dir, model_name, "eval", paste0("postpredict_insamp_facebook", ".jpg"))) + + +dat <- "y_F" +hat <- "F_hat" + +draws <- fit$draws(hat, format = "df") |> select(!starts_with(".")) +prob_lower <- 0.025 +prob_upper <- 0.975 + +hat_mean <- apply(draws, 2, mean) +hat_lower <- apply(draws, 2, quantile, probs = prob_lower) +hat_upper <- apply(draws, 2, quantile, probs = prob_upper) + +df_idx_F <- as.data.frame(md$idx_F) %>% + mutate(parameter = row_number()) + +plot_df <- data.frame( + observed = md[[dat]], + predicted = hat_mean, + lower = hat_lower, + upper = hat_upper) %>% +tibble::rownames_to_column(var = "parameter") %>% +mutate(parameter = as.integer(gsub("[^0-9]", "", parameter))) %>% +left_join(df_idx_F, by = c("parameter")) %>% +left_join(idx %>% dplyr::select(t, t_name, i, i_name, macroregion, s, s_name, a, a_name)) %>% + mutate(year = case_when( + t_name>="2022-01-01" & t_name<= "2022-12-31" ~ 2022, + t_name>="2023-01-01" & t_name<= "2023-12-31" ~ 2023, + t_name>="2024-01-01" & t_name<= "2024-12-31" ~ 2024 + ) %>% as.integer()) + + +p <- ggplot(plot_df, aes(x = observed, y = predicted, colour = s_name)) + + geom_point() + + geom_errorbar(aes(ymin = lower, ymax = upper), width = 0) + + geom_abline(slope = 1, intercept = 0, colour = "red", linetype = "dashed") + + facet_wrap(year~macroregion, scales = "free") + + labs( + title = "Posterior Predictive Check by Macroregion", + x = paste0("Observed (", dat, ")"), + y = paste0("Predicted (", hat, ")"), + colour = "s_name" + ) + + theme_minimal() +p +ggsave(file.path(out_dir, model_name, "eval", "postpredict_insamp_facebook1.jpg"), plot = p, width = 10, height = 8) + + +p <- ggplot(plot_df, aes(x = observed, y = predicted, colour = s_name)) + + geom_point() + + geom_errorbar(aes(ymin = lower, ymax = upper), width = 0) + + geom_abline(slope = 1, intercept = 0, colour = "red", linetype = "dashed") + + facet_wrap(a_name~macroregion, scales = "free") + + labs( + title = "Posterior Predictive Check by Macroregion", + x = paste0("Observed (", dat, ")"), + y = paste0("Predicted (", hat, ")"), + colour = "a_name" + ) + + theme_minimal() +p +ggsave(file.path(out_dir, model_name, "eval", "postpredict_insamp_facebook2.jpg"), plot = p, width = 10, height = 8) + + + +plot_postpred_fit( + dat = "y_G", + hat = "G_hat", + outfile = file.path(out_dir, model_name, "eval", paste0("postpredict_insamp_instagram", ".jpg"))) + + +dat <- "y_G" +hat <- "G_hat" + +draws <- fit$draws(hat, format = "df") |> select(!starts_with(".")) +prob_lower <- 0.025 +prob_upper <- 0.975 + +hat_mean <- apply(draws, 2, mean) +hat_lower <- apply(draws, 2, quantile, probs = prob_lower) +hat_upper <- apply(draws, 2, quantile, probs = prob_upper) + +df_idx_G <- as.data.frame(md$idx_G) %>% + mutate(parameter = row_number()) + +plot_df <- data.frame( + observed = md[[dat]], + predicted = hat_mean, + lower = hat_lower, + upper = hat_upper) %>% +tibble::rownames_to_column(var = "parameter") %>% +mutate(parameter = as.integer(gsub("[^0-9]", "", parameter))) %>% +left_join(df_idx_G, by = c("parameter")) %>% +left_join(idx %>% dplyr::select(t, t_name, i, i_name, macroregion, s, s_name, a, a_name)) %>% + mutate(year = case_when( + t_name>="2022-01-01" & t_name<= "2022-12-31" ~ 2022, + t_name>="2023-01-01" & t_name<= "2023-12-31" ~ 2023, + t_name>="2024-01-01" & t_name<= "2024-12-31" ~ 2024 + ) %>% as.integer()) + + +p <- ggplot(plot_df, aes(x = observed, y = predicted, colour = s_name)) + + geom_point() + + geom_errorbar(aes(ymin = lower, ymax = upper), width = 0) + + geom_abline(slope = 1, intercept = 0, colour = "red", linetype = "dashed") + + facet_wrap(year~macroregion, scales = "free") + + labs( + title = "Posterior Predictive Check by Macroregion", + x = paste0("Observed (", dat, ")"), + y = paste0("Predicted (", hat, ")"), + colour = "s_name" + ) + + theme_minimal() +p +ggsave(file.path(out_dir, model_name, "eval", "postpredict_insamp_instagram1.jpg"), plot = p, width = 10, height = 8) + + +p <- ggplot(plot_df, aes(x = observed, y = predicted, colour = s_name)) + + geom_point() + + geom_errorbar(aes(ymin = lower, ymax = upper), width = 0) + + geom_abline(slope = 1, intercept = 0, colour = "red", linetype = "dashed") + + facet_wrap(a_name~macroregion, scales = "free") + + labs( + title = "Posterior Predictive Check by Macroregion", + x = paste0("Observed (", dat, ")"), + y = paste0("Predicted (", hat, ")"), + colour = "a_name" + ) + + theme_minimal() +p +ggsave(file.path(out_dir, model_name, "eval", "postpredict_insamp_instagram2.jpg"), plot = p, width = 10, height = 8) + +##################################################### +# Time series plots +##################################################### +options(scipen = 999) + +params <- c("N_adult", "N_young", "r", "p_F", "p_G") + +param_labels <- c( + "N" = "Children population (N_young)", + "r" = "Growth rate (r)", + "p_F" = "Facebook penetration rate (p_F)", + "p_G" = "Instagram penetration rate (p_G)" +) + +draws <- fit$draws(params, format = "df") |> select(!starts_with(".")) +prob_lower <- 0.025 +prob_upper <- 0.975 + +params_mean <- apply(draws, 2, mean) +params_lower <- apply(draws, 2, quantile, probs = prob_lower) +params_upper <- apply(draws, 2, quantile, probs = prob_upper) + + +i_idx <- idx %>% + dplyr::select(i_key) %>% + distinct() %>% + filter(i_key %in% selected_locations_key) %>% + arrange(i_key) %>% + mutate(i = row_number()) + +idx1 <- expand.grid(t = 1:T, i = 1:I, a = 1:A, s = 1:S) %>% + left_join(i_idx) %>% + arrange(t, i, a, s) %>% + mutate(tias = row_number()) %>% + mutate(tias_adult = NA_integer_, + tias_young = NA_integer_) + +idx1$tias_adult[idx1$a != 7L] <- seq_len(sum(idx1$a != 7L)) +idx1$tias_young[idx1$a == 7L] <- seq_len(sum(idx1$a == 7L)) + +plot_df_adult <- tibble( + parameter = names(params_mean), + predicted = params_mean, + lower = params_lower, + upper = params_upper +) |> + mutate(param = sub("\\[.*", "", parameter), + tias_adult = sub(".*\\[([^]]+)]$", "\\1", parameter) |> as.integer(), + tias_young = NA_real_) %>% +filter(param == "N_adult"| param =="r" | param =="p_F"|param =="p_G") %>% +left_join(idx1 %>% dplyr::select(tias_adult, tias, t, i_key, a, s)) %>% +#dplyr::select(-tias_adult) %>% +dplyr::select(parameter, predicted, lower, upper, param, tias, t, i_key, a, s, tias_young, tias_adult) + + +plot_df_young <- tibble( + parameter = names(params_mean), + predicted = params_mean, + lower = params_lower, + upper = params_upper +) |> + mutate(param = sub("\\[.*", "", parameter), + tias_young = sub(".*\\[([^]]+)]$", "\\1", parameter) |> as.integer(), + tias_adult = NA_real_) %>% +filter(param %in% c("N_young")) %>% +left_join(idx1 %>% dplyr::select(tias_young, tias, t, i_key, a, s)) %>% +#dplyr::select(-tias_young) %>% +dplyr::select(parameter, predicted, lower, upper, param, tias, t, i_key, a, s, tias_young, tias_adult) + + +plot_df <- plot_df_adult %>% rbind(plot_df_young) %>% + mutate(param = ifelse(param=="N_young"|param=="N_adult", "N", param)) %>% + left_join(idx %>% distinct(t, i, s, t_name, i_name, s_name, i_key)) %>% + mutate(a_name = case_when(a == 7 ~ "00-14", + a == 6 ~ "15-19", + a == 1 ~ "20-29", + a == 2 ~ "30-39", + a == 3 ~ "40-49", + a == 4 ~ "50-59", + a == 5 ~ "60+")) %>% + arrange(t,i,a,s) + + +plot_df <- plot_df %>% + filter(!(param %in% c("p_F", "p_G") & a == 6L)) + +plot_df <- plot_df %>% + mutate(t_name = factor(t_name, levels = sort(unique(t_name)))) + + +out_dir_plots <- file.path(out_dir, model_name, "eval", "time_series_plots") +dir.create(out_dir_plots, recursive = TRUE, showWarnings = FALSE) + +for (name in unique(plot_df$i_name)) { + + df_i <- filter(plot_df, i_name == name) + plots <- vector("list", length(param_vec)) + + for (k in seq_along(param_vec)) { + + p_code <- param_vec[k] + df_p <- filter(df_i, param == p_code) + + if (p_code %in% c("p_F", "p_G")) { + df_p <- df_p %>% + mutate(a_name = factor(a_name, levels = c("20-29","30-39","40-49","50-59","60+"))) + } + + labs <- unique(df_p$t_name) + breaks_vec <- if (length(labs) >= 5) labs[seq(1, length(labs), 5)] else labs + + p <- ggplot(df_p, aes(t_name, predicted, group = a_name)) + + geom_ribbon(aes(ymin = lower, ymax = upper, fill = a_name), + alpha = .20, colour = NA) + + geom_line (aes(colour = a_name)) + + geom_point (aes(colour = a_name, shape = a_name), size = 1.8) + + scale_x_discrete(breaks = breaks_vec) + + labs(x = "Date", + y = param_labels[[p_code]] %||% p_code, + colour = "Age group", fill = "Age group", shape = "Age group") + + facet_wrap(~s_name, ncol = 2) + + theme_minimal(base_size = 11) + + theme(axis.text.x = element_text(angle = 90, vjust = .5, hjust = 1), + strip.text = element_text(size = 12)) + + if (p_code == "N") p <- p + ggtitle(name) + + plots[[k]] <- p + } + + combined_plot <- wrap_plots(plots, ncol = 1) + + ggsave( + filename = file.path(out_dir_plots, + paste0(gsub("[^A-Za-z0-9_]", "_", name), ".jpg")), + plot = combined_plot, + width = 14, + height = 12 + ) +} +##################################################### +# Evolution of population pyramids +##################################################### +n_steps <- 30 +times_to_plot <- unique(round(seq(1, md$T, length.out = n_steps))) + +draws_N <- fit$draws(variables = c("N_adult", "N_young"), format = "df") |> select(!starts_with(".")) +prob_lower <- 0.025 +prob_upper <- 0.975 + +params_mean <- apply(draws_N, 2, mean) +params_lower <- apply(draws_N, 2, quantile, probs = prob_lower) +params_upper <- apply(draws_N, 2, quantile, probs = prob_upper) + + +i_idx <- idx %>% + dplyr::select(i_key) %>% + distinct() %>% + filter(i_key %in% selected_locations_key) %>% + arrange(i_key) %>% + mutate(i = row_number()) + +idx1 <- expand.grid(t = 1:T, i = 1:I, a = 1:A, s = 1:S) %>% + left_join(i_idx) %>% + arrange(t, i, a, s) %>% + mutate(tias = row_number()) %>% + mutate(tias_adult = NA_integer_, + tias_young = NA_integer_) + +idx1$tias_adult[idx1$a != 7L] <- seq_len(sum(idx1$a != 7L)) +idx1$tias_young[idx1$a == 7L] <- seq_len(sum(idx1$a == 7L)) + + +plot_df_adult <- tibble( + parameter = names(params_mean), + predicted = params_mean, + lower = params_lower, + upper = params_upper +) |> + mutate(param = sub("\\[.*", "", parameter), + tias_adult = sub(".*\\[([^]]+)]$", "\\1", parameter) |> as.integer()) %>% +filter(param %in% c("N_adult", "r", "p_F", "p_G")) %>% +left_join(idx1 %>% dplyr::select(tias_adult, tias, t, i_key, a, s)) %>% +dplyr::select(-tias_adult) %>% +dplyr::select(parameter, predicted, lower, upper, param, tias, t, i_key, a, s) + + +plot_df_young <- tibble( + parameter = names(params_mean), + predicted = params_mean, + lower = params_lower, + upper = params_upper +) |> + mutate(param = sub("\\[.*", "", parameter), + tias_young = sub(".*\\[([^]]+)]$", "\\1", parameter) |> as.integer()) %>% +filter(param %in% c("N_young")) %>% +left_join(idx1 %>% dplyr::select(tias_young, tias, t, i_key, a, s)) %>% +dplyr::select(-tias_young) %>% +dplyr::select(parameter, predicted, lower, upper, param, tias, t, i_key, a, s) + + +plot_df <- plot_df_adult %>% rbind(plot_df_young) %>% + mutate(param = ifelse(param=="N_young"|param=="N_adult", "N", param), + pop = ifelse(s==1, -predicted, predicted), + lower1 = ifelse(s==1, -lower, lower), + upper1 = ifelse(s==1, -upper, upper)) %>% + left_join(idx %>% distinct(t, i, s, t_name, i_name, s_name, i_key)) %>% + mutate(a_name = case_when(a == 6 ~ "00-14", + a == 7 ~ "15-19", + a == 1 ~ "20-29", + a == 2 ~ "30-39", + a == 3 ~ "40-49", + a == 4 ~ "50-59", + a == 5 ~ "60+")) %>% + filter(t %in% times_to_plot) %>% + arrange(t,i,a,s) + + +oblasts <- unique(plot_df$i_name) + +for (ob in oblasts) { + df_sub <- plot_df %>% filter(i_name == ob) + + p <- ggplot(df_sub, aes(x = a_name)) + + geom_col(aes(y = pop/1000, fill = s_name), width = 0.8) + + geom_errorbar( + aes(ymin = lower1/1000, ymax = upper1/1000, colour = s_name), + width = 0.3, + position = position_identity() + ) + + coord_flip() + + scale_y_continuous( + labels = abs, + name = "Population (in thousands)" + ) + + facet_wrap(~ t_name) + + labs( + x = "Age group", + fill = "Sex", + colour = "Sex", + title = paste0("Population pyramids – ", ob) + ) + + theme_minimal() + + theme( + panel.grid.major.y = element_blank(), + axis.title.y = element_blank(), + strip.text.x = element_text(size = 8), + panel.background = element_rect(fill = "white", colour = NA), + plot.background = element_rect(fill = "white", colour = NA) + ) + + safe_name <- gsub("[^A-Za-z0-9_]", "_", ob) + ggsave( + filename = paste0(safe_name, ".png"), + plot = p, + path = file.path(out_dir, model_name, "eval", "population_pyramids"), + width = 12, + height = 4 + ) +} + + +#################################################### +# Trace plots for location-time-specific parameters +##################################################### +pars_model <- c("N_adult", "N_young", "log_r", "p_F", "p_G") + +t_values_to_plot <- seq(1, md$T, by = 20) + +grid_info_all <- idx1 %>% + left_join( + idx %>% filter(i_key %in% selected_locations_key) %>% distinct(t, s, t_name, i_key, i_name, s_name), + by = c("t", "i_key", "s")) %>% + mutate( + a_name = case_when( + a == 7 ~ "00-14", + a == 6 ~ "15-19", + a == 1 ~ "20-29", + a == 2 ~ "30-39", + a == 3 ~ "40-49", + a == 4 ~ "50-59", + a == 5 ~ "60+", + TRUE ~ as.character(a) + ) + ) %>% + distinct(t, i, a, s, t_name, i_name, s_name, tias_adult, tias_young, a_name) %>% + filter(t %in% t_values_to_plot) %>% + arrange(t, i, a, s) + +set.seed(123) + +n_t_keep <- 3 +n_a_keep <- 3 +n_s_keep <- 2 + +t_keep <- sample(t_values_to_plot, size = min(n_t_keep, length(t_values_to_plot))) +a_keep <- sample(sort(unique(grid_info_all$a)), size = min(n_a_keep, dplyr::n_distinct(grid_info_all$a))) +s_keep <- sample(sort(unique(grid_info_all$s)), size = min(n_s_keep, dplyr::n_distinct(grid_info_all$s))) + +grid_info_all <- grid_info_all %>% + filter(t %in% t_keep, a %in% a_keep, s %in% s_keep) %>% + arrange(t, i, a, s) + + +include_warmup <- FALSE +n_warmup <- fit$metadata()$iter_warmup +if (is.null(n_warmup) || is.na(n_warmup)) n_warmup <- 0 + +for (row_id in seq_len(nrow(grid_info_all))) { + + grid_info <- slice(grid_info_all, row_id) + + loc_folder <- gsub("\\s+", "_", grid_info$i_name) + folder_path <- file.path(out_dir, model_name, "eval", "trace_plots", loc_folder) + dir.create(folder_path, recursive = TRUE, showWarnings = FALSE) + + idx_adult <- grid_info$tias_adult + idx_young <- grid_info$tias_young + + idx_for_par <- function(p) if (p == "N_young") idx_young else idx_adult + + pars_batch <- vapply( + pars_model, + function(p) sprintf("%s[%d]", p, idx_for_par(p)), + character(1) + ) + + pars_batch <- pars_batch[!grepl("\\[NA\\]", pars_batch)] + if (length(pars_batch) == 0) next + + draws <- fit$draws( + variables = pars_batch, + inc_warmup = include_warmup + ) + + trace_plot <- bayesplot::mcmc_trace( + draws, + n_warmup = if (include_warmup) n_warmup else 0, + facet_args = list( + ncol = 4, + labeller = label_parsed + ) + ) + + plot_label <- paste( + gsub("\\s+", "_", grid_info$t_name), + loc_folder, + paste0("a", grid_info$a, "_", grid_info$a_name), + paste0("s", grid_info$s, "_", gsub("\\s+", "_", grid_info$s_name)), + sep = "__" + ) + + outfile <- file.path(folder_path, paste0("traceplot_", plot_label, ".jpg")) + + ggsave( + filename = outfile, + plot = trace_plot, + width = 12, + height = max(3, length(pars_batch) * 2), + dpi = 300 + ) +} + + +#################################################### +# Other trace plots for location–time–specific parameters +#################################################### +plot_traces <- function(fit, params, prefix, + out_dir, model_name, + ncol = 5, + n_warmup = 300, + inc_warmup = FALSE, + chunk_size = 10, + width = 12, + height = 8) { + +base_path <- file.path(out_dir, model_name, "eval", "trace_plots", prefix) +dir.create(base_path, recursive = TRUE, showWarnings = FALSE) + +all_vars <- posterior::variables(fit$draws()) +params <- intersect(params, all_vars) + +if (length(params) == 0L) { +message(sprintf("[plot_traces] No matching parameters for prefix '%s' — skipping.", prefix)) +return(invisible(NULL)) +} + +chunks <- split(params, ceiling(seq_along(params) / chunk_size)) + +n_warmup_plot <- if (isTRUE(inc_warmup)) n_warmup else 0 + +for (i in seq_along(chunks)) { +pars_chunk <- chunks[[i]] +draws <- fit$draws(variables = pars_chunk, inc_warmup = inc_warmup) + +plt <- bayesplot::mcmc_trace( +draws, +pars = pars_chunk, +n_warmup = n_warmup_plot, +facet_args = list(ncol = ncol, labeller = ggplot2::label_parsed) +) + +ggplot2::ggtitle(sprintf("%s: chunk %d / %d", prefix, i, length(chunks))) + +ggplot2::ggsave( +filename = sprintf("%s_chunk%02d.jpg", prefix, i), +plot = plt, +path = base_path, +width = width, +height = height +) +} + +invisible(NULL) +} + +T <- md$T +I <- md$I +A <- md$A +S <- md$S +K_r <- md$K_r +K_p <- md$K_p + +A_adult <- A - 1 +A_free <- A - 2 + +C_adult <- if (!is.null(md$C_adult)) md$C_adult else I * A_adult * S +N_adult_time <- T * C_adult + + +pars_global <- c( +"alpha_r", +sprintf("beta_r[%d]", seq_len(K_r)), +"sigma_r", +"sigma_delta_r", +"sigma_gamma_r_i", "sigma_gamma_r_a", "sigma_gamma_r_s", + +"alpha_p", +sprintf("beta_p[%d]", seq_len(K_p)), +"phi_p", +"sigma_delta_p", +"sigma_gamma_p_i", "sigma_gamma_p_a", "sigma_gamma_p_s", + +"sigma_phi_p_a", +"sigma_F", "sigma_G" +) + +plot_traces( +fit, pars_global, +prefix = "01_global", +chunk_size = length(pars_global), +out_dir = out_dir, +model_name = model_name +) + +pars_rw_innov <- c( +sprintf("delta_r_innov_raw[%d]", seq_len(T - 1)), +sprintf("delta_p_innov_raw[%d]", seq_len(T - 1)) +) + +plot_traces( +fit, pars_rw_innov, +prefix = "02_rw_innov_raw", +chunk_size = 10, +out_dir = out_dir, +model_name = model_name +) + +every_t_steps <- unique(pmin(seq(2, T, by = 30), T)) + +plot_traces( +fit, sprintf("delta_r[%d]", every_t_steps), +prefix = "03_delta_r_sample", +ncol = 3, +chunk_size = length(every_t_steps), +out_dir = out_dir, +model_name = model_name +) + +plot_traces( +fit, sprintf("delta_p[%d]", every_t_steps), +prefix = "04_delta_p_sample", +ncol = 3, +chunk_size = length(every_t_steps), +out_dir = out_dir, +model_name = model_name +) + +pars_gamma_r_raw <- c( +sprintf("gamma_r_i_raw[%d]", seq_len(I)), +sprintf("gamma_r_a_raw[%d]", seq_len(A_free)), +sprintf("gamma_r_s_raw[%d]", seq_len(S)) +) + +plot_traces( +fit, pars_gamma_r_raw, +prefix = "05_gamma_r_raw", +out_dir = out_dir, +model_name = model_name +) + +pars_gamma_p_raw <- c( +sprintf("gamma_p_i_raw[%d]", seq_len(I)), +sprintf("gamma_p_a_raw[%d]", seq_len(A_free)), +sprintf("gamma_p_s_raw[%d]", seq_len(S)) +) + +plot_traces( +fit, pars_gamma_p_raw, +prefix = "06_gamma_p_raw", +out_dir = out_dir, +model_name = model_name +) + +pars_gamma_r_expanded <- c( +sprintf("gamma_r_i[%d]", seq_len(I)), +sprintf("gamma_r_a[%d]", seq_len(A_adult)), +sprintf("gamma_r_s[%d]", seq_len(S)) +) + +plot_traces( +fit, pars_gamma_r_expanded, +prefix = "07_gamma_r_expanded", +out_dir = out_dir, +model_name = model_name +) + +pars_gamma_p_expanded <- c( +sprintf("gamma_p_i[%d]", seq_len(I)), +sprintf("gamma_p_a[%d]", seq_len(A_adult)), +sprintf("gamma_p_s[%d]", seq_len(S)) +) + +plot_traces( +fit, pars_gamma_p_expanded, +prefix = "08_gamma_p_expanded", +out_dir = out_dir, +model_name = model_name +) + +plot_traces( +fit, sprintf("phi_p_a_raw[%d]", seq_len(A_free)), +prefix = "09_phi_p_a_raw", +chunk_size = 10, +out_dir = out_dir, +model_name = model_name +) + +plot_traces( +fit, sprintf("phi_p_a[%d]", seq_len(A_adult)), +prefix = "10_phi_p_a_expanded", +chunk_size = 10, +out_dir = out_dir, +model_name = model_name +) + +plot_traces( +fit, sprintf("N_tot[%d]", every_t_steps), +prefix = "11_N_tot_sample", +ncol = 3, +chunk_size = length(every_t_steps), +out_dir = out_dir, +model_name = model_name +) + +set.seed(2025) +idx50 <- sample.int(N_adult_time, 50) + +plot_traces( +fit, sprintf("log_r_raw[%d]", idx50), +prefix = "12_log_r_raw_sample", +chunk_size = 10, +out_dir = out_dir, +model_name = model_name +) + +plot_traces( +fit, sprintf("log_r[%d]", idx50), +prefix = "13_log_r_sample", +chunk_size = 10, +out_dir = out_dir, +model_name = model_name +) + +plot_traces( +fit, sprintf("r[%d]", idx50), +prefix = "14_r_sample", +chunk_size = 10, +out_dir = out_dir, +model_name = model_name +) + +plot_traces( +fit, sprintf("p_F[%d]", idx50), +prefix = "15_pF_sample", +chunk_size = 10, +out_dir = out_dir, +model_name = model_name +) + +plot_traces( +fit, sprintf("p_G[%d]", idx50), +prefix = "16_pG_sample", +chunk_size = 10, +out_dir = out_dir, +model_name = model_name +) + +plot_traces( +fit, sprintf("nu_r_log[%d]", idx50), +prefix = "17_nu_r_log_sample", +chunk_size = 10, +out_dir = out_dir, +model_name = model_name +) + +plot_traces( +fit, sprintf("nu_p[%d]", idx50), +prefix = "18_nu_p_sample", +chunk_size = 10, +out_dir = out_dir, +model_name = model_name +) + + + +##################################################### +## Summary statistics +##################################################### +fit_summary <- fit$summary(.cores = ncores) +write.csv(fit_summary, file.path(out_dir, model_name, "eval", "fit_summary.csv"), row.names = FALSE) + +not_converged <- fit_summary %>% + filter(!is.na(rhat), round(rhat, 1) > 1.1) %>% + mutate(var_base = sub("\\[.*$", "", variable)) %>% + arrange(desc(rhat)) + +write.csv(not_converged, file.path(out_dir, model_name, "eval", "bad_pars.csv"), row.names = FALSE) + diff --git a/src/sandbox/multi-source_statistical/models_andrea/397_covs_model/397_covs_model.stan b/src/sandbox/multi-source_statistical/models_andrea/397_covs_model/397_covs_model.stan new file mode 100644 index 00000000..f88faeb8 --- /dev/null +++ b/src/sandbox/multi-source_statistical/models_andrea/397_covs_model/397_covs_model.stan @@ -0,0 +1,265 @@ +data { + int T; + int I; + int A; + int S; + int C_full; + int C_adult; + int K_r; + int K_p; + + vector[T] y_N_tot; + matrix[T * I * (A-1) * S, K_r] X_r; + matrix[T * I * (A-1) * S, K_p] X_p; + + int n_F; + int n_G; + array[n_F] real y_F; + array[n_G] real y_G; + array[n_F] int tias_F; + array[n_G] int tias_G; + + array[T * I * (A-1) * S] int tt; + array[T * I * (A-1) * S] int ii; + array[T * I * (A-1) * S] int aa; + array[T * I * (A-1) * S] int ss; + + vector[T * I * S] ratio_children_female; + + vector[C_adult] A0; + vector[C_adult] A1; + vector[C_adult] A2; + vector[I * S] Y0; + vector[I * S] Y1; + vector[I * S] Y2; + + array[T * (I * S)] int young_idx; + array[T * C_adult] int adult_idx; + + array[T, C_full] int slice_full; + array[T, C_adult] int slice_adult; + array[T, C_adult] int slice_adult_lag; + array[T, I * S] int slice_young; + + int N_fem_a; + array[T * I, 3] int idx_female_triple; +} + +parameters { + vector[T * C_adult] log_r_raw; + real sigma_r; + real alpha_r; + vector[K_r] beta_r; + + + real sigma_delta_r; + real sigma_gamma_r_i; + real sigma_gamma_r_a; + real sigma_gamma_r_s; + + vector[T-1] delta_r_innov_raw; + vector[I] gamma_r_i_raw; + vector[A-2] gamma_r_a_raw; + vector[S] gamma_r_s_raw; + + + real alpha_p; + vector[K_p] beta_p; + real phi_p; + + real sigma_delta_p; + real sigma_gamma_p_i; + real sigma_gamma_p_a; + real sigma_gamma_p_s; + + vector[T-1] delta_p_innov_raw; + vector[I] gamma_p_i_raw; + vector[A-2] gamma_p_a_raw; + vector[S] gamma_p_s_raw; + + real sigma_F; + real sigma_G; + + real sigma_phi_p_a; + + vector[A-2] phi_p_a_raw; +} + +transformed parameters { + vector[I] gamma_p_i = sigma_gamma_p_i * gamma_p_i_raw; + vector[S] gamma_p_s = sigma_gamma_p_s * gamma_p_s_raw; + + vector[A-2] gamma_p_a_free = sigma_gamma_p_a * gamma_p_a_raw; + vector[A-1] gamma_p_a; + + gamma_p_a[1] = gamma_p_a_free[1]; + gamma_p_a[6] = gamma_p_a_free[1]; + + gamma_p_a[2] = gamma_p_a_free[2]; + gamma_p_a[3] = gamma_p_a_free[3]; + gamma_p_a[4] = gamma_p_a_free[4]; + gamma_p_a[5] = gamma_p_a_free[5]; + + gamma_p_a -= mean(gamma_p_a); + gamma_p_i -= mean(gamma_p_i); + gamma_p_s -= mean(gamma_p_s); + + + vector[T] delta_p; + delta_p[1] = 0; + for (t in 2:T) + delta_p[t] = delta_p[t-1] + sigma_delta_p * delta_p_innov_raw[t-1]; + + + vector[T * I * (A-1) * S] nu_p = alpha_p + + delta_p[tt] + + gamma_p_i[ii] + + gamma_p_a[aa] + + gamma_p_s[ss] + + X_p * beta_p; + + vector[A-2] phi_p_a_free = sigma_phi_p_a * phi_p_a_raw; + vector[A-1] phi_p_a; + + phi_p_a[1] = phi_p_a_free[1]; + phi_p_a[6] = phi_p_a_free[1]; + + phi_p_a[2] = phi_p_a_free[2]; + phi_p_a[3] = phi_p_a_free[3]; + phi_p_a[4] = phi_p_a_free[4]; + phi_p_a[5] = phi_p_a_free[5]; + + phi_p_a -= mean(phi_p_a); + + vector[T * I * (A-1) * S] p_F = inv_logit(nu_p); + vector[T * I * (A-1) * S] p_G = inv_logit(nu_p + phi_p + phi_p_a[aa]); + + +vector[T] delta_r; +delta_r[1] = 0; +for (t in 2:T) + delta_r[t] = delta_r[t-1] + sigma_delta_r * delta_r_innov_raw[t-1]; + + + vector[I] gamma_r_i = sigma_gamma_r_i * gamma_r_i_raw; + vector[S] gamma_r_s = sigma_gamma_r_s * gamma_r_s_raw; + + vector[A-2] gamma_r_a_free = sigma_gamma_r_a * gamma_r_a_raw; + vector[A-1] gamma_r_a; + + gamma_r_a[1] = gamma_r_a_free[1]; + gamma_r_a[6] = gamma_r_a_free[1]; + + gamma_r_a[2] = gamma_r_a_free[2]; + gamma_r_a[3] = gamma_r_a_free[3]; + gamma_r_a[4] = gamma_r_a_free[4]; + gamma_r_a[5] = gamma_r_a_free[5]; + + gamma_r_a -= mean(gamma_r_a); + gamma_r_i -= mean(gamma_r_i); + gamma_r_s -= mean(gamma_r_s); + + + vector[T * I * (A-1) * S] nu_r_log = alpha_r + + delta_r[tt] + + gamma_r_i[ii] + + gamma_r_a[aa] + + gamma_r_s[ss] + + X_r * beta_r; + + vector[T * I * (A-1) * S] N_adult; + vector[T * I * S] N_young; + + N_adult[slice_adult[1]] = A0; + + vector[T * C_adult] log_r = nu_r_log + sigma_r * log_r_raw; + vector[T * C_adult] r = exp(log_r); + + for (t in 2:T) { + N_adult[slice_adult[t]] = N_adult[slice_adult_lag[t]] .* r[slice_adult[t]]; + } + + vector[T * I] N_female_repof = rep_vector(0, T * I); + for (a in 1:3) + N_female_repof += N_adult[idx_female_triple[, a]]; + + + vector[T * I * S] N_female_repof_sex = + to_vector( rep_matrix(N_female_repof', S) ); + + N_young = N_female_repof_sex .* ratio_children_female; + N_young[slice_young[1]] = Y0; + + vector[T * C_full] N_full; + N_full[young_idx] = N_young; + N_full[adult_idx] = N_adult; + + vector[T] N_tot; + for (t in 1:T) + N_tot[t] = sum(N_full[slice_full[t]]); +} + +model { + alpha_r ~ normal(0, 0.2); + beta_r ~ normal(0, 0.5); + + alpha_p ~ normal(0, 2); + beta_p ~ normal(0, 0.5); + phi_p ~ normal(0, 0.3); + + + sigma_delta_r ~ normal(0, 0.01); + sigma_gamma_r_i ~ normal(0, 0.05); + sigma_gamma_r_a ~ normal(0, 0.05); + sigma_gamma_r_s ~ normal(0, 0.05); + + delta_r_innov_raw ~ std_normal(); + gamma_r_i_raw ~ std_normal(); + gamma_r_a_raw ~ std_normal(); + gamma_r_s_raw ~ std_normal(); + + sigma_delta_p ~ normal(0, 0.01); + sigma_gamma_p_i ~ normal(0, 0.05); + sigma_gamma_p_a ~ normal(0, 0.05); + sigma_gamma_p_s ~ normal(0, 0.05); + + delta_p_innov_raw ~ std_normal(); + gamma_p_i_raw ~ std_normal(); + gamma_p_a_raw ~ std_normal(); + gamma_p_s_raw ~ std_normal(); + + sigma_phi_p_a ~ normal(0, 0.05); + phi_p_a_raw ~ std_normal(); + + sigma_F ~ exponential(5); + sigma_G ~ exponential(5); + + log_r_raw ~ std_normal(); + sigma_r ~ normal(0, 0.01); + + A1 ~ lognormal(log(N_adult[slice_adult[71]]), 0.02); + A2 ~ lognormal(log(N_adult[slice_adult[117]]), 0.02); + + Y1 ~ lognormal(log(fmax(N_young[slice_young[71]], 1e-6)), 0.02); + Y2 ~ lognormal(log(fmax(N_young[slice_young[117]], 1e-6)), 0.02); + + + for (t in 2:T) + if (t != 71 && t != 117) + N_tot[t] ~ lognormal(log(y_N_tot[t]), 0.05); + + + y_F ~ lognormal(log(N_adult[tias_F] .* p_F[tias_F]), sigma_F); + y_G ~ lognormal(log(N_adult[tias_G] .* p_G[tias_G]), sigma_G); +} + +generated quantities { + array[n_F] real F_hat; + array[n_G] real G_hat; + + for (i in 1:n_F) + F_hat[i] = lognormal_rng(log(N_adult[tias_F[i]] .* p_F[tias_F[i]]), sigma_F); + + for (j in 1:n_G) + G_hat[j] = lognormal_rng(log(N_adult[tias_G[j]] .* p_G[tias_G[j]]), sigma_G); +} diff --git a/src/sandbox/multi-source_statistical/models_andrea/397_covs_model/397_covs_model_config.R b/src/sandbox/multi-source_statistical/models_andrea/397_covs_model/397_covs_model_config.R new file mode 100644 index 00000000..5b9d7243 --- /dev/null +++ b/src/sandbox/multi-source_statistical/models_andrea/397_covs_model/397_covs_model_config.R @@ -0,0 +1,768 @@ +library(dplyr) +library(lubridate) +library(tidyr) + +env <- new.env() + +source(here::here(".env"), local = env) + +wd_dir <- file.path(here::here(), "wd") +dir.create(wd_dir, showWarnings = FALSE, recursive = TRUE) +setwd(wd_dir) + +repo_dir <- env$repo_dir +data_dir <- file.path(repo_dir, "data") +in_dir <- env$in_dir +out_dir <- file.path(env$out_dir) +dir.create(out_dir, showWarnings = FALSE, recursive = TRUE) + +country <- "UA" +model_name <- "397_covs_model" + +idx = read.csv(file.path(out_dir, paste0(tolower(country), "_master_index", ".csv"))) +idx_F = read.csv(file.path(out_dir, "population_proxy", "social_media_audience", paste0(tolower(country), "_facebook_audience", ".csv"))) +idx_G = read.csv(file.path(out_dir, "population_proxy", "social_media_audience", paste0(tolower(country), "_instagram_audience", ".csv"))) +covs = read.csv(file.path(out_dir, "covariates", "final", "ua_covariates_oblast.csv")) +codps = read.csv(file.path(data_dir, "cod-ps", "population_baseline.csv")) +codps_N1 = read.csv(file.path(in_dir, "cod-ps_2023", "DO_NOT_SHARE_UKR_ADM2_POP_2023.csv")) +codps_N2 = read.csv(file.path(in_dir, "cod-ps_2024", "DO_NOT_SHARE_UKR_ADM2_POP_2024_Sept_27.csv")) +outside_border = read.csv(file.path(out_dir, "population_proxy", "crossing_borders", "dat_refugees.csv")) +last_date = "2024-05-14" +process_drop_locations = c() +observation_drop_locations = c(3782,3788,3791,3797) #3782, 3791, 3797, 3788 +process_cov_select = read.csv(file.path(in_dir, "ua_r_covariates_oblast_selection.csv")) +observation_cov_select = read.csv(file.path(in_dir, "ua_p_covariates_oblast_selection.csv")) + +md <- list() + +seed <- round(runif(1, 1, 1e6)) +set.seed(seed) +md$seed <- seed + +combined_drop_locations <- c(process_drop_locations, observation_drop_locations) + +selected_locations <- as.integer(setdiff(seq_len(27), c(11, 5, 13, 19))) + +selected_locations_key <- as.integer(c( + 3800, 3801, 3781, 3804, 3802, 3803, 3783, 3790, 3787, + 3792, 3793, 3794, 3795, 3796, 3798, 3799, 3784, 3785, + 3786, 3778, 3780, 3779, 4290 +)) + +selected_locations_pcode <- c( + #"UA01", "UA14", "UA44", "UA85", + "UA71", "UA74", "UA73", "UA12", "UA26", "UA63", "UA65", "UA68", "UA35", + "UA32", "UA46", "UA48", "UA51", "UA53", "UA56", "UA59", "UA61", "UA05", + "UA07", "UA21", "UA23", "UA18", "UA80" +) + +i_idx <- idx %>% + select(i_key) %>% + distinct() %>% + filter(i_key %in% selected_locations_key) %>% + arrange(i_key) %>% + mutate(i = row_number()) + +md$idx <- idx %>% + filter(i %in% selected_locations, a != 7) %>% + select( + t, t_key, t_name, i_key, i_name, + a, a_key, a_name, s, s_key, s_name + ) %>% + left_join(i_idx, by = "i_key") %>% + arrange(t, i, a, s) %>% + mutate(tias = row_number()) %>% + group_by(a, s) %>% + mutate(ti = row_number()) %>% + ungroup() %>% + group_by(t, i) %>% + mutate(as = row_number()) %>% + ungroup() %>% + filter(t_name <= floor_date(as.Date(last_date), "week", week_start = 1)) %>% + select( + ti, as, tias, t, i, a, s, + t_key, t_name, i_key, i_name, + a_key, a_name, s_key, s_name + ) + +md$A <- 7 +md$S <- 2 +md$I <- 23 +md$T <- 117 +md$C_full <- md$I * md$A * md$S +md$C_adult <- md$I * (md$A - 1) * md$S +md$C_young <- md$I * md$S + +adult_map <- expand.grid( + t = seq_len(md$T), + i = seq_len(md$I), + a = 1:(md$A - 1), + s = seq_len(md$S) +) %>% + arrange(t, i, a, s) %>% + mutate(tias_adult = row_number()) %>% + select(t, i, a, s, tias_adult) + +idx_keys <- idx %>% + select(i, i_key, a, a_key, s, s_key) %>% + distinct() + +md$idx_F <- idx_F %>% + filter(i %in% selected_locations, a < 6) %>% + select(t, i, a, s, m, value) %>% + left_join(idx_keys, by = c("i", "a", "s")) %>% + select(-i) %>% + left_join(i_idx, by = "i_key") %>% + right_join(md$idx %>% select(t, i, a, s, ti, as, tias), by = c("t", "i", "a", "s")) %>% + left_join(adult_map, by = c("t", "i", "a", "s")) %>% + filter(!i_key %in% combined_drop_locations, value > 0, is.finite(value)) %>% + group_by(ti, as, tias, tias_adult, t, i, a, s) %>% + summarise(value = mean(value, na.rm = TRUE), .groups = "drop") %>% + select(ti, as, tias_full = tias, tias_adult, t, i, a, s, value) %>% + arrange(t, i, a, s) + +md$idx_G <- idx_G %>% + filter(i %in% selected_locations, a < 6) %>% + select(t, i, a, s, m, value) %>% + left_join(idx_keys, by = c("i", "a", "s")) %>% + select(-i) %>% + left_join(i_idx, by = "i_key") %>% + right_join(md$idx %>% select(t, i, a, s, ti, as, tias), by = c("t", "i", "a", "s")) %>% + left_join(adult_map, by = c("t", "i", "a", "s")) %>% + filter(!i_key %in% combined_drop_locations, value > 0, is.finite(value)) %>% + group_by(ti, as, tias, tias_adult, t, i, a, s) %>% + summarise(value = mean(value, na.rm = TRUE), .groups = "drop") %>% + select(ti, as, tias_full = tias, tias_adult, t, i, a, s, value) %>% + arrange(t, i, a, s) + +md$y_F <- md$idx_F$value +md$n_F <- length(md$y_F) +md$tias_F <- md$idx_F$tias_adult + +md$y_G <- md$idx_G$value +md$n_G <- length(md$y_G) +md$tias_G <- md$idx_G$tias_adult + +cols_to_pivot <- c( + "F_00_04", "F_05_09", "F_10_14", "F_15_19", "F_20_24", "F_25_29", "F_30_34", "F_35_39", + "F_40_44", "F_45_49", "F_50_54", "F_55_59", "F_60_64", "F_65_69", "F_70_74", "F_75_79", + "F_80Plus", + "M_00_04", "M_05_09", "M_10_14", "M_15_19", "M_20_24", "M_25_29", "M_30_34", "M_35_39", + "M_40_44", "M_45_49", "M_50_54", "M_55_59", "M_60_64", "M_65_69", "M_70_74", "M_75_79", + "M_80Plus" +) + +date_start <- "2022-02-25" +date_end <- last_date + +time_index_expanded <- tibble( + collection_date = seq(as.Date(date_start), as.Date(date_end), by = 1) +) |> + mutate( + t_name = floor_date(as.Date(collection_date), "week", week_start = 1), + t_key = stringr::str_replace_all(as.character(t_name), "-", "") |> as.integer(), + t = t_name |> as.character() |> haven::as_factor() |> as.integer() + ) + +time_index <- time_index_expanded |> + distinct(t_name, t_key, t) |> + arrange(t) + +############### +# COD-PS 2022 +############### +md$N010 <- codps %>% + select(ADM1_PCODE, fb_key, all_of(cols_to_pivot)) %>% + pivot_longer( + cols = all_of(cols_to_pivot), + names_to = "demog", + values_to = "pop0" + ) |> + separate(demog, into = c("sex", "age_min", "age_max"), sep = "_") |> + filter(sex != "T", age_min != "TL") |> + rename(pcode = ADM1_PCODE, i_key = fb_key) |> + filter(i_key %in% selected_locations_key) %>% + mutate( + age_group10y = case_when( + age_min == "00" & age_max == "04" ~ "00-14", + age_min == "05" & age_max == "09" ~ "00-14", + age_min == "10" & age_max == "14" ~ "00-14", + age_min == "15" & age_max == "19" ~ "15-19", + age_min == "20" & age_max == "24" ~ "20-29", + age_min == "25" & age_max == "29" ~ "20-29", + age_min == "30" & age_max == "34" ~ "30-39", + age_min == "35" & age_max == "39" ~ "30-39", + age_min == "40" & age_max == "44" ~ "40-49", + age_min == "45" & age_max == "49" ~ "40-49", + age_min == "50" & age_max == "54" ~ "50-59", + age_min == "55" & age_max == "59" ~ "50-59", + age_min == "60" & age_max == "64" ~ "60-999", + age_min == "65" & age_max == "69" ~ "60-999", + age_min == "70" & age_max == "74" ~ "60-999", + age_min == "75" & age_max == "79" ~ "60-999", + age_min == "80Plus" ~ "60-999" + ) + ) |> + group_by(i_key, pcode, age_group10y, sex) |> + summarise(value = sum(pop0), .groups = "drop") |> + mutate( + a = case_when( + age_group10y == "00-14" ~ 7L, + age_group10y == "15-19" ~ 6L, + age_group10y == "20-29" ~ 1L, + age_group10y == "30-39" ~ 2L, + age_group10y == "40-49" ~ 3L, + age_group10y == "50-59" ~ 4L, + age_group10y == "60-999" ~ 5L + ), + s = if_else(sex == "F", 2L, 1L) + ) |> + left_join(i_idx, by = "i_key") |> + arrange(i, a, s) |> + mutate(ias = row_number()) %>% + group_by(i, s) %>% + mutate(is = row_number()) %>% + ungroup() %>% + filter(!is.na(a)) %>% + select(ias, is, i, a, s, value) + +md$N0 <- md$N010$value +md$ias_N0 <- md$N010$ias + +md$Y0 <- md$N010 |> + arrange(i, a, s) |> + filter(a == 7) %>% + pull(value) + +md$A0 <- md$N010 |> + arrange(i, a, s) |> + filter(a != 7) %>% + pull(value) + +md$N01 <- codps |> + filter(fb_key %in% selected_locations_key) %>% + arrange(match(fb_key, i_idx$i_key)) |> + pull(T_TL) + +full_pop22 <- codps %>% + select(fb_key, all_of(cols_to_pivot)) %>% + pivot_longer( + cols = all_of(cols_to_pivot), + names_to = c("sex", "age_min", "age_max"), + names_sep = "_", + values_to = "pop0" + ) %>% + filter(sex != "T", age_min != "TL") + +total_pop_all22 <- sum(full_pop22$pop0, na.rm = TRUE) +total_pop_sel22 <- full_pop22 %>% + filter(fb_key %in% selected_locations_key) %>% + summarise(sum_sel = sum(pop0, na.rm = TRUE)) %>% + pull(sum_sel) + +prop_sel22 <- total_pop_sel22 / total_pop_all22 +prop_sel22 + +############## +# COD-PS 2023 +############## +md$N110 <- codps_N1 %>% + select(ADM1_PCODE, ADM2_PCODE, all_of(cols_to_pivot)) %>% + pivot_longer( + cols = all_of(cols_to_pivot), + names_to = "demog", + values_to = "pop1" + ) |> + group_by(ADM1_PCODE, demog) |> + summarise(pop1 = sum(pop1, na.rm = TRUE), .groups = "drop") |> + left_join(codps[, c("ADM1_PCODE", "fb_key")], by = "ADM1_PCODE") %>% + separate(demog, into = c("sex", "age_min", "age_max"), sep = "_") |> + filter(sex != "T", age_min != "TL") |> + rename(pcode = ADM1_PCODE, i_key = fb_key) |> + filter(i_key %in% selected_locations_key) %>% + mutate( + age_group10y = case_when( + age_min == "00" & age_max == "04" ~ "00-14", + age_min == "05" & age_max == "09" ~ "00-14", + age_min == "10" & age_max == "14" ~ "00-14", + age_min == "15" & age_max == "19" ~ "15-19", + age_min == "20" & age_max == "24" ~ "20-29", + age_min == "25" & age_max == "29" ~ "20-29", + age_min == "30" & age_max == "34" ~ "30-39", + age_min == "35" & age_max == "39" ~ "30-39", + age_min == "40" & age_max == "44" ~ "40-49", + age_min == "45" & age_max == "49" ~ "40-49", + age_min == "50" & age_max == "54" ~ "50-59", + age_min == "55" & age_max == "59" ~ "50-59", + age_min == "60" & age_max == "64" ~ "60-999", + age_min == "65" & age_max == "69" ~ "60-999", + age_min == "70" & age_max == "74" ~ "60-999", + age_min == "75" & age_max == "79" ~ "60-999", + age_min == "80Plus" ~ "60-999" + ) + ) %>% + group_by(i_key, pcode, age_group10y, sex) |> + summarise(value = sum(pop1), .groups = "drop") |> + mutate( + a = case_when( + age_group10y == "00-14" ~ 7L, + age_group10y == "15-19" ~ 6L, + age_group10y == "20-29" ~ 1L, + age_group10y == "30-39" ~ 2L, + age_group10y == "40-49" ~ 3L, + age_group10y == "50-59" ~ 4L, + age_group10y == "60-999" ~ 5L + ), + s = if_else(sex == "F", 2L, 1L) + ) |> + left_join(i_idx, by = "i_key") |> + arrange(i, a, s) |> + mutate(ias = row_number()) %>% + group_by(i, s) %>% + mutate(is = row_number()) %>% + ungroup() %>% + filter(!is.na(a)) %>% + select(ias, is, i, a, s, value) + +md$N1 <- md$N110$value +md$ias_N1 <- md$N110$ias + +md$Y1 <- md$N110 |> + arrange(i, a, s) |> + filter(a == 7) %>% + pull(value) + +md$A1 <- md$N110 |> + arrange(i, a, s) |> + filter(a != 7) %>% + pull(value) + +md$N11 <- codps_N1 %>% + select(ADM1_PCODE, ADM2_PCODE, all_of(cols_to_pivot)) %>% + pivot_longer( + cols = all_of(cols_to_pivot), + names_to = "demog", + values_to = "pop1" + ) |> + group_by(ADM1_PCODE) |> + summarise(pop1 = sum(pop1, na.rm = TRUE), .groups = "drop") |> + left_join(codps[, c("ADM1_PCODE", "fb_key")], by = "ADM1_PCODE") %>% + filter(fb_key %in% selected_locations_key) %>% + arrange(match(fb_key, i_idx$i_key)) %>% + pull(pop1) + +pop1_date <- as.Date("2023-07-01") +pop1_week <- floor_date(pop1_date, "week", week_start = 1) +t_pop1 <- time_index %>% + filter(t_name == pop1_week) %>% + pull(t) +t_pop1 + +full_pop23 <- codps_N1 %>% + select(ADM1_PCODE, ADM2_PCODE, all_of(cols_to_pivot)) %>% + pivot_longer( + cols = all_of(cols_to_pivot), + names_to = "demog", + values_to = "pop1" + ) %>% + group_by(ADM1_PCODE, demog) %>% + summarise(pop1 = sum(pop1, na.rm = TRUE), .groups = "drop") %>% + separate(demog, into = c("sex", "age_min", "age_max"), sep = "_") %>% + filter(sex != "T", age_min != "TL") + +total_pop_all23 <- sum(full_pop23$pop1, na.rm = TRUE) + +full_pop23_sel <- full_pop23 %>% + left_join(codps %>% select(ADM1_PCODE, fb_key), by = "ADM1_PCODE") %>% + filter(fb_key %in% selected_locations_key) + +total_pop_sel23 <- sum(full_pop23_sel$pop1, na.rm = TRUE) +prop_sel23 <- total_pop_sel23 / total_pop_all23 +prop_sel23 + +################ +# COD-PS 2024 +################ +md$N210 <- codps_N2 %>% + select(ADM1_PCODE, ADM2_PCODE, all_of(cols_to_pivot)) %>% + pivot_longer( + cols = all_of(cols_to_pivot), + names_to = "demog", + values_to = "pop2" + ) %>% + group_by(ADM1_PCODE, demog) %>% + summarise(pop2 = sum(pop2, na.rm = TRUE), .groups = "drop") %>% + left_join(codps[, c("ADM1_PCODE", "fb_key")], by = "ADM1_PCODE") %>% + separate(demog, into = c("sex", "age_min", "age_max"), sep = "_") %>% + filter(sex != "T", age_min != "TL") %>% + rename(pcode = ADM1_PCODE, i_key = fb_key) %>% + filter(i_key %in% selected_locations_key) %>% + mutate( + age_group10y = case_when( + age_min == "00" & age_max == "04" ~ "00-14", + age_min == "05" & age_max == "09" ~ "00-14", + age_min == "10" & age_max == "14" ~ "00-14", + age_min == "15" & age_max == "19" ~ "15-19", + age_min == "20" & age_max == "24" ~ "20-29", + age_min == "25" & age_max == "29" ~ "20-29", + age_min == "30" & age_max == "34" ~ "30-39", + age_min == "35" & age_max == "39" ~ "30-39", + age_min == "40" & age_max == "44" ~ "40-49", + age_min == "45" & age_max == "49" ~ "40-49", + age_min == "50" & age_max == "54" ~ "50-59", + age_min == "55" & age_max == "59" ~ "50-59", + age_min == "60" & age_max == "64" ~ "60-999", + age_min == "65" & age_max == "69" ~ "60-999", + age_min == "70" & age_max == "74" ~ "60-999", + age_min == "75" & age_max == "79" ~ "60-999", + age_min == "80Plus" ~ "60-999" + ) + ) %>% + group_by(i_key, pcode, age_group10y, sex) %>% + summarise(value = sum(pop2), .groups = "drop") %>% + mutate( + a = case_when( + age_group10y == "00-14" ~ 7L, + age_group10y == "15-19" ~ 6L, + age_group10y == "20-29" ~ 1L, + age_group10y == "30-39" ~ 2L, + age_group10y == "40-49" ~ 3L, + age_group10y == "50-59" ~ 4L, + age_group10y == "60-999" ~ 5L + ), + s = if_else(sex == "F", 2L, 1L) + ) %>% + left_join(i_idx, by = "i_key") %>% + arrange(i, a, s) %>% + mutate(ias = row_number()) %>% + group_by(i, s) %>% + mutate(is = row_number()) %>% + ungroup() %>% + filter(!is.na(a)) %>% + select(ias, is, i, a, s, value) + +md$N2 <- md$N210$value +md$ias_N2 <- md$N210$ias + +md$Y2 <- md$N210 |> + arrange(i, a, s) |> + filter(a == 7) %>% + pull(value) + +md$A2 <- md$N210 |> + arrange(i, a, s) |> + filter(a != 7) %>% + pull(value) + +md$N21 <- codps_N2 %>% + select(ADM1_PCODE, ADM2_PCODE, all_of(cols_to_pivot)) %>% + pivot_longer( + cols = all_of(cols_to_pivot), + names_to = "demog", + values_to = "pop2" + ) |> + group_by(ADM1_PCODE) |> + summarise(pop2 = sum(pop2, na.rm = TRUE), .groups = "drop") |> + left_join(codps[, c("ADM1_PCODE", "fb_key")], by = "ADM1_PCODE") %>% + arrange(match(fb_key, i_idx$i_key)) %>% + pull(pop2) + +pop2_date <- as.Date("2024-05-14") +pop2_week <- floor_date(pop2_date, "week", week_start = 1) +t_pop2 <- time_index %>% + filter(t_name == pop2_week) %>% + pull(t) +t_pop2 + +full_pop24 <- codps_N2 %>% + select(ADM1_PCODE, ADM2_PCODE, all_of(cols_to_pivot)) %>% + pivot_longer( + cols = all_of(cols_to_pivot), + names_to = "demog", + values_to = "pop2" + ) %>% + group_by(ADM1_PCODE, demog) %>% + summarise(pop2 = sum(pop2, na.rm = TRUE), .groups = "drop") %>% + separate(demog, into = c("sex", "age_min", "age_max"), sep = "_") %>% + filter(sex != "T", age_min != "TL") + +total_pop_all24 <- sum(full_pop24$pop2, na.rm = TRUE) + +full_pop24_sel <- full_pop24 %>% + left_join(codps %>% select(ADM1_PCODE, fb_key), by = "ADM1_PCODE") %>% + filter(fb_key %in% selected_locations_key) + +total_pop_sel24 <- sum(full_pop24_sel$pop2, na.rm = TRUE) +prop_sel24 <- total_pop_sel24 / total_pop_all24 +prop_sel24 + +###################################### +# Total population at every time step +###################################### +weekly_avg <- outside_border |> + mutate(week = floor_date(as.Date(date), "week", week_start = 1)) |> + group_by(week) |> + summarise(avg_value = mean(individuals, na.rm = TRUE), .groups = "drop") |> + filter( + week >= min(md$idx$t_name), + week <= max(md$idx$t_name) + ) + +border_df <- time_index %>% + select(t, t_name) %>% + left_join(weekly_avg %>% rename(t_name = week), by = "t_name") %>% + mutate(avg_value = if_else(t == 1, 0, avg_value)) + +md$y_N_tot <- border_df %>% + mutate( + baseline = case_when( + t < t_pop1 ~ sum(md$N01), + t == t_pop1 ~ sum(md$N11), + t < t_pop2 ~ sum(md$N11), + t == t_pop2 ~ sum(md$N21), + TRUE ~ sum(md$N21) + ), + prop = case_when( + t <= t_pop1 ~ prop_sel22, + t_pop1 < t & t <= t_pop2 ~ prop_sel23, + TRUE ~ prop_sel24 + ), + ytmp = baseline - (avg_value * prop), + y_N_tot = as.integer(ytmp) + ) %>% + pull(y_N_tot) + +rm(weekly_avg) + +children <- md$N010 %>% + filter(a == 7) %>% + mutate(t = 1) %>% + rename(children = value) %>% + rbind( + md$N110 %>% filter(a == 7) %>% mutate(t = t_pop1) %>% rename(children = value) + ) %>% + rbind( + md$N210 %>% filter(a == 7) %>% mutate(t = t_pop2) %>% rename(children = value) + ) + +female_repo <- md$N010 %>% + filter(a >= 1 & a <= 3, s == 2) %>% + group_by(i) %>% + summarise(female_repo = sum(value, na.rm = TRUE), .groups = "drop") %>% + mutate(t = 1) %>% + rbind( + md$N110 %>% + filter(a >= 1 & a <= 3, s == 2) %>% + group_by(i) %>% + summarise(female_repo = sum(value, na.rm = TRUE), .groups = "drop") %>% + mutate(t = t_pop1) + ) %>% + rbind( + md$N210 %>% + filter(a >= 1 & a <= 3, s == 2) %>% + group_by(i) %>% + summarise(female_repo = sum(value, na.rm = TRUE), .groups = "drop") %>% + mutate(t = t_pop2) + ) + +md$ratio_children_female <- expand.grid( + t = seq_len(md$T), + i = seq_len(md$I), + s = seq_len(md$S) +) %>% + left_join(children) %>% + left_join(female_repo) %>% + mutate(ratio_children_female = children / female_repo) %>% + group_by(i, s) %>% + arrange(t, .by_group = TRUE) %>% + mutate(ratio_children_female_im = zoo::na.approx(ratio_children_female, x = t, na.rm = FALSE)) %>% + ungroup() %>% + arrange(t, i, s) %>% + pull(ratio_children_female_im) + +table(!is.na(md$ratio_children_female)) + +t <- unique(md$idx$t) +i <- unique(md$idx$i) +a <- seq_len(md$A - 1) +s <- unique(md$idx$s) + +I <- md$I +A <- md$A +S <- md$S +T <- md$T + +C_full <- I * A * S +C_adult <- I * (A - 1) * S +C_young <- I * S +IS <- I * S + +adult_grid <- expand.grid(t = 1:T, i = 1:I, a = 1:(A - 1), s = 1:S) |> + dplyr::arrange(t, i, a, s) + +md$tt <- adult_grid$t +md$ii <- adult_grid$i +md$aa <- adult_grid$a +md$ss <- adult_grid$s + +md$adult_grid <- adult_grid[order(adult_grid$t, adult_grid$i, adult_grid$a, adult_grid$s), ] +md$female_core_mask <- as.integer(md$adult_grid$a %in% 1:7 & md$adult_grid$s == 2) +md$adult_grid$idx_tis <- (md$adult_grid$t - 1) * I * S + (md$adult_grid$i - 1) * S + md$adult_grid$s +md$tis_index <- md$adult_grid$idx_tis + +full_grid <- expand.grid(t = 1:T, i = 1:I, a = 1:A, s = 1:S) %>% + arrange(t, i, a, s) %>% + mutate(idx_full = row_number()) + +md$young_pos_full <- full_grid %>% + filter(t == 1, a == 7) %>% + arrange(i, a, s) %>% + pull(idx_full) %>% + as.integer() + +md$young_idx <- full_grid %>% + filter(a == 7) %>% + arrange(t, i, a, s) %>% + pull(idx_full) %>% + as.integer() + +md$adult_pos_full <- full_grid %>% + filter(t == 1, a != 7) %>% + arrange(i, a, s) %>% + pull(idx_full) %>% + as.integer() + +md$adult_idx <- full_grid %>% + filter(a != 7) %>% + arrange(t, i, a, s) %>% + pull(idx_full) %>% + as.integer() + +md$N_y <- T * length(md$young_pos_full) +md$N_a <- T * length(md$adult_pos_full) + +md$slice_full <- matrix(seq_len(T * C_full), nrow = T, byrow = TRUE) +md$slice_adult <- matrix(seq_len(T * C_adult), nrow = T, byrow = TRUE) +md$slice_adult_lag <- rbind(rep(1L, C_adult), md$slice_adult[-T, ]) + +md$slice_young <- matrix(seq_len(T * C_young), nrow = T, byrow = TRUE) +md$slice_young_lag <- rbind(rep(1L, C_young), md$slice_young[-T, ]) + +f_core <- (md$ss == 2L) & (md$aa %in% 1:3) +md$idx_female_a234 <- which(f_core) +md$idx_female_triple <- matrix(md$idx_female_a234, ncol = 3, byrow = TRUE) +md$N_fem_a <- length(md$idx_female_a234) + +idx_a1519 <- adult_map %>% filter(a == 6) %>% arrange(t, i, s) %>% pull(tias_adult) +idx_a2029 <- adult_map %>% filter(a == 1) %>% arrange(t, i, s) %>% pull(tias_adult) + +md$idx_a1519 <- as.integer(idx_a1519) +md$idx_a2029 <- as.integer(idx_a2029) +md$N_tie <- length(idx_a1519) +length(idx_a2029) + +standardise_cols <- function(df, vars) { + stats <- vector("list", length(vars)) + names(stats) <- vars + + for (v in vars) { + x <- df[[v]] + m <- mean(x, na.rm = TRUE) + s <- sd(x, na.rm = TRUE) + + if (!is.finite(s) || s <= 0) { + df[[v]] <- 0 + stats[[v]] <- c(mean = m, sd = 1) + } else { + df[[v]] <- (x - m) / s + stats[[v]] <- c(mean = m, sd = s) + } + } + + list(df = df, stats = stats) +} + +covs$time_std[is.na(covs$time_std)] <- "none" +covs$space_std[is.na(covs$space_std)] <- "none" + +rcov_select <- process_cov_select %>% + filter(select == 1) %>% + mutate(across(c(time_std, space_std), ~ replace(., is.na(.), "none"))) + +cov_names_r <- rcov_select %>% + transmute(nice_name = paste(covariate, sum_stat, time_std, space_std, sep = "_")) %>% + pull() + +idx <- expand.grid( + t = 1:md$T, + i = 1:md$I, + a = min(md$aa):max(md$aa), + s = 1:md$S +) %>% + arrange(t, i, a, s) + +Xr_df <- idx +for (k in seq_along(cov_names_r)) { + sel <- rcov_select[k, ] + col_real <- cov_names_r[k] + + cov_k <- covs %>% + filter( + covariate == sel$covariate, + sum_stat == sel$sum_stat, + time_std == sel$time_std, + space_std == sel$space_std + ) %>% + select(t, i, value_std) %>% + rename(!!col_real := value_std) + + Xr_df <- left_join(Xr_df, cov_k, by = c("t", "i")) +} + +Xr_df[is.na(Xr_df)] <- 0 + +std_r <- standardise_cols(Xr_df, cov_names_r) +Xr_df <- std_r$df +md$X_r_means <- sapply(std_r$stats, `[[`, "mean") +md$X_r_sds <- sapply(std_r$stats, `[[`, "sd") + +md$X_r <- as.matrix(Xr_df["acled_withfatalities_raw_none_none"]) +md$K_r <- ncol(md$X_r) + +pcov_select <- observation_cov_select %>% + filter(select == 1) %>% + mutate(across(c(time_std, space_std), ~ replace(., is.na(.), "none"))) + +cov_names_p <- pcov_select %>% + transmute(nice_name = paste(covariate, sum_stat, time_std, space_std, sep = "_")) %>% + pull() + +Xp_df <- idx +for (k in seq_along(cov_names_p)) { + sel <- pcov_select[k, ] + col_real <- cov_names_p[k] + + cov_k <- covs %>% + filter( + covariate == sel$covariate, + sum_stat == sel$sum_stat, + time_std == sel$time_std, + space_std == sel$space_std + ) %>% + select(t, i, value_std) %>% + rename(!!col_real := value_std) + + Xp_df <- left_join(Xp_df, cov_k, by = c("t", "i")) +} + +Xp_df[is.na(Xp_df)] <- 0 + +std_p <- standardise_cols(Xp_df, cov_names_p) +Xp_df <- std_p$df +md$X_p_means <- sapply(std_p$stats, `[[`, "mean") +md$X_p_sds <- sapply(std_p$stats, `[[`, "sd") + +md$X_p <- as.matrix(Xp_df["pwtt_sum_24week_none_none"]) +md$K_p <- ncol(md$X_p) + +saveRDS(md, file.path(out_dir, "modelling", model_name, "mcmc", paste0("md_", model_name, ".rds"))) diff --git a/src/sandbox/multi-source_statistical/models_andrea/397_covs_model/README.md b/src/sandbox/multi-source_statistical/models_andrea/397_covs_model/README.md new file mode 100644 index 00000000..137143d5 --- /dev/null +++ b/src/sandbox/multi-source_statistical/models_andrea/397_covs_model/README.md @@ -0,0 +1,54 @@ +# Beyond crisis nowcasting: A Bayesian framework for integrating multiple social media sources into population estimates + +This repository contains code and documentation to reproduce the modelling workflow for the paper: + +**Beyond crisis nowcasting: A Bayesian framework for integrating multiple social media sources into population estimates** + +Target journal: **Proceedings of the National Academy of Sciences of the United States of America (PNAS)** + +The case study is Ukraine (oblast level). The workflow fits a Bayesian model that integrates Facebook and Instagram audience data into weekly age and sex specific population estimates. + +## Data notes and access constraints: + +Some inputs used by the scripts are restricted and must not be committed to a public repository. In particular, the COD-PS 2023 and 2024. + +## Software requirements + +- R (recent version) +- CmdStan and the R interface `cmdstanr` +- R packages used by the modelling and evaluation scripts include: `cmdstanr`, `posterior`, `bayesplot`, `ggplot2`, `tidyverse`, `tidyr`, `readxl`, `patchwork`, `purrr`, and others used by helper scripts. + +## Configuration via .env + +The modelling scripts load a `.env` file, where key directories such as `repo_dir`, `in_dir`, and `out_dir` are defined. + +Create a `.env` file at the repository root. Use `example.env` as a template. + +## Workflow: + +**Step 1. Build population proxy inputs**: +Run these scripts (in order) from the folder `src/data_prep/pop_data/1_pop_data/`: + +- `10_create_master_index.R` +- `20_prepare_sma_data.R` +- `30_border_crossing.R` + +These scripts create outputs for: `out_dir/population_proxy/` + +**Step 2. Fit the model 397_covs_model**: +This workflow is currently model specific. Therefore, the Stan model file named `397_covs_model.stan` alongside all model specific files must be in the same folder. They are: + +- `397_covs_model_config.R`: Prepares the data list `md` used as input to Stan, including index construction, anchor populations, and derived ratios used to impute child counts from women of reproductive ages. +- `init_generator_fun_397.R`: Creates initial values for sampling with checks on expected dimensions and seeded initialisation by chain identifier. +- `1_mcmc_397.R`: Runs Bayesian sampling via `cmdstanr`, with threading enabled and multiple chains. +- `2_eval_397.R`: Produces evaluation outputs including trace plots, posterior predictive checks, time series plots, and population pyramid plots. + +If you run `1_mcmc_397.R` directly, it will source the model configuration script `397_covs_model_config.R`. + +## Outputs +The evaluation script expects these files to exist: +- `out_dir/modelling/397_covs_model/mcmc/md_397_covs_model.rds`: input data created by `397_covs_model_config.R`. +- `out_dir/modelling/397_covs_model/mcmc/fit_397_covs_model.rds`: fitted model create by `1_mcmc_397.R`. + +Bear in mind that the resulted MCMC file (`fit_397_covs_model.rds`) is around 6.3MB. The model takes nearly 8.4 hours to run with settings specified in `1_mcmc_397.R`. + diff --git a/src/sandbox/multi-source_statistical/models_andrea/397_covs_model/init_generator_fun_397.R b/src/sandbox/multi-source_statistical/models_andrea/397_covs_model/init_generator_fun_397.R new file mode 100644 index 00000000..99e52d1c --- /dev/null +++ b/src/sandbox/multi-source_statistical/models_andrea/397_covs_model/init_generator_fun_397.R @@ -0,0 +1,112 @@ +init_generator <- function(md, chain_id = 1) { + set.seed(2025 + chain_id) + + T <- md$T + I <- md$I + A <- md$A + S <- md$S + K_r <- md$K_r + K_p <- md$K_p + + if (is.null(md$C_adult)) { + stop("md$C_adult is NULL; please pass C_adult in md (= I*(A-1)*S).") + } + C_adult <- md$C_adult + N_cells <- T * C_adult + + A_free <- max(A - 2, 0) + Tm1 <- max(T - 1, 0) + + if (!is.null(md$tt) && length(md$tt) != N_cells) { + stop(sprintf("length(md$tt)=%d but T*C_adult=%d (T=%d, C_adult=%d).", + length(md$tt), N_cells, T, C_adult)) + } + if (!is.null(md$X_r) && nrow(md$X_r) != N_cells) { + stop(sprintf("nrow(md$X_r)=%d but T*C_adult=%d (T=%d, C_adult=%d).", + nrow(md$X_r), N_cells, T, C_adult)) + } + if (!is.null(md$X_p) && nrow(md$X_p) != N_cells) { + stop(sprintf("nrow(md$X_p)=%d but T*C_adult=%d (T=%d, C_adult=%d).", + nrow(md$X_p), N_cells, T, C_adult)) + } + + draw_sigma <- function(lo, hi) runif(1, lo, hi) + + beta_r <- if (K_r > 0) rep(0, K_r) else numeric(0) + beta_p <- if (K_p > 0) rep(0, K_p) else numeric(0) + + alpha_r <- 0 + alpha_p <- 0 + phi_p <- 0 + + sigma_F <- runif(1, 0.05, 0.20) + sigma_G <- runif(1, 0.05, 0.20) + + sigma_r <- runif(1, 0.005, 0.03) + + sigma_delta_r <- draw_sigma(0.005, 0.05) + sigma_gamma_r_i <- draw_sigma(0.005, 0.05) + sigma_gamma_r_a <- draw_sigma(0.005, 0.05) + sigma_gamma_r_s <- draw_sigma(0.005, 0.05) + + sigma_delta_p <- draw_sigma(0.005, 0.05) + sigma_gamma_p_i <- draw_sigma(0.005, 0.05) + sigma_gamma_p_a <- draw_sigma(0.005, 0.05) + sigma_gamma_p_s <- draw_sigma(0.005, 0.05) + + sigma_phi_p_a <- draw_sigma(0.005, 0.05) + + delta_r_innov_raw <- if (Tm1 > 0) rnorm(Tm1) else numeric(0) + delta_p_innov_raw <- if (Tm1 > 0) rnorm(Tm1) else numeric(0) + + gamma_r_i_raw <- if (I > 0) rnorm(I) else numeric(0) + gamma_p_i_raw <- if (I > 0) rnorm(I) else numeric(0) + + gamma_r_s_raw <- if (S > 0) rnorm(S) else numeric(0) + gamma_p_s_raw <- if (S > 0) rnorm(S) else numeric(0) + + gamma_r_a_raw <- if (A_free > 0) rnorm(A_free) else numeric(0) + gamma_p_a_raw <- if (A_free > 0) rnorm(A_free) else numeric(0) + + phi_p_a_raw <- if (A_free > 0) rnorm(A_free) else numeric(0) + + log_r_raw <- rnorm(N_cells, mean = 0, sd = 1) + log_r_raw <- pmin(pmax(log_r_raw, -3), 3) + + list( + log_r_raw = log_r_raw, + sigma_r = sigma_r, + alpha_r = alpha_r, + beta_r = beta_r, + + sigma_delta_r = sigma_delta_r, + sigma_gamma_r_i = sigma_gamma_r_i, + sigma_gamma_r_a = sigma_gamma_r_a, + sigma_gamma_r_s = sigma_gamma_r_s, + + delta_r_innov_raw = delta_r_innov_raw, + gamma_r_i_raw = gamma_r_i_raw, + gamma_r_a_raw = gamma_r_a_raw, + gamma_r_s_raw = gamma_r_s_raw, + + alpha_p = alpha_p, + beta_p = beta_p, + phi_p = phi_p, + + sigma_delta_p = sigma_delta_p, + sigma_gamma_p_i = sigma_gamma_p_i, + sigma_gamma_p_a = sigma_gamma_p_a, + sigma_gamma_p_s = sigma_gamma_p_s, + + delta_p_innov_raw = delta_p_innov_raw, + gamma_p_i_raw = gamma_p_i_raw, + gamma_p_a_raw = gamma_p_a_raw, + gamma_p_s_raw = gamma_p_s_raw, + + sigma_F = sigma_F, + sigma_G = sigma_G, + + sigma_phi_p_a = sigma_phi_p_a, + phi_p_a_raw = phi_p_a_raw + ) +}