---
title: "Assessing Canada’s Drug Shortage Problem"
output: 
  flexdashboard::flex_dashboard:
    orientation: rows
    source_code: embed
    self_contained: false
---

```{r setup, include=FALSE}
options(tidyverse.quiet=T)
library(flexdashboard)
library(tidyverse)
library(lubridate, warn.conflicts=F)
library(devtools, quietly=T)
library(rdrugshortages)     # devtools::install_github('pipitone/rdrugshortages')
library(kableExtra, warn.conflicts=F)
library(plotly, warn.conflicts=F)

# settings
DATA_DIR = 'data/'  # Load the Health Canada Drug Product Database
DPD_DOWNLOAD = getOption('dpd_download', default=F) # download the latest drug product database dataset
DSC_DOWNLOAD = getOption('dsc_download', default=F) # download latest DSC database
DSC_PATH     = 'data/dsc.csv'
DSC_START    = ymd("2018-01-01", tz="")
FAR_FUTURE   = ymd("2100-12-31", tz="")
FAR_PAST     = origin

theme_set(theme_minimal())
caption = list(labs(caption =  paste(
  paste(
    "source: drugshortagescanada.ca, last updated:", as_date(file.info(DSC_PATH)$ctime)), 
    "source: Health Canada Drug Product Database",
    "jon.pipitone.ca/research/drug-shortages", 
    sep="\n")))

# load data
if (DSC_DOWNLOAD) { 
  dsc = dsc_search()  %>%
    select_if(~ !is.list(.x)) %>%
    write_csv(DSC_PATH)
}
dsc.orig = read_csv(DSC_PATH, trim_ws = T, col_types = cols(
  status = col_factor(),
  drug.company.post_office_box = col_character(),
  hc_en_comments = col_character(), 
  hc_fr_comments = col_character(),
  shortage_reason.en_reason = col_factor()
  )) %>%
  select(-matches("^fr_|\\.fr_|_fr")) %>%
  mutate(
    reason = fct_recode(shortage_reason.en_reason,  
      "Demand increase" = "Demand increase for the drug.", 
      "Manufacturing distruption" = "Disruption of the manufacture of the drug.",
      "Ingredient shortage" = "Shortage of an active ingredient.", 
      "Ingredient shortage" = "Shortage of an inactive ingredient or component.",
      "Manufacturing practices" = "Requirements related to complying with good manufacturing practices.",
      "Other" = "Other (Please describe in comments)",
      "Shipping delay" = "Delay in shipping of the drug." 
    ) 
  ) %>%
  mutate(
    end_date_no_na = coalesce(actual_end_date, FAR_FUTURE))

## Parsual EML
# https://dx.doi.org/10.9778%2Fcmajo.20160122
eml = read_csv('data/cleanmeds.csv', trim_ws=T, col_types = cols()) %>% 
  select(category, medication, ATC)

## DPD
dpd = dpd_load(DATA_DIR, download=DPD_DOWNLOAD)

## Add the DPD ATC number to the DSC data
dsc = inner_join(dpd$drug, dpd$ther, by="DRUG_CODE") %>% 
  select(din = DRUG_IDENTIFICATION_NUMBER, dpd_atc_number = TC_ATC_NUMBER, TC_ATC) %>% 
  distinct(din, .keep_all = TRUE) %>%
  right_join(dsc.orig, by=c("din"))
dsc.shortages = dsc %>% filter(type.label == "shortage")
dsc.discontinuations = dsc %>% filter(type.label == "discontinuance")

this_year = interval(floor_date(now(), "year"), ceiling_date(now(), "year"))
week_start = floor_date(now(), "week")
this_week = interval(week_start, week_start + dweeks(1))
next_week = int_shift(this_week, by = dweeks(1))
```

Shortages
==========

Row {data-height=125}
--------------------------------

### Active shortages
```{r}
valueBox(dsc %>%  filter(status == "active_confirmed") %>% count() %>% pull(), 
         icon="fa-capsules", 
         color = "warning")
```

### New shortages expected this week
```{r}
.shortages_this_week =  dsc.shortages %>% 
  filter(
    status %in% c("active_confirmed", "anticipated_shortage"), 
    coalesce(actual_start_date, anticipated_start_date) %within% this_week)
valueBox(nrow(.shortages_this_week),icon="fa-needle")
```

