Common causes of unintentional injuries in Scotland

1 Overview

In this report I am going to analyse data on admissions to hospital and deaths in Scotland from unintentional injuries covering a ten-year period (2013 - 2022). The analysis will be centered on the task of exploring which types of injuries were the most common cause for admissions according to age groups and sex and on finding the death rates in admissions resulted for those main causes. Additionally, we will have a close look into the injury cause with the highest death ratio after admission.

2 Why

To identify the most common and deadly unintentional injuries in Scotland, so that prevention campaigns can be targeted to the highest-risk groups and causes and to provide with the appropriate resources to each health board.

3 Data Packages

We load the software packages we are going to use: tidyverse, janitor, lubridate, kableExtra, formatR, scales.

# Mount packages
library(tidyverse)
library(janitor)
library(lubridate)
library(kableExtra)
library(formatR)
library(scales)

4 About the data

To perform this report, two datasets has been retrieved from The Scottish Health and Social Care Open Data platform for their analysis:

  1. Admissions, Containing 391104 observations with information for 14 variables on emergency hospital admissions as a result of unintentional injuries and assaults.

  2. Deaths, Containing 1825202 observations with information for 14 variables on deaths as a result of an unintentional injury.

#read in .csv files with the data

raw_admissions <- read_csv("https://www.opendata.nhs.scot/dataset/b0135993-3d8a-4f3b-afcf-e01f4d52137c/resource/aee43295-2a13-48f6-bf05-92769ca7c6cf/download/ui_admissions_2023.csv")

raw_deaths <- read_csv("https://www.opendata.nhs.scot/dataset/b0135993-3d8a-4f3b-afcf-e01f4d52137c/resource/89807e07-fc5f-4b5e-a077-e4cf59491139/download/ui_deaths_-2023-.csv")

5 Exploratory Data Analysis

We started exploring the data by reading the Data Dictionary provided in the plataform and by opening both csv files in the R Studio environment, where we checked number of rows and columns, variable names and types of variables. We noticed a mixed naming style on variables and decided to use the clean_names() function right at the begging, creating two new objects: admissions and deaths that will be our main data feed for all the wrangling and analysis. After cleaning the names, we checked again our new created objects to see the new variable names. In this stage we also looked for variables susceptible to be used as the joining argument, and problems that we may encounter (eg: financial_year in admissions dataset is a character type variable where year in deaths dataset is numeric type).

6 cleaning names

#To have all variables with the same naming style
admissions <- raw_admissions %>%  clean_names()
deaths <- raw_deaths %>% clean_names()

7 categories

Possible values for all categorical variables were systematically explored with the distinct(), count() or unique() functions. One example included in the code below, where we can see injury_types names varies in both datasets. We also found duplication in the data due to aggregation of categories in an additional categories named as “All” (or similar). Although we noticed that the category referred as “All” was not always equal to the sums of the other, in any case, we considered we can safely exclude those observations from the analysis by subseting only the categories we want to explore.

# example of how columns were explored to find its unique categories
admissions %>%
  distinct(injury_type)
## # A tibble: 9 × 1
##   injury_type        
##   <chr>              
## 1 All Diagnoses      
## 2 RTA                
## 3 Poisoning          
## 4 Falls              
## 5 Struck by, against 
## 6 Crushing           
## 7 Scalds             
## 8 Accidental Exposure
## 9 Other

8 Plotting Unintended Injuries hospital admissions

For finding and visualizing the most common types injures by age group and sex we have use the data available for the full 10 years period. We have reduce the partitions for the age group aggregation into four groups, making its visualization and understanding more intuitive and less confusing: The wider groups created 0-14, 15-44, 45-65, 75+ can easily be related with common vulnerabilities, behaviors or risk exposures in children, young adults, middle age adults and older adults.

