Merge branch 'dev' into dev_update_postprocessing

Analysis/R/prepare_initial_values.R

@@ -297,55 +297,39 @@ if (config$obs_model == 3) {
   init.list <- purrr::map(
     1:nchain, function(x) {
       init <- list(
-        std_dev_w = runif(1, 1, 3),
-        rho = runif(1, .7, .8),
         # This assumes we have a fixed overdispersion parameter at admin level 1
         inv_od_param = abs(rnorm(stan_data$N_admin_lev - 1, 1, 1e-1))
       )
 
-      if (config$do_sd_w_mixture) {
+      if (config$use_intercept) {
         init <- append(
           init,
-          list(lambda = runif(1, .6, .8),
-               sigma_std_dev_w = abs(rnorm(2, .5, .1)))
-          )
-      } else {
-        init <- append(
-          init,
-          list(lambda = array(dim = 0),
-               sigma_std_dev_w = abs(rnorm(1, .5, .1)))
+          list(alpha = array(rnorm(1, -3, .5)))
         )
       }
 
-
-      if (config$time_effect) {
-        if (config$do_zerosum_cnst) {
-          n_eta <- stan_data$`T` - 1
-        } else {
-          n_eta <- stan_data$`T`
-        }
-
-        init <- append(
-          init,
-          list(eta_tilde = rnorm(n_eta, 0, .1))
-        )
-      }
-
-      # Spatial random effects
       if (config$spatial_effect) {
-        # Small values of spatial random effects
+
         init <- append(
           init,
-          list(w = rnorm(stan_data$smooth_grid_N, 0, .1))
-        )
+          list(std_dev_w = array(runif(1, 1, 3)),
+               rho = array(runif(1, .7, .8)),
+               w = rnorm(stan_data$smooth_grid_N, 0, .1))
+          )
 
-        if (config$use_intercept) {
+        if (config$do_sd_w_mixture) {
           init <- append(
             init,
-            list(alpha = rnorm(1, -3, .5))
+            list(lambda = array(runif(1, .6, .8)),
+                 sigma_std_dev_w = array(abs(rnorm(2, .5, .1))))
+          )
+        } else {
+          init <- append(
+            init,
+            list(lambda = array(dim = 0),
+                 sigma_std_dev_w = array(abs(rnorm(1, .5, .1))))
           )
         }
-
       } else {
 
         # Empty random effects
@@ -357,12 +341,20 @@ if (config$obs_model == 3) {
 
         if (!config$use_intercept) {
           warning("!! Runing a model with no space effect and no intercept.")
+        }
+      }
+
+      if (config$time_effect) {
+        if (config$do_zerosum_cnst) {
+          n_eta <- stan_data$`T` - 1
         } else {
-          init <- append(
-            init,
-            list(alpha = array(rnorm(1, 0, .1), dim = 1))
-          )
+          n_eta <- stan_data$`T`
         }
+
+        init <- append(
+          init,
+          list(eta_tilde = rnorm(n_eta, 0, .1))
+        )
       }
 
       init

Analysis/R/prepare_stan_input.R

@@ -185,8 +185,8 @@ prepare_stan_input <- function(
     dplyr::distinct(location_period_id, TL, TR) %>%
     dplyr::rename(locationPeriod_id = location_period_id) %>%
     dplyr::mutate(admin_lev = stringr::str_extract(locationPeriod_id, "ADM[0-9]{1}"),
-                  admin_lev = stringr::str_remove_all(admin_lev, "ADM") %>% as.integer())
-
+                  admin_lev = stringr::str_remove_all(admin_lev, "ADM") %>% as.integer()) %>% 
+    dplyr::arrange(locationPeriod_id)
 
   # Check that all fake observations appear in all time slices, drop if not
   fake_output_obs <- fake_output_obs %>% 
@@ -260,7 +260,7 @@ prepare_stan_input <- function(
       dplyr::pull(obs_id)
 
     # Get the corresponding mapping indices
-    output_adm0_locs <- purrr::map_dbl(output_adm0_obs, ~ which(stan_data$map_output_obs_loctime_obs == .))
+    output_adm0_locs <- purrr::map_dbl(output_adm0_obs, ~ stan_data$map_output_obs_loctime_loc[which(stan_data$map_output_obs_loctime_obs == .)])
 
     # Get the indices of the corresponding loctime/grid mapping
     output_adm0_ind <- purrr::map(output_adm0_locs, ~ which(stan_data$map_output_loc_grid_loc == .)) %>% 
@@ -269,7 +269,7 @@ prepare_stan_input <- function(
 
     # Check if all grid cells are covered
     if (nrow(sf_grid) != length(output_adm0_ind)) {
-      warning("Pop grid lengths do not match!")
+      warning("Pop grid lengths do not match! sf_grid:", nrow(sf_grid), " output cells:", length(output_adm0_ind))
     }
 
     # Get the grid cell ids
@@ -288,11 +288,12 @@ prepare_stan_input <- function(
 
     # Compute population adjustment factor
     adj_pop <- taxdat::compute_adjustment_UN_population(country = taxdat::get_country_isocode(config),
-                                                         pop = grid_pop,
-                                                         years = grid_years)
+                                                        pop = grid_pop,
+                                                        years = grid_years)
 
     # Probably a tidier way to do this but this should garantee the indexing is correct
-    stan_data$pop <- purrr::map_dbl(1:length(stan_data$pop), ~ stan_data$pop[.] * adj_pop[output_adm0_cells == .])
+    stan_data$pop[output_adm0_cells] <- purrr::map_dbl(1:length(output_adm0_cells), 
+                                                       ~ stan_data$pop[output_adm0_cells[.]] * adj_pop[.])
   }
 
   # ---- E. Pre-Aggregation Duplicates Removal in sf_cases ----
@@ -691,10 +692,10 @@ prepare_stan_input <- function(
   # We assume that the largest admin level (admin level 0 for national) has
   # an informative prior so as to produce little overdispersion. The over-dispersion for other
   # admin levels are allowed to have more prior support for larger amount of over-dispersion.
-  stan_data$mu_inv_od <- rep(0, stan_data$N_admin_lev)   # center at 0 (note that this is on the scale of 1/tau)
-  stan_data$sd_inv_od <- c(config$inv_od_sd_adm0,
-                           rep(config$inv_od_sd_nopool,
-                               stan_data$N_admin_lev - 1))
+  stan_data$mu_inv_od <- array(rep(0, stan_data$N_admin_lev))   # center at 0 (note that this is on the scale of 1/tau)
+  stan_data$sd_inv_od <- array(c(config$inv_od_sd_adm0,
+                                 rep(config$inv_od_sd_nopool,
+                                     stan_data$N_admin_lev - 1)))
 
   # Prior on the std_dev_w
   stan_data$mu_sd_w <- config$mu_sd_w