Figure2_Analysis.Rmd

---
title: "RampCodes_Figure2"
author: "Sarah Tennant & Matt Nolan"
date: "20/10/2021"
output: html_document
---

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


### ------------------------------------------------------------------------------------------ ### 

## Do ramping neurons encode speed, position or acceleration?

### ------------------------------------------------------------------------------------------ ### 



## Fit general linear mixed effect models (GLMER) using LME4 package to evaluate contributions of speed, acceleration and position to firing rate

Note here, there may be some overfitting with the mixed effect model. This is likely because some  coefficients likely do not vary with position. However, we can't address this by tailoring the model to each cell as we want to treat each cell the same way. The standardized coefficients should still be interpretable.
See discussion in ?isSingular and here: https://stats.stackexchange.com/questions/378939/dealing-with-singular-fit-in-mixed-models

spatial_firing$spikes_in_time is a list column. Each row corresponds to a single cell and contains a list of 6 data frames. These correspond to firing rate, position, acceleration, speed, trial number and trial type. The data cover a full behavioral session. Trial types are beaconed (0), non-beaconed and probe.
```{r, warning=FALSE}
# spatial_firing <- spatial_firing[1:5,]  %>%
#  mutate(o_mm_b = future_map2(spikes_in_time, 0, mm_fit_function),
#         #o_mm_p = map2(spikes_in_time, 1, mm_fit_function),
#         o_mm_p_b = future_map2(o_mm_b, session_id, mm_pvalues)) %>%
#         #o_mm_p_p = map(o_mm_p, session_id, mm_pvalues)) %>%
#  unnest_wider(o_mm_p_b, names_sep = "_", names_repair = "universal")
#  #unnest_wider(o_mm_p_p, names_sep = "_", names_repair = "universal")

df<-data.frame()
for (i in 1:nrow(spatial_firing)) {
  print(as.character(i))
  spatial_firing_i = spatial_firing[i,]
  o_mm_b <- mm_fit_function(spatial_firing_i$spikes_in_time[[1]], 0)
  o_mm_p_b <- mm_pvalues(o_mm_b, spatial_firing_i$session_id[1])
  o_b_mod_coefs <- std_coef(o_mm_b)

  spatial_firing_i$o_b_mod_coefs_pos[1] <- as.double(o_b_mod_coefs[1])
  spatial_firing_i$o_b_mod_coefs_speed[1] <- as.double(o_b_mod_coefs[2])
  spatial_firing_i$o_b_mod_coefs_accel[1] <- as.double(o_b_mod_coefs[3])
  spatial_firing_i$o_mm_p_b_pos[1] <- as.double(o_mm_p_b[1])
  spatial_firing_i$o_mm_p_b_speed[1] <- as.double(o_mm_p_b[2])
  spatial_firing_i$o_mm_p_b_accel[1] <- as.double(o_mm_p_b[3])
  df<-bind_rows(df, spatial_firing_i)
}
spatial_firing <- df
rm(df, spatial_firing_i)



# number with fits / NAs
table(is.na(spatial_firing$o_mm_p_b_pos))
```

# Adjust p values for multiple comparisons
```{r}
spatial_firing <- spatial_firing %>%
    mutate(o_mm_p_b_pos = p.adjust(o_mm_p_b_pos, "BH"),
         o_mm_p_b_speed = p.adjust(o_mm_p_b_speed, "BH"),
         o_mm_p_b_accel = p.adjust(o_mm_p_b_accel, "BH"))
```




```{r}
# Can load data here to save time. Important to do this after the step above, otherwise multiple comparison adjustment happens twice.
load_sf = 0
if (load_sf==1) {
  spatial_firing <- readRDS(file="data_in/SpatialFiring_with_Results.Rda")
}

```

### ----------------------------------------------------------------------------------------- ###

  
## Categorise neurons based on significant model coefficients

The idea here is to categorise each neuron according to whether it has coefficients for P, S or A that are 'significant' at a threshold of 0.01.
```{r}
spatial_firing <- spatial_firing  %>%
    mutate(final_model_o_b  = pmap(list(o_mm_p_b_pos, o_mm_p_b_speed, o_mm_p_b_accel, 0.01), coef_comparison))
```


### Plot model selection results performed in Figure1_Analysis.Rmd

Standardized coefficients are in (see above):
spatial_firing$o_b_mod_coefs_pos
spatial_firing$o_b_mod_coefs_speed
spatial_firing$o_b_mod_coefs_accel