#color palette for the plot (more contrasted that the default ggplot palette).
mycolors <- c( 
  "#6b83d7","#6ec8fa","#e24aae", "#b5c0ce", "#b1ff0b", "#38214e", "#ffc8f1", "#a94ad1"
           )
#subseting the data for the plot, selecting and changing some variable types.
admissions %>%
  separate(financial_year, into = c("year", NA), sep = "/") %>%
  transmute(
    year,
    sex = as.factor(sex),
    age_group = as.factor(age_group),
    injury_type = as.factor(injury_type),
    number_of_admissions
  ) %>%
  # filtering to eliminate aggregated data
  filter(age_group != "All" &
           sex != "All" &  injury_type != "All Diagnoses") %>%
  # Reducing age groups into 4 with a wider age range resulted in a better 
  # grasp of similar risky behaviors/exposure/vulnerabilities.
  mutate(
    age_group = case_when(
      age_group == "0-4 years" ~ "0-14 years",
      age_group == "5-9 years" ~ "0-14 years",
      age_group == "10-14 years" ~ "0-14 years",
      age_group == "15-24 years" ~ "15-44 years",
      age_group == "25-44 years" ~ "15-44 years",
      age_group == "45-64 years" ~ "45-74 years",
      age_group == "65-74 years" ~ "45-74 years",
      age_group == "75plus years" ~ "75+ years")) %>%
  group_by(injury_type, age_group, sex) %>%
  summarise(total = sum(number_of_admissions)) %>%
  ggplot(aes(sex, total, fill = injury_type)) +
  geom_col() + scale_fill_manual(values=mycolors) +
  facet_wrap(~  age_group, nrow = 1) +
  scale_y_continuous(labels = label_number(suffix = " K", scale = 1e-3)) +
  labs(
    title = "Hospital admissions by cause, age and sex",
    subtitle = "Unintentional Injuries in Scotland 2013-2022",
    caption = "Data source: Public Health Scotland",
    y = "Total number of admisions",
    x = " ",
    fill = "Type of injury" ) +
  theme_bw() +
  theme(
    plot.title = element_text(face= 'bold', size = 20),
    plot.subtitle = element_text(color = "steelblue", size = 14, margin = margin(3, 0, 20, 0)))

9 Plotting Death

For finding and visualizing the most common types injures that resulted in death

mycolors2 <- c( 
  "#6b83d7","#6ec8fa","#e24aae","#38214e", "#b5c0ce", "#b1ff0b", "#ffc8f1", "#a94ad1"
           )
#subseting the data for the plot, selecting and changing some variable types.
deaths %>%
  
  transmute(
    sex = as.factor(sex),
    age_group = as.factor(age_group),
    injury_type = as.factor(injury_type),
    numberof_deaths
  ) %>%
  # filtering to eliminate aggregated data
  filter(age_group != "All" &
           sex != "All" &  injury_type != "All") %>%
  # Reducing age groups into 4 with a wider age range resulted in a better 
  # grasp of similar risky behaviors/exposure/vulnerabilities.
  mutate(
    age_group = case_when(
      age_group == "0-4 years" ~ "0-14 years",
      age_group == "5-9 years" ~ "0-14 years",
      age_group == "10-14 years" ~ "0-14 years",
      age_group == "15-24 years" ~ "15-44 years",
      age_group == "25-44 years" ~ "15-44 years",
      age_group == "45-64 years" ~ "45-74 years",
      age_group == "65-74 years" ~ "45-74 years",
      age_group == "75plus years" ~ "75+ years")) %>%
  group_by(injury_type, age_group, sex) %>%
  summarise(total = sum(numberof_deaths)) %>%
  ggplot(aes(sex, total, fill = injury_type)) +
  geom_col() + scale_fill_manual(values=mycolors2) +
  facet_wrap(~  age_group, nrow = 1) +
  scale_y_continuous(labels = label_number(suffix = " K", scale = 1e-3)) +
  labs(
    title = "Deaths by cause, age and sex",
    subtitle = "Unintentional Injuries in Scotland 2013-2022",
    caption = "Data source: Public Health Scotland",
    y = "Total number of deaths",
    x = " ",
    fill = "Type of injury" ) +
  theme_bw() +
  theme(
    plot.title = element_text(face= 'bold', size = 20),
    plot.subtitle = element_text(color = "steelblue", size = 14, margin = margin(3, 0, 20, 0)))

