OAI_Clust

Introduction

Loading necessary libraries and dataset

library(tidyverse)
library(magrittr)
library(broom)
library(e1071)
library(clValid)
library(fpc)
library(gplots)
load("OAI_JSW_Data.RData")

Prepairing data to be analysed

feats <-
  OAI_JSW_Data %>% 
  filter(VISITID == 0, complete.cases(.)) %>% 
  select(-c(VISITID, IDSIDE, ID, SITE, BASEVISITDATE_C, FUVISITDATE, DaysFomStart, DAYSFROMBL,
            DaystoTKR, TimeToEvent, DeathVisit, BLKL2, WOMACTOTAL, 
            WOMAC2YTotalChange, `2YIncreaseSymptoms`, YearTOKL2), 
         -starts_with("Y2")) %>% 
  mutate_at(vars(BLKL, mJSN,lJSN, PainMedication, AnySurgery), as.integer) %>% 
  mutate_if(is.double, funs(scale))

Bootstrapping

Calculate optimal k s by:

1. Sample of 25% of data with replacement
2. Calculate Dunn index for said sample per clustering algorithm for k s between and 2 and a given value
3. Find k that maximizes the Dunn index per clustering algorithm
4. Repeat steps 1 to 3 x times
5. Get mean best k per clustering algorithm

res         <- clustest::find_k(dataset = feats, train_size = 0.25, iters = 28, max_k = 5)
train       <- res$train
test        <- res$test
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

Jaccard Stability Test

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)
}
jac %>% 
  map_df(get_stats, .id = "method") %>%  
  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))

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("WOMACPAIN", "BLKL", "TKREvent")))

|===========================================================================================|100% ~0 s remaining     

a_o_v %>% 
  group_by(dep_var, term) %>% 
  summarise(signif_percent = sum(p.value <= 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))

Print testing set ANOVA

Hierarchical Clustering Heatmap