drug_deaths_analysis.qmd

---
title: "Deaths associated with substance use"
subtitle: "Premature mortality and life-years lost due to drug and alcohol use"
author: "OHID"
date: today
number-sections: true
format:
  DHSC-pptx
execute:
  echo: false
  warning: false
  message: false

---

# Number of deaths associated with drug use

---

- The ONS classifies death related to drug poisoning according to ICD-10 codes.
- Certain ICD-10 codes classify a death as a "drug misuse death".[^1] 
- Each of these requires a specific substance (e.g. heroin) or substance category (e.g. opioids) to be indicated either in the ICD-10 code or on the death certificate.


- There are deaths each year where the ONS holds no information on the substance(s) involved (see table 1).

| Year of death registration | All drug poisonings | Number of deaths without substance information | Percentage without substance information |
|----------------------------|---------------------|-----------------------------------------------|------------------------------------------|
| 2023                       | 5448                | 1245                                          | 22.9                                     |
| 2022                       | 4907                | 1239                                          | 25.2                                     |
| 2021                       | 4859                | 1219                                          | 25.1                                     |
| 2020                       | 4561                | 1050                                          | 23.0                                     |


Table 1: Number of drug poisonings deaths without substance information.

------------------------------------------------------------------------

-   Some of these will be classified as related to drug misuse where an ICD-10 code indicates mental and behavioural disorders due to drug use (excluding alcohol and tobacco) without a specific substance (e.g. F19 "multiple drug use and use of other psychoactive substances").

-   Others, coded as accidental/intentional self-poisonings or self-poisonings of unknown intent, will not be classified as related to drug misuse unless a controlled drug under the Misuse of Drugs Act 1971 was mentioned on the death record.

-   The data linkage between ONS mortality and NDTMS allows some of those deaths to be identified indirectly as related to drug misuse where the person that died had had contact with the drug treatment.    

------------------------------------------------------------------------

```{r}
library(flextable)
source("R/drug_deaths_functions.R")
source("R/plotting_functions.R")
source("R/themes.R")
source("R/dhsc_colour_palette.R")

drug_poisoning_deaths_file <-
  "data/raw/ndtms_mortality_data.parquet"
deaths_in_treatment_file <-
  "data/raw/tx_deaths_la_2122_2223.parquet"

national_data <- 
  combine_national_data(
    poisoning_data = process_poisoning_data(file_path = drug_poisoning_deaths_file,
                                            date_of = "occurence",
                                            years = 2022),
    treatment_deaths_data = process_deaths_in_treatment(
      file_path = deaths_in_treatment_file,
      years = 2022,
      exclude_poisoning = TRUE,
      by_treatment_status = TRUE,
      by_death_cause = FALSE,
      exclude_alcohol_specific_deaths = TRUE
    )
  )

national_data <- 
  relabel_national_data(national_data = national_data)

p1 <- plot_national_data(plot_data = national_data)

p1 <- add_plot_annotations(plot = p1, data = national_data)

png(filename = "plots/plot_1.png", height = 18,  width = 26, units = "cm", res = 200)
p1
invisible(dev.off())

```

![](plots/plot_1.png){fig-alt="Bar chart showing the composition of deaths associated with drug use"}

------------------------------------------------------------------------

