################################################################################

### USER VARIABLES #############################################################

################################################################################

setwd('R:/Pop_Health/Farr_Luscombe/')

### Relating to constructed cohort #############################################

# the R data file containing the patient cohort

cohort.file <- '2_Cohort/COHORT_SCAD_makecohort2.rda'

# columns to remove

kill.cols <- c('anonpatid', 'pracid','year_of_birth')

# columns to make relative to indexdate

date.cols <- c('afterentry_acs','afterentry_mi_nos',

            'afterentry_nstemi','afterentry_stemi','afterentry_ua',

            'crd','date_entry','date_exit','date_mi_endpoint','deathdate','dob','dod',

            'dod_combined','earliest_chd','earliest_hf','earliest_mi','earliest_sa',

            'earliest_ua','frd','hes_end','hes_start','indexdate','praclcd',

            'pracuts','tod','recent_acs','recent_mi_nos','recent_nstemi',

            'recent_stemi','recent_ua','smokdate','endpoint_coronary_date',

            'endpoint_death_date'

)

### Relating to other files ####################################################

data.path <- '1a_ExtractedData'

hes.file <- file.path(data.path, 'hes_diag_epi.csv')

n.top.icd <- 100

procedures.file <- file.path(data.path, 'hes_procedures.csv')

n.top.procedures <- 100

tests.files <-

  file.path(

    data.path,

    paste0('test.part.', 0:3)

  )

# Number of tests to extract

n.top.tests <- 100

# Find tests as close as possible to the first value here, and no more than the

# second timepoint before it

test.timepoints <- c(0, 183) #days

clinical.files <-

  file.path(

    data.path,

    paste0('clinical.part.', 0:3)

  )

n.top.clinical.values <- 100

n.top.clinical.history <- 100

therapy.files <-

  file.path(

    data.path,

    paste0('therapy.part.', 0:7)

  )

n.top.therapy <- 100

################################################################################

### END USER VARIABLES #########################################################

################################################################################

### Load the patient cohort to act as a base ###################################

# load the standard COHORT variable, which is a data table called COHORT

load(cohort.file)

### Do a bunch of silly fixes on the data ######################################

# cabg_6mo is 1s and 0s, which should be TRUE and FALSE

COHORT$cabg_6mo <- as.logical(COHORT$cabg_6mo)

# long_nitrate is 1s and NAs, which should be TRUE and FALSE...

COHORT$long_nitrate <- as.logical(COHORT$long_nitrate)

COHORT$long_nitrate[is.na(COHORT$long_nitrate)] <- FALSE

# pci_6mo is 1s and 0s, which should be TRUE and FALSE

COHORT$pci_6mo <- as.logical(COHORT$pci_6mo)

### Append other data types ####################################################

# We'll be needing handymedical from here on in

source('Andrew/lib/handymedical.R', chdir = TRUE)

require(CALIBERlookups)

percentMissing <- function(x, sf = 3) {

  round(sum(is.na(x))/length(x), digits = sf)*100

}

### Hospital episodes statistics ###############################################

# read in the hospital data (HES)

hes.diag.epi <- readMedicalData(

  hes.file,

  c("anonpatid", "epistart", "epiend", "icd"),

  c("integer", "date", "date", "factor")

)

# remove non alphanumerics (trailing -s on some ICD codes)

hes.diag.epi$icd <- gsub('[^[:alnum:]]', '', hes.diag.epi$icd)

# Now, merge with the indexdate...we only want data before then, because looking

# after it is cheating, and using all data rather than just stuff before may

# distort our choice of variables as some variables may be very common after

# entering the study, but less so before... (This actually isn't much of an

# issue...only 9 variables differ. Still!)

hes.diag.epi <-

  merge(

    hes.diag.epi, COHORT[, c('anonpatid', 'indexdate')],

    by = 'anonpatid', all = TRUE

  )

hes.diag.epi$relativedate <- hes.diag.epi$indexdate - hes.diag.epi$epistart

# Now, remove all the negative relative dates, because they're in the past

hes.diag.epi <- hes.diag.epi[hes.diag.epi$relativedate >= 0, ]

# get aggregate statistics for each ICD code

hes.by.icd <- hes.diag.epi[, length(unique(anonpatid)), by = icd]

names(hes.by.icd) <- c('icd', 'n.pat')

# Take the top n by number of patients with that code

top.icd <-

  hes.by.icd$icd[order(hes.by.icd$n.pat, decreasing = TRUE)[1:n.top.icd]]

