ModifiedGap.Rmd

---
title: "ModifiedGapStatistic"
author: "Matt Nolan"
date: "28/02/2019"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

source("ModifiedGap_Functions.R")
```

## Modified gap statistic

Goal is to reimplement matlab code written by Hugh.

```{r}
K.max <- 8
```

Generate gap statistics and calculate significance threshold for simulated uniform distributions. Use these simulated data to calculate threshold values for detection of 'signficant' modes.
```{r}
# Parameters
n_sim <- 1000 # Use ≤ 20 for testing code. Use ≥ 500 for analyses.

# Find all unique dataset lengths ≥ 30.
exp_numbers <- data.sc %>% group_by(id) %>% summarise(n()) %>% select(-id) %>% unique()

# Generate simulated gap data. If n_sim is high this will take a long time to run.
# Will look for a loaded version of sim_tib, then for a saved verion.
# If it finds neither then will build sim_tib.
if(exists("sim_tib")){
  print("sim_tib already generated, remove it if you want to rebuild or reload it")
  } else {
    if (file.exists("sim_distributions.Rda")) {
      load("sim_distributions.Rda")
      } else {
        sim_tib <- tibble(cell_count = exp_numbers[[1]]) %>%
          mutate(sims = pmap(list(cell_count, K.max, n_sim), gap_uniform_sim))
      }
  }

# Write sim_tib to file.
save(sim_tib, file = "sim_distributions.Rda")

# Calculate cut off values from the data
sim_tib <- mutate(sim_tib,
                  threshold_vals = map2(sims, 0.99, calc_thresholds)
                  )
```

Calculate gap statistics, gap differences and K_est for all features for all animals.
```{r}
# Re-organise all-data by id
all_data_r <- data.sc %>%
  group_by(id) %>%
  nest()

# add column with numbers of observations for each group
all_data_r$n <- data.sc %>% group_by(id) %>% summarise(n()) %>% select(-id)
all_data_r <- filter(all_data_r, n >= 30)

# Add thresholds to all_data_r.
all_data_r$sim_vals = list("")
for (i in 1:length(all_data_r$id)) {
  all_data_r$sim_vals[[i]] <- filter(sim_tib, cell_count == as.numeric(all_data_r$n[i,]))
}

# Calculate gap statistics. Restricts analyses to n ≥ 30.
all_data_r <- filter(all_data_r, n >= 30) %>%
  mutate(clusGap = pmap(list(data, sim_vals, K.max), clusGap_AllData))

# Extract k_est and best_gap values for each mouse into a higher level column.
all_data_r <- mutate(all_data_r, K_est = map(clusGap, function(x) x[[4]]))
all_data_r <- mutate(all_data_r, Gap_best = map(clusGap, function(x) x[[5]]))


K_est_all_animals <- unnest(select(all_data_r, id, K_est))
K_est_all_animals$feature <- all_data_r$clusGap[[1]]$property
K_est_all_animals <- spread(K_est_all_animals, feature, K_est)

kableExtra::kable(K_est_all_animals) %>%
  kableExtra::kable_styling(bootstrap_options = "striped") %>%
  kableExtra::save_kable("Tables/K_est.html")
```

Check where clusGap is finding the thresholds. Seems like this needs to be passed to it.
    

Plot data and gap statistic analyses for a given cell and feature.
```{r}
(rheo_plot <- data_plot(all_data_r, "mouse_20140522", "rheo"))

(rheo_logW_plot <- logW_plot(all_data_r, "mouse_20140522", "rheo"))

(rheo_gap_plot <- gap_plot(all_data_r, "mouse_20140522", "rheo"))

(rheo_diff_plot <- diff_plot(all_data_r, "mouse_20140522", "rheo"))
```



Plots for resonance frequency.
Plot data,  gap statistic up to K.max .
```{r}
(resf_plot <- data_plot(all_data_r, "mouse_20140522", "resf"))

(resf_logW_plot <- logW_plot(all_data_r, "mouse_20140522", "resf"))

(resf_gap_plot <- gap_plot(all_data_r, "mouse_20140522", "resf"))

(resf_diff_plot <- diff_plot(all_data_r, "mouse_20140522", "resf"))
```

Plots for input resistance.
Plot data,  gap statistic up to K.max .
```{r}
(ir_plot <- data_plot(all_data_r, "mouse_20140522", "ir"))

(ir_logW_plot <- logW_plot(all_data_r, "mouse_20140522", "ir"))

(ir_gap_plot <- gap_plot(all_data_r, "mouse_20140522", "ir"))

(ir_diff_plot <- diff_plot(all_data_r, "mouse_20140522", "ir"))

cowplot::plot_grid(ir_plot, ir_logW_plot, ir_gap_plot, ir_diff_plot)
```

Gap analysis for simulated bimodal data with n = 30. 

```{r}
n = 30
N = 10
bimod <- function(n) {
  bimod <- tibble(y1 = rnorm(n, 0, 1),  y2 = rnorm(n, 5, 1), w = rbinom(n, 1, .5), x2 = w*y1 + (1-w)*y2)
}

bimod_examples <- tibble(n = rep(30, N)) %>%
  mutate(data = map(n, bimod)) 

# Add thresholds to bimod_examples.
bimod_examples$sim_vals = list("")
for (i in 1:N) {
  bimod_examples$sim_vals[[i]] <- filter(sim_tib, cell_count == as.numeric(all_data_r$n[i,]))
}

bimod_examples <- mutate(bimod_examples, clusGap = pmap(list(data, sim_vals, K.max), clusGap_SimData))
 


 ggplot(bimod, aes(x2)) +
   geom_histogram(bins = 10)
```





Calculate Gap statistic for data clustered on all features.
```{r}
example_data <- return_data(all_data_r, "mouse_20140522")

clus_out <- clusGap_Extra(select(example_data, vm:ahp))
all_k_plot <- ggplot(clus_out, aes(clus_num, gap)) +
                       geom_point()

tsne_all_k <- tsne::tsne(X = as.matrix(select(example_data, vm:ahp)))
colnames(tsne_all_k) <- c("R1", "R2")
tsne_all_k <- as_tibble(tsne_all_k)
kmeans_out <- kmeans(select(example_data, vm:ahp), centers = 7)
tsne_all_k$cluster_id <- as.factor(kmeans_out$cluster)
ggplot(tsne_all_k, aes(R1, R2, colour = cluster_id)) +
  geom_point()
```



Test clustMapDP????
```{r}
example_data <- return_data(all_data_r, "mouse_20140522")
example_data <- as.matrix(example_data$ir)
D = 1
N = length(example_data)
X <- matrix(example_data,D,N);
# Set up Normal-Wishart MAP-DP prior parameters
N0 <- 1;                           # Prior count (concentration parameter)
m0 <- rowMeans(X);                     # Normal-Wishart prior mean
a0 <- 10;                            # Normal-Wishart prior scale
c0 <- 10/N;                          # Normal-Wishart prior degrees of freedom
#B0 <- diag(1./(0.05*rowVars(X)));    # Normal-Wishart prior precision
# r <- clustMapDP(example_data,N0,m0,a0,c0,B0)
```