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


# Load the variable importances from the Cox model
cox.miss <-
  readRDS('../../output/caliber-replicate-with-missing-survreg-6-linear-age-surv-boot.rds')
cox.miss.vars <- bootStats(cox.miss, uncertainty = '95ci')
cox.miss.vars$var <- rownames(cox.miss.vars)
cox.miss.vars <- subset(cox.miss.vars, startsWith(var, 'vimp.c.index'))
cox.miss.vars$var <-
  substring(cox.miss.vars$var, nchar('vimp.c.index.') + 1)

# Load the variable importances from the random forest
rf.boot <- read.csv('../../output/rfsrc-cv-nsplit-try3-boot-all.csv')
rf.boot.vars <- bootStatsDf(rf.boot)
rf.boot.vars$var <- rownames(rf.boot.vars)
rf.boot.vars <- subset(rf.boot.vars, startsWith(var, 'vimp.c.index'))
rf.boot.vars$var <-
  substring(rf.boot.vars$var, nchar('vimp.c.index.') + 1)

var.imp.compare <- merge(cox.miss.vars, rf.boot.vars, by = c('var'))


# Plot a scatterplot of them
rf.vs.cox <-
  ggplot(
    var.imp.compare,
    aes(
      x = val.x, xmin = lower.x, xmax = upper.x,
      y = val.y, ymin = lower.y, ymax = upper.y
    )
  ) +
    geom_point() +
    geom_errorbar() +
    geom_errorbarh() +
    coord_cartesian(xlim = c(0, 0.03), ylim = c(0, 0.03))

print('Spearman correlation coefficient of variable importances:')
print(cor(var.imp.compare$val.x, var.imp.compare$val.y, method = 'spearman'))

# Load the variable importances from the big data model
cox.bigdata <- read.csv('../../output/cox-bigdata-varsellogrank-01-boot-all.csv')
cox.bigdata.vars <- bootStatsDf(cox.bigdata)
cox.bigdata.vars$var <- rownames(cox.bigdata.vars)
cox.bigdata.vars <- subset(cox.bigdata.vars, startsWith(var, 'vimp.c.index'))
cox.bigdata.vars$var <-
  substring(cox.bigdata.vars$var, nchar('vimp.c.index.') + 1)

cox.bigdata.vars <-
  cox.bigdata.vars[order(cox.bigdata.vars$val, decreasing = TRUE)[1:20], ]

cox.bigdata.vars <-
  cox.bigdata.vars[order(cox.bigdata.vars$val, decreasing = FALSE), ]

cox.bigdata.vars$description <- lookUpDescriptions(cox.bigdata.vars$var)

cat('c(', paste0("'", as.character(cox.bigdata.vars$description), "',"), ')', sep = '\n')

cox.bigdata.vars$description.manual <-
  factorOrderedLevels(
      c(
        'ALT',
        'PVD',
        'Hb',
        'Dementia',
        'Albumin',
        'Cardiac glycosides',
        'LV failure',
        'Home visit',
        'Oestrogens/HRT',
        'Chest pain',
        'Na',
        'WCC',
        'ALP',
        'Lymphocyte count',
        'Diabetes',
        'BMI ',
        'Weight',
        'Loop diuretics',
        'Smoking status',
        'Age'
      )
  )

# Plot a bar graph of them
cox.bigdata.plot <-
  ggplot(
    cox.bigdata.vars,
    aes(x = description.manual, y = val, ymin = lower, ymax = upper)
    ) +
    geom_bar(stat = 'identity') + 
    geom_errorbar(width = 0.25) +
    coord_flip() +
  theme(axis.title.y = element_blank(), axis.text.y = element_text(size = 10)) +
  ylim(0, 0.17)

# Random forest big data
rf.bigdata <- read.csv('../../output/rf-bigdata-varsellogrank-02-boot-all.csv')
rf.bigdata.vars <- bootStatsDf(rf.bigdata)
rf.bigdata.vars$var <- rownames(rf.bigdata.vars)
rf.bigdata.vars <- subset(rf.bigdata.vars, startsWith(var, 'vimp.c.index'))
rf.bigdata.vars$var <-
  substring(rf.bigdata.vars$var, nchar('vimp.c.index.') + 1)

rf.bigdata.vars <-
  rf.bigdata.vars[order(rf.bigdata.vars$val, decreasing = TRUE)[1:20], ]

rf.bigdata.vars <-
  rf.bigdata.vars[order(rf.bigdata.vars$val, decreasing = FALSE), ]

cat('c(', paste0("'", as.character(rf.bigdata.vars$description), "',"), ')', sep = '\n')

rf.bigdata.vars$description <- lookUpDescriptions(rf.bigdata.vars$var)

rf.bigdata.vars$description.manual <-
  factorOrderedLevels(
    c(
      'Fit note',
      'Stimulant laxatives',
      'Urea',
      'Hypertension',
      'Cardiac glycosides',
      'Beta2 agonists',
      'Telephone encounter',
      'Feet examination',
      'Creatinine',
      'Screening',
      'Osmotic laxatives',
      'ACE inhibitors',
      'Beta blockers',
      'Home visit',
      'Analgesics',
      'Blood pressure',
      'Chest pain',
      'Loop diuretics',
      'Smoking status',
      'Age'
    )
  )

