delta_gap_main.Rmd

---
title: "Modified Gap Statistic for cluster detection (uses delta gap)"
author: "Hugh Pastoll"
date: "31/03/2019"
output:
  html_document: default
  pdf_document: default
---

```{r setup, include=FALSE}
rm(list=ls())
knitr::opts_chunk$set(echo = TRUE)
library(tidyverse)
library(gridExtra)
source("delta_gap_clustering_functions.R")
#wd <- "/Users/hughpastoll/Documents/Research/delta_gap" # Specify working directory
# MN: have removed references to an additional wd so everything runs in the current working directory. Files relevant to this analysis are stored in the subfolder 'ModifiedGap'.
```

##Modified Gap Statistic
The Modified Gap Statistic (MGS) adds a mechanism to the original Gap Statistic (Tibshirani _et._ _al._ 2000) to determine whether a target dataset has a given number of clusters _with a specified false positive rate_. It achieves this by: 

* Computing the _decrease_ in average intra-cluster dispersion between incremental numbers of candidate clusters (the 'delta gaps'), and 
* Comparing the delta gaps for a target dataset against the average delta gap for a large number of reference datasets drawn from uniform distributions.  

If the difference between the target dataset delta gap and average reference delta gap for a given k is greater than would be expected a set proportion (1-'threshold_criterion') of the time if the data were _not_ clustered (i.e. drawn from a uniform distribution), the target dataset is considered to be clustered with k clusters (k~est~=k). 

##Parameters

###Scaling of inputs
The function 'scale_vec()' can be edited to scale data between 0 and 1 (currently active), or scale the data with zero mean, 1 stddev. Scaling method appears to make limited, if any, difference (at least with 1D data), but has not been systematically investigated. 

###General clustering parameters
This code uses the stats::kmeans algorithm to partition data into clusters. Up to 8 clusters per dataset are evaluated as it is considered unlikely that more clusters will be present in SC cells in mice given the number of clusters found in grid cells. 

###Threshold parameters
The thresholds (different for each k and n) are determined by the threshold criterion (0.01). The thresholds are fit for each k across n so that thresholds for any n value (within the range of tested n values, 20 - 100) can be used. The fits use non-linear least squares (a power function) which is very sensitive to the supplied initial fit parameter values. Make sure to change these according to the comments in the 'fit_thresholds' function if d.power changes.   

###Generating uniform reference data
The code generates sim_unif_n_sims >> 1000 simulated datasets drawn from a uniform distribution, calculates the average intra-cluster dispersions & delta gaps and saves the results to 'unif_sims_delta_gaps_and_dispersions.Rda'. The code will first look for this file and skip this step if it already exists. Note that if any simulation parameters change, the file should be deleted so it can be regenerated (this can take 10-15 mins if the number of simulations is around 10000). The simulated uniform data ranges from 20 observations per data to 100, in increments of 10. Note that intermediate numbers of observations (e.g. 36) are not simulated. 

###Generating multimodal test data
To evaluate the sensitivity of the modified gap statistic method, we synthesise multimodal datasets with different numbers of clusters, (k between 2 and 8), standard deviations separating the cluster centres and number of data points. We evaluate these datasets to see whether the clustered structure is detected (k_est > 1). This could be changed to detect the exact number of clusters (k_est = k).

