diff --git a/DESCRIPTION b/DESCRIPTION
index 42714e6..1caa4d0 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -35,6 +35,7 @@ Imports:
     ggsci,
     lubridate,
     rlang,
+    tidyverse,
     dplyr
 SystemRequirements: C++11
 Suggests: 
diff --git a/NAMESPACE b/NAMESPACE
index f19e3d0..cbccaa1 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -12,29 +12,28 @@ export(get_haplotype_coalescence2)
 export(get_haplotype_relatedness)
 export(get_sample_relatedness)
 export(life_table_Mali)
-export(plot_EIR_prevalence)
 export(plot_age_states)
 export(plot_daily_prevalence)
+export(plot_epi_data)
 export(plot_epi_distribution)
 export(plot_treatment_seeking)
 export(prune_transmission_record)
-export(retrieve_prev)
 export(sim_block_tree)
 export(sim_epi)
 export(sim_haplotype_tree)
 export(sim_vcf)
 export(simplegen_file)
 export(simplegen_project)
+import("ggplot2,")
+import("ggpubr,")
+import("tidyverse,")
 import(tidyr)
 importFrom(Rcpp,sourceCpp)
 importFrom(data.table,fread)
-importFrom(ggsci,scale_color_lancet)
 importFrom(grDevices,col2rgb)
 importFrom(grDevices,grey)
 importFrom(grDevices,rgb)
 importFrom(magrittr,"%>%")
-importFrom(rlang,.data)
-importFrom(stats,aggregate)
 importFrom(stats,dbeta)
 importFrom(stats,dgeom)
 importFrom(stats,dpois)
diff --git a/R/misc.R b/R/misc.R
index ffe2fa2..e296a60 100644
--- a/R/misc.R
+++ b/R/misc.R
@@ -55,74 +55,4 @@ pairwise_long <- function(n, include_diagonal = FALSE) {
   return(ret)
 }
 
-#-----------------------------
-#' @title Retrieve prevalence
-#' 
-#' @description Function to add detection layer to estimated prevalence
-#'
-#' @param data Model output from SIMPLEGEN sim_epi() function epi_output$daily_values
-#' @param case_detection Method of case detection, "Active" or "Passive"
-#' @param diagnosis Method of diagnosis, "Microscopy" or "PCR"
-#' @param sampling_time Sampling time, currently input should be a numeric day e.g. 100 would represent the 100th day of simulation.
-#' @param sampled TO DO - sample from bernoulli distribution to generate simulated numbers of cases detected on each day - done in SIMPLEGEN internally already
-#' @param deme deme of interest - currently assumes just a single deme analysed - i.e. can't pool results of multiple demes
-#'
-#' @importFrom stats aggregate
-#' @return returns named vector :  c("annual_EIR", "prevalence_true", "prevalence_detected")
-#' @export
 
-
-retrieve_prev <- function(data,
-                          case_detection,
-                          diagnosis,
-                          sampling_time,
-                          sampled = FALSE,
-                          deme = 1) {
-  
-  results <- vector(length=3)
-  names(results) <- c("annual_EIR", "prevalence_true", "prevalence_detected")
-  
-  # Calculating Annual EIR --------------------------------------------------
-  
-  # subset to desired rows and columns
-  df_wide <- data[, c("time", "deme", "EIR")]
-  df_wide <- df_wide[df_wide$deme == deme,]
-  
-  df_wide$year <- lubridate::year(lubridate::as_date(df_wide$time, origin = lubridate::origin))
-  
-  annual_EIR <- stats::aggregate(EIR ~ year + deme, data = df_wide, FUN = "sum")
-  
-  sample_year <- lubridate::year(lubridate::as_date(sampling_time, origin = lubridate::origin))
-  
-  results["annual_EIR"] <- annual_EIR[which(annual_EIR$year == sample_year),"EIR"]
-  
-  
-  # Calculating true prevalence ---------------------------------------------
-  data <- data[data$time == sampling_time,]
-  
-  results["prevalence_true"] <- sum(data$C,data$A) / data$H
-  
-  # Calculating observed prevalence -----------------------------------------
-  
-  if (case_detection == "Active") {
-    if (diagnosis == "PCR") {
-      results["prevalence_detected"] <- sum(data$A_detectable_PCR, data$C_detectable_PCR) / data$H
-    } else if (diagnosis == "Microscopy") {
-      results["prevalence_detected"] <-sum(data$A_detectable_microscopy,data$C_detectable_microscopy) / data$H
-    } else {
-      warning("diagnosis type not recognized, please enter Microscopy or PCR")
-    }
-  } else if (case_detection == "Passive") {
-    if (diagnosis == "PCR") {
-      results["prevalence_detected"] <- data$A_detectable_PCR / data$H
-    } else if (diagnosis == "Microscopy") {
-      results["prevalence_detected"] <- data$A_detectable_microscopy / data$H
-    } else {
-      warning("diagnosis type not recognized, please enter Microscopy or PCR")
-    }
-  } else {
-    warning("case_detection type not recognized, please enter Active or Passive")
-  }
-  
-  return(results)
-}
diff --git a/R/plot.R b/R/plot.R
index ad98bcd..7277142 100644
--- a/R/plot.R
+++ b/R/plot.R
@@ -304,51 +304,299 @@ plot_age_states <- function(project, sample_time = NULL, deme = 1, state = "S")
 
 
 #------------------------------------------------
