diff --git a/.github/CONTRIBUTING.md b/.github/CONTRIBUTING.md
new file mode 100644
index 00000000..2e5cdb42
--- /dev/null
+++ b/.github/CONTRIBUTING.md
@@ -0,0 +1,40 @@
+# Contributing to serofoi
+
+This outlines how to propose a change to serofoi.
+
+## Making changes
+
+If you want to make a change, it's a good idea to first file an issue and make sure someone from the team agrees that it’s needed.
+If you’ve found a bug, please file an issue that illustrates the bug with a minimal 
+[reprex](https://www.tidyverse.org/help/#reprex) (this will also help you write a unit test, if needed). See [bug report template](../.github/ISSUE_TEMPLATE/bug_report.md). If you have a feature request see [feature request](../.github/ISSUE_TEMPLATE/feature_request.md).
+
+### Pull request process
+
+See [pull request template](../.github/PULL_REQUEST_TEMPLATE/pull_request_template.md)
+
+*   Fork the package and clone onto your computer. If you haven't done this before, we recommend using `usethis::create_from_github("epiverse-trace/serofoi", fork = TRUE)`.
+
+*   Install all development dependencies with `devtools::install_dev_deps()`, and then make sure the package passes R CMD check by running `devtools::check()`. 
+    If R CMD check doesn't pass cleanly, it's a good idea to ask for help before continuing. 
+*   Create a Git branch for your pull request (PR). We recommend using `usethis::pr_init("brief-description-of-change")`.
+
+*   Make your changes, commit to git, and then create a PR by running `usethis::pr_push()`, and following the prompts in your browser.
+    The title of your PR should briefly describe the change.
+    The body of your PR should contain `Fixes #issue-number`.
+
+*  For user-facing changes, add a bullet to the top of `NEWS.md` (i.e. just below the first header). Follow the style described in <https://style.tidyverse.org/news.html>.
+
+### Code style
+
+*   New code should follow the tidyverse [style guide](https://style.tidyverse.org). 
+    You can use the [styler](https://CRAN.R-project.org/package=styler) package to apply these styles, but please don't restyle code that has nothing to do with your PR.  
+
+*  We use [roxygen2](https://cran.r-project.org/package=roxygen2), with [Markdown syntax](https://cran.r-project.org/web/packages/roxygen2/vignettes/rd-formatting.html), for documentation.  
+
+*  We use [testthat](https://cran.r-project.org/package=testthat) for unit tests. 
+   Contributions with test cases included are easier to accept.  
+
+## Code of Conduct
+
+Please note that the serofoi project is released with a
+[Contributor Code of Conduct](https://github.com/epiverse-trace/.github/blob/main/CODE_OF_CONDUCT.md). By contributing to this project you agree to abide by its terms.
diff --git a/.github/workflows/pkgdown.yaml b/.github/workflows/pkgdown.yaml
new file mode 100644
index 00000000..087f0b05
--- /dev/null
+++ b/.github/workflows/pkgdown.yaml
@@ -0,0 +1,46 @@
+# Workflow derived from https://github.com/r-lib/actions/tree/v2/examples
+# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help
+on:
+  push:
+    branches: [main, master]
+  pull_request:
+    branches: [main, master]
+  release:
+    types: [published]
+  workflow_dispatch:
+
+name: pkgdown
+
+jobs:
+  pkgdown:
+    runs-on: ubuntu-latest
+    # Only restrict concurrency for non-PR jobs
+    concurrency:
+      group: pkgdown-${{ github.event_name != 'pull_request' || github.run_id }}
+    env:
+      GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
+    steps:
+      - uses: actions/checkout@v3
+
+      - uses: r-lib/actions/setup-pandoc@v2
+
+      - uses: r-lib/actions/setup-r@v2
+        with:
+          use-public-rspm: true
+
+      - uses: r-lib/actions/setup-r-dependencies@v2
+        with:
+          extra-packages: any::pkgdown, local::.
+          needs: website
+
+      - name: Build site
+        run: pkgdown::build_site_github_pages(new_process = FALSE, install = FALSE)
+        shell: Rscript {0}
+
+      - name: Deploy to GitHub pages 🚀
+        if: github.event_name != 'pull_request'
+        uses: JamesIves/github-pages-deploy-action@v4.4.1
+        with:
+          clean: false
+          branch: gh-pages
+          folder: docs
diff --git a/.gitignore b/.gitignore
index 2af26a0e..f319541d 100644
--- a/.gitignore
+++ b/.gitignore
@@ -54,3 +54,4 @@ tests/testthat/Rplots.pdf
 t,
 tests/testthat/extdata/plots/actual/*.png
 tests/testthat/extdata/dataframes/actual/*.csv
+docs
diff --git a/DESCRIPTION b/DESCRIPTION
index 395ff896..cae4f40e 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -3,15 +3,33 @@ Type: Package
 Title: Estimates the Force-of-Infection of a given pathogen from population based sero-prevalence studies
 Version: 0.1.0
 Authors@R: 
-  c(
-    person("Zulma M. Cucunubá", role = "aut"),
-    person("Ben Lambert", role = "aut"),
-    person("Pierre Nouvellet", role = "aut"),
-    person("N. T. Domínguez", role = c("aut", "cre"), email = "ex-ntorres@javeriana.edu.co")
-  )
+    c(
+        person(
+            given = "Zulma M.",
+            family = "Cucunubá",
+            role = c("aut", "cre"),
+            email = "zulma.cucunuba@javeriana.edu.co",
+            comment = c(ORCID = "https://orcid.org/0000-0002-8165-3198")
+        ),
+        person(
+            given = "Nicolás",
+            family = "Torres",
+            role = c("aut")
+        ),
+        person(
+            given = "Ben",
+            family = "Lambert",
+            role = c("aut")
+        ),
+        person(
+            given = "Pierre",
+            family = "Nouvellet",
+            role = c("aut")
+        )
+    )
 Contributor: Miguel Gámez, Geraldine Gómez, Jaime A. Pavlich-Mariscal
 Maintainer: The package maintainer <zulma.cucunuba@javeriana.edu.co>
-Description: R package to estimate time-varying Force-of-Infection of a given pathogen from population based sero-prevalence studies on a bayesian framework.
+Description: R package to estimate time-varying Force-of-Infection of a given pathogen from population based sero-prevalence studies using a bayesian framework.
 License: MIT + file LICENSE
 Encoding: UTF-8
 LazyData: true
@@ -43,5 +61,4 @@ Suggests:
     knitr,
     rmarkdown
 VignetteBuilder: knitr
-Additional_repositories:
-    https://mc-stan.org/r-packages/
+URL: https://trace-lac.github.io/serofoi/
\ No newline at end of file
diff --git a/R/chagas2012.R b/R/chagas2012.R
new file mode 100644
index 00000000..e58bc111
--- /dev/null
+++ b/R/chagas2012.R
@@ -0,0 +1,15 @@
+#' Seroprevalence data on serofoi
+#'
+#' Data from a serological surveys
+#'
+#' @docType data
+#'
+#' @usage chagas2012
+#'
+#' @format An object of class \code{"cross"}; see \code{\link[qtl]{read.cross}}.
+#'
+#' @keywords datasets
+#'
+#' @examples
+#' chagas2012
+"chagas2012"
\ No newline at end of file
diff --git a/data/data_test.R b/R/chik2015.R
similarity index 83%
rename from data/data_test.R
rename to R/chik2015.R
index 0104ef0a..e72f7357 100644
--- a/data/data_test.R
+++ b/R/chik2015.R
@@ -4,12 +4,12 @@
 #'
 #' @docType data
 #'
-#' @usage data_test
+#' @usage chik2015
 #'
 #' @format An object of class \code{"cross"}; see \code{\link[qtl]{read.cross}}.
 #'
 #' @keywords datasets
 #'
 #' @examples
-#' data_test
-"serofoi"
+#' chik2015
+"chik2015"
\ No newline at end of file
diff --git a/R/modelling.R b/R/modelling.R
index e35d1c44..481e7def 100644
--- a/R/modelling.R
+++ b/R/modelling.R
@@ -1,7 +1,6 @@
 # TODO For some reason, the examples cannot access the mydata variable
-#' Run model
+#' Run the specified stan model for the force-of-infection
 #'
-#' Runs the specified stan model for the force-of-infection
 #' @param model_data A data frame containing the data from a seroprevalence survey.
 #' This data frame must contain the following columns:
 #' \tabular{ll}{
@@ -71,7 +70,7 @@ run_model <- function(model_data,
 #' Save or load model
 #'
 #' 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 
+#' 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 model_name Name of the selected model. Current version provides three options:
 #' \describe{
diff --git a/R/mydata.R b/R/mydata.R
new file mode 100644
index 00000000..9757e8ee
--- /dev/null
+++ b/R/mydata.R
@@ -0,0 +1,15 @@
+#' Seroprevalence data on serofoi
+#'
+#' Data from a serological surveys
+#'
+#' @docType data
+#'
+#' @usage mydata
+#'
+#' @format An object of class \code{"cross"}; see \code{\link[qtl]{read.cross}}.
+#'
+#' @keywords datasets
+#'
+#' @examples
+#' mydata
+"mydata"
\ No newline at end of file
diff --git a/R/veev2012.R b/R/veev2012.R
new file mode 100644
index 00000000..a19088c0
--- /dev/null
+++ b/R/veev2012.R
@@ -0,0 +1,15 @@
+#' Seroprevalence data on serofoi
+#'
+#' Data from a serological surveys
+#'
+#' @docType data
+#'
+#' @usage veev2012
+#'
+#' @format An object of class \code{"cross"}; see \code{\link[qtl]{read.cross}}.
+#'
+#' @keywords datasets
+#'
+#' @examples
+#' veev2012
+"veev2012"
\ No newline at end of file
diff --git a/README.Rmd b/README.Rmd
index 9a6bd179..eac385d1 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -1,5 +1,10 @@
 ---
 output: github_document
+editor_options: 
+  markdown: 
+    wrap: 72
+bibliography: references.bib
+link-citations: true
 ---
 
 ```{r, include = FALSE}
@@ -11,105 +16,199 @@ knitr::opts_chunk$set(
 )
 ```
 
-## *serofoi* <img src="man/figures/serofoi-logo.png" align="right" width="120" />
+## *serofoi*: force-of-infection from population based serosurveys with age-disagregated data <img src="man/figures/serofoi-logo.png" align="right" width="130"/>
 
-An R package to estimates the *Force-of-Infection* of a given pathogen from population
-based sero-prevalence studies on a Bayesian framework.
+<!-- badges: start -->
 
+[![License:
+MIT](https://img.shields.io/badge/License-MIT-yellow.svg)](https://opensource.org/licenses/MIT)
+[![R-CMD-check](https://github.com/TRACE-LAC/serofoi/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/TRACE-LAC/serofoi/actions/workflows/R-CMD-check.yaml)
+[![Codecov test
+coverage](https://codecov.io/gh/TRACE-LAC/serofoi/branch/bugfixes-jaime/graph/badge.svg)](https://github.com/TRACE-LAC/serofoi?branch=bugfixes-jaime)
+[![lifecycle-concept](https://raw.githubusercontent.com/reconverse/reconverse.github.io/master/images/badge-concept.svg)](https://www.reconverse.org/lifecycle.html#concept)
 
-<!-- badges: start -->
-[![License: MIT](https://img.shields.io/badge/License-MIT-yellow.svg)](https://opensource.org/licenses/MIT)
-[![R-CMD-check](https://github.com/epiverse-trace/readepi/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/epiverse-trace/readepi/actions/workflows/R-CMD-check.yaml)
-[![Codecov test coverage](https://codecov.io/gh/epiverse-trace/readepi/branch/main/graph/badge.svg)](https://app.codecov.io/gh/epiverse-trace/readepi?branch=main)
-[![lifecycle-concept](https://raw.githubusercontent.com/reconverse/reconverse.github.io/master/images/badge-concept.svg)](https://www.reconverse.org/lifecycle.html#concept) 
 <!-- badges: end -->
 
+***serofoi*** is an R package to estimates the *Force-of-Infection* of a
+given pathogen from age-disaggregated population based sero-prevalence
+studies, using a Bayesian framework.
+
+***serofoi*** implements methods outlined in [@cucunubá2017] and
+[@carrera2020]
+
+***serofoi*** relies on
+[`rstan`](https://mc-stan.org/users/interfaces/rstan)
+
+***serofoi*** is part of the [Epiverse
+Initiative](https://data.org/initiatives/epiverse/).
 
 ## Installation
+
 You can install the **development version** of `serofoi` from
 [GitHub](https://github.com/) with:
 
 ``` r
 # install.packages("remotes")
-remotes::install_github("TRACE-LAC/serofoi")
+# remotes::install_github("TRACE-LAC/serofoi")
+library(serofoi)
 ```
 
 ## Quick start
 
-These examples illustrate some of the current functionalities:
+```{r cleaning, include = FALSE, echo = TRUE}
+library(serofoi)
+rownames(mydata) <- NULL 
 
-```{r example}
-##library(serofoi)
+```
 
-remotes::install_github("TRACE-LAC/serofoi", ref = "dev-zulma", force = TRUE)
+The package provides an example dataset of the observed serosurvey data,
+`mydata`. This example is the basic entry for the package.
 
-library(serofoi)
+```{r ex, include = TRUE}
 
-# Prepare the data for using serofoi
+# Load example mydata data included with the package
+data(mydata)
+head(mydata, 5)
+
+```
+
+The function `prepare_data` will prepare the entry data for entering the
+modelling functions. The seroprevalence *prepared data* can be
+visualised with the `plot_seroprev` function. This function also plots
+the binomial confidence interval of the observed data.
+
+```{r data_test, include = TRUE, fig.height=6, fig.width=6, message=FALSE}
 
 data_test <- prepare_data(mydata)
 
+plot_seroprev(data_test, size_text = 15)
+
+```
+
+#### Current version of the package runs ***three*** different FoI models
+
+The `run_model` function allows specifying the Bayesian model from *R*,
+while running in the back from `rstan`. The number of iterations,
+thinning, and other parameters can be customised.
 
-# Current version of the package runs three different models of the FoI
+::: {.alert .alert-primary}
+NOTE: Running the *serofoi* models for the first time on your local
+computer make take a few minutes while the *rstan* code is compiled
+locally. Afterwards, no further local compilation is needed.
+:::
 
-# Constant Force-of-Infection with a binomial distribution
-model_0 <- run_model(model_data = data_test,
-                     model_name = "constant_foi_bi")
+#### Model 1. Constant Force-of-Infection (endemic model)
 
-# Time-varying Force-of-Infection with a prior normal-binomial distribution
+For the *endemic model* a small number of iterations is enough for
+achieving convergence, as it only fits one parameter (the constant FoI)
+from a binomial distribution.
+
+```{r model_1, include = TRUE, errors = FALSE, warning = FALSE, message = FALSE }
 model_1 <- run_model(model_data = data_test,
-                     model_name = "continuous_foi_normal_bi")
+                     model_name = "constant_foi_bi",
+                     n_iters = 500, 
+                     n_thin = 2)
 
-# Time-varying Force-of-Infection with a prior normal-log distribution
-model_2 <- run_model(model_data = data_test,
-                     model_name = "continuous_foi_normal_log")
+```
 
-# The following function can be used to visualize the sero-prevalence data with its corresponding binomial confidence interval before fitting to a model
-plot_seroprev(data_test, size_text = 15)
+#### Model 2. Time-varying Force-of-Infection (epidemic model)
+
+For the *epidemic model,* a larger number of iterations is required for
+achieving convergence, as it fits yearly FoI values from a binomial
+distribution. The number of iterations required may depend on the number
+of years, reflected by the difference between year of the serosurvey and
+the maximum age-class sampled.
 
-# For each model, you can use ploting functions such as
-plot_model(model_0, size_text = 6)
-plot_seroprev(model_0, size_text = 6)
-plot_foi(model_0, size_text = 6)
-plot_rhats(model_0, size_text = 6)
-extract_summary_model(model_0)
+```{r model_2, include = TRUE, errors = FALSE, warning = FALSE, message = FALSE }
+model_2 <- run_model(model_data = data_test,
+                     model_name = "continuous_foi_normal_bi",
+                     n_iters = 1500, 
+                     n_thin = 2)
+```
 
-# You can compare these three models based on convergence, elpd and p-values
-comp_table <- get_comparison_table(
-  model_objects_list = c(m0 = model_0,
-                         m1 = model_1,
-                         m2 = model_2))
+#### Model 3. Time-varying Force-of-Infection (fast epidemic model)
 
+For the *fast* *epidemic model,* a larger number of iterations is
+required for achieving convergence, compared to the previous models.
 
+```{r model_3, include = TRUE, errors = FALSE, warning = FALSE, message = FALSE }
+model_3 <- run_model(model_data = data_test,
+                     model_name = "continuous_foi_normal_log",
+                     n_iters = 1500, 
+                     n_thin = 2)
+```
 
+For each model, the plot_model function generate a vertical arrange of
+plots summarising the results of the model implementation plotting
+functions. Crucially, it shows the (expected) log-predictive density
+`elpd`, standard error `se`, and allows to check convergence based on
+`R-hat` convergence diagnostics.
 
+```{r plot_model, include = FALSE, echo = TRUE, errors = FALSE, warning = FALSE, message = FALSE, fig.height=12, fig.width=6}
 
+plot_model(model_2, size_text = 15)
 ```
 
-### Lifecycle
-This package is currently a *concept*, as defined by the [RECON software
-lifecycle](https://www.reconverse.org/lifecycle.html). This means that essential
-features and mechanisms are still being developed, and the package is not ready
-for use outside of the development team.
+Also, the `plot_seroprev_models_grid` allows a visual a comparison of
+the models based on the (expected) log-predictive density `elpd`,
+standard error `se`, and allows to check convergence based on `R-hat`
+convergence diagnostics.
+
+```{r, model_comp, include=FALSE, echo=TRUE, warning=FALSE, message=FALSE, fig.height=12, message=FALSE}
+p1 <- plot_model(model_1, size_text = 15)
+p2 <- plot_model(model_2, size_text = 15)
+p3 <- plot_model(model_3, size_text = 15)
 
+plot_seroprev_models_grid(p1, p2, p3, n_row = 1, n_col = 3)
+```
+
+For more detailed information and examples, please check the [online
+documentation](https://epiverse-trace.github.io/serofoi/articles) as
+package vignettes under Get Started.
 
 ### Contributions
-Contributions are welcome via [pull requests](https://github.com/TRACE-LAC/serofoi/pulls).
 
 Contributors to the project include:
 
-- [Zulma M. Cucunubá](https://github.com/zmcucunuba) (author)
-- [Nicolás Tórres] (author)
-- [Benjamin Lambert] (author)
-- [Pierre Nouvellet] (author)
+-   [Zulma M. Cucunubá](https://github.com/zmcucunuba) (author,
+    maintainer)
+
+-   [Nicolás Tórres](https://github.com/ntorresd) (author)
 
+-   [Benjamin Lambert](https://ben-lambert.com/about/) (author)
 
-- [Miguel Gamez](https://github.com/megamezl) (contributor)
-- [Geraldine Gómez](https://github.com/megamezl) (contributor)
-- [Jaime A. Pavlich-Mariscal](https://github.com/jpavlich) (contributor)
+-   [Pierre Nouvellet](https://github.com/pnouvellet) (author)
 
+-   [Miguel Gamez](https://github.com/megamezl) (contributor)
 
-### Code of Conduct
-Please note that the linelist project is released with a 
-[Contributor Code of Conduct](https://contributor-covenant.org/version/2/0/CODE_OF_CONDUCT.html).
+-   [Geraldine Gómez](https://github.com/megamezl) (contributor)
+
+-   [Jaime A. Pavlich-Mariscal](https://github.com/jpavlich)
+    (contributor)
+
+## Package vignettes
+
+More details on how to use ***serofoi*** can be found in the [online
+documentation as package
+vignettes](https://epiverse-trace.github.io/serofoi/), under "Get
+Started".
+
+## Help
+
+To report a bug please open an
+[issue](https://github.com/epiverse-trace/serofoi/issues/new/choose).
+
+## Contribute
+
+Contributions to ***serofoi*** are welcomed. Please follow the [package
+contributing
+guide](https://github.com/epiverse-trace/serofoi/blob/main/.github/CONTRIBUTING.md).
+
+## Code of conduct
+
+Please note that the ***serofoi*** project is released with a
+[Contributor Code of
+Conduct](https://github.com/epiverse-trace/.github/blob/main/CODE_OF_CONDUCT.md).
 By contributing to this project, you agree to abide by its terms.
+
+## References
diff --git a/README.md b/README.md
index 9f2f9269..5c9d4569 100644
--- a/README.md
+++ b/README.md
@@ -1,7 +1,5 @@
 
-## *serofoi*<img src="man/figures/serofoi-logo.png" align="right" width="120" />
-
-__serofoi version 0.1.1__ is an R package to estimate the *Force-of-Infection* of a given pathogen from population based sero-prevalence studies on a Bayesian framework.
+## *serofoi*: force-of-infection from population based serosurveys with age-disagregated data <img src="man/figures/serofoi-logo.png" align="right" width="130"/>
 
 <!-- badges: start -->
 
@@ -12,8 +10,22 @@ MIT](https://img.shields.io/badge/License-MIT-yellow.svg)](https://opensource.or
 coverage](https://codecov.io/gh/TRACE-LAC/serofoi/branch/bugfixes-jaime/graph/badge.svg)](https://github.com/TRACE-LAC/serofoi?branch=bugfixes-jaime)
 
 [![lifecycle-concept](https://raw.githubusercontent.com/reconverse/reconverse.github.io/master/images/badge-concept.svg)](https://www.reconverse.org/lifecycle.html#concept)
+
 <!-- badges: end -->
 
+***serofoi*** is an R package to estimates the *Force-of-Infection* of a
+given pathogen from age-disaggregated population based sero-prevalence
+studies, using a Bayesian framework.
+
+***serofoi*** implements methods outlined in ([Cucunubá et al.
+2017](#ref-cucunubá2017)) and ([Carrera et al. 2020](#ref-carrera2020))
+
+***serofoi*** relies on
+[`rstan`](https://mc-stan.org/users/interfaces/rstan)
+
+***serofoi*** is part of the [Epiverse
+Initiative](https://data.org/initiatives/epiverse/).
+
 ## Installation
 
 You can install the **development version** of `serofoi` from
@@ -21,113 +33,226 @@ You can install the **development version** of `serofoi` from
 
 ``` r
 # install.packages("remotes")
-remotes::install_github("TRACE-LAC/serofoi")
+# remotes::install_github("TRACE-LAC/serofoi")
 library(serofoi)
 ```
 
 ## Quick start
 
-These examples illustrate some of the current functionalities:
-
-The function `prepare_data()` helps the user to prepare the dataset for the use of `serofoi` package
+The package provides an example dataset of the observed serosurvey data,
+`mydata`. This example is the basic entry for the package.
 
 ``` r
-data_test <- prepare_data(mydata)
-head(data_test)
+head(mydata)
+#>       survey total counts age_min age_max year_init year_end tsur country
+#> 1 COL-035-18     2      0       1       1      2007     2007 2007     COL
+#> 2 COL-035-18     1      0       2       2      2007     2007 2007     COL
+#> 3 COL-035-18    13      2       4       4      2007     2007 2007     COL
+#> 4 COL-035-18    25      5       5       5      2007     2007 2007     COL
+#> 5 COL-035-18    17      0       6       6      2007     2007 2007     COL
+#> 6 COL-035-18    20      4       7       7      2007     2007 2007     COL
+#>      test antibody
+#> 1 ELISA .      IgG
+#> 2 ELISA .      IgG
+#> 3 ELISA .      IgG
+#> 4 ELISA .      IgG
+#> 5 ELISA .      IgG
+#> 6 ELISA .      IgG
 ```
 
-Current version of the package runs three different models of the FoI. The function `run_model()` allows to choose betwwen three different models.
-
+The function `prepare_data` will prepare the entry data for entering the
+modelling functions. The seroprevalence *prepared data* can be
+visualised with the `plot_seroprev` function. This function also plots
+the binomial confidence interval of the observed data.
 
-- Constant Force-of-Infection with a binomial distribution
 ``` r
-model_0 <- run_model(model_data = data_test,
-                     model_name = "constant_foi_bi")
-``` 
-
-- Time-varying Force-of-Infection with a prior normal-binomial distribution
-``` r
-model_1 <- run_model(model_data = data_test,
-                     model_name = "continuous_foi_normal_bi")
-```
 
-- Time-varying Force-of-Infection with a prior normal-log distribution
-``` r
-model_2 <- run_model(model_data = data_test,
-                     model_name = "continuous_foi_normal_log")
-``` 
+data_test <- prepare_data(mydata)
 
-Function `plot_seroprev()` can be used to visualize the sero-prevalence data with its corresponding binomial confidence interval before fitting to a model:
-```
 plot_seroprev(data_test, size_text = 15)
 ```
-![](man/figures/plot_seroprev_example.png)
 
-For each model, there are three plotting functions:
-```
-plot_seroprev_fitted(model_0)
-``` 
-![](man/figures/plot_seroprev_fitted_example.png)
+<img src="man/figures/README-data_test-1.png" width="100%" />
 
+#### Current version of the package runs ***three*** different FoI models
 
-```
-plot_foi(model_0)
-``` 
-![](man/figures/plot_foi_example.png)
+The `run_model` function allows specifying the Bayesian model from *R*,
+while running in the back from `rstan`. The number of iterations,
+thinning, and other parameters can be customised.
 
+<div class="alert alert-primary">
 
-```
-plot_rhats(model_0)
-``` 
-![](man/figures/plot_rhats_example.png)
+NOTE: Running the *serofoi* models for the first time on your local
+computer make take a few minutes while the *rstan* code is compiled
+locally. Afterwards, no further compilation is needed.
+
+</div>
 
+#### Model 1. Constant Force-of-Infection (endemic model)
+
+For the *endemic model* a small number of iterations is enough for
+achieving convergence, as it only fits one parameter (the constant FoI)
+from a binomial distribution.
 
-The three plots can be obtained at once with the function `plot_model()`
 ``` r
-plot_model(model_0)
+model_1 <- run_model(model_data = data_test,
+                     model_name = "constant_foi_bi",
+                     n_iters = 500, 
+                     n_thin = 2)
+#> [1] "serofoi model constant_foi_bi finished running ------"
+#>            [,1]             
+#> model_name "constant_foi_bi"
+#> dataset    "COL-035-18"     
+#> country    "COL"            
+#> year       "2007"           
+#> test       "ELISA ."        
+#> antibody   "IgG"            
+#> n_sample   "212"            
+#> n_agec     "27"             
+#> n_iter     "500"            
+#> elpd       "-30.85"         
+#> se         "4.9"            
+#> converged  "Yes"
 ```
-![](man/figures/plot_model_example.png)
 
+#### Model 2. Time-varying Force-of-Infection (epidemic model)
 
-Finally, the package provides a funcion `get_comparison_table()`  allows to compare these three models based on convergence, elpd and p-values
+For the *epidemic model,* a larger number of iterations is required for
+achieving convergence, as it fits yearly FoI values from a binomial
+distribution. The number of iterations required may depend on the number
+of years, reflected by the difference between year of the serosurvey and
+the maximum age-class sampled.
 
 ``` r
-comp_table <- get_comparison_table(
-  model_objects_list = c(m0 = model_0,
-                         m1 = model_1,
-                         m2 = model_2))
-                         
-``` 
+model_2 <- run_model(model_data = data_test,
+                     model_name = "continuous_foi_normal_bi",
+                     n_iters = 1500, 
+                     n_thin = 2)
+#> [1] "serofoi model continuous_foi_normal_bi finished running ------"
+#>            [,1]                      
+#> model_name "continuous_foi_normal_bi"
+#> dataset    "COL-035-18"              
+#> country    "COL"                     
+#> year       "2007"                    
+#> test       "ELISA ."                 
+#> antibody   "IgG"                     
+#> n_sample   "212"                     
+#> n_agec     "27"                      
+#> n_iter     "1500"                    
+#> elpd       "-29.82"                  
+#> se         "5.28"                    
+#> converged  "Yes"
+```
+
+#### Model 3. Time-varying Force-of-Infection (fast epidemic model)
 
+For the *fast* *epidemic model,* a larger number of iterations is
+required for achieving convergence, compared to the previous models.
 
+``` r
+model_3 <- run_model(model_data = data_test,
+                     model_name = "continuous_foi_normal_log",
+                     n_iters = 1500, 
+                     n_thin = 2)
+#> [1] "serofoi model continuous_foi_normal_log finished running ------"
+#>            [,1]                       
+#> model_name "continuous_foi_normal_log"
+#> dataset    "COL-035-18"               
+#> country    "COL"                      
+#> year       "2007"                     
+#> test       "ELISA ."                  
+#> antibody   "IgG"                      
+#> n_sample   "212"                      
+#> n_agec     "27"                       
+#> n_iter     "1500"                     
+#> elpd       "-29.73"                   
+#> se         "5.53"                     
+#> converged  "Yes"
+```
 
+For each model, the plot_model function generate a vertical arrange of
+plots summarising the results of the model implementation plotting
+functions. Crucially, it shows the (expected) log-predictive density
+`elpd`, standard error `se`, and allows to check convergence based on
+`R-hat` convergence diagnostics.
 
-### Lifecycle
+Also, the `plot_models_list` allows a visual a comparison of the models
+based on the (expected) log-predictive density `elpd`, standard error
+`se`, and allows to check convergence based on `R-hat` convergence
+diagnostics.
 
-This package is currently a *concept*, as defined by the [RECON software
-lifecycle](https://www.reconverse.org/lifecycle.html). This means that
-essential features and mechanisms are still being developed, and the
-package is not ready for use outside of the development team.
+For more detailed information and examples, please check the [online
+documentation](https://epiverse-trace.github.io/serofoi/articles) as
+package vignettes under Get Started.
 
 ### Contributions
 
 Contributors to the project include:
 
-- [Zulma M. Cucunubá](https://github.com/zmcucunuba) (author)
+- [Zulma M. Cucunubá](https://github.com/zmcucunuba) (author,
+  maintainer)
+
 - [Nicolás Tórres](https://github.com/ntorresd) (author)
-- Benjamin Lambert (author)
-- Pierre Nouvellet (author)
+
+- [Benjamin Lambert](https://ben-lambert.com/about/) (author)
+
+- [Pierre Nouvellet](https://github.com/pnouvellet) (author)
+
 - [Miguel Gamez](https://github.com/megamezl) (contributor)
-- [Jaime Pavlich-Mariscal](https://github.com/jpavlich) (contributor)
-- [Geraldine Gómez](https://github.com/GeraldineGomez) (contributor)
 
-Contributions are welcome via [pull
-requests](https://github.com/TRACE-LAC/serofoi/pulls).
+- [Geraldine Gómez](https://github.com/megamezl) (contributor)
+
+- [Jaime A. Pavlich-Mariscal](https://github.com/jpavlich) (contributor)
+
+## Package vignettes
+
+More details on how to use ***serofoi*** can be found in the [online
+documentation as package
+vignettes](https://epiverse-trace.github.io/serofoi/), under “Get
+Started”.
+
+## Help
 
+To report a bug please open an
+[issue](https://github.com/epiverse-trace/serofoi/issues/new/choose).
 
-### Code of Conduct
+## Contribute
 
-Please note that the linelist project is released with a [Contributor
-Code of
-Conduct](https://contributor-covenant.org/version/2/0/CODE_OF_CONDUCT.html).
+Contributions to ***serofoi*** are welcomed. Please follow the [package
+contributing
+guide](https://github.com/epiverse-trace/serofoi/blob/main/.github/CONTRIBUTING.md).
+
+## Code of conduct
+
+Please note that the ***serofoi*** project is released with a
+[Contributor Code of
+Conduct](https://github.com/epiverse-trace/.github/blob/main/CODE_OF_CONDUCT.md).
 By contributing to this project, you agree to abide by its terms.
+
+## References
+
+<div id="refs" class="references csl-bib-body hanging-indent">
+
+<div id="ref-carrera2020" class="csl-entry">
+
+Carrera, Jean-Paul, Zulma M. Cucunubá, Karen Neira, Ben Lambert, Yaneth
+Pittí, Jesus Liscano, Jorge L. Garzón, et al. 2020. “Endemic and
+Epidemic Human Alphavirus Infections in Eastern Panama: An Analysis of
+Population-Based Cross-Sectional Surveys.” *The American Journal of
+Tropical Medicine and Hygiene* 103 (6): 2429–37.
+<https://doi.org/10.4269/ajtmh.20-0408>.
+
+</div>
+
+<div id="ref-cucunubá2017" class="csl-entry">
+
+Cucunubá, Zulma M, Pierre Nouvellet, Lesong Conteh, Mauricio Javier
+Vera, Victor Manuel Angulo, Juan Carlos Dib, Gabriel Jaime Parra -Henao,
+and María Gloria Basáñez. 2017. “Modelling Historical Changes in the
+Force-of-Infection of Chagas Disease to Inform Control and Elimination
+Programmes: Application in Colombia.” *BMJ Global Health* 2 (3):
+e000345. <https://doi.org/10.1136/bmjgh-2017-000345>.
+
+</div>
+
+</div>
diff --git a/_pkgdown.yml b/_pkgdown.yml
new file mode 100644
index 00000000..89edf49d
--- /dev/null
+++ b/_pkgdown.yml
@@ -0,0 +1,4 @@
+url: https://trace-lac.github.io/serofoi/
+template:
+  bootstrap: 5
+
diff --git a/data/chagas2012.RDS b/data/chagas2012.RDS
new file mode 100644
index 00000000..885b2727
Binary files /dev/null and b/data/chagas2012.RDS differ
diff --git a/data/chagas2012.RData b/data/chagas2012.RData
new file mode 100644
index 00000000..352fe7b7
Binary files /dev/null and b/data/chagas2012.RData differ
diff --git a/data/chik2015.RDS b/data/chik2015.RDS
new file mode 100644
index 00000000..390300e3
Binary files /dev/null and b/data/chik2015.RDS differ
diff --git a/data/chik2015.RData b/data/chik2015.RData
new file mode 100644
index 00000000..af87681d
Binary files /dev/null and b/data/chik2015.RData differ
diff --git a/data/data.RDS b/data/data.RDS
deleted file mode 100644
index 78f058b9..00000000
Binary files a/data/data.RDS and /dev/null differ
diff --git a/data/data.RData b/data/data.RData
deleted file mode 100644
index d6cbc1f0..00000000
Binary files a/data/data.RData and /dev/null differ
diff --git a/data/fake_data_1569526695.RDS b/data/fake_data_1569526695.RDS
new file mode 100644
index 00000000..9d0ad0b3
Binary files /dev/null and b/data/fake_data_1569526695.RDS differ
diff --git a/data/mydata.RData b/data/mydata.RData
new file mode 100644
index 00000000..2e8a51fe
Binary files /dev/null and b/data/mydata.RData differ
diff --git a/data/veev2012.RDS b/data/veev2012.RDS
new file mode 100644
index 00000000..c021ca47
Binary files /dev/null and b/data/veev2012.RDS differ
diff --git a/data/veev2012.RData b/data/veev2012.RData
new file mode 100644
index 00000000..8d33994f
Binary files /dev/null and b/data/veev2012.RData differ
diff --git a/man/chagas2012.Rd b/man/chagas2012.Rd
new file mode 100644
index 00000000..39d0832a
--- /dev/null
+++ b/man/chagas2012.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/chagas2012.R
+\docType{data}
+\name{chagas2012}
+\alias{chagas2012}
+\title{Seroprevalence data on serofoi}
+\format{
+An object of class \code{"cross"}; see \code{\link[qtl]{read.cross}}.
+}
+\usage{
+chagas2012
+}
+\description{
+Data from a serological surveys
+}
+\examples{
+chagas2012
+}
+\keyword{datasets}
diff --git a/man/chik2015.Rd b/man/chik2015.Rd
new file mode 100644
index 00000000..044418ba
--- /dev/null
+++ b/man/chik2015.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/chik2015.R
+\docType{data}
+\name{chik2015}
+\alias{chik2015}
+\title{Seroprevalence data on serofoi}
+\format{
+An object of class \code{"cross"}; see \code{\link[qtl]{read.cross}}.
+}
+\usage{
+chik2015
+}
+\description{
+Data from a serological surveys
+}
+\examples{
+chik2015
+}
+\keyword{datasets}
diff --git a/man/figures/README-data_test-1.png b/man/figures/README-data_test-1.png
new file mode 100644
index 00000000..b01c9a91
Binary files /dev/null and b/man/figures/README-data_test-1.png differ
diff --git a/man/figures/README-ex-1.png b/man/figures/README-ex-1.png
new file mode 100644
index 00000000..5ab8b89a
Binary files /dev/null and b/man/figures/README-ex-1.png differ
diff --git a/man/figures/README-model_comp-1.png b/man/figures/README-model_comp-1.png
new file mode 100644
index 00000000..cf9fdebe
Binary files /dev/null and b/man/figures/README-model_comp-1.png differ
diff --git a/man/figures/README-plot_model-1.png b/man/figures/README-plot_model-1.png
new file mode 100644
index 00000000..268d4f26
Binary files /dev/null and b/man/figures/README-plot_model-1.png differ
diff --git a/man/mydata.Rd b/man/mydata.Rd
new file mode 100644
index 00000000..b50dfdd0
--- /dev/null
+++ b/man/mydata.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/mydata.R
+\docType{data}
+\name{mydata}
+\alias{mydata}
+\title{Seroprevalence data on serofoi}
+\format{
+An object of class \code{"cross"}; see \code{\link[qtl]{read.cross}}.
+}
+\usage{
+mydata
+}
+\description{
+Data from a serological surveys
+}
+\examples{
+mydata
+}
+\keyword{datasets}
diff --git a/man/run_model.Rd b/man/run_model.Rd
index 93b01e2f..0e26bede 100644
--- a/man/run_model.Rd
+++ b/man/run_model.Rd
@@ -2,7 +2,7 @@
 % Please edit documentation in R/modelling.R
 \name{run_model}
 \alias{run_model}
-\title{Run model}
+\title{Run the specified stan model for the force-of-infection}
 \usage{
 run_model(
   model_data,
@@ -61,7 +61,7 @@ For further details refer to the \code{control} parameter in \link[rstan]{sampli
 \code{model_object}. An object containing relevant information about the implementation of the model. For further details refer to \link{fit_model}.
 }
 \description{
-Runs the specified stan model for the force-of-infection
+Run the specified stan model for the force-of-infection
 }
 \examples{
 \dontrun{
diff --git a/man/save_or_load_model.Rd b/man/save_or_load_model.Rd
index d64c0419..777ca6e7 100644
--- a/man/save_or_load_model.Rd
+++ b/man/save_or_load_model.Rd
@@ -19,7 +19,7 @@ save_or_load_model(model_name = "constant_foi_bi")
 }
 \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 
+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{
diff --git a/man/serofoi-package.Rd b/man/serofoi-package.Rd
index 07ea8f2a..6c1380ea 100644
--- a/man/serofoi-package.Rd
+++ b/man/serofoi-package.Rd
@@ -6,14 +6,21 @@
 \alias{serofoi-package}
 \title{serofoi: Estimates the Force-of-Infection of a given pathogen from population based sero-prevalence studies}
 \description{
-R package to estimate time-varying Force-of-Infection of a given pathogen from population based sero-prevalence studies on a bayesian framework.
+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}: N. T. Domínguez \email{ex-ntorres@javeriana.edu.co}
+\strong{Maintainer}: Zulma M. Cucunubá \email{zulma.cucunuba@javeriana.edu.co} (\href{https://orcid.org/0000-0002-8165-3198}{ORCID})
 
 Authors:
 \itemize{
-  \item Zulma M. Cucunubá
+  \item Nicolás Torres
   \item Ben Lambert
   \item Pierre Nouvellet
 }
diff --git a/man/veev2012.Rd b/man/veev2012.Rd
new file mode 100644
index 00000000..c095cc7f
--- /dev/null
+++ b/man/veev2012.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/veev2012.R
+\docType{data}
+\name{veev2012}
+\alias{veev2012}
+\title{Seroprevalence data on serofoi}
+\format{
+An object of class \code{"cross"}; see \code{\link[qtl]{read.cross}}.
+}
+\usage{
+veev2012
+}
+\description{
+Data from a serological surveys
+}
+\examples{
+veev2012
+}
+\keyword{datasets}
diff --git a/references.bib b/references.bib
new file mode 100644
index 00000000..fc3a76e1
--- /dev/null
+++ b/references.bib
@@ -0,0 +1,30 @@
+
+@article{cucunubá2017,
+	title = {Modelling historical changes in the force-of-infection of Chagas disease to inform control and elimination programmes: application in Colombia},
+	author = {{Cucunubá}, Zulma M and Nouvellet, Pierre and Conteh, Lesong and Vera, Mauricio Javier and Angulo, Victor Manuel and Dib, Juan Carlos and Parra -Henao, Gabriel Jaime and {Basáñez}, {María Gloria}},
+	year = {2017},
+	month = {09},
+	date = {2017-09},
+	journal = {BMJ Global Health},
+	pages = {e000345},
+	volume = {2},
+	number = {3},
+	doi = {10.1136/bmjgh-2017-000345},
+	url = {http://dx.doi.org/10.1136/bmjgh-2017-000345},
+	langid = {en}
+}
+
+@article{carrera2020,
+	title = {Endemic and Epidemic Human Alphavirus Infections in Eastern Panama: An Analysis of Population-Based Cross-Sectional Surveys},
+	author = {Carrera, Jean-Paul and {Cucunubá}, Zulma M. and Neira, Karen and Lambert, Ben and {Pittí}, Yaneth and Liscano, Jesus and {Garzón}, Jorge L. and Beltran, Davis and Collado-Mariscal, Luisa and Saenz, Lisseth and Sosa, {Néstor} and Rodriguez-Guzman, Luis D. and {González}, Publio and Lezcano, {Andrés G.} and {Pereyra-Elías}, {Reneé} and Valderrama, Anayansi and Weaver, Scott C. and Vittor, Amy Y. and {Armién}, Blas and Pascale, Juan-Miguel and Donnelly, Christl A.},
+	year = {2020},
+	month = {12},
+	date = {2020-12-02},
+	journal = {The American Journal of Tropical Medicine and Hygiene},
+	pages = {2429--2437},
+	volume = {103},
+	number = {6},
+	doi = {10.4269/ajtmh.20-0408},
+	url = {http://dx.doi.org/10.4269/ajtmh.20-0408},
+	langid = {en}
+}
diff --git a/test/testing_vignette.R b/test/testing_vignette.R
new file mode 100644
index 00000000..17b62975
--- /dev/null
+++ b/test/testing_vignette.R
@@ -0,0 +1,15 @@
+library(remotes)
+remotes::install_github("TRACE-LAC/serofoi",
+                        ref = "dev",
+                        force = TRUE)
+
+
+library(serofoi)
+library(usethis)
+library(roxygen2)
+
+devtools::document()
+usethis::use_pkgdown_github_pages()
+pkgdown::build_site()
+use_github_pages(branch = "dev-webrd", path = "/", cname = NA)
+
diff --git a/tests/test_sim_data.R b/tests/test_sim_data.R
new file mode 100644
index 00000000..ff83827c
--- /dev/null
+++ b/tests/test_sim_data.R
@@ -0,0 +1,134 @@
+rm(list=ls())
+
+
+
+# Generate exposure matrix
+
+generate_fake_data <- function(n, foi, name_example = 'fake', grouping = FALSE) {
+  # dat <- readRDS('data/chik_serop_IND 1965(S095).RDS')
+  dat <- data.frame(birth_year=c(2000:2049)) %>%
+    mutate(tSur=rep(2050, length(birth_year)),
+           age_mean_f=2050-birth_year) %>%
+    mutate(tsur=tSur)
+  yexpo <- make_yexpo(dat)
+  yexpo <- yexpo[-length(yexpo)]
+  RealYexpo <- (min(dat$birth_year):dat$tsur[1])[-1]
+  ExposureMatrix <- get_exposure_matrix(dat, yexpo)
+  Nobs <- nrow(dat)
+
+  fGenerateProbs <- function(ExposureMatrix, foi){
+    p <- map_dbl(1:nrow(ExposureMatrix), ~1-exp(-dot(ExposureMatrix[., ], foi)))
+    return(p)
+  }
+
+  fGenerateFakeData <- function(foi, sample_size, ExposureMatrix){
+    p <- fGenerateProbs(ExposureMatrix, foi)
+    counts <- map_int(p, ~rbinom(1, sample_size, .))
+    return(counts)
+  }
+
+
+  dat$counts <- fGenerateFakeData(foi, n, ExposureMatrix)
+  dat$total <- n
+  dat$survey <- name_example
+  dat$country <- 'None'
+
+  dat <- dat %>% mutate(age_group = 'NA', age = age_mean_f)
+  dat$age_group[dat$age >0 & dat$age <5] <-   '01-04'
+  dat$age_group[dat$age >4 & dat$age <10] <-  '05-09'
+  dat$age_group[dat$age >9 & dat$age <15] <-  '10-14'
+  dat$age_group[dat$age >14 & dat$age <20] <- '15-19'
+  dat$age_group[dat$age >19 & dat$age <25] <- '20-24'
+  dat$age_group[dat$age >24 & dat$age <30] <- '25-29'
+  dat$age_group[dat$age >29 & dat$age <35] <- '30-34'
+  dat$age_group[dat$age >34 & dat$age <40] <- '35-39'
+  dat$age_group[dat$age >39 & dat$age <45] <- '40-44'
+  dat$age_group[dat$age >44 & dat$age <51] <- '45-50'
+
+  dat <- dat%>% arrange(age)
+
+  if (grouping == TRUE) {
+
+    datg <- dat %>% group_by(age_group) %>%
+      dplyr::summarise(total = sum(total), counts = sum(counts)) %>%
+      mutate(tsur = dat$tsur[1], tSur = dat$tSur[1]) %>%
+      mutate(age_min = as.numeric(substr(age_group, 1, 2)),
+             age_max = as.numeric(substr(age_group, 4, 5))) %>%
+      mutate(age_mean_f = floor((age_min + age_max)/2)) %>%
+      mutate (age = age_mean_f, age_mean = age_mean_f,
+              survey = name_example,
+              country = 'None') %>%
+      mutate(birth_year = tsur - age_mean_f)
+
+    dat <- datg[, names(dat)]
+
+  }
+
+  dat <- dat%>% arrange(age)
+  conf_int <- data.frame(Hmisc::binconf(x=dat$counts,n = dat$total, alpha = 0.05, method ='exact'))
+  dat <- cbind(dat, conf_int)
+  dat <- dat %>% rename(prev_obs_lower = Lower,
+                        prev_obs = PointEst,
+                        prev_obs_upper = Upper)
+
+
+  return(dat)
+}
+
+
+no_transm <- 0.0000000001
+small_outbreak <- 0.5
+big_outbreak   <- 1.5
+
+foiA <- rep(0.02, 50)
+foiB <- c(rep(0.2, 25), rep(0.1, 10), rep(no_transm, 15)) # interruption
+foiC <- seq(no_transm, 0.05, length.out = 50) # interruption
+foiD <- c(rep(no_transm, 32), rep(big_outbreak, 3), rep(no_transm, 15)) # 1 epidemics
+foiE <- c(rep(no_transm, 37), rep(small_outbreak, 3), rep(no_transm, 10)) # 1 epidemics
+foiF <- c(rep(no_transm, 11), rep(big_outbreak, 2), rep(no_transm, 23), rep(small_outbreak, 3), rep(no_transm, 11)) # 1 epidemics
+
+genetate_4_datasets <- function (foi, name_foi)  {
+  datfoiA1 <- generate_fake_data(n = 5, foi,  paste0(name_foi, '_n05_sing')) %>% mutate(simulation = name_foi)
+  datfoiA2 <- generate_fake_data(n = 5, foi,  paste0(name_foi, '_n05_group'), grouping = TRUE) %>% mutate(simulation = name_foi)
+  datfoiA3 <- generate_fake_data(n = 10, foi, paste0(name_foi, '_n10_sing')) %>% mutate(simulation = name_foi)
+  datfoiA4 <- generate_fake_data(n = 10, foi, paste0(name_foi, '_n10_group'), grouping = TRUE)%>% mutate(simulation = name_foi)
+  datasets <- list(foi1 = datfoiA1,
+                   foi2 = datfoiA2,
+                   foi3 = datfoiA3,
+                   foi4 = datfoiA4)
+  return(datasets)
+}
+
+foiA_ds <- genetate_4_datasets(foiA, 'foiA')
+foiB_ds <- genetate_4_datasets(foiB, 'foiB')
+foiC_ds <- genetate_4_datasets(foiC, 'foiC')
+foiD_ds <- genetate_4_datasets(foiD, 'foiD')
+foiE_ds <- genetate_4_datasets(foiE, 'foiE')
+foiF_ds <- genetate_4_datasets(foiF, 'foiF')
+
+dat0 <- rbind (foiA_ds$foi1, foiA_ds$foi2, foiA_ds$foi3, foiA_ds$foi4,
+               foiB_ds$foi1, foiB_ds$foi2, foiB_ds$foi3, foiB_ds$foi4,
+               foiC_ds$foi1, foiC_ds$foi2, foiC_ds$foi3, foiC_ds$foi4,
+               foiD_ds$foi1, foiD_ds$foi2, foiD_ds$foi3, foiD_ds$foi4,
+               foiE_ds$foi1, foiE_ds$foi2, foiE_ds$foi3, foiE_ds$foi4,
+               foiF_ds$foi1, foiF_ds$foi2, foiF_ds$foi3, foiF_ds$foi4)
+
+dat0 <- dat0 %>% mutate(Antibody = 'IgG', Test = 'Fake')
+# Here we assume the test used is IgG
+dat0 %>%
+  ggplot(aes(x=age, y=prev_obs)) +
+  geom_point() +
+  facet_wrap(~survey, nrow = 2)
+
+
+fake_data <- list  (dat0 = dat0,
+                    foiA = foiA,
+                    foiB = foiB,
+                    foiC = foiC,
+                    foiD = foiD,
+                    foiE = foiE,
+                    foiF = foiF)
+
+
+(name_file <- paste0('data/fake_data_', round(as.numeric(Sys.time()), 0), '.RDS'))
+saveRDS(fake_data, name_file)
diff --git a/tests/testepidemics.R b/tests/testepidemics.R
new file mode 100644
index 00000000..139597f9
--- /dev/null
+++ b/tests/testepidemics.R
@@ -0,0 +1,2 @@
+
+
diff --git a/vignettes/references.bib b/vignettes/references.bib
new file mode 100644
index 00000000..185bcf95
--- /dev/null
+++ b/vignettes/references.bib
@@ -0,0 +1,60 @@
+
+@article{Cucunubá2017,
+	title = {Modelling historical changes in the force-of-infection of Chagas disease to inform control and elimination programmes: application in Colombia},
+	author = {{Cucunubá}, Zulma M and Nouvellet, Pierre and Conteh, Lesong and Vera, Mauricio Javier and Angulo, Victor Manuel and Dib, Juan Carlos and Parra -Henao, Gabriel Jaime and {Basáñez}, {María Gloria}},
+	year = {2017},
+	month = {09},
+	date = {2017-09},
+	journal = {BMJ Global Health},
+	pages = {e000345},
+	volume = {2},
+	number = {3},
+	doi = {10.1136/bmjgh-2017-000345},
+	url = {http://dx.doi.org/10.1136/bmjgh-2017-000345},
+	langid = {en}
+}
+
+@article{Carrera2020,
+	title = {Endemic and Epidemic Human Alphavirus Infections in Eastern Panama: An Analysis of Population-Based Cross-Sectional Surveys},
+	author = {Carrera, Jean-Paul and {Cucunubá}, Zulma M. and Neira, Karen and Lambert, Ben and {Pittí}, Yaneth and Liscano, Jesus and {Garzón}, Jorge L. and Beltran, Davis and Collado-Mariscal, Luisa and Saenz, Lisseth and Sosa, {Néstor} and Rodriguez-Guzman, Luis D. and {González}, Publio and Lezcano, {Andrés G.} and {Pereyra-Elías}, {Reneé} and Valderrama, Anayansi and Weaver, Scott C. and Vittor, Amy Y. and {Armién}, Blas and Pascale, Juan-Miguel and Donnelly, Christl A.},
+	year = {2020},
+	month = {12},
+	date = {2020-12-02},
+	journal = {The American Journal of Tropical Medicine and Hygiene},
+	pages = {2429--2437},
+	volume = {103},
+	number = {6},
+	doi = {10.4269/ajtmh.20-0408},
+	url = {http://dx.doi.org/10.4269/ajtmh.20-0408},
+	langid = {en}
+}
+
+@article{cunha2017,
+	title = {Seroprevalence of Chikungunya Virus in a Rural Community in Brazil},
+	author = {Cunha, Rivaldo V. and Trinta, Karen S. and Montalbano, Camila A. and Sucupira, Michel V. F. and de Lima, Maricelia M. and Marques, Erenilde and Romanholi, Izilyanne H. and Croda, Julio},
+	editor = {Barrera, Roberto},
+	year = {2017},
+	month = {01},
+	date = {2017-01-20},
+	journal = {PLOS Neglected Tropical Diseases},
+	pages = {e0005319},
+	volume = {11},
+	number = {1},
+	doi = {10.1371/journal.pntd.0005319},
+	url = {http://dx.doi.org/10.1371/journal.pntd.0005319},
+	langid = {en}
+}
+
+@article{dias2018,
+	title = {Seroprevalence of Chikungunya Virus after Its Emergence in Brazil},
+	author = {Dias, Juarez P. and Costa, {Maria da Conceição N.} and Campos, Gubio Soares and {Paixão}, Enny S. and Natividade, Marcio S. and Barreto, Florisneide R. and Itaparica, Martha Suely C. and Goes, Cristina and Oliveira, Francisca L.S. and Santana, Eloisa B. and Silva, Neusa S.J. and Brito, Carlos A.A. and Rodrigues, Laura C. and Sardi, Silvia Inez and Saavedra, Ramon C. and Teixeira, {Maria Glória}},
+	year = {2018},
+	month = {04},
+	date = {2018-04},
+	journal = {Emerging Infectious Diseases},
+	pages = {617--624},
+	volume = {24},
+	number = {4},
+	doi = {10.3201/eid2404.171370},
+	url = {http://dx.doi.org/10.3201/eid2404.171370}
+}
diff --git a/vignettes/serofoi.Rmd b/vignettes/serofoi.Rmd
deleted file mode 100644
index 16334100..00000000
--- a/vignettes/serofoi.Rmd
+++ /dev/null
@@ -1,19 +0,0 @@
----
-title: "serofoi"
-output: rmarkdown::html_vignette
-vignette: >
-  %\VignetteIndexEntry{serofoi}
-  %\VignetteEngine{knitr::rmarkdown}
-  %\VignetteEncoding{UTF-8}
----
-
-```{r, include = FALSE}
-knitr::opts_chunk$set(
-  collapse = TRUE,
-  comment = "#>"
-)
-```
-
-```{r setup}
-library(serofoi)
-```
diff --git a/vignettes/use_cases.Rmd b/vignettes/use_cases.Rmd
new file mode 100644
index 00000000..4fc67b56
--- /dev/null
+++ b/vignettes/use_cases.Rmd
@@ -0,0 +1,235 @@
+---
+title: "serofoi"
+output: rmarkdown::html_vignette
+bibliography: references.bib
+link-citations: true
+vignette: >
+  %\VignetteIndexEntry{serofoi}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>"
+)
+```
+
+```{r cleaning, include = FALSE, echo = TRUE, message=FALSE}
+library(serofoi)
+library(tidyverse)
+rownames(mydata) <- NULL 
+data_test <- prepare_data(mydata)
+
+
+```
+
+# Epidemiological scenarios for the *serofoi* models  
+
+**serofoi** is an R package that allows estimations of the Force-of-Infection from a population based serosurvey data.
+
+## **Sero-survey data**
+
+Surveys had to meet all of the following [inclusion criteria]{.underline}:
+
+1.  Be population-based (not hospital based)
+
+2.  Specify individuals' age or age group
+
+3.  Indicate diagnostic test(s) used. The current version of serofoi only applies to *IgG* antibodies)
+
+4.  Identify the location and date (year) of sample collection
+
+## Model assumptions
+
+Current version of **serofoi** includes the following assumptions on the underlying biological process:
+
+-   The Force-of-infection (FoI) is estimated using a catalytic model
+
+-   There is no sero-reversion (from positive to negative). It means IgG antibodies are life-long duration with no waning immunity. This may not be the case for several pathogens. This feature is planned for future versions of *serofoi***.**
+
+-   There is no age-dependency. This may not be the case for several pathogens. This feature is planned for future versions of *serofoi***.**
+
+-   There is no impact from migration processes in the sampled population
+
+-   There are no differences in the mortality rate of infected versus susceptible individuals
+
+::: {.alert .alert-primary}
+NOTE: Running the *serofoi* models for the first time on your local computer make take a few minutes while the *rstan* code is compiled locally. Afterwards, no further local compilation is needed.
+:::
+
+# **Constant vs Time-varying** FoI models
+
+## Constant Force-of-Infection (endemic model)
+
+For *constant* *endemic model*, the rate of infection acquisition―rate of seroconversion―is constant over time, and infection (sero)prevalence will increase monotonically with age as cumulative exposure increases.
+
+Denoting $n(a,t)$ the number of seropositive individuals of age $a$ at time $t$ , $N$ the serosurvey sample size, and $P(a,t)$ the underlying seroprevalence at age $a$ at time $t$, we assumed that the number of seropositive subjects follows a Binomial distribution $n(a,t)$ \~ $B(N,P(a,t))$.
+
+For the force-of-infection (FoI) that is constant over time, denoted $\lambda$, we modelled seroprevalence for age $a$ in year $t$ (i.e. the time when the serosurvey occurred) as:
+
+$$ P(a,t) = 1−exp(−\lambda{a})$$
+
+The corresponding code for the constant-FoI model with binomial distribution is:
+
+```{r model_1, include = FALSE, echo=TRUE, errors = FALSE, warning = FALSE, message = FALSE }
+model_1 <- run_model(model_data = data_test,
+                     model_name = "constant_foi_bi",
+                     n_iters = 500, 
+                     n_thin = 2)
+
+```
+
+## Time-varying Force-of-Infection
+
+For the time-varying FoI model, the FoI values are estimated from a binomial distribution. The sero-prevalence as age a can be expressed as:
+
+$$
+𝑃(𝑎,𝑡)=1−exp(−∑_{i=t-a+1}^{t}𝜆{i})
+$$
+
+Therefore, a serosurvey completed at time $t$ and including ages $a$ from {$age_{min}, age_{max}$} is informative on exposure (and FoI) between {$t, age_{max}$} and $t$.
+
+Currently, two models within *serofoi* allow time-varying FoI, with corresponding code:
+
+```{r model_2, include = FALSE, echo = TRUE, errors = FALSE, warning = FALSE, message = FALSE }
+model_2 <- run_model(model_data = data_test,
+                     model_name = "continuous_foi_normal_bi",
+                     n_iters = 1500, 
+                     n_thin = 2)
+
+model_3 <- run_model(model_data = data_test,
+                     model_name = "continuous_foi_normal_log",
+                     n_iters = 1500, 
+                     n_thin = 2)
+```
+
+## Fitting process
+
+Models are implemented in **RStan** using the No-U-Turn sampler, a type of Hamiltonian Monte Carlo sampling. After convergence, the best-fitting models can be selected based on the expected log predictive density for an out-of-sample data point (elpd).
+
+## Case study 1. Chagas disease (endemic disease)
+
+Based on the data and analysis shown in [@Cucunubá2017], we use one of the datasets for measuring the sero-prevalence of IgG antibodies against *Trypanosoma cruzi* infection in rural area of Colombia in 2012. The dataset is part of the serofoi package as `chagas2012`.
+
+```{r data_chagas, include = TRUE, errors = FALSE, warning = FALSE, message = FALSE }
+chagas2012 <- readRDS("../data/chagas2012.RDS") # This can be removed once data is incorporated
+head(chagas2012, 5)
+
+```
+
+After preparing the data we can run two models (constant and time varying) available on serofoi.
+
+```{r three models, include = TRUE, errors = FALSE, warning = FALSE, message=FALSE }
+
+chagas2012p <- prepare_data(chagas2012)
+
+m1_cha <- run_model(model_data = chagas2012p,
+                     model_name = "constant_foi_bi",
+                     n_iters = 500, 
+                     n_thin = 2)
+
+m2_cha <- run_model(model_data = chagas2012p,
+                     model_name = "continuous_foi_normal_bi",
+                     n_iters = 1500, 
+                     n_thin = 2)
+
+```
+
+Now, we can plot the results of the two models to compare (Figure 1). As shown in @Cucunubá2017, interventions for Chagas control have been ongoing from the 1980s in Colombia having heterogeneous impact depending on the type of population. For this serosurvey, which is from tradittional indigenous rural area, the serofoi models are able to detect a modest still relevant slow decreasing trend consistent (model 2) as slightly better supported than the constant model. Notice that model 3 does not converge despite a high number of iterations.
+
+```{r, model_comp_cha, include=TRUE, echo=TRUE, fig.width = 8, fig.height=10, warning=FALSE, message=FALSE, message=FALSE}
+p1_cha <- plot_model(m1_cha, size_text = 12, max_lambda = 0.02)
+p2_cha <- plot_model(m2_cha, size_text = 12, max_lambda = 0.02)
+
+plot_seroprev_models_grid(p1_cha, p2_cha, n_row = 1, n_col = 2)
+```
+
+## Case study 2. Hidden Alphaviruses epidemics in Panama
+
+As shown in [@Carrera2020], hidden epidemic and endemic transmission of alphaviruses in Eastern Panama have been around for decades. From this paper we use a dataset measuring IgG antibodies againts Venezuelan Equine Encephalitis Virus (VEEV) in a rural village in Panamá in 2017. This dataset, `veev2017` is included in serofoi.
+
+```{r data_veev, include = TRUE, errors = FALSE, warning = FALSE, message = FALSE }
+
+veev2012 <- readRDS("../data/veev2012.RDS")
+head(veev2012, 5)
+
+```
+
+```{r models_veee, include = TRUE, errors = FALSE, warning = FALSE, message=FALSE }
+
+veev2012p <- prepare_data(veev2012)
+
+m1_veev <- run_model(model_data = veev2012p,
+                     model_name = "constant_foi_bi",
+                     n_iters = 500, 
+                     n_thin = 2)
+
+m2_veev <- run_model(model_data = veev2012p,
+                     model_name = "continuous_foi_normal_bi",
+                     n_iters = 500, 
+                     n_thin = 2)
+
+m3_veev <- run_model(model_data = veev2012p,
+                     model_name = "continuous_foi_normal_log",
+                     n_iters = 500, 
+                     n_thin = 2)
+```
+
+Now, we can plot the results of the three models to compare. On Figure 3, we can observe a large increase in the estimated FOI. As suggested in @Carrera2020, an important increase in the transmission of VEEV in this region is inferred. Although the three fitted models converge well, we see much larger support for a more sudden increase in recent years with highest FoI values at 0.25 (consistent with an epidemic).
+
+```{r, model_comp_veev, include=TRUE, echo=TRUE, fig.width = 10, fig.height=10, warning=FALSE, message=FALSE, message=FALSE}
+p1_veev <- plot_model(m1_veev, size_text = 10, max_lambda = 0.6)
+p2_veev <- plot_model(m2_veev, size_text = 10, max_lambda = 0.6)
+p3_veev <- plot_model(m3_veev, size_text = 10, max_lambda = 0.6)
+
+plot_seroprev_models_grid(p1_veev, p2_veev, p3_veev, n_row = 1, n_col = 3)
+
+```
+
+Figure 3. Results of fitted models, Force-of-Infection estimates and convergence.
+
+## Case study 3. A known Large Chikungunya epidemic in Brazil
+
+Chikungunya outbreaks ocurred rapidly after the introduction of the virus to Brazil in 2013-2014. We use the dataset from [@dias2018] that conducts a population-based study through household interviews and serologic surveys (measuring IgG antibodies against Chikungunya virus) in Bahia, Brazil during October-December 2015, right after a large CHIKV epidemic that occurred in that area.
+
+```{r data_chik, include = TRUE, errors = FALSE, warning = FALSE, message = FALSE }
+
+chik2015 <- readRDS("../data/chik2015.RDS")
+head(chik2015, 5)
+
+```
+
+```{r models_chik, include = TRUE, errors = FALSE, warning = FALSE, message=FALSE }
+
+chik2015p <- prepare_data(chik2015)
+
+mod1_chik <- run_model(model_data = chik2015p,
+                     model_name = "constant_foi_bi",
+                     n_iters = 1000, 
+                     n_thin = 2)
+
+mod2_chik <- run_model(model_data = chik2015p,
+                     model_name = "continuous_foi_normal_bi",
+                     n_iters = 1500, 
+                     n_thin = 2)
+
+mod3_chik <- run_model(model_data = chik2015p,
+                     model_name = "continuous_foi_normal_log",
+                     n_iters = 1500, 
+                     n_thin = 2)
+```
+
+In Figure 4, we can observe the comparison between the three serofoi models. Here serofoi shows strong statistical support for a sudden increase in the transmission of CHIKV close to the year of the serosurvey (2015). The exact year is not possible to estimate, mainly given the data used is largely aggregated by 20-years age groups. Despite that, the results are consistent with the empirical evidence shown by [@dias2018] with both interviews, and IgM testing.
+
+```{r, plot_chik, include=TRUE, echo=TRUE, fig.width = 10, fig.height=10, warning=FALSE, message=FALSE, message=FALSE}
+p1_chik <- plot_model(mod1_chik, size_text = 10, max_lambda = 0.08)
+p2_chik <- plot_model(mod2_chik, size_text = 10, max_lambda = 0.08)
+p3_chik <- plot_model(mod3_chik, size_text = 10, max_lambda = 0.08)
+
+plot_seroprev_models_grid(p1_chik, p2_chik, p3_chik, n_row = 1, n_col = 3)
+
+```
+
+## References