```{r, out.width=60, out.height=30, dpi = 300}
drug_misuse_age_data <- function(){
 df1 <- 
process_poisoning_data(file_path = drug_poisoning_deaths_file, by = "age", years = 2022) |> 
  mutate(death_cause = "Drug poisoning") |> 
  select(-additional_poisoning_deaths) |> 
  select(-dod_year) |> 
  rename("age" = ageinyrs)

 df2 <- 
process_deaths_in_treatment(file_path = deaths_in_treatment_file, by = "age", by_treatment_status = TRUE, by_death_cause = TRUE, years = 2022) |> 
  filter(treatment_status != "Died one or more years following discharge") |> 
  group_by(age, death_cause) |> 
  summarise(count = sum(count))
 
 bind_rows(df1, df2) |> 
   mutate(death_category = "Deaths associated with drug use")
 
}

ons_leading_mortality_causes <- 
read_csv("data/processed/ons_leading_mortality_causes.csv") |> 
  filter(age_group != "All ages")


cut_age_groups <- 
  function(x) {
    cut(
      x,
      breaks = c(19, 34, 49, 64, 79, Inf),
      labels = c(
        "20 to 34 years",
        "35 to 49 years",
        "50 to 64 years",
        "65 to 79 years",
        "80 years and over"
      ),
      right = TRUE
    )
  }

deaths_related_to_drug_misuse <- 
drug_misuse_age_data() |>
  group_by(age, death_category) |>
  summarise(count = sum(count)) |>
  filter(age > 19) |> 
  rowwise() |>
  mutate(age_group = cut_age_groups(age)) |> 
  group_by(age_group, death_category) |> 
  summarise(count = sum(count))


ons_leading_mortality_causes <- 
  read_csv("data/processed/ons_leading_mortality_causes.csv") |> 
  filter(age_group != "All ages") |> 
  select(age_group, leading_cause, deaths) |> 
  filter(age_group %in% c(
    "20 to 34 years",
    "35 to 49 years",
    "50 to 64 years",
    "65 to 79 years",
    "80 years and over"
  )) |> 
  rename("death_category" = leading_cause) |> 
  rename("count" = deaths)

plot <- 
bind_rows(deaths_related_to_drug_misuse, ons_leading_mortality_causes) |> 
  filter(age_group != "80 years and over") |> 
  arrange(age_group, count) |> 
  mutate(death_category = as_factor(death_category)) |> 
  ggplot(aes(x = death_category, y = count, fill = death_category)) + 
  geom_col(colour = "black") + 
  facet_wrap(~age_group, scales = "free") + 
  coord_flip() + 
  scale_y_continuous(labels = scales::comma) + 
  scale_fill_manual(values = c("Deaths associated with drug use" = "red")) + 
  my_theme +
  theme(legend.position = "none", plot.title.position = "plot") + 
  labs(
    title = "Deaths associated with drug misuse",
    subtitle = "Compared with leading causes of death by age group",
    y = "Number of deaths (2022)",
    x = "Cause of death"
    )

plot 

```

------------------------------------------------------------------------
## Causes of non-poisoning deaths by treatment status
```{r}
ft <- 
 create_cause_of_death_table(
  process_deaths_in_treatment(
    file_path = deaths_in_treatment_file,
    years = 2022,
    exclude_poisoning = TRUE,
    by_treatment_status = TRUE,
    by_death_cause = TRUE
  )
) 

ft
 
```

------------------------------------------------------------------------

