12_multiple_disease.Rmd

---
title: "Diffusion scores on several diseases"
author: "Sergio Picart-Armada"
date: "October 9, 2017"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE)
```

# Getting started

```{r}
# Data handling
library(plyr)
library(dplyr)
library(tidyr)
library(magrittr)

# ML libraries
library(kernlab)
library(caret)
library(mlr)
library(parallelMap)
parallelStartSocket(12)

library(igraph)
library(ggplot2)
library(ggsci)
library(diffuStats)

library(EGAD)
library(RANKS)
library(COSNet)

library(multcomp)

# have all config variables in a different env
config <- new.env(parent = globalenv())
source("03_config.R", local = config)

# load dataset and kernel
load(config$graph_alldisease)
load(config$file_kernel3)
# adjacency matrix with 1's in the 
A <- igraph::as_adj(g_filter, sparse = TRUE, attr = "weight") %>% as.matrix
# diag(A) <- 1

# Load MashUp features
nm_mashup <- readLines(config$file_mashup_names)
df_mashup <- data.table::fread(config$file_mashup_features) %>% 
  t %>% 
  as.data.frame %>% 
  set_rownames(nm_mashup)

# the dataset
# x: training, y: response
df_disease <- mutate(
  g_filter$dataset, 
  disease = as.factor(disease.efo_info.label)) %>% 
  rename(drugs = known_drug_binary, genetic = known_gene_binary) %>% 
  mutate(validation = drugs)
n <- vcount(g_filter)
```

# Reference streams

```{r}
# data frame with input scores
df_input <- plyr::ddply(
  reshape2::melt(df_disease, measure.vars = c("drugs", "genetic"), 
                 variable.name = "input_type", value.name = "input"),
  c("disease", "input_type"), 
  function(disease) {
    # browser()
    nm <- V(g_filter)$name
    x <- nm %in% (filter(disease, input == 1)$STRING_id)
    val <- nm %in% (filter(disease, validation == 1)$STRING_id)
    
    # genetic_original <- setNames(numeric(length(nm)), nm)
    # genetic_original[as.character(disease$STRING_id)] <- 
    #   disease$association_score.datatypes.genetic_association
    
    data.frame(
      STRING_id = nm, 
      score = x*1, 
      validation = val*1
      # genetic_original = genetic_original
    )
  }, 
  .progress = "text"
)

# Save the drugs data as a long data frame
df_drugs <- filter(df_input, input_type == "drugs") %>%
  dplyr::select(disease, STRING_id, validation) %>% 
  rename(drugs = validation)
# pairs disease-gene
write.csv(filter(df_drugs, drugs == 1L) %>% 
            dplyr::select(-drugs), 
          file = paste0(config$dir_data10, "/disease_drugs.csv"), 
          row.names = FALSE)
# genes in several diseases
df_ndisease <- plyr::ddply(
  df_drugs, "STRING_id", 
  function(df) c(n_disease = sum(df$drugs))) %>% 
  dplyr::filter(n_disease >= 10) %>%
  dplyr::arrange(desc(n_disease))
write.csv(df_ndisease, 
          file = paste0(config$dir_data10, "/disease_count_genes.csv"), 
          row.names = FALSE)

# Small tests
# a <- subset(df_input, input_type == "drugs")
# b <- subset(df_input, input_type == "genetic")
# all(a$score == a$validation)
# all(a$score == b$validation)