10 Death ratio after admission for each cause

In this piece of code we first created the table “total_admissions” with all admissions grouped by injury type. We had to change some names of the categories in order to match those on the second dataset. Then, we created a second table “total_deaths” with number of deaths by injury type. Both newly created tables were joined into a new one called admissions_deaths. Death rates in admissions were calculated for each injury cause with the mutate() function. The resulting ranking table is shown below.

#create a table with number of admission per type of injury
total_admissions <- admissions %>%
  separate(financial_year, into = c("year", NA), sep = "/") %>% 
  select(sex, age_group, injury_type, number_of_admissions) %>%
  #categories in both datasets need to have same names
  mutate(injury_type = str_replace(injury_type, "RTA", "Land transport accidents"), 
         injury_type = str_replace(injury_type, "Struck by, against", "Struck by,against"),
         injury_type = str_replace(injury_type, "Accidental Exposure", "Accidental exposure")
         ) %>%
  #filtering to eliminate duplicated data
  filter(age_group != "All" & sex != "All" & injury_type != "All Diagnoses") %>%
  group_by(injury_type) %>%
  summarise(total_admissions = sum(number_of_admissions, na.rm=TRUE))


#create a table with number of death per type of injury
total_deaths <- deaths %>%
  select(sex, age_group, injury_type, numberof_deaths) %>%
  #filtering to eliminate duplicated data
  filter(age_group != "All" &  sex != "All" &  injury_type != "All") %>%
  group_by(injury_type) %>%
  summarise(total_deaths = sum(numberof_deaths, na.rm=TRUE))


#join both tables and calculate rate of death in admissions 
#per type of injury 
admissions_deaths <- total_admissions %>%
  left_join(total_deaths, by = c("injury_type")) %>%
  mutate(death_ratio = total_deaths/total_admissions) %>%
  #ordering the values in descending order by death ratio
  arrange(desc(death_ratio))

#to display the table limiting decimals to 3 
kable(admissions_deaths, 
      caption = "Death rates in Scotland for injury types. 2013-2022", 
      digits = 3) %>%
  kable_styling(latex_options = c("HOLD_position"), font_size = 12)
Death rates in Scotland for injury types. 2013-2022
injury_type total_admissions total_deaths death_ratio
Poisoning 132096 38808 0.294
Land transport accidents 116100 7112 0.061
Accidental exposure 135080 4048 0.030
Falls 1440776 36548 0.025
Other 314676 6436 0.020
Scalds 18992 56 0.003
Crushing 44268 56 0.001
Struck by,against 90436 112 0.001

11 Visualizing death rate by injury type

admissions_deaths %>%
  ggplot(aes(x = reorder(injury_type, death_ratio), y = death_ratio)) +
  geom_col(color="red", fill='pink') + 
  coord_flip() +
  labs(
    title = "Death ratio by Injury type",
    subtitle = "Scotland 2013-2022",
    caption = "Data source: Public Health Scotland",
    y = "Deaths/Admissions ratio",
    x = "",
    fill = "total_deaths" ) +
  geom_text(aes(label = round(death_ratio, 3)), hjust = -0.1, size = 3, color='red')

