diff --git a/lib/functions/create_tables.R b/lib/functions/create_tables.R
new file mode 100644
index 0000000..910fbb0
--- /dev/null
+++ b/lib/functions/create_tables.R
@@ -0,0 +1,131 @@
+# Function to create clinical pathways table
+create_clinical_pathways_table <- function(title) {
+  data <- tibble(
+    Condition = c(
+      "Uncomplicated Urinary Tract Infection",
+      "Shingles",
+      "Impetigo",
+      "Infected Insect Bites",
+      "Acute Sore Throat",
+      "Acute Sinusitis",
+      "Acute Otitis Media"
+    ),
+    Age = c(
+      "16 to 64 years",
+      "18 years and over",
+      "1 year and over",
+      "1 year and over",
+      "5 years and over",
+      "12 years and over",
+      "1 to 17 years"
+    ),
+    Sex = c(
+      "Female",
+      "Any",
+      "Any",
+      "Any",
+      "Any",
+      "Any",
+      "Any"
+    ),
+    Exclusions = c(
+      "Pregnant individuals, urinary catheter, recurrent UTI (2 episodes in last 6 months, or 3 episodes in last 12 months)",
+      "Pregnant individuals",
+      "Bullous impetigo, recurrent impetigo (2 or more episodes in the same year), pregnant individuals under 16 years",
+      "Pregnant individuals under 16 years",
+      "Pregnant individuals under 16 years",
+      "Immunosuppressed individuals, chronic sinusitis (symptoms lasting more than 12 weeks), pregnant individuals under 16 years",
+      "Recurrent acute otitis media (3 or more episodes in 6 months or four or more episodes in 12 months), pregnant individuals under 16 years"
+    )
+  )
+
+  data %>%
+    gt() %>%
+    tab_header(
+      title = title
+      # subtitle = "Inclusion and exclusion criteria for clinical pathway/conditions"
+    ) %>%
+    cols_label(
+      Condition = "Condition",
+      Age = "Age Range",
+      Sex = "Sex",
+      Exclusions = "Exclusions"
+    ) %>%
+    tab_options(
+      table.font.size = "medium",
+      heading.title.font.size = "large",
+      heading.subtitle.font.size = "small"
+    ) %>%
+    tab_style(
+      style = cell_text(weight = "bold"),
+      locations = cells_column_labels(columns = everything())
+    )
+}
+
+# Function to create pharmacy first service codes table
+create_pf_service_codes_table <- function(title) {
+  data <- tibble(
+    codelist = c(
+      "Community Pharmacist (CP) Consultation Service for minor illness (procedure)",
+      "Pharmacy First service (qualifier value)"
+    ),
+    code = c(
+      "1577041000000109",
+      "983341000000102"
+    )
+  )
+
+  data %>%
+    gt() %>%
+    tab_header(
+      title = title,
+      # subtitle = "Codelist descriptions and their respective SNOMED codes"
+    ) %>%
+    cols_label(
+      codelist = md("**Codelist Description**"),
+      code = md("**SNOMED Code**")
+    ) %>%
+    tab_options(
+      table.font.size = "medium",
+      heading.title.font.size = "large",
+      heading.subtitle.font.size = "small"
+    )
+}
+
+create_clinical_conditions_codes_table <- function(title) {
+  data <- tibble(
+    condition = c(
+      "Acute otitis media",
+      "Herpes zoster",
+      "Acute sinusitis",
+      "Impetigo",
+      "Infected insect bite",
+      "Acute pharyngitis",
+      "Uncomplicated urinary tract infection"
+    ),
+    code = c(
+      "3110003",
+      "4740000",
+      "15805002",
+      "48277006",
+      "262550002",
+      "363746003",
+      "1090711000000102"
+    )
+  )
+  data %>%
+    gt() %>%
+    tab_header(
+      title = title
+      # subtitle = "Clinical conditions and their corresponding SNOMED codes"
+    ) %>%
+    cols_label(
+      condition = md("**Clinical Condition**"),
+      code = md("**SNOMED Code**")
+    ) %>%
+    tab_options(
+      table.font.size = "medium",
+      heading.title.font.size = "large",
+      heading.subtitle.font.size = "small"
+    )
+}
diff --git a/lib/functions/load_opensafely_outputs.R b/lib/functions/load_opensafely_outputs.R
index 432ce74..cf60f1a 100644
--- a/lib/functions/load_opensafely_outputs.R
+++ b/lib/functions/load_opensafely_outputs.R
@@ -36,53 +36,9 @@ df_measures <- tidy_measures(
   pf_measures_groupby_dict = pf_measures_groupby_dict
 )
 
-df_measures$ethnicity <- factor(
-  df_measures$ethnicity,
-  levels = c(
-    "White",
-    "Mixed",
-    "Asian or Asian British",
-    "Black or Black British",
-    "Chinese or Other Ethnic Groups",
-    "Missing"
-  ),
-  ordered = TRUE
-)
-
-df_measures$age_band <- factor(
-  df_measures$age_band,
-  levels = c(
-    "0-19",
-    "20-39",
-    "40-59",
-    "60-79",
-    "80+",
-    "Missing"
-  ),
-  ordered = TRUE
-)
-
-df_measures$region <- factor(
-  df_measures$region,
-  levels = c(
-    "East",
-    "East Midlands",
-    "London",
-    "North East",
-    "North West",
-    "South East",
-    "South West",
-    "West Midlands",
-    "Yorkshire and The Humber",
-    "Missing"
-  ),
-  ordered = TRUE
-)
-
-df_measures <- df_measures %>%
-  mutate(sex = factor(sex,
-    levels = c("female", "male"),
-    labels = c("Female", "Male")
-  ))
+# str(df_measures$ethnicity)
+# str(df_measures$age_band)
+# str(df_measures$region)
+# str(df_measures$sex)
 
 df_measures$age_band[is.na(df_measures$age_band)] <- "Missing"
