Drug utilization in Sweden and Denmark

Part II

In this post (Part II), I will look into the Danish national trends on of sedatives, antidepressants, anti-psychotics and anxiolytics utilization in women. Then I’d like to compare Danish national trends with Swedish national trends from Part I.

library(tidyverse)
library(magrittr)
library(wesanderson)

Dataset

First things first - uploading data.

link_list_dk <- list(
  "1996_atc_code_data.txt" = "https://medstat.dk/da/download/file/MTk5Nl9hdGNfY29kZV9kYXRhLnR4dA==",
  "1997_atc_code_data.txt"  = "https://medstat.dk/da/download/file/MTk5N19hdGNfY29kZV9kYXRhLnR4dA==",
  "1998_atc_code_data.txt" = "https://medstat.dk/da/download/file/MTk5OF9hdGNfY29kZV9kYXRhLnR4dA==",
  "1999_atc_code_data.txt" = "https://medstat.dk/da/download/file/MTk5OV9hdGNfY29kZV9kYXRhLnR4dA==",
  "2000_atc_code_data.txt" = "https://medstat.dk/da/download/file/MjAwMF9hdGNfY29kZV9kYXRhLnR4dA==",
  "2001_atc_code_data.txt" = "https://medstat.dk/da/download/file/MjAwMV9hdGNfY29kZV9kYXRhLnR4dA==",
  "2002_atc_code_data.txt" = "https://medstat.dk/da/download/file/MjAwMl9hdGNfY29kZV9kYXRhLnR4dA==",
  "2003_atc_code_data.txt" = "https://medstat.dk/da/download/file/MjAwM19hdGNfY29kZV9kYXRhLnR4dA==",
  "2004_atc_code_data.txt" = "https://medstat.dk/da/download/file/MjAwNF9hdGNfY29kZV9kYXRhLnR4dA==",
  "2006_atc_code_data.txt" = "https://medstat.dk/da/download/file/MjAwNl9hdGNfY29kZV9kYXRhLnR4dA==",
  "2007_atc_code_data.txt" = "https://medstat.dk/da/download/file/MjAwN19hdGNfY29kZV9kYXRhLnR4dA==",
  "2008_atc_code_data.txt" = "https://medstat.dk/da/download/file/MjAwOF9hdGNfY29kZV9kYXRhLnR4dA==",
  "2009_atc_code_data.txt" = "https://medstat.dk/da/download/file/MjAwOV9hdGNfY29kZV9kYXRhLnR4dA==",
  "2010_atc_code_data.txt" = "https://medstat.dk/da/download/file/MjAxMF9hdGNfY29kZV9kYXRhLnR4dA==",
  "2011_atc_code_data.txt" = "https://medstat.dk/da/download/file/MjAxMV9hdGNfY29kZV9kYXRhLnR4dA==",
  "2012_atc_code_data.txt" = "https://medstat.dk/da/download/file/MjAxMl9hdGNfY29kZV9kYXRhLnR4dA==",
  "2013_atc_code_data.txt" = "https://medstat.dk/da/download/file/MjAxM19hdGNfY29kZV9kYXRhLnR4dA==",
  "2014_atc_code_data.txt" = "https://medstat.dk/da/download/file/MjAxNF9hdGNfY29kZV9kYXRhLnR4dA==",
  "2015_atc_code_data.txt" = "https://medstat.dk/da/download/file/MjAxNV9hdGNfY29kZV9kYXRhLnR4dA==",
  "2016_atc_code_data.txt" = "https://medstat.dk/da/download/file/MjAxNl9hdGNfY29kZV9kYXRhLnR4dA==",
  "2017_atc_code_data.txt" = "https://medstat.dk/da/download/file/MjAxN19hdGNfY29kZV9kYXRhLnR4dA==",
  "2018_atc_code_data.txt" = "https://medstat.dk/da/download/file/MjAxOF9hdGNfY29kZV9kYXRhLnR4dA==",
  "2019_atc_code_data.txt" = "https://medstat.dk/da/download/file/MjAxOV9hdGNfY29kZV9kYXRhLnR4dA==",
  "atc_code_text.txt" = "https://medstat.dk/da/download/file/YXRjX2NvZGVfdGV4dC50eHQ=",
  "atc_groups.txt" = "https://medstat.dk/da/download/file/YXRjX2dyb3Vwcy50eHQ=",
  "population_data.txt" = "https://medstat.dk/da/download/file/cG9wdWxhdGlvbl9kYXRhLnR4dA=="
)

