Atopo Performance vs Fitness Model Analysis

Intro

The purpose of this analysis is to find a correlation between the boolean models fitness to a steady state activity profile and their performance in terms of the number of True Positive (TP) synergies predicted and/or the overall MCC score (Matthews Correlation Coefficient score). We want to show that a closer fitness to the steady state suggests more predictive models, corroborating thus our proof of concept of using an ensemble-based approach where models are trained towards a specific steady state signaling pattern for drug combination predictions.

The boolean model datasets we will use are in total \(9\): one for each cell line of interest (8 cell lines) where the models were fitted to a specific steady state in each case and one for the so-called random models which were generated randomly in the sense that were fitted only to a proliferation state (simulations were done using the DrugLogics software modules Gitsbe and Drabme).

Each boolean model dataset constitues of:

• The model predictions file which has for each model the prediction for each drug combination tested (0 = no synergy predicted, 1 = synergy predicted, NA = couldn’t find stable states in either the drug combination inhibited model or in any of the two single-drug inhibited models)
• The models stable state (one per model). A fitness score for each model can easily be calculated then by matching the model’s stable state (which is something inherent in the boolean’s model structure, a unique fixpoint attractor) with the steady state of interest, node per node. A higher fitness score would mean a better match of a model’s stable state to the cell line derived steady state (a perfect match would result in a fitness of 1).
• The models link operators which is a representation of the boolean equations of each model. Each boolean equation is in the form: **Target *= (Activator OR Activator OR…) AND NOT (Inhibitor OR Inhibitor OR…)** and the difference between the models can be found in the link operator (1 = ‘OR NOT’, 0 = ‘AND NOT’, or absent) which has been mutated (changed) through the genetic algorithm in Gitsbe. Note that the equations that do not have link operators are the same for every model and are thus discarded.
• The observed synergies file which lists the drug combinations that were observed as synergistic for each cell line.
• The steady state file which lists the network nodes (protein, gene, complexes names, etc) and their activity value (0 or 1, representing an inhibited or active node respectively). This input is provided per cell line and not for the random models since they are just trained to a profileration state.

Input

Loading libraries:

library(DT)
library(ggpubr)
library(emba)
library(usefun)
library(nnet)
library(pscl)
library(ComplexHeatmap)
library(circlize)
library(dplyr)
library(tibble)
library(Ckmeans.1d.dp)
library(RColorBrewer)

First we load the cell-specific input data:

# Cell Lines
cell.lines = c("A498", "AGS", "DU145", "colo205", "SW620", "SF295", "UACC62", "MDA-MB-468")

cell.line.dirs = sapply(cell.lines, function(cell.line) {
  paste0(getwd(), "/", cell.line)
})

# Model predictions
model.predictions.files = sapply(cell.line.dirs, function(cell.line.dir) {
  paste0(cell.line.dir, "/model_predictions")
})

model.predictions.per.cell.line = lapply(model.predictions.files, 
  function(file) {
    get_model_predictions(file)
  }
)

# Observed synergies
observed.synergies.files = sapply(cell.line.dirs, function(cell.line.dir) {
  paste0(cell.line.dir, "/observed_synergies")
})

observed.synergies.per.cell.line = lapply(observed.synergies.files, 
  function(file) {
    get_observed_synergies(file)
  }
)

# Models Stable State (1 per model)
models.stable.state.files = sapply(cell.line.dirs, function(cell.line.dir) {
  paste0(cell.line.dir, "/models_stable_state")
})

models.stable.state.per.cell.line = lapply(models.stable.state.files,
  function(file) {
    as.matrix(read.table(file, check.names = FALSE))
  }
)

# Models Link Operators
models.link.operator.files = sapply(cell.line.dirs, function(cell.line.dir) {
  paste0(cell.line.dir, "/models_equations")
})

models.link.operators.per.cell.line = lapply(models.link.operator.files,
  function(file) {
    as.matrix(read.table(file, check.names = FALSE))
  }
)

# the node names used in our analysis
node.names = colnames(models.stable.state.per.cell.line[[1]])