```{r, echo=TRUE}
# All parameters

# General clustering parameters.
k.max <- 8 # Maximum number of clusters

# d.power is the  power to raise Euclidian distances to. Default is 1, Tibshirani used 2.
# NB: changing this requires the threshold fit initial parameter guess values to be changed according to the comments in the 'fit_thresholds' function. 
d.power <- 2 

iter.max <- 100 # Maximum iterations the kmeans algorithm will run. Default is 10
nstart <- 50 # The number of times the kmeans algorithm initialises. Total stability needs > 200

# Threshold parameters
threshold_criterion=0.01 # Specify false positive rate for each evaluated k

# Parameters for uniform sampled simulated datasets. For calculating thresholds and calculated 
# reference dispersions.
sim_unif_dataset_sizes <- seq(from=20, to=100, by=10) # A wide range is useful for stable fits
sim_unif_n_sims <- 20000 # >10000 provides sufficienty good fits

# Parameters for example multimodal simulation illustration of clustering with different k
mm_sims <- list()
mm_sims$n_ex_modes <- 5 # Number of clusters
mm_sims$ex_mode_sep_sd <- c(3.5, 4, 4.5, 5, 5.5, 6) # Separations between cluster centers

# Parameters for evaluating cluster detection sensitivity 
mm_sims$n <- 5000 # Simulated datasets per combination of factors. 5000 provides stable results.
mm_sims$sep_sd_max <- 6 # Range of separations of standard deviations to evaluate
mm_sims$sep_sd_min <- 3
mm_sims$sep_sd_incr <- 0.5

mm_sims$n_data_max <- 70 # Range of multimodal dataset sizes to evaluate.
mm_sims$n_data_min <- 20
mm_sims$n_data_incr <- 10

mm_sims$k_max <- 8 # Range of numbers of clusters (k) to evaluate
mm_sims$k_min <- 2
mm_sims$k_incr <- 1

# Other variables (help with import, categorisation and plotting)
mouse_ind <- 8 # Mouse to visualise clustering for. Takes integer values.

props_subthresh <- c("vm", "ir", "sag", "tau", "resf", "resmag")
props_suprathresh <- c("spkthr", "spkmax", "spkhlf", "rheo", "ahp", "fi")
property_names <- c(props_subthresh, props_suprathresh)
props_sub_labels <- c("vm"="Vm", "ir"="IR", "sag"="Sag", "tau"="Mem. tau", 
                    "resf"="Res. freq.", "resmag"="Res. mag.")
props_supra_labels <- c("spkthr"="Spk. thresh.", "spkmax"="Spk. max.",
                        "spkhlf"="Spk. halfwidth", "rheo"="Rheobase", "ahp"="AHP", "fi"="FI")
property_labels <- c(props_sub_labels, props_supra_labels)

figs <- list() # List to hold plots
```
##Establish uniform reference data

###Uniform simulations

```{r, echo=TRUE}
# Generate simulated delta gap data. If sim_unif_n_sims is high this will take a long time 
# to run. Will look for a loaded version of slopes_store, then for a saved version.
# If it finds neither then will build slopes_store (and dispersion_store).
if(exists("delta_gaps_store")){
  print("delta_gaps_store already generated, remove it if you want to rebuild or reload it")
  if (NROW(delta_gaps_store) != sim_unif_n_sims){
    warning("Number of simulated uniform datasets differs from parameter value")
  }
} else {
  if (file.exists("ModifiedGap/unif_sims_delta_gaps_and_dispersions.Rda")) {
    load("ModifiedGap/unif_sims_delta_gaps_and_dispersions.Rda")
    print("Loaded simulated delta gaps and dispersions")
    if (NROW(delta_gaps_store) != sim_unif_n_sims){
      warning("Number of simulated uniform datasets differs from parameter value")
    }
  } else {
    # Create unif_sims_delta_gaps_and_dispersions.Rda
    gen_unif_sims(sim_unif_n_sims, sim_unif_dataset_sizes, k.max=k.max, d.power=d.power) 
    # Load 'delta_gaps_store' and 'dispersion_store' variables
    load("ModifiedGap/unif_sims_delta_gaps_and_dispersions.Rda") 
  }
}
```

###Fitting the delta gap thresholds and dispersions
It's important to check after each regeneration of simulated data that the threshold fits are appropriate. With less than 10000 runs, the thresholds can be quite noisy and it's possible that the fits won't converge. The form of the fitted equations is $y = (a/n^{b}) + c$, where $n$ is the number of data points in the evaluated dataset and $a, b, c$ are the parameters to be fitted.

In addition, to speed up evaluation of the reference dispersion to calculate the modified gap statistic k_est, it's also possible to fit the dispersions obtained from the simulated data and reuse this information so that multiple bootstrap samples don't have to be drawn for every evaluated dataset. Dispersion fits use a different fitting equation: $y = a*log(n) + b$, where $a, b$ are the parameters to be fitted. 