# select streams other than drug and overall
df_streams <- dplyr::select(
  df_disease, disease, STRING_id, 
  contains("association_score"), -contains("drug"), 
  -contains("overall")) %>% 
  gather(key = stream, value = score, 
         contains("association_score"), factor_key = TRUE) %>%
  ddply(
    c("disease", "stream"), 
    function(df_dis) {
      nm <- V(g_filter)$name
      x <- nm %in% df_dis$STRING_id
      
      original <- setNames(numeric(length(nm)), nm)
      original[as.character(df_dis$STRING_id)] <- df_dis$score
      
      data.frame(STRING_id = nm, score = original)
    }, 
    .progress = "text"
)
```

```{r}
df_descriptive <- plyr::ddply(
  dplyr::select(df_input, -validation) %>% 
    reshape2::dcast(STRING_id+disease~input_type, fun.aggregate = NULL, value.var = "score"), 
  "disease", 
  function(df) {
    tibble(
      n_genetic = sum(df$genetic), 
      n_drug = sum(df$drugs), 
      overlap = sum(df$genetic*df$drugs), 
      p_value = fisher.test(x = df$genetic, y = df$drugs)$p.value
    ) 
  }, 
  .id = "disease") %>% 
  mutate(fdr = p.adjust(p_value, method = "fdr"))

write.csv(df_descriptive, 
          file = paste0(config$dir_data10, "/descriptive_diseases.csv"), 
          row.names = FALSE)
df_descriptive

# check for NAs
stopifnot(all(!is.na(df_input)))

# group_by(df_input, disease, input_type) %>% 
#   select(score) %>% 
#   summarise_all(c("min", "mean", "max"))

# cross-validation parameters
k <- 3
times <- 25
```

# Performance measures

```{r}
# performance measures
top_k <- function(k) {
  function(actual, predicted) {
    inds <- head(order(predicted, decreasing = TRUE), k)
    sum(actual[inds])
  }
}

list_metrics <- list(
  auroc = metric_fun(curve = "ROC"), 
  partial_auroc_0.9 = metric_fun(curve = "ROC", c(0, .1)), 
  partial_auroc_0.95 = metric_fun(curve = "ROC", c(0, .05)), 
  auprc = metric_fun(curve = "PRC"),
  top_20_hits = top_k(20), 
  top_100_hits = top_k(100)
) 
```

# Defining functions for ML-based methods

```{r}
# C-SVM Bagging wrapper 
# Assumption: number of positives << number of negatives 
# Aggregation is at the level of "decision" values!
# Because if it is on the predicted class, it 
# does not work well, i.e. too many times the SVM 
# predicts only positives (or negatives), although the
# raking of the "decision" values can be meaningful
# (bootstrap is on the negatives!)
# ind_train, ind_test: numeric or character vectors with the 
# ids of the training and the testing samples
# ytrain binary vector with training labels (+:1, -:0)
# K graph kernel matrix
# B number of bootstrapping iterations
# ... further arguments for ksvm
bag_svm <- function(ind_train, ind_test, K, ytrain, B = 30, ...) {
  stopifnot(length(ind_train) == length(ytrain))
  
  # SVM bagging
  ind_pos <- which(ytrain == 1L)
  ind_neg <- which(ytrain == 0L)
  npos <- length(ind_pos)
  nneg <- length(ind_neg)
  
  yt <- as.factor(ytrain)
  mat_bag <- plyr::ldply(
    1:B, 
    function(rep) {
      # browser()
      ind_bag_neg <- sample(ind_neg, npos, replace = TRUE)
      ind_bag_all <- c(ind_pos, ind_bag_neg)
      # this one can contain names:
      ind_bag_orig <- ind_train[ind_bag_all]
      
      # train svms
      svm_mod <- ksvm(
        as.kernelMatrix(K[ind_bag_orig, ind_bag_orig]), 
        yt[ind_bag_all], 
        kernel = "matrix", 
        ...
      )
      # find support vectors
      svm_vec <- ind_bag_orig[SVindex(svm_mod)]
      
      # predict using precomputed kernel
      predict(svm_mod, 
              as.kernelMatrix(K[ind_test, svm_vec]), 
              type = "decision") %>% as.vector
    }, 
    .progress = "text"
    # this last line will work if ind_test are rownames and if they are 
    # numeric values
  ) %>% colMeans %>% setNames(rownames(K[ind_test, 1, drop = FALSE]))
}