First organise the data.
```{r}
# Distributions of fit categories for all neurons with ramping activity before the reward zone
(all_neurons_groups <- spatial_firing %>% 
  filter(lm_group_b == "Positive" | lm_group_b == "Negative") %>%
   select(final_model_o_b) %>%
   make_coeffs_table())

# All neurons in the group
sum(all_neurons_groups$num)
# Positional neurons
sum(all_neurons_groups[grep("P", all_neurons_groups$ramp_id),]$num)

# For positive slope neurons
pos_neurons_groups <- spatial_firing %>%
  filter(lm_group_b == "Positive") %>%
  select(final_model_o_b) %>%
  make_coeffs_table() %>%
  mutate(ramp_type = "Positive")

# For negative slope neurons
neg_neurons_groups <- spatial_firing %>%
  filter(lm_group_b == "Negative") %>%
  select(final_model_o_b) %>%
  make_coeffs_table() %>%
  mutate(ramp_type = "Negative")

# For unclassified slope neurons
uc_neurons_groups <- spatial_firing %>%
  filter(lm_group_b == "Unclassified") %>%
  select(final_model_o_b) %>%
  make_coeffs_table() %>%
  mutate(ramp_type = "Unclassified")

# Combine
mixed_ramps_2 <- bind_rows(pos_neurons_groups, neg_neurons_groups, uc_neurons_groups)
```

Plot model results.
```{r}
# plot data
level_order <- c("P", "S", "A", "PS", "PA","SA", "PSA", "None")
ggplot(mixed_ramps_2, aes(x= factor(ramp_type), y = perc, fill=factor(ramp_id, level = level_order))) +
  geom_bar(stat="identity",width = 0.9, alpha = .7) +
  labs(y = "Percent of neurons", x="") +
  scale_fill_manual(values=c("firebrick1","gold", "dodgerblue2", "darkorange", "darkorchid1", "chartreuse3", "darkslategray", "grey78")) +
  geom_text(aes(label = num), hjust = 0.7, vjust = 0.2, size = 5, position = position_stack(vjust= 0.5)) +
  #scale_fill_brewer(palette= "RdYlBu") +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 50, vjust=0.65)) +
  theme(axis.text.x = element_text(size=17),
        axis.text.y = element_text(size=16),
        legend.position="bottom", 
        legend.title = element_blank(),
        text = element_text(size=16), 
        legend.text=element_text(size=16), 
        axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)))


if (save_figures == 1) {
  ggsave(file = "plots/LMERProportions_beaconed_update.png", width = 4, height = 5.5)
}
```




### plot distribution of coefficients for all cells
_split by positive and negative sloping neurons_

make unique id for each neuron (session_id + cluster_id)
```{r}
# run on all cells
spatial_firing <- spatial_firing  %>%
  mutate(unique_id2 = pmap(list(session_id, cluster_id, sep = "_"), paste))

```

1. Reorganize data
```{r}
spatial_firing_coefs <- spatial_firing %>% 
  select(session_id, cluster_id, unique_id, o_b_mod_coefs_pos, o_b_mod_coefs_speed, o_b_mod_coefs_accel, final_model_o_b, lm_group_b) %>%
  unnest(cols = c(unique_id, o_b_mod_coefs_pos, o_b_mod_coefs_speed, o_b_mod_coefs_accel,
    final_model_o_b, lm_group_b)) %>%
  as.tibble() %>%
  gather("coef_type", "coef", o_b_mod_coefs_pos, o_b_mod_coefs_speed, o_b_mod_coefs_accel, factor_key = TRUE) %>%
  mutate(coef_type = sub("o_b_mod_coefs_pos", "P", coef_type),
         coef_type = sub("o_b_mod_coefs_speed", "S", coef_type),
         coef_type = sub("o_b_mod_coefs_accel", "A", coef_type)) %>%
  rename(group = "final_model_o_b")
  
spatial_firing_save <- select(spatial_firing, session_id, cluster_id, unique_id, o_b_mod_coefs_pos, o_b_mod_coefs_speed, o_b_mod_coefs_accel, final_model_o_b, lm_group_b, lm_group_nb, lm_group_p, lm_group_b_h, lm_group_nb_h, lm_group_p_h, asr_b_o_rewarded_fit_r.squared, asr_b_o_rewarded_fit_p.value, asr_b_o_rewarded_fit_intercept, asr_b_o_rewarded_fit_slope, asr_nb_o_rewarded_fit_r.squared, asr_nb_o_rewarded_fit_p.value, asr_nb_o_rewarded_fit_intercept, asr_nb_o_rewarded_fit_slope, asr_p_o_rewarded_fit_r.squared, asr_p_o_rewarded_fit_p.value, asr_p_o_rewarded_fit_intercept, asr_p_o_rewarded_fit_slope, asr_b_h_rewarded_fit_r.squared, asr_b_h_rewarded_fit_p.value, asr_b_h_rewarded_fit_intercept, asr_b_h_rewarded_fit_slope, asr_nb_h_rewarded_fit_r.squared, asr_nb_h_rewarded_fit_p.value, asr_nb_h_rewarded_fit_intercept, asr_nb_h_rewarded_fit_slope, asr_p_h_rewarded_fit_r.squared, asr_p_h_rewarded_fit_p.value, asr_p_h_rewarded_fit_intercept, asr_p_h_rewarded_fit_slope) %>%
  unnest(cols = c(unique_id, o_b_mod_coefs_pos, o_b_mod_coefs_speed, o_b_mod_coefs_accel,
    final_model_o_b, lm_group_b, lm_group_nb, lm_group_p, lm_group_b_h, lm_group_nb_h, lm_group_p_h)) %>%
  as.tibble()
spatial_firing_save$o_b_mod_coefs_pos <- as.numeric(spatial_firing_save$o_b_mod_coefs_pos)
spatial_firing_save$o_b_mod_coefs_speed <- as.numeric(spatial_firing_save$o_b_mod_coefs_speed)
spatial_firing_save$o_b_mod_coefs_accel <- as.numeric(spatial_firing_save$o_b_mod_coefs_accel)
spatial_firing_save$cluster_id <- as.numeric(spatial_firing_save$cluster_id)
```