# ATC codes are stored in files
seq_atc <- c(1:23)
atc_code_data_list <- link_list_dk[seq_atc]

I will compile data on population size, which I may use to compute denominators for drug utilization per capita.

# age structure data
# parse with column type character
age_structure <- read_delim(link_list_dk[["population_data.txt"]], delim = ";", col_names = c(paste0("V", c(1:7))), col_types = cols(V1 = col_character(), V2 = col_character(), V3 = col_character(), V4 = col_character(), V5 = col_character(), V6 = col_character(), V7 = col_character())) %>%
  # rename and keep columns
  select(year = V1,
         region_text = V2,
         region = V3,
         gender = V4,
         age = V5,
         denominator_per_year = V6) %>% 
  # human-reading friendly label on sex categories
  mutate(
      gender_text = case_when(
      gender == "0" ~ "both sexes",
      gender == "1" ~ "men",
      gender == "2" ~ "women",
      T ~ as.character(gender)
    )
  ) %>% 
  # make numeric variables
  mutate_at(vars(year, age, denominator_per_year), as.numeric) %>% 
  # keep only data for the age categories I will work with and only on women
  filter(age %in% 20:54, gender == "2") %>% 
  arrange(year, age, region, gender)

age_structure %>% slice(1:10)
## # A tibble: 10 x 7
##     year region_text region gender   age denominator_per_year gender_text
##    <dbl> <chr>       <chr>  <chr>  <dbl>                <dbl> <chr>      
##  1  1996 DK          0      2         20                36363 women      
##  2  1996 DK          0      2         21                36366 women      
##  3  1996 DK          0      2         22                36498 women      
##  4  1996 DK          0      2         23                38452 women      
##  5  1996 DK          0      2         24                37873 women      
##  6  1996 DK          0      2         25                36383 women      
##  7  1996 DK          0      2         26                36115 women      
##  8  1996 DK          0      2         27                37514 women      
##  9  1996 DK          0      2         28                40669 women      
## 10  1996 DK          0      2         29                43619 women

Danish data has a handy file with all drug classes in English, which I intend to use as labels later on.

eng_drug <- read_delim(link_list_dk[["atc_groups.txt"]], delim = ";", col_names = c(paste0("V", c(1:6))), col_types = cols(V1 = col_character(), V2 = col_character(), V3 = col_character(), V4 = col_character(), V5 = col_character(), V6 = col_character())) %>% 
  # keep drug classes labels in English
  filter(V5 == "1") %>% 
  select(ATC = V1,
         drug = V2,
         unit_dk = V4)

eng_drug %>% slice(1:10)
## # A tibble: 10 x 3
##    ATC     drug                                    unit_dk                  
##    <chr>   <chr>                                   <chr>                    
##  1 A       Alimentary tract and metabolism         No common unit for volume
##  2 A01     Stomatological preparations             No common unit for volume
##  3 A01A    Stomatological preparations             No common unit for volume
##  4 A01AA   Caries prophylactic agents              Packages                 
##  5 A01AA01 Sodium fluoride                         Packages                 
##  6 A01AB   Antiinfectives for local oral treatment No common unit for volume
##  7 A01AB02 Hydrogenperoxide                        <NA>                     
##  8 A01AB03 Chlorhexidine                           DDD                      
##  9 A01AB04 Amphotericin                            DDD                      
## 10 A01AB09 Miconazole                              DDD

I keep only ATC codes I am interested in. I use $ in the end of each regular expression (regex) to make sure that I capture the widest ATC code category, which already encompasses all the subcategories. Hypnotics and sedatives - N05C already aggregates N05CA, N05CB, N05CC, N05CD, N05CE, N05CF, N05CH ATC codes.

regex_antidepress <- "^N06A$"
regex_antipsych <- "^N05A$"
regex_anxiolyt <- "^N05B$"
regex_sedat <- "^N05C$"

# combine all ATC codes in one regex string
all_regex <- paste(regex_antidepress, regex_antipsych, regex_anxiolyt, regex_sedat, sep = "|")
# ATC data
atc_data <- map(atc_code_data_list, ~read_delim(file = .x, delim = ";", trim_ws = T, col_names = c(paste0("V", c(1:14))), col_types = cols(V1 = col_character(), V2 = col_character(), V3 = col_character(), V4 = col_character(), V5 = col_character(),V6 = col_character(), V7 = col_character(), V8 = col_character(), V9 = col_character(), V10 = col_character(), V11 = col_character(), V12 = col_character(), V13 = col_character(), V14 = col_character()))) %>% 
  # bind data from all years
  bind_rows() 

