Test dev ggmcmc

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: serodynamics
 Title: What the Package Does (One Line, Title Case)
-Version: 0.0.0.9020
+Version: 0.0.0.9022
 Authors@R: c(
     person("Peter", "Teunis", , "p.teunis@emory.edu", role = c("aut", "cph"),
            comment = "Author of the method and original code."),
@@ -20,7 +20,7 @@ Imports:
     cli,
     coda,
     dplyr,
-    ggmcmc,
+    ggmcmc (>= 1.5.1.1),
     ggplot2,
     purrr,
     rlang,
@@ -60,8 +60,9 @@ Config/needs/test:
     readr,
     rlist
 Config/needs/check: commonmark, xml2
-Remotes: 
-    ucd-serg/serocalculator
+Remotes:  
+    ucd-serg/serocalculator,
+    xfim/ggmcmc
 Depends: 
     R (>= 4.1.0)
 LazyData: true

NAMESPACE

@@ -8,6 +8,7 @@ export(get_biomarker_names_var)
 export(initsfunction)
 export(load_data)
 export(plot_jags_dens)
+export(plot_jags_effect)
 export(postprocess_jags_output)
 export(prep_data)
 export(prep_priors)
@@ -19,6 +20,7 @@ importFrom(dplyr,.data)
 importFrom(dplyr,all_of)
 importFrom(dplyr,filter)
 importFrom(dplyr,mutate)
+importFrom(dplyr,n)
 importFrom(dplyr,reframe)
 importFrom(dplyr,rename)
 importFrom(dplyr,rowwise)

NEWS.md

@@ -4,6 +4,8 @@
 
 ## New features
 
+* Diagnostic function to produce effective sample size plots with
+stratification (#66)
 * Diagnostic function to produce density plots with stratification (#27)
 * Added SEES data set data folder and documentation (#41)
 * Fixing SEES data and added jags_post for SEES (#63)

R/plot_jags_effectivedx.R

@@ -0,0 +1,83 @@
+
+#' @title Plot Effective Sample Size Diagnostics
+#' @author Sam Schildhauer
+#' @description
+#'  plot_jags_effect() takes a [list] output from [serodynamics::run_mod()]
+#'  to create summary diagnostics for each chain run in the mcmc estimation.
+#'  Defaults will produce every combination of antigen/antibody, parameters,
+#'  and stratifications, unless otherwise specified. At least 2 chains are 
+#'  required to run function. 
+#'  Antigen/antibody combinations and stratifications will vary by analysis.
+#'  The antibody dynamic curve includes the following parameters:
+#'  - y0 = baseline antibody concentration
+#'  - y1 = peak antibody concentration
+#'  - t1 = time to peak
+#'  - r = shape parameter
+#'  - alpha = decay rate
+#' @param data A [list] outputted from run_mod().
+#' @param iso Specify [character] string to produce plots of only a
+#' specific antigen/antibody combination, entered with quotes. Default outputs
+#' all antigen/antibody combinations.
+#' @param param Specify [character] string to produce plots of only a
+#' specific parameter, entered with quotes. Options include:
+#' - `alpha` = posterior estimate of decay rate
+#' - `r` = posterior estimate of shape parameter
+#' - `t1` = posterior estimate of time to peak
+#' - `y0` = posterior estimate of baseline antibody concentration
+#' - `y1` = posterior estimate of peak antibody concentration
+#' @param strat Specify [character] string to produce plots of specific
+#' stratification entered in quotes.
+#' @return A [list] of [ggplot2::ggplot] objects showing the 
+#' proportion of effective samples taken/total samples taken for all parameter
+#' iso combinations. The estimate with the highest proportion of effective
+#' samples taken will be listed first.
+#' @export
+#' @example inst/examples/examples-plot_jags_effectivedx.R
+
+
+plot_jags_effect <- function(data,
+                             iso = unique(data$curve_params$Iso_type),
+                             param = unique(data$curve_params$Parameter_sub),
+                             strat = unique(data$curve_params$Stratification)) {
+  visualize_jags <- data[["curve_params"]]
+  attributes_jags <- data[["attributes"]]
+
+  eff_strat_list <- list()
+  for (i in strat) {
+
+    visualize_jags_sub <- visualize_jags |>
+      dplyr::filter(.data$Stratification == i)
+
+    # Creating open list to store ggplots
+    eff_out <- list()
+    # Looping through the isos
+    for (j in iso) {
+      visualize_jags_plot <- visualize_jags_sub |>
+        dplyr::filter(.data$Iso_type == j)
+
+      # Will not loop through parameters, as we may want each to show on the
+      # same plot by default.
+      visualize_jags_plot <- visualize_jags_plot |>
+        dplyr::filter(.data$Parameter_sub %in% param)
+
+      visualize_jags_plot <- visualize_jags_plot |>
+        # Changing parameter name to reflect the input
+        dplyr::mutate(Parameter = paste0("iso = ", j, ", parameter = ",
+                                         .data$Parameter_sub, ", strat = ",
+                                         i))
+      # Assigning attributes, which are needed to run ggs_density
+      attributes(visualize_jags_plot) <- c(attributes(visualize_jags_plot),
+                                           attributes_jags)
+      # Creating density plot
+      eff <- ggmcmc::ggs_effective(visualize_jags_plot) +
+        theme_bw()
+      eff_out[[j]] <- eff
+    }
+    eff_strat_list[[i]] <- eff_out
+  }
+  #Printing only one plot if only one exists.
+  if (length(eff_strat_list) == 1) {
+    eff_strat_list <- eff_strat_list[[1]][[iso]]
+  } 
+  eff_strat_list
+}

R/serodynamics-package.R

@@ -6,6 +6,7 @@
 #' @importFrom dplyr all_of
 #' @importFrom dplyr filter
 #' @importFrom dplyr mutate
+#' @importFrom dplyr n
 #' @importFrom dplyr reframe
 #' @importFrom dplyr rename
 #' @importFrom dplyr rowwise

inst/WORDLIST

@@ -9,7 +9,6 @@ bayesian
 behaviour
 biomarker
 biomarkers
-bldculres
 dobson
 geoms
 ggproto
@@ -18,7 +17,6 @@ iso
 isotype
 isotypes
 mcmc
-nepal
 newperson
 nmc
 params

inst/examples/examples-plot_jags_effectivedx.R

@@ -0,0 +1,10 @@
+
+if (!is.element(runjags::findjags(), c("", NULL))) {
+
+  data <- serodynamics::nepal_sees_jags_post
+
+  plot_jags_effect(data = data,
+                   iso = "HlyE_IgA",
+                   param = "alpha",
+                   strat = "typhi")
+}

pkgdown/_pkgdown.yml

@@ -32,6 +32,7 @@ reference:
 - title: Model diagnostics
   contents:
   - plot_jags_dens
+  - plot_jags_effect
 - title: Postprocess JAGS output
   contents:
   - postprocess_jags_output

tests/testthat/test-plot_jags_effectivedx.R

@@ -0,0 +1,20 @@
+
+test_that(
+  desc = "results are consistent with ggplot output showing ESS",
+  code = {
+    skip_if(getRversion() < "4.4.1") # 4.3.3 had issues
+
+    data <- serodynamics::nepal_sees_jags_post |>
+      suppressWarnings()
+
+    # Testing for any errors:
+    results <- plot_jags_effect(data) |> expect_no_error()
+
+    # Test to ensure output is a list object:
+    expect_true(is.list(results))
+
+    # Test to ensure that a piece of the list is a ggplot object:
+    results$typhi$HlyE_IgA |> 
+      vdiffr::expect_doppelganger(title = "typhoid_plot_ESS")
+  }
+)