# And now, let's put how far in the past the patient was first diagnosed with

# each of these things into the table...

# First, discard all the rows corresponding to non-top ICDs

hes.diag.epi.top <- hes.diag.epi[hes.diag.epi$icd %in% top.icd, ]

# Then, keep only the earliest instance of each per patient, because those are

# the ones we care about

hes.diag.epi.top.earliest <- 

  hes.diag.epi.top[, min(relativedate), by = c('anonpatid', 'icd')]

# Per ICD code, add a new column to the cohort and put in numbers

hes.diag.epi.top.earliest <-

  dcast(

    data = hes.diag.epi.top.earliest,

    formula = anonpatid ~ icd,

    value.var = "V1"

  )

# Add a prefix to the names to keep track

names(hes.diag.epi.top.earliest)[2:(n.top.icd + 1)] <-

  paste0('hes.icd.', names(hes.diag.epi.top.earliest)[2:(n.top.icd + 1)])

# Keep the names for when we need to make them relative to the indexdate when

# anonymising later

names.hes.diag.icd <- names(hes.diag.epi.top.earliest)[2:(n.top.icd + 1)]

# Now, merge with the original cohort

COHORT <-

  merge(

    COHORT,

    hes.diag.epi.top.earliest,

    by = c('anonpatid'),

    all = TRUE

  )

hes.icd.summary <-

  data.frame(

    percent.missing = 

      sort(

        sapply(

          COHORT[, names(COHORT)[startsWith(names(COHORT), 'hes.icd.')], with = FALSE],

          percentMissing

        )

      )

  )

hes.icd.summary$code <- substring(rownames(hes.icd.summary), 9)

hes.icd.summary <- merge(hes.icd.summary, CALIBER_DICT[, c('code', 'term')], by = 'code')

print(hes.icd.summary[order(hes.icd.summary$percent.missing),])

### Hospital procedures data ###################################################

# read in the hospital data (HES)

hes.procedures <- readMedicalData(

  procedures.file,

  c("anonpatid", "opcs", "evdate"),

  c("integer", "factor", "date")

)

# Now, merge with the indexdate...we only want data before then, because looking

# after it is cheating, and using all data rather than just stuff before may

# distort our choice of variables as some variables may be very common after

# entering the study, but less so before... (14 differ...)

hes.procedures <-

  merge(

    hes.procedures, COHORT[, c('anonpatid', 'indexdate')],

    by = 'anonpatid', all = TRUE

  )

hes.procedures$relativedate <- hes.procedures$indexdate - hes.procedures$evdate

# Now, remove all the negative relative dates, because they're in the past

hes.procedures <- hes.procedures[relativedate >= 0]

# get aggregate statistics for each OPCS code

hes.by.opcs <- hes.procedures[, length(unique(anonpatid)), by = opcs]

names(hes.by.opcs) <- c('opcs', 'n.pat')

# Take the top n by number of patients with that code

top.opcs <-

  hes.by.opcs$opcs[order(hes.by.opcs$n.pat, decreasing = TRUE)[1:n.top.procedures]]

# And now, let's put how far in the past the patient was first diagnosed with

# each of these things into the table...

# First, discard all the rows corresponding to non-top OPCS codes

hes.procedures.top <- hes.procedures[hes.procedures$opcs %in% top.opcs, ]

# Then, keep only the earliest instance of each per patient, because those are

# the ones we care about

hes.procedures.top.earliest <-

  hes.procedures.top[, min(relativedate), by = c('anonpatid', 'opcs')]

# Per OPCS code, add a new column to the cohort and put in numbers

hes.procedures.top.earliest <-

  dcast(

    data = hes.procedures.top.earliest,

    formula = anonpatid ~ opcs,

    value.var = "V1"

  )

# Add a prefix to the names to keep track

names(hes.procedures.top.earliest)[2:(n.top.procedures + 1)] <-

  paste0('hes.opcs.', names(hes.procedures.top.earliest)[2:(n.top.procedures + 1)])

# Keep the names for when we need to make them relative to the indexdate when

# anonymising later

names.procedures.opcs <- names(hes.procedures.top.earliest)[2:(n.top.icd + 1)]

# Now, merge with the original cohort

COHORT <-

  merge(

    COHORT,

    hes.procedures.top.earliest,

    by = c('anonpatid'),

    all = TRUE

  )

# Summarise by percent missing