12 Findings

  1. Having a fall was the most common reason for hospital admission for all age groups and sexes.

  2. For under 75 years old, there have been more hospital admissions for males than for females, being this difference greater for the age group 15-44 years. This may be due to males more prone to engage themselves in risky activities and behaviors than their females counterparts.

  3. For the 75+ years group, there have been more admissions of females than males. Probably due to the higher proportion of female versus males in the total population within this age bracket.

  4. Males between 15 to 44 years have the higher total number of admissions for poisoning, accidental exposure, other injuries, traffic accident and Struct by.

  5. Among all the unintentional injury causes registered in admission, poisoning have the higher death ratio (0.294), calculated for the 10 years data. Followed by transport accidents (death rate = 0.061) and falls (death rate = 0.030).

13 Futher exploration based on results

As poisoning was by far the type of injury with the highest death rate in admissions among all injury types, further exploration focused in this cause were performed. We decided to explore how admissions, deaths and death rates differ across all the NHS health boards.

14 Poisoning death rates on each NHS Scottish Health Board

We have calculated total admissions for poisoning for each health board form the data in the admissions dataset, and total deaths for poisoning for each health board from the data in the deaths dataset. Then, joined both together.

#create a table with total poison adissions by health board
poison_admissions_hb <- admissions %>%
  separate(financial_year, into = c("year", NA), sep = "/") %>%
  select ( year, hbr, sex, age_group, injury_type, number_of_admissions) %>%
  #filtering to eliminate aggregated data and subseting poisoning
  filter(age_group != "All" & sex != "All" & hbr != "S92000003" & 
           injury_type == "Poisoning") %>%
  group_by(hbr)%>%
  summarise(total_poison_admissions = sum(number_of_admissions, na.rm=TRUE))

#create a table with number of deaths due to poisoning in each Health Board
poison_deaths_hb <- deaths %>%
  select (year, hbr, sex, age_group, injury_type, numberof_deaths) %>%
  #filtering to eliminate aggregated data and subseting poisoning
  filter(age_group != "All" & sex != "All" & hbr != "S92000003" & 
           injury_type == "Poisoning") %>%
  group_by(hbr)%>%
  summarise(total_poison_deaths = sum(numberof_deaths, na.rm=TRUE))

#create a table by joining both tables 
#and add column with rate of death calculation per Heath Board
poison_deaths_rates_hb <- poison_admissions_hb %>%
  left_join(poison_deaths_hb, by = c("hbr")) %>%
  mutate(death_ratio = total_poison_deaths/total_poison_admissions)

#to display the table limiting decimals to 3 
kable(poison_deaths_rates_hb, 
      caption = "Poisoning data by HB 2013-2022", 
      digits = 3) %>%
  kable_styling(latex_options = c("HOLD_position"), font_size = 12)
Poisoning data by HB 2013-2022
hbr total_poison_admissions total_poison_deaths death_ratio
S08000015 4890 1566 0.320
S08000016 866 242 0.279
S08000017 1444 468 0.324
S08000019 3286 1030 0.313
S08000020 5460 1548 0.284
S08000022 2978 690 0.232
S08000024 8770 2518 0.287
S08000025 92 18 0.196
S08000026 196 36 0.184
S08000028 592 36 0.061
S08000029 5918 1100 0.186
S08000030 3478 1518 0.436
S08000031 19844 6122 0.309
S08000032 8234 2512 0.305

15 Adding the NHS Health Boards names to our table

Not necessary for our map, but if you are no familiar with the code and want to add the names of the Health boards, you can do it by joining our data with the following one also from Public Health Scotland in which all 14 Health Boards are listed with their corresponding name and the country code for Scotland.

HBR <- read_csv('https://www.opendata.nhs.scot/dataset/9f942fdb-e59e-44f5-b534-d6e17229cc7b/resource/652ff726-e676-4a20-abda-435b98dd7bdc/download/hb14_hb19.csv')

hbr_names_poison <- HBR %>%
inner_join(poison_deaths_rates_hb, by = c("HB" = "hbr")) %>%
  select(HB, HBName, total_poison_admissions, total_poison_deaths, death_ratio) %>%
  arrange(desc(death_ratio))