# nu-SVM RBF and random forest wrapper, with data downsampling
# Assumption: number of positives << number of negatives 
# ind_train, ind_test: numeric or character vectors with the 
# ids of the training and the testing samples
# ytrain binary vector with training labels (+:1, -:0)
# df_features MashUp features as a data frame
mlr_svm_rf <- function(ind_train, ind_test, df_features, ytrain) {
  # levels are 0, 1 (in this order)
  yclass <- as.factor(ytrain)
  mean_y <- mean(ytrain)
  # how many more negatives are there?
  ratio <- (1 - mean_y)/mean_y
  
  # create task
  # shared by all learners
  # positives are coded as "1" and negatives as "0"
  tsk <- makeClassifTask(
    id = "diffu", 
    data = data.frame(class = yclass, df_features[ind_train, ]), 
    target = "class", 
    positive = "1")
  tsk <- undersample(tsk, 1/ratio)

  ###### SVM LEARNER ###### 
  num_ps <- makeParamSet(
    # makeNumericParam("C", lower = -5, upper = 5, trafo = function(x) 10^x),
    makeNumericParam("nu", lower = .1, upper = .9),
    makeNumericParam("sigma", lower = -6, upper = 2, trafo = function(x) 10^x)
  )
  ctrl <- makeTuneControlGrid(resolution = 5L, tune.threshold = FALSE)
  
  # define learner
  lrn <- makeLearner("classif.ksvm", predict.type = "prob")
  rdesc <- makeResampleDesc("CV", iters = 3L, stratify = TRUE)
  
  # Grid search in parallel
  res <- tuneParams(
    lrn, 
    task = tsk, 
    resampling = rdesc,
    par.set = num_ps, 
    measures = list(auc),
    control = ctrl)
  
  # Fit optimal params
  lrn.optim <- setHyperPars(lrn, par.vals = res$x)
  m <- train(lrn.optim, tsk)
  m
  
  # predict
  pred_svm <- predict(m, newdata = df_features[ind_test, ])
  
  ###### RandomForest LEARNER ###### 
  num_ps <- makeParamSet(
    makeIntegerParam("ntree", lower = 10, upper = 500), 
    # makeIntegerParam("mtry", lower = 10, upper = 50), 
    makeIntegerParam("nodesize", lower = 1, upper = 5)
  )
  ctrl <- makeTuneControlGrid(resolution = 3L, tune.threshold = TRUE)
  
  # define learner
  lrn <- makeLearner("classif.randomForest", predict.type = "prob")
  rdesc <- makeResampleDesc("CV", iters = 3L, stratify = TRUE)
  
  # Grid search in parallel
  res <- tuneParams(
    lrn, 
    task = tsk, 
    resampling = rdesc,
    par.set = num_ps, 
    measures = list(auc),
    control = ctrl)
  
  # Fit optimal params
  lrn.optim <- setHyperPars(lrn, par.vals = res$x)
  m <- train(lrn.optim, tsk)
  m
  
  # Test set
  pred_rf <- predict(m, newdata = df_features[ind_test, ])
  
  list(
    svm = setNames(pred_svm$data$prob.1, rownames(pred_svm$data)), 
    rf = setNames(pred_rf$data$prob.1, rownames(pred_rf$data))
  )
}
```


# Computing the diffusion scores through cross-validation

```{r}
# reproducibility
set.seed(1)

# A centrality measure
pr <- page.rank(g_filter)$vector

# keep appending on csv file
# write file and header (careful, overwrites prior file)
header_csv <- c("disease", "input_type", "split_cv", 
                "method", names(list_metrics))
write(header_csv, file = config$file_csv_perf12, 
      ncolumns = length(header_csv), sep = "\t")