# Steady States
steady.state.files = sapply(cell.line.dirs, function(cell.line.dir) {
  paste0(cell.line.dir, "/steady_state")
})

steady.state.per.cell.line = lapply(steady.state.files,
  function(file) {
    ss.df = read.table(file, sep = "\t", stringsAsFactors = FALSE)
    steady.state = ss.df[,2]
    names(steady.state) = ss.df[,1]
    
    # change value to NA for nodes for which there was no activity found (dash)
    steady.state[steady.state == "-"] = NA
    
    # keep only the nodes that are included in the analysis
    steady.state = prune_and_reorder_vector(steady.state, node.names)
    
    # return an integer vector since the activity values are binarized (0,1)
    return(sapply(steady.state, as.integer))
  }
)

The random model input data:

random.dir = paste0(getwd(), "/random")
random.model.predictions = get_model_predictions(paste0(random.dir, "/model_predictions"))

random.models.stable.state = as.matrix(
  read.table(file = paste0(random.dir, "/models_stable_state"), check.names = FALSE)
)

random.models.link.operator =
  as.matrix(read.table(file = paste0(random.dir, "/models_equations"), check.names = FALSE))

Model Analysis

In order to find the number of true positive (TP) predicted synergies, MCC scores and fitness scores for each of the models in each of the 9 datasets, we use functions from the emba R package.

Cell-specific

We find the MCC, TP and fitness values for each model per cell line (note that each model’s stable state in a specific cell line is matched against the steady state from that cell line):

models.tp.per.cell.line = list()
models.mcc.per.cell.line = list()
models.fitness.per.cell.line = list()

for (cell.line in cell.lines) {
  model.predictions = model.predictions.per.cell.line[[cell.line]]
  observed.synergies = observed.synergies.per.cell.line[[cell.line]]
  number.of.drug.comb.tested = ncol(model.predictions.per.cell.line[[cell.line]])
  
  # Split model.predictions to positive (observed) and negative (non-observed) results
  observed.model.predictions =
    get_observed_model_predictions(model.predictions, observed.synergies)
  unobserved.model.predictions =
    get_unobserved_model_predictions(model.predictions, observed.synergies)
  
  # Count the predictions of the observed synergies per model (TP)
  models.tp.per.cell.line[[cell.line]] = calculate_models_synergies_tp(observed.model.predictions)
  
  # Calculate Matthews Correlation Coefficient (MCC) for every model
  models.mcc.per.cell.line[[cell.line]] = 
    calculate_models_mcc(observed.model.predictions,
                         unobserved.model.predictions,
                         number.of.drug.comb.tested)
  
  # Fitness per model contrasted to steady state from cell line
  models.fitness.per.cell.line[[cell.line]] = 
    apply(models.stable.state.per.cell.line[[cell.line]], 1, get_percentage_of_matches, 
          steady.state.per.cell.line[[cell.line]])
}

Random Models

Next, we find the MCC, TP and fitness values for each random model per cell line (note that each random model’s stable state in a specific cell line is matched against the steady state from that cell line and that the random models’ stable state data does not change per cell line, i.e. same random.models.stable.state object):

random.models.mcc.per.cell.line = list()
random.models.tp.per.cell.line = list()
random.models.fitness.per.cell.line = list()

for (cell.line in cell.lines) {
  observed.synergies = observed.synergies.per.cell.line[[cell.line]]
  number.of.drug.comb.tested = ncol(random.model.predictions)
  
  # Split model.predictions to positive (observed) and negative (non-observed) results
  observed.model.predictions =
    get_observed_model_predictions(random.model.predictions, observed.synergies)
  unobserved.model.predictions =
    get_unobserved_model_predictions(random.model.predictions, observed.synergies)
  
   # Count the predictions of the observed synergies per model (TP)
  random.models.tp.per.cell.line[[cell.line]] = 
    calculate_models_synergies_tp(observed.model.predictions)
  
  # Calculate Matthews Correlation Coefficient (MCC) for every model
  random.models.mcc.per.cell.line[[cell.line]] = 
    calculate_models_mcc(observed.model.predictions,
                         unobserved.model.predictions,
                         number.of.drug.comb.tested)
  
  # Fitness per model contrasted to steady state from cell line
  random.models.fitness.per.cell.line[[cell.line]] = 
    apply(random.models.stable.state, 1, get_percentage_of_matches, 
          steady.state.per.cell.line[[cell.line]])
}