3. save to csv file _this is for matching the coefficient values to plots of instantaneous rates from python_
```{r}
write.table(spatial_firing_save, "data_out/all_results_coefficients.csv", quote=FALSE, sep="\t")

```

_Plot coefficients for all cells_

5a. Split by positive slope _all cells_
```{r}

(pos_coef_plot <- spatial_firing_coefs %>%
  filter(lm_group_b == "Positive") %>%
  standard_plot())

if (save_figures == 1) {
   ggsave(file = "plots/CoefficientValues_positivecells.png", width = 2.5, height = 2.5)
}

one.way_p <- aov(coef ~ as.factor(coef_type) + Error(as.factor(unique_id)), data = filter(spatial_firing_coefs, lm_group_b == "Positive"))
summary(one.way_p)

(t_tests_p <- spatial_firing %>%
  filter(lm_group_b == "Positive") %>%
  summarise(ttest = list(t.test(o_b_mod_coefs_pos, o_b_mod_coefs_accel, paired = TRUE)$p.value,
                         t.test(o_b_mod_coefs_pos, o_b_mod_coefs_speed, paired = TRUE)$p.value,
            t.test(o_b_mod_coefs_accel, o_b_mod_coefs_speed, paired = TRUE)$p.value)) %>%
  unlist() %>%
  as.double() %>%
  p.adjust("bonferroni"))
```

5b. Split by negative slope _all cells_
```{r}
(neg_coef_plot <- spatial_firing_coefs %>%
  filter(lm_group_b == "Negative" ) %>%
    standard_plot())

if (save_figures == 1) {
  ggsave(file = "plots/CoefficientValues_negativecells.png", width = 2.5, height = 2.5)
}

one.way_n <- aov(coef ~ coef_type + Error(as.factor(unique_id)), data = filter(spatial_firing_coefs, lm_group_b == "Negative"))
summary(one.way_n)

(t_tests_n <- spatial_firing %>%
  filter(lm_group_b == "Negative") %>%
  summarise(ttest = list(t.test(o_b_mod_coefs_pos, o_b_mod_coefs_accel, paired = TRUE)$p.value,
                         t.test(o_b_mod_coefs_pos, o_b_mod_coefs_speed, paired = TRUE)$p.value,
            t.test(o_b_mod_coefs_accel, o_b_mod_coefs_speed, paired = TRUE)$p.value)) %>%
  unlist() %>%
  as.double() %>%
  p.adjust("bonferroni"))
```

_Plot coefficients for conjunctive position cells_

5c. Split by positive slope _PS cells_
```{r}
(pos_PS_coef_plot <- spatial_firing_coefs %>%
  filter(lm_group_b == "Positive", group == "PS") %>%
    standard_plot())

if (save_figures == 1) {
  ggsave(file = "plots/CoefficientValues_PositivePScells.png", width = 2.5, height = 2.5)
}

one.way_p_ps <- aov(coef ~ coef_type, data = filter(spatial_firing_coefs, lm_group_b == "Positive", group == "PS"))
summary(one.way_p_ps)

tukey.test_n <- TukeyHSD(one.way_p_ps)
tukey.test_n$coef_type
```

5d. Split by negative slope _PS cells_
```{r}
(neg_PS_coef_plot <- spatial_firing_coefs %>%
  filter(lm_group_b == "Negative", group == "PS") %>%
    standard_plot())

if (save_figures == 1) {
  ggsave(file = "plots/CoefficientValues_NegativePScells.png", width = 2.5, height = 2.5)
}

one.way_n_ps <- aov(coef ~ coef_type, data = filter(spatial_firing_coefs, lm_group_b == "Negative", group == "PS"))
summary(one.way_n_ps)

tukey.test_n <- TukeyHSD(one.way_n_ps)
tukey.test_n$coef_type
```


