2024-09-03 update : adding spatial vignette. (#143)

* 2024-09-03 update : adding spatial vignette.

* 2024-09-03 update : fixing libraries in spatial vignette.

* 2024-09-04 update : cont. spatial vignette

* 2024-09-06 update : pushing latest.

* 2024-09-09 update : editing vignette.

* 2024-09-10 update : merging spatial-main.

* 2024-09-10 update : updating vignette.

* Kj tweaks to 123 (#159)

* modify seed and revert to old R(t) for testing, add inits

* fix length of norm vec aux site

* add identity matrix as initial cholesky decomp matrix

* tweak initialization

* tweak plots

* make prior on phi lower to approx iid case more

* 2024-09-11 update : testing spatial vignette.

---------

Co-authored-by: cbernalz <[email protected]>

* 2024-09-11 update : cont. vignette.

* 2024-09-13 update : vignette comp.

* solution to initialization, pass in 2 x2 matrix for L omega if corr structure switch !=2 (#165)

* 2024-09-16 update : resolving KJ requested changes.

* 2024-09-17 update : resolving KJ requested changes.

* Update _pkgdown.yml (#167)

* Update _pkgdown.yml

* add mathcal R(t) to getting started vignette

* 2024-09-17 update : merging KJ branch.

* 2024-09-17 update : resolving KJ requested changes.

---------

Co-authored-by: Kaitlyn Johnson <[email protected]>

DESCRIPTION

@@ -86,7 +86,10 @@ Imports:
     MASS,
     cramer,
     expm,
-    checkmate
+    checkmate,
+    patchwork,
+    scoringutils,
+    ggh4x
 Remotes:
     stan-dev/cmdstanr
 VignetteBuilder:

NAMESPACE

@@ -37,11 +37,12 @@ export(get_subpop_data)
 export(get_ww_data_indices)
 export(get_ww_data_sizes)
 export(get_ww_values)
-export(independence_corr_func)
+export(independence_corr_func_r)
 export(indicate_ww_exclusions)
 export(parameter_diagnostics)
 export(preprocess_count_data)
 export(preprocess_ww_data)
+export(rand_corr_matrix_func)
 export(spatial_rt_process)
 export(to_simplex)
 export(validate_paramlist)

R/generate_simulated_data.R

@@ -68,7 +68,7 @@
 #' defaulted to presets for exponential decay correlation function
 #' @param phi_rt Coefficient for AR(1) temporal correlation on subpopulation
 #' deviations
-#' @param sigma_generalized  Generalized variance of the spatial epsilon
+#' @param sigma_sqrd_generalized  Generalized variance of the spatial epsilon
 #' (determinant of the covariance matrix).
 #' @param scaling_factor Scaling factor for aux site
 #' @param aux_site_bool Boolean to use the aux site framework with
@@ -116,7 +116,7 @@
 #'     l = 1
 #'   ),
 #'   phi_rt = 0.6,
-#'   sigma_generalized = 0.05^6,
+#'   sigma_sqrd_generalized = 0.05^6,
 #'   scaling_factor = 1.1,
 #'   aux_site_bool = TRUE,
 #'   init_stat = TRUE
@@ -182,8 +182,8 @@ generate_simulated_data <- function(r_in_weeks = # nolint
                                       l = 1
                                     ),
                                     phi_rt = 0.6,
-                                    sigma_generalized = 0.05^4,
-                                    scaling_factor = 1.1,
+                                    sigma_sqrd_generalized = 0.005^4,
+                                    scaling_factor = 1,
                                     aux_site_bool = TRUE,
                                     init_stat = TRUE) {
   # Some quick checks to make sure the inputs work as expected-----------------
@@ -368,27 +368,14 @@ generate_simulated_data <- function(r_in_weeks = # nolint
       corr_fun_params$dist_matrix
     )
   }
-  sigma_matrix <- (sigma_generalized^(1 / n_sites)) * matrix_normalization(
-    corr_function(corr_fun_params)
+  corr_matrix <- corr_function(corr_fun_params)
+  sigma_matrix <- (sigma_sqrd_generalized^(1 / n_sites)) * matrix_normalization(
+    matrx = corr_matrix
   )
   spatial_deviation_noise_matrix <- spatial_deviation_noise_matrix_rng(
     sigma_matrix,
     n_weeks
   )
-  if (!use_spatial_corr) {
-    spatial_deviation_init <- mvrnorm(
-      n = 1,
-      mu = rep(0, n_sites + 1),
-      Sigma = sigma_matrix
-    )
-  } else {
-    spatial_deviation_init <- mvrnorm(
-      n = 1,
-      mu = rep(0, n_sites),
-      Sigma = sigma_matrix
-    )
-  }
-
   log_r_site <- construct_spatial_rt(
     log_state_rt = log_r_state_week,
     spatial_deviation_ar_coeff = phi_rt,
@@ -401,11 +388,12 @@ generate_simulated_data <- function(r_in_weeks = # nolint
       mean = 0,
       sd = 1
     )
+    stdev_generalized <- sqrt(sigma_sqrd_generalized^(1 / n_sites))
     log_r_site_aux <- construct_aux_rt(
       log_state_rt = log_r_state_week,
       state_deviation_ar_coeff = phi_rt,
       scaling_factor = scaling_factor,
-      sigma_eps = sigma_generalized^(1 / n_sites),
+      sigma_eps = stdev_generalized,
       z = state_deviation_noise_vec,
       init_stat = init_stat
     )
@@ -598,7 +586,8 @@ generate_simulated_data <- function(r_in_weeks = # nolint
     hosp_data = hosp_data,
     hosp_data_eval = hosp_data_eval,
     rt_site_data = r_site,
-    rt_global_data = rt
+    rt_global_data = rt,
+    corr_matrix = corr_matrix
   )
 
   return(example_data)

R/get_stan_data.R

@@ -375,7 +375,7 @@ get_stan_data <- function(input_count_data,
   # If user does / doesn't want spatial comps.
   # We can add an extra step here for when spatial desired and dist_matrix
   #   not given.
-  if (corr_structure_switch == 0) {
+  if (corr_structure_switch == 0 || is.null(dist_matrix)) {
     # This dist_matrix will not be used, only needed for stan data specs.
     dist_matrix <- matrix(
       0,

R/independence_corr_func.R → R/independence_corr_func_r.R

@@ -5,7 +5,7 @@
 #'
 #' @return Correlation matrix of diagonal ones.
 #'
-independence_corr_func <- function(
+independence_corr_func_r <- function(
     corr_function_params = list(
       num_sites = NULL
     )) {
@@ -21,5 +21,5 @@ independence_corr_func <- function(
       is.numeric(corr_function_params$num_sites)
   )
 
-  return(diag(x = 1, nrow = corr_function_params$num_sites))
+  return(diag(x = 1, nrow = as.integer(corr_function_params$num_sites)))
 }

R/initialization.R

@@ -94,7 +94,28 @@ get_inits_for_one_chain <- function(stan_data, stdev = 0.01) {
     hosp_wday_effect = to_simplex(abs(
       stats::rnorm(7, 1 / 7, stdev)
     )),
-    infection_feedback = abs(stats::rnorm(1, 500, 20))
+    infection_feedback = abs(stats::rnorm(1, 500, 20)),
+    # Spatial inits
+    log_sigma_generalized = stats::rnorm(1, log(0.05^(n_subpops - 1)), 0.5),
+    log_phi = stats::rnorm(1, log(0.25), 0.1),
+    log_scaling_factor = stats::rnorm(1, log(1), 0.1),
+    non_cent_spatial_dev_ns_mat = matrix(
+      stats::rnorm((n_subpops - 1) * n_weeks,
+        mean = 0,
+        sd = stdev
+      ),
+      (n_subpops - 1),
+      n_weeks
+    ),
+    norm_vec_aux_site = stats::rnorm(n_weeks, 0, stdev),
+    # Initialize the cholesky decomposition matrix if inferring
+    # unstructured correlation matrix
+    L_Omega = as.matrix(diag(2))
   )
+
+  if (stan_data$corr_structure_switch == 2) {
+    init_list$L_Omega <- diag((n_subpops - 1))
+  }
+
   return(init_list)
 }

R/rand_corr_matrix_func.R

@@ -0,0 +1,49 @@
+#' Generates a random correlation matrix.
+#'
+#' @param corr_function_params Structured as follows : {n_sites, seed}.
+#' Ensure that this list follows this structure. Number of sites n_sites, and
+#' seed desired for rng.  Defaults are set to NULL.
+#' @export
+#'
+#' @return Correlation matrix.  A symmetric matrix with off diagonal values
+#' sampled randomly from -1 to 1, using a uniform distribution.
+#'
+rand_corr_matrix_func <- function(
+    corr_function_params = list(
+      n_sites = NULL,
+      seed = NULL
+    )) {
+  # error managing
+  stopifnot(
+    "corr_function_params : NOT STRUCTURED CORRECTLY!!!\n
+             List names should be : {n_sites, seed}" =
+      names(corr_function_params) == c("n_sites", "seed")
+  )
+  # presets
+  n_sites <- corr_function_params$n_sites
+  seed <- corr_function_params$seed
+  set.seed(seed)
+  random_matrix <- matrix(
+    data = runif(
+      n_sites * n_sites,
+      min = -1, max = 1
+    ),
+    nrow = n_sites
+  )
+  sym_matrix <- (random_matrix + t(random_matrix)) / 2
+  eigen_decomp <- eigen(
+    sym_matrix
+  )
+  positive_semi_definite_matrix <- eigen_decomp$vectors %*%
+    diag(
+      pmax(
+        eigen_decomp$values,
+        0
+      )
+    ) %*%
+    t(eigen_decomp$vectors)
+  diag(positive_semi_definite_matrix) <- 1
+
+
+  return(round(positive_semi_definite_matrix, 4))
+}

_pkgdown.yml

@@ -1,3 +1,4 @@
 url: https://cdcgov.github.io/ww-inference-model/
 template:
   bootstrap: 5
+  math-rendering: katex

inst/extdata/example_params.toml

@@ -3,13 +3,13 @@ uot = 50
 
 [infection_process]
 # spatial params
-log_phi_mu_prior = -1.609438 # log(0.2)
-log_phi_sd_prior = 0.2
+log_phi_mu_prior = -1.386294 # log(0.25)
+log_phi_sd_prior = 0.75
 l = 1
-log_sigma_generalized_mu_prior = -11.98293 # log(0.05^4)
-log_sigma_generalized_sd_prior = 0.2
-log_scaling_factor_mu_prior = 0.09531018 # log(1.1)
-log_scaling_factor_sd_prior = 0.15
+log_sigma_generalized_mu_prior = -5.298317 # log(0.05^4)
+log_sigma_generalized_sd_prior = 0.5
+log_scaling_factor_mu_prior = 0 # log(1.1)
+log_scaling_factor_sd_prior = 0.5
 
 r_prior_mean = 1
 r_prior_sd = 1

inst/stan/wwinference.stan

@@ -135,7 +135,6 @@ parameters {
   real<lower=0> sigma_rt; // magnitude of site level variation from state level
   real<lower=0, upper=1> autoreg_rt_site;
   real<lower=0, upper=1> autoreg_p_hosp;
-  matrix[n_subpops, n_weeks] error_site; // matrix of subpopulations
   real<lower=0,upper=1> i_first_obs_over_n; // per capita
   // infection incidence on the day of the first observed infection
   vector[n_subpops] eta_i_first_obs; // z-score on logit scale of site
@@ -170,7 +169,7 @@ parameters {
   real log_scaling_factor;
   matrix[n_subpops-1,n_weeks] non_cent_spatial_dev_ns_mat;
   vector[n_weeks] norm_vec_aux_site;
-  cholesky_factor_corr[corr_structure_switch == 2 ? n_subpops-1 : 0] L_Omega;
+  cholesky_factor_corr[corr_structure_switch == 2 ? n_subpops-1 : 2] L_Omega;
   //----------------------------------------------------------------------------
 }
 //
@@ -259,7 +258,7 @@ transformed parameters {
     log_r_mu_t_in_weeks,
     autoreg_rt_site,
     scaling_factor,
-    sigma_generalized,
+    sqrt(sigma_generalized),
     norm_vec_aux_site,
     0
   );