diff --git a/lib/functions/plot_measures.R b/lib/functions/plot_measures.R
index 4171c0c..b39d0d5 100644
--- a/lib/functions/plot_measures.R
+++ b/lib/functions/plot_measures.R
@@ -29,14 +29,20 @@ plot_measures <- function(
     facet_wrap = FALSE,
     facet_var = NULL,
     colour_var = NULL,
+    shape_var = NULL,
+    colour_palette = NULL,
+    y_scale = NULL,
+    scale_measure = NULL,
+    shapes = NULL,
+    date_breaks = "1 month",
     legend_position = "bottom") {
   # Test if all columns expected in output from generate measures exist
-  expected_names <- c("measure", "interval_start", "interval_end", "ratio", "numerator", "denominator")
-  missing_columns <- setdiff(expected_names, colnames(data))
+  # expected_names <- c("measure", "interval_start", "interval_end", "ratio", "numerator", "denominator")
+  # missing_columns <- setdiff(expected_names, colnames(data))
 
-  if (length(missing_columns) > 0) {
-    stop("Data does not have expected column(s): ", paste(missing_columns, collapse = ", "), call. = FALSE)
-  }
+  # if (length(missing_columns) > 0) {
+  #   stop("Data does not have expected column(s): ", paste(missing_columns, collapse = ", "), call. = FALSE)
+  # }
 
   plot_tmp <- ggplot(
     data,
@@ -44,11 +50,13 @@ plot_measures <- function(
       x = {{ select_interval_date }},
       y = {{ select_value }},
       colour = {{ colour_var }},
-      group = {{ colour_var }}
+      group = {{ colour_var }},
+      shape = {{ colour_var }},
+      fill = {{ colour_var }}
     )
   ) +
-    geom_point() +
-    geom_line(alpha = .5) +
+    geom_point(size = 2) +
+    geom_line(alpha = .3) +
     geom_vline(
       xintercept = lubridate::as_date("2024-02-01"),
       linetype = "dotted",
@@ -56,35 +64,63 @@ plot_measures <- function(
       linewidth = .7
     ) +
     scale_x_date(
-      date_breaks = "1 month",
+      date_breaks = {{ date_breaks }},
       labels = scales::label_date_short()
     ) +
     guides(
-      color = guide_legend(nrow = guide_nrow)
+      color = guide_legend(nrow = guide_nrow),
+      shape = guide_legend(nrow = guide_nrow)
     ) +
     labs(
       title = title,
       x = x_label,
       y = y_label,
       colour = guide_label,
+      shape = NULL,
+      fill = NULL
     ) +
     theme(
       legend.position = legend_position,
-      plot.title = element_text(hjust = 0.5)
+      plot.title = element_text(hjust = 0.5),
+      text = element_text(size = 14)
     )
 
+  # Change colour based on specified colour palette
+  if (!is.null(colour_palette)) {
+    if (length(colour_palette) == 1 && colour_palette == "plasma") {
+      plot_tmp <- plot_tmp + scale_colour_viridis_d(option = "plasma", end = .75) +
+        geom_line(size = 0.5) +
+        geom_point(size = 2.5)
+    } else {
+      plot_tmp <- plot_tmp + scale_colour_manual(values = colour_palette)
+    }
+  } else {
+    plot_tmp <- plot_tmp + scale_colour_viridis_d(end = .75)
+  }
+
+  if (!is.null(shapes) && shapes == "condition_shapes") {
+    plot_tmp <- plot_tmp + scale_shape_manual(values = condition_shapes)
+  }
+
   # Automatically change y scale depending selected value
-  if (rlang::as_label(enquo(select_value)) %in% c("numerator", "denominator")) {
+  scale_label <- rlang::as_label(enquo(scale_measure))
+  if (is.null(scale_measure)) {
     plot_tmp <- plot_tmp + scale_y_continuous(
       limits = c(0, NA),
       labels = scales::label_number()
     )
-  } else {
+  } else if (scale_measure == "rate") {
     plot_tmp <- plot_tmp + scale_y_continuous(
       limits = c(0, NA),
-      # scale = 1000 to calculate rate per 1000 people
       labels = scales::label_number(scale = 1000)
     )
+  } else if (scale_measure == "percent") {
+    plot_tmp <- plot_tmp + scale_y_continuous(labels = scales::percent)
+  } else {
+    plot_tmp <- plot_tmp + scale_y_continuous(
+      limits = c(0, NA),
+      labels = scales::label_number()
+    )
   }
 
   # Add facets if requested
@@ -94,12 +130,51 @@ plot_measures <- function(
     plot_tmp <- plot_tmp +
       facet_wrap(vars({{ facet_var }}), ncol = 2)
   }
+  # Add y_scale to add option for free_y
+  if (!is.null(y_scale) && y_scale == "free_y") {
+    plot_tmp <- plot_tmp +
+      facet_wrap(~source, scales = "free_y")
+  }
 
   plot_tmp
 }
 
+set_patchwork_theme <- function(patchwork_figure) {
+  patchwork_figure +
+    plot_annotation(tag_levels = "A") +
+    plot_layout(guides = "collect", widths = c(2, 1)) &
+    theme(
+      legend.position = "bottom",
+      text = element_text(size = 15),
+      strip.background = element_rect(size = 0),
+      strip.text.x = element_text(size = 13, face = "bold")
+    )
+}
+
+save_figure <- function(figure, width = 10, height = 6) {
+  # this uses the 'figure' argument as a string to later generate a filename
+  figure_name <- deparse(substitute(figure))
+  ggsave(
+    filename = here("released_output", "results", "figures", paste(figure_name, "png",sep = ".")),
+    figure,
+    width = width, height = height
+  )
+}
+
 # Colour palettes
 gradient_palette <- c("#001F4D", "#0056B3", "#007BFF", "#66B3E2", "#A4D8E1", "grey")
 region_palette <- c("red", "navy", "#018701", "#ffa600ca", "purple", "brown", "#f4a5b2", "cyan", "green", "grey")
 ethnicity_palette <- c("#42db0188", "#0056B3", "#ff0000c2", "#a52a2a5a", "purple", "grey")
 sex_palette <- c("red", "blue")
+dark2_palette <- RColorBrewer::brewer.pal(n = 8, name = "Dark2")
+
+# Custom shapes
+condition_shapes <- c(
+  "Acute Sinusitis" = 15,
+  "Infected Insect Bite" = 19,
+  "UTI" = 4,
+  "Acute Otitis Media" = 23,
+  "Acute Pharyngitis" = 3,
+  "Herpes Zoster" = 17,
+  "Impetigo" = 8
+)
diff --git a/lib/functions/tidy_measures.R b/lib/functions/tidy_measures.R
index 44ab4e6..6d1ff6f 100644
--- a/lib/functions/tidy_measures.R
+++ b/lib/functions/tidy_measures.R
@@ -33,6 +33,37 @@ pf_measures_groupby_dict <- list(
   ethnicity = "Ethnicity"
 )
 
+pf_measures_ethnicity_list <- list(
+  "White",
+  "Mixed",
+  "Asian or Asian British",
+  "Black or Black British",
+  "Chinese or Other Ethnic Groups",
+  "Missing"
+)
+
+pf_measures_age_list <- list(
+  "0-19",
+  "20-39",
+  "40-59",
+  "60-79",
+  "80+",
+  "Missing"
+)
+
+pf_measures_region_list <- list(
+    "East",
+    "East Midlands",
+    "London",
+    "North East",
+    "North West",
+    "South East",
+    "South West",
+    "West Midlands",
+    "Yorkshire and The Humber",
+    "Missing"
+)
+
 #' Tidy measures data
 #'
 #' Creates a tidier dataframe of measures data.
@@ -49,14 +80,22 @@ pf_measures_groupby_dict <- list(
 #' @return A dataframe
 tidy_measures <- function(data, pf_measures_name_dict, pf_measures_name_mapping, pf_measures_groupby_dict) {
   data_tmp <- data %>%
+  # Separate 'measure' column into 'summary_stat_measure' and 'group_by'
+  # Separate 'summary_stat_measure' into 'summary_stat' and 'measure'
     separate(measure, into = c("summary_stat_measure", "group_by"), sep = "_by_") %>%
     separate(summary_stat_measure, into = c("summary_stat", "measure"), sep = "_", extra = "merge")
 
+  # Modify columns based on recoding and factor levels
   data_tmp <- data_tmp %>%
     mutate(
+      # Recode 'measure' to be more readable
       measure_desc = recode(factor(measure), !!!pf_measures_name_mapping),
       measure = recode(factor(measure), !!!pf_measures_name_dict),
-      group_by = recode(factor(group_by), !!!pf_measures_groupby_dict)
+      group_by = recode(factor(group_by), !!!pf_measures_groupby_dict),
+      ethnicity = factor(ethnicity, levels = pf_measures_ethnicity_list, labels = pf_measures_ethnicity_list),
+      age_band = factor(age_band, levels = pf_measures_age_list, labels = pf_measures_age_list),
+      region = factor(region, levels = pf_measures_region_list, labels = pf_measures_region_list),
+      sex = factor(sex, levels = c("female", "male"), labels = c("Female", "Male"))
     )
 
   data_tmp
diff --git a/reports/create_figures.Rmd b/reports/create_figures.Rmd
new file mode 100644
index 0000000..818d6a1
--- /dev/null
+++ b/reports/create_figures.Rmd
@@ -0,0 +1,795 @@
+---
+title: "Pharmacy First"
+output:
+  html_document:
+    toc: true
+    toc_depth: 4
+  pdf_document: default
+date: "`r format(Sys.time(), '%d %B, %Y')`"
+---
+
+```{r setup, include=FALSE}
+knitr::opts_chunk$set(echo = FALSE)
+library(tidyverse)
+library(here)
+library(readr)
+library(gt)
+library(patchwork)
+```
+
+```{r load-data, message=FALSE, warning=FALSE}
+# Load functions
+source(here("lib", "functions", "tidy_measures.R"))
+source(here("lib", "functions", "plot_measures.R"))
+
+# Load validation data:
+# - df_bsa_medication_validation: date, pharmacy_advanced_service, bnf_paragraph, count
+# - df_bsa_consultation_validation: date, consultation_type, source, count_method, count
+source(here("lib", "functions", "load_validation_data.R"))
+
+# Load opensafely ouputs:
+# - df_measures: measure, interval_start, interval_end, ratio numerator, denominator, age_band, sex,imd, region, ethnicity
+# - df_descriptive_stats: measure, interval_start, interval_end, ratio numerator, denominator
+# - df_pfmed: measure, interval_start, interval_end, ratio, numerator, denominator, dmd_code
+# - df_condition_provider: measure, interval_start, interval_end, ratio, numerator, denominator, pf_status, imd
+source(here("lib", "functions", "load_opensafely_outputs.R"))
+```
+
+```{r, message=FALSE, warning=FALSE}
+# Create figure for total count of Pharmacy First consultations for each code (3 codes)
+df_measures_selected <- df_measures %>%
+  filter(measure_desc == "clinical_service") %>%
+  filter(is.na(group_by)) |>
+  select(measure, interval_start, numerator) |>
+  mutate(measure = factor(measure,
+    levels = c("Consultation Service", "Pharmacy First Consultation"),
+    labels = c(
+      "Consultation Service for minor illness (1577041000000109)",
+      "Pharmacy First service (983341000000102)"
+    )
+  ))
+
+fig_pf_individual_consultations_count <- plot_measures(
+  df_measures_selected,
+  select_value = numerator,
+  select_interval_date = interval_start,
+  legend_position = "bottom",
+  facet_wrap = FALSE,
+  facet_var = measure,
+  y_label = "Total Count",
+  colour_var = measure,
+  guide_nrow = 1,
+)
+
+save_figure(fig_pf_individual_consultations_count)
+
+# Create figure for total count of Pharmacy First Consultations (GROUPED)
+df_measures_selected <- df_measures_selected |>
+  group_by(interval_start) |>
+  mutate(
+    pf_consultation_total = sum(numerator, na.rm = TRUE),
+    data_desc = "Pharmacy First Consultation"
+  )
+
+fig_pf_grouped_consultations_count <- plot_measures(
+  df_measures_selected,
+  select_value = pf_consultation_total,
+  select_interval_date = interval_start,
+  legend_position = "bottom",
+  facet_wrap = FALSE,
+  facet_var = data_desc,
+  y_label = "Total Count",
+  colour_var = data_desc,
+  guide_nrow = 1,
+)
+
+save_figure(fig_pf_grouped_consultations_count)
+
+```
+
+```{r, message=FALSE, warning=FALSE, fig.height=10, fig.width=8}
+# Create figure for total count of Pharmacy First grouped conditions (no breakdowns)
+df_measures_selected <- df_measures %>%
+  filter(measure_desc == "clinical_condition") %>%
+  filter(is.na(group_by))
+
+# Create visualisation
+fig_pf_grouped_conditions_count <- plot_measures(
+  df_measures_selected,
+  select_value = numerator,
+  select_interval_date = interval_end,
+  guide_nrow = 1,
+  facet_wrap = TRUE,
+  facet_var = measure,
+  title = "Pharmacy First Conditions",
+  y_label = "Number of codes for PF conditions",
+)
+
+save_figure(fig_pf_grouped_conditions_count)
+```
+
+### Breakdown by age
+
+```{r, message=FALSE, warning=FALSE, fig.height=4, fig.width=8}
+# Create figure for total count of PF consultations by age
+df_measures_selected <- df_measures %>%
+  filter(measure_desc == "pharmacy_first_services") %>%
+  filter(group_by == "Age band")
+
+# Create visualisation
+fig_pf_consultations_by_age_count <- plot_measures(
+  df_measures_selected,
+  select_value = numerator,
+  select_interval_date = interval_end,
+  colour_var = age_band,
+  guide_nrow = 1,
+  facet_wrap = TRUE,
+  facet_var = measure,
+  title = "Pharmacy First Consultations",
+  y_label = "Number of codes for PF consultations",
+  colour_palette = gradient_palette
+)
+
+save_figure(fig_pf_consultations_by_age_count)
+
+# Create figure for rate of PF consultations by age
+df_measures_selected <- df_measures %>%
+  filter(measure_desc == "pharmacy_first_services") %>%
+  filter(group_by == "Age band")
+
+# Create visualisation
+fig_pf_consultations_by_age_rate <- plot_measures(
+  df_measures_selected,
+  select_value = ratio,
+  select_interval_date = interval_end,
+  colour_var = age_band,
+  guide_nrow = 1,
+  facet_wrap = TRUE,
+  facet_var = measure,
+  scale_measure = "rate",
+  title = "Rate of Pharmacy First Consultations per 1000 people",
+  y_label = "Number of codes for FP consultations",
+  colour_palette = gradient_palette
+)
+
+save_figure(fig_pf_consultations_by_age_rate)
+```
+
+```{r, message=FALSE, warning=FALSE, fig.height=8, fig.width=8}
+# Create figure for total count of PF clinical conditions broken down by age
+df_measures_selected <- df_measures %>%
+  filter(measure_desc == "clinical_condition") %>%
+  filter(group_by == "Age band")
+
+# Create visualisation
+fig_pf_conditions_by_age_count <- plot_measures(
+  df_measures_selected,
+  select_value = numerator,
+  select_interval_date = interval_end,
+  colour_var = age_band,
+  guide_nrow = 1,
+  facet_wrap = TRUE,
+  facet_var = measure,
+  title = "Pharmacy First Conditions",
+  y_label = "Number of codes for PF conditions",
+  colour_palette = gradient_palette
+)
+
+save_figure(fig_pf_conditions_by_age_count)
+
+# Create figure for rate of PF clinical conditions broken down by age
+df_measures_selected <- df_measures %>%
+  filter(measure_desc == "clinical_condition") %>%
+  filter(group_by == "Age band")
+
+# Create visualisation
+fig_pf_conditions_by_age_rate <- plot_measures(
+  df_measures_selected,
+  select_value = ratio,
+  select_interval_date = interval_end,
+  colour_var = age_band,
+  guide_nrow = 1,
+  facet_wrap = TRUE,
+  facet_var = measure,
+  scale_measure = "rate",
+  title = "Rate of Pharmacy First Conditions per 1000 people",
+  y_label = "Number of codes for PF conditions",
+  colour_palette = gradient_palette
+)
+
+save_figure(fig_pf_conditions_by_age_rate)
+
+```
+
+### Breakdown by sex
+
+```{r, message=FALSE, warning=FALSE, fig.height=4, fig.width=8}
+# Create figure for total count of PF consultations broken down by sex
+df_measures_selected <- df_measures %>%
+  filter(measure_desc == "pharmacy_first_services") %>%
+  filter(group_by == "Sex")
+
+# Create visualisation
+fig_pf_consultations_by_sex_count <- plot_measures(
+  df_measures_selected,
+  select_value = numerator,
+  select_interval_date = interval_end,
+  colour_var = sex,
+  guide_nrow = 1,
+  facet_wrap = TRUE,
+  facet_var = measure,
+  title = "Pharmacy First Consultations",
+  y_label = "Number of codes for FP consultations",
+  colour_palette = sex_palette
+)
+
+save_figure(fig_pf_consultations_by_sex_count)
+
+# Create figure for rate of PF consultations broken down by sex
+df_measures_selected <- df_measures %>%
+  filter(measure_desc == "pharmacy_first_services") %>%
+  filter(group_by == "Sex")
+
+# Create visualisation
+fig_pf_consultations_by_sex_rate <- plot_measures(
+  df_measures_selected,
+  select_value = ratio,
+  select_interval_date = interval_end,
+  colour_var = sex,
+  guide_nrow = 1,
+  facet_wrap = TRUE,
+  facet_var = measure,
+  scale_measure = "rate",
+  title = "Rate of Pharmacy First Consultations per 1000 people",
+  y_label = "Number of codes for FP consultations",
+  colour_palette = sex_palette
+)
+
+save_figure(fig_pf_consultations_by_sex_rate)
+```
+
+```{r, message=FALSE, warning=FALSE, fig.height=8, fig.width=8}
+# Create figure for total count of PF clinical conditions by sex
+df_measures_selected <- df_measures %>%
+  filter(measure_desc == "clinical_condition") %>%
+  filter(group_by == "Sex")
+
+# Create visualisation
+fig_pf_conditions_by_sex_count <- plot_measures(
+  df_measures_selected,
+  select_value = numerator,
+  select_interval_date = interval_end,
+  colour_var = sex,
+  guide_nrow = 1,
+  facet_wrap = TRUE,
+  facet_var = measure,
+  title = "Pharmacy First Conditions",
+  y_label = "Number of codes for PF conditions",
+  colour_palette = sex_palette
+)
+
+save_figure(fig_pf_conditions_by_sex_count)
+
+# Create figure for rate of PF clinical conditions by sex
+df_measures_selected <- df_measures %>%
+  filter(measure_desc == "clinical_condition") %>%
+  filter(group_by == "Sex")
+
+# Create visualisation
+fig_pf_conditions_by_sex_rate <- plot_measures(
+  df_measures_selected,
+  select_value = ratio,
+  select_interval_date = interval_end,
+  colour_var = sex,
+  guide_nrow = 1,
+  facet_wrap = TRUE,
+  facet_var = measure,
+  scale_measure = "rate",
+  title = "Rate of Pharmacy First Conditions per 1000 people",
+  y_label = "Number of codes for PF conditions",
+) + scale_color_manual(values = sex_palette)
+
+save_figure(fig_pf_conditions_by_sex_rate)
+```
+
+### Breakdown by IMD
+
+```{r, message=FALSE, warning=FALSE, fig.height=4, fig.width=8}
+# Create figure for total count of PF consultations by IMD
+df_measures_selected <- df_measures %>%
+  filter(measure_desc == "pharmacy_first_services") %>%
+  filter(group_by == "IMD")
+
+# Create visualisation
+fig_pf_consultations_by_imd_count <- plot_measures(
+  df_measures_selected,
+  select_value = numerator,
+  select_interval_date = interval_end,
+  colour_var = imd,
+  guide_nrow = 1,
+  facet_wrap = TRUE,
+  facet_var = measure,
+  title = "Pharmacy First Consultations",
+  y_label = "Number of codes for FP consultations",
+  colour_palette = gradient_palette
+)
+
+save_figure(fig_pf_consultations_by_imd_count)
+
+# Create figure for rate of PF consultations by IMD
+df_measures_selected <- df_measures %>%
+  filter(measure_desc == "pharmacy_first_services") %>%
+  filter(group_by == "IMD")
+
+# Create visualisation
+fig_pf_consultations_by_imd_rate <- plot_measures(
+  df_measures_selected,
+  select_value = ratio,
+  select_interval_date = interval_end,
+  colour_var = imd,
+  guide_nrow = 1,
+  facet_wrap = TRUE,
+  facet_var = measure,
+  scale_measure = "rate",
+  title = "Rate of Pharmacy First Consultations per 1000 people",
+  y_label = "Number of codes for FP consultations",
+  colour_palette = gradient_palette
+)
+
+save_figure(fig_pf_consultations_by_imd_rate)
+```
+
+```{r, message=FALSE, warning=FALSE, fig.height=8, fig.width=8}
+# Create figure for total count of PF clinical conditions by IMD
+df_measures_selected <- df_measures %>%
+  filter(measure_desc == "clinical_condition") %>%
+  filter(group_by == "IMD")
+
+# Create visualisation
+fig_pf_conditions_by_imd_count <- plot_measures(
+  df_measures_selected,
+  select_value = numerator,
+  select_interval_date = interval_end,
+  colour_var = imd,
+  guide_nrow = 1,
+  facet_wrap = TRUE,
+  facet_var = measure,
+  title = "Pharmacy First Conditions",
+  y_label = "Number of codes for PF conditions",
+  colour_palette = gradient_palette
+)
+
+save_figure(fig_pf_conditions_by_imd_count)
+
+# Create figure for rate of PF clinical conditions by IMD
+df_measures_selected <- df_measures %>%
+  filter(measure_desc == "clinical_condition") %>%
+  filter(group_by == "IMD")
+
+# Create visualisation
+fig_pf_conditions_by_imd_rate <- plot_measures(
+  df_measures_selected,
+  select_value = ratio,
+  select_interval_date = interval_end,
+  colour_var = imd,
+  guide_nrow = 1,
+  facet_wrap = TRUE,
+  facet_var = measure,
+  scale_measure = "rate",
+  title = "Rate of Pharmacy First Conditions per 1000 people",
+  y_label = "Number of codes for PF conditions",
+  colour_palette = gradient_palette
+)
+
+save_figure(fig_pf_conditions_by_imd_rate)
+
+```
+
+### Breakdown by region
+
+```{r, message=FALSE, warning=FALSE, fig.height=4, fig.width=8}
+# Create figure for total count of PF consultations by region
+df_measures_selected <- df_measures %>%
+  filter(measure_desc == "pharmacy_first_services") %>%
+  filter(group_by == "Region")
+
+# Create visualisation
+fig_pf_consultations_by_region_count <- plot_measures(
+  df_measures_selected,
+  select_value = numerator,
+  select_interval_date = interval_end,
+  colour_var = region,
+  guide_nrow = 2,
+  facet_wrap = TRUE,
+  facet_var = measure,
+  title = "Pharmacy First Consultations",
+  y_label = "Number of codes for FP consultations",
+  colour_palette = region_palette
+)
+
+save_figure(fig_pf_consultations_by_region_count)
+
+# Create figure for rate of PF consultations by IMD
+df_measures_selected <- df_measures %>%
+  filter(measure_desc == "pharmacy_first_services") %>%
+  filter(group_by == "Region")
+
+# Create visualisation
+fig_pf_consultations_by_region_rate <- plot_measures(
+  df_measures_selected,
+  select_value = ratio,
+  select_interval_date = interval_end,
+  colour_var = region,
+  guide_nrow = 2,
+  facet_wrap = TRUE,
+  facet_var = measure,
+  scale_measure = "rate",
+  title = "Rate of Pharmacy First Consultations per 1000 people",
+  y_label = "Number of codes for FP consultations",
+  colour_palette = region_palette
+)
+
+save_figure(fig_pf_consultations_by_region_rate)
+
+```
+
+```{r, message=FALSE, warning=FALSE, fig.height=8, fig.width=8}
+# Create figure for total count of PF clinical conditions by region
+df_measures_selected <- df_measures %>%
+  filter(measure_desc == "clinical_condition") %>%
+  filter(group_by == "Region")
+
+# Create visualisation
+fig_pf_conditions_by_region_count <- plot_measures(
+  df_measures_selected,
+  select_value = numerator,
+  select_interval_date = interval_end,
+  colour_var = region,
+  guide_nrow = 2,
+  facet_wrap = TRUE,
+  facet_var = measure,
+  title = "Pharmacy First Conditions",
+  y_label = "Number of codes for PF conditions",
+  colour_palette = region_palette
+)
+
+save_figure(fig_pf_conditions_by_region_count)
+
+# Create figure for rate of PF clinical conditions by region
+df_measures_selected <- df_measures %>%
+  filter(measure_desc == "clinical_condition") %>%
+  filter(group_by == "Region")
+
+# Create visualisation
+fig_pf_conditions_by_region_rate <- plot_measures(
+  df_measures_selected,
+  select_value = ratio,
+  select_interval_date = interval_end,
+  colour_var = region,
+  guide_nrow = 2,
+  facet_wrap = TRUE,
+  facet_var = measure,
+  scale_measure = "rate",
+  title = "Rate of Pharmacy First Conditions per 1000 people",
+  y_label = "Number of codes for PF conditions",
+  colour_palette = region_palette
+)
+
+save_figure(fig_pf_conditions_by_region_rate)
+
+```
+
+### Breakdown by ethnicity
+
+```{r, message=FALSE, warning=FALSE, fig.height=4, fig.width=8}
+# Create figure for total count of PF consultations by ethnicity
+df_measures_selected <- df_measures %>%
+  filter(measure_desc == "pharmacy_first_services") %>%
+  filter(group_by == "Ethnicity")
+
+# Create visualisation
+fig_pf_consultations_by_ethnicity_count <- plot_measures(
+  df_measures_selected,
+  select_value = numerator,
+  select_interval_date = interval_end,
+  colour_var = ethnicity,
+  guide_nrow = 2,
+  facet_wrap = TRUE,
+  facet_var = measure,
+  title = "Pharmacy First Consultations",
+  y_label = "Number of codes for FP consultations",
+  colour_palette = ethnicity_palette
+)
+
+save_figure(fig_pf_consultations_by_ethnicity_count)
+
+# Create figure for rate of PF consultations by ethnicity
+df_measures_selected <- df_measures %>%
+  filter(measure_desc == "pharmacy_first_services") %>%
+  filter(group_by == "Ethnicity")
+
+# Create visualisation
+fig_pf_consultations_by_ethnicity_rate <- plot_measures(
+  df_measures_selected,
+  select_value = ratio,
+  select_interval_date = interval_end,
+  colour_var = ethnicity,
+  guide_nrow = 2,
+  facet_wrap = TRUE,
+  facet_var = measure,
+  scale_measure = "rate",
+  title = "Rate of Pharmacy First Consultations per 1000 people",
+  y_label = "Number of codes for FP consultations",
+  colour_palette = ethnicity_palette
+)
+
+save_figure(fig_pf_consultations_by_ethnicity_rate)
+```
+
+```{r, message=FALSE, warning=FALSE, fig.height=8, fig.width=8}
+# Create figure for total count of PF clinical conditions by ethnicity
+df_measures_selected <- df_measures %>%
+  filter(measure_desc == "clinical_condition") %>%
+  filter(group_by == "Ethnicity")
+
+# Create visualisation
+fig_pf_conditions_by_ethnicity_count <- plot_measures(
+  df_measures_selected,
+  select_value = numerator,
+  select_interval_date = interval_end,
+  colour_var = ethnicity,
+  guide_nrow = 2,
+  facet_wrap = TRUE,
+  facet_var = measure,
+  title = "Pharmacy First Conditions",
+  y_label = "Number of codes for PF conditions",
+  colour_palette = ethnicity_palette
+)
+
+save_figure(fig_pf_conditions_by_ethnicity_count)
+
+# Create figure for rate of PF clinical conditions by ethnicity
+df_measures_selected <- df_measures %>%
+  filter(measure_desc == "clinical_condition") %>%
+  filter(group_by == "Ethnicity")
+
+# Create visualisation
+fig_pf_conditions_by_ethnicity_rate <- plot_measures(
+  df_measures_selected,
+  select_value = ratio,
+  select_interval_date = interval_end,
+  colour_var = ethnicity,
+  guide_nrow = 2,
+  facet_wrap = TRUE,
+  facet_var = measure,
+  scale_measure = "rate",
+  title = "Rate of Pharmacy First Conditions per 1000 people",
+  y_label = "Number of codes for PF conditions",
+  colour_palette = ethnicity_palette
+)
+
+save_figure(fig_pf_conditions_by_ethnicity_rate)
+
+```
+```{r, message=FALSE, warning=FALSE, echo = FALSE}
+# Create figure to compare OS and BSA counts for PF clinical conditions
+
+# OpenSAFELY data for clinical conditions into a tidy df
+df_opensafely_validation <- df_measures %>%
+  filter(measure_desc == "clinical_condition") %>%
+  # filter(interval_start >= as.Date("2024-02-01") & interval_start <= as.Date("2024-07-30")) %>%
+  filter(is.na(group_by)) %>%
+  select(date = interval_start, consultation_type = measure, count = numerator) %>%
+  mutate(
+    source = "opensafely",
+    count_method = "opensafely_tpp"
+  ) |>
+  filter(date >= "2024-01-01") %>%
+  relocate(date, consultation_type, source, count_method, count)
+
+# Combining rows from OS and BSA validation dataframes
+df_validation_condition <- bind_rows(df_opensafely_validation, df_bsa_consultation_validation)
+
+# Line graph comparing clinical condition counts of BSA and OS data
+df_validation_condition_counts <- df_validation_condition %>%
+  filter(count_method %in% c("opensafely_tpp", "count_40pct")) %>%
+  filter(date >= "2024-01-01") %>%
+  mutate(source = factor(source,
+    levels = c("opensafely", "nhs_bsa"),
+    labels = c("OpenSAFELY-TPP", "NHS BSA (40%)")
+  ))
+
+# Create visualisation
+fig_validation_condition_count <- plot_measures(
+  df_validation_condition_counts,
+  select_value = count,
+  select_interval_date = date,
+  colour_var = consultation_type,
+  guide_nrow = 2,
+  facet_wrap = TRUE,
+  facet_var = source,
+  y_label = "Count",
+  y_scale = "free_y",
+  shapes = "condition_shapes",
+  colour_palette = "plasma",
+  date_breaks = "2 month"
+)
+
+# Another plot visualising the percentage
+df_validation_condition_pct <- df_validation_condition %>%
+  filter(count_method %in% c("opensafely_tpp", "count_40pct")) %>%
+  filter(date >= "2024-01-01") %>%
+  pivot_wider(names_from = c(source, count_method), values_from = count) %>%
+  mutate(source = "Percentage of NHS BSA (40%) in OpenSAFELY")
+
+fig_validation_condition_pct <- plot_measures(
+  df_validation_condition_pct,
+  select_value = opensafely_opensafely_tpp / nhs_bsa_count_40pct,
+  select_interval_date = date,
+  colour_var = consultation_type,
+  guide_nrow = 2,
+  facet_wrap = TRUE,
+  facet_var = source,
+  scale_measure = "percent",
+  y_label = "Percent",
+  y_scale = "free_y",
+  shapes = "condition_shapes",
+  colour_palette = "plasma",
+  date_breaks = "2 month"
+)
+
+fig_validation_condition_comparison <- (fig_validation_condition_count + fig_validation_condition_pct) %>%
+set_patchwork_theme()
+save_figure(fig_validation_condition_comparison, width = 15)
+
+```
+
+```{r, message=FALSE, warning=FALSE, echo = FALSE, fig.width=8}
+# Create figure to show & of PF Med, Condition and both with linked PF consultations
+
+df_descriptive_stats <- df_descriptive_stats %>%
+  mutate(
+    measure = factor(measure,
+      levels = c("pf_with_pfmed", "pf_with_pfcondition", "pf_with_pfmed_and_pfcondition"),
+      labels = c("PF Med", "PF Condition", "PF Med & PF Condition")
+    )
+  )
+
+fig_pf_descriptive_stats <- plot_measures(
+  df_descriptive_stats,
+  select_value = ratio,
+  select_interval_date = interval_start,
+  colour_var = measure,
+  guide_nrow = 2,
+  facet_wrap = FALSE,
+  facet_var = measure,
+  scale_measure = "percent",
+  y_label = "Percent",
+  colour_palette = dark2_palette,
+  date_breaks = "1 month"
+)
+
+save_figure(fig_pf_descriptive_stats)
+```
+
+```{r, message=FALSE, warning=FALSE, echo = FALSE}
+# Create figure to compare OS and BSA counts for PF medication
+
+df_bsa_medication_validation_sum <- df_bsa_medication_validation %>%
+  group_by(date) %>%
+  summarise(count = sum(count) * 0.4) %>%
+  mutate(
+    source = "nhs_bsa",
+    count_method = "count_40pct"
+  )
+range(df_pfmed$interval_start)
+df_opensafely_pfmed_sum <- df_pfmed %>%
+  rename(date = interval_start) %>%
+  group_by(date) %>%
+  summarise(count = sum(numerator)) %>%
+  mutate(
+    source = "opensafely_tpp",
+    count_method = "opensafely_tpp"
+  )
+
+df_validation_med_counts <- bind_rows(df_opensafely_pfmed_sum, df_bsa_medication_validation_sum) |>
+  filter(date >= "2024-01-01" & date <= "2024-07-01")
+
+df_validation_med_counts <- df_validation_med_counts %>%
+  mutate(
+    source = factor(source, levels = c("opensafely_tpp", "nhs_bsa"), labels = c("OpenSAFELY-TPP", "NHS BSA")),
+    count_method = factor(count_method, levels = c("opensafely_tpp", "count_40pct"), labels = c("OpenSAFELY-TPP", "NHS BSA (40%)"))
+  )
+
+fig_validation_med_count <- plot_measures(
+  df_validation_med_counts,
+  select_value = count,
+  select_interval_date = date,
+  colour_var = count_method,
+  guide_nrow = 1,
+  facet_wrap = TRUE,
+  facet_var = source,
+  y_scale = "free_y",
+  y_label = "Count",
+  colour_palette = "plasma",
+  date_breaks = "1 month"
+)
+
+# Another plot visualising the percentage
+df_validation_med_pct <- df_validation_med_counts %>%
+  filter(count_method %in% c("OpenSAFELY-TPP", "NHS BSA (40%)")) %>%
+  pivot_wider(names_from = c(source, count_method), values_from = count) %>%
+  mutate(source = "Percentage of NHS BSA (40%) in OpenSAFELY")
+
+fig_validation_med_pct <- plot_measures(
+  df_validation_med_pct,
+  select_value = `OpenSAFELY-TPP_OpenSAFELY-TPP` / `NHS BSA_NHS BSA (40%)`,
+  select_interval_date = date,
+  colour_var = source,
+  guide_nrow = 1,
+  facet_wrap = TRUE,
+  facet_var = source,
+  scale_measure = "percent",
+  y_scale = "free_y",
+  y_label = "Count",
+  colour_palette = "plasma",
+  date_breaks = "1 month"
+)
+
+fig_validation_med_comparison <- (fig_validation_med_count + fig_validation_med_pct) %>% 
+set_patchwork_theme()
+save_figure(fig_validation_med_comparison, width = 15)
+
+```
+```{r, message=FALSE, warning=FALSE, echo = FALSE, fig.width=8}
+# Create figure to compare clinical events linked to PF consultation and not linked
+
+df_condition_provider_grouped <- df_condition_provider %>%
+  group_by(measure, interval_start, pf_status) %>%
+  summarise(
+    count = sum(numerator)
+  ) %>%
+  mutate(
+    measure = factor(measure,
+      levels = c(
+        "count_acute_sinusitis_total",
+        "count_infected_insect_bite_total",
+        "count_uncomplicated_urinary_tract_infection_total",
+        "count_acute_otitis_media_total",
+        "count_acute_pharyngitis_total",
+        "count_herpes_zoster_total",
+        "count_impetigo_total"
+      ),
+      labels = c(
+        "Acute Sinusitis",
+        "Infected Insect Bite",
+        "UTI",
+        "Acute Otitis Media",
+        "Acute Pharyngitis",
+        "Herpes Zoster",
+        "Impetigo"
+      )
+    ),
+    pf_status = factor(pf_status,
+      levels = c(TRUE, FALSE),
+      labels = c("Linked to Pharmacy First consultation", "Not linked to Pharmacy First consultation")
+    )
+  )
+
+fig_pf_condition_provider_count <- plot_measures(
+  df_condition_provider_grouped,
+  select_value = count,
+  select_interval_date = interval_start,
+  colour_var = pf_status,
+  guide_nrow = 1,
+  facet_wrap = TRUE,
+  facet_var = measure,
+  y_label = "Count",
+  date_breaks = "6 month",
+  colour_palette = "plasma"
+)
+
+save_figure(fig_pf_condition_provider_count)
+```
+
+# References
diff --git a/reports/create_tables.Rmd b/reports/create_tables.Rmd
new file mode 100644
index 0000000..38b795f
--- /dev/null
+++ b/reports/create_tables.Rmd
@@ -0,0 +1,96 @@
+```{r setup, include=FALSE}
+knitr::opts_chunk$set(echo = FALSE)
+library(tidyverse)
+library(here)
+library(readr)
+library(gt)
+library(patchwork)
+```
+
+```{r echo=FALSE, message=FALSE}
+source(here("lib", "functions", "create_tables.R"))
+```
+
+```{r echo=FALSE, message=FALSE}
+# Create clinical pathways table
+tab_clinical_pathways <- create_clinical_pathways_table("Table 1. Pharmacy First population criteria")
+gtsave(tab_clinical_pathways, filename = here("released_output", "results", "tables", "tab_pf_condition_criteria.png"))
+```
+
+```{r echo=FALSE}
+# Create pharmacy first service codes dataframe
+tab_pf_service_codes <- create_pf_service_codes_table("Table 2. Pharmacy First consultation codes")
+gtsave(tab_pf_service_codes, filename = here("released_output", "results", "tables", "tab_pf_service_codelist.png"))
+```
+
+```{r echo=FALSE, message=FALSE}
+tab_pf_condition_codes <- create_clinical_conditions_codes_table("Table 3. Pharmacy First condition codes")
+gtsave(tab_pf_condition_codes, filename = here("released_output", "results", "tables", "tab_pf_condition_codes.png"))
+```
+
+```{r, message=FALSE, warning=FALSE, echo = FALSE}
+
+df_consultation_med_counts <- read_csv(
+  here("released_output", "measures", "consultation_med_counts_measures.csv"),
+  col_types = cols(dmd_code = col_character())
+)
+
+codelist_file_names <- list.files(
+  here("codelists"),
+  pattern = "\\.csv$",
+  full.names = FALSE
+)
+
+pf_med_code_desc <- here("codelists", codelist_file_names) |>
+  map(~ read_csv(.x, col_types = cols(
+    code = col_character()
+  ))) |>
+  bind_rows() |>
+  select(code, term) |>
+  distinct()
+
+top10_nonpf_med_code_desc  <- tribble(
+  ~code, ~term,
+  "37388111000001102", "Macrobid 100mg modified-release capsules",
+  "531611000001106", "Difflam 0.15% spray",
+  "3697711000001104", "EarCalm 2% spray",
+  "42533911000001101", "Coryen 27.5micrograms/dose nasal spray",
+  "623911000001105", "Hydrocortisone 1% cream",
+  "4530711000001104", "Covonia Sore Throat 0.2%/0.05% oromucosal spray menthol",
+  "42100111000001106", "Benzydamine 0.15% oromucosal spray sugar free",
+  "29311000001104", "Piriton 2mg/5ml syrup",
+  "17290311000001107", "Benzydamine 0.15% oromucosal spray sugar free",
+  "4648111000001108", "Robitussin Chesty Cough 100mg/5ml oral solution"
+)
+
+combined_med_code_desc <- pf_med_code_desc |>
+  dplyr::bind_rows(top10_nonpf_med_code_desc)
+
+df_pf_med_counts <- df_consultation_med_counts |>
+  select(numerator, code = dmd_code, pharmacy_first_med) |>
+  left_join(combined_med_code_desc, by = "code") |>
+  filter(numerator > 0) |>
+  select(-code) %>%
+  group_by(pharmacy_first_med, term) |>
+  summarise(count = sum(numerator, na.rm = TRUE)) |>
+  filter(!is.na(term)) %>%
+  ungroup() |>
+  group_by(pharmacy_first_med) |>
+  mutate(ratio_by_group = count / sum(count, na.na.rm = TRUE)) |>
+  slice_max(order_by = ratio_by_group, n = 5) |>
+  ungroup()
+
+df_pf_and_non_pf_med_counts <- df_pf_med_counts %>%
+  arrange(!pharmacy_first_med) %>%
+  mutate(pharmacy_first_med = factor(pharmacy_first_med,
+  levels = c(FALSE, TRUE),
+  labels = c(("Medication not included in codelists"), "Medication included in codelists"))) %>%
+  group_by(pharmacy_first_med)
+  
+gt_top_meds(df_pf_and_non_pf_med_counts)
+
+gtsave(
+  tab_pf_med_counts,
+  here("released_output", "results", "tables", "tab_pf_med_counts.png"),
+)
+```
\ No newline at end of file
diff --git a/reports/pharmacy_first_report.Rmd b/reports/pharmacy_first_report.Rmd
deleted file mode 100644
index da3e951..0000000
--- a/reports/pharmacy_first_report.Rmd
+++ /dev/null
@@ -1,1151 +0,0 @@
----
-title: "Pharmacy First"
-output:
-  html_document:
-    toc: true
-    toc_depth: 4
-  pdf_document: default
-date: "`r format(Sys.time(), '%d %B, %Y')`"
-bibliography: references.bib
-link-citations: true
----
-
-```{r setup, include=FALSE}
-knitr::opts_chunk$set(echo = FALSE)
-library(tidyverse)
-library(here)
-library(readr)
-library(gt)
-library(patchwork)
-```
-
-```{r load-data, message=FALSE, warning=FALSE}
-# Load functions
-source(here("lib", "functions", "tidy_measures.R"))
-source(here("lib", "functions", "plot_measures.R"))
-
-# Load validation data:
-# - df_bsa_medication_validation: date, pharmacy_advanced_service, bnf_paragraph, count
-# - df_bsa_consultation_validation: date, consultation_type, source, count_method, count
-source(here("lib", "functions", "load_validation_data.R"))
-
-# Load opensafely ouputs:
-# - df_measures: measure, interval_start, interval_end, ratio numerator, denominator, age_band, sex,imd, region, ethnicity
-# - df_descriptive_stats: measure, interval_start, interval_end, ratio numerator, denominator
-# - df_pfmed: measure, interval_start, interval_end, ratio, numerator, denominator, dmd_code
-# - df_condition_provider: measure, interval_start, interval_end, ratio, numerator, denominator, pf_status, imd
-source(here("lib", "functions", "load_opensafely_outputs.R"))
-```
-
-# Background
-
-Add background here.
-
-# Methods
-
-## Data source
-
-Primary care records managed by the GP software provider TPP were accessed through OpenSAFELY (https://opensafely.org).
-OpenSAFELY-TPP covers approximately 40% of the population of England, for a detailed description of the representativeness see @Andrews2022.
-
-Between 1st November 2023 and 1st September 2024, we identified individuals aged between 0 and 120 years who were registered at a TPP practice each month across the study period. 
-From this cohort, we counted the number of patients with at least one SNOMED CT code indicating a consultation using the Pharmacy First Service. 
-Patients with more than one instance of the same clinical code within a month were counted only once.
-
-Our study relied on the correct entry of relevant Pharmacy First codes in patients' GP records. 
-The Pharmacy First service uses [GP Connect - Update Record](https://digital.nhs.uk/services/gp-connect/gp-connect-in-your-organisation/gp-connect-update-record) to update a patient's GP record with consultation information from community pharmacies. 
-Following the launch of the Pharmacy First service, there has been a [Gradual roll-out of GP Connect - Update Record](https://cpe.org.uk/our-news/gp-connect-update-record-rollout-and-flow-of-information/) across approved community pharmacy IT system suppliers.
-
-### Population
-
-The eligible patient population for each clinical condition associated with Pharmacy First is detailed in the table below.
-Currently, we have not yet applied any of the inclusion or exlusion citeria to restrict the population used in this report.
-This will initially help us to understand the underlying data.
-
-```{r echo=FALSE, message=FALSE}
-# Create clinical pathways dataframe
-clinical_pathways_table <- data.frame(
-  Condition = c(
-    "Uncomplicated Urinary Tract Infection",
-    "Shingles",
-    "Impetigo",
-    "Infected Insect Bites",
-    "Acute Sore Throat",
-    "Acute Sinusitis",
-    "Acute Otitis Media"
-  ),
-  Age = c(
-    "16 to 64 years",
-    "18 years and over",
-    "1 year and over",
-    "1 year and over",
-    "5 years and over",
-    "12 years and over",
-    "1 to 17 years"
-  ),
-  Sex = c(
-    "Female",
-    "Any",
-    "Any",
-    "Any",
-    "Any",
-    "Any",
-    "Any"
-  ),
-  Exclusions = c(
-    "Pregnant individuals, urinary catheter, recurrent UTI (2 episodes in last 6 months, or 3 episodes in last 12 months)",
-    "Pregnant individuals",
-    "Bullous impetigo, recurrent impetigo (2 or more episodes in the same year), pregnant individuals under 16 years",
-    "Pregnant individuals under 16 years",
-    "Pregnant individuals under 16 years",
-    "Immunosuppressed individuals, chronic sinusitis (symptoms lasting more than 12 weeks), pregnant individuals under 16 years",
-    "Recurrent acute otitis media (3 or more episodes in 6 months or four or more episodes in 12 months), pregnant individuals under 16 years"
-  )
-)
-
-# Create clinical pathways table
-clinical_pathways_table %>%
-  gt() %>%
-  tab_header(
-    title = "Table 1. Pharmacy First population criteria"
-    # subtitle = "Inclusion and exclusion criteria for clinical pathway/conditions"
-  ) %>%
-  cols_label(
-    Condition = "Condition",
-    Age = "Age Range",
-    Sex = "Sex",
-    Exclusions = "Exclusions"
-  ) %>%
-  tab_options(
-    table.font.size = "medium",
-    heading.title.font.size = "large",
-    heading.subtitle.font.size = "small"
-  ) %>%
-  tab_style(
-    style = cell_text(weight = "bold"),
-    locations = cells_column_labels(columns = everything())
-  )
-```
-
-### Codelists
-
-We used the following codelists to identify Pharmacy First consultations, conditions, and demographic breakdowns.
-
-#### Pharmacy First consultation codes
-
-The following two SNOMED codes were used to identify Pharmacy First consultations.
-For clarity, we combined these codes for the presentation of the results.
-
-```{r echo=FALSE}
-# Create pharmacy first service codes dataframe
-pharmacy_first_table <- data.frame(
-  codelist = c(
-    "Community Pharmacist (CP) Consultation Service for minor illness (procedure)",
-    "Pharmacy First service (qualifier value)"
-  ),
-  code = c(
-    "1577041000000109",
-    "983341000000102"
-  )
-)
-
-# Create pharmacy first service codes table
-pharmacy_first_table %>%
-  gt() %>%
-  tab_header(
-    title = "Table 2. Pharmacy First consultation codes",
-    # subtitle = "Codelist descriptions and their respective SNOMED codes"
-  ) %>%
-  cols_label(
-    codelist = md("**Codelist Description**"),
-    code = md("**SNOMED Code**")
-  ) %>%
-  tab_options(
-    table.font.size = "medium",
-    heading.title.font.size = "large",
-    heading.subtitle.font.size = "small"
-  )
-```
-
-#### Pharmacy First condition codes
-
-To categorise clinical events related to Pharmacy First services used the Pharmacy First [Clinical Pathways Codelist](https://www.opencodelists.org/codelist/opensafely/pharmacy-first-clinical-pathway-conditions/7ec97762/#full-list).
-
-```{r echo=FALSE, message=FALSE}
-clinical_codes_table <- data.frame(
-  condition = c(
-    "Acute otitis media",
-    "Herpes zoster",
-    "Acute sinusitis",
-    "Impetigo",
-    "Infected insect bite",
-    "Acute pharyngitis",
-    "Uncomplicated urinary tract infection"
-  ),
-  code = c(
-    "3110003",
-    "4740000",
-    "15805002",
-    "48277006",
-    "262550002",
-    "363746003",
-    "1090711000000102"
-  )
-)
-
-clinical_codes_table %>%
-  gt() %>%
-  tab_header(
-    title = "Table 3. Pharmacy First condition codes"
-    # subtitle = "Clinical conditions and their corresponding SNOMED codes"
-  ) %>%
-  cols_label(
-    condition = md("**Clinical Condition**"),
-    code = md("**SNOMED Code**")
-  ) %>%
-  tab_options(
-    table.font.size = "medium",
-    heading.title.font.size = "large",
-    heading.subtitle.font.size = "small"
-  )
-```
-
-#### Pregnancy Codelist
-
-The [Pregnancy Codelist](https://www.opencodelists.org/codelist/nhsd-primary-care-domain-refsets/preg_cod/20200812/#full-list) was used to identify patients who were pregnant during each month.
-
-#### Ethnicity Codelist
-
-We used the [Ethnicity Codelist](https://www.opencodelists.org/codelist/opensafely/ethnicity-snomed-0removed/2e641f61/) identify ethnicity in Electronic Health Records.
-To ensure comprehensive ethnicity data, we supplemented missing ethnicity values with data from the Secondary Uses Service (SUS).
-
-# Results
-
-### Total population
-
-```{r, message=FALSE, warning=FALSE, fig.height=6, fig.width=10}
-# Select measures and breakdown
-df_measures_selected <- df_measures %>%
-  filter(measure_desc == "clinical_service") %>%
-  filter(is.na(group_by)) |>
-  select(measure, interval_start, numerator) |>
-  mutate(measure = factor(measure,
-    levels = c("Consultation Service", "Pharmacy First Consultation"),
-    labels = c(
-      "Consultation Service for minor illness (1577041000000109)",
-      "Pharmacy First service (983341000000102)"
-    )
-  ))
-
-fig_pf_individual_consultations_count <- df_measures_selected |>
-  select(measure, interval_start, numerator) |>
-  ggplot(aes(
-    x = interval_start,
-    y = numerator,
-    colour = measure,
-    shape = measure,
-  )) +
-  geom_point(size = 2) +
-  geom_line(alpha = .3) +
-  labs(
-    title = NULL,
-    x = NULL,
-    y = "Total count",
-    colour = NULL,
-    shape = NULL
-  ) +
-  scale_y_continuous(
-    labels = scales::label_number(),
-  ) +
-  theme(legend.position = "bottom") +
-  guides(
-    colour = guide_legend(ncol = 2),
-    shape = guide_legend(ncol = 2)
-  ) +
-  scale_x_date(
-    date_breaks = "1 month",
-    labels = scales::label_date_short()
-  ) +
-  geom_vline(
-    xintercept = lubridate::as_date(c(
-      "2024-01-31"
-    )),
-    linetype = "dotted",
-    colour = "orange",
-    size = .7
-  ) +
-  scale_colour_viridis_d(end = .75) +
-  theme(
-    text = element_text(size = 14)
-  )
-
-
-ggsave(
-  here("released_output", "results", "figures", "fig_pf_individual_consultations_count.png"),
-  fig_pf_individual_consultations_count,
-  width = 10, height = 6
-)
-
-fig_pf_individual_consultations_count
-
-fig_pf_grouped_consultations_count <- df_measures_selected |>
-  group_by(interval_start) |>
-  mutate(
-    pf_consultation_total = sum(numerator, na.rm = TRUE),
-    data_desc = "Pharmacy First Consultation"
-  ) |>
-  ggplot(aes(
-    x = interval_start,
-    y = pf_consultation_total,
-    colour = data_desc,
-    shape = data_desc,
-  )) +
-  geom_point(size = 2) +
-  geom_line(alpha = .3) +
-  labs(
-    title = NULL,
-    x = NULL,
-    y = "Total count",
-    colour = NULL,
-    shape = NULL
-  ) +
-  scale_y_continuous(
-    labels = scales::label_number(),
-  ) +
-  theme(legend.position = "bottom") +
-  guides(
-    colour = guide_legend(ncol = 2),
-    shape = guide_legend(ncol = 2)
-  ) +
-  scale_x_date(
-    date_breaks = "1 month",
-    labels = scales::label_date_short()
-  ) +
-  geom_vline(
-    xintercept = lubridate::as_date(c(
-      "2024-01-31"
-    )),
-    linetype = "dotted",
-    colour = "orange",
-    size = .7
-  ) +
-  scale_colour_viridis_d(end = .75) +
-  theme(
-    text = element_text(size = 14)
-  )
-
-
-ggsave(
-  here("released_output", "results", "figures", "fig_pf_grouped_consultations_count.png"),
-  fig_pf_grouped_consultations_count,
-  width = 10, height = 6
-)
-```
-
-```{r, message=FALSE, warning=FALSE, fig.height=10, fig.width=8}
-# Select measures and breakdown
-df_measures_selected <- df_measures %>%
-  filter(measure_desc == "clinical_condition") %>%
-  filter(is.na(group_by))
-
-# Create visualisation
-plot_measures(
-  df_measures_selected,
-  select_value = numerator,
-  select_interval_date = interval_end,
-  guide_nrow = 1,
-  facet_wrap = TRUE,
-  facet_var = measure,
-  title = "Pharmacy First Conditions",
-  y_label = "Number of codes for PF conditions",
-)
-```
-
-### Breakdown by age
-
-```{r, message=FALSE, warning=FALSE, fig.height=4, fig.width=8}
-# Select measures and breakdown
-df_measures_selected <- df_measures %>%
-  filter(measure_desc == "pharmacy_first_services") %>%
-  filter(group_by == "Age band")
-
-# Create visualisation
-plot_measures(
-  df_measures_selected,
-  select_value = numerator,
-  select_interval_date = interval_end,
-  colour_var = age_band,
-  guide_nrow = 1,
-  facet_wrap = TRUE,
-  facet_var = measure,
-  title = "Pharmacy First Consultations",
-  y_label = "Number of codes for FP consultations",
-) + scale_color_manual(values = gradient_palette)
-
-# Select measures and breakdown
-df_measures_selected <- df_measures %>%
-  filter(measure_desc == "pharmacy_first_services") %>%
-  filter(group_by == "Age band")
-
-# Create visualisation
-plot_measures(
-  df_measures_selected,
-  select_value = ratio,
-  select_interval_date = interval_end,
-  colour_var = age_band,
-  guide_nrow = 1,
-  facet_wrap = TRUE,
-  facet_var = measure,
-  title = "Rate of Pharmacy First Consultations per 1000 people",
-  y_label = "Number of codes for FP consultations",
-) + scale_color_manual(values = gradient_palette)
-```
-
-```{r, message=FALSE, warning=FALSE, fig.height=8, fig.width=8}
-# Select measures and breakdown
-df_measures_selected <- df_measures %>%
-  filter(measure_desc == "clinical_condition") %>%
-  filter(group_by == "Age band")
-
-# Create visualisation
-plot_measures(
-  df_measures_selected,
-  select_value = numerator,
-  select_interval_date = interval_end,
-  colour_var = age_band,
-  guide_nrow = 1,
-  facet_wrap = TRUE,
-  facet_var = measure,
-  title = "Pharmacy First Conditions",
-  y_label = "Number of codes for PF conditions"
-) + scale_color_manual(values = gradient_palette)
-
-# Select measures and breakdown
-df_measures_selected <- df_measures %>%
-  filter(measure_desc == "clinical_condition") %>%
-  filter(group_by == "Age band")
-
-# Create visualisation
-plot_measures(
-  df_measures_selected,
-  select_value = ratio,
-  select_interval_date = interval_end,
-  colour_var = age_band,
-  guide_nrow = 1,
-  facet_wrap = TRUE,
-  facet_var = measure,
-  title = "Rate of Pharmacy First Conditions per 1000 people",
-  y_label = "Number of codes for PF conditions"
-) + scale_color_manual(values = gradient_palette)
-```
-
-### Breakdown by sex
-
-```{r, message=FALSE, warning=FALSE, fig.height=4, fig.width=8}
-# Select measures and breakdown
-df_measures_selected <- df_measures %>%
-  filter(measure_desc == "pharmacy_first_services") %>%
-  filter(group_by == "Sex")
-
-# Create visualisation
-plot_measures(
-  df_measures_selected,
-  select_value = numerator,
-  select_interval_date = interval_end,
-  colour_var = sex,
-  guide_nrow = 1,
-  facet_wrap = TRUE,
-  facet_var = measure,
-  title = "Pharmacy First Consultations",
-  y_label = "Number of codes for FP consultations",
-) + scale_color_manual(values = sex_palette)
-
-# Select measures and breakdown
-df_measures_selected <- df_measures %>%
-  filter(measure_desc == "pharmacy_first_services") %>%
-  filter(group_by == "Sex")
-
-# Create visualisation
-plot_measures(
-  df_measures_selected,
-  select_value = ratio,
-  select_interval_date = interval_end,
-  colour_var = sex,
-  guide_nrow = 1,
-  facet_wrap = TRUE,
-  facet_var = measure,
-  title = "Rate of Pharmacy First Consultations per 1000 people",
-  y_label = "Number of codes for FP consultations",
-) + scale_color_manual(values = sex_palette)
-```
-
-```{r, message=FALSE, warning=FALSE, fig.height=8, fig.width=8}
-# Select measures and breakdown
-df_measures_selected <- df_measures %>%
-  filter(measure_desc == "clinical_condition") %>%
-  filter(group_by == "Sex")
-
-# Create visualisation
-plot_measures(
-  df_measures_selected,
-  select_value = numerator,
-  select_interval_date = interval_end,
-  colour_var = sex,
-  guide_nrow = 1,
-  facet_wrap = TRUE,
-  facet_var = measure,
-  title = "Pharmacy First Conditions",
-  y_label = "Number of codes for PF conditions"
-) + scale_color_manual(values = sex_palette)
-
-# Select measures and breakdown
-df_measures_selected <- df_measures %>%
-  filter(measure_desc == "clinical_condition") %>%
-  filter(group_by == "Sex")
-
-# Create visualisation
-plot_measures(
-  df_measures_selected,
-  select_value = ratio,
-  select_interval_date = interval_end,
-  colour_var = sex,
-  guide_nrow = 1,
-  facet_wrap = TRUE,
-  facet_var = measure,
-  title = "Rate of Pharmacy First Conditions per 1000 people",
-  y_label = "Number of codes for PF conditions"
-) + scale_color_manual(values = sex_palette)
-```
-
-### Breakdown by IMD
-
-```{r, message=FALSE, warning=FALSE, fig.height=4, fig.width=8}
-# Select measures and breakdown
-df_measures_selected <- df_measures %>%
-  filter(measure_desc == "pharmacy_first_services") %>%
-  filter(group_by == "IMD")
-
-# Create visualisation
-plot_measures(
-  df_measures_selected,
-  select_value = numerator,
-  select_interval_date = interval_end,
-  colour_var = imd,
-  guide_nrow = 1,
-  facet_wrap = TRUE,
-  facet_var = measure,
-  title = "Pharmacy First Consultations",
-  y_label = "Number of codes for FP consultations",
-) + scale_color_manual(values = gradient_palette)
-
-# Select measures and breakdown
-df_measures_selected <- df_measures %>%
-  filter(measure_desc == "pharmacy_first_services") %>%
-  filter(group_by == "IMD")
-
-# Create visualisation
-plot_measures(
-  df_measures_selected,
-  select_value = ratio,
-  select_interval_date = interval_end,
-  colour_var = imd,
-  guide_nrow = 1,
-  facet_wrap = TRUE,
-  facet_var = measure,
-  title = "Rate of Pharmacy First Consultations per 1000 people",
-  y_label = "Number of codes for FP consultations",
-) + scale_color_manual(values = gradient_palette)
-```
-
-```{r, message=FALSE, warning=FALSE, fig.height=8, fig.width=8}
-# Select measures and breakdown
-df_measures_selected <- df_measures %>%
-  filter(measure_desc == "clinical_condition") %>%
-  filter(group_by == "IMD")
-
-# Create visualisation
-plot_measures(
-  df_measures_selected,
-  select_value = numerator,
-  select_interval_date = interval_end,
-  colour_var = imd,
-  guide_nrow = 1,
-  facet_wrap = TRUE,
-  facet_var = measure,
-  title = "Pharmacy First Conditions",
-  y_label = "Number of codes for PF conditions"
-) + scale_color_manual(values = gradient_palette)
-
-# Select measures and breakdown
-df_measures_selected <- df_measures %>%
-  filter(measure_desc == "clinical_condition") %>%
-  filter(group_by == "IMD")
-
-# Create visualisation
-plot_measures(
-  df_measures_selected,
-  select_value = ratio,
-  select_interval_date = interval_end,
-  colour_var = imd,
-  guide_nrow = 1,
-  facet_wrap = TRUE,
-  facet_var = measure,
-  title = "Rate of Pharmacy First Conditions per 1000 people",
-  y_label = "Number of codes for PF conditions"
-) + scale_color_manual(values = gradient_palette)
-```
-
-### Breakdown by region
-
-```{r, message=FALSE, warning=FALSE, fig.height=4, fig.width=8}
-# Select measures and breakdown
-df_measures_selected <- df_measures %>%
-  filter(measure_desc == "pharmacy_first_services") %>%
-  filter(group_by == "Region")
-
-# Create visualisation
-plot_measures(
-  df_measures_selected,
-  select_value = numerator,
-  select_interval_date = interval_end,
-  colour_var = region,
-  guide_nrow = 2,
-  facet_wrap = TRUE,
-  facet_var = measure,
-  title = "Pharmacy First Consultations",
-  y_label = "Number of codes for FP consultations",
-) + scale_color_manual(values = region_palette)
-
-# Select measures and breakdown
-df_measures_selected <- df_measures %>%
-  filter(measure_desc == "pharmacy_first_services") %>%
-  filter(group_by == "Region")
-
-# Create visualisation
-plot_measures(
-  df_measures_selected,
-  select_value = ratio,
-  select_interval_date = interval_end,
-  colour_var = region,
-  guide_nrow = 2,
-  facet_wrap = TRUE,
-  facet_var = measure,
-  title = "Rate of Pharmacy First Consultations per 1000 people",
-  y_label = "Number of codes for FP consultations",
-) + scale_color_manual(values = region_palette)
-```
-
-```{r, message=FALSE, warning=FALSE, fig.height=8, fig.width=8}
-# Select measures and breakdown
-df_measures_selected <- df_measures %>%
-  filter(measure_desc == "clinical_condition") %>%
-  filter(group_by == "Region")
-
-# Create visualisation
-plot_measures(
-  df_measures_selected,
-  select_value = numerator,
-  select_interval_date = interval_end,
-  colour_var = region,
-  guide_nrow = 2,
-  facet_wrap = TRUE,
-  facet_var = measure,
-  title = "Pharmacy First Conditions",
-  y_label = "Number of codes for PF conditions"
-) + scale_color_manual(values = region_palette)
-
-# Select measures and breakdown
-df_measures_selected <- df_measures %>%
-  filter(measure_desc == "clinical_condition") %>%
-  filter(group_by == "Region")
-
-# Create visualisation
-plot_measures(
-  df_measures_selected,
-  select_value = ratio,
-  select_interval_date = interval_end,
-  colour_var = region,
-  guide_nrow = 2,
-  facet_wrap = TRUE,
-  facet_var = measure,
-  title = "Rate of Pharmacy First Conditions per 1000 people",
-  y_label = "Number of codes for PF conditions"
-) + scale_color_manual(values = region_palette)
-```
-
-### Breakdown by ethnicity
-
-```{r, message=FALSE, warning=FALSE, fig.height=4, fig.width=8}
-# Select measures and breakdown
-df_measures_selected <- df_measures %>%
-  filter(measure_desc == "pharmacy_first_services") %>%
-  filter(group_by == "Ethnicity")
-
-# Create visualisation
-plot_measures(
-  df_measures_selected,
-  select_value = numerator,
-  select_interval_date = interval_end,
-  colour_var = ethnicity,
-  guide_nrow = 2,
-  facet_wrap = TRUE,
-  facet_var = measure,
-  title = "Pharmacy First Consultations",
-  y_label = "Number of codes for FP consultations",
-) + scale_color_manual(values = ethnicity_palette)
-
-# Select measures and breakdown
-df_measures_selected <- df_measures %>%
-  filter(measure_desc == "pharmacy_first_services") %>%
-  filter(group_by == "Ethnicity")
-
-# Create visualisation
-plot_measures(
-  df_measures_selected,
-  select_value = ratio,
-  select_interval_date = interval_end,
-  colour_var = ethnicity,
-  guide_nrow = 2,
-  facet_wrap = TRUE,
-  facet_var = measure,
-  title = "Rate of Pharmacy First Consultations per 1000 people",
-  y_label = "Number of codes for FP consultations",
-) + scale_color_manual(values = ethnicity_palette)
-```
-
-```{r, message=FALSE, warning=FALSE, fig.height=8, fig.width=8}
-# Select measures and breakdown
-df_measures_selected <- df_measures %>%
-  filter(measure_desc == "clinical_condition") %>%
-  filter(group_by == "Ethnicity")
-
-# Create visualisation
-plot_measures(
-  df_measures_selected,
-  select_value = numerator,
-  select_interval_date = interval_end,
-  colour_var = ethnicity,
-  guide_nrow = 2,
-  facet_wrap = TRUE,
-  facet_var = measure,
-  title = "Pharmacy First Conditions",
-  y_label = "Number of codes for PF conditions"
-) + scale_color_manual(values = ethnicity_palette)
-
-# Select measures and breakdown
-df_measures_selected <- df_measures %>%
-  filter(measure_desc == "clinical_condition") %>%
-  filter(group_by == "Ethnicity")
-
-# Create visualisation
-plot_measures(
-  df_measures_selected,
-  select_value = ratio,
-  select_interval_date = interval_end,
-  colour_var = ethnicity,
-  guide_nrow = 2,
-  facet_wrap = TRUE,
-  facet_var = measure,
-  title = "Rate of Pharmacy First Conditions per 1000 people",
-  y_label = "Number of codes for PF conditions"
-) + scale_color_manual(values = ethnicity_palette)
-```
-```{r, message=FALSE, warning=FALSE, echo = FALSE}
-# OpenSAFELY data for clinical conditions into a tidy df
-df_opensafely_validation <- df_measures %>%
-  filter(measure_desc == "clinical_condition") %>%
-  # filter(interval_start >= as.Date("2024-02-01") & interval_start <= as.Date("2024-07-30")) %>%
-  filter(is.na(group_by)) %>%
-  select(date = interval_start, consultation_type = measure, count = numerator) %>%
-  mutate(
-    source = "opensafely",
-    count_method = "opensafely_tpp"
-  ) |>
-  filter(date >= "2024-01-01") %>%
-  relocate(date, consultation_type, source, count_method, count)
-
-# Combining rows from OS and BSA validation dataframes
-df_validation_condition_counts <- bind_rows(df_opensafely_validation, df_bsa_consultation_validation)
-
-# Line graph comparing clinical condition counts of BSA and OS data
-fig_validation_condition_count <- df_validation_condition_counts %>%
-  filter(count_method %in% c("opensafely_tpp", "count_40pct")) %>%
-  mutate(source = factor(source,
-    levels = c("opensafely", "nhs_bsa"),
-    labels = c("OpenSAFELY-TPP", "NHS BSA (40%)")
-  )) %>%
-  ggplot(
-    aes(
-      x = date,
-      y = count,
-      shape = consultation_type,
-      color = consultation_type,
-      fill = consultation_type,
-      group = consultation_type
-    )
-  ) +
-  geom_point(size = 2.5) +
-  geom_line(size = 0.5) +
-  facet_wrap(~source, scales = "free_y") +
-  scale_x_date(
-    labels = scales::label_date_short()
-  ) +
-  labs(x = NULL, y = "Count", colour = NULL, shape = NULL, fill = NULL) +
-  scale_color_viridis_d(
-    option = "plasma",
-    end = 0.9
-  ) +
-  scale_fill_viridis_d(
-    option = "plasma",
-    end = 0.9
-  ) +
-  scale_shape_manual(
-    values = c(
-      "Acute Sinusitis" = 15,
-      "Infected Insect Bite" = 19,
-      "UTI" = 4,
-      "Acute Otitis Media" = 23,
-      "Acute Pharyngitis" = 3,
-      "Herpes Zoster" = 17,
-      "Impetigo" = 8
-    )
-  ) +
-  theme(
-    text = element_text(size = 14)
-  ) +
-  geom_vline(
-    xintercept = lubridate::as_date("2024-02-01"),
-    linetype = "dotted",
-    colour = "orange",
-    linewidth = .7
-  ) +
-  scale_y_continuous(labels = scales::number)
-
-# Another plot visualising the percentage
-fig_validation_condition_pct <- df_validation_condition_counts %>%
-  filter(count_method %in% c("opensafely_tpp", "count_40pct")) %>%
-  pivot_wider(names_from = c(source, count_method), values_from = count) %>%
-  mutate(source = "Percentage of NHS BSA (40%) in OpenSAFELY") %>%
-  ggplot(
-    aes(
-      x = date,
-      y = opensafely_opensafely_tpp / nhs_bsa_count_40pct,
-      shape = consultation_type,
-      color = consultation_type,
-      fill = consultation_type,
-      group = consultation_type
-    )
-  ) +
-  geom_point(size = 2.5) +
-  geom_line(size = 0.5) +
-  facet_wrap(~source, scales = "free_y") +
-  scale_x_date(
-    labels = scales::label_date_short()
-  ) +
-  labs(x = NULL, y = "Percent", colour = NULL, shape = NULL, fill = NULL) +
-  scale_color_viridis_d(
-    option = "plasma",
-    end = 0.9
-  ) +
-  scale_fill_viridis_d(
-    option = "plasma",
-    end = 0.9
-  ) +
-  scale_shape_manual(
-    values = c(
-      "Acute Sinusitis" = 15,
-      "Infected Insect Bite" = 19,
-      "UTI" = 4,
-      "Acute Otitis Media" = 23,
-      "Acute Pharyngitis" = 3,
-      "Herpes Zoster" = 17,
-      "Impetigo" = 8
-    )
-  ) +
-  theme(
-    text = element_text(size = 14)
-  ) +
-  geom_vline(
-    xintercept = lubridate::as_date("2024-02-01"),
-    linetype = "dotted",
-    colour = "orange",
-    linewidth = .7
-  ) +
-  scale_y_continuous(labels = scales::percent)
-
-# Combining the plots with patchwork
-fig_validation_condition_count_pct <- (fig_validation_condition_count + fig_validation_condition_pct) +
-  plot_annotation(tag_levels = "A") +
-  plot_layout(guides = "collect", widths = c(2, 1)) &
-  theme(
-    legend.position = "bottom",
-    text = element_text(size = 15),
-    strip.background = element_rect(size = 0),
-    strip.text.x = element_text(size = 13, face = "bold")
-  )
-
-fig_validation_condition_count_pct
-
-ggsave(
-  here("released_output", "results", "figures", "fig_validation_condition_count_pct.png"),
-  fig_validation_condition_count_pct,
-  width = 15, height = 6
-)
-```
-
-```{r, message=FALSE, warning=FALSE, echo = FALSE, fig.width=8}
-# Line graph comparing clinical condition counts of BSA and OS data
-fig_pf_descriptive_stats <- df_descriptive_stats %>%
-  mutate(
-    measure = factor(measure,
-      levels = c("pf_with_pfmed", "pf_with_pfcondition", "pf_with_pfmed_and_pfcondition"),
-      labels = c("PF Med", "PF Condition", "PF Med & PF Condition")
-    )
-  ) |>
-  ggplot(aes(
-    x = interval_start,
-    y = ratio,
-    colour = measure,
-    shape = measure,
-  )) +
-  geom_point(size = 2.5) +
-  geom_line(size = 0.5) +
-  labs(
-    x = NULL,
-    y = NULL,
-    shape = "PF consultation linked to:",
-    colour = "PF consultation linked to:"
-  ) +
-  scale_x_date(
-    labels = scales::label_date_short(), breaks = "month"
-  ) +
-  scale_y_continuous(
-    labels = scales::percent,
-  ) +
-  theme(
-    text = element_text(size = 14)
-  ) +
-  geom_vline(
-    xintercept = lubridate::as_date("2024-02-01"),
-    linetype = "dotted",
-    colour = "orange",
-    linewidth = .7
-  ) +
-  scale_colour_brewer(palette = "Dark2")
-
-
-fig_pf_descriptive_stats
-
-ggsave(
-  here("released_output", "results", "figures", "fig_pf_descriptive_stats.png"),
-  fig_pf_descriptive_stats,
-  width = 10, height = 6
-)
-```
-
-```{r, message=FALSE, warning=FALSE, echo = FALSE}
-# Validation of pharmacy first medication counts figure
-# OS data - waiting on released output
-
-df_bsa_medication_validation_sum <- df_bsa_medication_validation %>%
-  group_by(date) %>%
-  summarise(count = sum(count) * 0.4) %>%
-  mutate(
-    source = "nhs_bsa",
-    count_method = "count_40pct"
-  )
-range(df_pfmed$interval_start)
-df_opensafely_pfmed_sum <- df_pfmed %>%
-  rename(date = interval_start) %>%
-  group_by(date) %>%
-  summarise(count = sum(numerator)) %>%
-  mutate(
-    source = "opensafely_tpp",
-    count_method = "opensafely_tpp"
-  )
-
-df_validation_med_counts <- bind_rows(df_opensafely_pfmed_sum, df_bsa_medication_validation_sum) |>
-  filter(date >= "2024-01-01" & date <= "2024-07-01")
-
-fig_validation_med_count <- df_validation_med_counts |>
-  mutate(
-    source = factor(source, levels = c("opensafely_tpp", "nhs_bsa"), labels = c("OpenSAFELY-TPP", "NHS BSA")),
-    count_method = factor(count_method, levels = c("opensafely_tpp", "count_40pct"), labels = c("OpenSAFELY-TPP", "NHS BSA (40%)"))
-  ) |>
-  ggplot(aes(
-    x = date,
-    y = count,
-    colour = count_method,
-    shape = count_method
-  )) +
-  geom_point(size = 2) +
-  facet_wrap(~count_method, scales = "free_y") +
-  geom_line(size = 0.5) +
-  labs(
-    x = NULL,
-    y = "Count",
-    colour = NULL,
-    shape = NULL,
-  ) +
-  scale_x_date(
-    labels = scales::label_date_short(), breaks = "month"
-  ) +
-  theme(
-    text = element_text(size = 14)
-  ) +
-  geom_vline(
-    xintercept = lubridate::as_date("2024-02-01"),
-    linetype = "dotted",
-    colour = "orange",
-    linewidth = .7
-  ) +
-  scale_colour_brewer(palette = "Dark2")
-
-# Another plot visualising the percentage
-fig_validation_med_pct <- df_validation_med_counts %>%
-  filter(count_method %in% c("opensafely_tpp", "count_40pct")) %>%
-  pivot_wider(names_from = c(source, count_method), values_from = count) %>%
-  mutate(source = "Percentage of NHS BSA (40%) in OpenSAFELY") %>%
-  ggplot(
-    aes(
-      x = date,
-      y = opensafely_tpp_opensafely_tpp / nhs_bsa_count_40pct,
-      shape = source,
-      color = source,
-      fill = source,
-      group = source
-    )
-  ) +
-  geom_point(size = 2.5) +
-  geom_line(size = 0.5) +
-  facet_wrap(~source, scales = "free_y") +
-  scale_x_date(
-    labels = scales::label_date_short(), breaks = "month"
-  ) +
-  labs(x = NULL, y = "Percent", colour = NULL, shape = NULL, fill = NULL) +
-  scale_color_viridis_d(
-    option = "plasma",
-    end = 0.9
-  ) +
-  scale_fill_viridis_d(
-    option = "plasma",
-    end = 0.9
-  ) +
-  scale_shape_manual(
-    values = c("Percentage of NHS BSA (40%) in OpenSAFELY" = 15)
-  ) +
-  theme(
-    text = element_text(size = 14)
-  ) +
-  geom_vline(
-    xintercept = lubridate::as_date("2024-02-01"),
-    linetype = "dotted",
-    colour = "orange",
-    linewidth = .7
-  ) +
-  scale_y_continuous(labels = scales::percent)
-
-# Combining the plots with patchwork
-fig_validation_medication_count_pct <- (fig_validation_med_count + fig_validation_med_pct) +
-  plot_annotation(tag_levels = "A") +
-  plot_layout(guides = "collect", widths = c(2, 1)) &
-  theme(
-    legend.position = "bottom",
-    text = element_text(size = 15),
-    strip.background = element_rect(size = 0),
-    strip.text.x = element_text(size = 13, face = "bold")
-  )
-
-fig_validation_medication_count_pct
-
-ggsave(
-  here("released_output", "results", "figures", "fig_validation_medication_count_pct.png"),
-  fig_validation_medication_count_pct,
-  width = 15, height = 6
-)
-```
-```{r, message=FALSE, warning=FALSE, echo = FALSE, fig.width=8}
-# GP vs PF provider graph
-
-df_condition_provider_grouped <- df_condition_provider %>%
-  group_by(measure, interval_start, pf_status) %>%
-  summarise(
-    count = sum(numerator)
-  ) %>%
-  mutate(
-    measure = factor(measure,
-      levels = c(
-        "count_acute_sinusitis_total",
-        "count_infected_insect_bite_total",
-        "count_uncomplicated_urinary_tract_infection_total",
-        "count_acute_otitis_media_total",
-        "count_acute_pharyngitis_total",
-        "count_herpes_zoster_total",
-        "count_impetigo_total"
-      ),
-      labels = c(
-        "Acute Sinusitis",
-        "Infected Insect Bite",
-        "UTI",
-        "Acute Otitis Media",
-        "Acute Pharyngitis",
-        "Herpes Zoster",
-        "Impetigo"
-      )
-    ),
-    pf_status = factor(pf_status,
-      levels = c(TRUE, FALSE),
-      labels = c("Linked to Pharmacy First consultation", "Not linked to Pharmacy First consultation")
-    )
-  )
-
-fig_pf_condition_provider_count <- ggplot(
-  df_condition_provider_grouped,
-  aes(
-    x = interval_start,
-    y = count,
-    colour = pf_status,
-    shape = pf_status
-  )
-) +
-  geom_point(size = 1.5) +
-  geom_line(size = 0.5) +
-  facet_wrap(~measure, scales = "free_y") +
-  labs(
-    x = NULL, y = "Count", color = NULL, shape = NULL
-  ) +
-  theme(
-    plot.title = element_text(hjust = 0.5),
-    legend.position = "bottom",
-    axis.title.x = element_blank()
-  ) +
-  scale_x_date(
-    labels = scales::label_date_short()
-  ) +
-  geom_vline(
-    xintercept = lubridate::as_date("2024-02-01"),
-    linetype = "dotted",
-    colour = "orange",
-    linewidth = .7
-  ) +
-  scale_color_viridis_d(
-    option = "plasma",
-    end = 0.75
-  ) +
-  theme(
-    legend.position = "bottom",
-    text = element_text(size = 14),
-    strip.background = element_rect(size = 0),
-    # strip.text.x = element_text(size = 13, face = "bold")
-  )
-
-fig_pf_condition_provider_count
-
-ggsave(
-  here("released_output", "results", "figures", "fig_pf_condition_provider_count.png"),
-  fig_pf_condition_provider_count,
-  width = 13, height = 8
-)
-```
-
-# References