Choose best dataset

We now want to find the best dataset & cell line for our subsequent analysis - that is to show the performance vs fitness correlation. The argument here is that we want to choose a boolean model dataset that has a large enough fitness value range combined with a large TP and/or MCC value range, since with smaller value ranges it would be harder to distinguish the difference of the estimated distributions of the fitness scores belonging to different performance classes (i.e. models’ fitnesses that belong to different classification groups with the metric being either the number of TPs or the MCC score).

The next summary statistics tables and box-plots will help us determine exactly which cell line and dataset to use:

cell.specific.model.data = matrix(data = NA, nrow = length(cell.lines), ncol = 11)
rownames(cell.specific.model.data) = cell.lines
colnames(cell.specific.model.data) = c("fitness range", "Min fitness", 
  "Max fitness", "Mean fitness", "Median fitness", "MCC range", "Min MCC", 
  "Max MCC", "Mean MCC", "Median MCC", "Max TPR")

for (cell.line in cell.lines) {
  models.fitness = models.fitness.per.cell.line[[cell.line]]
  models.mcc = models.mcc.per.cell.line[[cell.line]]
  models.tp = models.tp.per.cell.line[[cell.line]]
  
  fit.summary = unclass(summary(models.fitness))
  mcc.summary = unclass(summary(models.mcc))
  max.tpr = max(models.tp) / length(observed.synergies.per.cell.line[[cell.line]])
  
  cell.specific.model.data[cell.line, "fitness range"] = fit.summary["Max."] - fit.summary["Min."]
  cell.specific.model.data[cell.line, "Min fitness"] = fit.summary["Min."]
  cell.specific.model.data[cell.line, "Max fitness"] = fit.summary["Max."]
  cell.specific.model.data[cell.line, "Mean fitness"] = fit.summary["Mean"]
  cell.specific.model.data[cell.line, "Median fitness"] = fit.summary["Median"]
  
  cell.specific.model.data[cell.line, "MCC range"] = mcc.summary["Max."] - mcc.summary["Min."]
  cell.specific.model.data[cell.line, "Min MCC"] = mcc.summary["Min."]
  cell.specific.model.data[cell.line, "Max MCC"] = mcc.summary["Max."]
  cell.specific.model.data[cell.line, "Mean MCC"] = mcc.summary["Mean"]
  cell.specific.model.data[cell.line, "Median MCC"] = mcc.summary["Median"]
  
  cell.specific.model.data[cell.line, "Max TPR"] = max.tpr
}

# color columns
fit.breaks = quantile(cell.specific.model.data[,"fitness range"], probs = seq(.05, .95, .05), na.rm = TRUE)
fit.colors = round(seq(255, 40, length.out = length(fit.breaks) + 1), 0) %>%
  {paste0("rgb(255,", ., ",", ., ")")} # red
mcc.breaks = quantile(cell.specific.model.data[,"MCC range"], probs = seq(.05, .95, .05), na.rm = TRUE)
mcc.colors = round(seq(255, 40, length.out = length(mcc.breaks) + 1), 0) %>%
  {paste0("rgb(", ., ",255,", ., ")")} # green

caption.title = "Table 1: Fitness, MCC scores and TP (True positives) for the Cell-specific model predictions across 8 Cell Lines"
datatable(data = cell.specific.model.data, 
          options = list(dom = "t"), # just show the table
          caption = htmltools::tags$caption(caption.title, style="color:#dd4814; font-size: 18px")) %>% 
  formatRound(1:11, digits = 3) %>%
  formatStyle(columns = c("fitness range"), backgroundColor = styleInterval(fit.breaks, fit.colors)) %>%
  formatStyle(columns = c("MCC range"), backgroundColor = styleInterval(mcc.breaks, mcc.colors))