### New discontinuations this week
```{r}
.active_disconts = dsc.discontinuations  %>% 
  filter(
    status == "discontinued", 
    discontinuation_date %within% this_week)
valueBox(nrow(.active_disconts), icon="fa-ban")
```

### Shortages of Essential Medicines this year
```{r}
.emls_with_shortages = dsc.shortages %>%
  filter(actual_start_date %within% this_year) %>%
  right_join(eml, by=c("dpd_atc_number" = "ATC"))  %>%
  filter(!is.na(din))  %>%
  group_by(category, medication, dpd_atc_number)
valueBox(nrow(.emls_with_shortages), icon="fa-ambulance")
```

Row  
----------------------------------
### Drugs in shortage
```{r}
weeks = tibble(
  week.start = seq(as_date(floor_date(DSC_START, "week")), 
                   as_date(floor_date(now(), "week")), by = "1 week"), 
  i = 1
)
p = dsc.shortages %>%
  mutate(i = 1) %>%
  full_join(weeks, by = "i") %>%
  filter(week.start %within% interval(actual_start_date, end_date_no_na))%>%
  group_by(week.start) %>%
  summarize(n = n()) %>%
  mutate(week = week.start) %>%
  ggplot( aes(x = week, y = n)) + 
  geom_line() + 
  scale_x_date(date_minor_breaks = "1 month", 
               date_labels="%b %Y") + 
  labs(x = NULL, y = NULL)
ggplotly(p)
```


### Causes of drug shortages this year {.no-padding}
```{r, fig.height=4, asp=0.6}
p = dsc.shortages %>%
  transmute(started = as_date(floor_date(actual_start_date, "month")), 
            reason = reason) %>%  
  filter(started >= ymd('2019-10-01'), started < now()) %>%
  count(started, reason) %>%
  group_by(started) %>%
  mutate(percent = round(n / sum(n) * 100, 1)) %>%
  transmute(month = started, percent, reason) %>%
  ggplot(aes(x = month, y = percent, fill=reason)) + 
  geom_bar(stat='identity') + 
  scale_x_date(date_breaks="1 month", 
               date_labels="%b %Y") +
  labs(x=NULL, y="% of total") + 
  theme(legend.position="bottom")
ggplotly(p) %>% layout(legend = list(orientation = 'h', x=0, y=-.1))
```


Row
-----------------------------------------------------------------------

### New shortages per week since 2018
```{r}
p = dsc.shortages %>%
  mutate(week = as_date(floor_date(actual_start_date, "week"))) %>%
  filter(week >= ymd('2018-01-01'), week < now()) %>%
  count(week) %>%
  ggplot(aes(x=week, y=n)) + 
  geom_bar(stat='identity') + 
  scale_x_date(date_minor_breaks = "1 month", 
               date_labels = "%b %Y") +  
  labs(x=NULL, y=NULL)
 ggplotly(p) 
```

### Shortage durations
```{r}
p = dsc.shortages %>% 
  filter(actual_start_date >= ymd("2018-01-01")) %>%
  mutate(duration = (actual_end_date - coalesce(actual_start_date, estimated_end_date))/dweeks(1), 
         month_str=floor_date(actual_start_date, unit="year")) %>%
  filter(!is.na(duration), duration > 0) %>%
  ggplot(aes(x=duration)) +
  geom_histogram(binwidth=4) +
  facet_wrap(~factor(month_str), nrow=1, labeller = as_labeller(function(x) format.Date(x, "%Y"))) + 
  labs(x="Duration (weeks)", y="# of shortages")
ggplotly(p)
```

COVID-19 {data-orientation=columns}
=======================================

