BC_Clust

Loading libraries and data

library(tidyverse)
library(magrittr)
library(fpc)
library(e1071)
library(gplots)
load("DataOnco/SJ_BC_Data.RData")

Getting dataset ready

feats <- 
  SJ_BC_Data %>% 
  select(-(ImageID:HER2P)) %>% 
  select_if(funs(!anyNA(.))) %>%
  mutate(Oncotype = SJ_BC_Data$Oncotype) %>% 
  filter(!is.na(Oncotype)) %>% 
  mutate_if(is.numeric, scale)

Getting optimal k

res   <- clustest::partition(dataset = feats, train_size = 0.25)
train <- res$train
test  <- res$test

res         <- clustest::find_k(dataset = train, iters = 30, max_k = 6)
k_stats     <- res$k_stats
best_k      <- res$best_k 
all_best_ks <- res$all_best_ks
c_samples   <- res$c_samples

remove(res)

Plotting optimal k s

k_stats %>% 
  ggplot(aes(x = method, y = mean_best_k)) +
    geom_col(aes(fill = method)) +
    guides(fill = F) +
    geom_errorbar(aes(ymin = ci_min, ymax = ci_max))

best_k

Fuzzy c-Means Clustering Bootstrap Interface

cmeansCBI <- function(data, k) {
  result <- cmeans(x = data, center = k)
  
  list(result = result, 
       partition = result$cluster,
       clusterlist = map(1:k, ~. == result$cluster),
       clustermethod = "fcmeans",
       nc = k)
}

Jaccard evaluation of all algorithms

jac <- list(
  KMeans  = clusterboot(train, B = 30, clustermethod = kmeansCBI, k = best_k$KMeans)$bootresult,
  FCMeans = clusterboot(train, B = 30, clustermethod = cmeansCBI, k = best_k$FCMeans)$bootresult,
  HCWard  = clusterboot(train, B = 30, clustermethod = hclustCBI, k = best_k$HCWard, 
                        method = "ward.D")$bootresult,
  HCSing  = clusterboot(train, B = 30, clustermethod = hclustCBI, k = best_k$HCSing, 
                        method = "single")$bootresult,
  HCComp  = clusterboot(train, B = 30, clustermethod = hclustCBI, k = best_k$HCComp, 
                        method = "complete")$bootresult)

Plotting stability tests summary

get_stats <- function(jac_mat){
  jac_mat %>% 
    t() %>% 
    as_tibble() %>% 
    summarise_all(funs(ci_min = t.test(.)$conf.int[1],
                       ci_max = t.test(.)$conf.int[2],
                       mean_jaccard = mean(.))) %>% 
    gather() %>% 
    separate(key, into = c("cluster", "stat"), extra = "merge") %>% 
    spread(stat, value)
}
mean_jac <- 
  jac %>% 
  map_df(get_stats, .id = "method") %>% 
  mutate_at(vars(ci_max, ci_min), funs(if_else(is.na(.), 0, .))) 
mean_jac %>% 
  ggplot(aes(x = cluster, y = mean_jaccard)) +
    facet_wrap(~method) +
    geom_col(aes(fill = cluster)) +
    guides(fill = F) +
    geom_errorbar(aes(ymin = ci_min, ymax = ci_max))

mean_jac %>% 
  group_by(method) %>% 
  summarise(ci_min = t.test(mean_jaccard)$conf.int[1], 
            ci_max = t.test(mean_jaccard)$conf.int[2],
            mean_jaccard = mean(mean_jaccard)) %>% 
  ggplot(aes(x = method, y = mean_jaccard)) +
    geom_col(aes(fill = method)) +
    guides(fill = F) +
    geom_errorbar(aes(ymin = ci_min, ymax = ci_max))

ANOVA

Plot training sets ANOVA

a_o_v <- 
  c_samples %>% 
  group_by(sample) %>% 
  do(clustest::var_clust_aov(., methods =  c("KMeans", "FCMeans", "HCWard", "HCSing", "HCComp"), 
                                dep_vars = c("Oncotype", "lesion")))
a_o_v %>% 
  mutate(p.value = log10(p.value)) %>% 
  group_by(dep_var, term) %>% 
  summarise(signif_percent = sum(p.value <= log10(0.05))/n(),
            ci_min = -t.test(p.value)$conf.int[1],
            ci_max = -t.test(p.value)$conf.int[2],
            mean_p_value = -mean(p.value)) %T>%
  print() %>% 
  ggplot(aes(x = term, y = mean_p_value, fill = term)) +
    facet_wrap(~dep_var, ncol = 2) +
    geom_col() +
    guides(fill = F) +
    geom_errorbar(aes(ymin = ci_min, ymax = ci_max)) +
    ylab("-mean log10(p-value)") +
    xlab("method")

Print testing set ANOVA

d_test <- dist(test)
cls <- map(c("ward.D", "single", "complete"), ~hclust(d_test, .))
names(cls) <- c("HCWard", "HCSing", "HCComp")
clustered_test <- 
  test %>% 
  bind_cols(
    KMeans  = kmeans(., best_k$KMeans)$cluster,
    FCMeans = cmeans(., best_k$FCMeans)$cluster,
    HCWard  = 
      cls$HCWard %>% 
      cutree(best_k$HCWard),
    HCSing  = 
      cls$HCSing %>% 
      cutree(best_k$HCSing),
    HCComp  = 
      cls$HCComp %>% 
      cutree(best_k$HCComp)
  )
clustered_test %>%  
  clustest::var_clust_aov(c("KMeans", "FCMeans", "HCWard", "HCSing", "HCComp"),
                          c("Oncotype", "lesion")) %>% 
  filter(p.value <= 0.05)

Hierarchical Clustering Heatmap

plotfun <- function(method){
  
  colorvec <- rainbow(best_k[[method]])[clustered_test[[method]]]
  
  cindx <- get('colInd', parent.frame(2))
  marg  <- get('margins', parent.frame(2))
  
  par(mar = c(0.25, 0, 0, marg[2]), cex = 0.66)
  image(cbind(seq_along(colorvec)), col = colorvec[cindx], axes = F)
  mtext(method, side = 4, line = 1, cex=0.4,las=1, col="black")
}
hmfun <- function(data, method){
  data %>% 
    t() %>% 
    heatmap.2(dendrogram = "col", Colv = as.dendrogram(cls[[method]]), trace = "none", 
            scale = "none", breaks = seq(-2, 2, length.out = 101), 
            col = colorRampPalette(colors = c("cyan","blue","black","red","yellow")),
            density.info = "density", cexRow = 0.4, extrafun = function(x) plotfun(method = method), 
            lmat = rbind(c(4, 3), 
                         c(0, 5), 
                         c(2, 1)), 
            lwid = c(1.2, 4), 
            lhei = c(1.5, 0.15, 4))
}
test %T>% 
  hmfun(method = "HCWard") %>% 
  hmfun(method = "HCComp")