kable(hbr_names_poison, 
      caption = "Poisoning data for each health board (2013-2022)", 
      digits = 3) %>%
  kable_styling(latex_options = c("HOLD_position"), font_size = 12)
Poisoning data for each health board (2013-2022)
HB HBName total_poison_admissions total_poison_deaths death_ratio
S08000030 NHS Tayside 3478 1518 0.436
S08000017 NHS Dumfries and Galloway 1444 468 0.324
S08000015 NHS Ayrshire and Arran 4890 1566 0.320
S08000019 NHS Forth Valley 3286 1030 0.313
S08000031 NHS Greater Glasgow and Clyde 19844 6122 0.309
S08000032 NHS Lanarkshire 8234 2512 0.305
S08000024 NHS Lothian 8770 2518 0.287
S08000020 NHS Grampian 5460 1548 0.284
S08000016 NHS Borders 866 242 0.279
S08000022 NHS Highland 2978 690 0.232
S08000025 NHS Orkney 92 18 0.196
S08000029 NHS Fife 5918 1100 0.186
S08000026 NHS Shetland 196 36 0.184
S08000028 NHS Western Isles 592 36 0.061

16 Spatial data

For this additional analysis we have used geographical spatial data for the Scottish Health Boards, a ESRI Shape file spatial data defining the boundaries of NHS Health Boards in Scotland, Available open source from the Spatial Data Metadata Portal, Scotland’s catalogue of spatial data.

17 Packages for dealing with spatial data

We load additional libraries for reading a dealing with this type of data

# Load packages

library(sp)
library(sf)
library(gridExtra)
library(latticeExtra)

18 Joining the spatial data with our health data

In this chunk of code we read the .shp file containing the vector data with the health board boundaries and join our newly created table poison_deaths_rates_hb to it.

#read in .shp file 
# you need to have all spatial data files in your working directory write your own file path
path <- "C:/Users/Casa/Desktop/IRM/Injuries/SG_NHS_HealthBoards_2019.shp"
scotland_hb <- st_read(path)
## Reading layer `SG_NHS_HealthBoards_2019' from data source 
##   `C:\Users\Casa\Desktop\IRM\Injuries\SG_NHS_HealthBoards_2019.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 14 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 5512.998 ymin: 530250.8 xmax: 470332 ymax: 1220302
## Projected CRS: OSGB36 / British National Grid
#join the spatial data with the poisoning data
join_data <- scotland_hb %>%
inner_join(poison_deaths_rates_hb, by = c("HBCode" = "hbr"))

19 Building and plotting the chloropeth maps

We create the 3 map plots

map3 <- ggplot(join_data, aes(fill = death_ratio)) + 
  geom_sf(size = 0.3, color = "#1f1b39") + 
  scale_fill_viridis_c(option = "viridis", direction = -1) +
    labs(
    title = "Death ratio for poisoning",
    subtitle = "period 2013-2022") +
  coord_sf() +
  theme_void()

map2 <- ggplot(join_data, aes(fill = total_poison_deaths)) + 
  geom_sf(size = 0.3, color = "#1f1b39") + 
  scale_fill_viridis_c(option = "viridis", direction = -1) +
    labs(
    title = "Deaths for poisoning ",
    subtitle = "on 2013-2022") +
  coord_sf() +
  theme_void()

map1 <- ggplot(join_data, aes(fill = total_poison_admissions)) + 
  geom_sf(size = 0.3, color = "#1f1b39") + 
  scale_fill_viridis_c(option = "viridis", direction = -1) +
    labs(
    title = "Admissions for poisoning",
    subtitle = "on 2013-2022") +
  coord_sf() +
  theme_void()
map1

map2

map3

20 Plotting the maps all together in one figure

We can also print the 3 maps plots together as one figure for comparing the visualization. For that we can use the cowplot package.

#library containing the function plot_grid()
library(cowplot)

#plotting the 3 maps together
plot_grid(map1, map2, map3, align = "h")