```{r, echo=TRUE}
# Get the delta gap thresholds (and the delta gaps themselves for plotting)
delta_gaps_thresholds <- get_delta_gaps_thresh_vals(delta_gaps_store, 
                                                    threshold_criterion=threshold_criterion)
# Fit delta_slope_thresh_vals
threshold_fit_results <- fit_thresholds(delta_gaps_thresholds, sim_unif_dataset_sizes)

# Find the average dispersion for each k across all simulations (prevents having to re-simulate
# several uniform distributions every time the modified gap statistic is evaluated)
mean_dispersions <- get_mean_dispersions(dispersion_store)
# Fit mean_dispersions with log function
dispersion_fit_results <- fit_dispersions(mean_dispersions, sim_unif_dataset_sizes)

```


```{r, echo=FALSE}
# Prepare tibbles for plotting threshold fits
thresholds <- as_tibble(delta_gaps_thresholds$delta_gaps_thresh_vals[[1]])
thresholds <- add_column(thresholds, n_data_pts=sim_unif_dataset_sizes)
threshold_fits <- threshold_fit_results$threshold_fits[[1]]
threshold_fits <- add_column(threshold_fits, n_data_pts=sim_unif_dataset_sizes)

chart_data <- thresholds %>% gather(key="method", value="threshold", -n_data_pts)
fit_data <- threshold_fits %>% gather(key="method", value="threshold", -n_data_pts)

figs$p_thr <- ggplot() +
  geom_point(data=chart_data, mapping=aes(x=n_data_pts, y=threshold, colour=method)) +
  geom_line(data=fit_data, mapping=aes(x=n_data_pts, y=threshold, colour=method)) +
  labs(x="Dataset size", y="Delta gap threshold", title="Delta gap threshold fits") +
  scale_color_discrete(name="Delta gap k", labels=c(2:8))

# Prepare tibbles for plotting dispersion fits
dispersions_for_plotting <- as_tibble(mean_dispersions)
dispersions_for_plotting <- add_column(dispersions_for_plotting, n_data_pts=sim_unif_dataset_sizes)
dispersion_fits <- dispersion_fit_results$dispersion_fits[[1]]
dispersion_fits <- add_column(dispersion_fits, n_data_pts=sim_unif_dataset_sizes)

chart_data <- dispersions_for_plotting %>% gather(key="method", value="threshold", -n_data_pts)
fit_data <- dispersion_fits %>% gather(key="method", value="threshold", -n_data_pts)

figs$p_disp <- ggplot() +
  geom_point(data=chart_data, mapping=aes(x=n_data_pts, y=threshold, colour=method)) +
  geom_line(data=fit_data, mapping=aes(x=n_data_pts, y=threshold, colour=method)) +
  labs(x="Dataset size", y="Log dispersion", title="Dispersion fits") +
  scale_color_discrete(name="Dispersion k", labels=c(1:8))

figs$pg <- grid.arrange(figs$p_thr, figs$p_disp, ncol=2, heights=unit(2.5, "in"))

ggsave("ModifiedGap/figure_fits.pdf", plot=figs$pg)
```

##Evaluating SC data inputs for clustering
Apply the modified gap statistic method (or 'delta gap' whatever name we decide to give it), to real and simulated data. 