random.model.data = matrix(data = NA, nrow = length(cell.lines), ncol = 11)
rownames(random.model.data) = cell.lines
colnames(random.model.data) = c("fitness range", "Min fitness", 
  "Max fitness", "Mean fitness", "Median fitness", "MCC range", "Min MCC", 
  "Max MCC", "Mean MCC", "Median MCC", "Max TPR")

for (cell.line in cell.lines) {
  models.fitness = random.models.fitness.per.cell.line[[cell.line]]
  models.mcc = random.models.mcc.per.cell.line[[cell.line]]
  models.tp = random.models.tp.per.cell.line[[cell.line]]
  
  fit.summary = unclass(summary(models.fitness))
  mcc.summary = unclass(summary(models.mcc))
  max.tpr = max(models.tp) / length(observed.synergies.per.cell.line[[cell.line]])
  
  random.model.data[cell.line, "fitness range"] = fit.summary["Max."] - fit.summary["Min."]
  random.model.data[cell.line, "Min fitness"] = fit.summary["Min."]
  random.model.data[cell.line, "Max fitness"] = fit.summary["Max."]
  random.model.data[cell.line, "Mean fitness"] = fit.summary["Mean"]
  random.model.data[cell.line, "Median fitness"] = fit.summary["Median"]
  
  random.model.data[cell.line, "MCC range"] = mcc.summary["Max."] - mcc.summary["Min."]
  random.model.data[cell.line, "Min MCC"] = mcc.summary["Min."]
  random.model.data[cell.line, "Max MCC"] = mcc.summary["Max."]
  random.model.data[cell.line, "Mean MCC"] = mcc.summary["Mean"]
  random.model.data[cell.line, "Median MCC"] = mcc.summary["Median"]
  
  random.model.data[cell.line, "Max TPR"] = max.tpr
}

# color columns
fit.breaks = quantile(random.model.data[,"fitness range"], probs = seq(.05, .95, .05), na.rm = TRUE)
fit.colors = round(seq(255, 40, length.out = length(fit.breaks) + 1), 0) %>%
  {paste0("rgb(255,", ., ",", ., ")")} # red
mcc.breaks = quantile(random.model.data[,"MCC range"], probs = seq(.05, .95, .05), na.rm = TRUE)
mcc.colors = round(seq(255, 40, length.out = length(mcc.breaks) + 1), 0) %>%
  {paste0("rgb(", ., ",255,", ., ")")} # green

caption.title = "Table 2: Fitness, MCC scores and TP (True positives) for the random model predictions across 8 Cell Lines"
datatable(data = random.model.data, 
          options = list(dom = "t"), # just show the table
          caption = htmltools::tags$caption(caption.title, style="color:#dd4814; font-size: 18px")) %>% 
  formatRound(1:11, digits = 3) %>%
  formatStyle(columns = c("fitness range"), backgroundColor = styleInterval(fit.breaks, fit.colors)) %>%
  formatStyle(columns = c("MCC range"), backgroundColor = styleInterval(mcc.breaks, mcc.colors))

The below box plots compare the MCC values and fitness scores across all cell lines between the cell-specific models (trained to steady state) and the random ones (trained to proliferation):

num.of.models = nrow(random.models.stable.state)
data.list = list()

for (cell.line in cell.lines) {
  cell.line.vec = as.data.frame(rep(cell.line, num.of.models), stringsAsFactors = FALSE)
  models.mcc.cell.specific = remove_rownames(as.data.frame(models.mcc.per.cell.line[[cell.line]]))
  models.mcc.random        = remove_rownames(as.data.frame(random.models.mcc.per.cell.line[[cell.line]]))
  
  models.fitness.cell.specific = remove_rownames(as.data.frame(models.fitness.per.cell.line[[cell.line]]))
  models.fitness.random        = remove_rownames(as.data.frame(random.models.fitness.per.cell.line[[cell.line]]))
  
  data.list[[cell.line]] = bind_cols(cell.line.vec, models.mcc.cell.specific, 
    models.mcc.random, models.fitness.cell.specific, models.fitness.random)
}

