Introduction

Is there an association between antibiotic prescription rate and cold weather? How does this vary by health board?

Antibiotics are among the most widely prescribed classes of medications, with 31.2% of the population being prescribed annually in the Scotland [1]. Consequently, understanding the factors that influence prescription is paramount to healthcare. One benefit of improving our understanding of prescription trends is the ability to combat antimicrobial resistance, which is consistently rising [2]. In this report, we are looking at the antibiotic dosage rate, over 2 June-December periods, and observing if there is any correlation between dosage rate in health boards, and the minimum average temperature.

To investigate this, we will be using prescription data from NHS Scotland, population data from Public Health Scotland, and temperature data from the MET office. We will start by using static graphs to show preliminary findings. Secondly, an interactive app to show temporal/spatial relationships in a dynamically.

Our Packages

Packages = c(
  "tidyverse", "janitor", "sf", "terra", "exactextractr", "shiny", "here", "reactable", "rmdformats", "ggiraph", "scales", "reactablefmtr") # creates list of required packages
Missing = Packages[!Packages %in% installed.packages()[, "Package"]]
if (length(Missing) > 0) install.packages(Missing) # checks to see if they are already installed, in not it will install missing packages

invisible(lapply(Packages, library, character.only = TRUE)) #opens library of all packages

Reading in Our Data

Temperature Data

First, let’s load our temperature data. We can’t directly download our data from the CEDA website due to the login required to download (and the API, which bypasses this, is not yet functional).

To access the data, you can create a CEDA account and download the data and put it in your temperature data folder: HadUK-grid 1km grid, v1.3.1.ceda We’ll be using data from the HadUK gridded data set, which provides a comprehensive overview of temperatures trends with high resolution. From this data set, we can construct a highly precise record of average minimum temperatures.

Years = c(2023, 2024) # we'll save this value, so we can rename our layers later easily using rep()

if (!dir.exists(here("Temperature_Data"))) {dir.create(here("Temperature_Data"))}

Temp_Data_Files = list.files(path = here("Temperature_Data"), pattern = "\\.nc$", full.names = TRUE)
All_Temp_Data_Rast = lapply(Temp_Data_Files, rast)
All_Temp_Data = c(All_Temp_Data_Rast[[1]], All_Temp_Data_Rast[[2]])

names(All_Temp_Data) = paste0(rep(Years, each = 12), "_", month.name)

Antibiotics Data

We can download our data directly from the Website at:
Public Health Scotland, Prescriptions in the Community We’ll use paid_quantity as our antibiotic metric. Unlike number_of_paid_items, which counts both single and multiple orders as 1. paid_quantity consistently reflects the actual medication volumes, making it more appropriate.

# First you'll need a folder called "Antibiotics_Data", the script below will check and create one if its not present.
if(!dir.exists(here("Antibiotics_Data"))) {dir.create(here("Antibiotics_Data"))}
# Download the data from the link above. we'll be working with June-December in 2023 and 2024, in a .CSV format, once downloaded, put them into your antibiotics folder.
Antibiotics_Data = list.files(path = here("Antibiotics_Data"), pattern = "\\.csv$", full.names = TRUE)

Antibiotics_Data_List = map(Antibiotics_Data, ~ read_csv(.x) %>%
  clean_names() %>%
  select(hbt, bnf_item_code, paid_quantity, paid_date_month) %>%
  filter(startsWith(
            bnf_item_code, "05010"),
            hbt != "SB0806"))
Antibiotics_Data_Combo = bind_rows(Antibiotics_Data_List)
# First we will collate our data into a list, so we can process it all with one command, and add more data if required easily, here's the first couple, expand the code below to see the full thing.

Population Data

In order to accurately compare between more populated and less populated regions, we’re going to use the metric ‘Antibiotic doses per 1000’. We’ll need to divide our Antibiotic data by population. We’ll be using data from Public Health Scotland again to provide population statistics. If you’d like to manually download the data, it can be found here:
Public Health Scotland, Health Board (2019) Population Estimates

Pop = read_csv("https://www.opendata.nhs.scot/dataset/7f010430-6ce1-4813-b25c-f7f335bdc4dc/resource/27a72cc8-d6d8-430c-8b4f-3109a9ceadb1/download/hb2019_pop_est_03102025.csv") %>%
   clean_names() %>%
   filter(year %in% c("2023", "2024") &
     sex == "All" &
     hb != "S92000003") %>%
    select(hb, all_ages, year) # R can get confused to which select() you want to use, so adding dplyr:: before select() helps R pick the correct one