# This took a bit less than 6 days
# The metrics are written in each iteration, providing 
# access to partial results while running
df_perf <- plyr::ddply(
  df_input, 
  # subset(df_input, input_type == "genetic"), 
  c("disease", "input_type"), 
  function(df_in) {
    # browser()
    name_disease <- as.character(df_in$disease[1])
    name_input_type <- as.character(df_in$input_type[1])
    
    df_streams_disease <- filter(df_streams, disease == name_disease)
    
    # Stratified split
    # x: input
    # y: validation
    # Both are named vectors
    x <- setNames(df_in$score, df_in$STRING_id)
    y <- setNames(df_in$validation, df_in$STRING_id)
    
    # split the dataset, stratified CV on the validation labels
    # this returns the index of the training instances
    list_split_cv <- caret::createMultiFolds(y = y, k = k, times = times)
    
    # diffusion scores
    list_perf <- plyr::ldply(
      setNames(names(list_split_cv), names(list_split_cv)), 
      function(split_cv_name) {
        # browser()
        ######## define training and validation ########
        split_cv_train <- list_split_cv[[split_cv_name]]
        # train vectors, with three formats
        # diffuStats + EGAD: positive 1, negative 0, unabelled NULL
        vec_diffustats <- x[split_cv_train]
        # COSnet: positive 1, negative -1, unlabelled 0
        vec_cosnet <- ifelse(x == 1, 1, -1)
        vec_cosnet[-split_cv_train] <- 0
        # RANKS: which(1) - but only for training fold!
        vec_ranks <- which(vec_cosnet == 1)
        # EGAD: positive 1, negative/unlabelled 0
        vec_egad <- ifelse(vec_cosnet == 1, 1, 0)
        # debug
        # table(vec_diffustats)
        # table(vec_cosnet)
        # length(vec_ranks)
        # table(vec_egad)
        
        # validation labels
        vec_val <- y[-split_cv_train]
        # for safety, we will index all the results using the 
        # names of the validation genes, making sure the order 
        # is kept...
        names_train <- names(vec_diffustats)
        names_val <- names(vec_val)

        ######## diffusion-based approaches ########
        # diffuStats
        list_scores <- plyr::llply(
          setNames(config$list_methods, config$list_methods),
          function(method) {
            diffuStats::diffuse(
              K = K, scores = vec_diffustats, 
              method = method, n.perm = 1e3)[names_val] 
          }
        )
        
        # personalized PageRank
        list_scores$ppr <- page.rank(
          g_filter, personalized = vec_egad)$vector[names_val]
        
        # EGAD (gba)
        list_scores$EGAD <- EGAD::predictions(
          genes.labels = vec_egad, 
          network = A
        )[, 1]
        list_scores$EGAD <- list_scores$EGAD[names_val]
        
        # RANKS (wsld + knn) kernelized scores
        list_scores$wsld <- RANKS::WSLD.score(
          RW = K, x = 1:nrow(K), x.pos = vec_ranks, d = config$wsld_d) %>% 
          setNames(rownames(K))
        list_scores$wsld <- list_scores$wsld[names_val]
        list_scores$knn <- RANKS::KNN.score(
          RW = K, x = 1:nrow(K), x.pos = vec_ranks, k = config$knn_k) %>% 
          setNames(rownames(K))
        list_scores$knn <- list_scores$knn[names_val]
        
        ######## other machine learning approaches ########
        # based in prodige1: SVM
        list_scores$bagsvm <- bag_svm(
          ind_train = names_train, ind_test = names_val, 
          K = K, ytrain = vec_diffustats, B = 30, 
          C = 1, type = "C-svc", scaled = FALSE)
        # based in MashUp: SVM and RandomForest
        list_mashup <- mlr_svm_rf(
          ind_train = names_train, ind_test = names_val, 
          df_features = df_mashup, ytrain = vec_diffustats)
        list_scores <- c(list_scores, list_mashup)
        
        # COSNet: neural net
        list_scores$COSNet <- COSNet::COSNet(
          W = A, labeling = vec_cosnet, cost = config$cosnet_cost
        )$scores
        list_scores$COSNet <- list_scores$COSNet[names_val]
        
        ######## reference scores ########
        # random
        list_scores$random <- setNames(
          sample(length(names_val)), names_val
        )
        # network properties that ignore input
        list_scores$randomraw <- diffuStats::diffuse(
          K = K, 
          scores = setNames(sample(vec_diffustats), 
                            names(vec_diffustats)), 
          method = "raw")[names_val] 
        list_scores$pr <- pr[names_val]
        
        # reference: other data streams
        list_streams <- plyr::dlply(
          df_streams_disease, 
          "stream", 
          function(df) {
            stream_scores <- setNames(df$score, df$STRING_id)
            stream_scores[names_val]
          }
        )
        
        ######## performance metrics ########
        # compute metrics 
        df_metrics <- plyr::ldply(
          c(list_scores, list_streams), 
          function(scores) {
            perf_eval(
              prediction = scores, 
              validation = vec_val, 
              metric = list_metrics
            )
          }, 
          .id = "method"
        ) 
        
        # Append to file, so that we can keep track on how the file grows 
        # and how advanced the process is
        df_append <- data.frame(
          disease = name_disease, 
          input_type = name_input_type, 
          split_cv = split_cv_name, 
          df_metrics
        )
        write.table(df_append, config$file_csv_perf12, sep = "\t",
                    append = TRUE, row.names = FALSE, col.names = FALSE)
        
        # return the metrics
        df_metrics
      }, 
      .id = "split_cv", 
      .progress = "text"
    )
  }
)
```

# Prediction on the whole dataset

```{r, eval=FALSE}
# We have estimated the performance already.
# This gives the actual predictions using current disease data
# Now this predicts based on drugs and genetic scores
df_predict <- plyr::ddply(
  df_input,
  # subset(df_input, input_type == "drugs"),
  c("input_type", "disease"), 
  function(df_in) {
    set.seed(1)
    # browser()

    # x: input
    x <- setNames(df_in$score, df_in$STRING_id)
    
    # browser()
    ######## define input vector in all formats ########
    # diffuStats + EGAD: positive 1, negative 0, unabelled NULL
    vec_diffustats <- x
    # COSnet: positive 1, negative -1, unlabelled 0
    # THIS IS DIFFERENT THAN IN THE CV!
    vec_cosnet <- ifelse(x == 1, 1, 0)
    # RANKS: which(1) - but only for training fold!
    vec_ranks <- which(vec_cosnet == 1)
    # EGAD: positive 1, negative/unlabelled 0
    vec_egad <- ifelse(vec_cosnet == 1, 1, 0)
    # debug
    # table(vec_diffustats)
    # table(vec_cosnet)
    # length(vec_ranks)
    # table(vec_egad)
    
    # for safety, we will index all the results using the 
    # names...
    names_train <- names(vec_diffustats)

    ######## diffusion-based approaches ########
    # diffuStats
    list_scores <- plyr::llply(
      setNames(config$list_methods, config$list_methods),
      function(method) {
        diffuStats::diffuse(
          K = K, scores = vec_diffustats, 
          method = method, n.perm = 1e3)[names_train] 
      }
    )
    
    # personalized PageRank
    list_scores$ppr <- page.rank(
      g_filter, personalized = vec_egad)$vector[names_train]
    
    # EGAD (gba)
    list_scores$EGAD <- EGAD::predictions(
      genes.labels = vec_egad, 
      network = A
    )[, 1]
    list_scores$EGAD <- list_scores$EGAD[names_train]
    
    # RANKS (wsld + knn) kernelized scores
    list_scores$wsld <- RANKS::WSLD.score(
      RW = K, x = 1:nrow(K), x.pos = vec_ranks, d = config$wsld_d) %>% 
      setNames(rownames(K))
    list_scores$wsld <- list_scores$wsld[names_train]
    list_scores$knn <- RANKS::KNN.score(
      RW = K, x = 1:nrow(K), x.pos = vec_ranks, k = config$knn_k) %>% 
      setNames(rownames(K))
    list_scores$knn <- list_scores$knn[names_train]
    
    ######## other machine learning approaches ########
    # based in prodige1: SVM
    list_scores$bagsvm <- bag_svm(
      ind_train = names_train, ind_test = names_train, 
      K = K, ytrain = vec_diffustats, B = 30, 
      C = 1, type = "C-svc", scaled = FALSE)
    # based in MashUp: SVM and RandomForest
    list_mashup <- mlr_svm_rf(
      ind_train = names_train, ind_test = names_train, 
      df_features = df_mashup, ytrain = vec_diffustats)
    list_scores <- c(list_scores, list_mashup)
    
    # COSNet: neural net
    scores_COSNet <- COSNet::COSNet(
      W = A, labeling = vec_cosnet, cost = config$cosnet_cost
    )$scores
    list_scores$COSNet <- ifelse(
      names_train %in% names(scores_COSNet), 
      scores_COSNet[names_train], 
      Inf
    ) %>% setNames(names_train)
    
    # return as a long data frame
    plyr::ldply(
      list_scores, 
      function(x) {
        data.frame(STRING_id = names(x), score = x, stringsAsFactors = TRUE)
      }, 
      .id = "method")
  }, 
  .progress = "text"
)
save(df_predict, 
     file = paste0(config$dir_performance12, "/all_predictions.RData"), 
     compress = "xz")