```{r}
.dpd_status_ordered = dpd$status %>%
  group_by(DRUG_CODE) %>%
  arrange(HISTORY_DATE, .by_group = T) %>%
  mutate(.curr = 1:n(), .prev = .curr-1, .last = n()) %>%
  ungroup()

# complex join so that LHS has FROM status and RHS has TO status
# Also includes pair when LHS is the most recent status and then RHS is ignored
dpd_status_intervals = inner_join(.dpd_status_ordered, .dpd_status_ordered, by="DRUG_CODE") %>%
  select(-starts_with("CURRENT"), -starts_with("LOT"), -starts_with("EXP")) %>%
  filter(.prev.y == .curr.x | (.curr.x == .last.x & .curr.x == .curr.y)) %>%
  transmute(
    DRUG_CODE = DRUG_CODE, 
    start_date = HISTORY_DATE.x,
    end_date = if_else(.curr.x == .last.x, as_date(FAR_FUTURE), HISTORY_DATE.y), 
    status = STATUS.x)

# Map ATC code to NAME
atc_codes = dpd$ther %>% transmute(dpd_atc_number = TC_ATC_NUMBER, atc_name = TC_ATC) %>% distinct()

# date range to show 
FROM = floor_date(ymd("2019-11-01"), "week")
TO = as_date(now())
BY = "1 week"  # warning: if the FROM-TO range is large (e.g. > 20 years), change BY to something larger, e.g. 3 months

# now, for each month, let's compute the total number of drugs marketed, for each ATC code
dpd.months = tibble(month.date = seq(FROM, TO, by=BY),  n = 1)

dpd_num_marketed_by_atc = dpd_status_intervals %>% mutate(n=1) %>%
  inner_join(dpd.months, by="n") %>%
  filter(month.date %within% interval(start_date, end_date)) %>%
  select(DRUG_CODE, status, month.date) %>%
  inner_join(select(dpd$ther, DRUG_CODE, TC_ATC_NUMBER), by = "DRUG_CODE") %>%
  group_by(month.date, TC_ATC_NUMBER) %>%
  summarize(marketed = sum(status == "MARKETED"))

dsc_num_shortages_by_atc = dsc.shortages %>%
  mutate(n=1) %>%
  inner_join(dpd.months, by="n") %>%
  filter(month.date %within% interval(actual_start_date, end_date_no_na)) %>%
  distinct(month.date, dpd_atc_number, din) %>%
  group_by(month.date, dpd_atc_number) %>%
  summarize(shortages = n())


atc_plot = function(atc) {
  shortages = dsc_num_shortages_by_atc %>%
    filter(str_starts(dpd_atc_number, atc))
  
  marketed = dpd_num_marketed_by_atc  %>%
    filter(str_starts(TC_ATC_NUMBER, atc))
  
  right_join(shortages, marketed, 
             by = c("month.date" = "month.date", "dpd_atc_number" = "TC_ATC_NUMBER")) %>%
    ungroup() %>%
    mutate(shortages = replace_na(shortages, 0)) %>%
    mutate(available = marketed - shortages) %>%
    group_by(month.date) %>%
    summarize(Available = sum(available), Marketed = sum(marketed)) %>%
    pivot_longer(cols = Available:Marketed, names_to = "status", values_to = "count") %>%
    ggplot(aes(x = week, y = count, linetype = fct_rev(status))) +
    geom_line()  + 
    labs(x=NULL, y = "# of DINs", 
         linetype = NULL) + 
    expand_limits(y=0)  + 
    theme(legend.position="bottom")
}

integer_breaks <- function(n = 5, ...) { 
  function(x) { 
    breaks = floor(pretty(x, n, ...))
    names(breaks) = attr(breaks, "labels")
    breaks
  }
}


atc_group_plot = function(atcs) {
  atc = paste(atcs, collapse="|")
  shortages = dsc_num_shortages_by_atc %>%
    filter(str_starts(dpd_atc_number, atc))
  
  marketed = dpd_num_marketed_by_atc  %>%
    filter(str_starts(TC_ATC_NUMBER, atc))
  
  right_join(shortages, marketed, 
             by = c("month.date" = "month.date", "dpd_atc_number" = "TC_ATC_NUMBER")) %>%
    ungroup() %>%
    mutate(shortages = replace_na(shortages, 0)) %>%
    mutate(available = marketed - shortages) %>%
    group_by(month.date, dpd_atc_number) %>%
    summarize(Available = sum(available), Marketed = sum(marketed)) %>%
    ungroup() %>%
    group_by(dpd_atc_number) %>%
    filter(sum(Marketed-Available) > 0, Available > 0) %>%
    mutate(percent_available = Available/max(Available)*100) %>%
    mutate(week = month.date) %>%
    left_join(atc_codes, by="dpd_atc_number") %>%
    mutate(name_and_atc = paste0(atc_name, ' [', dpd_atc_number, ']')) %>%
    ggplot(aes(x = week, y = Available, fill=name_and_atc)) +
    geom_area()  + 
    scale_x_date(date_breaks="2 month", 
                 date_labels="%b %Y", 
                 minor_breaks=NULL) + 
    scale_y_continuous(breaks=integer_breaks(3)) + 
    facet_wrap(~name_and_atc, ncol=2, strip.position = "right", scales="free_y") + 
    theme(legend.position="none", 
          strip.text.y = element_text(angle=0)
          )
}
```