5e. Split by positive slope _PSA cells_
```{r}
(pos_PSA_coef_plot <- spatial_firing_coefs %>%
  filter(lm_group_b == "Positive", group == "PSA") %>%
    standard_plot())

if (save_figures == 1) {
  ggsave(file = "plots/CoefficientValues_PositivePSAcells.png", width = 2.5, height = 2.5)
}

one.way_p_psa <- aov(coef ~ coef_type, data = filter(spatial_firing_coefs, lm_group_b == "Positive", group == "PSA"))
summary(one.way_p_psa)

tukey.test_p_psa <- TukeyHSD(one.way_p_psa)
tukey.test_p_psa$coef_type
```

5f. Split by negative slope _PSA cells_
```{r}
(neg_PSA_coef_plot <- spatial_firing_coefs %>%
  filter(lm_group_b == "Negative", group == "PSA") %>%
    standard_plot())

if (save_figures == 1) {
  ggsave(file = "plots/CoefficientValues_NegativePSAcells.png", width = 2.5, height = 2.5)
}

one.way_n_psa <- aov(coef ~ coef_type, data = filter(spatial_firing_coefs, lm_group_b == "Negative", group == "PSA"))
summary(one.way_n_psa)

tukey.test_n_psa <- TukeyHSD(one.way_n_psa)
tukey.test_n_psa$coef_type
```

_Plot coefficients for just position cells_

5g. Plot coefficients for cells classified in Figure 2 as having a positive slope and in Figure 3 as having significant fit coefficients for position only. _P cells_
```{r}
(pos_P_coef_plot <- spatial_firing_coefs %>%
  filter(lm_group_b == "Positive", group == "P") %>%
    standard_plot())

if (save_figures == 1) {
  ggsave(file = "plots/CoefficientValues_PositivePcells.png", width = 2.5, height = 2.5)
}

one.way_p_p <- aov(coef ~ coef_type, data = filter(spatial_firing_coefs, lm_group_b == "Positive", group == "P"))
summary(one.way_p_p)

tukey.test_p_p <- TukeyHSD(one.way_p_p)
tukey.test_p_p$coef_type
```


5h. Plot coefficients for cells classified in Figure 2 as having a negative slope and in Figure 3 as having significant fit coefficients for position only. _P cells_
```{r}
(neg_P_coef_plot <- spatial_firing_coefs %>%
  filter(lm_group_b == "Negative", group == "P") %>%
    standard_plot())

if (save_figures == 1) {
  ggsave(file = "plots/CoefficientValues_NegativePcells.png", width = 2.5, height = 2.5)
}

one.way_n_p <- aov(coef ~ coef_type, data = filter(spatial_firing_coefs, lm_group_b == "Negative", group == "P"))
summary(one.way_n_p)

tukey.test_n_p <- TukeyHSD(one.way_n_p)
tukey.test_n_p$coef_type
```


_Plot coefficients for just speed cells_
5i. Plot coefficients for cells classified in Figure 1 as having a positive slope and in Figure 2 as having significant fit coefficients for position only. _SA cells_

```{r}
(pos_SA_coef_plot <- spatial_firing_coefs %>%
  filter(lm_group_b == "Positive", group == "SA") %>%
    standard_plot())

if (save_figures == 1) {
  ggsave(file = "plots/CoefficientValues_PositiveSAcells.png", width = 2.5, height = 2.5)
}
```


5j. Plot coefficients for cells classified in Figure 1 as having a negative slope and in Figure 2 as having significant fit coefficients for position only. _SA cells_

```{r}
(neg_SA_coef_plot <- spatial_firing_coefs %>%
  filter(lm_group_b == "Negative", group == "SA") %>%
    standard_plot())

if (save_figures == 1) {
  ggsave(file = "plots/CoefficientValues_NegativeSAcells.png", width = 2.5, height = 2.5)
}
```


5k. Plot coefficients for cells classified in Figure 1 as being unclassified but that have 'significant' positional coefficients in the mixed effect model.

```{r}
(UC_PSA_coef_plot <- spatial_firing_coefs %>%
  filter(lm_group_b == "Unclassified", group == "P" | group == "PS" | group == "PA" | group == "PSA") %>%
    standard_plot())


if (save_figures == 1) {
  ggsave(file = "plots/CoefficientValues_UCPSAcells.png", width = 2.5, height = 2.5)
}
```





```{r}
# Save results so code can be run later without repeating model fitting.
if (save_results == 1) {
  saveRDS(spatial_firing, "data_in/SpatialFiring_with_Results.Rda")
}
```