Calinski_Harabasz_omega <-  function(max_k = 5, distance_matrix = NULL, hierarchical_cluster_solution = NULL) {
  
  require(fpc)
  require(dplyr)
  
  #VRC values #Calinski et al.
  VRC <- rep(2:(max_k+1)) %>%
    map(~ cluster.stats(d = distance_matrix, cutree(hierarchical_cluster_solution, k=.))$ch) %>%
    unlist() %>% as_tibble()
  
  names(VRC)[1] <- "VRC"
  
  #cluster numbers
  cluster_number_stats <- tibble(.rows=max_k)
  cluster_number_stats$clusters <- rep(2:(max_k+1))
  cluster_number_stats <- cbind(cluster_number_stats, VRC)
  
  
  #add Omega VRC
  cluster_number_stats$VRC_omega <- NA
  
  i=2
  while(i<=(nrow(cluster_number_stats)-1)) {
    cluster_number_stats$VRC_omega[i] <- (cluster_number_stats$VRC[i+1]-cluster_number_stats$VRC[i])-(cluster_number_stats$VRC[i]-cluster_number_stats$VRC[i-1])
    i <- i+1
  }
  
  
  return(cluster_number_stats)
  
}