Too few medicines to treat neuropathic pain on national essential medicines lists
Peter Kamerman
11 February 2017Map: Inclusion of medicines to treat neuropathic pain on national essential medicines lists
Interactive plot: Choose the drug class or IMF income category using the radio buttons in the top-right corner, and click on the country of interest to get more information.
############################################################
# #
# Load raw data #
# #
############################################################
access_data <- read_csv('./_data/2017-02-10-eml-drug-map/2015-neml-data.csv')
############################################################
# #
# Prepare data #
# #
############################################################
# Create filtering objects
drug_class <- c('Anticonvulsant', 'SNRI', 'TCA', 'Opioid', 'Topical')
drug_name <- c('Gabapentin', 'Pregabalin', 'Duloxetine', 'Venlafaxine',
'Amitriptyline', 'Clomipramine', 'Desipramine', 'Imipramine',
'Nortriptyline', 'Morphine', 'Oxycodone', 'Tramadol',
'Capsaicin', 'Lidocaine')
first_line <- c('Gabapentin', 'Pregabalin', 'Duloxetine', 'Venlafaxine',
'Amitriptyline', 'Clomipramine', 'Desipramine', 'Imipramine',
'Nortriptyline')
simplify_tca <- 'Amitriptyline|Clomipramine|Desipramine|Imipramine|Nortriptyline'
second_line <- c('Tramadol', 'Capsaicin', 'Lidocaine')
third_line <- c('Morphine', 'Oxycodone')
# Remove extraneous drugs/countries
access_data <- access_data %>%
# Select required columns
select(Country, Class, Drug, Listed) %>%
# Get rid of Sweden (OECD High Income)
filter(Country != 'Sweden') %>%
# Make lower case
rename(country = Country,
class = Class,
drug = Drug,
listed = Listed) %>%
## filter by drug name
filter(drug %in% drug_name)
# Get country names from access_data df for expanding first_line df (see below)
countries_access <- access_data$country
#-- Extract first-line drugs per country -------------------------------------#
first_line <- access_data %>%
filter(drug %in% first_line) %>%
mutate(grade = rep('first_line', nrow(.))) %>%
filter(listed == 'Yes') %>%
# Simplify TCAs
mutate(drug = str_replace(drug,
pattern = simplify_tca,
replacement = 'Tricyclic')) %>%
select(country, grade, class, drug) %>%
group_by(country, grade, class, drug) %>%
summarise(n = n()) %>%
# Flesh-out countries so that they all have all second-line drugs
complete(country = unique(countries_access), drug) %>%
unique(.) %>%
# Convert NA to 0
mutate(n = str_replace_na(n)) %>%
mutate(n = as.numeric(str_replace(n,
pattern = 'NA',
replacement = '0'))) %>%
# Recode all values > 0 as 1 (i.e., yes or no)
mutate(n = ifelse(test = n > 0, yes = 1, no = 0)) %>%
# Spread and flatten across drug classes
spread(drug, n) %>%
ungroup() %>%
select(-class) %>%
mutate(Tricyclic = str_replace_na(Tricyclic),
Gabapentin = str_replace_na(Gabapentin),
Pregabalin = str_replace_na(Pregabalin),
Duloxetine = str_replace_na(Duloxetine),
Venlafaxine = str_replace_na(Venlafaxine)) %>%
mutate(Tricyclic = as.numeric(str_replace(Tricyclic,
pattern = 'NA',
replacement = '0')),
Gabapentin = as.numeric(str_replace(Gabapentin,
pattern = 'NA',
replacement = '0')),
Pregabalin = as.numeric(str_replace(Pregabalin,
pattern = 'NA',
replacement = '0')),
Duloxetine = as.numeric(str_replace(Duloxetine,
pattern = 'NA',
replacement = '0')),
Venlafaxine = as.numeric(str_replace(Venlafaxine,
pattern = 'NA',
replacement = '0'))) %>%
# Calculate totals per drug per country
# (used to reduce rows from 345 to 115 after expansion of '0' rows
# caused by previous steps)
group_by(grade, country) %>%
summarise(tricyclic = sum(Tricyclic),
gabapentin = sum(Gabapentin),
pregabalin = sum(Pregabalin),
duloxetine = sum(Duloxetine),
venlafaxine = sum(Venlafaxine)) %>%
# Calculate class totals
mutate(class_tca = ifelse(tricyclic > 0,
yes = 1, no = 0),
class_alpha2delta = ifelse(gabapentin + pregabalin > 0,
yes = 1, no = 0),
class_snri = ifelse(duloxetine + venlafaxine > 0,
yes = 1, no = 0)) %>%
# Calculate classes per country
mutate(classes_first = class_tca + class_alpha2delta + class_snri) %>%
# Convert 0 and 1 to 'No' and 'Yes' for the drug columns (for popup)
mutate(tricyclic = ifelse(tricyclic == 1,
yes = 'yes', no = 'no'),
gabapentin = ifelse(gabapentin == 1,
yes = 'yes', no = 'no'),
pregabalin = ifelse(pregabalin == 1,
yes = 'yes', no = 'no'),
duloxetine = ifelse(duloxetine == 1,
yes = 'yes', no = 'no'),
venlafaxine = ifelse(venlafaxine == 1,
yes = 'yes', no = 'no')) %>%
# Rename for merging to SPDF
rename(name = country) %>%
# Reorder columns
select(name, grade, tricyclic, gabapentin, pregabalin, duloxetine,
venlafaxine, class_tca, class_alpha2delta,
class_snri, classes_first) %>%
# Get rid of strange attr that have cropped-up
as.data.frame()
#-- Extract second-line drugs per country ------------------------------------#
second_line <- access_data %>%
filter(drug %in% second_line) %>%
mutate(grade = rep('second_line', nrow(.))) %>%
filter(listed == 'Yes') %>%
select(country, grade, class, drug) %>%
group_by(country, grade, class, drug) %>%
summarise(n = n()) %>%
# Flesh-out countries so that they all have all second-line drugs
complete(country = countries_access, drug) %>%
unique(.) %>%
# Convert NA to 0
mutate(n = str_replace_na(n)) %>%
mutate(n = as.numeric(str_replace(n,
pattern = 'NA',
replacement = '0'))) %>%
# Recode all values > 0 as 1 (i.e., yes or no)
mutate(n = ifelse(test = n > 0, yes = 1, no = 0)) %>%
# Spread and flatten across drug classes
spread(drug, n) %>%
ungroup() %>%
select(-class) %>%
mutate(Lidocaine = str_replace_na(Lidocaine),
Tramadol = str_replace_na(Tramadol)) %>%
mutate(Lidocaine = as.numeric(str_replace(Lidocaine,
pattern = 'NA',
replacement = '0')),
Tramadol = as.numeric(str_replace(Tramadol,
pattern = 'NA',
replacement = '0'))) %>%
# Calculate totals per drug per country
group_by(grade, country) %>%
summarise(lidocaine = sum(Lidocaine),
tramadol = sum(Tramadol)) %>%
# Calculate class totals
mutate(class_opioid2 = ifelse(tramadol > 0,
yes = 1, no = 0),
class_topical = ifelse(lidocaine > 0,
yes = 1, no = 0)) %>%
# Calculate classes per country
mutate(classes_second = class_opioid2 + class_topical) %>%
# Convert 0 and 1 to 'No' and 'Yes' for the drug columns (for popup)
mutate(tramadol = ifelse(tramadol == 1,
yes = 'yes', no = 'no'),
lidocaine = ifelse(lidocaine == 1,
yes = 'yes', no = 'no')) %>%
# Make sure the is a zero for the classes_total (for plotting)
complete(classes_second = c(0, 1, 2)) %>%
# Rename for merging to SPDF
rename(name = country) %>%
# Reorder columns
select(name, grade, tramadol, lidocaine, class_opioid2,
class_topical, classes_second) %>%
# Get rid of strange attr that have cropped-up
as.data.frame()
#-- Extract third-line drugs per country -------------------------------------#
third_line <- access_data %>%
filter(drug %in% third_line) %>%
mutate(grade = rep('third_line', nrow(.))) %>%
filter(listed == 'Yes') %>%
select(country, grade, class, drug) %>%
group_by(country, grade, class, drug) %>%
summarise(n = n()) %>%
# Flesh-out countries so that they all have all second-line drugs
complete(country = countries_access, drug) %>%
unique(.) %>%
# Convert NA to 0
mutate(n = str_replace_na(n)) %>%
mutate(n = as.numeric(str_replace(n,
pattern = 'NA',
replacement = '0'))) %>%
# Recode all values > 0 as 1 (i.e., yes or no)
mutate(n = ifelse(test = n > 0, yes = 1, no = 0)) %>%
# Spread and flatten across drug classes
spread(drug, n) %>%
ungroup() %>%
select(-class) %>%
mutate(Morphine = str_replace_na(Morphine), # None had Botox
Oxycodone = str_replace_na(Oxycodone)) %>%
mutate(Morphine = as.numeric(str_replace(Morphine,
pattern = 'NA',
replacement = '0')),
Oxycodone = as.numeric(str_replace(Oxycodone,
pattern = 'NA',
replacement = '0'))) %>%
# Calculate totals per drug per country
group_by(grade, country) %>%
summarise(morphine = sum(Morphine),
oxycodone = sum(Oxycodone)) %>%
# Calculate class totals
mutate(classes_third = ifelse(morphine > 0 | oxycodone > 0,
yes = 1, no = 0)) %>%
# Convert 0 and 1 to 'No' and 'Yes' for the drug columns (for popup)
mutate(morphine = ifelse(morphine == 1,
yes = 'yes', no = 'no'),
oxycodone = ifelse(oxycodone == 1,
yes = 'yes', no = 'no')) %>%
# Make sure the is a zero for the classes_total (for plotting)
complete(classes_third = c(0, 1, 2)) %>%
# Get rid of strange NA row that gets added when using 'complete'
filter(!is.na(country)) %>%
# Rename for merging to SPDF
rename(name = country) %>%
# Reorder columns
select(name, grade, morphine, oxycodone, classes_third) %>%
# Get rid of strange attr that have cropped-up
as.data.frame()
############################################################
# #
# Make geospacial object #
# #
############################################################
# Load geojson
geo_data <- geojson_read(
'./_data/2017-02-10-eml-drug-map/admin0_countries_50m.geojson',
what = 'sp')
# Other countries
## Filter out countries with data
geo_other <- geo_data
geo_other <- geo_other[!(geo_data$name %in% countries_access), 'name']
# Simplify SPDF to avoid NAs
## Filter geo_data by countries in first_line df
geo_data <- geo_data[geo_data$name %in% countries_access,
c('name', 'income_grp', 'gdp_md_est')]
## Filter *_line dfs by countries in geo_data
countries_geo <- geo_data$name
first_line <- first_line[first_line$name %in% countries_geo, ]
second_line <- second_line[second_line$name %in% countries_geo, ]
third_line <- third_line[third_line$name %in% countries_geo, ]
## Factor 'classes_*'
first_line$factor_first <- factor(first_line$classes_first,
levels = c('0', '1', '2', '3'),
ordered = TRUE)
second_line$factor_second <- factor(second_line$classes_second,
levels = c('0', '1', '2'),
ordered = TRUE)
third_line$factor_third <- factor(third_line$classes_third,
levels = c('0', '1', '2'),
ordered = TRUE)
## Clean-up 'income grp'
imf <- data.frame(name = geo_data$name, imf = geo_data$income_grp) %>%
mutate(imf = str_replace(imf,
pattern = '1. High income: OECD',
replacement = 'High income')) %>%
mutate(imf = str_replace(imf,
pattern = '2. High income: nonOECD',
replacement = 'High income')) %>%
mutate(imf = str_replace(imf,
pattern = '3. Upper middle income',
replacement = 'Upper-middle income')) %>%
mutate(imf = str_replace(imf,
pattern = '4. Lower middle income',
replacement = 'Lower-middle income')) %>%
mutate(imf = str_replace(imf,
pattern = '5. Low income',
replacement = 'Low income')) %>%
mutate(imf = factor(imf,
levels = c('Low income', 'Lower-middle income',
'Upper-middle income', 'High income'),
ordered = TRUE))
# Merge SPDF with *_line dfs
geo_data <- merge(x = geo_data, y = first_line,
by.x = 'name', by.y = 'name') %>%
merge(x = ., y = second_line,
by.x = 'name', by.y = 'name') %>%
merge(x = ., y = third_line,
by.x = 'name', by.y = 'name') %>%
merge(x = ., y = imf,
by.x = 'name', by.y = 'name')
# Remove 'income_grp' column
geo_data@data <- geo_data@data[-2]
# Palette (from Color Brewer single hue sequential palette)
pal_first <- colorFactor(c('#feedde','#fdbe85','#fd8d3c','#d94701'),
geo_data$factors_first, ordered = TRUE)
pal_second <- colorFactor(c('#efedf5','#bcbddc','#756bb1'),
geo_data$factor_second, ordered = TRUE)
pal_third <- colorFactor(c('#e5f5e0','#a1d99b','#31a354'),
geo_data$factor_third, ordered = TRUE)
pal_income <- colorFactor(c('#eff3ff','#bdd7e7','#6baed6','#2171b5'),
geo_data$imf, ordered = TRUE)
# Popups
popup_first <- paste0('<b>Country: </b>', geo_data$name,
' (', geo_data$imf, ')',
'<br><b>No. 1st-line classes on NEML: </b>',
geo_data$classes_first,
'<br><em>TCA (any): </em>',
geo_data$tricyclic,
'<br><em>α2δ-ligand (gabapentin): </em>',
geo_data$gabapentin,
'<br><em>α2δ-ligand (pregabalin): </em>',
geo_data$pregabalin,
'<br><em>SNRI (duloxetine): </em>',
geo_data$duloxetine,
'<br><em>SNRI (venlafaxine): </em>',
geo_data$venlafaxine)
popup_second <- paste0('<b>Country: </b>', geo_data$name,
' (', geo_data$imf, ')',
'<br><b>No. 2nd-line classes on NEML: </b>',
geo_data$classes_second,
'<br><em>Weak Opioid (Tramadol): </em>',
geo_data$tramadol,
'<br><em>Topical agent (5% lidocaine): </em>',
geo_data$lidocaine,
'<br><em>Topical agent (8% capsaicin): </em> no')
popup_third <- paste0('<b>Country: </b>', geo_data$name,
' (', geo_data$imf, ')',
'<br><b>No. 3rd-line classes on NEML: </b>',
geo_data$classes_third,
'<br><em>Strong opioid (morphine): </em>',
geo_data$morphine,
'<br><em>Strong opioid (oxycodone): </em>',
geo_data$oxycodone,
'<br><em>Other (Botulinium toxin A): </em> no')
popup_NA <- paste0('<b>Country: </b>', geo_other$name,
'<br><b>No data</b>')
############################################################
# #
# Plot #
# #
############################################################
leaflet(data = geo_data,
width = '100%',
height = 600) %>%
setView(lat = 10, lng = 0, zoom = 2) %>%
# Add other countries
addPolygons(data = geo_other, weight = 1, color = '#777777',
opacity = 1, fillColor = '#777777',
fillOpacity = 1, popup = popup_NA) %>%
# Add analysed countries
addPolygons(weight = 1, color = '#000000',
opacity = 1, fillColor = ~ pal_first(factor_first),
fillOpacity = 1, popup = popup_first,
group = 'First-line drug classes') %>%
addPolygons(weight = 1, color = '#000000',
opacity = 1,fillColor = ~ pal_second(factor_second),
fillOpacity = 1, popup = popup_second,
group = 'Second-line drug classes') %>%
addPolygons(weight = 1, color = '#000000',
opacity = 1, fillColor = ~ pal_third(factor_third),
fillOpacity = 1, popup = popup_third,
group = 'Third-line drug classes') %>%
addPolygons(weight = 1, color = '#000000',
opacity = 1, fillColor = ~ pal_income(imf),
fillOpacity = 1, group = 'IMF income group') %>%
addLayersControl(baseGroups = c('First-line drug classes',
'Second-line drug classes',
'Third-line drug classes',
'IMF income group'),
options = layersControlOptions(collapsed = FALSE)) %>%
# Added new class to leafletfix.css - .info-small with smaller font size
addLegend("bottomright",
pal = pal_income,
values = ~ imf,
title = "International Monetary Fund (IMF) <br>income category",
opacity = 1) %>%
addLegend("bottomleft",
pal = pal_third,
values = ~ factor_third,
title = "No. of 3rd-line classes on NEML <br><em>(strong opioid, other)<sup>§</sup></em>",
opacity = 1) %>%
addLegend("bottomleft",
pal = pal_second,
values = ~ factor_second,
title = "No. of 2nd-line classes on NEML <br><em>(weak opioid, topical agents)<sup>‡</sup></em>",
opacity = 1) %>%
addLegend("bottomleft",
pal = pal_first, values = ~ factor_first,
title = "No. of 1st-line classes on NEML <br><em>(TCA, α2δ-ligand, SNRI)<sup>†</sup></em>",
opacity = 1)