Map data

We can download our map data directly from the Scottish Government website: spatialdata.gov.scot

 Scotland_Map = st_read("/vsizip/vsicurl/https://maps.gov.scot/ATOM/shapefiles/SG_NHS_HealthBoards_2019.zip/SG_NHS_HealthBoards_2019.shp", quiet = TRUE) %>% # /vsizip/vsicurl allows our data to open a zip file, and download it over the internet
   st_simplify(dTolerance = 1000) %>% # st_simplify makes our map data less data-intensive, and helps it run quicker
   st_transform(crs(All_Temp_Data)) # st_transform makes sure that our map data will work with our other data

Data Wrangling and Joining

Temperature Data

Temp_df = exact_extract(All_Temp_Data, Scotland_Map, "mean",
  full_colnames = FALSE,
  append_cols = "HBName") %>%
  select(-matches("(January|February|March|April|May)"))
# exact_extract adds "mean." to our column names. You can't prevent this. We have to remove it, so we can reliably process our data downstream.
names(Temp_df) = gsub("mean.", "", names(Temp_df))

Temp_df_Long = Temp_df %>%
  pivot_longer(
    cols = -c(HBName), # we select the columns we want to turn into rows (all except HBName which we dont want to change)
    names_to = "Year_Month", # new column name for month
    values_to = "Min_Temperature") %>% # new column name for values
  separate("Year_Month", into = c("Year", "Month")) %>%
  mutate(Year = as.numeric(Year)) # we encode our year as a numeric, to ensure it can join other data later.

Antibiotics Data

Date_Extracted_Antibiotics = Antibiotics_Data_Combo %>%
  mutate(
    Year = as.numeric(substr(paid_date_month, 1, 4)), # first 4 chars = year
    Month = month.name[as.numeric(substr(paid_date_month, 5, 6))]) # turn month into long form, for easy reading.
Sum_Antibiotics_Data = Date_Extracted_Antibiotics %>% 
  group_by(hbt,Year,Month) %>% 
  summarise(total_quanity = sum(paid_quantity, na.rm = TRUE), .groups = "keep")%>% 
  full_join(Pop, by = join_by(hbt == hb, Year == year))
# full  join ensures we keep all the data in both data sets, even when there isn't a match. This means if there's incomplete sets, we can identify them by using anyNA(Tally_Antibiotics_Data) 

Antibiotics_Dose_PC = Sum_Antibiotics_Data %>% 
  mutate(Dose_per_1000 = total_quanity/(all_ages/1000)) %>% 
  subset(select = -c(all_ages, total_quanity))
 Antibiotics_Temp_Map = Scotland_Map %>% 
  full_join(Antibiotics_Dose_PC, join_by(HBCode == hbt))%>% 
             right_join(Temp_df_Long) %>% 
  relocate(c("Shape_Leng", "Shape_Area"), .before = geometry) # we will move our map data out of the way, so the data we want to review are together for improved readability
    
  Antibiotics_Temp_Table = Antibiotics_Temp_Map %>% 
    st_drop_geometry() %>% # we have to remove geometry data, reactable doesn't display data if it has geometry attached.
  select(-HBCode, -Shape_Leng, -Shape_Area) # removing excess/redundant information helps improve readability

Tables of our data

 dropdown_filter = function(values, name) { # drop down provides filtering options for our data, once we put it into a table
   tags$select(
     onchange = sprintf(
       "Reactable.setFilter('antibiotics-table', '%s', event.target.value || undefined)", name),
     tags$option(value = "", "All"),
     lapply(sort(unique(values)), tags$option),
     style = "width: 100%; height: 30px;")}

reactable( # we'll use reactable() to make an interactive table
  Antibiotics_Temp_Table,
  filterable = TRUE,
  pagination = FALSE, # This lets us have a scroll bar rather than pages
  height = "400px",
  defaultColDef = colDef(filterInput = dropdown_filter),
    columns = list(
    Dose_per_1000 = colDef(
                    filterable = FALSE,
                    name = "Doses per 1000"),
   Min_Temperature = colDef(
       filterable = FALSE,
       name = " Average Minimum Temperature (ºC)"),
   HBName = colDef(
     name = "Health Board"),
     elementId = "antibiotics-table")) %>% 
  add_subtitle("Table 1: Antibiotic dosage rate vs average minimum temperature (ºC) in 2023 and 2024, by date and healthboard",
               align = "center")

Table 1: Antibiotic dosage rate vs average minimum temperature (ºC) in 2023 and 2024, by date and healthboard