data = bind_rows(data.list)
colnames(data) = c("cell.line", "MCC cell-specific", "MCC random", 
                   "fitness cell-specific", "fitness random")

# Note the cell-specific models have NaN MCC values
ggboxplot(data, x = "cell.line", y = c("MCC cell-specific", "MCC random"),
          palette = brewer.pal(3, "Set1"), merge = "asis",
          xlab = "Cell Lines", ylab = "MCC values", add = "point", add.params = list(size = 0.5))

ggboxplot(data, x = "cell.line", y = c("fitness cell-specific", "fitness random"),
          palette = brewer.pal(3, "Set1"), merge = "asis",
          xlab = "Cell Lines", ylab = "Fitness values", add = "point", add.params = list(size = 0.5))

So, from the two tables and the two box plots above we conclude that:

• In general, the random models offer a larger range of fitness values when their stable states are matched against the steady state of each particular cell line. This is something we expected since these models weren’t fitted to a specific cell-line steady state (meaning that they weren’t chosen from the Gitsbe module as the 3 best from each simulation that match that steady state as best as possible) but rather to a more generic state of proliferation. Thus, they represent a set of models with larger variation in terms of structure (boolean model equations) compared to the cell-specific generated ones.
• The cell-specific models have always a larger maximum fitness and median value compared to the random ones for each respective cell line.
• In each cell line (with the expection of SW620) there are always cell-specific models that show better performance than the random ones (have a higher MCC value)

All in all, we will use the random models data, contrasted to the steady state of the A498 cell line (see Table 2). This dataset has the largest fitness and MCC value range (second largest TPR value as well) combined in both tables above.

best.cell.line = "A498"

fit = random.models.fitness.per.cell.line[[best.cell.line]]
tp  = random.models.tp.per.cell.line[[best.cell.line]]
mcc = random.models.mcc.per.cell.line[[best.cell.line]]

Statistical Analysis

Data preprocessing

Firstly, we filter the data by finding the unique models - those that have strictly different boolean equations. Then, we take a random sample out of these (while stabilizing the seed number for reproducibility purposes):

# For reproducibility
set.seed(0)
sample.size = 1000

unique.models = rownames(unique(random.models.link.operator))
unique.models.sample = sample(unique.models, size = sample.size)

fit.unique = fit[names(fit) %in% unique.models.sample]
tp.unique  = tp[names(tp) %in% unique.models.sample]
mcc.unique = mcc[names(mcc) %in% unique.models.sample]

df = as.data.frame(cbind(fit.unique, mcc.unique, tp.unique))

Note that the fitness and MCC score are continuous variables while the number of true positives (TP) is discrete.

Correlation Plots

Then, we check if the our data is normally distributed (using the Shapiro-Wilk test for normality and the Q-Q plots):

ggqqplot(data = df, x = "fit.unique", ylab = "Fitness")

ggqqplot(data = df, x = "mcc.unique", ylab = "MMC")

shapiro.test(x = sample(fit.unique))

    Shapiro-Wilk normality test

data:  sample(fit.unique)
W = 0.99266, p-value = 7.375e-05

shapiro.test(x = sample(mcc.unique))

    Shapiro-Wilk normality test

data:  sample(mcc.unique)
W = 0.91827, p-value < 2.2e-16

From the above results we observe that both the fitness and MCC scores are surely not normally distributed (with statistical significance). Thus, we will use non-parametric correlation scores, namely the Spearman and Kendall rank-based correlation tests, to check the correlation between the two continuous variables (models fitness values and their corresponding MCC score):

ggscatter(df, x = "mcc.unique", y = "fit.unique", 
          title = "Fitness vs Performance (MCC) - Spearman Correlation",
          add = "reg.line", conf.int = TRUE, 
          cor.coef = TRUE, cor.coeff.args = list(method = "spearman", label.x.npc = 0.7, label.y.npc = 1),
          xlab = "MCC scores", ylab = "Fitness Values")