Column  {data-width=150}
------------------------------------- 

### COVID-19 Indicator Medicines
We are tracking the following list of medications (and groups of medications) to
assess for COVID-related shortages:

```{r}
# To add: propofol, rocuronium, cisatracurium, ketamine, midazolam, norepinephrine, phenylephrine, fentanyl
covid_meds = tribble(
  ~ATC,      ~Type,             ~Description,
  "N01A",    "Anesthetic use",  "Anesthetics (includes opioid, ketamines)",  
  "M03AC",   "Anesthetic use",  "Muscle relaxants (e.g. paralytics)", 
  "M03AB01", "Anesthetic use",  "Succinylcholine (suxamethonium)",
  "N05CD",   "Anesthetic use",  "Benzodiazepines", 
  "C01CA",   "Anesthetic use",  "Adrenergics",  
  "N05CM18", "Anesthetic use",  "Dexmedetomidine", 
  "H01BA",   "Anesthetic use",  "Vasopressin and analogues",
  "R03AC",   "Other",           "Respiratory Beta Agonists",  
  "R01AA04", "Other",           "Phenylephrine (nasal)",
  "C03CA01", "Other",           "Furosemide",
  "C01BD01", "Other",	          "Amiodarone",
  "J05AR10", "Other",           "Kaletra (anti-viral)", 
  "P01BA02", "Other",           "Hydroxychloroquine", 
  "N02BE01", "Other",           "Acetaminophen (Paracetamol)",  
)
   
kable(covid_meds %>% arrange(Type, ATC) %>% select(-Type)) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"), 
    full_width=T) 
```

The following charts shows the availability of COVID-19 indicator medications since November 2019. 
Only medications (by ATC code) that have had a shortages will be shown. 

Column {.tabset}
-----------------------------------------------------

### Anesthetic Use {data-padding=30}
```{r}
p = atc_group_plot(filter(covid_meds, Type == "Anesthetic use")$ATC)
p1 = p + 
    geom_vline(xintercept=as.numeric(ymd("2020-01-15")), linetype=1) + 
    geom_vline(xintercept=as.numeric(ymd("2020-03-09")), linetype=1, colour="red") + 
    labs(x=NULL, y = "# DIN available") 
p1_ly = ggplotly(p1) %>% 
  layout(annotations = list(x = 1, y = -0.06, xref="paper", yref="paper", 
                            xanchor="right", yanchor="bottom", xshift=0, yshift=0, 
                            font=list(size=10), showarrow=F,
                            text="Black line indicates first COVID-19 case in canada, red line indicates first death"))  
p1_ly
```

### Other Use {data-padding=30}
```{r}
p = atc_group_plot(filter(covid_meds, Type == "Other")$ATC)
p1 = p + 
    geom_vline(xintercept=as.numeric(ymd("2020-01-15")), linetype=1) + 
    geom_vline(xintercept=as.numeric(ymd("2020-03-09")), linetype=1, colour="red") + 
    labs(x=NULL, y = "# DIN available") 
p1_ly = ggplotly(p1) %>% 
  layout(annotations = list(x = 1, y = -0.06, xref="paper", yref="paper", 
                            xanchor="right", yanchor="bottom", xshift=0, yshift=0, 
                            font=list(size=10), showarrow=F,
                            text="Black line indicates first COVID-19 case in canada, red line indicates first death"))  
p1_ly
```