We can see that August is consistently the least cold month, accounting for 8/10 of the highest minimum temperatures. 2023 also saw higher minimum temperatures, accounting for 8/10 of the highest minimum temperatures. 7/10 of the 10 coldest months occurred in their Grampians, Tayside, or the Highlands. We see that Orkney has the 4 lowest monthly dosage rates. The western central belt saw the highest dosage rates, with Lanarkshire and Ayrshire and Arran consistently seeing high dosage rates, making up 6/10 of the highest dosage rates across both years.

Preliminary Graphs, and Looking for Associations

ggplot(Antibiotics_Temp_Table, aes(x = Min_Temperature, y = Dose_per_1000)) +
  geom_point(size = 1,
             color = "blue") +
  facet_wrap(~ Year) +
  geom_smooth(method = "lm",
              color = "red") +  # geom_smoth() adds a trend line, method = "lm" make it straight, 
  labs(x = "Average Minimum Temperature (ºC)", y = " Antibiotic Doses per 1000", #data and graph labels
       title = "Scatter plot of Antibiotic Doses per 1000 vs Average Minimum Temperature") +
  theme_minimal() # formatting to it fits our page nicely
Fig. 1 Scatter plot of antibiotic use vs average minimum temperature (ºC) in 2023 and 2024. Red = regression line; Shaded = Standard error

Fig. 1 Scatter plot of antibiotic use vs average minimum temperature (ºC) in 2023 and 2024. Red = regression line; Shaded = Standard error

From Fig 1, we can see a correlation between average minimum temperature and the antibiotic dosage rate, as indicated by the downward trend of the regression lines. However, we cannot determine how different health boards are affected individually, and without the trend line, broader patterns would be harder to detect. Plotting this data as a choropleth would reveal spatial relationships and allow us to observe temporal trends among and between health boards.

Interactive choropleths

We’ll use Shiny, and shiny.IO to make and host an app which will let us browse the full range of choropleths and identify any spatial-temporal trends in the data. You can find the introduction to Shiny, and Shiny.IO below:
Shiny
Shiny.IO

From the choropleths, we can see an inverse relationship between antibiotic prescription and minimum temperature - as it gets less cold, less antibiotics are prescribed. This cooperates with the correlation we observed in Fig. 1. We see that Tayside is consistently the coldest region (although not consistently the most prescribed). The western central belt appears to regularly see the highest antibiotic dosage rates, despite being one of the least cold regions.

Whilst we have observed a correlation. One important aspect to consider is the reported over prescription by doctors in the winter months [3]. This could artificially inflate antibiotic numbers, especially in areas with higher elderly populations. our data corroborates with this idea, with 8 out of the 10 highest dosage rates being from health boards with larger elderly demographics[4].

if(!dir.exists(here("App_Data"))) {dir.create(here("App_Data"))}
if(!file.exists(here("App_Data/app.R"))) {dir.create(here("App_Data/app.R"))}
saveRDS(Antibiotics_Temp_Map, here("App_Data/Antibiotics_Temp_Map.rds"))
#Although downloading the full tidyverse is convenient, loading only the needed packages is better. Shiny.IO has limited memory.
library(shiny) #this chunk should be loaded into your new app.R file
library(scales)
library(sf)
library(dplyr)
library(ggplot2)
library(ggiraph)
# If deploying/uploading, library(rsconnect) / install.packages(rsconnect)
Antibiotics_Temp_Map = readRDS(here("App_Data/Antibiotics_Temp_Map.rds"))
#remove here(App_Data/) when uploading to Shiny.IO
UI = fluidPage( # fluid page tells R that this is an interactive element
  style = "padding: 0; margin: 0;",
  titlePanel("Fig 2: Antibiotic Doses per 1000 / Average Minimum Temperature (ºC), by Health board"),
  sidebarLayout( #the sidebar will be our filter drop down area
    sidebarPanel( width = 3, style = "padding: 0; margin: 0;",
      selectInput("Year", "Select Year",  # 1st drop down box,for year
                  choices = c(2023, 2024)), #options for drop down box
      selectInput("Month", "Select Month", # 2nd drop down box, for month
                  choices = month.name[6:12])) ,# options for drop down box
 mainPanel(
      splitLayout(
        cellWidths = c("50%", "50%"), # splits the row 50/50
        cellArgs = list(style = "padding: 0;"), #removes padding
        girafeOutput("Min_Temp_Map", height = "100%", width = "100%"),
        girafeOutput("Antibiotic_Map", height = "100%", width = "100%") )))) 
