diff --git a/R/estimate_infections.R b/R/estimate_infections.R
index b7524ce23..34b29ec19 100644
--- a/R/estimate_infections.R
+++ b/R/estimate_infections.R
@@ -131,8 +131,10 @@
 #'
 #' 
 #' # Plot output
-#' report_plots(summarised_estimates = snapshot$summarised,
-#'                       reported = snapshot_cases)     
+#' plots <- report_plots(summarised_estimates = snapshot$summarised,
+#'                       reported = snapshot_cases)
+#'
+#' plots$summary    
 #' 
 #' ## Run model with stationary Rt assumption (likely to provide biased real-time estimates)
 #' stat <- estimate_infections(reported_cases, family = "negbin",
@@ -144,8 +146,10 @@
 #'
 #' 
 #' # Plot output
-#' report_plots(summarised_estimates = stat$summarised,
-#'              reported = reported_cases)
+#' plots <- report_plots(summarised_estimates = stat$summarised,
+#'                       reported = reported_cases)
+#'
+#' plots$summary
 #'        
 #' # Run model with fixed Rt assumption 
 #' fixed <- estimate_infections(reported_cases, family = "negbin",
@@ -157,8 +161,10 @@
 #'
 #' 
 #' # Plot output
-#' report_plots(summarised_estimates = fixed$summarised,
-#'              reported = reported_cases)
+#' plots <- report_plots(summarised_estimates = fixed$summarised,
+#'                       reported = reported_cases)
+#'                       
+#' plots$summary
 #' 
 #' # Run model with breakpoints                                                                      
 #' bkp <- estimate_infections(reported_cases, family = "negbin",
@@ -170,8 +176,10 @@
 #'
 #' 
 #' # Plot output
-#' report_plots(summarised_estimates = bkp$summarised,
-#'              reported = reported_cases)
+#' plots <- report_plots(summarised_estimates = bkp$summarised,
+#'                       reported = reported_cases)
+#'
+#' plots$summary
 #'              
 #' # Run model with breakpoints but with constrained non-linear change over time 
 #' # This formulation may increase the apparent effect of the breakpoint but needs to be tested using
@@ -187,9 +195,11 @@
 #'
 #' 
 #' # Plot output
-#' report_plots(summarised_estimates = cbkp$summarised,
+#' plots <- report_plots(summarised_estimates = cbkp$summarised,
 #'              reported = reported_cases)
 #'              
+#' plots$summary
+#' 
 #' # Pull out breakpoint summary
 #' cbkp$summarised[variable == "breakpoints"]
 #' 
@@ -205,9 +215,11 @@
 #'
 #' 
 #' # Plot output
-#' report_plots(summarised_estimates = fbkp$summarised,
+#' plots <- report_plots(summarised_estimates = fbkp$summarised,
 #'              reported = reported_cases)
 #'              
+#' plots$summary
+#' 
 #' # Pull out breakpoint summary
 #' fbkp$summarised[variable == "breakpoints"]
 #' 
diff --git a/man/estimate_infections.Rd b/man/estimate_infections.Rd
index da3b55b6d..9018cc592 100644
--- a/man/estimate_infections.Rd
+++ b/man/estimate_infections.Rd
@@ -181,8 +181,10 @@ snapshot <- estimate_infections(snapshot_cases, family = "negbin",
 
 
 # Plot output
-report_plots(summarised_estimates = snapshot$summarised,
-                      reported = snapshot_cases)     
+plots <- report_plots(summarised_estimates = snapshot$summarised,
+                      reported = snapshot_cases)
+
+plots$summary    
 
 ## Run model with stationary Rt assumption (likely to provide biased real-time estimates)
 stat <- estimate_infections(reported_cases, family = "negbin",
@@ -194,8 +196,10 @@ stat <- estimate_infections(reported_cases, family = "negbin",
 
 
 # Plot output
-report_plots(summarised_estimates = stat$summarised,
-             reported = reported_cases)
+plots <- report_plots(summarised_estimates = stat$summarised,
+                      reported = reported_cases)
+
+plots$summary
        
 # Run model with fixed Rt assumption 
 fixed <- estimate_infections(reported_cases, family = "negbin",
@@ -207,8 +211,10 @@ fixed <- estimate_infections(reported_cases, family = "negbin",
 
 
 # Plot output
-report_plots(summarised_estimates = fixed$summarised,
-             reported = reported_cases)
+plots <- report_plots(summarised_estimates = fixed$summarised,
+                      reported = reported_cases)
+                      
+plots$summary
 
 # Run model with breakpoints                                                                      
 bkp <- estimate_infections(reported_cases, family = "negbin",
@@ -220,8 +226,10 @@ bkp <- estimate_infections(reported_cases, family = "negbin",
 
 
 # Plot output
-report_plots(summarised_estimates = bkp$summarised,
-             reported = reported_cases)
+plots <- report_plots(summarised_estimates = bkp$summarised,
+                      reported = reported_cases)
+
+plots$summary
              
 # Run model with breakpoints but with constrained non-linear change over time 
 # This formulation may increase the apparent effect of the breakpoint but needs to be tested using
@@ -237,9 +245,11 @@ cbkp <- estimate_infections(reported_cases, family = "negbin",
 
 
 # Plot output
-report_plots(summarised_estimates = cbkp$summarised,
+plots <- report_plots(summarised_estimates = cbkp$summarised,
              reported = reported_cases)
              
+plots$summary
+
 # Pull out breakpoint summary
 cbkp$summarised[variable == "breakpoints"]
 
@@ -255,9 +265,11 @@ fbkp <- estimate_infections(reported_cases, family = "negbin",
 
 
 # Plot output
-report_plots(summarised_estimates = fbkp$summarised,
+plots <- report_plots(summarised_estimates = fbkp$summarised,
              reported = reported_cases)
              
+plots$summary
+
 # Pull out breakpoint summary
 fbkp$summarised[variable == "breakpoints"]