```


# Plotting performance metrics

```{r}
# column for repetition only
df_perf$rep_cv <- gsub("(Fold\\d+\\.)(Rep.+)", "\\2", df_perf$split_cv)
dim(df_perf)

df_perf_means <- df_perf %>% 
  select(-split_cv) %>%
  group_by(rep_cv, method, input_type, disease) %>%
  summarise_all(funs(mean))
  
dim(df_perf_means)

# obtain data frame to plot
df_plot <- reshape2::melt(
  df_perf_means, 
  id.vars = c("rep_cv", "method", "input_type", "disease"))

levels_plot <- c(
  "pr", 
  "randomraw", 
  "random", 
  "EGAD", 
  "association_score.datatypes.affected_pathway", 
  "association_score.datatypes.animal_model", 
  "association_score.datatypes.genetic_association", 
  "association_score.datatypes.literature", 
  "association_score.datatypes.rna_expression", 
  "association_score.datatypes.somatic_mutation", 
  "ppr", 
  "raw", 
  "gm", 
  "mc", 
  "z", 
  "knn", 
  "wsld", 
  "COSNet", 
  "bagsvm", 
  "rf", 
  "svm")
df_plot$method <- factor(df_plot$method, levels_plot)

# save variables for a posterior analysis
save(df_perf, df_perf_means, df_plot, 
     file = paste0(config$dir_performance12, "/all_diseases.RData"), 
     compress = "xz")