hes.opcs.summary <-

  data.frame(

    percent.missing = 

      sort(

        sapply(

          COHORT[, names(COHORT)[startsWith(names(COHORT), 'hes.opcs.')], with = FALSE],

          percentMissing

        )

      )

  )

hes.opcs.summary$code <- substring(rownames(hes.opcs.summary), 10)

hes.opcs.summary <- merge(hes.opcs.summary, CALIBER_DICT[, c('code', 'term')], by = 'code')

print(hes.opcs.summary[order(hes.opcs.summary$percent.missing),])

### Test results ###############################################################

# read in the test data

tests.data <- readMedicalData(

  tests.files,

  # data1 is operator (eg =, > etc), data2 is value, data 3 is unit of measure

  c("anonpatid", "eventdate", "enttype", "data1", "data2", "data3"),

  c("integer", "date", "integer", "integer", "numeric", "integer")

)

# First, discard those where operator is not =, because > and < etc will

# introduce complexity, and drop data1 since it's now useless

tests.data <- tests.data[data1 == 3]

tests.data$data1 <- NULL

# Now, let's subtract the indexdate from every test so we can choose the ones

# closest to the desired dates...

tests.data <-

  merge(

    tests.data, COHORT[, c('anonpatid', 'indexdate')],

    by = 'anonpatid', all = TRUE

  )

tests.data$relativedate <- tests.data$indexdate - tests.data$eventdate

# Only keep positive values; negative ones are in the future which is cheating!

tests.data <- tests.data[relativedate >= 0]

# Only 89001 have any test results from before the indexdate!

# Discard any test results from times greater than the longest time ago to check

tests.data <- tests.data[relativedate < test.timepoints[2]]

# And this drops to 57972 in the preceding six months

# Per patient and per test, keep the smallest relativedate value only

tests.data <-

  tests.data[, .SD[which.min(relativedate)], by = list(anonpatid, enttype)]

# aggregate by test type (enttype, which covers a range of Read codes which can

# mean the same test) and unit of measure (so we can choose the tests with the

# units with the best coverage)

# We do this after all the preprocessing (ie only looking at the most recent

# value per patient before but not too far before the indexdate) because

# otherwise a lot of the top tests change. Presumably 

tests.by.test <-

  tests.data[, length(unique(anonpatid)), by = c('enttype', 'data3')]

names(tests.by.test) <- c('enttype', 'data3', 'n.pat')

top.tests <- tests.by.test$enttype[order(tests.by.test$n.pat, decreasing = TRUE)[1:n.top.tests]]

top.units <- tests.by.test$data3[order(tests.by.test$n.pat, decreasing = TRUE)[1:n.top.tests]]

# Now, discard all the rows corresponding to non-top tests

tests.data <-

  tests.data[

    # Needs to be exact match of enttype and variable combination

    paste0(enttype, '!!!', data3) %in% paste0(top.tests, '!!!', top.units)

  ]

# Make a column per test

tests.wide <-

  dcast(

    data = tests.data,

    formula = anonpatid ~ enttype + data3,

    value.var = "data2"

  )

# Add a prefix to the names to keep track

# (There may not be n.top.tests columns here as some get lost during the paring

# down processes above...)

names(tests.wide)[-1] <-

  paste0('tests.enttype.data3.', names(tests.wide)[-1])

# Now, merge with the original cohort

COHORT <-

  merge(

    COHORT,

    tests.wide,

    by = c('anonpatid'),

    all = TRUE

  )

### GP diagnosis data ##########################################################

# read in the GP data (clinical)

clinical.data <- readMedicalData(

  clinical.files,

  c("anonpatid", "eventdate", "medcode", "enttype", "data1", "data2", "data3", "data4", "data5", "data6", "data7"),

  c("integer", "date", "integer", "integer", "numeric", "numeric", "numeric", "numeric", "numeric", "numeric", "numeric")

)

# Now, merge with the indexdate...we only want data before then, because looking

# after it is cheating, and using all data rather than just stuff before may

# distort our choice of variables as some variables may be very common after

# entering the study, but less so before... (14 differ...)

clinical.data <-

  merge(

    clinical.data, COHORT[, c('anonpatid', 'indexdate')],

    by = 'anonpatid', all = TRUE

  )

clinical.data$relativedate <- clinical.data$indexdate - clinical.data$eventdate

# Now, remove all the negative relative dates, because they're in the past

clinical.data <- clinical.data[relativedate >= 0]

# Find out which enttypes have associated data values, to distinguish purely

