Dev plot foi (#41)

* remove unused test datasets

* add option to plot additional plot data to plot_foi and plot_seromodel methods

* add test for the plot_foi method of the visualization module

* update documentation for plot_foi and plot_seromodel. Add new simulated dataset serodata_simD and add test simulated data serodata_simD.R

* add test test_plot_foi

* change simulated foi data plot in plot_foi from scattered to line

tests/testthat/extdata/data_test.R → R/serodata_simD.R

@@ -4,14 +4,12 @@
 #'
 #' @docType data
 #'
-#' @usage data_test
+#' @usage serodata_simD
 #'
 #' @format An object of class \code{"cross"}; see \code{\link[qtl]{read.cross}}.
 #'
 #' @keywords datasets
 #'
 #' @examples
-#' \dontrun{
-#' data_test
-#' }
-"data_test"
+#' serodata_simD
+"serodata_simD"

R/visualization.R

@@ -139,9 +139,9 @@ plot_seroprev_fitted <- function(seromodel_object,
 #' }
 #' @export
 plot_foi <- function(seromodel_object,
-                     lambda_sim = NA,
                      max_lambda = NA,
-                     size_text = 25) {
+                     size_text = 25,
+                     foi_sim = NULL) {
   if (is.character(seromodel_object$fit) == FALSE) {
     if (class(seromodel_object$fit@sim$samples) != "NULL") {
       foi <- rstan::extract(seromodel_object$fit,
@@ -151,18 +151,6 @@ plot_foi <- function(seromodel_object,
       #-------- This bit is to get the actual length of the foi data
       foi_data <- seromodel_object$foi_cent_est
 
-      if (!is.na(lambda_sim)) {
-        lambda_mod_length <- NROW(foi_data)
-        lambda_sim_length <- length(lambda_sim)
-
-        if (lambda_mod_length < lambda_sim_length) {
-          remove_x_values <- lambda_sim_length - lambda_mod_length
-          lambda_sim <- lambda_sim[-c(1:remove_x_values)]
-        }
-
-        foi_data$simulated <- lambda_sim
-      }
-
       #--------
       foi_data$medianv[1] <- NA
       foi_data$lower[1] <- NA
@@ -186,6 +174,19 @@ plot_foi <- function(seromodel_object,
         ggplot2::coord_cartesian(ylim = c(0, max_lambda)) +
         ggplot2::ylab("Force-of-Infection") +
         ggplot2::xlab("Year")
+      #TODO Add warning for foi_sim of different length than exposure years
+      if (!is.null(foi_sim)){
+        foi_data_length <- nrow(foi_data)
+        foi_sim_length <- length(foi_sim)
+        remove_x_values <- foi_sim_length - foi_data_length
+
+        foi_sim_data <- data.frame(year = foi_data$year, 
+                                  foi_sim = foi_sim[-c(1:remove_x_values)])    
+        foi_plot <- foi_plot + 
+          ggplot2::geom_line(data = foi_sim_data, ggplot2::aes(x = year, y = foi_sim),
+                            colour = "#b30909",
+                            size = size_text/8)
+      }
     }
   } else {
     print("model did not run")
@@ -296,19 +297,19 @@ plot_rhats <- function(seromodel_object,
 #' }
 #' @export
 plot_seromodel <- function(seromodel_object,
-                       lambda_sim = NA,
-                       max_lambda = NA,
-                       size_text = 25) {
+                          max_lambda = NA,
+                          size_text = 25,
+                          foi_sim = NULL) {
   if (is.character(seromodel_object$fit) == FALSE) {
     if (class(seromodel_object$fit@sim$samples) != "NULL") {
       prev_plot <- plot_seroprev_fitted(seromodel_object = seromodel_object,
                                  size_text = size_text)
 
       foi_plot <- plot_foi(
         seromodel_object = seromodel_object,
-        lambda_sim = lambda_sim,
         max_lambda = max_lambda,
-        size_text = size_text
+        size_text = size_text,
+        foi_sim = foi_sim
       )
 
       rhats_plot <- plot_rhats(seromodel_object = seromodel_object,

tests/testthat/test_plot_foi.R

@@ -0,0 +1,51 @@
+test_that("individual models", {
+  # So far we are skipping tests on these platforms until
+  # we find an efficient way to update rstan testthat snapshots on all of them
+  skip_on_os(c("windows", "mac"))
+  skip_on_ci()
+  source("testing_utils.R")
+  set.seed(1234) # For reproducibility
+
+  library(devtools)
+  library(dplyr)
+  library(vdiffr)
+
+  #----- Read and prepare data
+  data("serodata_simD")
+  simdata <- serodata_simD %>% prepare_serodata()
+  no_transm <- 0.0000000001
+  big_outbreak   <- 1.5
+  foi_sim <- c(rep(no_transm, 32), rep(big_outbreak, 3), rep(no_transm, 15)) # 1 epidemics
+
+
+  #----- Run models
+  models_to_run <- c("constant",
+                     "tv_normal",
+                     "tv_normal_log")
+  models_list <- lapply(models_to_run,
+                        run_seromodel,
+                        serodata = simdata,
+                        n_iters = 1000)
+
+  #----- Results visualisation
+  size_text <- 6
+  max_lambda <- 1.55
+  constant_plot <- plot_seromodel(models_list[[1]],
+                                  size_text = size_text,
+                                  max_lambda = max_lambda,
+                                  foi_sim = foi_sim)
+  tv_normal_plot <- plot_seromodel(models_list[[2]],
+                                   size_text = size_text,
+                                   max_lambda = max_lambda,
+                                   foi_sim = foi_sim)
+  tv_normal_log_plot <- plot_seromodel(models_list[[3]],
+                                       size_text = size_text,
+                                       max_lambda = max_lambda,
+                                       foi_sim = foi_sim)
+  plot_arrange <- cowplot::plot_grid(constant_plot,
+                                     tv_normal_plot,
+                                     tv_normal_log_plot,
+                                     ncol = 3, labels = "AUTO")
+  vdiffr::expect_doppelganger("plot_arrange_simdata_foi", plot_arrange)
+
+})