write.csv(
  x = df_perf_means, 
  file = paste0(config$dir_performance12, "/all_diseases.csv"), 
  row.names = FALSE)

```

```{r}
# plot the results
plyr::d_ply(
  df_plot, 
  "input_type", 
  function(df_type) {
    # browser()
    input_t <- df_type$input_type[1] %>% as.character
    gg <- ggplot(df_type, aes(x = method, y = value)) + 
      geom_boxplot(aes(fill = method), outlier.size = .3, lwd = .2) +
      # scale_fill_manual(guide = FALSE) +
      # scale_fill_brewer(palette = "Set3", guide = FALSE) +
      theme_bw() + 
      facet_grid(variable~disease, scales = "free_y") + 
      xlab("Method") + 
      ylab("Performance") + 
      ggtitle(paste0("Input: ", input_t, ", ", k, "-fold (repeated x", times, ") CV"), 
              subtitle = "Measures averaged per fold") + 
      theme(aspect.ratio = 1, 
            axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1, size = 6.5), 
            legend.position = "none")
    
    ggsave(filename = paste0(config$dir_performance12, "/boxplot_", input_t, ".png"), 
           device = NULL, plot = gg, width = 140, height = 50, units = "cm", 
           limitsize = FALSE)
  }
)
```

# Logistic models on AUROC and AUPRC

```{r}
old_maxpr <- getOption("max.print")
old_width <- getOption("width")