ggscatter(df, x = "mcc.unique", y = "fit.unique", 
          title = "Fitness vs Performance (MCC) - Kendall Correlation",
          add = "reg.line", conf.int = TRUE, 
          cor.coef = TRUE, cor.coeff.args = list(method = "kendall", label.x.npc = 0.7, label.y.npc = 1),
          xlab = "MCC scores", ylab = "Fitness Values")

From the correlation plots above, we observe a weak/small positive correlation between performance and fitness to the steady state.

To assess the correlation between the number of TP predictions of the models (discrete variable) and the models fitness (continuous variable), we construct a predictor of the categorical variable from the continuous variable: if the resulting classifier has a high degree of fit we can conclude the two variables share a relationship and are indeed correlated. Since there are more than 2 TP classes (values) in the dataset, we will use Multinomial Logistic Regression and fit log-linear models via neural networks:

model.classifier = multinom(data = df, formula = tp.unique ~ fit.unique)

# weights:  18 (10 variable)
initial  value 1791.759469 
iter  10 value 1361.771381
iter  20 value 1356.825440
iter  30 value 1356.792829
iter  40 value 1356.765176
final  value 1356.764853 
converged

pseudo.r2.measures = pR2(model.classifier)

fitting null model for pseudo-r2
# weights:  12 (5 variable)
initial  value 1791.759469 
iter  10 value 1415.698062
final  value 1415.682090 
converged

pseudo.r2.measures["McFadden"]

  McFadden 
0.04161756 

To measure the goodness-of-fit for our model, there are several measures proposed for logistic regression. We emphasize on the McFaddens’s pseudo-\(R^2\) measure for our classifier, for which a value of \(0.2-0.4\) would indicate an excellent fit [source]. Since we found less, we can also assume that there is only a weak/small positive correlation between the number of TPs and the fitness score of the models.

Next, we will next proceed with a more elaborate analysis, where the models will be split to different performance classes (TP or MCC score-derived) and the statistical correlation between the individual groups will be tested (with regards to their fitness scores).

TP-class vs fitness

First, some box and density plots to see practically what and where is the difference between the models fitness values belonging to different TP-classes:

# Box plots
ggboxplot(df, x = "tp.unique", y = "fit.unique", color = "tp.unique",
          palette = usefun:::colors.100[1:length(unique(tp.unique))],
          xlab = "True Positives (TP)", ylab = "Fitness values")

# Density Plots
densities = list()
for (tp.num in sort(unique(tp.unique))) {
  x = df %>%
    filter(tp.unique == tp.num) %>%
    select_at(.vars = c("fit.unique"))
  den = density(x$fit.unique)
  densities[[paste0(tp.num, " (", den$n, ")")]] = den
}

make_multiple_density_plot(densities, legend.title = "TP classes (#models)",
        title = "Density Estimation", x.axis.label = "Fitness score")

As we can see from the above plots, there is positive correlation between the classes that predicted 0,1 and 2 TP synergies as well as between the 0 TP class and the 3,4,5-TP classes.

Mathematically, we show that the group distributions are indeed different using the Kruskal-Wallis Rank Sum Test:

# Hypothesis testing: Are the location parameters of the distribution of x the same in each group?
kruskal.test(x = df[,"fit.unique"], g = df[, "tp.unique"])

    Kruskal-Wallis rank sum test

data:  df[, "fit.unique"] and df[, "tp.unique"]
Kruskal-Wallis chi-squared = 108.06, df = 5, p-value < 2.2e-16

As the p-value is less than the significance level \(0.05\), we can conclude that there are significant differences between the different groups of fitness values. To see exactly which pair of groups are different we perform pairwise Wilcoxon rank sum tests and we draw a heatmap of the each test’s p-value:

res = pairwise.wilcox.test(x = df[,"fit.unique"], g = df[, "tp.unique"], p.adjust.method = "BH")
p.value.mat = res$p.value