Search
======

Row 
------------------------------------- 
### Searchable list of shortage reports since January, 2023
```{r}
dsc.shortages %>%
  mutate(start_date = date(coalesce(actual_start_date, anticipated_start_date)),) %>%
  filter(start_date %within% interval(ymd("2023-01-01"), now()+dweeks(4))) %>% 
  arrange(start_date) %>%
  mutate(drug_name = if_else(str_length(en_drug_brand_name) > 30, paste0(substring(en_drug_brand_name, 0, 30), '...'), en_drug_brand_name), 
         company_name = str_to_title(company_name), 
         duration = round((start_date %--% estimated_end_date)/days(1)),
         status = fct_recode(status, 
                             "Active" = "active_confirmed", 
                             "Anticipated" = "anticipated_shortage", 
                             "Avoided" = "avoided_shortage", 
                             "Resolved" = "resolved"),
         link = text_spec(id, link=paste0("https://www.drugshortagescanada.ca/shortage/", id))) %>%
  select(dpd_atc_number, drug_name, drug_strength, company_name, status, start_date, duration, reason, link) %>%
  DT::datatable(escape=F, 
                filter="top",
                rownames=T,
                options = list(bPaginate = FALSE, 
                               autoWidth = T,
                               rowid = F, 
                               search = list(regex = T), 
                               order=list(list(6, "desc"), list(5,"asc"))), 
                colnames = c("ATC", "Drug Name", "Dose", "Company", "Status", "Start", "Duration (days)", "Reason", "View report")) %>%
  DT::formatStyle("start_date", "white-space"="nowrap") %>%
  DT::formatStyle("drug_name", "white-space"="nowrap") %>%
  DT::formatStyle("reason", "white-space"="nowrap") %>%
  DT::formatStyle('status', color = DT::styleEqual(c('Anticipated', 'Active'), c('blue', 'green')))
```

About
=====

Column
------

### About this website

Hello. I am Jon Pipitone, resident in Psychiatry at Queen's University, Kingston,
Ontario, Canada. This dashboard is my attempt to show interesting trends and
summaries of Canadian drug shortage reports. Code for this website can be found
on [github](https://github.com/pipitone/drug-shortages-dashboard). This work is
licensed under a <a rel="license"
href="http://creativecommons.org/licenses/by/4.0/">Creative Commons Attribution
4.0 International License</a>.

Much of this work has been done in collaboration with [Dr. Jacalyn
Duffin](https://jacalynduffin.ca/), who has a her own site dedicated to
analysis of the canadian drug shortages: https://canadadrugshortage.com. 

We published a peer-reviwed paper in 2018 on this issue:

* Donelle, Jessy, Jacalyn Duffin, Jonathan Pipitone, and Brian White-Guay.
“Assessing Canada’s Drug Shortage Problem.” CD Howe Institute Commentary 515, 2018. https://doi.org/10.2139/ssrn.3192558

I also keep [a blog](https://jon.pipitone.ca/blog) and I sometimes discuss drug
shortages.

Questions, comments, requests? Please send email to drugshortages@pipitone.ca.

### FAQ 

1. Why are there different totals listed under "active shortages" and on the "drugs in shortage" graph?

    In short, it's because there is missing or incorrect data in some of the
shortage reports that means they get excluded from historical counts. I've
written [a blog post about it](https://jon.pipitone.ca/blog/2020-04-18-how-many-drug-shortages-are-there/).

1. How often is the dashboard updated? 

    Daily. It was last updated `r now()`. 

Column  {data-width=100}
----------------------------

### Data sources

##### [drugshortagescanada.ca](https://www.drugshortagescanada.ca)
Database of shortage and discontinuation reports.  Active from March 2017 to
present. Manufacturers are required to report.

Last updated: `r now()`

##### [Drug Product Database](https://www.canada.ca/en/health-canada/services/drugs-health-products/drug-products/drug-product-database.html)
Database of historical records of drug products in Canada. Maintained by Health Canada.

Last updated: `r now()`

##### [CLEANMeds Essential Medicines List](https://cleanmeds.ca)