###Evaluate the Stellate Cell recording data
```{r, echo=TRUE}
# Evaluate the stellate cell (sc) data for k_est
fname <- "Data/datatable.txt"
sc_data <- read_tsv(fname, col_types = cols(.default="d", mlpos="c", hemi="c", 
                                            id="c", housing="c", expr="c", patchdir="c"))

sc_data_long <- sc_data %>% filter(totalcells >= 30) %>% select(property_names, id, dvloc) %>% gather(key="property", value="measurement", property_names, na.rm=TRUE)

# Initialise list to store tibbles for cluster evaluation of each combination of animal 
# and measured property 
sc_cluster_data_list <- list()
list_ind <- 0L # To track the list index inside the double loop

for (m in unique(sc_data_long$id)){
  for (p in property_names){
    # The data for one measured property for a single animal
    data <- sc_data_long %>% filter(id==m, property==p)
    
    # Use the fitted parameters to obtain thresholds and dispersions. Note that where there 
    # are measurement NAs in the original data, they will have already been removed, so the 
    # fits are for the right dataset size. 
    thresholds <- get_thresh_from_params(data, k.max,
                                         threshold_fit_results$thresh_params[[1]])
    dispersions <- get_dispersions_from_params(data, k.max,
                                               dispersion_fit_results$dispersion_params[[1]])
    
    # The clustering evaluation outputs
    cluster_out <- get_k_est(data$measurement, kmeans, iter.max=iter.max, nstart=nstart,
                                K.max=k.max, B=NULL, d.power=d.power, thresholds=thresholds,
                                dispersions=dispersions)
    
    # Construct tibble that combines cluster evaluation output with other data in long format
    list_ind <- list_ind + 1 # Increment list index counter
    dimnames(cluster_out$cluster)[[2]] <- as.integer(1:k.max) # to set 'as.tibble()' behaviour
    sc_cluster_data_list[[list_ind]] <- as.tibble(cluster_out$cluster) %>% bind_cols(data) %>% 
      mutate(k_est=rep(cluster_out$k_est, nrow(cluster_out$cluster))) %>%
      gather(key="k_eval", value="cluster", dimnames(cluster_out$cluster)[[2]]) 
  }
}
# Aggregated data in long format. Most granular level is single measurement cluster membership
# for a specified k. This is the main input to further SC data analyses. 
sc_cluster_data_long <- bind_rows(sc_cluster_data_list) %>% 
  mutate(k_eval=as.integer(k_eval), cluster=as.factor(cluster)) %>% # Change class for plotting
  mutate(prop_type=case_when(property %in% props_subthresh ~ "subthr", TRUE ~ "suprathr"))
```

```{r, echo=FALSE}
figs$sc_cl <- sc_cluster_data_long %>% filter(id==unique(id)[mouse_ind], 
                                              prop_type=="subthr") %>% 
  ggplot() + geom_point(mapping=aes(x=dvloc, y=measurement, color=cluster), size=0.5) +
  scale_color_discrete(name="Cluster") +
  facet_grid(property + k_est ~ k_eval, scales="free_y", switch="y",
             labeller=label_value) +
  labs(y="Values (mixed units - not shown)", x="Dorsoventral location (mm)") +
  theme_bw() + theme(axis.text.y = element_blank(), axis.ticks.y = element_blank(),
                     panel.grid = element_blank()) +
  scale_x_continuous(labels=c("0", "2.5"), breaks=c(0, 2500), limits=c(0, 2500))

figs$sc_cl
ggsave("ModifiedGap/figure_sc_clustering_example.pdf", plot=figs$sc_cl)

# Plot k_est counts for subthreshold and suprathreshold properties
figs$sb <- sc_cluster_data_long %>% filter(prop_type=="subthr") %>% distinct(id, property, k_est) %>% 
  ggplot() + geom_bar(aes(x=k_est)) + 
  facet_wrap(vars(property), ncol=6, labeller=as_labeller(props_sub_labels)) +
  labs(y="Count", x=" ") + 
  scale_x_continuous(labels=c(rep(" ", k.max)), breaks=1:k.max,
                     limits=c(0.5, k.max+0.5))

figs$sp <- sc_cluster_data_long %>% filter(prop_type=="suprathr") %>% 
  distinct(id, property, k_est) %>% 
  ggplot() + geom_bar(aes(x=k_est)) + 
  facet_wrap(vars(property), ncol=6, labeller=as_labeller(props_supra_labels))+
  labs(x=bquote(~"Estimated number of clusters, "~k[est]), y="Count") +
  scale_x_continuous(labels=c("1", rep("", k.max-2), toString(k.max)), breaks=1:k.max,
                     limits=c(0.5, k.max+0.5))

figs$sc_k <- grid.arrange(figs$sb, figs$sp, ncol=1, nrow=2)
figs$sc_k
ggsave("ModifiedGap/figure_sc_data_k_est.pdf", plot=figs$sc_k)
```