# tells R which plots we want to display (we define what they are in the server section)

Server = function(input, output, session) {
  
Filter = reactive({ # reactive tells R that it needs to update everytime a new selection is made
    Antibiotics_Temp_Map %>%
      filter(Year == input$Year, #tells R that it has to filter based on the input of the drop downs.
             Month == input$Month) })
  
output$Antibiotic_Map = renderGirafe({
    Antibiotics_Data = Filter() #Filtered() acts as a placeholder for our filtered data
    ggplot_Antibiotics = ggplot(Antibiotics_Data) +
  geom_sf_interactive(
        aes(fill = Dose_per_1000,
            tooltip = paste0("HB :", HBName, "<br>",
                             "Dose per 1000: ", round(Dose_per_1000, 2)),
            data_id = HBName),
            color = "white") + # we'll make the outline of the map white 
  scale_fill_gradient(
        low = "#F5E3FA",
        high = "#A50026",
        limits = c(1200, 3200), # we set a uniform limit to all data, we we can compare month by month easily
        labels = label_comma()) +
  labs(
        title = paste("Antibiotic Doses per 1000 —",input$Month, input$Year),
        subtitle = "By Scottish Health Boards",
        fill = "Dose Per 1000") +
  guides(fill = guide_colorbar(
        barwidth = 10,   # width of the legend bar
        barheight = 1)) +
  theme_minimal() +
  theme(legend.position = "bottom") # move legend to the bottom to prevent ambiguity 
    
  girafe(ggobj = ggplot_Antibiotics,
           width_svg = 5,   # width
           height_svg = 7) })

output$Min_Temp_Map = renderGirafe({
    Min_Temp = Filter()
    ggplot_Temp = ggplot(data = Min_Temp) +
  geom_sf_interactive(aes(fill = Min_Temperature, 
                            tooltip = paste0("HB :", HBName,"<br>", 
                                        "Min temp :", round(Min_Temperature, 2)), 
                                         data_id = HBName),
                            color = "white") +
  scale_fill_gradient2(
      low = "#313695",
      mid = "#F5E3FA",
      high = "#A50026",
      midpoint = 5.5,
      limits = c(-0.75, 11.75),
      labels = label_number()) +
  labs(
      title = paste("Minimum Temperature —", input$Month, input$Year),
      subtitle = "By Scottish Health Boards",
      fill = "Minimum Temperature (ºC)") +
  guides(fill = guide_colorbar(
      barwidth = 10,
      barheight = 1)) +
  theme_minimal() +
  theme(legend.position = "bottom")
  girafe(ggobj = ggplot_Temp,
         width_svg = 5,   #width
         height_svg = 7) })}
shinyApp(UI, Server) # tells R to run the app

#To upload your app to Shiny.IO, use rsconnect with the token provided when creating your account. You don’t need to include it in the app itself.
# Copy and paste this, outside of the chunk on its own:
# <iframe 
# src="[link given when deploying app" 
# data-external="1"
# style="border: none"
# width="100%" 
# height="700px"
# scrolling="no"
# seamless = TRUE>
# </iframe>

Next steps

This gives us a meaningful insight into the trends that could affect antibiotic prescription. With these findings, healthcare providers could use weather forecasts to help predict pathogenesis in the community and prepare for seasonal demand efficiently. By sourcing antibiotics more efficiently, healthcare providers can save money by preventing over-ordering and disposing of unused expired medications.

This is not the only issue faced by seasonal prescriptions, as mentioned previously, over-prescription can promote antibiotic resistance. Future studies should also look at the bacterial:non-bacterial infection rate and compare it with the antibiotic dosage rate. In regions where antibiotics are incorrectly prescribed, we would be able to spot that bacterial:non-bacterial infection rates are low with a disproportionately high antibiotic dosage rate and could provide significant public benefits by helping curb antibiotic resistance development.

References

[1] “Scottish One Health Antimicrobial Use and Antimicrobial Resistance in 2023” - Antimicrobial Resistance and Healthcare Associated Infection Scotland, 19 November 2024, Link
[2] “Policy paper - Confronting antimicrobial resistance 2024 to 2029” - Department of Health and Social Care, 8 May 2024 Link
[3]‘Analysis of seasonal variation of antibiotic prescribing for respiratory tract diagnoses in primary care practices’, Serletti, L. et al. 5 Septemeber 2023, Antimicrobial Stewardship & Healthcare Epidemiology, 3(1), p. e147. Link.
[4] “Population: estimates, by NHS Board”, The Scottish Public Health Observatory, 13 April 2023 Link