col_fun = colorRamp2(breaks = c(0, 0.05, 0.5, 1), c("green", "white", "orange", "red"))
ht = Heatmap(matrix = p.value.mat, cluster_rows = FALSE, cluster_columns = FALSE,
  na_col = "white", name = "p-value", col = col_fun, column_names_rot = 0,
  row_title = "TP", row_title_rot = 0, row_names_side = "left",
  column_title = "P-values (Pairwise Wilcoxon tests) - TP-classified fitnesses", column_title_gp = gpar(fontsize = 20),
  cell_fun = function(j, i, x, y, width, height, fill) {
    if (!is.na(p.value.mat[i,j]))
        grid.text(sprintf("%.6f", p.value.mat[i, j]), x, y, gp = gpar(fontsize = 16))
})
draw(ht)

As we can see above there is significant difference between groups of fitness values belonging to models that predicted \(0\), \(1\) and \(2\) TP synergies, but not between these and the larger performant groups (\(TP=3,4,5\)) - with the exception of the group with \(TP=2\).

Note that the number of true positive predictions is an one-dimensional metric and by excluding the TN, FP and FNs we have a less-informant (and arguably incorrect) picture of the models performance classification. That’s why we will proceed with the more balanced MCC score for classifing the models fitnesses to different performance groups.

MCC-class vs fitness

Firstly, we perform a univariate k-means clustering to split the models MCC values to different classes and plot the data histogram. The number of MCC classes to split the data can be arbitrarily chosen and so we chose one less than the maximum number of TPs predicted:

num.of.mcc.classes = 5
mcc.class.ids = 1:num.of.mcc.classes

# find the clusters
res = Ckmeans.1d.dp(x = df[,"mcc.unique"], k = num.of.mcc.classes)
mcc.class.id = res$cluster
df = cbind(df, mcc.class.id)

plot_mcc_classes_hist(df[,"mcc.unique"], df[,"mcc.class.id"],
                      num.of.mcc.classes, mcc.class.ids)

Then we show some box and density plots to see practically what and where is the difference between the models fitness values belonging to different MCC-classes:

# Box plots
# if you want to add p-values on the boxplot
#mcc.class.id.cmps = list(c("1", "2"), c("1", "3"), c("3", "5"))
ggboxplot(df, x = "mcc.class.id", y = "fit.unique", color = "mcc.class.id",
          palette = usefun:::colors.100[1:num.of.mcc.classes],
          xlab = "MCC class", ylab = "Fitness values")

# + stat_compare_means(comparisons = mcc.class.id.cmps) + stat_compare_means()

# Density Plots
densities = list()
for (id in mcc.class.ids) {
  x = df %>%
    filter(mcc.class.id == id) %>%
    select_at(.vars = c("fit.unique"))
  den = density(x$fit.unique)
  densities[[paste0(id, " (", res$size[id], ")")]] = den
}

make_multiple_density_plot(densities, legend.title = "MCC classes (#models)",
        title = "Density Estimation", x.axis.label = "Fitness score")

As we can see from the above plots, there is positive correlation between the distribution of MCC scores in each class and the respective fitness scores, though not between the 3rd class and the 4th or 5th (it’s negative).

Note also that the 5th class has a lot less models than the rest. Mathematically, we show that the group distributions are indeed different using the Kruskal-Wallis Rank Sum Test:

# Hypothesis testing: Are the location parameters of the distribution of x the same in each group?
kruskal.test(x = df[,"fit.unique"], g = df[, "mcc.class.id"])

    Kruskal-Wallis rank sum test

data:  df[, "fit.unique"] and df[, "mcc.class.id"]
Kruskal-Wallis chi-squared = 135.52, df = 4, p-value < 2.2e-16

As the p-value is less than the significance level \(0.05\), we can conclude that there are significant differences between the different groups of fitness values.

Next, to show that most of the groups are statistically different, we perform pairwise Wilcoxon rank sum tests and draw a heatmap of the each test’s p-value:

pair.res = pairwise.wilcox.test(x = df[,"fit.unique"], g = df[, "mcc.class.id"], p.adjust.method = "BH")
p.value.mat = pair.res$p.value
# same: compare_means(fit.unique ~ mcc.class.id, data = df)