### Evaluate uniformly distributed data
Check the results of using uniformly sampled data that has the same size as the real sc data. This is mostly a sanity check - we do not systematically evaluate uniform data inputs across a wide range of n. 
```{r, echo=FALSE}
# Quantify the false positive rates for uniform distributed datasetsdatasets

# Replace actual property measurement values with uniformly sampled data
unif_data_long <- sc_data_long %>% mutate(measurement=runif(NROW(sc_data_long)))

# Initialise list to store tibbles for cluster evaluation of each combination of animal 
# and measured property 
unif_cluster_data_list <- list()
list_ind <- 0L # To track the list index inside the double loop

for (m in unique(unif_data_long$id)){
  for (p in property_names){
    # The data for one measured property for a single animal
    data <- unif_data_long %>% filter(id==m, property==p)
    
    # Use the fitted parameters to obtain thresholds and dispersions. Note that where there 
    # are measurement NAs in the original data, they will have already been removed, so the 
    # fits are for the right dataset size. 
    thresholds <- get_thresh_from_params(data, k.max,
                                         threshold_fit_results$thresh_params[[1]])
    dispersions <- get_dispersions_from_params(data, k.max,
                                               dispersion_fit_results$dispersion_params[[1]])
    
    # The clustering evaluation outputs
    cluster_out <- get_k_est(data$measurement, kmeans, iter.max=iter.max, nstart=nstart,
                                K.max=k.max, B=NULL, d.power=d.power, thresholds=thresholds,
                                dispersions=dispersions)
    
    # Construct tibble that combines cluster evaluation output with other data in long format
    list_ind <- list_ind + 1 # Increment list index counter
    dimnames(cluster_out$cluster)[[2]] <- as.integer(1:k.max) # to set 'as.tibble()' behaviour
    unif_cluster_data_list[[list_ind]] <- as.tibble(cluster_out$cluster) %>% bind_cols(data) %>% 
      mutate(k_est=rep(cluster_out$k_est, nrow(cluster_out$cluster))) %>%
      gather(key="k_eval", value="cluster", dimnames(cluster_out$cluster)[[2]]) 
  }
}
# Aggregated data in long format. Most granular level is single measurement cluster membership
# for a specified k. This is the main input to further uniform data analyses. 
unif_cluster_data_long <- bind_rows(unif_cluster_data_list) %>% 
  mutate(k_eval=as.integer(k_eval), cluster=as.factor(cluster)) %>% # Change class for plotting
  mutate(prop_type=case_when(property %in% props_subthresh ~ "subthr", TRUE ~ "suprathr"))
```

```{r, echo=FALSE}
figs$unif_cl <- unif_cluster_data_long %>% filter(id==unique(id)[mouse_ind], 
                                                  prop_type=="subthr") %>% 
  ggplot() + geom_point(mapping=aes(x=dvloc, y=measurement, color=cluster), size=0.5) +
  scale_color_discrete(name="Cluster") +
  facet_grid(property + k_est ~ k_eval, scales="free_y", switch="y",
             labeller=label_value) +
  labs(y="Simulated values", x="Simulated location (mm)") +
  theme_bw() + theme(axis.text.y = element_blank(), axis.ticks.y = element_blank(),
                     panel.grid = element_blank()) +
  scale_x_continuous(labels=c("0", "2.5"), breaks=c(0, 2500), limits=c(0, 2500))

figs$unif_cl
ggsave("ModifiedGap/figure_unif_clustering_example.pdf", plot=figs$unif_cl)

# Plot k_est counts for subthreshold and suprathreshold properties
figs$unif_k <- unif_cluster_data_long %>% 
  distinct(id, property, k_est) %>% 
  ggplot() + geom_bar(aes(x=k_est)) + 
  facet_wrap(vars(property), ncol=6, labeller=as_labeller(property_names)) + # Give NA labels
  labs(x=bquote(~"Estimated number of clusters, "~k[est]), y="Count") +
  scale_x_continuous(labels=c(rep(" ", k.max)), breaks=1:k.max,
                     limits=c(0.5, k.max+0.5))

figs$unif_k
ggsave("ModifiedGap/figure_unif_data_k_est.pdf", plot=figs$unif_k)
```