# binary variables (like medical history, family history etc) from test results

# and so on which have data1, data2, etc values.

# From a quick scan, anything with any data at all has a data1 value, so we can

# use that as a proxy.

clinical.by.data1 <- clinical.data[, sum(!is.na(data1)), by = enttype]

# Split into two data tables, the ones with associated data, and without

clinical.history <-

  clinical.data[

    enttype %in% clinical.by.data1$enttype[clinical.by.data1$V1 == 0]

  ]

clinical.values <-

  clinical.data[

    enttype %in% clinical.by.data1$enttype[clinical.by.data1$V1 > 0]

    ]

### Clinical history

# get aggregate statistics for each medcode

clinical.history.by.medcode <-

  clinical.history[, length(unique(anonpatid)), by = medcode]

names(clinical.history.by.medcode) <- c('medcode', 'n.pat')

# Print a table of the leading medcodes we've found

print(

  merge(

    clinical.history.by.medcode[order(n.pat, decreasing = TRUE)[1:100],],

    CALIBER_DICT[, c('medcode', 'term')], by = 'medcode'

  )

)

# Take the top n by number of patients with that code

top.history <-

  clinical.history.by.medcode$medcode[

    order(

      clinical.history.by.medcode$n.pat, decreasing = TRUE

    )[1:n.top.clinical.history]

  ]

# And now, let's put how far in the past the patient was first diagnosed with

# each of these things into the table...

# First, discard all the rows corresponding to non-top OPCS codes

clinical.history <- clinical.history[clinical.history$medcode %in% top.history, ]

# Then, keep only the earliest instance of each per patient, because those are

# the ones we care about

clinical.history <-

  clinical.history[, min(relativedate), by = c('anonpatid', 'medcode')]

# Per medcode, add a new column to the cohort and put in numbers

clinical.history <-

  dcast(

    data = clinical.history,

    formula = anonpatid ~ medcode,

    value.var = "V1"

  )

# Add a prefix to the names to keep track

names(clinical.history)[2:(n.top.clinical.history + 1)] <-

  paste0('clinical.history.', names(clinical.history)[2:(n.top.clinical.history + 1)])

# Now, merge with the original cohort

COHORT <-

  merge(

    COHORT,

    clinical.history,

    by = c('anonpatid'),

    all = TRUE

  )

### Clinical values

# Next, let's melt the values data by data1, data2 etc. Each potentially

# contains a separate measurement (eg for entity type 13, which is weight, data1

# is the weight, data2 is the weight centile [always blank in this dataset!] and

# data3 is BMI)

clinical.values <-

  melt(

    clinical.values,

    id.vars = c("anonpatid", "relativedate", "medcode", "enttype"),

    measure.vars = paste0("data", 1:7)

  )

# Because there are lots of data types, the value column gets coerced to double.

# This is fine for now, because it's all numeric, and whilst some values are

# categorical represented as integers, random forests don't care about that.

# Remove all NA values

clinical.values <- clinical.values[!is.na(value)]

# Remove all values measured too long ago

clinical.values <- clinical.values[relativedate <= 183]

# aggregate by enttype and which data column the test came from

clinical.values.by.type <-

  clinical.values[, length(unique(anonpatid)), by = c('enttype', 'variable')]

names(clinical.values.by.type) <- c('enttype', 'variable', 'n.pat')

top.clinical.enttypes <-

  clinical.values.by.type$enttype[order(clinical.values.by.type$n.pat, decreasing = TRUE)[1:n.top.clinical.values]]

top.clinical.dataN <-

  clinical.values.by.type$variable[order(clinical.values.by.type$n.pat, decreasing = TRUE)[1:n.top.clinical.values]]

# Now, discard all the rows corresponding to non-top tests

clinical.values <-

  clinical.values[

    # Needs to be exact match of enttype and variable combination

    paste0(enttype, '!!!', variable) %in%

      paste0(top.clinical.enttypes, '!!!', top.clinical.dataN),

    ]

# Per patient and per test, keep the smallest relativedate value only

clinical.values.most.recent <-

  clinical.values[, .SD[which.min(relativedate)], by = c('anonpatid', 'enttype', 'variable')]

# Make a column per test

clinical.values.most.recent.wide <-

  dcast(

    data = clinical.values.most.recent,

    formula = anonpatid ~ enttype + variable,

    value.var = "value"

  )

# Add a prefix to the names to keep track

# (There may not be n.top.tests columns here as some get lost during the paring