```{r, fig.height=20, fig.width=40, fig.dpi=200}

library(tidyverse)
library(afcolours)
library(openxlsx)
library(curl)

# Getting all cause deaths data, initially to see which age groups to use for
# drug/alcohol deaths; then to compare.

df <-
  read_csv("data/processed/drug_deaths_national.csv")


if (!file.exists("data/raw/ons_deaths_data_2023.xlsx")) {
  url <-
    "https://www.ons.gov.uk/file?uri=/peoplepopulationandcommunity/birthsdeathsandmarriages/deaths/datasets/deathsregisteredinenglandandwalesseriesdrreferencetables/2023/annualdeathregistrations2023.xlsx"
  
  curl_download(url = url,
                destfile = "data/raw/ons_deaths_data_2023.xlsx")
  
}


ons_all_deaths_age <-
  read.xlsx(
    xlsxFile = "data/raw/ons_deaths_data_2023.xlsx",
    sheet = 6,
    sep.names = "_",
    startRow = 5
  ) |>
  janitor::clean_names()

ons_all_deaths_age <-
  ons_all_deaths_age |>
  as_tibble()

period_age_all_cause <-
  # Deaths from all causes by year and age only.
  ons_all_deaths_age |>
  filter(
    year_of_registration %in% c(2022, 2023),
    sex != "All people",
    age_years != "All ages",
    marital_status == "All marital statuses"
  ) |>
  group_by(year_of_registration, age_years) |>
  summarise(number_of_deaths = sum(number_of_deaths))

period_age_all_cause <-
  period_age_all_cause |>
  mutate(age_years = as.numeric(str_remove(age_years, "[^\\d]+"))) |>
  # Exclude ages outside the range of the drug-related deaths figures
  filter(age_years >= 18) |>
  filter(age_years <= 100)


period_age_drd <-
  # All deaths we've categorized as drug-related, by year and age
  df |>
  group_by(period, age) |>
  summarise(count = sum(count))

colnames(period_age_all_cause) <-
  c("period", "age", "deaths_all_causes")

colnames(period_age_drd) <-
  c("period", "age", "deaths_related_to_drug_misuse")


data <-
  left_join(period_age_all_cause, period_age_drd) |>
  mutate(deaths_related_to_drug_misuse = zoo::na.fill(deaths_related_to_drug_misuse, 0))

cut_age_groups <- 
function(x) {
  cut(
    x,
    breaks = c(17, 24, 34, 44, 54, 64, 74, 84, Inf),
    labels = c(
      "18-24",
      "25-34",
      "35-44",
      "45-54",
      "55-64",
      "65-74",
      "75-84",
      "85+"
    ),
    right = TRUE
  )
}


p1 <- 
data |> 
  filter(age < 55) |> 
  rowwise() |> 
  mutate(age_group = cut_age_groups(age)) |> 
  group_by(period, age_group) |> 
  summarise(across(where(is.numeric), sum)) |> 
  group_by(age_group) |> 
  summarise(across(where(is.numeric), mean)) |>  
  mutate(deaths_all_causes = deaths_all_causes - deaths_related_to_drug_misuse) |> 
  pivot_longer(cols = deaths_all_causes:deaths_related_to_drug_misuse) |> 
  mutate(name = str_replace(name, "all_causes", "all_other_causes")) |> 
  mutate(name = snakecase::to_sentence_case(name)) |> 
  ggplot(aes(x = age_group, y = value)) + 
  geom_col(aes(fill = name), position = "fill", width = 0.5, colour = "black", alpha = 0.8) + 
  theme_bw() + 
  scale_y_continuous(labels = scales::percent) +
  scale_fill_manual(values = af_colours(n = 4)) + 
  theme(
    text = element_text(size =50),
    legend.position = "bottom",
    legend.justification = "left"
  ) + 
  labs(fill = NULL,
       y = NULL,
       x = "Age group",
       title = "Deaths related to drug misuse as a proportion of all deaths, by age group",
       subtitle = "Average of 2022 and 2023 data"
       )


p1

```

------------------------------------------------------------------------

```{r, fig.height=15, fig.width=30}
p2 <- 
data |> 
  filter(age < 55) |> 
  rowwise() |> 
  mutate(age_group = cut_age_groups(age)) |> 
  group_by(period, age_group) |> 
  summarise(across(where(is.numeric), sum)) |> 
  group_by(age_group) |> 
  summarise(across(where(is.numeric), mean)) |>  
  mutate(deaths_all_causes = deaths_all_causes - deaths_related_to_drug_misuse) |> 
  pivot_longer(cols = deaths_all_causes:deaths_related_to_drug_misuse) |> 
  mutate(name = str_replace(name, "all_causes", "all_other_causes")) |> 
  mutate(name = snakecase::to_sentence_case(name)) |> 
  ggplot(aes(x = age_group, y = value)) + 
  geom_col(aes(fill = name), position = "identity", width = 0.5, colour = "black", alpha = 0.8) + 
  theme_bw() + 
  scale_y_continuous(labels = scales::comma) +
  scale_fill_manual(values = af_colours(n = 4)) + 
  theme(
    text = element_text(size = 50),
    legend.position = "bottom",
    legend.justification = "left"
  ) + 
  labs(fill = NULL,
       y = NULL,
       x = "Age group",
       title = "Count of deaths related to drug misuse deaths by all other causes\nby age group",
       subtitle = "Average of 2022 and 2023 data"
  )
p2
```