# age categories: in 6 year bands and in each year
age_cat_in_6 <- atc_data %>% distinct(V6) %>% filter(row_number() %in% 2:7) %>% pull()
age_cat_in_6
## [1] "00-17" "18-24" "25-44" "45-64" "65-79" "80+"
age_in_years <- atc_data %>% distinct(V6) %>% filter(row_number() %in% 8:103) %>% pull()
age_in_years
##  [1] "000" "001" "002" "003" "004" "005" "006" "007" "008" "009" "010" "011"
## [13] "012" "013" "014" "015" "016" "017" "018" "019" "020" "021" "022" "023"
## [25] "024" "025" "026" "027" "028" "029" "030" "031" "032" "033" "034" "035"
## [37] "036" "037" "038" "039" "040" "041" "042" "043" "044" "045" "046" "047"
## [49] "048" "049" "050" "051" "052" "053" "054" "055" "056" "057" "058" "059"
## [61] "060" "061" "062" "063" "064" "065" "066" "067" "068" "069" "070" "071"
## [73] "072" "073" "074" "075" "076" "077" "078" "079" "080" "081" "082" "083"
## [85] "084" "085" "086" "087" "088" "089" "090" "091" "092" "093" "094" "95+"
# clean dataset
atc_data %<>% 
  # rename and keep columns
  rename(
    ATC = V1,
    year = V2,
    sector = V3,
    region = V4,
    gender = V5, # 0 - both; 1 - male; 2 - female
    age = V6,
    number_of_people = V7,
    patients_per_1000_inhabitants = V8,
    turnover = V9,
    regional_grant_paid = V10,
    quantity_sold_1000_units = V11,
    quantity_sold_units_per_unit_1000_inhabitants_per_day = V12,
    percentage_of_sales_in_the_primary_sector = V13
  ) %>% 
  # clean columns names and set-up labels
  mutate(
    year = as.numeric(year),
    region_text = case_when(
      region == "0" ~ "DK",
      region == "1" ~ "Region Hovedstaden",
      region == "2" ~ "Region Midtjylland",
      region == "3" ~ "Region Nordjylland",
      region == "4" ~ "Region Sjælland",
      region == "5" ~ "Region Syddanmark",
      T ~ NA_character_
    ),
    gender_text = case_when(
      gender == "0" ~ "both sexes",
      gender == "1" ~ "men",
      gender == "2" ~ "women",
      T ~ as.character(gender)
    ),
    age_cats = case_when(
      age %in% "A" ~ "A",
      age %in% age_cat_in_6 ~ "6 years",
      age %in% age_in_years ~ "yearly categories",
      age %in% "T" ~ "T",
      T ~ NA_character_
    )
  ) %>% 
  mutate_at(vars(turnover, regional_grant_paid, quantity_sold_1000_units, quantity_sold_units_per_unit_1000_inhabitants_per_day, number_of_people, patients_per_1000_inhabitants), as.numeric) %>% 
  select(-V14) %>%
  filter(
  age_cats == "yearly categories",
  gender == "2",
  region == "0",
  sector =="0"
) %>% 
  filter(str_detect(ATC, all_regex)) %>% 
  # deal with non-numeric age and re-code as numeric
  mutate(
    age = case_when(
      age == "95+" ~ "95",
      T ~ as.character(age)
    ),
    age = as.numeric(age)
    ) %>% 
  filter(age %in% 20:54) %>% 
  mutate(
    age_cat = case_when(
      age %in% 20:24 ~ "20-24",
      age %in% 25:34 ~ "25-34",
      age %in% 35:44 ~ "35-44",
      age %in% 45:54 ~ "45-54",
      T ~ NA_character_
    )
) %>% 
  left_join(age_structure) %>% 
  group_by(ATC, year, age_cat) %>% 
  mutate(
    numerator = sum(number_of_people),
    denominator = sum(denominator_per_year),
    patients_per_1000_inhabitants = numerator / denominator * 1000
  ) %>% 
  filter(row_number() == 1) %>% 
  ungroup() %>% 
  select(year, ATC, gender, age_cat, patients_per_1000_inhabitants) %>% 
  mutate(
    country = "Denmark"
  )
