diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml
index fbfc2e01..1e3a7e52 100644
--- a/.github/workflows/R-CMD-check.yaml
+++ b/.github/workflows/R-CMD-check.yaml
@@ -2,7 +2,7 @@
 # Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help
 on:
   push:
-    branches: [main, master]
+    branches: [main, master, dev]
   pull_request:
     branches: [main, master]
 
@@ -18,11 +18,11 @@ jobs:
       fail-fast: false
       matrix:
         config:
-          # - {os: macos-latest,   r: 'release'}
-          # - {os: windows-latest, r: 'release'}
-          # - {os: ubuntu-latest,   r: 'devel', http-user-agent: 'release'}
+          - {os: macos-latest,   r: 'release'}
+          - {os: windows-latest, r: 'release'}
+          - {os: ubuntu-latest,   r: 'devel', http-user-agent: 'release'}
           - {os: ubuntu-latest,   r: 'release'}
-          # - {os: ubuntu-latest,   r: 'oldrel-1'}
+          - {os: ubuntu-latest,   r: 'oldrel-1'}
 
     env:
       GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
diff --git a/.gitignore b/.gitignore
index aa3ed4e1..b5fb862f 100644
--- a/.gitignore
+++ b/.gitignore
@@ -40,6 +40,9 @@ vignettes/*.pdf
 .Rhistory
 
 
+# RStudio Connect folder
+rsconnect/
+
 # Mac File
 .DS_Store
 
@@ -53,5 +56,41 @@ _snaps/
 !tests/testthat/_snaps
 tests/testthat/Rplots.pdf
 t,
-tests/testthat/extdata/plots/actual/*.png
-tests/testthat/extdata/dataframes/actual/*.csv
\ No newline at end of file
+
+
+# Prerequisites
+*.d
+
+# Compiled Object files
+*.slo
+*.lo
+*.o
+*.obj
+
+# Precompiled Headers
+*.gch
+*.pch
+
+# Compiled Dynamic libraries
+*.so
+*.dylib
+*.dll
+
+# Fortran module files
+*.mod
+*.smod
+
+# Compiled Static libraries
+*.lai
+*.la
+*.a
+*.lib
+
+# Executables
+*.exe
+*.out
+*.app
+
+inst/stan/*.rds
+
+src/
diff --git a/DESCRIPTION b/DESCRIPTION
index 99062b0c..3ceccaa9 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -13,7 +13,7 @@ Authors@R:
         ),
         person(
             given = "Nicolás",
-            family = "T-Domínguez",
+            family = "T. Domínguez",
             role = c("aut"),
             email = "ex-ntorres@javeriana.edu.co",
             comment = c(ORCID = "https://orcid.org/0009-0002-8484-1298")
@@ -38,38 +38,47 @@ LazyData: true
 RoxygenNote: 7.2.3
 Depends: 
     R (>= 3.5.0)
-Imports:
-    config,
-    rstan (>= 2.21.1),
-    StanHeaders,
-    tidyverse,
-    reshape2,
-    bayesplot,
-    loo,
-    Hmisc,
+Imports: 
+    methods,
+    Rcpp (>= 0.12.0),
+    RcppParallel (>= 5.0.1),
+    rstan (>= 2.18.1),
+    rstantools (>= 2.3.1),
     dplyr,
-    gsubfn,
-    usethis,
-    testthat (>= 3.0.0),
-    vdiffr (>= 1.0.0),
+    bayesplot,
+    config,
+    cowplot,
     devtools,
-    methods,
-    Rcpp,
     ggplot2,
-    BH,
-    RcppEigen,
-    RcppParallel,
+    gsubfn,
+    Hmisc,
+    jsonlite,
+    loo,
     purrr,
-    cowplot
+    qtl,
+    reshape2,
+    tidyverse,
+    usethis,
+    vdiffr (>= 1.0.0)    
+LinkingTo: 
+    BH (>= 1.66.0),
+    Rcpp (>= 0.12.0),
+    RcppEigen (>= 0.3.3.3.0),
+    RcppParallel (>= 5.0.1),
+    rstan (>= 2.18.1),
+    StanHeaders (>= 2.18.0)
 Suggests: 
     knitr,
-    rmarkdown
+    rmarkdown,
+    testthat (>= 3.0.0)
+Config/testthat/edition: 3
 Config/Needs/website:
     epiverse-trace/epiversetheme
 VignetteBuilder: knitr
 URL: https://trace-lac.github.io/serofoi/
 Additional_repositories:
     https://mc-stan.org/r-packages/
-Config/testthat/edition: 3
 Remotes:
     tidyverse/purrr
+Biarch: true
+SystemRequirements: GNU make
diff --git a/NAMESPACE b/NAMESPACE
index ea0f3e9d..f5ffc57b 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -15,7 +15,6 @@ export(plot_seroprev_fitted)
 export(prepare_bin_data)
 export(prepare_serodata)
 export(run_seromodel)
-export(save_or_load_model)
 import(Rcpp)
 import(dplyr)
 import(methods)
@@ -24,3 +23,4 @@ importFrom(graphics,text)
 importFrom(rstan,sampling)
 importFrom(stats,quantile)
 importFrom(utils,read.table)
+useDynLib(serofoi, .registration = TRUE)
diff --git a/R/chagas2012.R b/R/chagas2012.R
index e58bc111..ac00f646 100644
--- a/R/chagas2012.R
+++ b/R/chagas2012.R
@@ -1,3 +1,4 @@
+# TODO Check if we really need to have the package `qtl` installed. Otherwise remove all entries of the form `see \code{\link[qtl]...`
 #' Seroprevalence data on serofoi
 #'
 #' Data from a serological surveys
diff --git a/R/modelling.R b/R/modelling.R
index efac7f11..8558f109 100644
--- a/R/modelling.R
+++ b/R/modelling.R
@@ -1,6 +1,9 @@
-#' Function that runs the specified stan model for the Force-of-Infection and estimates de seroprevalence based on the result of the fit 
-#' 
-#' This function runs the specified model for the Force-of-Infection \code{foi_model} using the data froma seroprevalence survey 
+  # TODO Complete @param documentation
+
+
+#' Function that runs the specified stan model for the Force-of-Infection and estimates de seroprevalence based on the result of the fit
+#'
+#' This function runs the specified model for the Force-of-Infection \code{foi_model} using the data froma seroprevalence survey
 #' \code{serodata} as the input data. See \link{fit_seromodel} for further details.
 #'
 #' @param serodata A data frame containing the data from a seroprevalence survey.
@@ -36,6 +39,7 @@
 #' For further details refer to the \code{control} parameter in \link[rstan]{sampling} or \href{https://mc-stan.org/rstanarm/reference/adapt_delta.html}{here}.
 #' @param m_treed Maximum tree depth for the binary tree used in the NUTS stan sampler. For further details refer to the \code{control} parameter in \link[rstan]{sampling}.
 #' @param decades Number of decades covered by the survey data.
+#' @param print_summary TBD
 #' @return \code{seromodel_object}. An object containing relevant information about the implementation of the model. For further details refer to \link{fit_seromodel}.
 #' @examples
 #' \dontrun{
@@ -69,46 +73,11 @@ run_seromodel <- function(serodata,
   return(seromodel_object)
 }
 
-# TODO The warning 'recompiling to avoid crashing R session' still appears when the function is run for a second time.
-#' Function used to determine whether the stan model corresponding to the specified serological model has been already compiled or not
-#'
-#' This function determines whether the corresponding .RDS file of the selected model exists or not.
-#' In case the .RDS file exists, it is read and returned; otherwise, the object model is created through the
-#' \link[rstan]{stan_model} function, saved as an .RDS file and returned as the output of the function.
-#' @param foi_model Name of the selected model. Current version provides three options:
-#' \describe{
-#' \item{\code{"constant"}}{Runs a constant model}
-#' \item{\code{"tv_normal"}}{Runs a normal model}
-#' \item{\code{"tv_normal_log"}}{Runs a normal logarithmic model}
-#' }
-#' @return \code{model}. The rstan model object corresponding to the selected model.
-#' @examples
-#' \dontrun{
-#' model <- save_or_load_model(foi_model="constant")
-#' }
-#' 
-#' @export
-
-save_or_load_model <- function(foi_model = "constant") {
-  base_path <- config::get("stan_models_base_path",
-    file = system.file("config.yml", package = "serofoi", mustWork = TRUE))
-  rds_path <- system.file(base_path, paste(foi_model, ".rds", sep = ""), package = getPackageName())
-  if (!file.exists(rds_path)) {
-    message(sprintf("\nNo rds file found for model %s. Compiling stan model...", foi_model))
-  }
-  stan_path <- system.file(base_path, paste(foi_model, ".stan", sep = ""), package = getPackageName())
-
-  model <- rstan::stan_model(stan_path, auto_write = TRUE)
-
-  return(model)
-}
-
-
 #' Function that fits the selected model to the specified seroprevalence survey data
-#' 
-#' This function fits the specified model \code{foi_model} to the serological survey data \code{serodata} 
-#' by means of the \link[rstan]{sampling} method. The function determines whether the corresponding stan model 
-#' object needs to be compiled by means of the function \link{save_or_load_model}. 
+#'
+#' This function fits the specified model \code{foi_model} to the serological survey data \code{serodata}
+#' by means of the \link[rstan]{sampling} method. The function determines whether the corresponding stan model
+#' object needs to be compiled by rstan.
 #' @param serodata A data frame containing the data from a seroprevalence survey. For further details refer to \link{run_seromodel}.
 #' @param foi_model Name of the selected model. Current version provides three options:
 #' \describe{
@@ -160,7 +129,7 @@ fit_seromodel <- function(serodata,
                           m_treed = 10,
                           decades = 0) {
   # TODO Add a warning because there are exceptions where a minimal amount of iterations is needed
-  model <- save_or_load_model(foi_model)
+  model <- stanmodels[[foi_model]]
   exposure_ages <- get_exposure_ages(serodata)
   exposure_years <- (min(serodata$birth_year):serodata$tsur[1])[-1]
   exposure_matrix <- get_exposure_matrix(serodata)
@@ -280,7 +249,7 @@ fit_seromodel <- function(serodata,
 
 #' Function that generates an atomic vector containing the corresponding exposition years of a serological survey
 #'
-#' This function generates an atomic vector containing the exposition years corresponding to the specified serological survey data \code{serodata}. 
+#' This function generates an atomic vector containing the exposition years corresponding to the specified serological survey data \code{serodata}.
 #' The exposition years to the disease for each individual corresponds to the time from birth to the moment of the survey.
 #' @param serodata A data frame containing the data from a seroprevalence survey. This data frame must contain the year of birth for each individual (birth_year) and the time of the survey (tsur). birth_year can be constructed by means of the \link{prepare_serodata} function.
 #' @return \code{exposure_ages}. An atomic vector with the numeration of the exposition years in serodata
@@ -321,15 +290,15 @@ get_exposure_matrix <- function(serodata) {
 
 #' Method to extact a summary of the specified serological model object
 #'
-#' This method extracts a summary corresponding to a serological model object that contains information about the original serological 
-#' survey data used to fit the model, such as the year when the survey took place, the type of test taken and the corresponding antibody, 
-#' as well as information about the convergence of the model, like the expected log pointwise predictive density \code{elpd} and its 
+#' This method extracts a summary corresponding to a serological model object that contains information about the original serological
+#' survey data used to fit the model, such as the year when the survey took place, the type of test taken and the corresponding antibody,
+#' as well as information about the convergence of the model, like the expected log pointwise predictive density \code{elpd} and its
 #' corresponding standar deviation.
-#' @param seromodel_object \code{seromodel_object}. An object containing relevant information about the implementation of the model. 
+#' @param seromodel_object \code{seromodel_object}. An object containing relevant information about the implementation of the model.
 #' Refer to \link{fit_seromodel} for further details.
 #' @return \code{model_summary}. Object with a summary of \code{seromodel_object} containing the following:
 #' \tabular{ll}{
-#' \code{foi_model} \tab Name of the selected model. For further details refer to \link{save_or_load_model}. \cr \tab \cr
+#' \code{foi_model} \tab Name of the selected model. \cr \tab \cr
 #' \code{data_set} \tab Seroprevalence survey label.\cr \tab \cr
 #' \code{country} \tab Name of the country were the survey was conducted in. \cr \tab \cr
 #' \code{year} \tab Year in which the survey was conducted. \cr \tab \cr
@@ -385,14 +354,15 @@ extract_seromodel_summary <- function(seromodel_object) {
   return(model_summary)
 }
 
-
-#' Function that generates an object containing the confidence interval based on a 
+# TODO Complete @param documentation
+#' Function that generates an object containing the confidence interval based on a
 #' Force-of-Infection fitting
 #'
 #' This function computes the corresponding binomial confidence intervals for the obtained prevalence based on a fitting
-#' of the Force-of-Infection \code{foi} for plotting an analysis purposes. 
-#' @param serodata A data frame containing the data from a seroprevalence survey. For further details refer to \link{run_seromodel}.
+#' of the Force-of-Infection \code{foi} for plotting an analysis purposes.
 #' @param foi Object containing the information of the force of infection. It is obtained from \code{rstan::extract(seromodel_object$fit, "foi", inc_warmup = FALSE)[[1]]}.
+#' @param serodata A data frame containing the data from a seroprevalence survey. For further details refer to \link{run_seromodel}.
+#' @param bin_data TBD
 #' @return \code{prev_final}. The expanded prevalence data. This is used for plotting purposes in the \code{visualization} module.
 #' @examples
 #' \dontrun{
@@ -404,8 +374,10 @@ extract_seromodel_summary <- function(seromodel_object) {
 #' }
 #' @export
 get_prev_expanded <- function(foi,
-                              serodata) {
+                              serodata,
+                              bin_data = FALSE) {
   dim_foi <- dim(foi)[2]
+  # TODO: check whether this conditional is necessary
   if (dim_foi < 80) {
     oldest_year <- 80 - dim_foi + 1
     foin <- matrix(NA, nrow = dim(foi)[1], 80)
@@ -438,7 +410,7 @@ get_prev_expanded <- function(foi,
   medianv <- apply(prev_pn, 2, function(x) quantile(x, 0.5))
 
   predicted_prev <- data.frame(
-    age = 1:80,
+    age = 1:age_max,
     predicted_prev = medianv,
     predicted_prev_lower = lower,
     predicted_prev_upper = upper
@@ -460,22 +432,23 @@ get_prev_expanded <- function(foi,
                 observed_prev,
                 by = "age",
                 all.x = TRUE) %>% dplyr::mutate(survey = serodata$survey[1])
-
-  # I added this here for those cases when binned is prefered for plotting
-  if (serodata$age_max[1] - serodata$age_min[1] < 3) {
-  xx <- prepare_bin_data(serodata)
-    prev_final <-
-      base::merge(prev_expanded, xx, by = "age", all.x = TRUE)
-  } else {
-    prev_final <- prev_expanded %>% dplyr::mutate(
-      cut_ages = "original",
-      bin_size = .data$sample_by_age,
-      bin_pos = .data$positives,
-      p_obs_bin = .data$prev_obs,
-      p_obs_bin_l = .data$prev_obs_lower,
-      p_obs_bin_u = .data$prev_obs_upper
-    )
+  if (bin_data) {
+    # I added this here for those cases when binned is prefered for plotting
+    if (serodata$age_max[1] - serodata$age_min[1] < 3) {
+    xx <- prepare_bin_data(serodata)
+      prev_expanded <-
+        base::merge(prev_expanded, xx, by = "age", all.x = TRUE)
+    } else {
+      prev_expanded <- prev_expanded %>% dplyr::mutate(
+        cut_ages = "original",
+        bin_size = .data$sample_by_age,
+        bin_pos = .data$positives,
+        p_obs_bin = .data$prev_obs,
+        p_obs_bin_l = .data$prev_obs_lower,
+        p_obs_bin_u = .data$prev_obs_upper
+      )
+    }
   }
 
-  return(prev_final)
+  return(prev_expanded)
 }
diff --git a/R/serofoi-package.R b/R/serofoi-package.R
new file mode 100644
index 00000000..defa92f7
--- /dev/null
+++ b/R/serofoi-package.R
@@ -0,0 +1,21 @@
+#' The 'serofoi' package.
+#'
+#' @description A DESCRIPTION OF THE PACKAGE
+#'
+#' @docType package
+#' @name serofoi-package
+#' @aliases serofoi
+#' @useDynLib serofoi, .registration = TRUE
+#' @import dplyr
+#' @importFrom dplyr %>%
+
+#' @import methods
+#' @import Rcpp
+#' @importFrom rstan sampling
+#' @importFrom graphics text
+#' @importFrom utils read.table
+#' @importFrom stats quantile
+#' @references
+#' Stan Development Team (NA). RStan: the R interface to Stan. R package version 2.26.22. https://mc-stan.org
+#'
+NULL
diff --git a/R/serofoi_package.R b/R/serofoi_package.R
deleted file mode 100644
index de07169d..00000000
--- a/R/serofoi_package.R
+++ /dev/null
@@ -1,14 +0,0 @@
-#' @keywords internal
-"_PACKAGE"
-
-## usethis namespace: start
-#' @import dplyr
-#' @importFrom dplyr %>%
-#' @import Rcpp
-#' @import methods
-#' @importFrom graphics text
-#' @importFrom utils read.table
-#' @importFrom rstan sampling
-#' @importFrom stats quantile
-## usethis namespace: end
-NULL
diff --git a/R/seroprevalence_data.R b/R/seroprevalence_data.R
index 29352522..f4168cfe 100644
--- a/R/seroprevalence_data.R
+++ b/R/seroprevalence_data.R
@@ -1,3 +1,6 @@
+# TODO Complete @param documentation
+
+
 #' Function that prepares the data from a serological survey for modelling
 #'
 #' This function adds the necessary additional variables to the given dataset \code{serodata} corresponding to a serological survey.
@@ -14,7 +17,8 @@
 #' \code{test} \tab The type of test taken \cr \tab \cr
 #' \code{antibody} \tab antibody \cr \tab \cr
 #' }
-#' @param alpha probability of a type I error. For further details refer to \link{Hmisc::binconf}.
+#' @param alpha probability of a type I error. For further details refer to \link[Hmisc]{binconf}.
+#' @param add_age_mean_f TBD
 #' @return serodata with additional columns necessary for the analysis. These columns are:
 #' \tabular{ll}{
 #' \code{age_mean_f} \tab Floor value of the average between age_min and age_max \cr \tab \cr
diff --git a/R/stanmodels.R b/R/stanmodels.R
new file mode 100644
index 00000000..4d385063
--- /dev/null
+++ b/R/stanmodels.R
@@ -0,0 +1,27 @@
+# Generated by rstantools.  Do not edit by hand.
+
+# names of stan models
+stanmodels <- c("constant", "tv_normal", "tv_normal_log")
+
+# load each stan module
+Rcpp::loadModule("stan_fit4constant_mod", what = TRUE)
+Rcpp::loadModule("stan_fit4tv_normal_mod", what = TRUE)
+Rcpp::loadModule("stan_fit4tv_normal_log_mod", what = TRUE)
+
+# instantiate each stanmodel object
+stanmodels <- sapply(stanmodels, function(model_name) {
+  # create C++ code for stan model
+  stan_file <- if(dir.exists("stan")) "stan" else file.path("inst", "stan")
+  stan_file <- file.path(stan_file, paste0(model_name, ".stan"))
+  stanfit <- rstan::stanc_builder(stan_file,
+                                  allow_undefined = TRUE,
+                                  obfuscate_model_name = FALSE)
+  stanfit$model_cpp <- list(model_cppname = stanfit$model_name,
+                            model_cppcode = stanfit$cppcode)
+  # create stanmodel object
+  methods::new(Class = "stanmodel",
+               model_name = stanfit$model_name,
+               model_code = stanfit$model_code,
+               model_cpp = stanfit$model_cpp,
+               mk_cppmodule = function(x) get(paste0("rstantools_model_", model_name)))
+})
diff --git a/R/visualization.R b/R/visualisation.R
similarity index 97%
rename from R/visualization.R
rename to R/visualisation.R
index 35147840..26f7fd81 100644
--- a/R/visualization.R
+++ b/R/visualisation.R
@@ -67,7 +67,7 @@ plot_seroprev_fitted <- function(seromodel_object,
     if  (class(seromodel_object$fit@sim$samples)  != "NULL" ) {
 
       foi <- rstan::extract(seromodel_object$fit, "foi", inc_warmup = FALSE)[[1]]
-      prev_expanded <- get_prev_expanded(foi, serodata = seromodel_object$serodata)
+      prev_expanded <- get_prev_expanded(foi, serodata = seromodel_object$serodata, bin_data = TRUE)
       prev_plot <-
         ggplot2::ggplot(prev_expanded) +
         ggplot2::geom_ribbon(
@@ -119,6 +119,8 @@ plot_seroprev_fitted <- function(seromodel_object,
   return(prev_plot)
 }
 
+# TODO Complete @param documentation
+
 #' Function that generates a Force-of-Infection plot corresponding to the specified fitted serological model
 #'
 #' This function generates a Force-of-Infection plot from the results obtained by fitting a serological model.
@@ -126,6 +128,8 @@ plot_seroprev_fitted <- function(seromodel_object,
 #' The x axis corresponds to the decades covered by the survey the y axis to the Force-of-Infection.
 #' @param seromodel_object Object containing the results of fitting a model by means of \link{run_seromodel}.
 #' @param size_text Text size use in the theme of the graph returned by the function.
+#' @param max_lambda TBD
+#' @param foi_sim TBD
 #' @return A ggplot2 object containing the Force-of-infection-vs-time including the corresponding confidence interval.
 #' @examples
 #' \dontrun{
@@ -278,12 +282,15 @@ plot_rhats <- function(seromodel_object,
   return(rhats_plot)
 }
 
+# TODO Complete @param documentation
 
 #' Function that generates a vertical arrange of plots showing a summary of a model, the estimated seroprevalence,
 #' the Force-of-Infection fit and the R-hat estimates plots.
 #'
 #' @param seromodel_object Object containing the results of fitting a model by means of \link{run_seromodel}.
 #' @param size_text Text size use in the theme of the graph returned by the function.
+#' @param max_lambda TBD
+#' @param foi_sim TBD
 #' @return A ggplot object with a vertical arrange containing the seropositivity, force of infection, and convergence plots.
 #' @examples
 #' \dontrun{
@@ -362,12 +369,13 @@ plot_seromodel <- function(seromodel_object,
   return(plot_arrange)
 }
 
-
+# TODO Improve documentation of @return. 
+# TODO Give more details about the generated plot
 #' Function that generates a plot for a given table
 #'
 #' @param info the information that will be contained in the table
 #' @param size_text Text size of the graph returned by the function
-#' @return p, a variable that will be used in the \link{visualisation} module
+#' @return p the plot for the given table
 #' @examples
 #' \dontrun{
 #'  data_test <- prepare_serodata(serodata)
diff --git a/README.Rmd b/README.Rmd
index 42ea15f7..e15802fb 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -92,7 +92,7 @@ More details on how to use ***serofoi*** can be found in the
 [online documentation](https://epiverse-trace.github.io/serofoi/) as package vignettes, under 
 [**Get Started**](https://epiverse-trace.github.io/serofoi/articles/serofoi.html), 
 [**An Introduction to FoI Models**](https://epiverse-trace.github.io/serofoi/articles/foi_models.html) and 
-[**Real-life Use Cases for _serofoi_**](https://epiverse-trace.github.io/serofoi/articles/use_cases.html)
+[**Real-life Use Cases for serofoi**](https://epiverse-trace.github.io/serofoi/articles/use_cases.html)
 
 ## Help
 
diff --git a/configure b/configure
new file mode 100755
index 00000000..0304fc54
--- /dev/null
+++ b/configure
@@ -0,0 +1,5 @@
+#! /bin/sh
+
+# Generated by rstantools.  Do not edit by hand.
+
+"${R_HOME}/bin/Rscript" -e "rstantools::rstan_config()"
diff --git a/configure.win b/configure.win
new file mode 100755
index 00000000..5e2dceb8
--- /dev/null
+++ b/configure.win
@@ -0,0 +1,5 @@
+#! /bin/sh
+
+# Generated by rstantools.  Do not edit by hand.
+
+"${R_HOME}/bin${R_ARCH_BIN}/Rscript.exe" -e "rstantools::rstan_config()"
diff --git a/inst/config.yml b/inst/config.yml
deleted file mode 100644
index 9fc90915..00000000
--- a/inst/config.yml
+++ /dev/null
@@ -1,2 +0,0 @@
-default:
-  stan_models_base_path: "extdata/stanmodels"
diff --git a/inst/desktop.ini b/inst/desktop.ini
deleted file mode 100644
index ab17096e..00000000
--- a/inst/desktop.ini
+++ /dev/null
@@ -1,4 +0,0 @@
-[ViewState]
-Mode=
-Vid=
-FolderType=Documents
diff --git a/data/chagas2012.RDS b/inst/extdata/chagas2012.RDS
similarity index 100%
rename from data/chagas2012.RDS
rename to inst/extdata/chagas2012.RDS
diff --git a/data/chik2015.RDS b/inst/extdata/chik2015.RDS
similarity index 100%
rename from data/chik2015.RDS
rename to inst/extdata/chik2015.RDS
diff --git a/data/serodata.RDS b/inst/extdata/serodata.RDS
similarity index 100%
rename from data/serodata.RDS
rename to inst/extdata/serodata.RDS
diff --git a/data/simdata_constant.RDS b/inst/extdata/simdata_constant.RDS
similarity index 100%
rename from data/simdata_constant.RDS
rename to inst/extdata/simdata_constant.RDS
diff --git a/data/simdata_large_epi.RDS b/inst/extdata/simdata_large_epi.RDS
similarity index 100%
rename from data/simdata_large_epi.RDS
rename to inst/extdata/simdata_large_epi.RDS
diff --git a/data/simdata_sw_dec.RDS b/inst/extdata/simdata_sw_dec.RDS
similarity index 100%
rename from data/simdata_sw_dec.RDS
rename to inst/extdata/simdata_sw_dec.RDS
diff --git a/data/veev2012.RDS b/inst/extdata/veev2012.RDS
similarity index 100%
rename from data/veev2012.RDS
rename to inst/extdata/veev2012.RDS
diff --git a/inst/include/stan_meta_header.hpp b/inst/include/stan_meta_header.hpp
new file mode 100644
index 00000000..3b914da2
--- /dev/null
+++ b/inst/include/stan_meta_header.hpp
@@ -0,0 +1 @@
+// Insert all #include<foo.hpp> statements here
diff --git a/inst/extdata/stanmodels/constant.stan b/inst/stan/constant.stan
similarity index 100%
rename from inst/extdata/stanmodels/constant.stan
rename to inst/stan/constant.stan
diff --git a/inst/extdata/stanmodels/tv_normal.stan b/inst/stan/tv_normal.stan
similarity index 99%
rename from inst/extdata/stanmodels/tv_normal.stan
rename to inst/stan/tv_normal.stan
index 9e679455..5e7a80d4 100644
--- a/inst/extdata/stanmodels/tv_normal.stan
+++ b/inst/stan/tv_normal.stan
@@ -42,4 +42,4 @@ generated quantities{
     P_sim[i] = Npos_sim[i] / Ntotal[i];
     logLikelihood[i] = binomial_lpmf(Npos[i] | Ntotal[i], P[i]);
   }
-}
\ No newline at end of file
+}
diff --git a/inst/extdata/stanmodels/tv_normal_log.stan b/inst/stan/tv_normal_log.stan
similarity index 99%
rename from inst/extdata/stanmodels/tv_normal_log.stan
rename to inst/stan/tv_normal_log.stan
index 27219287..c65884e3 100644
--- a/inst/extdata/stanmodels/tv_normal_log.stan
+++ b/inst/stan/tv_normal_log.stan
@@ -48,4 +48,4 @@ generated quantities{
     P_sim[i] = Npos_sim[i] / Ntotal[i];
     logLikelihood[i] = binomial_lpmf(Npos[i] | Ntotal[i], P[i]);
   }
-}
\ No newline at end of file
+}
diff --git a/man/extract_seromodel_summary.Rd b/man/extract_seromodel_summary.Rd
index ff9e916e..cd037799 100644
--- a/man/extract_seromodel_summary.Rd
+++ b/man/extract_seromodel_summary.Rd
@@ -7,13 +7,13 @@
 extract_seromodel_summary(seromodel_object)
 }
 \arguments{
-\item{seromodel_object}{\code{seromodel_object}. An object containing relevant information about the implementation of the model. 
+\item{seromodel_object}{\code{seromodel_object}. An object containing relevant information about the implementation of the model.
 Refer to \link{fit_seromodel} for further details.}
 }
 \value{
 \code{model_summary}. Object with a summary of \code{seromodel_object} containing the following:
 \tabular{ll}{
-\code{foi_model} \tab Name of the selected model. For further details refer to \link{save_or_load_model}. \cr \tab \cr
+\code{foi_model} \tab Name of the selected model. \cr \tab \cr
 \code{data_set} \tab Seroprevalence survey label.\cr \tab \cr
 \code{country} \tab Name of the country were the survey was conducted in. \cr \tab \cr
 \code{year} \tab Year in which the survey was conducted. \cr \tab \cr
@@ -28,9 +28,9 @@ Refer to \link{fit_seromodel} for further details.}
 }
 }
 \description{
-This method extracts a summary corresponding to a serological model object that contains information about the original serological 
-survey data used to fit the model, such as the year when the survey took place, the type of test taken and the corresponding antibody, 
-as well as information about the convergence of the model, like the expected log pointwise predictive density \code{elpd} and its 
+This method extracts a summary corresponding to a serological model object that contains information about the original serological
+survey data used to fit the model, such as the year when the survey took place, the type of test taken and the corresponding antibody,
+as well as information about the convergence of the model, like the expected log pointwise predictive density \code{elpd} and its
 corresponding standar deviation.
 }
 \examples{
diff --git a/man/fit_seromodel.Rd b/man/fit_seromodel.Rd
index 61fc06d8..67908578 100644
--- a/man/fit_seromodel.Rd
+++ b/man/fit_seromodel.Rd
@@ -58,9 +58,9 @@ This object is used as an input for the \link[rstan]{sampling} function \cr \tab
 }
 }
 \description{
-This function fits the specified model \code{foi_model} to the serological survey data \code{serodata} 
-by means of the \link[rstan]{sampling} method. The function determines whether the corresponding stan model 
-object needs to be compiled by means of the function \link{save_or_load_model}.
+This function fits the specified model \code{foi_model} to the serological survey data \code{serodata}
+by means of the \link[rstan]{sampling} method. The function determines whether the corresponding stan model
+object needs to be compiled by rstan.
 }
 \examples{
 \dontrun{
diff --git a/man/get_exposure_ages.Rd b/man/get_exposure_ages.Rd
index 96125fd3..aec51f06 100644
--- a/man/get_exposure_ages.Rd
+++ b/man/get_exposure_ages.Rd
@@ -13,7 +13,7 @@ get_exposure_ages(serodata)
 \code{exposure_ages}. An atomic vector with the numeration of the exposition years in serodata
 }
 \description{
-This function generates an atomic vector containing the exposition years corresponding to the specified serological survey data \code{serodata}. 
+This function generates an atomic vector containing the exposition years corresponding to the specified serological survey data \code{serodata}.
 The exposition years to the disease for each individual corresponds to the time from birth to the moment of the survey.
 }
 \examples{
diff --git a/man/get_prev_expanded.Rd b/man/get_prev_expanded.Rd
index 3bd09927..a95f4899 100644
--- a/man/get_prev_expanded.Rd
+++ b/man/get_prev_expanded.Rd
@@ -2,15 +2,17 @@
 % Please edit documentation in R/modelling.R
 \name{get_prev_expanded}
 \alias{get_prev_expanded}
-\title{Function that generates an object containing the confidence interval based on a 
+\title{Function that generates an object containing the confidence interval based on a
 Force-of-Infection fitting}
 \usage{
-get_prev_expanded(foi, serodata)
+get_prev_expanded(foi, serodata, bin_data = FALSE)
 }
 \arguments{
 \item{foi}{Object containing the information of the force of infection. It is obtained from \code{rstan::extract(seromodel_object$fit, "foi", inc_warmup = FALSE)[[1]]}.}
 
 \item{serodata}{A data frame containing the data from a seroprevalence survey. For further details refer to \link{run_seromodel}.}
+
+\item{bin_data}{TBD}
 }
 \value{
 \code{prev_final}. The expanded prevalence data. This is used for plotting purposes in the \code{visualization} module.
@@ -23,7 +25,7 @@ of the Force-of-Infection \code{foi} for plotting an analysis purposes.
 \dontrun{
 serodata <- prepare_serodata(serodata)
 seromodel_object <- run_seromodel(serodata = serodata,
-                                  foi_model = "constant")
+                          foi_model = "constant")
 foi <- rstan::extract(seromodel_object$fit, "foi")[[1]]
 get_prev_expanded <- function(foi, serodata)
 }
diff --git a/man/plot_foi.Rd b/man/plot_foi.Rd
index 82f7d417..6e376c70 100644
--- a/man/plot_foi.Rd
+++ b/man/plot_foi.Rd
@@ -1,5 +1,5 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/visualization.R
+% Please edit documentation in R/visualisation.R
 \name{plot_foi}
 \alias{plot_foi}
 \title{Function that generates a Force-of-Infection plot corresponding to the specified fitted serological model}
@@ -9,7 +9,11 @@ plot_foi(seromodel_object, max_lambda = NA, size_text = 25, foi_sim = NULL)
 \arguments{
 \item{seromodel_object}{Object containing the results of fitting a model by means of \link{run_seromodel}.}
 
+\item{max_lambda}{TBD}
+
 \item{size_text}{Text size use in the theme of the graph returned by the function.}
+
+\item{foi_sim}{TBD}
 }
 \value{
 A ggplot2 object containing the Force-of-infection-vs-time including the corresponding confidence interval.
diff --git a/man/plot_info_table.Rd b/man/plot_info_table.Rd
index 27145fca..72ddc003 100644
--- a/man/plot_info_table.Rd
+++ b/man/plot_info_table.Rd
@@ -1,5 +1,5 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/visualization.R
+% Please edit documentation in R/visualisation.R
 \name{plot_info_table}
 \alias{plot_info_table}
 \title{Function that generates a plot for a given table}
@@ -12,7 +12,7 @@ plot_info_table(info, size_text)
 \item{size_text}{Text size of the graph returned by the function}
 }
 \value{
-p, a variable that will be used in the \link{visualisation} module
+p the plot for the given table
 }
 \description{
 Function that generates a plot for a given table
diff --git a/man/plot_rhats.Rd b/man/plot_rhats.Rd
index 1840cb91..02b73f00 100644
--- a/man/plot_rhats.Rd
+++ b/man/plot_rhats.Rd
@@ -1,5 +1,5 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/visualization.R
+% Please edit documentation in R/visualisation.R
 \name{plot_rhats}
 \alias{plot_rhats}
 \title{Function that generates a plot of the R-hat estimates of the specified fitted serological model}
diff --git a/man/plot_seromodel.Rd b/man/plot_seromodel.Rd
index 78e10010..6a36f736 100644
--- a/man/plot_seromodel.Rd
+++ b/man/plot_seromodel.Rd
@@ -1,5 +1,5 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/visualization.R
+% Please edit documentation in R/visualisation.R
 \name{plot_seromodel}
 \alias{plot_seromodel}
 \title{Function that generates a vertical arrange of plots showing a summary of a model, the estimated seroprevalence,
@@ -15,7 +15,11 @@ plot_seromodel(
 \arguments{
 \item{seromodel_object}{Object containing the results of fitting a model by means of \link{run_seromodel}.}
 
+\item{max_lambda}{TBD}
+
 \item{size_text}{Text size use in the theme of the graph returned by the function.}
+
+\item{foi_sim}{TBD}
 }
 \value{
 A ggplot object with a vertical arrange containing the seropositivity, force of infection, and convergence plots.
diff --git a/man/plot_seroprev.Rd b/man/plot_seroprev.Rd
index 7a5de2d0..8948decb 100644
--- a/man/plot_seroprev.Rd
+++ b/man/plot_seroprev.Rd
@@ -1,5 +1,5 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/visualization.R
+% Please edit documentation in R/visualisation.R
 \name{plot_seroprev}
 \alias{plot_seroprev}
 \title{Function that generates the sero-positivity plot from a raw serological survey dataset}
diff --git a/man/plot_seroprev_fitted.Rd b/man/plot_seroprev_fitted.Rd
index 6c8d4366..18740584 100644
--- a/man/plot_seroprev_fitted.Rd
+++ b/man/plot_seroprev_fitted.Rd
@@ -1,5 +1,5 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/visualization.R
+% Please edit documentation in R/visualisation.R
 \name{plot_seroprev_fitted}
 \alias{plot_seroprev_fitted}
 \title{Function that generates a seropositivity plot corresponding to the specified fitted serological model}
diff --git a/man/prepare_serodata.Rd b/man/prepare_serodata.Rd
index 5b03d6e2..33b4a1de 100644
--- a/man/prepare_serodata.Rd
+++ b/man/prepare_serodata.Rd
@@ -21,7 +21,9 @@ This data frame must contain the following columns:
 \code{antibody} \tab antibody \cr \tab \cr
 }}
 
-\item{alpha}{probability of a type I error. For further details refer to \link{Hmisc::binconf}.}
+\item{alpha}{probability of a type I error. For further details refer to \link[Hmisc]{binconf}.}
+
+\item{add_age_mean_f}{TBD}
 }
 \value{
 serodata with additional columns necessary for the analysis. These columns are:
diff --git a/man/run_seromodel.Rd b/man/run_seromodel.Rd
index 12d78585..becb135a 100644
--- a/man/run_seromodel.Rd
+++ b/man/run_seromodel.Rd
@@ -55,12 +55,14 @@ For further details refer to the \code{control} parameter in \link[rstan]{sampli
 \item{m_treed}{Maximum tree depth for the binary tree used in the NUTS stan sampler. For further details refer to the \code{control} parameter in \link[rstan]{sampling}.}
 
 \item{decades}{Number of decades covered by the survey data.}
+
+\item{print_summary}{TBD}
 }
 \value{
 \code{seromodel_object}. An object containing relevant information about the implementation of the model. For further details refer to \link{fit_seromodel}.
 }
 \description{
-This function runs the specified model for the Force-of-Infection \code{foi_model} using the data froma seroprevalence survey 
+This function runs the specified model for the Force-of-Infection \code{foi_model} using the data froma seroprevalence survey
 \code{serodata} as the input data. See \link{fit_seromodel} for further details.
 }
 \examples{
diff --git a/man/save_or_load_model.Rd b/man/save_or_load_model.Rd
deleted file mode 100644
index 762da384..00000000
--- a/man/save_or_load_model.Rd
+++ /dev/null
@@ -1,30 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/modelling.R
-\name{save_or_load_model}
-\alias{save_or_load_model}
-\title{Function used to determine whether the stan model corresponding to the specified serological model has been already compiled or not}
-\usage{
-save_or_load_model(foi_model = "constant")
-}
-\arguments{
-\item{foi_model}{Name of the selected model. Current version provides three options:
-\describe{
-\item{\code{"constant"}}{Runs a constant model}
-\item{\code{"tv_normal"}}{Runs a normal model}
-\item{\code{"tv_normal_log"}}{Runs a normal logarithmic model}
-}}
-}
-\value{
-\code{model}. The rstan model object corresponding to the selected model.
-}
-\description{
-This function determines whether the corresponding .RDS file of the selected model exists or not.
-In case the .RDS file exists, it is read and returned; otherwise, the object model is created through the
-\link[rstan]{stan_model} function, saved as an .RDS file and returned as the output of the function.
-}
-\examples{
-\dontrun{
-model <- save_or_load_model(foi_model="constant")
-}
-
-}
diff --git a/man/serofoi-package.Rd b/man/serofoi-package.Rd
index 6c1380ea..1ecdd986 100644
--- a/man/serofoi-package.Rd
+++ b/man/serofoi-package.Rd
@@ -1,29 +1,13 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/serofoi_package.R
+% Please edit documentation in R/serofoi-package.R
 \docType{package}
 \name{serofoi-package}
-\alias{serofoi}
 \alias{serofoi-package}
-\title{serofoi: Estimates the Force-of-Infection of a given pathogen from population based sero-prevalence studies}
+\alias{serofoi}
+\title{The 'serofoi' package.}
 \description{
-R package to estimate time-varying Force-of-Infection of a given pathogen from population based sero-prevalence studies using a bayesian framework.
-}
-\seealso{
-Useful links:
-\itemize{
-  \item \url{https://trace-lac.github.io/serofoi/}
-}
-
-}
-\author{
-\strong{Maintainer}: Zulma M. Cucunubá \email{zulma.cucunuba@javeriana.edu.co} (\href{https://orcid.org/0000-0002-8165-3198}{ORCID})
-
-Authors:
-\itemize{
-  \item Nicolás Torres
-  \item Ben Lambert
-  \item Pierre Nouvellet
+A DESCRIPTION OF THE PACKAGE
 }
-
+\references{
+Stan Development Team (NA). RStan: the R interface to Stan. R package version 2.26.22. https://mc-stan.org
 }
-\keyword{internal}
diff --git a/tests/testthat/extdata/haiti_ssa_sample.RDS b/tests/testthat/extdata/haiti_ssa_sample.RDS
new file mode 100644
index 00000000..04f9fa75
Binary files /dev/null and b/tests/testthat/extdata/haiti_ssa_sample.RDS differ
diff --git a/tests/testthat/extdata/prev_expanded_constant.RDS b/tests/testthat/extdata/prev_expanded_constant.RDS
new file mode 100644
index 00000000..02286383
Binary files /dev/null and b/tests/testthat/extdata/prev_expanded_constant.RDS differ
diff --git a/tests/testthat/extdata/prev_expanded_tv_normal.RDS b/tests/testthat/extdata/prev_expanded_tv_normal.RDS
new file mode 100644
index 00000000..41948d0b
Binary files /dev/null and b/tests/testthat/extdata/prev_expanded_tv_normal.RDS differ
diff --git a/tests/testthat/extdata/prev_expanded_tv_normal_log.RDS b/tests/testthat/extdata/prev_expanded_tv_normal_log.RDS
new file mode 100644
index 00000000..71b7344b
Binary files /dev/null and b/tests/testthat/extdata/prev_expanded_tv_normal_log.RDS differ
diff --git a/tests/testthat/extdata/simdata_foiD_n05_group.RDS b/tests/testthat/extdata/simdata_foiD_n05_group.RDS
deleted file mode 100644
index 14e493e9..00000000
Binary files a/tests/testthat/extdata/simdata_foiD_n05_group.RDS and /dev/null differ
diff --git a/tests/testthat/models_serialization.R b/tests/testthat/models_serialization.R
new file mode 100644
index 00000000..6306fac0
--- /dev/null
+++ b/tests/testthat/models_serialization.R
@@ -0,0 +1,33 @@
+library(devtools)
+library(dplyr)
+library(serofoi)
+library(testthat)
+
+set.seed(1234) # For reproducibility
+
+#----- Read and prepare data
+data("simdata_large_epi")
+simdata <- simdata_large_epi %>% 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)
+
+model_constant_json <- jsonlite::serializeJSON(models_list[[1]])
+write_json(model_constant_json, testthat::test_path("extdata", "model_constant.json"))
+
+model_tv_normal_json <- jsonlite::serializeJSON(models_list[[2]])
+write_json(model_tv_normal_json, testthat::test_path("extdata", "model_tv_normal.json"))
+
+model_tv_normal_log_json <- jsonlite::serializeJSON(models_list[[3]])
+write_json(model_tv_normal_json, testthat::test_path("extdata", "model_tv_normal_log.json"))
diff --git a/tests/testthat/test_constant_model.R b/tests/testthat/test_constant_model.R
deleted file mode 100644
index 19883d17..00000000
--- a/tests/testthat/test_constant_model.R
+++ /dev/null
@@ -1,35 +0,0 @@
-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")
-  data_test <- serodata %>% prepare_serodata(alpha = 0.05)
-
-  #----- Generate plots for the constant model
-
-  model_name <- "constant"
-  model <- run_seromodel(serodata = data_test,
-                         foi_model = model_name,
-                         n_iters = 1000)
-  model_plot <- plot_seromodel(model, size_text = 6)
-  vdiffr::expect_doppelganger(paste0(model_name, "_model_plot"), model_plot)
-
-  model_seroprev_plot <- plot_seroprev_fitted(model, size_text = 15)
-  vdiffr::expect_doppelganger(paste0(model_name, "_sp_fitted_plot"), model_seroprev_plot)
-
-  model_foi_plot <- plot_foi(model, size_text = 15)
-  vdiffr::expect_doppelganger(paste0(model_name, "_foi_plot"), model_foi_plot)
-
-  model_rhats_plot <- plot_rhats(model, size_text = 15)
-  vdiffr::expect_doppelganger(paste0(model_name, "_rhats_plot"), model_rhats_plot)
-
-})
diff --git a/tests/testthat/test_issue_47.R b/tests/testthat/test_issue_47.R
new file mode 100644
index 00000000..e46b2740
--- /dev/null
+++ b/tests/testthat/test_issue_47.R
@@ -0,0 +1,21 @@
+test_that("issue 47", {
+  skip_on_os(c("windows", "mac"))
+  source("testing_utils.R")
+
+  library(devtools)
+  library(dplyr)
+
+  # Load data
+  ## This dataset is already prepared
+  data_path <- testthat::test_path("extdata", "haiti_ssa_sample.RDS")
+  data_issue <- readRDS(data_path)
+
+  # Error reproduction
+  model_test <- run_seromodel(data_issue, foi_model = "tv_normal", print_summary = FALSE)
+  foi <- rstan::extract(model_test$fit, "foi", inc_warmup = FALSE)[[1]]
+  age_max <- max(data_issue$age_mean_f)
+  prev_expanded <- get_prev_expanded(foi, serodata = data_issue)
+
+  # Test
+  expect_length(prev_expanded$age, n = age_max)
+})
diff --git a/tests/testthat/test_modelling.R b/tests/testthat/test_modelling.R
new file mode 100644
index 00000000..76494b19
--- /dev/null
+++ b/tests/testthat/test_modelling.R
@@ -0,0 +1,61 @@
+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"))
+  source("testing_utils.R")
+  set.seed(1234) # For reproducibility
+
+  library(devtools)
+  library(dplyr)
+  library(vdiffr)
+
+  #----- Read and prepare data
+  data("serodata")
+  data_test <- serodata %>% prepare_serodata(alpha = 0.05)
+
+  data_constant_path <- testthat::test_path("extdata", "prev_expanded_constant.RDS")
+  data_tv_normal_path <- testthat::test_path("extdata", "prev_expanded_tv_normal.RDS")
+  data_tv_normal_log_path <- testthat::test_path("extdata", "prev_expanded_tv_normal_log.RDS")
+
+  prev_expanded_tv_normal_log <- readRDS(data_constant_path)
+
+  #----- Test for the constant model
+
+  model_name <- "constant"
+  model_object <- run_seromodel(serodata = data_test,
+                                foi_model = model_name,
+                                n_iters = 1000,
+                                print_summary = FALSE)
+
+  foi <- rstan::extract(model_object$fit, "foi", inc_warmup = FALSE)[[1]]
+  prev_expanded <- get_prev_expanded(foi, serodata = model_object$serodata)
+  prev_expanded_constant <- readRDS(data_constant_path)
+
+  testthat::expect_equal(prev_expanded, prev_expanded_constant, tolerance = TRUE)
+
+  #----- Test for the tv_normal model
+
+  model_name <- "tv_normal"
+  model_object <- run_seromodel(serodata = data_test,
+                                foi_model = model_name,
+                                n_iters = 1000)
+
+  foi <- rstan::extract(model_object$fit, "foi", inc_warmup = FALSE)[[1]]
+  prev_expanded <- get_prev_expanded(foi, serodata = model_object$serodata)
+  prev_expanded_tv_normal <- readRDS(data_tv_normal_path)
+  testthat::expect_equal(prev_expanded, prev_expanded_tv_normal, tolerance = TRUE)
+
+  #----- Test for the tv_normal_log model
+
+  model_name <- "tv_normal_log"
+  model_object <- run_seromodel(serodata = data_test,
+                                foi_model = model_name,
+                                n_iters = 1000)
+
+  foi <- rstan::extract(model_object$fit, "foi", inc_warmup = FALSE)[[1]]
+  prev_expanded <- get_prev_expanded(foi, serodata = model_object$serodata)
+  prev_expanded_tv_normal <- readRDS(data_tv_normal_path)
+  testthat::expect_equal(prev_expanded, prev_expanded_tv_normal_log, tolerance = TRUE)
+
+})
diff --git a/tests/testthat/test_plot_foi.R b/tests/testthat/test_plot_foi.R
deleted file mode 100644
index ca3146b3..00000000
--- a/tests/testthat/test_plot_foi.R
+++ /dev/null
@@ -1,51 +0,0 @@
-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)
-
-})
diff --git a/tests/testthat/test_tv_normal_log_model.R b/tests/testthat/test_tv_normal_log_model.R
deleted file mode 100644
index 466485d2..00000000
--- a/tests/testthat/test_tv_normal_log_model.R
+++ /dev/null
@@ -1,39 +0,0 @@
-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")
-  data_test <- serodata %>% prepare_serodata(alpha = 0.05)
-
-  #----- Plot raw data
-  data_test_plot <- plot_seroprev(data_test, size_text = 15)
-  vdiffr::expect_doppelganger("serodata_plot", data_test_plot)
-
-  #----- Generate plots for the constant model
-
-  model_name <- "tv_normal_log"
-  model <- run_seromodel(serodata = data_test,
-                         foi_model = model_name,
-                         n_iters = 1000)
-  model_plot <- plot_seromodel(model, size_text = 6)
-  vdiffr::expect_doppelganger(paste0(model_name, "_model_plot"), model_plot)
-
-  model_seroprev_plot <- plot_seroprev_fitted(model, size_text = 15)
-  vdiffr::expect_doppelganger(paste0(model_name, "_sp_fitted_plot"), model_seroprev_plot)
-
-  model_foi_plot <- plot_foi(model, size_text = 15)
-  vdiffr::expect_doppelganger(paste0(model_name, "_foi_plot"), model_foi_plot)
-
-  model_rhats_plot <- plot_rhats(model, size_text = 15)
-  vdiffr::expect_doppelganger(paste0(model_name, "_rhats_plot"), model_rhats_plot)
-
-})
diff --git a/tests/testthat/test_tv_normal_model.R b/tests/testthat/test_tv_normal_model.R
deleted file mode 100644
index 0e5d94c2..00000000
--- a/tests/testthat/test_tv_normal_model.R
+++ /dev/null
@@ -1,39 +0,0 @@
-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")
-  data_test <- serodata %>% prepare_serodata(alpha = 0.05)
-
-  #----- Plot raw data
-  data_test_plot <- plot_seroprev(data_test, size_text = 15)
-  vdiffr::expect_doppelganger("serodata_plot", data_test_plot)
-
-  #----- Generate plots for the constant model
-
-  model_name <- "tv_normal"
-  model <- run_seromodel(serodata = data_test,
-                         foi_model = model_name,
-                         n_iters = 1000)
-  model_plot <- plot_seromodel(model, size_text = 6)
-  vdiffr::expect_doppelganger(paste0(model_name, "_model_plot"), model_plot)
-
-  model_seroprev_plot <- plot_seroprev_fitted(model, size_text = 15)
-  vdiffr::expect_doppelganger(paste0(model_name, "_sp_fitted_plot"), model_seroprev_plot)
-
-  model_foi_plot <- plot_foi(model, size_text = 15)
-  vdiffr::expect_doppelganger(paste0(model_name, "_foi_plot"), model_foi_plot)
-
-  model_rhats_plot <- plot_rhats(model, size_text = 15)
-  vdiffr::expect_doppelganger(paste0(model_name, "_rhats_plot"), model_rhats_plot)
-
-})
diff --git a/tests/testthat/test_visualisation.R b/tests/testthat/test_visualisation.R
new file mode 100644
index 00000000..52ef09fa
--- /dev/null
+++ b/tests/testthat/test_visualisation.R
@@ -0,0 +1,146 @@
+# Test for the function plot_seroprev_fitted
+
+library(testthat)
+
+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"))
+  source("testing_utils.R")
+  set.seed(1234) # For reproducibility
+
+  library(devtools)
+  library(dplyr)
+  library(vdiffr)
+  library(jsonlite)
+
+  data("simdata_large_epi")
+  simdata <- simdata_large_epi %>% 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
+
+
+  #----- Results visualisation
+  size_text <- 6
+  max_lambda <- 1.55
+
+  model_constant_json <- jsonlite::fromJSON(testthat::test_path("extdata", "model_constant.json"))
+  model_constant <- jsonlite::unserializeJSON(model_constant_json)
+  constant_plot <- plot_seromodel(model_constant,
+    size_text = size_text,
+    max_lambda = max_lambda,
+    foi_sim = foi_sim
+  )
+
+  model_tv_normal_json <- fromJSON(testthat::test_path("extdata", "model_tv_normal.json"))
+  model_tv_normal <- jsonlite::unserializeJSON(model_tv_normal_json)
+  tv_normal_plot <- plot_seromodel(model_tv_normal,
+    size_text = size_text,
+    max_lambda = max_lambda,
+    foi_sim = foi_sim
+  )
+
+  model_tv_normal_log_json <- fromJSON(testthat::test_path("extdata", "model_tv_normal_log.json"))
+  model_tv_normal_log <- jsonlite::unserializeJSON(model_tv_normal_log_json)
+  tv_normal_log_plot <- plot_seromodel(model_tv_normal_log,
+    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)
+})
+
+
+# #Test for the function plot_rhats
+#
+# library(testthat)
+#
+# # Define a test context
+# context("Testing plot_rhats function")
+#
+# # Create a mock seromodel_object
+# mock_seromodel_object <- list(fit = "model did not run")
+#
+# # Define a test case for the else statement
+# test_that("plot_rhats function works for else statement", {
+#
+#   # Call the function with the mock seromodel_object
+#   rhats_plot <- plot_rhats(mock_seromodel_object)
+#
+#   # Expect the output to be a ggplot object
+#   expect_is(rhats_plot, "ggplot")
+#
+#   # Expect the plot to have a single point
+#   expect_equal(length(rhats_plot$layers[[1]]$data), 1)
+#
+#   # Expect the plot to have a single label
+#   expect_equal(length(rhats_plot$layers[[2]]$data), 1)
+#
+#   # Expect the label to be "errors"
+#   expect_equal(rhats_plot$layers[[2]]$label, "errors")
+# })
+#
+#
+# #Test for the function plot_seromodel
+#
+# library(testthat)
+#
+# # Test for exception in else
+# test_that("plot_seromodel prints an error message and returns an empty plot object when a model cannot be fitted", {
+#   # Create a seromodel object with fit as a feature
+#   seromodel_object <- list(fit = "no_fit", model = "my_model")
+#   # Run the function and check that it returns an empty plot object
+#   expect_silent(plot_seromodel(seromodel_object))
+#   expect_equal(length(plot_seromodel(seromodel_object)$grobs), 5)
+# })
+#
+# # We create a helper function that returns an unwrapped seromodel object
+# create_dummy_seromodel <- function() {
+#   # empty object
+#   seromodel <- list()
+#   seromodel$fit <- "dummy"
+#   seromodel$serodata <- data.frame(age = c(0, 10, 20, 30, 40, 50, 60),
+#                                    p_obs_bin = c(0.2, 0.4, 0.6, 0.8, 0.9, 0.95, 0.99),
+#                                    bin_size = c(5, 10, 20, 30, 40, 50, 60))
+#   return(seromodel)
+# }
+#
+# # Unit tests for plot_seroprev_fitted()
+# test_that("plot_seroprev_fitted() returns an empty plot for an unfitted seromodel object", {
+#   # We create an unadjusted seromodel object
+#   seromodel <- create_dummy_seromodel()
+#   # We call the function plot_seroprev_fitted()
+#   plot <- plot_seroprev_fitted(seromodel)
+#   # We verify that the plot object is an empty ggplot
+#   expect_true(class(plot) == "ggplot")
+#   expect_true(length(plot$layers) == 0)
+# })
+#
+#
+# #Test for the function plot_foi
+#
+# library(testthat)
+#
+# # Define the test context
+# context("Test of the function plot_foi")
+#
+# # Create a test to check the else block
+# test_that("The plot_foi function should output an empty plot when the model is not running", {
+#
+#   # Create an empty object that simulates the output of the model that failed
+#   empty_model <- list(fit = "failure")
+#
+#   # Run the plot_foi function with the empty model
+#   plot <- plot_foi(empty_model)
+#
+#   # Check if the output is an empty graph
+#   expect_identical(ggplot2::ggplot(), plot)
+# })
+#
diff --git a/vignettes/foi_models.Rmd b/vignettes/foi_models.Rmd
index 5f9a82ad..2c53d6a9 100644
--- a/vignettes/foi_models.Rmd
+++ b/vignettes/foi_models.Rmd
@@ -4,7 +4,7 @@ output: rmarkdown::html_vignette
 bibliography: references.bib
 link-citations: true
 vignette: >
-  %\VignetteIndexEntry{FoI Models}
+  %\VignetteIndexEntry{An introduction to Force-of-Infection (FoI) Models}
   %\VignetteEngine{knitr::rmarkdown}
   %\VignetteEncoding{UTF-8}
 ---
diff --git a/vignettes/serofoi.Rmd b/vignettes/serofoi.Rmd
index 95158746..082e86ef 100644
--- a/vignettes/serofoi.Rmd
+++ b/vignettes/serofoi.Rmd
@@ -2,7 +2,7 @@
 title: "An introduction to serofoi"
 output: rmarkdown::html_vignette
 vignette: >
-  %\VignetteIndexEntry{serofoi}
+  %\VignetteIndexEntry{An introduction to serofoi}
   %\VignetteEngine{knitr::rmarkdown}
   %\VignetteEncoding{UTF-8}
 ---
diff --git a/vignettes/use_cases.Rmd b/vignettes/use_cases.Rmd
index a0813c89..361fd1f2 100644
--- a/vignettes/use_cases.Rmd
+++ b/vignettes/use_cases.Rmd
@@ -1,10 +1,10 @@
 ---
-title: "Real-life use cases for _serofoi_"
+title: "Real-life use cases for serofoi"
 output: rmarkdown::html_vignette
 bibliography: references.bib
 link-citations: true
 vignette: >
-  %\VignetteIndexEntry{FoI Models}
+  %\VignetteIndexEntry{Real-life use cases for serofoi}
   %\VignetteEngine{knitr::rmarkdown}
   %\VignetteEncoding{UTF-8}
 ---