##Evaluate multimodal (simulated clustered) data

###Illustrate k evaluation for multimodal data
First, generate example multimodal data with the same size as the experimental data (and with k=mm_sims$n_ex_modes. This allows visualising how clustering evaluations compare between experimental and simulated multimodal data. 

Changing the SC mouse index (in the Parameters section) will apply the relevant set of dorsoventral locations to the simulated data for plotting. 
```{r, echo=FALSE}
# This chunk is a bit of a hack! It produces multimodal data in place of the actual SC data so that the same clustering evaluation code can be applied to the multimodal data. 

# Keep SC data for now, replace measurement values inside loop
mm_data_long <- sc_data_long %>% filter(property %in% props_subthresh) # Convenience: only want 6 values

# Initialise list to store tibbles for cluster evaluation of each combination of animal 
# and measured property 
mm_cluster_data_list <- list()
list_ind <- 0L # To track the list index inside the double loop

for (m in unique(mm_data_long$id)){
  mm_sims$ex_sd_ind <- 0L
  for (p in props_subthresh){
    # The data for one measured property for a single animal
    data <- mm_data_long %>% filter(id==m, property==p)
    
    # Use the fitted parameters to obtain thresholds and dispersions. Note that where there 
    # are measurement NAs in the original data, they will have already been removed, so the 
    # fits are for the right dataset size. 
    thresholds <- get_thresh_from_params(data, k.max,
                                         threshold_fit_results$thresh_params[[1]])
    dispersions <- get_dispersions_from_params(data, k.max,
                                               dispersion_fit_results$dispersion_params[[1]])
    
    # we are mapping each sd onto a different subthresh prop.
    mm_sims$ex_sd_ind <- mm_sims$ex_sd_ind + 1
    
    # The clustering evaluation outputs. Instead of sending the real data for evaluation use the
    # synthesised multimodal data.
    data$measurement <- get_multimodal_data_sample(NROW(data), n_modes=mm_sims$n_ex_modes,
                                           mm_sims$ex_mode_sep_sd[mm_sims$ex_sd_ind])
    
    cluster_out <- get_k_est(data$measurement, kmeans, iter.max=iter.max, nstart=nstart,
                                K.max=k.max, B=NULL, d.power=d.power, thresholds=thresholds,
                                dispersions=dispersions)
    
    # Construct tibble that combines cluster evaluation output with other data in long format
    list_ind <- list_ind + 1 # Increment list index counter
    dimnames(cluster_out$cluster)[[2]] <- as.integer(1:k.max) # to set 'as.tibble()' behaviour
    mm_cluster_data_list[[list_ind]] <- as.tibble(cluster_out$cluster) %>% bind_cols(data) %>% 
      mutate(k_est=rep(cluster_out$k_est, nrow(cluster_out$cluster))) %>%
      gather(key="k_eval", value="cluster", dimnames(cluster_out$cluster)[[2]]) 
  }
}
# Aggregated data in long format. Most granular level is single measurement cluster membership
# for a specified k. This is the main input to further multimodal data analyses. 
mm_cluster_data_long <- bind_rows(mm_cluster_data_list) %>% 
  mutate(k_eval=as.integer(k_eval), cluster=as.factor(cluster)) %>% # Change class for plotting
  mutate(prop_type=case_when(property %in% props_subthresh ~ "subthr", TRUE ~ "suprathr")) %>% 
  mutate(sd=case_when(property=="vm" ~ "3.5", property=="ir" ~ "4.0", 
                          property=="sag" ~ "4.5", property=="tau" ~ "5.0",
                          property=="resf" ~ "5.5", property=="resmag" ~ "6.0")) %>% 
  select(-property)
```

```{r, echo=FALSE}
figs$mm_cl <- mm_cluster_data_long %>% 
  filter(id==unique(id)[mouse_ind], prop_type=="subthr") %>% 
  ggplot() + geom_point(mapping=aes(x=dvloc, y=measurement, color=cluster), size=0.5) +
  scale_color_discrete(name="Cluster") +
  facet_grid(sd + k_est ~ k_eval, scales="free_y", switch="y",
             labeller=label_value) +
  labs(y="Standard deviations + Estimated number of clusters", x="Simulated location (mm)") +
  theme_bw() + theme(axis.text.y = element_blank(), axis.ticks.y = element_blank(),
                     panel.grid = element_blank()) +
  scale_x_continuous(labels=c("0", "2.5"), breaks=c(0, 2500), limits=c(0, 2500))

figs$mm_cl
ggsave("ModifiedGap/figure_mm_clustering_example.pdf", plot=figs$mm_cl)
```


###Quantify false positive rates and detection power for multimodal distributions
Use different combinations of k, n and standard deviation separation between cluster centres to generate a wide range of multimodal datasets to evaluate for clustering. 
```{r, echo=TRUE}
# Get estimates for k for simulated clustered (multimodal) data. 
# If n_mm_sims is high this will take a long time to run.
# Will look for a loaded version of k_estimates_mm, then for a saved version.
# If it finds neither then will build k_estimates_mm.
if(exists("k_estimates_mm")){
  print("k_estimates_mm already exists, remove it if you want to rebuild or reload it")
  } else {
    if (file.exists("ModifiedGap/k_estimates_mm.Rda")) {
      load("ModifiedGap/k_estimates_mm.Rda")
      print("Loaded multimodal data k estimates. Delete file if you want to rebuild")
      } else {
        # Evaluate k_est for each combination of factors
        mm_sims$sep_sd_vec <- seq(from=mm_sims$sep_sd_min, 
                             to=mm_sims$sep_sd_max, by=mm_sims$sep_sd_incr)
        mm_sims$n_data_vec <- seq(from=mm_sims$n_data_min, 
                             to=mm_sims$n_data_max, by=mm_sims$n_data_incr)
        mm_sims$k_vec <- seq(from=mm_sims$k_min, to=mm_sims$k_max, by=mm_sims$k_incr)
        
        k_estimates_mm <- get_mm_sim_k_ests(mm_sims$n, mm_sims$k_vec, mm_sims$n_data_vec, 
                                  mm_sims$sep_sd_vec, threshold_fit_results,
                                  dispersion_fit_results)
      }
  }
```

```{r, echo=FALSE}
# Which datasets have been evaluated as clustered (k_est > 1)?
k_est_clustered_mm <- k_estimates_mm > 1

# Average across n_sims dimension to get the proportion of datasets that are clustered. 
k_est_proportions_mm <- apply(k_est_clustered_mm, c(2, 3, 4), mean)

# Change n_data to integer for correct ordering when plotting
k_est_proportions_mm_chart_data <- as.tibble(as.table(k_est_proportions_mm)) %>%
  mutate(n_data=as.integer(n_data))

figs$mm1 <- ggplot(data=k_est_proportions_mm_chart_data, 
                   aes(x=k, y=n, group=sd_sep, colour=sd_sep)) +
  geom_line() + geom_point() + facet_wrap(vars(n_data), labeller=label_bquote(n==.(n_data))) + 
  scale_y_continuous(breaks=seq(0, 1, 0.2), labels=c(0, rep("", 4), 1)) +
  labs(x="Number of simulated clusters (k)", y=bquote(p(k[est]>1)),
       title="Effect of dataset size on cluster detection", 
       colour=bquote(sigma~"s"))
figs$mm1
ggsave("ModifiedGap/figure_mm_detection_by_n_fits.pdf", plot=figs$mm1)

figs$mm2 <- ggplot(data=k_est_proportions_mm_chart_data, aes(x=n_data, y=n, group=sd_sep, colour=sd_sep)) +
  geom_line() + geom_point() + facet_wrap(vars(k), labeller=label_bquote(k==.(k))) +
  scale_y_continuous(breaks=seq(0, 1, 0.2), labels=c(0, rep("", 4), 1)) +
  labs(x="Number of data points in simulated datasets", y=bquote(p(k[est]>1)),
       title="Effect of number of clusters on cluster detection", 
       colour=bquote(sigma~"s"))
figs$mm2
ggsave("ModifiedGap/figure_mm_detection_by_k_fits.pdf", plot=figs$mm2)
```