------------------------------------------------------------------------

# Years of life lost due to substance use

---

- Years of life lost (YLL) is a measure of the impact of premature mortality, helpfully defined by Public Health England [here](https://fingertips.phe.org.uk/static-reports/health-profile-for-england/definitions-regional.html#years-of-life-lost-yll).

- Chudasama et al. (2022) investigated five methods for estimating YLL.[^2] The first two methods are feasible with the available data for YLL from drug use and alcohol-specific deaths. Only the drug-related YLL could be segmented by geographical estimates of deprivation.

- All five methods are detailed in the supplementary PDF [here](https://ars.els-cdn.com/content/image/1-s2.0-S0895435622001639-mmc1.pdf).

---

The crude expected years of life lost is:

$$
YLL = (D_x)(e_{x}^s)
$$

Where \(D_x\) is the number of deaths and \(e_{x}^s\) is the standard age of death from the external life expectancy.

- Total YLL associated with drug use was **215,148**. 

- For comparison, Heald et al. (2024)[^3] estimated the impact of obesity on YLL in England at 791,689 in 2019.

- Deaths at ages between 35 and 54 accounted for 60% of the YLL.

![](plots/yll_plot_age_group.png)

---
## YLL method II

"Extend the basic formula and incorporate the discounting and age weighting requires a range of assumptions."

$$

Y_x = d_x \left[\frac{KCe^{r}(n^{a}x)}{(r + \beta)^2}(e^{z}[-(r+\beta)^{2}(e^{s}_{x} + a_{x}) - 1]) \right]

$$


The discount rate ($r$) is based on the concept that future mortality outcomes
have less value than the present mortality outcomes. The discounting rate is a
continuous and exponential function, therefore the higher the discount rate,
the lower the impact is in the future. The usual time for the discounting rate
is set at 3% per year. Age weighting is used because it enables the YLL
equation to value differently at various ages, i.e. the value of years lived
by a young adult could be noted to be worth more than that of a very young or
an old individual. So, the parameters were estimated at $C=0.1658$ and $\beta=0.04$
to account for this, which was calibrated and chosen from the 1990 GBD
study.[1, 2]


Aragón, T.J., Lichtensztajn, D.Y., Katcher, B.S., Reiter, R. and Katz, M.H., 2007. Calculating Expected Years of Life Lost to Rank the Leading Causes of Premature Death in San Francisco. San Francisco Department of Public Health.
### References

[^1]: The criteria for this classification are described in `Box 2` of the `Definition` tab of the latest release of *Deaths related to drug poisoning, England and Wales*. Available [here](https://www.ons.gov.uk/peoplepopulationandcommunity/birthsdeathsandmarriages/deaths/datasets/deathsrelatedtodrugpoisoningenglandandwalesreferencetable).

[^2]: Chudasama, Y.V., Khunti, K., Gillies, C.L., Dhalwani, N.N., Davies, M.J., Yates, T., & Zaccardi, F. (2022). Estimates of years of life lost depended on the method used: tutorial and comparative investigation. *Journal of Clinical Epidemiology*, 150, pp. 42–50. Available at: [https://doi.org/10.1016/j.jclinepi.2022.06.012](https://doi.org/10.1016/j.jclinepi.2022.06.012) [Accessed 6 Nov. 2024].

[^3]: Heald, A., Stedman, M., Fryer, A.A., Davies, M.B., Rutter, M.K., Gibson, J.M., & Whyte, M. (2024). Counting the lifetime cost of obesity: Analysis based on national England data. *Diabetes, Obesity and Metabolism*, 26(4), pp. 1464–1478.