-#' @title Plot EIR versus prevalence
+#' Plotting SIMPLEGEN outputs against epidemiological data
 #'
-#' @description Plots relationship between annual EIR and prevalence
+#' @param model_out list of SIMPLEGEN model outputs from running sim_epi function
+#' @param dataname  Data to compare SIMPLEGEN output to. "micro_vs_PCR" for plotting
+#'  the relationship between prevalence by microscopy and prevalence by PCR.
+#'  "prevalence_EIR" is the relationship between Annual EIR and prevalence,
+#'  "prevalence_incidence" plots the relationship between Prevalence in 2-10 year olds
+#'   against incidence in 0-5 year olds, and the modelled relationship between these 
+#'   from Griffin et a. 2014
+#'   
+#'   TODO - Need to add Griffin et al. 2014  model fit data and example list of SIMPLEGEN
+#'   model outputs 
 #'
-#' @param data dataset produced by the retrieve_prev function
-#' @param plot_studies logical argument whether to plot data from previous
-#'   empirical studies (currently Hay et al., 2005)
-#' @param scale_x linear or log scale for x axis
 #'
-#' @importFrom ggsci scale_color_lancet
-#' @importFrom rlang .data
-#' @export
-#' 
-plot_EIR_prevalence <- function(data, plot_studies = TRUE, scale_x = "linear"){
+#' @export 
+#' @import tidyverse, ggpubr, ggplot2, 
+#'
+#' @examples
+#'
+#' plot_epi_data(model_out, "micro_vs_PCR")
+#' plot_epi_data(model_out, "prevalence_EIR")
+#' plot_epi_data(model_out, "prevalence_incidence")
+
+plot_epi_data <- function(model_out, dataname) {
+  
+  #check input
+  assert_class(model_out[[1]], "simplegen_project")
+  assert_class(model_out, "list")
+  assert_in(dataname, c("micro_vs_PCR", "prevalence_EIR", "prevalence_incidence"))
+  
+  
+  # Microscopy detection plot
+  if (dataname == "micro_vs_PCR") {
+    out <- lapply(
+      FUN = function(x) {
+        prev <- model_out[[x]]$epi_output$surveys %>%
+          filter(measure == "prevalence" &
+                   diagnostic == "microscopy") %>%
+          select(value) %>%
+          rename(prev_micro = value)
+        
+        incid <- model_out[[x]]$epi_output$surveys %>%
+          filter(measure == "prevalence" & diagnostic == "PCR") %>%
+          select(value) %>%
+          rename(prev_PCR = value)
+        
+        res <- cbind(prev, incid)
+        return(res)
+      },
+      X = 1:length(model_out)
+    )
+    
+    
+    data("PCR_micro_full_whittaker2021")
+    SIMPLEGEN_dat <- do.call(rbind.data.frame, out)
+    
+    
+    # error if nothing to plot
+    if (nrow(SIMPLEGEN_dat) == 0) {
+      stop("SIMPLEGEN output doesn't exist or incompatible with function")
+    }
+    
+    # error if nothing to plot
+    if (nrow(PCR_micro_full_whittaker2021) == 0) {
+      stop("data missing or formatted incorrectly")
+    }
+    
+    p <- ggplot2::ggplot() +
+      ggplot2::geom_point(
+        data = PCR_micro_full_whittaker2021,
+        aes(x = PCR_Prev, y = Micro_Prev, col = "Whittaker et al. 2021"),
+        size = 2
+      ) +
+      ggplot2:: geom_point(
+        data = SIMPLEGEN_dat,
+        aes(
+          x = prev_PCR / 100,
+          y = prev_micro / 100,
+          col = "SIMPLEGEN"
+        ),
+        size = 2
+      ) +
+      ggplot2::geom_abline(intercept = 0, lwd = 1) +
+      ggplot2::theme_bw() +
+      ggpubr::labs_pubr(base_size = 14) +
+      ggplot2::xlab("Prevalence by PCR") +
+      ggplot2::ylab("Prevalence by microscopy") +
+      ggplot2::labs(colour = "Source")
+    p
+    set_palette(p, "npg")
+    
+    return(p)
+  }
   
-  # load data (this is apparently not good practice!)
-  #data("EIRprev_hay2005")
-  #data("EIRprev_beier1999")
   
-  EIRprev <- tidyr::gather(as.data.frame(EIRprev), key = "detection_type", value = "prevalence", -1)
+  # Microscopy detection plot
+  if (dataname == "prevalence_EIR") {
+    # check ages for reported prevalence
+    
+    data("EIR_prev_hay2005")
+    #data("EIR_prev_beier1999")
+    
+    out <- lapply(function(x) {
+      # get prevalence
+      prev <- model_out[[x]]$epi_output$surveys %>%
+        filter(measure == 'prevalence' &
+                 diagnostic == "microscopy") %>%
+        select(value) %>%
+        rename(prev = value)
+      
+      
+      # get EIR
+      EIR <- model_out[[x]]$epi_output$daily %>%
+        filter(
+          measure == "EIR" &
+            time == unique(model_out[[x]]$epi_output$surveys$reporting_time)
+        ) %>%
+        select(value) %>%
+        rename(EIR = value)
+      
+      res <- data.frame(EIR = EIR,
+                        prev = mean(prev$prev, na.rm = TRUE))
+      return(res)
+      
+    }, X = 1:length(model_out))
+    
+    SIMPLEGEN_dat <- do.call(rbind.data.frame, out)
+    
+    # error if nothing to plot
+    if (nrow(SIMPLEGEN_dat) == 0) {
+      stop("SIMPLEGEN output doesn't exist or incompatible with function")
+    }
+    
+    # error if nothing to plot
+    if (nrow(EIR_prev_hay2005) == 0) {
+      stop("data missing or formatted incorrectly")
+    }
+    
+    
+    p <- ggplot() +
+      geom_point(
+        data = EIR_prev_hay2005,
+        aes(x = annual_EIR, y = prevalence, col = "Hay et al. 2005"),
+        size = 2
+      ) +
+      geom_point(data =  SIMPLEGEN_dat,
+                 aes(x = EIR, y = prev / 100, col = "SIMPLEGEN"),
+                 size = 2) +
+      theme_bw() +
+      labs_pubr(base_size = 14) +
+      xlim(c(0, 750)) +
+      xlab("Annual EIR") +
+      ylab("Prevalence by Microscopy") +
+      labs(colour = "Source")
+    p
+    set_palette(p, "npg")
+    return(p)
+  }
   
-  # remove zeros to avoid issues when logging
-  EIRprev <- EIRprev[is.finite(log(EIRprev[,"annual_EIR"])),]
   
-  if (plot_studies) {
-    EIRprev_hay2005$detection_type <- rep("Hay et al., 2005", length.out = nrow( EIRprev_hay2005))
-    EIRprev_beier1999$detection_type <- rep("Beier et al., 1999", length.out = nrow( EIRprev_beier1999))
-    EIRprev <- rbind( EIRprev, EIRprev_hay2005)
-    EIRprev <- rbind(EIRprev,EIRprev_beier1999)
-  }
-  if (scale_x == "linear") {
-    p <- ggplot2::ggplot(data = as.data.frame(EIRprev), ggplot2::aes(.data$annual_EIR, .data$prevalence, colour = .data$detection_type)) +
-      ggplot2::theme_bw() + ggsci::scale_color_lancet() +
-      ggplot2::geom_point() +
-      ggplot2::labs(title = "Annual EIR and Prevalence", x = "Annual EIR", y = "Prevalence") +
-      ggplot2::coord_cartesian(ylim = c(0,1), xlim = c(0,400))
-  } else if (scale_x == "log") {
-    p <- ggplot2::ggplot(as.data.frame(EIRprev), ggplot2::aes(.data$annual_EIR, .data$prevalence, colour = .data$detection_type)) +
+  # Prev-incidence relationship by age
+  if (dataname == "prevalence_incidence") {
+    # check ages for reported prevalence
+    
+    data("prev_inc_griffin2014")
+    
+    # Griffin Model data must be added to package
+    # data("griffinmod")
+    # griffin_mod$site_name <- "griffin_model"
+    # griffin_mod <- griffin_mod[, c(3, 1, 2)]
+    # colnames(griffin_mod) <- c("site_name", "Prev_2_10", "incid_0_5")
+    # griffin_mod$Prev_2_10 <- as.numeric(griffin_mod$Prev_2_10) / 100
+    
+    # adding data used to fit griffin et al
+    prev_inc_griffin2014$value <-
+      prev_inc_griffin2014$numer / prev_inc_griffin2014$denom
+    
+    # separate out prevalence and incidence
+    prev <- prev_inc_griffin2014 %>% filter(type == "prevalence")
+    incid <- prev_inc_griffin2014 %>% filter(type == "incidence")
+    
+    # for prevalence, active or passive detection
+    passive_det <- prev %>% filter(case_detection == "PCD")
+    active_det <- prev %>% filter(case_detection != "PCD")
+    
+    # filter to ages 2-10
+    active_prev <- prev_inc_griffin2014 %>%
+      filter(site_index %in% unique(active_det$site_index)) %>%
+      filter(age1 < 11) %>%
+      filter(age0 > 1) %>%
+      filter(type == 'prevalence')
+    
+    
+    # filter to age 0 - 5
+    incid <- prev_inc_griffin2014 %>%
+      filter(age1 < 6) %>%
+      filter(type == 'incidence')
+    
+    dat_reshaped <- as.data.frame(matrix(ncol = 3, nrow = 3))
+    colnames(dat_reshaped) <-
+      c("site_name", "Prev_0_10", "incid_0_5")
+    
+    for (i in 1:length(unique(active_prev$site_name))) {
+      site_n <- unique(active_prev$site_name)[i]
+      i_tmp <- incid %>% filter(site_name == site_n)
+      ap_tmp <- active_prev %>% filter(site_name  == site_n)
+      
+      dat_reshaped[i, "site_name"] <- site_n
+      dat_reshaped[i, "Prev_0_10"] <-
+        sum(ap_tmp$numer) / sum(ap_tmp$denom)
+      dat_reshaped[i, "incid_0_5"] <-
+        sum(i_tmp$numer) / sum(i_tmp$denom)
+      
+    }
+    
+    
+    
+    out <- lapply(function(x) {
+      prev210 <- model_out[[x]]$epi_output$surveys %>%
+        filter(measure == 'prevalence' &
+                 age_max == 10 & diagnostic == "PCR") %>%
+        select(value) %>%
+        rename(prev210 = value)
+      
+      
+      incid05 <- model_out[[x]]$epi_output$surveys %>%
+        filter(measure == 'incidence' & age_max == 5) %>%
+        select(value) %>%
+        rename(incid05 = value)
+      
+      
+      # incid05 <- model_out[[x]]$epi_output$daily %>%
+      #   filter(
+      #     measure == 'incidence' &
+      #       time <  model_out[[x]]$epi_output$surveys$reporting_time &
+      #       time >  (
+      #         model_out[[x]]$epi_output$surveys$reporting_time[1] - model_out[[x]]$epi_output$surveys$reporting_interval[1]
+      #       )
+      #   ) %>%
+      #   #select(value)%>%
+      #   summarize(incid05 = sum(value, na.rm = TRUE)/sum(denom%>%
+      #   as.data.frame()
+      # colnames(incid05)<- c("incid05")
+      
+      
+      res <- cbind(prev210, incid05)
+      return(res)
+      
+    }, X = 1:length(model_out))
+    
+    # combine model output and griffin et al 2014 data
+    SIMPLEGEN_dat <- do.call(rbind.data.frame, out)
+    
+    
+    SIMPLEGEN_dat$site_name <-
+      rep('SIMPLEGEN', length.out = nrow(SIMPLEGEN_dat))
+    dat_reshaped <-
+      dat_reshaped %>% mutate(site_name = gsub('_.*', '', site_name))
+    
+    # error if nothing to plot
+    if (nrow(SIMPLEGEN_dat) == 0) {
+      stop("SIMPLEGEN output doesn't exist or incompatible with function")
+    }
+    
+    # error if nothing to plot
+    if (nrow(dat_reshaped) == 0) {
+      stop("data missing or formatted incorrectly")
+    }
+    p <-
+      ggplot2::ggplot(dat_reshaped,
+                      ggplot2::aes(x =  Prev_0_10, y = incid_0_5, col = site_name)) +
       ggplot2::theme_bw() +
-      ggplot2::geom_point() + ggplot2::scale_x_log10() + ggsci::scale_color_lancet() +
-      ggplot2::labs(title = "Annual EIR (log scale) and Prevalence", x = "Annual EIR", y = "Prevalence") +
-      ggplot2::coord_cartesian(ylim = c(0,1))
-  } else {
-    warning("scale_x must be log or linear")
+      # ggplot2::theme(legend.position = "none") +
+      ggplot2::geom_point(size = 2) +
+      labs_pubr(base_size = 14) +
+      ggplot2::xlab("Prevalence in 2 - 10 year olds") +
+      ggplot2::ylab("Incidence in 0 - 5 year olds") +
+      ggplot2::ggtitle("Prevalence in 2-10 year olds plotted against incidence in 0-5 year olds")
+    # p <-
+    #   p + geom_line(data = griffin_mod,
+    #                 aes(x = Prev_2_10, y = incid_0_5),
+    #                 size = 1)
+    p <- p + geom_point(data = SIMPLEGEN_dat,
+                        aes(x = prev210 / 100, y = incid05),
+                        size = 2)
+    
+    set_palette(p, "npg")
+    return(p)
+    
   }
-  print(p)
+  
 }
 
+
+
+
diff --git a/man/plot_EIR_prevalence.Rd b/man/plot_EIR_prevalence.Rd
deleted file mode 100644
index 859031f..0000000
--- a/man/plot_EIR_prevalence.Rd
+++ /dev/null
@@ -1,19 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/plot.R
-\name{plot_EIR_prevalence}
-\alias{plot_EIR_prevalence}
-\title{Plot EIR versus prevalence}
-\usage{
-plot_EIR_prevalence(data, plot_studies = TRUE, scale_x = "linear")
-}
-\arguments{
-\item{data}{dataset produced by the retrieve_prev function}
-
-\item{plot_studies}{logical argument whether to plot data from previous
-empirical studies (currently Hay et al., 2005)}
-
-\item{scale_x}{linear or log scale for x axis}
-}
-\description{
-Plots relationship between annual EIR and prevalence
-}
diff --git a/man/plot_epi_data.Rd b/man/plot_epi_data.Rd
new file mode 100644
index 0000000..d2de17b
--- /dev/null
+++ b/man/plot_epi_data.Rd
@@ -0,0 +1,30 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plot.R
+\name{plot_epi_data}
+\alias{plot_epi_data}
+\title{Plotting SIMPLEGEN outputs against epidemiological data}
+\usage{
+plot_epi_data(model_out, dataname)
+}
+\arguments{
+\item{model_out}{list of SIMPLEGEN model outputs from running sim_epi function}
+
+\item{dataname}{Data to compare SIMPLEGEN output to. "micro_vs_PCR" for plotting
+the relationship between prevalence by microscopy and prevalence by PCR.
+"prevalence_EIR" is the relationship between Annual EIR and prevalence,
+"prevalence_incidence" plots the relationship between Prevalence in 2-10 year olds
+ against incidence in 0-5 year olds, and the modelled relationship between these 
+ from Griffin et a. 2014
+ 
+ TODO - Need to add Griffin et al. 2014  model fit data and example list of SIMPLEGEN
+ model outputs}
+}
+\description{
+Plotting SIMPLEGEN outputs against epidemiological data
+}
+\examples{
+
+plot_epi_data(model_out, "micro_vs_PCR")
+plot_epi_data(model_out, "prevalence_EIR")
+plot_epi_data(model_out, "prevalence_incidence")
+}
diff --git a/man/retrieve_prev.Rd b/man/retrieve_prev.Rd
deleted file mode 100644
index 07a4225..0000000
--- a/man/retrieve_prev.Rd
+++ /dev/null
@@ -1,34 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/misc.R
-\name{retrieve_prev}
-\alias{retrieve_prev}
-\title{Retrieve prevalence}
-\usage{
-retrieve_prev(
-  data,
-  case_detection,
-  diagnosis,
-  sampling_time,
-  sampled = FALSE,
-  deme = 1
-)
-}
-\arguments{
-\item{data}{Model output from SIMPLEGEN sim_epi() function epi_output$daily_values}
-
-\item{case_detection}{Method of case detection, "Active" or "Passive"}
-
-\item{diagnosis}{Method of diagnosis, "Microscopy" or "PCR"}
-
-\item{sampling_time}{Sampling time, currently input should be a numeric day e.g. 100 would represent the 100th day of simulation.}
-
-\item{sampled}{TO DO - sample from bernoulli distribution to generate simulated numbers of cases detected on each day - done in SIMPLEGEN internally already}
-
-\item{deme}{deme of interest - currently assumes just a single deme analysed - i.e. can't pool results of multiple demes}
-}
-\value{
-returns named vector :  c("annual_EIR", "prevalence_true", "prevalence_detected")
-}
-\description{
-Function to add detection layer to estimated prevalence
-}