# Plot a bar graph of them
rf.bigdata.plot <-
  ggplot(
    rf.bigdata.vars,
    aes(x = description.manual, y = val, ymin = lower, ymax = upper)
  ) +
    geom_bar(stat = 'identity') + 
    geom_errorbar(width = 0.25) +
    coord_flip() +
  theme(axis.title.y = element_blank(), axis.text.y = element_text(size = 10)) +
  ylim(0, 0.17)



# Elastic net Cox model
elastic.bigdata <- read.csv('../../output/cox-discrete-elasticnet-01-boot-all.csv')
elastic.bigdata.vars <- bootStatsDf(elastic.bigdata)
elastic.bigdata.vars$var <- rownames(elastic.bigdata.vars)
elastic.bigdata.vars <- subset(elastic.bigdata.vars, startsWith(var, 'vimp.c.index'))
elastic.bigdata.vars$var <-
  substring(elastic.bigdata.vars$var, nchar('vimp.c.index.') + 1)

elastic.bigdata.vars <-
  elastic.bigdata.vars[order(elastic.bigdata.vars$val, decreasing = TRUE)[1:20], ]

elastic.bigdata.vars <-
  elastic.bigdata.vars[order(elastic.bigdata.vars$val, decreasing = FALSE), ]

cat('c(', paste0("'", as.character(elastic.bigdata.vars$description), "',"), ')', sep = '\n')

elastic.bigdata.vars$description <- factorOrderedLevels(lookUpDescriptions(elastic.bigdata.vars$var))

elastic.bigdata.vars$description.manual <-
  factorOrderedLevels(
    c(
      'Biguanides',
      'CKD',
      'Osmotic laxatives',
      'MCV',
      'IMD score',
      'Dementia',
      'Home visit',
      'Sulphonylureas',
      'Insulin',
      'LV failure',
      'Cardiac glycosides',
      'Telephone encounter',
      'Records held date',
      'Chest pain',
      'Diabetes',
      'Smoking status',
      'Loop diuretics',
      'Gender',
      'Type 2 diabetes',
      'Age'
      )
  )

# Plot a bar graph of them
elastic.bigdata.plot <-
  ggplot(
    elastic.bigdata.vars,
    aes(x = description.manual, y = val, ymin = lower, ymax = upper)
  ) +
  geom_bar(stat = 'identity') + 
  geom_errorbar(width = 0.25) +
  coord_flip() +
  theme(axis.title.y = element_blank(), axis.text.y = element_text(size = 10)) +
  ylim(0, 0.17)


# Combine for output
plot_grid(
  rf.bigdata.plot, cox.bigdata.plot, elastic.bigdata.plot,
  labels = c('A', 'B', 'C'),
  ncol = 3
)

ggsave(
  '../../output/variable-importances.pdf',
  width = 24,
  height = 8,
  units = 'cm',
  useDingbats = FALSE
)


# Print a helpful list of overlapping predictors
print('RF vs Cox')
print(
  paste(
    rf.bigdata.vars$description.manual[
      rf.bigdata.vars$description.manual %in%
        cox.bigdata.vars$description.manual
      ], collapse = ', ')
)
print('Cox vs elastic net')
print(
  paste(
    elastic.bigdata.vars$description.manual[
      elastic.bigdata.vars$description.manual %in%
        cox.bigdata.vars$description.manual
      ], collapse = ', ')
)
print('RF vs elastic net not Cox')
print(
  paste(
    elastic.bigdata.vars$description.manual[
      elastic.bigdata.vars$description.manual %in%
        rf.bigdata.vars$description.manual & !(
          elastic.bigdata.vars$description.manual %in%
          cox.bigdata.vars$description.manual)
      ], collapse = ', ')
)
print('Cox vs elastic net not RF')
print(
  paste(
    elastic.bigdata.vars$description.manual[
      elastic.bigdata.vars$description.manual %in%
        cox.bigdata.vars$description.manual & !(
          elastic.bigdata.vars$description.manual %in%
            rf.bigdata.vars$description.manual)
      ], collapse = ', ')
)
# Should be none or the graph will be challenging to draw!

# Aborted idea to draw a rank-change chart which is too messy to be useful...
elastic.bigdata.vars$model <- 'enet'
rf.bigdata.vars$model <- 'rf'
cox.bigdata.vars$model <- 'cox'

all.models <- rbind(cox.bigdata.vars, rf.bigdata.vars, elastic.bigdata.vars)

all.models$model <- factor(all.models$model, levels = c('rf', 'cox', 'enet'))

ggplot(
  all.models,
  aes(
    x = model, y = log(val), ymin = log(lower), ymax = log(upper),
    label = description.manual, group = description.manual
  )) +
  geom_text(position = position_dodge(width = 0.7)) +
  geom_errorbar(width = 0.1, position = position_dodge(width = 0.7)) +
  geom_point(position = position_dodge(width = 0.7))