NM-MRI
02_medication
02_medication
clinical data pulled: 2020-11-16
code written: 2021-11-05
last ran: 2021-11-17
Description. This script performs all analyses related to psychotropic medication. Note that medication data is pulled from the screening visit (not baseline), which is further from scan, as the medication history assessment was not administered at baseline until midway through the study.
Load libraries and data
#clear environment
rm(list = ls())
#list required libraries
packages <- c(
'tidyverse',
'reshape2', #wrangling
'kableExtra',
'ggrepel'
)
#load required libraries
lapply(packages, require, character.only = TRUE)
#read in medication data
meds <- read.csv('../clinical/SenDep_MedicationData_2021.04.15.csv')
#adjust medication IDs, so match participants' list
meds$id <- substr(meds$id, 11, 16)
#read in cleaned participant demographic/clinical/cognition/lc data
df <- read.csv(dir('../clinical', full.names=T, pattern="^df_2021")) #48Wrangle data
#keep only participants we need
meds <- meds[meds$id %in% df$id,] #139
#keep medication only from screening arm -- closest to scan
meds <- meds[meds$redcap_event_name == 'screening_arm_1',] #48
#for clarity, keep only the data we need
meds <- meds[, grepl('^id$|pmeds_meds_yn|pmeds_drug|pmeds_continue', names(meds))]
#recode if taking medication or not
meds$pmeds_meds_yn <- ifelse(meds$pmeds_meds_yn == 0, 'no', 'yes') # if take medication or not
#omit data (rows) if not taking pmed
meds <- meds[!meds$pmeds_meds_yn == 'no',] #only 38 -- meaning some don't have any meds
#melt -- easier to manipulate
meds <- reshape2::melt(meds, id.vars=c('id','pmeds_meds_yn'))
#delete rows with a blank 'value'
meds <- meds[!meds$value=="", ]
#delete rows that are all NA
meds <- meds[rowSums(is.na(meds)) != ncol(meds), ]
#rename drug started / ended variables
meds$variable <- gsub('started$', 'started_month', meds$variable)
meds$variable <- gsub('started2$', 'started_year', meds$variable)
#get drug number in a new column
meds$drug_number <- as.numeric(stringr::str_extract(meds$variable, "[[:digit:]]+"))
#remove drug number from the variables
meds$variable <- gsub('[[:digit:]]+', '', meds$variable)
#also rename 'continuing variables, for clarity
meds$variable <- gsub('pmeds_drug_continuing', 'pmeds_drug_continuingSameDose', meds$variable)
#remove 'continue' variable we don't need
meds <- meds[!meds$variable == 'pmeds_continue',]
#create a mutate condition
mutate_cond <- function(.data, condition, ..., envir = parent.frame()) {
condition <- eval(substitute(condition), .data, envir)
.data[condition, ] <- .data[condition, ] %>% mutate(...)
.data
}Recode medication variables
#recode drug name
meds <- meds %>%
mutate_cond(variable == 'pmeds_drug_name', value = recode(value,
'2' = 'Irbesartan',
'6' = 'Fosinopril',
'9'= 'Perindopril',
'16' = 'Losartan',
'18' = 'Telmisartan',
'30' = 'Propranolol',
'39' = 'Spironolactone',
'49' = 'Amlodipine',
'60' = 'Chlorthalidone',
'70' = 'Diltiazem',
'80' = 'Hydrochlorothiazide',
'94' = 'Pravastatin',
'96' = 'Rosuvastatin',
'202' = 'Dabigatran',
'302' = 'Hydroxyzine',
'311' = 'Bimatoprost',
'312' = 'Budesonide/Formoterol',
'313'= 'Budesonide/Formoterol',
'502' = 'Calcium Carbonate',
'504' = 'Esomeprazole',
'509' = 'Pantoprazole',
'511'= 'Rabeprazole',
'937' = 'Levothyroxine Sodium',
'940' = 'Metformin',
'1000' = 'OTHER',
'1256' = 'Acetaminophen',
'1259' = 'ASA',
'1267' = 'Meloxicam',
'1269' = 'Naproxen',
'1300' = 'B Complex',
'1303' = 'Vitamin C',
'1304' = 'Vitamin D',
'1308' = 'Calcium',
'1310' = 'Fish oil',
'1316' = 'Iron Supplement',
'1317' = 'Magnesium',
'1318' = 'Melatonin',
'1319' = 'Multi-Vitamin',
'1321' = 'Omega 3',
'2006' = 'Pregabalin',
'2015' = 'Sertraline',
'2019' = 'Venlafaxine',
'2027' = 'Escitalopram',
'2030' = 'Paroxetine',
'2033' = 'Mirtazapine',
'2043' = 'Olanzapine',
'2048' = 'Lorazepam',
'2055' = 'Zolpidem',
'2056' = 'Zopiclone',
'2060' = 'Amitryptiline',
'2068' = 'Bupropion',
'2074'= 'Citalopram',
'2090' = 'Methylphenidate'
))
#recode drug indication
meds <- meds %>%
mutate_cond(variable == 'pmeds_drug_indication', value = recode(value,
'6' = 'High Blood Pressure',
'7' = 'High Cholesterol',
'8' = 'Hypertension',
'200' = 'Iron Deficiency',
'400'= 'Asthma',
'500' = 'Acid Reflux',
'800'= 'Arthritis',
'802' = 'Muscle and Joint Pain',
'803'= 'Osteoporosis',
'895' = 'Menopause',
'896'= 'Low Testosterone',
'897'= 'Hypothyroidism',
'899' = 'Diabetes Type II',
'1000' = 'OTHER',
'1001' = 'UNKNOWN',
'1200' = 'Allergies',
'1203' = 'Inflammation',
'1206' = 'Pain Relief',
'1300' = 'Supplement',
'2000' = 'ADHD',
'2002' = 'Anxiety',
'2003' = 'Depression',
'2005' = 'Migraines',
'2012' = 'Sleep',
'2013' = 'Tremor'
))
#recode drug class
meds <- meds %>%
mutate_cond(variable == 'pmeds_drug_class', value = recode(value,
'1' = 'ACE Inhibitor',
'2' = 'Angiotensin Receptor Blocker',
'5' = 'Beta Blocker',
'7' = 'Calcium Channel Blocker',
'8' = 'Diuretic',
'14'= 'Statin',
'300' = 'Antihistamine',
'301' = 'Bronchodilator',
'504' = 'Proton Pump Inhibitor',
'806' = 'Steroid',
'902' = 'Hormonal Replacement Therapy',
'905' = 'Monoclonal Antibody',
'906' = 'Oral Anti-Hyperglycemic Agent',
'908' = 'Thyroid Agent',
'1000' = 'OTHER',
'1001' = 'UNKNOWN',
'1201' = 'AntiInflammatory',
'1203' = 'Analgesic',
'1204' = 'NSAID',
'1300' = 'Supplement',
'2000' = 'Anticholinergic Agent',
'2001' = 'Anticonvulsant',
'2002' = 'Antidepressant',
'2004' = 'Atypical Antipsychotic',
'2005' = 'Benzodiazepine',
'2009' = 'Non-Benzodiazepine Hypnotic Agent',
'2010' = 'Stimulant'
))
#recode drug frequency
meds <- meds %>%
mutate_cond(variable == 'pmeds_drug_frequency', value = recode(value,
'1' = 'QD',
'2' = 'QAM',
'3' = 'QHS',
'4' = 'BID',
'5' = 'TID',
'7' = 'Q1wk',
'8' = 'Q2wk',
'12' = 'PRN',
'1000' = 'OTHER',
'1001' = 'UNKNOWN'
))
#recode start month
meds <- meds %>%
mutate_cond(variable == 'pmeds_drug_started_month', value = recode(value,
'1' = 'winter',
'2' = 'spring',
'3' = 'summer',
'4' = 'fall',
'5' = 'January',
'6' = 'Feburary',
'7' = 'March',
'8' = 'April',
'9' = 'May',
'11' = 'July',
'13' = 'September',
'14' = 'October',
'15' = 'November',
'16' = 'December',
'1000' = 'UNKNOWN'
))
#recode drug units
meds <- meds %>%
mutate_cond(variable == 'pmeds_drug_units', value = recode(value,
'1' = 'Milligram',
'2' = 'Drop',
'3' = 'Puff',
'4' = 'Spray/Squirt',
'6' = 'Tablet',
'8' = 'Microgram',
'9' = 'Gram',
'19' = 'International Units',
'1001' = 'UNKNOWN'
))
#recode continuing same dose
meds <- meds %>%
mutate_cond(variable == 'pmeds_drug_continuingSameDose', value = recode(value,
'1' = 'Yes',
'2' = 'No'
))
#move back data into wide format
meds <- dcast(meds, id + pmeds_meds_yn + drug_number ~ variable, value.var='value')
#move 'other' columns back into main column
meds$pmeds_drug_name <- ifelse(meds$pmeds_drug_name == 'OTHER', meds$pmeds_drug_name_other, meds$pmeds_drug_name)
meds$pmeds_drug_class <- ifelse(!is.na(meds$pmeds_drug_class_other), meds$pmeds_drug_class_other, meds$pmeds_drug_class)
meds$pmeds_drug_frequency <- ifelse(!is.na(meds$pmeds_drug_freq_othersp), meds$pmeds_drug_freq_othersp, meds$pmeds_drug_frequency)
meds$pmeds_drug_indication <- ifelse(!is.na(meds$pmeds_drug_indication_other), meds$pmeds_drug_indication_other, meds$pmeds_drug_indication)
#delete the 'other' columns
meds <- meds %>% select(-contains("_other"))Identify drugs that are not relevant
#now, list drugs not relevant (not antidepressant or beta-blocker)
omit_indication <- c('Acid Reflux',
'Allergies',
'Anti-Diarrhea',
'Arthritis',
'Diabetes Type II',
'Digestive problems (Pancreas)',
'Diuretic',
'Erectile Dysfunction',
'High Cholesterol', #all statins
'Hypothyroidism',
'Immune System',
'Inflammation',
'Iron Deficiency',
'Low Testosterone',
'Menopause',
'Migraines',
'Muscle and Joint Pain',
'Osteoporosis',
'Pain Relief',
'Sleep',
'Supplement')
#omit the drugs with those indications
meds <- meds[!meds$pmeds_drug_indication %in% omit_indication, ] #51
#manually review what's left
table(meds$pmeds_drug_class)
#next, omit from remaining drugs on the basis of class
omit_class <- c(
'ACE Inhibitor', #these are not beta-blockers
'Angiotensin Receptor Blocker', #these overlap with ACE inhibitors...
'Anticholinergic Agent',
'Benzodiazepine',
'Calcium Channel Blocker',
'Diuretic',
'NSAID',
'Supplement',
'Thiazide', #urine,
'Thiazide-like diurectic/ACE inhibitor')
#omit the drugs
meds <- meds[!meds$pmeds_drug_class %in% omit_class, ] #22
#also omit the 'Bimatoprost' (glaucoma)
meds <- meds[!meds$pmeds_drug_name == 'Bimatoprost', ] #21
#now that we have the medications, get them into some standard format
meds <- subset(meds, select=-c(
pmeds_meds_yn,
pmeds_drug_continuingSameDose,
pmeds_drug_started_month,
pmeds_drug_started_year,
drug_number
))Manual adjustments, as required
#manual adjustments, as required- first, the asthma medications to puffs -- looks like inconsistent
meds$pmeds_drug_strength <- ifelse(meds$id == 'SEN039' & meds$pmeds_drug_name == 'umeclidinium and vilanterol', 1, meds$pmeds_drug_strength)
meds$pmeds_drug_units<- ifelse(meds$id == 'SEN039' & meds$pmeds_drug_name == 'umeclidinium and vilanterol', 'Puff', meds$pmeds_drug_units)
meds$pmeds_drug_strength <- ifelse(meds$id == 'SEN011' & meds$pmeds_drug_name == 'Budesonide/Formoterol', 2, meds$pmeds_drug_strength)
meds$pmeds_drug_units<- ifelse(meds$id == 'SEN011' & meds$pmeds_drug_name == 'Budesonide/Formoterol', 'Puff', meds$pmeds_drug_units)
meds$pmeds_drug_frequency <- ifelse(meds$id == 'SEN011' & meds$pmeds_drug_name == 'Budesonide/Formoterol', 'BID', meds$pmeds_drug_frequency)
#fix class error
meds$pmeds_drug_class <- ifelse(meds$pmeds_drug_name == 'Bricanyl', 'Bronchodilator', meds$pmeds_drug_class)
meds$pmeds_drug_class <- ifelse(meds$pmeds_drug_name == 'Budesonide/Formoterol', 'Long-acting beta agonist', meds$pmeds_drug_class)
#fix drug name
meds$pmeds_drug_name<- ifelse(meds$pmeds_drug_name == 'umeclidinium and vilanterol', 'Umeclidinium/Vilanterol', meds$pmeds_drug_name)Standardize drug strength, frequency, type
#make sure that drug strength is numeric
meds$pmeds_drug_strength <- as.numeric(meds$pmeds_drug_strength)
#also make sure that the units are standardized to QD (once a day)
meds$pmeds_drug_strength <- ifelse(meds$pmeds_drug_frequency == 'BID', meds$pmeds_drug_strength * 2, meds$pmeds_drug_strength)
meds$pmeds_drug_frequency <- ifelse(meds$pmeds_drug_frequency == 'BID', 'QD', meds$pmeds_drug_frequency)
#check unique antidepressants in sample
unique(meds$pmeds_drug_name[meds$pmeds_drug_class == 'Antidepressant'])
#recode
meds <- meds %>%
mutate(antidepressantType = case_when(
#SSRI
pmeds_drug_name == 'Sertraline' ~ "SSRI",
pmeds_drug_name == 'Escitalopram' ~ "SSRI",
pmeds_drug_name == 'Paroxetine' ~ "SSRI",
pmeds_drug_name == 'Citalopram' ~ "SSRI",
#NDRI
pmeds_drug_name == 'Bupropion' ~ "NDRI",
#Tricyclic
pmeds_drug_name == 'Amitryptiline' ~ "Tricyclic",
#NaSSA
pmeds_drug_name == 'Mirtazapine' ~ "NaSSA",
#SNRI
pmeds_drug_name == 'Venlafaxine' ~ "SNRI"))Medication tables
[1] All depression medication.
The following table summarises all of the depression medications taken by participants in the study. As expected, only LLD participants are taking depression medications. In total, 14 depression medications (8 unique) are taken by 13 unique participants (i.e., one participant is taking 2 depression medications).
To ensure that antidepressant drug is not confounding LC NM-MRI contrast, we opted to run alternate analyses with participants taking noradrenergic antidepressants excluded. Thus, we will retain participants taking only SSRIs, but omit those taking any NaSSA (noradrenergic and specific serotonergic antidepressant), NDRI (norepinephrine and dopamine reuptake inhibitors), SNRI (serotonin and norepinephrine reuptake inhibitors), and Tricyclic medications (increase levels of norepinephrine and serotonin). Thus, we will run an alternate analysis with n=4 unique LLD participants excluded on the basis of their depression medication. The participants in question are SEN014, SEN020, SEN057, SEN061.
#add Diagnosis information to table
meds <- merge(df[, c('id', 'Diagnosis')], meds, by='id')
#table
df_medDepressionTable <- meds %>%
filter(pmeds_drug_class == 'Antidepressant') %>%
select(c(id, Diagnosis, pmeds_drug_name, antidepressantType, pmeds_drug_strength))
#rename columns for table
df_medDepressionTable <- df_medDepressionTable %>%
dplyr::rename(
'drug name' = pmeds_drug_name,
'antidepressant type' = antidepressantType,
'dose (mg/day)' = pmeds_drug_strength
)
#put into table (sortable)
DT::datatable(df_medDepressionTable,
rownames = FALSE,
options= list(pageLength = nrow(df_medDepressionTable)))[2] All beta agonists/antagonists.
The following table summarises noradrenergic medications taken for non-depression purposes, i.e., beta agonists and antagonists. Thus, this table includes both HC and LLD. Note that I have omitted ACE Inhibitors and Angiotensin Receptor Blockers, though it is possible that these may impact norepinephrine. We may run an alternate analysis excluding the n=5 unique participants taking a total of 6 beta-blocking medications (5 unique). The participants in question are SEN002, SEN010, SEN011, SEN039, SEN053. Note that the dose unit for the Asthma medications is ‘PUFFS’, not mg/day.
df_medBetaTable <- meds %>%
filter(pmeds_drug_class != 'Antidepressant',
pmeds_drug_class != 'Atypical Antipsychotic') %>%
select(c(id, Diagnosis, pmeds_drug_name, pmeds_drug_indication, pmeds_drug_strength))
#rename columns for table
df_medBetaTable <- df_medBetaTable %>%
dplyr::rename(
'drug name' = pmeds_drug_name,
'indication' = pmeds_drug_indication,
'dose (mg/day)' = pmeds_drug_strength
)
#put into table (sortable)
DT::datatable(df_medBetaTable,
rownames = FALSE,
options= list(pageLength = nrow(df_medBetaTable)))[3] Additional drug-based exclusions.
We may also run an alternate analysis excluding one participant who is taking an anti-psychotic. That participant is SEN061.
#identify the antipsychotic
df_medOtherTable <- meds %>%
filter(pmeds_drug_class == 'Atypical Antipsychotic') %>%
select(c(id, Diagnosis, pmeds_drug_name, pmeds_drug_indication, pmeds_drug_strength))
#rename columns for table
df_medOtherTable <- df_medOtherTable %>%
dplyr::rename(
'drug name' = pmeds_drug_name,
'indication' = pmeds_drug_indication,
'dose (mg/day)' = pmeds_drug_strength
)
#put into table (sortable)
DT::datatable(df_medOtherTable, rownames = FALSE)#participants to omit
dep <- unique(df_medDepressionTable$id[df_medDepressionTable$`antidepressant type` != 'SSRI']) #4
beta <- unique(df_medBetaTable$id) #5
other <- unique(df_medOtherTable$id) #1
#figure out who shared
beta %in% dep
other %in% dep #SEN061
other %in% beta
#bind together and get unique
participantsOmit <- unique(c(dep, beta, other))
#write out participant ID, and reason for exclusion, to a txt file
participantsOmit_txt <- file('../clinical/participantOmit.txt')
writeLines(paste(unique(dep), 'depression med'), participantsOmit_txt)
writeLines(paste(unique(beta), 'beta-blocker, non-depression med'), participantsOmit_txt)
writeLines(paste(unique(other), 'antipsychotic'), participantsOmit_txt)
close(participantsOmit_txt)Summary. The comprehensive list of the n=9 unique participants to consider excluding (7 LLD, 2 HC), in alternate analysis are: SEN014, SEN020, SEN057, SEN061, SEN002, SEN010, SEN011, SEN039, SEN053.
Medication and LC outliers
We also wanted to understand if the participants taking noradrenergic medications are more likely to show outlying LC NM-MRI values.
#pull imaging data out of df
lc <- df[, grep('^id$|Diagnosis|avg_max', names(df))]
#function to identify outliers
is_outlier <- function(x) {
return(x < quantile(x, 0.25) - 1.5 * IQR(x) | x > quantile(x, 0.75) + 1.5 * IQR(x))
}
#melt data for ggplot
lc <- reshape2::melt(lc, id.var=c('id', 'Diagnosis'))
#indicate if outlier
contingency <- lc %>%
group_by(variable) %>%
mutate(outlier = ifelse(is_outlier(value), value, as.numeric(NA)))
#get ready for plot
lc_outlier <- contingency %>%
drop_na(outlier)
#make a separate df containing only participants taking noradrenergic medications
lc_noradrenergic <- lc[lc$id %in% participantsOmit,]
#plot
ggplot(lc, aes(x=id, y=value)) +
geom_point(aes(colour=Diagnosis)) +
theme_bw() +
theme(legend.position = 'top',
axis.title = element_blank(),
panel.grid.major.x= element_blank(),
axis.text.x = element_blank(),
axis.ticks = element_blank()) + #text(angle=45, hjust=1),
facet_wrap(~variable) +
#now, add labels re: noradrengic medication particiapnts
geom_label_repel(
data=lc_noradrenergic,
aes(label=id),
box.padding=.5,
alpha=.5) +
#overlay black points re: medication
geom_point(data=lc_outlier, aes(x=id, y=value), size=3) 
Proportion test.
We see that we have n=10 outlying values across the 6 segments (4 HC and 6 LLD). We can summarize the number of outlying LC values as a function of medication taken, in the following basic summary table:
##
## not outlier outlier Sum
## noradrenergic med 112 14 126
## other med 476 70 546
## Sum 588 84 672We can also review the column-wise contingencies, i.e., the percentage of participants with outlying LC values that are/are not taking noradrenergic medication:
##
## not outlier outlier
## noradrenergic med 19.04762 16.66667
## other med 80.95238 83.33333Lastly, we can perform a statistical test to determine if the row and column values are independent, i.e., if medication type is independent of outlying LC values. Because we have fewer than 5 observations in some bins, we perform Fisher’s exact test as an alternative to Pearson’s Chi-squared. We find that those taking noradrenergic medication are not more likely to have outlying LC values. Thus, we see that participants taking nonadrengic medications are not more likely to have outlying LC NM-MRI values.
##
## Fisher's Exact Test for Count Data
##
## data: contingency$noradrenergic and contingency$LC_outlier
## p-value = 0.6569
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
## 0.6279986 2.3456840
## sample estimates:
## odds ratio
## 1.176192Summary table
The data in the table below summarises the depression-related and beta agonist/antagonist medication taken by participants (n=number of participants, percent=percent of LLD sample, mean/min/max refer to mg/QD.)
Depression medication
#make sure that dosage information is numeric
df_medDepressionTable$`dose (mg/day)` <- as.numeric(df_medDepressionTable$`dose (mg/day)`)
df_medDepressionTable %>%
group_by(`drug name`) %>%
summarise(n=n(),
percent = round(n()/(sum(df$Diagnosis == 'LLD')) * 100,2),
mean = mean(`dose (mg/day)`, na.rm=TRUE),
min = min(`dose (mg/day)`),
max = max(`dose (mg/day)`)) %>%
kable() %>%
kable_styling()| drug name | n | percent | mean | min | max |
|---|---|---|---|---|---|
| Amitryptiline | 1 | 4 | 20 | 20 | 20 |
| Bupropion | 2 | 8 | 300 | 150 | 450 |
| Citalopram | 1 | 4 | 15 | 15 | 15 |
| Escitalopram | 2 | 8 | 20 | 20 | 20 |
| Mirtazapine | 1 | 4 | 15 | 15 | 15 |
| Paroxetine | 1 | 4 | 20 | 20 | 20 |
| Sertraline | 5 | 20 | 78 | 50 | 150 |
| Venlafaxine | 1 | 4 | 75 | 75 | 75 |
Beta agonists/antagonists
| drug name | n | percent | mean | min | max |
|---|---|---|---|---|---|
| Bricanyl | 1 | 4 | 1 | 1 | 1 |
| Budesonide/Formoterol | 2 | 8 | 4 | 4 | 4 |
| Methylphenidate | 1 | 4 | 10 | 10 | 10 |
| Propranolol | 1 | 4 | 20 | 20 | 20 |
| Umeclidinium/Vilanterol | 1 | 4 | 1 | 1 | 1 |