col_fun = colorRamp2(breaks = c(0, 0.05, 0.5, 1), c("green", "white", "orange", "red"))
ht = Heatmap(matrix = p.value.mat, cluster_rows = FALSE, cluster_columns = FALSE,
  na_col = "white", name = "p-value", col = col_fun, column_names_rot = 0,
  row_title = "MCC class id", row_title_rot = 90, row_names_side = "left",
  column_title = "P-values (Pairwise Wilcoxon tests) - MCC-classified fitnesses", column_title_gp = gpar(fontsize = 20),
  cell_fun = function(j, i, x, y, width, height, fill) {
    if (!is.na(p.value.mat[i,j]))
        grid.text(sprintf("%.6f", p.value.mat[i, j]), x, y, gp = gpar(fontsize = 20))
})
draw(ht)

R session info

xfun::session_info()

R version 3.6.1 (2019-07-05)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04.3 LTS, RStudio 1.2.5001

Locale:
  LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
  LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
  LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
  LC_PAPER=en_US.UTF-8       LC_NAME=C                 
  LC_ADDRESS=C               LC_TELEPHONE=C            
  LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

Package version:
  assertthat_0.2.1     backports_1.1.5      base64enc_0.1-3     
  BH_1.69.0.1          bibtex_0.4.2         bookdown_0.14       
  circlize_0.4.8       Ckmeans.1d.dp_4.3.0  cli_1.1.0           
  clue_0.3-57          cluster_2.1.0        codetools_0.2-16    
  colorspace_1.4-1     compiler_3.6.1       ComplexHeatmap_2.0.0
  cowplot_1.0.0        crayon_1.3.4         crosstalk_1.0.0     
  datasets_3.6.1       digest_0.6.21        dplyr_0.8.3         
  DT_0.9               ellipsis_0.3.0       emba_0.1.2          
  evaluate_0.14        fansi_0.4.0          fastmap_1.0.1       
  gbRd_0.4-11          GetoptLong_0.1.7     ggplot2_3.2.1       
  ggpubr_0.2.3         ggrepel_0.8.1        ggsci_2.9           
  ggsignif_0.6.0       GlobalOptions_0.1.1  glue_1.3.1          
  graphics_3.6.1       grDevices_3.6.1      grid_3.6.1          
  gridExtra_2.3        gtable_0.3.0         highr_0.8           
  htmltools_0.4.0      htmlwidgets_1.5.1    httpuv_1.5.2        
  igraph_1.2.4.1       jsonlite_1.6         knitr_1.25          
  labeling_0.3         later_1.0.0          lattice_0.20.38     
  lazyeval_0.2.2       lifecycle_0.1.0      magrittr_1.5        
  markdown_1.1         MASS_7.3-51.4        Matrix_1.2.17       
  methods_3.6.1        mgcv_1.8.29          mime_0.7            
  munsell_0.5.0        nlme_3.1.141         nnet_7.3-12         
  packrat_0.5.0        parallel_3.6.1       pillar_1.4.2        
  pkgconfig_2.0.3      plogr_0.2.0          plyr_1.8.4          
  png_0.1-7            polynom_1.4.0        promises_1.1.0      
  pscl_1.5.2           purrr_0.3.2          R6_2.4.0            
  RColorBrewer_1.1-2   Rcpp_1.0.2           Rdpack_0.11-0       
  reshape2_1.4.3       rje_1.10.10          rjson_0.2.20        
  rlang_0.4.0          rmarkdown_1.16       rstudioapi_0.10     
  scales_1.0.0         shape_1.4.4          shiny_1.4.0         
  sourcetools_0.1.7    splines_3.6.1        stats_3.6.1         
  stringi_1.4.3        stringr_1.4.0        tibble_2.1.3        
  tidyr_1.0.0          tidyselect_0.2.5     tinytex_0.16        
  tools_3.6.1          usefun_0.4.3         utf8_1.1.4          
  utils_3.6.1          vctrs_0.2.0          viridisLite_0.3.0   
  visNetwork_2.0.8     withr_2.1.2          xfun_0.10           
  xtable_1.8-4         yaml_2.2.0           zeallot_0.1.0