# WATCHOUT: is input_type genetic included or not?
list_mod_proportions <- plyr::dlply(
  filter(df_plot, grepl("au", variable) & input_type == "drugs"), 
  "variable", 
  function(df) {
    options(max.print = 2000, width = 1000)
    
    var <- df$variable[1]
    
    if(var == "partial_auroc_0.9") df$value <- df$value/.1
    if(var == "partial_auroc_0.95") df$value <- df$value/.05
    mod <- glm(value ~ method + disease, 
        data = df, 
        family = quasibinomial(link = "logit"))
    writeLines(
      capture.output(summary(mod)), 
      con = paste0(config$dir_performance12, "/logistic_", var, ".txt")
    )
    
    # compute contrasts between methods
    mod_glht <- multcomp::glht(
      mod, 
      linfct = multcomp::mcp(method = "Tukey"))
    mod_confint <- summary(mod_glht)
    writeLines(
      capture.output(print(mod_confint)), 
      con = paste0(config$dir_performance12, "/logistic_", var, "_contrasts.txt")
    )

    # Plot contrasts and estimates
    cols <- names(mod_confint$test$coefficients) %>% 
      strsplit(split = " - ") %>% plyr::ldply(
      function(x) c(row = x[1], col = x[2])
    ) 
    cols$row <- factor(cols$row, levels_plot)
    cols$col <- factor(cols$col, levels_plot)
    cols$estimate <- mod_confint$test$coefficients 
    cols$pvalue <- mod_confint$test$pvalues
    
    e.cor <- reshape2::acast(cols, row~col, value.var = "estimate", fill = 0)
    p.cor <- reshape2::acast(cols, row~col, value.var = "pvalue", fill = 0)
    
    png(filename = paste0(config$dir_performance12, "/logistic_", var, "_contrasts.png"), 
        width = 1200, height = 1200)
    corrplot::corrplot(
      corr = e.cor, method = "color", type = "lower", 
       addCoef.col = "black", 
      p.mat = p.cor, diag = TRUE, is.corr = FALSE)
    dev.off()
    
    list(model = mod, contrasts = mod_confint)
  }, 
  .progress = "text"
)
```

```{r}
list_mod_counts <- plyr::dlply(
  filter(df_plot, grepl("top", variable) & input_type == "drugs"), 
  "variable", 
  function(df) {
    options(max.print = 2000, width = 1000)
    var <- df$variable[1]
    
    # model
    mod <- glm(value ~ method + disease, 
        data = df, 
        family = quasipoisson(link = "log"))
    writeLines(
      capture.output(summary(mod)), 
      con = paste0(config$dir_performance12, "/poisson_", var, ".txt")
    )
    
    # compute contrasts between methods
    mod_glht <- multcomp::glht(
      mod, 
      linfct = multcomp::mcp(method = "Tukey"))
    mod_confint <- summary(mod_glht)
    writeLines(
      capture.output(print(mod_confint)), 
      con = paste0(config$dir_performance12, "/poisson_", var, "_contrasts.txt")
    )

    # Plot contrasts and estimates
    cols <- names(mod_confint$test$coefficients) %>% 
      strsplit(split = " - ") %>% plyr::ldply(
      function(x) c(row = x[1], col = x[2])
    ) 
    cols$row <- factor(cols$row, levels_plot)
    cols$col <- factor(cols$col, levels_plot)
    cols$estimate <- mod_confint$test$coefficients 
    cols$pvalue <- mod_confint$test$pvalues
    
    e.cor <- reshape2::acast(cols, row~col, value.var = "estimate", fill = 0)
    p.cor <- reshape2::acast(cols, row~col, value.var = "pvalue", fill = 0)
    
    png(filename = paste0(config$dir_performance12, "/poisson_", var, "_contrasts.png"), 
        width = 1200, height = 1200)
    corrplot::corrplot(
      corr = e.cor, method = "color", type = "lower", 
       addCoef.col = "black", 
      p.mat = p.cor, diag = TRUE, is.corr = FALSE)
    dev.off()
    
    list(model = mod, contrasts = mod_confint)
  }, 
  .progress = "text"
)
```

# Reproducibility

```{r}
out <- capture.output(sessionInfo())
writeLines(out, con = paste0(config$dir_metadata, "/12_sessionInfo_disease.txt"))
```