diff --git a/NAMESPACE b/NAMESPACE
index 4ba1e15d..aa29ccd8 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -24,6 +24,7 @@ export(set_carrying_capacity)
 export(set_clinical_treatment)
 export(set_demography)
 export(set_drugs)
+export(set_epi_outputs)
 export(set_equilibrium)
 export(set_mass_pev)
 export(set_mda)
diff --git a/R/output_parameters.R b/R/output_parameters.R
new file mode 100644
index 00000000..7e752ecf
--- /dev/null
+++ b/R/output_parameters.R
@@ -0,0 +1,39 @@
+#' @title Parameterise age grouped output rendering
+#'
+#' @details this function produces discrete and contiguous age groups, inclusive of the lower
+#' age limit and exclusive of the upper age limit: e.g., list(c(0, 10, 100), c(200, 250) will produce
+#' three age groups: 0-9, 10-99 and 200-249 in days.
+#' @param parameters the model parameters
+#' @param age_group age breaks for population size outputs; default = NULL
+#' @param incidence age breaks for incidence outputs (D+Tr+A); default = NULL
+#' @param clinical_incidence age breaks for clinical incidence outputs (symptomatic); default = c(0, 1825)
+#' @param severe_incidence age breaks for severe incidence outputs; default = NULL
+#' @param prevalence age breaks for clinical prevalence outputs (pcr and lm detectable infections); default = c(730, 3650)
+#' @export
+#'
+set_epi_outputs <- function(parameters,
+                            age_group = NULL,
+                            incidence = NULL,
+                            clinical_incidence = NULL,
+                            severe_incidence = NULL,
+                            prevalence = NULL,
+                            ica = NULL,
+                            icm = NULL,
+                            iva = NULL,
+                            ivm = NULL,
+                            id = NULL,
+                            ib = NULL
+){
+  
+  parent_formals <- names(formals())
+  parent_formals <- parent_formals[which(parent_formals != "parameters")]
+  outputs <- parent_formals[!unlist(lapply(parent_formals, function(x){is.null(get(x))}))]
+  
+  for (output in outputs) {
+    if(!is.list(get(output))){stop("Each age input must be a list of vectors")}
+    parameters[[paste0(output, "_rendering_min_ages")]] <- unlist(lapply(get(output), function(x){x[-length(x)]}))
+    parameters[[paste0(output, "_rendering_max_ages")]] <- unlist(lapply(get(output), function(x){x[-1]-1}))
+  }
+  
+  parameters
+}
diff --git a/R/processes.R b/R/processes.R
index 12f8c3db..985d2019 100644
--- a/R/processes.R
+++ b/R/processes.R
@@ -181,6 +181,9 @@ create_processes <- function(
   # =========
   # Rendering
   # =========
+  
+  imm_var_names <- c('ica', 'icm', 'id', 'ib', 'iva', 'ivm')
+  
   processes <- c(
     processes,
     individual::categorical_count_renderer_process(
@@ -190,10 +193,10 @@ create_processes <- function(
     ),
     create_variable_mean_renderer_process(
       renderer,
-      c('ica', 'icm', 'ib', 'id', 'iva', 'ivm'),
-      variables[c('ica', 'icm', 'ib', 'id', 'iva', 'ivm')]
+      imm_var_names,
+      variables[imm_var_names]
     ),
-    create_prevelance_renderer(
+    create_prevalence_renderer(
       variables$state,
       variables$birth,
       variables$id,
@@ -205,6 +208,13 @@ create_processes <- function(
       parameters,
       renderer
     ),
+    create_age_variable_mean_renderer_process(
+      imm_var_names[paste0(imm_var_names,"_rendering_min_ages") %in% names(parameters)],
+      variables[imm_var_names[paste0(imm_var_names,"_rendering_min_ages") %in% names(parameters)]],
+      variables$birth,
+      parameters,
+      renderer
+    ),
     create_compartmental_rendering_process(renderer, solvers, parameters)
   )
 
diff --git a/R/render.R b/R/render.R
index 60cb5b73..820f15c0 100644
--- a/R/render.R
+++ b/R/render.R
@@ -14,7 +14,7 @@ in_age_range <- function(birth, timestep, lower, upper) {
 #' @param renderer model renderer
 #' 
 #' @noRd
-create_prevelance_renderer <- function(
+create_prevalence_renderer <- function(
   state,
   birth,
   immunity,
@@ -119,6 +119,30 @@ create_variable_mean_renderer_process <- function(
   }
 }
 
+create_age_variable_mean_renderer_process <- function(
+    names,
+    variables,
+    birth,
+    parameters,
+    renderer
+) {
+  function(timestep) {
+    for (i in seq_along(variables)) {
+      for (j in seq_along(parameters[[paste0(names[[i]],"_rendering_min_ages")]])) {
+        lower <- parameters[[paste0(names[[i]],"_rendering_min_ages")]][[j]]
+        upper <- parameters[[paste0(names[[i]],"_rendering_max_ages")]][[j]]
+        in_age <- in_age_range(birth, timestep, lower, upper)
+        renderer$render(paste0('n_', lower, '_', upper), in_age$size(), timestep)
+        renderer$render(
+          paste0(names[[i]], '_mean_', lower, '_', upper),
+          mean(variables[[i]]$get_values(index = in_age)),
+          timestep
+        )
+      }
+    }
+  }
+}
+
 create_vector_count_renderer_individual <- function(
   mosquito_state,
   species,
diff --git a/man/set_epi_outputs.Rd b/man/set_epi_outputs.Rd
new file mode 100644
index 00000000..7e2b07ef
--- /dev/null
+++ b/man/set_epi_outputs.Rd
@@ -0,0 +1,42 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/output_parameters.R
+\name{set_epi_outputs}
+\alias{set_epi_outputs}
+\title{Parameterise age grouped output rendering}
+\usage{
+set_epi_outputs(
+  parameters,
+  age_group = NULL,
+  incidence = NULL,
+  clinical_incidence = NULL,
+  severe_incidence = NULL,
+  prevalence = NULL,
+  ica = NULL,
+  icm = NULL,
+  iva = NULL,
+  ivm = NULL,
+  id = NULL,
+  ib = NULL
+)
+}
+\arguments{
+\item{parameters}{the model parameters}
+
+\item{age_group}{age breaks for population size outputs; default = NULL}
+
+\item{incidence}{age breaks for incidence outputs (D+Tr+A); default = NULL}
+
+\item{clinical_incidence}{age breaks for clinical incidence outputs (symptomatic); default = c(0, 1825)}
+
+\item{severe_incidence}{age breaks for severe incidence outputs; default = NULL}
+
+\item{prevalence}{age breaks for clinical prevalence outputs (pcr and lm detectable infections); default = c(730, 3650)}
+}
+\description{
+Parameterise age grouped output rendering
+}
+\details{
+this function produces discrete and contiguous age groups, inclusive of the lower
+age limit and exclusive of the upper age limit: e.g., list(c(0, 10, 100), c(200, 250) will produce
+three age groups: 0-9, 10-99 and 200-249 in days.
+}
diff --git a/tests/testthat/test-output.R b/tests/testthat/test-output.R
new file mode 100644
index 00000000..ee221d76
--- /dev/null
+++ b/tests/testthat/test-output.R
@@ -0,0 +1,49 @@
+test_that('Test age parameter function works', {
+  parameters <- get_parameters()
+  age_limits <- c(0,1,2,3)*365
+  parameters <- set_epi_outputs(parameters,
+                                age_group = list(age_limits),
+                                incidence = list(age_limits+1),
+                                clinical_incidence = list(age_limits+2),
+                                severe_incidence = list(age_limits+3),
+                                prevalence = list(age_limits+4),
+                                ica = list(age_limits+5),
+                                icm = list(age_limits+6),
+                                id = list(age_limits+7),
+                                ib = list(age_limits+8),
+                                iva = list(age_limits+9),
+                                ivm = list(age_limits+10)
+  )
+  
+  sim <- run_simulation(timesteps = 1, parameters)
+  
+  expect_true(
+    all(
+      paste0(rep(c("n_age",
+                   "n", "n_inc", "p_inc",
+                   "n","n_inc_clinical","p_inc_clinical",
+                   "n","n_inc_severe","p_inc_severe",
+                   "n","n_detect","p_detect",
+                   "ica_mean", "icm_mean","id_mean","ib_mean","iva_mean","ivm_mean"), each = 3),"_",
+             age_limits[-4]+rep(c(0,rep(c(1:4), each = 3),5:10), each = 3),"_",age_limits[-1]-1+rep(c(0,rep(c(1:4), each = 3),5:10), each = 3)
+      ) %in%
+        names(sim)))
+  
+})
+
+test_that('Test age parameter multiple sequences function works', {
+  parameters <- get_parameters()
+  age_limits <- c(0,1,2,3)*365
+  parameters <- set_epi_outputs(parameters,
+                                age_group = list(age_limits, age_limits+100))
+  
+  sim <- run_simulation(timesteps = 1, parameters)
+  
+  expect_true(
+    all(
+      paste0(rep(c("n_age"), each = 3),"_",
+             c(age_limits[-4], age_limits[-4]+100), "_", c(age_limits[-1], age_limits[-1]+100)-1
+      ) %in%
+        names(sim)))
+  
+})
diff --git a/tests/testthat/test-render.R b/tests/testthat/test-render.R
index 920f09bb..07f256c3 100644
--- a/tests/testthat/test-render.R
+++ b/tests/testthat/test-render.R
@@ -11,7 +11,7 @@ test_that('that default rendering works', {
   immunity <- individual::DoubleVariable$new(rep(1, 4))
 
   renderer <- mock_render(1)
-  process <- create_prevelance_renderer(
+  process <- create_prevalence_renderer(
     state,
     birth,
     immunity,
@@ -62,7 +62,7 @@ test_that('that default rendering works when no one is in the age range', {
   immunity <- individual::DoubleVariable$new(rep(1, 4))
 
   renderer <- mock_render(1)
-  process <- create_prevelance_renderer(
+  process <- create_prevalence_renderer(
     state,
     birth,
     immunity,