From 5001aaee8a9797ecae7cafbbb7ebbcf04d7fd3b7 Mon Sep 17 00:00:00 2001 From: cbernalz Date: Mon, 15 Jul 2024 14:18:42 -0700 Subject: [PATCH 1/7] 2024-07-15 update : added spatial unadj rt proc. func. --- R/spatial_unadj_rt_process.R | 45 ++++++++++++++++++++++++++++++++++++ 1 file changed, 45 insertions(+) create mode 100644 R/spatial_unadj_rt_process.R diff --git a/R/spatial_unadj_rt_process.R b/R/spatial_unadj_rt_process.R new file mode 100644 index 0000000..bcaea30 --- /dev/null +++ b/R/spatial_unadj_rt_process.R @@ -0,0 +1,45 @@ +#' Generates the subpopulation unadjusted Rt trajectories from proposal 2. +#' +#' @param state_rt "State" level unadjusted Rt, on log scale. +#' @param corr_function Function for the correlation structure desired. +#' @param corr_function_params List of parameters of desired correlation function, should contain dist matrix. +#' @param phi_rt Coefficient for AR(1) temporal correlation on subpopulation deviations. +#' @param sigma_eps Parameter for construction of covariance matrix of spatial epsilon. +#' +#' @return A matrix for unadjusted Rt where rows are subpopulations and columns are time points. +#' +spatial_unadj_rt_process = function( log_state_Rt, corr_function, corr_function_params, phi_rt, sigma_eps, pre_seed = 1 ) { + + set.seed(pre_seed) + + # correlation matrix constr. + omega_matrix_eps = corr_function( corr_function_params ) + sigma_matrix_eps = sigma_eps^2 * omega_matrix_eps + + # presets + n_subpopulations = nrow( sigma_matrix_eps ) + n_time = length( log_state_Rt ) + log_site_rt = matrix( data = 0, nrow = n_subpopulations, ncol = n_time ) + delta = matrix( data = 0, nrow = n_subpopulations, ncol = n_time ) + + # delta constr. + delta[,1] = MASS::mvrnorm( n = 1, + mu = matrix( data = 0, nrow = 1, ncol = n_subpopulations ), + Sigma = sigma_matrix_eps ) + for (t_i in 2:n_time) { + eps_vec = MASS::mvrnorm( n = 1, + mu = matrix( data = 0, nrow = 1, ncol = n_subpopulations ), + Sigma = sigma_matrix_eps ) + delta[,t_i] = phi_Rt * delta[,t_i - 1] + eps_vec + } + + # Subpopulation unadjusted Rt constr. + for (t_i in 1:n_time) { + log_site_rt[,t_i] = (log_state_rt[t_i] * matrix( data = 1, + nrow = n_subpopulations, + ncol = 1 )) + + delta[,t_i] + } + + return( log_site_rt ) +} \ No newline at end of file From 13eb48f3d662c854b34757c8bc7c71fc113e780c Mon Sep 17 00:00:00 2001 From: cbernalz Date: Mon, 15 Jul 2024 14:36:41 -0700 Subject: [PATCH 2/7] 2024-07-15 update : added corr funcs. and changed gen. sim. --- R/exponential_decay_corr_func.R | 36 +++++++++++++++++++++++ R/generate_simulated_data.R | 51 ++++++++++++++++++++++++++------- R/independence_corr_func.R | 18 ++++++++++++ 3 files changed, 94 insertions(+), 11 deletions(-) create mode 100644 R/exponential_decay_corr_func.R create mode 100644 R/independence_corr_func.R diff --git a/R/exponential_decay_corr_func.R b/R/exponential_decay_corr_func.R new file mode 100644 index 0000000..421907a --- /dev/null +++ b/R/exponential_decay_corr_func.R @@ -0,0 +1,36 @@ +#' Generates the correlation matrix for an exponential decay correlation function. +#' +#' @param corr_function_params Sturcuted as follows : {dist.matrix, phi, l}. +#' +#' @return Correlation matrix. +#' +exponential_decay_corr_func = function( corr_function_params = list( dist_matrix = NULL, phi = NULL, l = NULL ) ) { + + # error managing + stopifnot( 'corr_function_params : NOT STRUCTURED CORRECTLY!!!\n + List names should be : {dist_matrix, phi, l}' = + names( corr_function_params ) == c( 'dist_matrix', 'phi', 'l' ) ) + stopifnot( 'corr_function_params : DOES NOT CONTAIN VALID VALUES!!!\n + List should be {matrix, int, int}' = + is.matrix( corr_function_params$dist_matrix ) & + is.numeric( corr_function_params$phi ) & + is.numeric( corr_function_params$l ) ) + stopifnot( 'corr_function_params : DIST.MATRIX NOT SQUARE!!!\n + Make sure dist_matrix is square' = + nrow( corr_function_params$dist_matrix ) == ncol( corr_function_params$dist_matrix ) ) + + # presets + dist_matrix = corr_function_params$dist_matrix + phi = corr_function_params$phi + l = corr_function_params$l + corr_matrix = matrix( data = 0, nrow = nrow( dist_matrix ), ncol = ncol( dist_matrix ) ) + + for (i in 1:nrow( dist_matrix ) ) { + for (j in 1:ncol( dist_matrix ) ) { + d_ij = dist.matrix[i,j] + corr_matrix[i,j] = exp( -( d_ij / phi )^l ) + } + } + + return( corr_matrix ) +} \ No newline at end of file diff --git a/R/generate_simulated_data.R b/R/generate_simulated_data.R index 43c3fd4..f6e304d 100644 --- a/R/generate_simulated_data.R +++ b/R/generate_simulated_data.R @@ -50,6 +50,12 @@ #' LOD across sites #' @param input_params_path path to the toml file with the parameters to use #' to generate the simulated data +#' @param corr_function Correlation function for spatial site correlations +#' @param corr_fun_params Parameters required for correlation function +#' @param phi_rt Coefficient for AR(1) temporal correlation on subpopulation deviations. +#' @param sigma_eps Parameter for construction of covariance matrix of spatial epsilon. +#' @param scaling_factor Scaling factor for aux site. +#' @param aux_site_bool Boolean to use the aux site framework with scaling factor. #' #' @return a list containing three dataframes. hosp_data is a dataframe #' containing the number of daily hospital admissions by day for a theoretical @@ -108,7 +114,13 @@ generate_simulated_data <- function(r_in_weeks = # nolint fs::path_package("extdata", "example_params.toml", package = "wwinference" - )) { + ), + corr_function = independence_corr_func, + corr_fun_params = NULL, + phi_rt = 0.6, + sigma_eps = sqrt(0.02), + scaling_factor = NULL, + aux_site_bool = FALSE) { # Some quick checks to make sure the inputs work as expected stopifnot( "weekly R(t) passed in isn't long enough" = @@ -122,6 +134,9 @@ generate_simulated_data <- function(r_in_weeks = # nolint "Site and lab indices don't align" = length(site) == length(lab) ) + + # presetting corr_fun_params + if(is.null(corr_fun_params)){ corr_fun_params = list(num_sites = n_sites)} # Get pop fractions of each subpop. There will n_sites + 1 subpops @@ -266,19 +281,9 @@ generate_simulated_data <- function(r_in_weeks = # nolint r_site <- matrix(nrow = n_sites + 1, ncol = (ot + ht)) # Generate site-level R(t) log_r_state_week <- log(unadj_r_weeks) - log_r_site <- matrix(nrow = n_sites + 1, ncol = n_weeks) initial_growth_site <- vector(length = n_sites + 1) log_i0_over_n_site <- vector(length = n_sites + 1) for (i in 1:(n_sites + 1)) { - # This creates each R(t) vector for each subpopulation, by sampling - # from a normal distribution centered on the state R(t). - # In the model, this is an AR(1) process based on the previous deviation - log_r_site[i, ] <- rnorm( - n = n_weeks, - mean = log_r_state_week, - sd = 0.05 - ) # sigma_rt - # Generate deviations in the initial growth rate and initial incidence initial_growth_site[i] <- rnorm( n = 1, mean = initial_growth, @@ -290,6 +295,30 @@ generate_simulated_data <- function(r_in_weeks = # nolint sd = 0.5 ) } + log_r_site <- spatial_rt_process(log_r_state_week, corr_function, corr_fun_params, phi_rt, sigma_eps) + # Auxiliary Site + if(aux_site_bool){ + log_r_site_aux <- matrix(data = 0, + nrow = 1, + ncol = ncol(log_r_site)) + delta <- matrix(data = 0, + nrow = 1, + ncol = ncol(log_r_site)) + delta[1] <- rnorm(n = 1, + mean = 0, + sd = sqrt(scaling_factor) * sigma_eps) + for (i in 2:ncol(log_r_site)) { + eps_temp <- rnorm(n = 1, + mean = 0, + sd = sqrt(scaling_factor) * sigma_eps) + delta[i] <- phi_rt * delta[i - 1] + eps_temp + } + for (i in 1:ncol(log_r_site)) { + log_r_site_aux[i] <- log_r_state_week[i] + delta[i] + } + log_r_site = rbind(log_r_site, + log_r_site_aux) + } new_i_over_n <- rep(0, (uot + ot + ht)) # State infections diff --git a/R/independence_corr_func.R b/R/independence_corr_func.R new file mode 100644 index 0000000..c8e5dc3 --- /dev/null +++ b/R/independence_corr_func.R @@ -0,0 +1,18 @@ +#' Generates the correlation matrix for a non-spatial process. +#' +#' @param corr_function_params NULL +#' +#' @return Correlation matrix of diagonal ones. +#' +independence_corr_func = function( corr_function_params = list( num_sites = NULL ) ) { + + # error managing + stopifnot( 'corr.function.params : NOT STRUCTURED CORRECTLY!!! + List names should be : {num_sites}' = + names( corr_function_params ) == c( 'num_sites' ) ) + stopifnot( 'corr.function.params : DOES NOT CONTAIN VALID VALUES!!! + List should be {int}' = + is.numeric( corr_function_params$num_sites ) ) + + return( diag(x = 1, nrow = corr_function_params$num_sites) ) +} \ No newline at end of file From b084cf244fc2809c33033e0ad3a25856a5aa4603 Mon Sep 17 00:00:00 2001 From: cbernalz Date: Mon, 15 Jul 2024 14:45:00 -0700 Subject: [PATCH 3/7] 2024-07-15 update : removing seeding. --- R/spatial_unadj_rt_process.R | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/R/spatial_unadj_rt_process.R b/R/spatial_unadj_rt_process.R index bcaea30..a356be3 100644 --- a/R/spatial_unadj_rt_process.R +++ b/R/spatial_unadj_rt_process.R @@ -8,9 +8,7 @@ #' #' @return A matrix for unadjusted Rt where rows are subpopulations and columns are time points. #' -spatial_unadj_rt_process = function( log_state_Rt, corr_function, corr_function_params, phi_rt, sigma_eps, pre_seed = 1 ) { - - set.seed(pre_seed) +spatial_unadj_rt_process = function(log_state_Rt, corr_function, corr_function_params, phi_rt, sigma_eps) { # correlation matrix constr. omega_matrix_eps = corr_function( corr_function_params ) From 786b4b927d5a4246a064c5562998ae972ee28e5a Mon Sep 17 00:00:00 2001 From: cbernalz Date: Mon, 15 Jul 2024 15:18:19 -0700 Subject: [PATCH 4/7] 2024-07-15 update : fixing to snake_script. --- R/exponential_decay_corr_func.R | 56 ++++++++++++++++++--------------- R/generate_simulated_data.R | 44 ++++++++++++++++---------- R/independence_corr_func.R | 23 ++++++++------ R/spatial_unadj_rt_process.R | 45 ++++++++++++++------------ 4 files changed, 97 insertions(+), 71 deletions(-) diff --git a/R/exponential_decay_corr_func.R b/R/exponential_decay_corr_func.R index 421907a..257806f 100644 --- a/R/exponential_decay_corr_func.R +++ b/R/exponential_decay_corr_func.R @@ -4,33 +4,39 @@ #' #' @return Correlation matrix. #' -exponential_decay_corr_func = function( corr_function_params = list( dist_matrix = NULL, phi = NULL, l = NULL ) ) { +exponential_decay_corr_func <- function(corr_function_params = list(dist_matrix = NULL, phi = NULL, l = NULL)) { # error managing - stopifnot( 'corr_function_params : NOT STRUCTURED CORRECTLY!!!\n - List names should be : {dist_matrix, phi, l}' = - names( corr_function_params ) == c( 'dist_matrix', 'phi', 'l' ) ) - stopifnot( 'corr_function_params : DOES NOT CONTAIN VALID VALUES!!!\n - List should be {matrix, int, int}' = - is.matrix( corr_function_params$dist_matrix ) & - is.numeric( corr_function_params$phi ) & - is.numeric( corr_function_params$l ) ) - stopifnot( 'corr_function_params : DIST.MATRIX NOT SQUARE!!!\n - Make sure dist_matrix is square' = - nrow( corr_function_params$dist_matrix ) == ncol( corr_function_params$dist_matrix ) ) - + stopifnot( + "corr_function_params : NOT STRUCTURED CORRECTLY!!!\n + List names should be : {dist_matrix, phi, l}" = + names(corr_function_params) == c("dist_matrix", "phi", "l") + ) + stopifnot( + "corr_function_params : DOES NOT CONTAIN VALID VALUES!!!\n + List should be {matrix, int, int}" = + is.matrix(corr_function_params$dist_matrix) & + is.numeric(corr_function_params$phi) & + is.numeric(corr_function_params$l) + ) + stopifnot( + "corr_function_params : DIST.MATRIX NOT SQUARE!!!\n + Make sure dist_matrix is square" = + nrow(corr_function_params$dist_matrix) == ncol(corr_function_params$dist_matrix) + ) + # presets - dist_matrix = corr_function_params$dist_matrix - phi = corr_function_params$phi - l = corr_function_params$l - corr_matrix = matrix( data = 0, nrow = nrow( dist_matrix ), ncol = ncol( dist_matrix ) ) - - for (i in 1:nrow( dist_matrix ) ) { - for (j in 1:ncol( dist_matrix ) ) { - d_ij = dist.matrix[i,j] - corr_matrix[i,j] = exp( -( d_ij / phi )^l ) - } - } + dist_matrix <- corr_function_params$dist_matrix + phi <- corr_function_params$phi + l <- corr_function_params$l + corr_matrix <- matrix(data = 0, nrow = nrow(dist_matrix), ncol = ncol(dist_matrix)) + + for (i in 1:nrow(dist_matrix)) { + for (j in 1:ncol(dist_matrix)) { + d_ij <- dist_matrix[i, j] + corr_matrix[i, j] <- exp(-(d_ij / phi)^l) + } + } - return( corr_matrix ) + return(corr_matrix) } \ No newline at end of file diff --git a/R/generate_simulated_data.R b/R/generate_simulated_data.R index f6e304d..3b282c2 100644 --- a/R/generate_simulated_data.R +++ b/R/generate_simulated_data.R @@ -136,7 +136,9 @@ generate_simulated_data <- function(r_in_weeks = # nolint ) # presetting corr_fun_params - if(is.null(corr_fun_params)){ corr_fun_params = list(num_sites = n_sites)} + if(is.null(corr_fun_params)){ + corr_fun_params = list(num_sites = n_sites) + } # Get pop fractions of each subpop. There will n_sites + 1 subpops @@ -297,27 +299,37 @@ generate_simulated_data <- function(r_in_weeks = # nolint } log_r_site <- spatial_rt_process(log_r_state_week, corr_function, corr_fun_params, phi_rt, sigma_eps) # Auxiliary Site - if(aux_site_bool){ - log_r_site_aux <- matrix(data = 0, - nrow = 1, - ncol = ncol(log_r_site)) - delta <- matrix(data = 0, - nrow = 1, - ncol = ncol(log_r_site)) - delta[1] <- rnorm(n = 1, - mean = 0, - sd = sqrt(scaling_factor) * sigma_eps) + if (aux_site_bool) { + log_r_site_aux <- matrix( + data = 0, + nrow = 1, + ncol = ncol(log_r_site) + ) + delta <- matrix( + data = 0, + nrow = 1, + ncol = ncol(log_r_site) + ) + delta[1] <- rnorm( + n = 1, + mean = 0, + sd = sqrt(scaling_factor) * sigma_eps + ) for (i in 2:ncol(log_r_site)) { - eps_temp <- rnorm(n = 1, - mean = 0, - sd = sqrt(scaling_factor) * sigma_eps) + eps_temp <- rnorm( + n = 1, + mean = 0, + sd = sqrt(scaling_factor) * sigma_eps + ) delta[i] <- phi_rt * delta[i - 1] + eps_temp } for (i in 1:ncol(log_r_site)) { log_r_site_aux[i] <- log_r_state_week[i] + delta[i] } - log_r_site = rbind(log_r_site, - log_r_site_aux) + log_r_site <- rbind( + log_r_site, + log_r_site_aux + ) } diff --git a/R/independence_corr_func.R b/R/independence_corr_func.R index c8e5dc3..ae3fe12 100644 --- a/R/independence_corr_func.R +++ b/R/independence_corr_func.R @@ -4,15 +4,18 @@ #' #' @return Correlation matrix of diagonal ones. #' -independence_corr_func = function( corr_function_params = list( num_sites = NULL ) ) { - +independence_corr_func <- function(corr_function_params = list(num_sites = NULL)) { # error managing - stopifnot( 'corr.function.params : NOT STRUCTURED CORRECTLY!!! - List names should be : {num_sites}' = - names( corr_function_params ) == c( 'num_sites' ) ) - stopifnot( 'corr.function.params : DOES NOT CONTAIN VALID VALUES!!! - List should be {int}' = - is.numeric( corr_function_params$num_sites ) ) - - return( diag(x = 1, nrow = corr_function_params$num_sites) ) + stopifnot( + "corr.function.params : NOT STRUCTURED CORRECTLY!!! + List names should be : {num_sites}" = + names(corr_function_params) == c("num_sites") + ) + stopifnot( + "corr.function.params : DOES NOT CONTAIN VALID VALUES!!! + List should be {int}" = + is.numeric(corr_function_params$num_sites) + ) + + return(diag(x = 1, nrow = corr_function_params$num_sites)) } \ No newline at end of file diff --git a/R/spatial_unadj_rt_process.R b/R/spatial_unadj_rt_process.R index a356be3..f4c2b1c 100644 --- a/R/spatial_unadj_rt_process.R +++ b/R/spatial_unadj_rt_process.R @@ -8,36 +8,41 @@ #' #' @return A matrix for unadjusted Rt where rows are subpopulations and columns are time points. #' -spatial_unadj_rt_process = function(log_state_Rt, corr_function, corr_function_params, phi_rt, sigma_eps) { - +spatial_unadj_rt_process <- function(log_state_Rt, corr_function, corr_function_params, phi_rt, sigma_eps) { # correlation matrix constr. - omega_matrix_eps = corr_function( corr_function_params ) - sigma_matrix_eps = sigma_eps^2 * omega_matrix_eps + omega_matrix_eps <- corr_function(corr_function_params) + sigma_matrix_eps <- sigma_eps^2 * omega_matrix_eps # presets - n_subpopulations = nrow( sigma_matrix_eps ) - n_time = length( log_state_Rt ) - log_site_rt = matrix( data = 0, nrow = n_subpopulations, ncol = n_time ) - delta = matrix( data = 0, nrow = n_subpopulations, ncol = n_time ) + n_subpopulations <- nrow(sigma_matrix_eps) + n_time <- length(log_state_Rt) + log_site_rt <- matrix(data = 0, nrow = n_subpopulations, ncol = n_time) + delta <- matrix(data = 0, nrow = n_subpopulations, ncol = n_time) # delta constr. - delta[,1] = MASS::mvrnorm( n = 1, - mu = matrix( data = 0, nrow = 1, ncol = n_subpopulations ), - Sigma = sigma_matrix_eps ) + delta[, 1] <- MASS::mvrnorm( + n = 1, + mu = matrix(data = 0, nrow = 1, ncol = n_subpopulations), + Sigma = sigma_matrix_eps + ) for (t_i in 2:n_time) { - eps_vec = MASS::mvrnorm( n = 1, - mu = matrix( data = 0, nrow = 1, ncol = n_subpopulations ), - Sigma = sigma_matrix_eps ) - delta[,t_i] = phi_Rt * delta[,t_i - 1] + eps_vec + eps_vec <- MASS::mvrnorm( + n = 1, + mu = matrix(data = 0, nrow = 1, ncol = n_subpopulations), + Sigma = sigma_matrix_eps + ) + delta[, t_i] <- phi_Rt * delta[, t_i - 1] + eps_vec } # Subpopulation unadjusted Rt constr. for (t_i in 1:n_time) { - log_site_rt[,t_i] = (log_state_rt[t_i] * matrix( data = 1, - nrow = n_subpopulations, - ncol = 1 )) + - delta[,t_i] + log_site_rt[, t_i] <- (log_state_rt[t_i] * matrix( + data = 1, + nrow = n_subpopulations, + ncol = 1 + )) + + delta[, t_i] } - return( log_site_rt ) + return(log_site_rt) } \ No newline at end of file From d74add56945b32a03586ab141069c0617da66863 Mon Sep 17 00:00:00 2001 From: cbernalz Date: Tue, 16 Jul 2024 08:11:38 -0700 Subject: [PATCH 5/7] 2025-07-16 update : fixing style errors. --- R/exponential_decay_corr_func.R | 46 +++++++++++++++++++-------------- R/generate_simulated_data.R | 28 +++++++++++++------- R/independence_corr_func.R | 9 ++++--- R/spatial_unadj_rt_process.R | 32 ++++++++++++++--------- 4 files changed, 72 insertions(+), 43 deletions(-) diff --git a/R/exponential_decay_corr_func.R b/R/exponential_decay_corr_func.R index 257806f..64b7e76 100644 --- a/R/exponential_decay_corr_func.R +++ b/R/exponential_decay_corr_func.R @@ -1,42 +1,50 @@ -#' Generates the correlation matrix for an exponential decay correlation function. +#' Generates the correlation matrix for an exponential +#' decay correlation function. #' #' @param corr_function_params Sturcuted as follows : {dist.matrix, phi, l}. #' #' @return Correlation matrix. -#' -exponential_decay_corr_func <- function(corr_function_params = list(dist_matrix = NULL, phi = NULL, l = NULL)) { - +#' +exponential_decay_corr_func <- function( + corr_function_params = list( + dist_matrix = NULL, + phi = NULL, + l = NULL + )) { # error managing - stopifnot( + stopifnot( "corr_function_params : NOT STRUCTURED CORRECTLY!!!\n List names should be : {dist_matrix, phi, l}" = names(corr_function_params) == c("dist_matrix", "phi", "l") ) + # presets + dist_matrix <- corr_function_params$dist_matrix + phi <- corr_function_params$phi + l <- corr_function_params$l stopifnot( "corr_function_params : DOES NOT CONTAIN VALID VALUES!!!\n List should be {matrix, int, int}" = - is.matrix(corr_function_params$dist_matrix) & - is.numeric(corr_function_params$phi) & - is.numeric(corr_function_params$l) + is.matrix(dist_matrix) & is.numeric(phi) & is.numeric(l) ) stopifnot( "corr_function_params : DIST.MATRIX NOT SQUARE!!!\n Make sure dist_matrix is square" = - nrow(corr_function_params$dist_matrix) == ncol(corr_function_params$dist_matrix) + nrow(dist_matrix) == ncol(dist_matrix) ) # presets - dist_matrix <- corr_function_params$dist_matrix - phi <- corr_function_params$phi - l <- corr_function_params$l - corr_matrix <- matrix(data = 0, nrow = nrow(dist_matrix), ncol = ncol(dist_matrix)) + corr_matrix <- matrix( + data = 0, + nrow = nrow(dist_matrix), + ncol = ncol(dist_matrix) + ) - for (i in 1:nrow(dist_matrix)) { - for (j in 1:ncol(dist_matrix)) { + for (i in seq_len(nrow(dist_matrix))) { + for (j in seq_len(ncol(dist_matrix))) { d_ij <- dist_matrix[i, j] corr_matrix[i, j] <- exp(-(d_ij / phi)^l) - } - } - + } + } + return(corr_matrix) -} \ No newline at end of file +} diff --git a/R/generate_simulated_data.R b/R/generate_simulated_data.R index 3b282c2..9ef31c6 100644 --- a/R/generate_simulated_data.R +++ b/R/generate_simulated_data.R @@ -52,10 +52,13 @@ #' to generate the simulated data #' @param corr_function Correlation function for spatial site correlations #' @param corr_fun_params Parameters required for correlation function -#' @param phi_rt Coefficient for AR(1) temporal correlation on subpopulation deviations. -#' @param sigma_eps Parameter for construction of covariance matrix of spatial epsilon. +#' @param phi_rt Coefficient for AR(1) temporal correlation on subpopulation +#' deviations. +#' @param sigma_eps Parameter for construction of covariance matrix of spatial +#' epsilon. #' @param scaling_factor Scaling factor for aux site. -#' @param aux_site_bool Boolean to use the aux site framework with scaling factor. +#' @param aux_site_bool Boolean to use the aux site framework with +#' scaling factor. #' #' @return a list containing three dataframes. hosp_data is a dataframe #' containing the number of daily hospital admissions by day for a theoretical @@ -134,10 +137,10 @@ generate_simulated_data <- function(r_in_weeks = # nolint "Site and lab indices don't align" = length(site) == length(lab) ) - + # presetting corr_fun_params - if(is.null(corr_fun_params)){ - corr_fun_params = list(num_sites = n_sites) + if (is.null(corr_fun_params)) { + corr_fun_params <- list(num_sites = n_sites) } @@ -297,7 +300,13 @@ generate_simulated_data <- function(r_in_weeks = # nolint sd = 0.5 ) } - log_r_site <- spatial_rt_process(log_r_state_week, corr_function, corr_fun_params, phi_rt, sigma_eps) + log_r_site <- spatial_rt_process( + log_r_state_week, + corr_function, + corr_fun_params, + phi_rt, + sigma_eps + ) # Auxiliary Site if (aux_site_bool) { log_r_site_aux <- matrix( @@ -315,7 +324,8 @@ generate_simulated_data <- function(r_in_weeks = # nolint mean = 0, sd = sqrt(scaling_factor) * sigma_eps ) - for (i in 2:ncol(log_r_site)) { + n_t + for (i in 2:n_t) { eps_temp <- rnorm( n = 1, mean = 0, @@ -323,7 +333,7 @@ generate_simulated_data <- function(r_in_weeks = # nolint ) delta[i] <- phi_rt * delta[i - 1] + eps_temp } - for (i in 1:ncol(log_r_site)) { + for (i in 1:n_t) { log_r_site_aux[i] <- log_r_state_week[i] + delta[i] } log_r_site <- rbind( diff --git a/R/independence_corr_func.R b/R/independence_corr_func.R index ae3fe12..a93df89 100644 --- a/R/independence_corr_func.R +++ b/R/independence_corr_func.R @@ -3,8 +3,11 @@ #' @param corr_function_params NULL #' #' @return Correlation matrix of diagonal ones. -#' -independence_corr_func <- function(corr_function_params = list(num_sites = NULL)) { +#' +independence_corr_func <- function( + corr_function_params = list( + num_sites = NULL + )) { # error managing stopifnot( "corr.function.params : NOT STRUCTURED CORRECTLY!!! @@ -18,4 +21,4 @@ independence_corr_func <- function(corr_function_params = list(num_sites = NULL) ) return(diag(x = 1, nrow = corr_function_params$num_sites)) -} \ No newline at end of file +} diff --git a/R/spatial_unadj_rt_process.R b/R/spatial_unadj_rt_process.R index f4c2b1c..b3eb354 100644 --- a/R/spatial_unadj_rt_process.R +++ b/R/spatial_unadj_rt_process.R @@ -2,23 +2,31 @@ #' #' @param state_rt "State" level unadjusted Rt, on log scale. #' @param corr_function Function for the correlation structure desired. -#' @param corr_function_params List of parameters of desired correlation function, should contain dist matrix. -#' @param phi_rt Coefficient for AR(1) temporal correlation on subpopulation deviations. -#' @param sigma_eps Parameter for construction of covariance matrix of spatial epsilon. +#' @param corr_function_params List of parameters of desired correlation +#' function, should contain dist matrix. +#' @param phi_rt Coefficient for AR(1) temporal correlation on +#' subpopulation deviations. +#' @param sigma_eps Parameter for construction of covariance matrix +#' of spatial epsilon. #' -#' @return A matrix for unadjusted Rt where rows are subpopulations and columns are time points. -#' -spatial_unadj_rt_process <- function(log_state_Rt, corr_function, corr_function_params, phi_rt, sigma_eps) { +#' @return A matrix for unadjusted Rt where rows are subpopulations +#' and columns are time points. +#' +spatial_unadj_rt_process <- function(log_state_rt, + corr_function, + corr_function_params, + phi_rt, + sigma_eps) { # correlation matrix constr. omega_matrix_eps <- corr_function(corr_function_params) sigma_matrix_eps <- sigma_eps^2 * omega_matrix_eps - + # presets n_subpopulations <- nrow(sigma_matrix_eps) - n_time <- length(log_state_Rt) + n_time <- length(log_state_rt) log_site_rt <- matrix(data = 0, nrow = n_subpopulations, ncol = n_time) delta <- matrix(data = 0, nrow = n_subpopulations, ncol = n_time) - + # delta constr. delta[, 1] <- MASS::mvrnorm( n = 1, @@ -33,7 +41,7 @@ spatial_unadj_rt_process <- function(log_state_Rt, corr_function, corr_function_ ) delta[, t_i] <- phi_Rt * delta[, t_i - 1] + eps_vec } - + # Subpopulation unadjusted Rt constr. for (t_i in 1:n_time) { log_site_rt[, t_i] <- (log_state_rt[t_i] * matrix( @@ -43,6 +51,6 @@ spatial_unadj_rt_process <- function(log_state_Rt, corr_function, corr_function_ )) + delta[, t_i] } - + return(log_site_rt) -} \ No newline at end of file +} From c446eb3e9dc805ed207a193bd914edc619c11e80 Mon Sep 17 00:00:00 2001 From: cbernalz Date: Tue, 16 Jul 2024 08:55:17 -0700 Subject: [PATCH 6/7] 2025-07-16 update : updating DESCRIPTION. --- DESCRIPTION | 2 +- NAMESPACE | 1 + R/spatial_unadj_rt_process.R | 4 ++-- R/wwinference-package.R | 1 + man/exponential_decay_corr_func.Rd | 21 ++++++++++++++++++ man/generate_simulated_data.Rd | 23 +++++++++++++++++++- man/independence_corr_func.Rd | 17 +++++++++++++++ man/spatial_unadj_rt_process.Rd | 35 ++++++++++++++++++++++++++++++ man/wwinference-package.Rd | 3 ++- 9 files changed, 102 insertions(+), 5 deletions(-) create mode 100644 man/exponential_decay_corr_func.Rd create mode 100644 man/independence_corr_func.Rd create mode 100644 man/spatial_unadj_rt_process.Rd diff --git a/DESCRIPTION b/DESCRIPTION index 0233493..4b8746b 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -62,7 +62,7 @@ Depends: SystemRequirements: CmdStan (>=2.35.0) Encoding: UTF-8 Roxygen: list(markdown = TRUE) -RoxygenNote: 7.3.1 +RoxygenNote: 7.3.2 Suggests: testthat (>= 3.0.0), bookdown, diff --git a/NAMESPACE b/NAMESPACE index 310bc70..c3b5bfa 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -32,6 +32,7 @@ export(preprocess_ww_data) export(to_simplex) export(validate_paramlist) export(wwinference) +importFrom(MASS,mvrnorm) importFrom(cmdstanr,cmdstan_model) importFrom(dplyr,arrange) importFrom(dplyr,as_tibble) diff --git a/R/spatial_unadj_rt_process.R b/R/spatial_unadj_rt_process.R index b3eb354..e55de2a 100644 --- a/R/spatial_unadj_rt_process.R +++ b/R/spatial_unadj_rt_process.R @@ -28,13 +28,13 @@ spatial_unadj_rt_process <- function(log_state_rt, delta <- matrix(data = 0, nrow = n_subpopulations, ncol = n_time) # delta constr. - delta[, 1] <- MASS::mvrnorm( + delta[, 1] <- mvrnorm( n = 1, mu = matrix(data = 0, nrow = 1, ncol = n_subpopulations), Sigma = sigma_matrix_eps ) for (t_i in 2:n_time) { - eps_vec <- MASS::mvrnorm( + eps_vec <- mvrnorm( n = 1, mu = matrix(data = 0, nrow = 1, ncol = n_subpopulations), Sigma = sigma_matrix_eps diff --git a/R/wwinference-package.R b/R/wwinference-package.R index 7283a14..d0479b5 100644 --- a/R/wwinference-package.R +++ b/R/wwinference-package.R @@ -15,4 +15,5 @@ #' @importFrom rlang sym #' @importFrom stats dnbinom dweibull ecdf plogis qlogis rlnorm rnbinom rnorm #' rt sd time +#' @importFrom MASS mvrnorm NULL diff --git a/man/exponential_decay_corr_func.Rd b/man/exponential_decay_corr_func.Rd new file mode 100644 index 0000000..9d3f025 --- /dev/null +++ b/man/exponential_decay_corr_func.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/exponential_decay_corr_func.R +\name{exponential_decay_corr_func} +\alias{exponential_decay_corr_func} +\title{Generates the correlation matrix for an exponential +decay correlation function.} +\usage{ +exponential_decay_corr_func( + corr_function_params = list(dist_matrix = NULL, phi = NULL, l = NULL) +) +} +\arguments{ +\item{corr_function_params}{Sturcuted as follows : {dist.matrix, phi, l}.} +} +\value{ +Correlation matrix. +} +\description{ +Generates the correlation matrix for an exponential +decay correlation function. +} diff --git a/man/generate_simulated_data.Rd b/man/generate_simulated_data.Rd index 2cd6158..1f779e4 100644 --- a/man/generate_simulated_data.Rd +++ b/man/generate_simulated_data.Rd @@ -27,7 +27,13 @@ generate_simulated_data( mean_log_lod = 3.8, sd_log_lod = 0.2, input_params_path = fs::path_package("extdata", "example_params.toml", package = - "wwinference") + "wwinference"), + corr_function = independence_corr_func, + corr_fun_params = NULL, + phi_rt = 0.6, + sigma_eps = sqrt(0.02), + scaling_factor = NULL, + aux_site_bool = FALSE ) } \arguments{ @@ -98,6 +104,21 @@ LOD across sites} \item{input_params_path}{path to the toml file with the parameters to use to generate the simulated data} + +\item{corr_function}{Correlation function for spatial site correlations} + +\item{corr_fun_params}{Parameters required for correlation function} + +\item{phi_rt}{Coefficient for AR(1) temporal correlation on subpopulation +deviations.} + +\item{sigma_eps}{Parameter for construction of covariance matrix of spatial +epsilon.} + +\item{scaling_factor}{Scaling factor for aux site.} + +\item{aux_site_bool}{Boolean to use the aux site framework with +scaling factor.} } \value{ a list containing three dataframes. hosp_data is a dataframe diff --git a/man/independence_corr_func.Rd b/man/independence_corr_func.Rd new file mode 100644 index 0000000..7725812 --- /dev/null +++ b/man/independence_corr_func.Rd @@ -0,0 +1,17 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/independence_corr_func.R +\name{independence_corr_func} +\alias{independence_corr_func} +\title{Generates the correlation matrix for a non-spatial process.} +\usage{ +independence_corr_func(corr_function_params = list(num_sites = NULL)) +} +\arguments{ +\item{corr_function_params}{NULL} +} +\value{ +Correlation matrix of diagonal ones. +} +\description{ +Generates the correlation matrix for a non-spatial process. +} diff --git a/man/spatial_unadj_rt_process.Rd b/man/spatial_unadj_rt_process.Rd new file mode 100644 index 0000000..0447c1b --- /dev/null +++ b/man/spatial_unadj_rt_process.Rd @@ -0,0 +1,35 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/spatial_unadj_rt_process.R +\name{spatial_unadj_rt_process} +\alias{spatial_unadj_rt_process} +\title{Generates the subpopulation unadjusted Rt trajectories from proposal 2.} +\usage{ +spatial_unadj_rt_process( + log_state_rt, + corr_function, + corr_function_params, + phi_rt, + sigma_eps +) +} +\arguments{ +\item{corr_function}{Function for the correlation structure desired.} + +\item{corr_function_params}{List of parameters of desired correlation +function, should contain dist matrix.} + +\item{phi_rt}{Coefficient for AR(1) temporal correlation on +subpopulation deviations.} + +\item{sigma_eps}{Parameter for construction of covariance matrix +of spatial epsilon.} + +\item{state_rt}{"State" level unadjusted Rt, on log scale.} +} +\value{ +A matrix for unadjusted Rt where rows are subpopulations +and columns are time points. +} +\description{ +Generates the subpopulation unadjusted Rt trajectories from proposal 2. +} diff --git a/man/wwinference-package.Rd b/man/wwinference-package.Rd index 7520008..2103b37 100644 --- a/man/wwinference-package.Rd +++ b/man/wwinference-package.Rd @@ -17,10 +17,11 @@ Useful links: } \author{ -\strong{Maintainer}: Kaitlyn Johnson \email{uox1@cdc.gov} (\href{https://orcid.org/0000-0001-8011-0012}{ORCID}) +\strong{Maintainer}: Christian Bernal Zelaya \email{xuk0@cdc.gov} Authors: \itemize{ + \item Kaitlyn Johnson \email{uox1@cdc.gov} (\href{https://orcid.org/0000-0001-8011-0012}{ORCID}) \item Sam Abbott \email{contact@samabbott.co.uk} (\href{https://orcid.org/0000-0001-8057-8037}{ORCID}) \item Zachary Susswein \email{utb2@cdc.gov} \item Andrew Magee \email{rzg0@cdc.gov} From 96ab5aec4237c87cd48c720df6f8737619d71706 Mon Sep 17 00:00:00 2001 From: cbernalz Date: Tue, 16 Jul 2024 11:26:02 -0700 Subject: [PATCH 7/7] 2025-07-16 update : updating DESCRIPTION. --- DESCRIPTION | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 4b8746b..c8bac52 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -86,7 +86,8 @@ Imports: rlang, scales, ggplot2, - posterior + posterior, + MASS Remotes: stan-dev/cmdstanr VignetteBuilder: