Merge pull request #33 from sempwn/32-add-option-in-plot_kit_use-to-f…

…ilter-for-given-regions

32 add option in plot kit use to filter for given regions

R/utils-plot.R

@@ -6,10 +6,14 @@
 #' @param ... named list of [stanfit] objects
 #' @param data data used for model fitting. Can also include `p_use` column
 #' which can be used to plot true values if derived from simulated data.
+#' @param regions_to_plot Optional list to filter which regions are
+#' plotted
 #' @return [ggplot2] object
 #' @export
 #' @family plots
-plot_kit_use <- function(..., data = NULL) {
+plot_kit_use <- function(...,
+                         data = NULL,
+                         regions_to_plot = NULL) {
   combined_plot_data <- combine_model_fits(..., data = data)
 
   region <- model <- times <- sim_p <- p_use <- region_name <- NULL
@@ -20,6 +24,11 @@ plot_kit_use <- function(..., data = NULL) {
     combined_plot_data[["p_use"]] <- NA_real_
   }
 
+  if (!is.null(regions_to_plot)) {
+    combined_plot_data <- combined_plot_data %>%
+      dplyr::filter(region_name %in% regions_to_plot)
+  }
+
   combined_plot_data <- combined_plot_data %>%
     dplyr::group_by(region_name, model, times) %>%
     dplyr::summarise(

R/utils-table.R

@@ -45,6 +45,7 @@ kit_summary_table <- function(fit, ..., data = NULL,
 
   i <- sim_p <- times <- estimate <- Distributed <- Reported_Distributed <- NULL
   sim_used <- Reported_Used <- Orders <- NULL
+  .chain <- .draw <- .iteration <- value <- NULL
 
   percent_func <- scales::percent_format(accuracy = accuracy)
   comma_func <- scales::comma_format(accuracy = accuracy)

README.Rmd

@@ -63,7 +63,7 @@ mcmc_pairs(fit, pars = c("sigma","mu0","zeta"),
            off_diag_args = list(size = 1, alpha = 0.5))
 ```
 
-An overall summary of the model output can also be provided as a datframe
+An overall summary of the model output can also be provided as a data frame
 
 ```{r, warning=FALSE}
 kit_summary_table(fit, data = d)

vignettes/Introduction.Rmd

@@ -26,9 +26,11 @@ options(mc.cores = 2)
 ## basic example code
 d <- generate_model_data()
 # note iter should be at least 2000 to generate a reasonable posterior sample
-fit <- est_naloxone(d,iter=100)
-mcmc_pairs(fit, pars = c("sigma","mu0","zeta"),
-           off_diag_args = list(size = 1, alpha = 0.5))
+fit <- est_naloxone(d, iter = 100)
+mcmc_pairs(fit,
+  pars = c("sigma", "mu0", "zeta"),
+  off_diag_args = list(size = 1, alpha = 0.5)
+)
 ```
 
 Return median and 90th percentiles for posterior samples
@@ -39,11 +41,21 @@ library(posterior)
 summarise_draws(fit, default_summary_measures())
 ```
 
-We can compare two different model fits of the probability of prior kit use
-to simulated data using the `plot_kit_use` function
+We can plot the posterior predictive distribution of the probability of prior 
+kit use and compare to simulated data using the `plot_kit_use` function
 
 ```{r plot draws}
-plot_kit_use(model = fit,data=d)
+plot_kit_use(model = fit, data = d)
+```
+
+Alternatively, if we only want to plot certain regions we can include that
+as an optional argument in the `plot_kit_use` function
+
+```{r plot draws filter regions}
+plot_kit_use(
+  model = fit, data = d,
+  regions_to_plot = c("2")
+)
 ```
 
 We can compare the posterior predictive check to the prior predictive check
@@ -53,18 +65,21 @@ parameter to `FALSE`
 ```{r run prior}
 # note iter should be at least 2000 to generate a reasonable posterior sample
 prior_fit <- est_naloxone(d,
-                          run_estimation = FALSE,
-                          iter=100)
-mcmc_pairs(prior_fit, pars = c("sigma","mu0","zeta"),
-           off_diag_args = list(size = 1, alpha = 0.5))
+  run_estimation = FALSE,
+  iter = 100
+)
+mcmc_pairs(prior_fit,
+  pars = c("sigma", "mu0", "zeta"),
+  off_diag_args = list(size = 1, alpha = 0.5)
+)
 ```
 
 
 We can compare the prior and posterior predictive distribution of the 
 probability of prior kit use using the `plot_kit_use` function
 
 ```{r plot compare prior posterior}
-plot_kit_use(prior = prior_fit,posterior = fit, data=d)
+plot_kit_use(prior = prior_fit, posterior = fit, data = d)
 ```
 
 
@@ -86,10 +101,13 @@ naloxone use in time.
 # note iter should be at least 2000 to generate a reasonable posterior sample
 # note use `rw_type` to specify order of random walk
 rw2_fit <- est_naloxone(d,
-                    rw_type = 2,
-                    iter=100)
-mcmc_pairs(rw2_fit, pars = c("sigma","mu0","zeta"),
-           off_diag_args = list(size = 1, alpha = 0.5))
+  rw_type = 2,
+  iter = 100
+)
+mcmc_pairs(rw2_fit,
+  pars = c("sigma", "mu0", "zeta"),
+  off_diag_args = list(size = 1, alpha = 0.5)
+)
 ```
 
 
@@ -103,7 +121,7 @@ We can compare two different model fits of the probability of prior kit use
 to simulated data using the `plot_kit_use` function
 
 ```{r plot draws comparison}
-plot_kit_use(rw_1 = fit,rw_2 = rw2_fit,data=d)
+plot_kit_use(rw_1 = fit, rw_2 = rw2_fit, data = d)
 ```
 
 
@@ -118,25 +136,28 @@ distribution data once every three months
 ```{r data generation lower frequency}
 ## basic example code
 d_missing <- generate_model_data(reporting_freq = 3)
-ggplot(aes(x=times,y=Reported_Used,color=as.factor(regions)),data=d_missing) +
+ggplot(aes(x = times, y = Reported_Used, color = as.factor(regions)),
+  data = d_missing
+) +
   geom_point()
-
 ```
 
 
 
 ```{r missing data inference}
 missing_fit <- est_naloxone(d_missing,
-                    rw_type = 2,
-                    iter=100)
-mcmc_pairs(missing_fit, pars = c("sigma","mu0","zeta"),
-           off_diag_args = list(size = 1, alpha = 0.5))
-
+  rw_type = 2,
+  iter = 100
+)
+mcmc_pairs(missing_fit,
+  pars = c("sigma", "mu0", "zeta"),
+  off_diag_args = list(size = 1, alpha = 0.5)
+)
 ```
 
 
 ```{r plot draws with missing data}
-plot_kit_use(model = missing_fit, data=d_missing)
+plot_kit_use(model = missing_fit, data = d_missing)
 ```
 
 ## Change scale of data
@@ -153,22 +174,22 @@ To convert the delay distribution from month into week set $c = 1/4$ as one over
 # monthly reporting delay distribution
 psi_vec <- c(0.7, 0.2, 0.1)
 # convert to weeks using interpolation
-weekly_psi_vec <- rep(psi_vec,4) / sum(psi_vec)
+weekly_psi_vec <- rep(psi_vec, 4) / sum(psi_vec)
 
 # properties of order to distributed delay distribution in months
 max_delays <- 3
 delay_alpha <- 2
 delay_beta <- 1
 
 # convert to weeks
-weekly_max_delays <- max_delays*4
+weekly_max_delays <- max_delays * 4
 weekly_delay_alpha <- delay_alpha
-weekly_delay_beta <- 0.25 * delay_beta 
+weekly_delay_beta <- 0.25 * delay_beta
 
 result <- est_naloxone(d,
- psi_vec = weekly_psi_vec,
- max_delays = weekly_max_delays,
- delay_alpha = weekly_delay_alpha,
- delay_beta = weekly_delay_beta)
-
+  psi_vec = weekly_psi_vec,
+  max_delays = weekly_max_delays,
+  delay_alpha = weekly_delay_alpha,
+  delay_beta = weekly_delay_beta
+)
 ```