## Joining, by = c("year", "region", "gender", "age", "region_text", "gender_text")
atc_data %>% slice(1:10)
## # A tibble: 10 x 6
##     year ATC   gender age_cat patients_per_1000_inhabitants country
##    <dbl> <chr> <chr>  <chr>                           <dbl> <chr>  
##  1  1999 N05A  2      20-24                            5.60 Denmark
##  2  1999 N05A  2      25-34                            9.71 Denmark
##  3  1999 N05A  2      35-44                           19.9  Denmark
##  4  1999 N05A  2      45-54                           29.4  Denmark
##  5  1999 N05B  2      20-24                           12.5  Denmark
##  6  1999 N05B  2      25-34                           26.3  Denmark
##  7  1999 N05B  2      35-44                           61.7  Denmark
##  8  1999 N05B  2      45-54                          106.   Denmark
##  9  1999 N05C  2      20-24                            7.86 Denmark
## 10  1999 N05C  2      25-34                           18.0  Denmark

Combining Danish and Swedish data

# combine Swedish data from previous post and Danish data
data <- data_se %>% 
  select(year, ATC, gender, age_cat, patients_per_1000_inhabitants) %>% 
  mutate(
    country = "Sweden"
  ) %>% 
  full_join(atc_data) %>% 
  # join with ATC classes labels in English
  left_join(eng_drug)
## Joining, by = c("year", "ATC", "gender", "age_cat", "patients_per_1000_inhabitants", "country")
## Joining, by = "ATC"

To avoid copy-pasting the same ggplot code, I use a wrapper around the plotting function.

plot_utilization <- function(data, drug_regex, title){
  data %>% 
  filter(str_detect(ATC, {{ drug_regex }})) %>% 
  group_by(country) %>% 
  ggplot(aes(x = year, y = patients_per_1000_inhabitants, color = country)) +
  geom_path(position = position_dodge(width = 1)) +
  facet_grid(rows = vars(drug, ATC), cols = vars(age_cat), scales = "free", drop = T) +
  theme_light(base_size = 14) +
  scale_x_continuous(breaks = c(seq(1999, 2019, 10))) +
  scale_color_manual(values = wes_palette(name = "GrandBudapest1", type = "discrete")) +
  theme(plot.caption = element_text(hjust = 0, size = 10),
        legend.position = "bottom",
        panel.spacing = unit(0.8, "cm")) +
  labs(y = "Female patients\nper 1,000 women in the population", title = paste0(title, " utilization in Swedish and Danish women"), subtitle = "aged between 20 and 54 years", caption = "dataviz by Elena Dudukina Twitter @evpatora\nSources: sdb.socialstyrelsen.se; medstat.dk")
}

list_regex <- list(regex_antidepress, regex_antipsych, regex_anxiolyt, regex_sedat)
list_title <- list("Antidepressants", "Antipsychotics", "Anxiolytics", "Sedatives")
list_plots <- map2(list_regex, list_title, ~plot_utilization(data = data, drug_regex = .x, title = .y))

Results

list_plots[[1]]
### Antidepressants Among Danish women aged between 20 and 54 years, antidepressants utilizatioin was steadily decreasing, while utilization among Swedish women was increasing. In Denmark, 56 women aged between 20-24 and 72 women aged between 35-45 per 1000 women in the population received antidepressants prescription.
list_plots[[2]]

Antipsychotics

Antipsychotics utilization was increasing over the years in both Denmark and Sweden. The rates of antipsychotics utilization were higher among Danish women when compared to Swedish women in the same age categories.

list_plots[[3]]

Anxiolytics

Anxiolytics utilization rates were higher among Swedish women. In Denmark, anxiolytics utilization was decreasing since 1999 until 2019 for women in all age categories between 20 and 54 years.

list_plots[[4]]

Sedatives

The rates of sedatives utilization were higher for Swedish women when compared with Danish women of the same age groups. Utilization of sedatives among Danish women was below 50 women per 1000 women of the same age group in the population in all investigated age categories. Among Swedish women, the rate of sedatives utilization was below 50 women per 1000 women in the population only for women between 20 and 24 years of age.

Take home message

  • Countries similar in many ways may exhibit different data patterns
  • Antidepressants, antipsychotics, anxiolytics, and sedatives utilization among Danish and Swedish women is far from identical

Thanks for tuning in 🤓

All feedback is appreciated ✌️

Elena Dudukina
Elena Dudukina
Clinical Specialist and PhD student

I am interested in women’s health, reproductive epidemiology, pharmacoepidemiology, causal inference, directed acyclic graphs, and R stats.

Related