# down processes above...)

names(clinical.values.most.recent.wide)[-1] <-

  paste0('clinical.values.', names(clinical.values.most.recent.wide)[-1])

# Now, merge with the original cohort

COHORT <-

  merge(

    COHORT,

    clinical.values.most.recent.wide,

    by = c('anonpatid'),

    all = TRUE

  )

### Therapy ####################################################################

# read in the therapy data

therapy.data <- readMedicalData(

  therapy.files,

  c("anonpatid", "eventdate", "bnfcode"),

  c("integer", "date", "integer")

)

# The other option than bnfcode is prodcode, which refers to specific products

# rather than BNF categories. There are far more of these so, assuming the BNF

# classification is somewhat rational, I'm going to go with that first to

# reduce data sparsity.

# Now, merge with the indexdate...we only want data before then, because looking

# after it is cheating, and using all data rather than just stuff before may

# distort our choice of variables as some variables may be very common after

# entering the study, but less so before... (14 differ...)

therapy.data <-

  merge(

    therapy.data, COHORT[, c('anonpatid', 'indexdate')],

    by = 'anonpatid', all.x = TRUE

  )

therapy.data$relativedate <- therapy.data$indexdate - therapy.data$eventdate

# Now, remove all the negative relative dates, because they're in the past

therapy.data <- therapy.data[relativedate >= 0]

# And remove all data which is too far into the past

therapy.data <- therapy.data[relativedate < 366]

# get aggregate statistics for each BNF code

therapy.by.bnf <- therapy.data[, length(unique(anonpatid)), by = bnfcode]

names(therapy.by.bnf) <- c('bnfcode', 'n.pat')

# Take the top n by number of patients with that code

top.bnf <-

  therapy.by.bnf$bnfcode[order(therapy.by.bnf$n.pat, decreasing = TRUE)[1:n.top.therapy]]

# Discard all the rows corresponding to non-top BNF codes

therapy.data <- therapy.data[therapy.data$bnfcode %in% top.bnf, ]

# Aggregate by number of prescriptions per patient

therapy.data <- therapy.data[, .N, by = list(anonpatid, bnfcode)]

# Per BNF code, add a new column to the cohort and put in numbers

therapy.wide <-

  dcast(

    data = therapy.data,

    formula = anonpatid ~ bnfcode,

    value.var = "N"

  )

# Add a prefix to the names to keep track

names(therapy.wide)[2:(n.top.therapy + 1)] <-

  paste0('bnf.', names(therapy.wide)[2:(n.top.therapy + 1)])

# And merge into the overall cohort

COHORT <-

  merge(

    COHORT,

    therapy.wide,

    by = c('anonpatid'),

    all = TRUE

  )

### Anonymisation steps ########################################################

# delete the columns with obvious privacy issues

COHORT[, (kill.cols) := NULL]

# make the remaining columns relative to the indexdate

indexdate <- as.Date(COHORT$indexdate)

# date.cols specified

for(date.col in c(date.cols)) {

    COHORT[[date.col]] <-

      # Do indexdate minus, so this is days in the past

      as.numeric(indexdate - as.Date(COHORT[[date.col]]))

    # ...and therefore negative values are in the future, hence cheating

    COHORT[[date.col]][COHORT[[date.col]] < 0] <- NA

}

# make age an integer

COHORT$age <- round(COHORT$age)

# make pracregion and ethnicity into categories

lookup_pracregion <-

  sample(unique(COHORT$pracregion),length(unique(COHORT$pracregion)))

lookup_ethnicity <-

  sample(unique(COHORT$hes_ethnicity),length(unique(COHORT$hes_ethnicity)))

write.csv(lookup_pracregion, 'lookup_pracregion.csv')

write.csv(lookup_ethnicity, 'lookup_ethnicity.csv')

COHORT$pracregion <-

  as.integer(factor(COHORT$pracregion, levels = lookup_pracregion))

COHORT$hes_ethnicity <-

  as.integer(factor(COHORT$hes_ethnicity, levels = lookup_ethnicity))

# make IMD score into deciles-ish

COHORT$imd_score <- round(COHORT$imd_score/10)

write.csv(COHORT, 'cohort-datadriven.csv')

fake.df <- data.frame(id = 1:10)

for (colname in names(COHORT)) {

    fake.df[,colname] <- sample(COHORT[!is.na(COHORT[, colname]), colname], 10)

}

write.csv(fake.df, 'cohort-sample.csv')