diff --git a/DESCRIPTION b/DESCRIPTION
index ec596a5..35c7d8a 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -30,7 +30,9 @@ Imports:
     etm,
     Rdpack
 Suggests:
+    kableExtra,
     knitr,
+    rmarkdown,
     testthat (>= 2.0),
     roxytypes,
     roxylint
diff --git a/vignettes/savvyr_example.Rmd b/vignettes/savvyr_example.Rmd
new file mode 100644
index 0000000..c8aa6d5
--- /dev/null
+++ b/vignettes/savvyr_example.Rmd
@@ -0,0 +1,114 @@
+---
+title: "Estimation of AE probabilities with savvyr"
+package: savvyr
+bibliography: "../inst/REFERENCES.bib"
+output:
+  rmarkdown::html_vignette:
+          toc: true
+vignette: |
+  %\VignetteEncoding{UTF-8}
+  %\VignetteIndexEntry{Estimation of AE probabilities with savvyr}
+  %\VignetteEngine{knitr::rmarkdown}
+editor_options:
+  chunk_output_type: console
+---
+
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>"
+)
+library(savvyr)
+library(kableExtra)
+```
+
+
+
+
+# Example using dummy data
+
+We generate the dataset $S1$ in @stegherr_meta_analytic_2021 using the parameter
+values for Arm A.  
+First we define sample size and range of censoring times. Then we set the hazard
+of the three event types (adverse event, death/hard competing event and soft
+competing event). After the dataset has been generated we set $\tau$ as the
+maximum event time.
+
+```{r, include=TRUE, echo=TRUE}
+n <- 200
+
+min_cens <- 0
+max_cens <- 1000
+
+haz_ae <- 0.00265
+haz_death <- 0.00151
+haz_soft <- 0.00227
+
+set.seed(2020)
+dat1 <- generate_data(n, cens = c(min_cens, max_cens), haz_ae, haz_death, haz_soft)
+
+tau <- max(dat1[, "time_to_event"])
+```
+
+The structure of the dataset looks as follows:
+
+```{r, include=TRUE, echo=TRUE}
+kable(head(dat1, 10), align = c("crcr"))
+```
+
+For this dataset we then compute all the estimators used in the comparisons 
+in @stegherr_survival_2021 and @stegherr_estimating_2021.  
+We start with the estimators that do not account for competing events (incidence
+proportion, incidence density, Inverse Kaplan Meier), then incidence proportion 
+accounting for competing events and Aalen-Johansen (both first with death only
+as hard competing event, then using all competing events):
+
+```{r, include=TRUE, echo=TRUE}
+ip <- inc_prop(dat1, tau)
+id <- prop_trans_inc_dens(dat1, tau)
+km <- one_minus_kaplan_meier(dat1, tau)
+
+idce_2 <- prop_trans_inc_dens_ce(dat1, ce = 2, tau)
+aj_2 <- aalen_johansen(dat1, ce = 2, tau)
+
+idce_3 <- prop_trans_inc_dens_ce(dat1, ce = 3, tau)
+aj_3 <- aalen_johansen(dat1, ce = 3, tau)
+```
+
+The AE risks look as follows:
+
+```{r, include=TRUE, echo=TRUE}
+tab <- rbind(ip, id, km, idce_2, aj_2[1:2], idce_3, aj_3[1:2])
+colnames(tab) <- c(
+  "estimated AE probability",
+  "variance of estimation"
+)
+rownames(tab) <- c(
+  "incidence proportion",
+  "probability transform incidence density ignoring competing event",
+  "1 - Kaplan-Meier", "probability transform incidence density (death only)",
+  "Aalen-Johansen (death only), AE risk", "probability transform incidence density (all CEs)",
+  "Aalen-Johansen (all CEs), AE risk"
+)
+
+kable(tab, digits = c(3, 5))
+```
+
+Finally, the estimated probabilities of competing events based on the
+Aalen-Johansen estimators:
+
+```{r, include=TRUE, echo=TRUE}
+tab <- rbind(aj_2[3:4], aj_3[3:4])
+colnames(tab) <- c(
+  "estimated probability",
+  "variance of estimation"
+)
+rownames(tab) <- c(
+  "Aalen-Johansen (death only), CE risk",
+  "Aalen-Johansen (all CEs), CE risk"
+)
+
+kable(tab, digits = c(3, 5))
+```
+
+# References