new functions and bug fixes

functions estRodis_init_params_model_one and estRodis_init_params_model_two added, models renamed and some bugs fixed

NAMESPACE

@@ -1,7 +1,9 @@
 # Generated by roxygen2: do not edit by hand
 
-export(estRodis_estimate_parameters_1)
-export(estRodis_estimate_parameters_2)
+export(estRodis_estimate_parameters_one)
+export(estRodis_estimate_parameters_two)
+export(estRodis_init_params_model_one)
+export(estRodis_init_params_model_two)
 export(estRodis_plot_transmission_chain)
 export(estRodis_plot_transmission_chain_mutation)
 export(estRodis_scaled_beta_distribution_pdf)

R/estRodis_estimate_parameters_1.R → R/estRodis_estimate_parameters_one.R

@@ -1,7 +1,7 @@
 
-#' @title Estimate parameters.
+#' @title Estimate parameters using Model 1
 #'
-#' @description `estRodis_estimate_parameters_1()` estimates the effective reproduction number,
+#' @description `estRodis_estimate_parameters_one()` estimates the effective reproduction number,
 #' the dispersion parameter, the probability of a case undergoing a mutation and
 #' the probability that a case is confirmed by a test based on the size distribution
 #' of identical sequence clusters.
@@ -36,7 +36,7 @@
 #' @param open_progress Parameter that is passed on to rstan::sampling, see manual of rstan::sampling for more details.
 #' @param show_messages Parameter that is passed on to rstan::sampling, see manual of rstan::sampling for more details.
 #'
-#' @details The core of the function `estRodis_estimate_parameters_1()` is a mathematical model
+#' @details The core of the function `estRodis_estimate_parameters_one()` is a mathematical model
 #' of the size distribution of sequence clusters, in which viral transmission,
 #' the mutation of the virus, and incomplete case-detection are integrated.
 #' Parameters are estimated by a Bayesian inference model implemented in Stan.
@@ -57,36 +57,36 @@
 #'
 #' options(mc.cores = parallelly::availableCores())
 #'
-#' estRodis_estimate_parameters_1(
+#' estRodis_estimate_parameters_one(
 #'   clusters_size = simulated_clusters$size,
 #'   clusters_freq = simulated_clusters$frequency,
 #'   sequencing_proba = 0.44)
 
-estRodis_estimate_parameters_1 <- function(clusters_size,
-                                           clusters_freq,
-                                           prior_r = c(10, 10),
-                                           prior_k = c(5, 10),
-                                           mean_generation_interval  = 5.2,
-                                           prior_number_yearly_mutations = c(14, 0.5),
-                                           prior_testing = c(1, 3, 0.05, 1),
-                                           sequencing_proba = 1,
-                                           pars = NA,
-                                           chains = 4,
-                                           iter = 2000,
-                                           warmup = floor(iter/2),
-                                           thin = 1,
-                                           seed = sample.int(.Machine$integer.max, 1),
-                                           init = 'random',
-                                           check_data = TRUE,
-                                           sample_file = NULL,
-                                           diagnostic_file = NULL,
-                                           verbose = FALSE,
-                                           algorithm = c("NUTS", "HMC", "Fixed_param"),
-                                           control = NULL,
-                                           include = TRUE,
-                                           cores = getOption("mc.cores", 1L),
-                                           open_progress = interactive() && !isatty(stdout()) && !identical(Sys.getenv("RSTUDIO"), "1"),
-                                           show_messages = TRUE) {
+estRodis_estimate_parameters_one <- function(clusters_size,
+                                             clusters_freq,
+                                             prior_r = c(10, 10),
+                                             prior_k = c(5, 10),
+                                             mean_generation_interval  = 5.2,
+                                             prior_number_yearly_mutations = c(14, 0.5),
+                                             prior_testing = c(1, 3, 0.05, 1),
+                                             sequencing_proba = 1,
+                                             pars = NA,
+                                             chains = 4,
+                                             iter = 2000,
+                                             warmup = floor(iter/2),
+                                             thin = 1,
+                                             seed = sample.int(.Machine$integer.max, 1),
+                                             init = 'random',
+                                             check_data = TRUE,
+                                             sample_file = NULL,
+                                             diagnostic_file = NULL,
+                                             verbose = FALSE,
+                                             algorithm = c("NUTS", "HMC", "Fixed_param"),
+                                             control = NULL,
+                                             include = TRUE,
+                                             cores = getOption("mc.cores", 1L),
+                                             open_progress = interactive() && !isatty(stdout()) && !identical(Sys.getenv("RSTUDIO"), "1"),
+                                             show_messages = TRUE) {
 
   standata <- list(M = length(clusters_size),
                    clusters_size = clusters_size,
@@ -98,7 +98,7 @@ estRodis_estimate_parameters_1 <- function(clusters_size,
                    prior_testing = prior_testing,
                    sequencing_proba = sequencing_proba)
 
-  out <- rstan::sampling(object = stanmodels$estRodis_stan_model_estimate_parameters_1,
+  out <- rstan::sampling(object = stanmodels$estRodis_stan_model_estimate_parameters_one,
                          data = standata,
                          pars = pars,
                          chains = chains,

R/estRodis_estimate_parameters_2.R → R/estRodis_estimate_parameters_two.R

@@ -1,7 +1,7 @@
 
-#' @title Estimate parameters.
+#' @title Estimate parameters using Model 2
 #'
-#' @description `estRodis_estimate_parameters_2()` estimates the effective reproduction number,
+#' @description `estRodis_estimate_parameters_two()` estimates the effective reproduction number,
 #' the dispersion parameter, the probability of a case undergoing a mutation and
 #' the probability that a case is confirmed by a test based on the size distribution
 #' of identical sequence clusters.
@@ -35,7 +35,7 @@
 #' @param open_progress Parameter that is passed on to rstan::sampling, see manual of rstan::sampling for more details.
 #' @param show_messages Parameter that is passed on to rstan::sampling, see manual of rstan::sampling for more details.
 #'
-#' @details The core of the function `estRodis_estimate_parameters_2()` is a mathematical model
+#' @details The core of the function `estRodis_estimate_parameters_two()` is a mathematical model
 #' of the size distribution of sequence clusters, in which viral transmission,
 #' the mutation of the virus, and incomplete case-detection are integrated.
 #' Parameters are estimated by a Bayesian inference model implemented in Stan.
@@ -55,37 +55,37 @@
 #'
 #' options(mc.cores = parallelly::availableCores())
 #'
-#' estRodis_estimate_parameters_2(
+#' estRodis_estimate_parameters_two(
 #'   clusters_size = simulated_clusters$size,
 #'   clusters_freq = simulated_clusters$frequency,
 #'   testing_proba = 0.55,
 #'   sequencing_proba = 0.44)
 
-estRodis_estimate_parameters_2 <- function(clusters_size,
-                                           clusters_freq,
-                                           prior_r = c(10, 10),
-                                           prior_k = c(5, 10),
-                                           mean_generation_interval  = 5.2,
-                                           prior_number_yearly_mutations = c(14, 0.5),
-                                           testing_proba = 1,
-                                           sequencing_proba = 1,
-                                           pars = NA,
-                                           chains = 4,
-                                           iter = 2000,
-                                           warmup = floor(iter/2),
-                                           thin = 1,
-                                           seed = sample.int(.Machine$integer.max, 1),
-                                           init = 'random',
-                                           check_data = TRUE,
-                                           sample_file = NULL,
-                                           diagnostic_file = NULL,
-                                           verbose = FALSE,
-                                           algorithm = c("NUTS", "HMC", "Fixed_param"),
-                                           control = NULL,
-                                           include = TRUE,
-                                           cores = getOption("mc.cores", 1L),
-                                           open_progress = interactive() && !isatty(stdout()) && !identical(Sys.getenv("RSTUDIO"), "1"),
-                                           show_messages = TRUE) {
+estRodis_estimate_parameters_two <- function(clusters_size,
+                                             clusters_freq,
+                                             prior_r = c(10, 10),
+                                             prior_k = c(5, 10),
+                                             mean_generation_interval  = 5.2,
+                                             prior_number_yearly_mutations = c(14, 0.5),
+                                             testing_proba = 1,
+                                             sequencing_proba = 1,
+                                             pars = NA,
+                                             chains = 4,
+                                             iter = 2000,
+                                             warmup = floor(iter/2),
+                                             thin = 1,
+                                             seed = sample.int(.Machine$integer.max, 1),
+                                             init = 'random',
+                                             check_data = TRUE,
+                                             sample_file = NULL,
+                                             diagnostic_file = NULL,
+                                             verbose = FALSE,
+                                             algorithm = c("NUTS", "HMC", "Fixed_param"),
+                                             control = NULL,
+                                             include = TRUE,
+                                             cores = getOption("mc.cores", 1L),
+                                             open_progress = interactive() && !isatty(stdout()) && !identical(Sys.getenv("RSTUDIO"), "1"),
+                                             show_messages = TRUE) {
 
   standata <- list(M = length(clusters_size),
                    clusters_size = clusters_size,
@@ -97,7 +97,7 @@ estRodis_estimate_parameters_2 <- function(clusters_size,
                    testing_proba = testing_proba,
                    sequencing_proba = sequencing_proba)
 
-  out <- rstan::sampling(object = stanmodels$estRodis_stan_model_estimate_parameters_2,
+  out <- rstan::sampling(object = stanmodels$estRodis_stan_model_estimate_parameters_two,
                          data = standata,
                          pars = pars,
                          chains = chains,

R/estRodis_init_params_model_one.R

@@ -0,0 +1,34 @@
+
+#' @title Initial values for Model 1
+#'
+#' @description Sample from prior distributions of R, k, the number of yearly mutations and the testing probability.
+#'
+#' @param prior_r Parameters for prior distribution of R (gamma)
+#' @param prior_k Parameters for prior distribution of k (gamma)
+#' @param prior_number_yearly_mutations Parameters for prior distribution of the number of yearly mutations (normal)
+#' @param prior_testing Parameters for prior distribution of testing probability (scaled beta)
+#'
+#' @return A named list containing samples of the prior distributions of R, k, the number of yearly mutations and the testing probability.
+#'
+#' @export
+#'
+#' @examples
+#' estRodis_init_params_model_one()
+
+estRodis_init_params_model_one <- function(prior_r = c(10, 10), prior_k = c(5, 10), prior_number_yearly_mutations = c(14, 0.5), prior_testing = c(1, 3, 0.05, 1)) {
+
+  distribution <- unlist(lapply(X = seq(prior_testing[3], prior_testing[4], 0.00001),
+                                FUN = function(x) estRodis_scaled_beta_distribution_pdf(prior_testing[1], prior_testing[2], prior_testing[3], prior_testing[4], x)))
+
+  # initial value for R sampled from gamma distribution(prior_r[1], prior_r[2])
+  # initial value for k sampled from gamma distribution(prior_k[1], prior_k[2])
+  # initial value for number_yearly_mutations sampled from normal distribution(prior_number_yearly_mutations[1], prior_number_yearly_mutations[2])
+  # initial value for testing_proba sampled from scaled beta distribution(prior_testing[1], prior_testing[2]) on the interval (prior_testing[3], prior_testing[4])
+  result <- list(R = stats::rgamma(1, prior_r[1], prior_r[2]),
+                 k = stats::rgamma(1, prior_k[1], prior_k[2]),
+                 number_yearly_mutations = stats::rnorm(1, prior_number_yearly_mutations[1], prior_number_yearly_mutations[2]),
+                 testing_proba =  sample(seq(prior_testing[3], prior_testing[4], 0.00001), 1, replace = TRUE, prob = distribution))
+
+  return(result)
+
+}

R/estRodis_init_params_model_two.R

@@ -0,0 +1,28 @@
+
+#' @title Initial values for Model 2
+#'
+#' @description Sample from prior distributions of R, k, the number of yearly mutations and the testing probability.
+#'
+#' @param prior_r Parameters for prior distribution of R (gamma)
+#' @param prior_k Parameters for prior distribution of k (gamma)
+#' @param prior_number_yearly_mutations Parameters for prior distribution of the number of yearly mutations (normal)
+#'
+#' @return A named list containing samples of the prior distributions of R, k and the number of yearly mutations.
+#'
+#' @export
+#'
+#' @examples
+#' estRodis_init_params_model_two()
+
+estRodis_init_params_model_two <- function(prior_r = c(10, 10), prior_k = c(5, 10), prior_number_yearly_mutations = c(14, 0.5)) {
+
+  # initial value for R sampled from gamma distribution(prior_r[1], prior_r[2])
+  # initial value for k sampled from gamma distribution(prior_k[1], prior_k[2])
+  # initial value for number_yearly_mutations sampled from normal distribution(prior_number_yearly_mutations[1], prior_number_yearly_mutations[2])
+  result <- list(R = stats::rgamma(1, prior_r[1], prior_r[2]),
+                 k = stats::rgamma(1, prior_k[1], prior_k[2]),
+                 number_yearly_mutations = stats::rnorm(1, prior_number_yearly_mutations[1], prior_number_yearly_mutations[2]))
+
+  return(result)
+
+}

R/estRodis_plot_transmission_chain.R

@@ -1,5 +1,4 @@
 
-
 #' @title Plot transmission chain
 #'
 #' @description A plot of a transmission chain is created.
@@ -8,6 +7,7 @@
 #' @param transmission_chain_edges Data frame containing information about the edges of the transmission chain
 #' @param max_generation Maximum generation of nodes and edges that are included into the plot
 #' @param colour_nodes_edges Colour for nodes and edges
+#' @param style_plot Positioning of transmission chain on canvas
 #'
 #' @details The data frame `transmission_chain_nodes` has to contain at least the columns `node_key` (id of node)
 #' and `generation` (generation of nodes).
@@ -17,22 +17,37 @@
 #' `to` (index of infectee node in `transmission_chain_nodes`)
 #' and `generation` (generation of infectee node).
 #'
+#' Using the inputs `max_generation` and `style_plot`,
+#' the graphical presentation of the transmission chain can be influenced.
+#'
+#' `max_generation` determines the last generation of nodes and edges
+#' that are included into the plot.
+#'
+#' Setting `style_plot = fixed`, nodes and edges up to (and including)
+#' generation `max_generation` are drawn,
+#' while nodes and edges of later generations are made invisible.
+#'
+#' Setting `style_plot = flexible`, nodes and edges up to (and including)
+#' generation `max_generation` are drawn,
+#' while nodes and edges of later generations are cut away.
+#'
 #' @return Plot of the transmission chain defined by nodes, edges and max_generation in the form of a ggplot2 object.
 #' @export
-#'
-#' @examples transmission_chain <- estRodis_simulate_transmission_chain()
-#' transmission_chain_nodes <- transmission_chain$nodes |> dplyr::select(c("node_key", "generation"))
-#' transmission_chain_edges <- transmission_chain$edges |> dplyr::select(c("from", "to", "generation"))
-#' plot_transmission_chain <- estRodis_plot_transmission_chain(
-#' transmission_chain_nodes = transmission_chain_nodes,
-#' transmission_chain_edges = transmission_chain_edges,
-#' max_generation = min(max(nodes$generation), 5))
+
 estRodis_plot_transmission_chain <- function(transmission_chain_nodes,
                                              transmission_chain_edges,
                                              max_generation = max(transmission_chain_nodes |> dplyr::pull("generation")),
                                              colour_nodes_edges = "navyblue",
                                              style_plot = "flexible") {
 
+  # #' @examples transmission_chain <- estRodis_simulate_transmission_chain()
+  # #' transmission_chain_nodes <- transmission_chain$nodes |> dplyr::select(c("node_key", "generation"))
+  # #' transmission_chain_edges <- transmission_chain$edges |> dplyr::select(c("from", "to", "generation"))
+  # #' plot_transmission_chain <- estRodis_plot_transmission_chain(
+  # #' transmission_chain_nodes = transmission_chain_nodes,
+  # #' transmission_chain_edges = transmission_chain_edges,
+  # #' max_generation = min(max(nodes$generation), 5))
+
   # determine nodes and edges that will be added to the plot
   if (style_plot == "fixed") {