Dataset_Stats.Rmd

---
title: "DataSet_Stats"
output: html_document
---

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

Returns basic features of the data reported at the start of the results section.

### Check identity and number of mice
```{r}
# Cohorts
unique(as.character(spatial_firing$cohort))

# Mouse numbers within each cohort
unique(spatial_firing$Mouse)

# Miceidentified by cohort and number
unique(paste(spatial_firing$Mouse, as.character(spatial_firing$cohort)))
# or
unique(spatial_firing$unique_mouse)

# Number of mice
length(unique(paste(spatial_firing$Mouse, as.character(spatial_firing$cohort))))
```





### Check the number of sessions and the number per mouse
```{r}
# Number of sessions
length(unique(spatial_firing$session_id))

# Numbers of neurons per session
(sessions <- spatial_firing  %>%
  group_by(unique_mouse, session_id, .drop=FALSE) %>%
  summarise(n = n())) 

# Number of sessions per mouse
table(sessions$unique_mouse)

# Mean
mean(table(sessions$unique_mouse))
# Range
range(table(sessions$unique_mouse))


```


### Check the number of neurons and the number per mouse

```{r}
# Total number of neurons
dim(spatial_firing)[[1]]

# Show numbers of neurons per mouse
(neuron_numbers <- spatial_firing  %>%
  select(Mouse, cohort) %>%
  group_by(Mouse, cohort, .drop=FALSE) %>%
  summarise(n = n())) 
  
# Numbers of neurons per mouse
spatial_firing  %>%
  group_by(unique_mouse, .drop=FALSE) %>%
  summarise(n = n())

# mean
mean(neuron_numbers$n)

#range
range(neuron_numbers$n)
```


### Number of trials per session
```{r}
# Trial numbers for all sessions
trials_numbers <- spatial_firing  %>%
    select(session_id, cohort, max_trial_number) %>%
  mutate(max_trial_number = unlist(max_trial_number)) %>%
  distinct(session_id, .keep_all = TRUE)

# Mean
mean(trials_numbers$max_trial_number)
# SEM
sqrt(var(trials_numbers$max_trial_number)/length(trials_numbers$max_trial_number))
# Range
range(trials_numbers$max_trial_number)
```




### ----------------------------------------------------------------------
### Code below here should be deleted or incorporated into specific figures



First we might want to know how many sessions for each mouse

```{r}
position_ramps <- spatial_firing %>%
  filter(cue_group_o == "Positive-PI" | cue_group_o == "Negative-PI")  %>% 
  filter(final_model_o_b == "P" | final_model_o_b == "PS" | final_model_o_b == "PA" | final_model_o_b == "PSA")
  
```

```{r}
data <- spatial_firing  %>%
  select(Mouse, Day, cohort, session_id) %>%
  filter(Mouse != "1124") %>%
  dplyr::group_by(Mouse, cohort, .drop=FALSE) %>%
  dplyr::summarise(n = n()) 
  
data <- data[!duplicated(data$session_id), ]

data <- data  %>%
  select(session_id)
```

```{r}
write_csv2(data, "mouse_session_data_harry.csv")

```

```{r}
data <- position_ramps  %>%
  select(Mouse, cohort, session_id) %>%
  dplyr::group_by(Mouse, cohort, .drop=FALSE) %>%
  dplyr::summarise(n = n()) %>%
  dplyr::mutate(proportion = round(n / sum(n)*100, 2)) 
  #group_by(Mouse, Day) 
```


```{r}
data <- spatial_firing  %>%
  select(Mouse, cohort, session_id, track_category, track_category_p) %>%
  filter(Mouse != "1124") %>%
  dplyr::group_by(Mouse, cohort, track_category, track_category_p, .drop=FALSE) %>%
  dplyr::summarise(n = n()) %>%
  dplyr::group_by(track_category, track_category_p, .drop=FALSE) %>%
  dplyr::mutate(proportion = round(n / sum(n)*100, 2)) 
  #group_by(Mouse, Day) 

data$track_category <- as.character(data$track_category)
data$track_category_p <- as.character(data$track_category_p)

```

```{r}
write_csv2(data, "position_session_data_harry.csv")

```


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

### How many neurons passed criteria in the linear model ? - either outbound or homebound


```{r}
unclassified_neurons_total <- subset(spatial_firing, lm_group_b != "Unclassified" )
unclassified_neurons <- subset(unclassified_neurons_total, lm_group_b_h == "Unclassified")

```

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


### How many neurons passed criteria in the linear model ? 

1. get numbers of cells for each lm group (positive/negative/unclassified)
```{r}
start <- nrow(subset(spatial_firing, lm_group_b == "Negative" | lm_group_b == "Positive"))
reward <- nrow(subset(spatial_firing, lm_group_b == "Positive"))
nonslope <- nrow(subset(spatial_firing, lm_group_b == "Unclassified"))
```

2. find and plot proportions for ramp types in dataset
```{r}
start <- nrow(subset(spatial_firing, lm_group_b == "Negative"))/nrow(spatial_firing)*100
reward <- nrow(subset(spatial_firing, lm_group_b == "Positive"))/nrow(spatial_firing)*100
nonslope <- nrow(subset(spatial_firing, lm_group_b == "Unclassified"))/nrow(spatial_firing)*100
```



# How many numbers on average for each animal?

```{r}
data <- spatial_firing %>%
  dplyr::select(track_category, Mouse, cohort, Day) %>%
  filter(Mouse != "1124") %>%
  dplyr::group_by(Mouse, cohort, track_category, .drop=FALSE) %>%
  dplyr::summarise(n = n()) %>%
  dplyr::mutate(proportion = round(n / sum(n)*100, 2))  %>%
  filter(track_category != "None" & track_category != "NAs") %>%
  dplyr::group_by(Mouse, cohort, .drop=FALSE) %>%
  dplyr::summarise(avg = sum(proportion), sd = std.error(proportion))() %>%
  dplyr::ungroup() %>%
  summarise(mean = mean(avg), sd = std.error(avg), min=min(avg), max=max(avg))

data <- data %>%
  select(proportion,lm_group_b, cohort, Mouse) %>%
  filter(lm_group_b != "Unclassified" | lm_group_b != "NaN") %>%
  summarise(avg = mean(proportion), sd = std.error(proportion), min=min(proportion), max=max(proportion))


data$lm_group_b <- as.character(data$lm_group_b)
```


```{r}
write_csv2(data, "mouse_session_data_proportions.csv")

```

```{r}
data <- spatial_firing %>%
  filter(cue_group_o == "Positive-PI" | cue_group_o == "Negative-PI") %>%
  select(track_category,Mouse, cohort) %>%
  group_by(Mouse, track_category, cohort) %>%
  summarise(n = n()) %>%
  summarise(proportion = n / sum(n)) %>%
  filter(Mouse != "1124") %>%
  mutate(lmgroup = c("Negative", "Unclassified", "Positive")) %>%
  select(proportion,lmgroup) %>%
  group_by(lmgroup) %>%
  summarise(avg = mean(proportion), sd = std.error(proportion), min=min(proportion), max=max(proportion))
```


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

### How much of the shuffled dataset is past criteria? 
1. Extract shuffled slopes and rsquared values. 
```{r}
shuff_slopes <- tibble(slopes = unlist(spatial_firing$shuffle_results_b_o_slope), 
                       r2 = unlist(spatial_firing$shuffle_results_b_o_r.squared), 
                       pval = unlist(spatial_firing$shuffle_results_b_o_p.value), 
                       min_slope = rep(spatial_firing$shuffle_min_slope_b_o, times = 999), 
                       max_slope = rep(spatial_firing$shuffle_max_slope_b_o, times = 999))
```

2. Function to classify shuffled cells based on shuffled distribution
```{r}
shuff_slopes <- shuff_slopes %>%
  mutate(
    shuff_lm_group_b = pmap(
      list(
        min_slope,
        max_slope,
        slopes,
        pval
      ),
      compare_slopes
    )
  )
```

3. Calculate proportion of cells in the shuffled datasets that pass criteria
```{r}
total_shuffles <- nrow(shuff_slopes)
shuff_ramps <- nrow(subset(shuff_slopes,shuff_lm_group_b == "Positive" | shuff_lm_group_b == "Negative" ))
non_shuff_ramps <- nrow(subset(shuff_slopes,shuff_lm_group_b == "Unclassified" | shuff_lm_group_b == "NA" ))
non_shuff_ramps_percentage <- nrow(subset(shuff_slopes,shuff_lm_group_b == "Unclassified" | shuff_lm_group_b == "NA"))/nrow(shuff_slopes)*100

```


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

### correlate linear model result (i.e. positive/negative/unclassified) with ramp score

1. Extract ramp score 
_since its a list of three (outbound/homebound/all) we extract the first one_
```{r}
extract_ramp_score <- function(df){
  dx <- df[[1]]
  return(dx)
}

spatial_firing <- spatial_firing %>%
  mutate(start_ramp_score = map(ramp_score, extract_ramp_score))

```

2. find apsolute ramp scores for -ve and +ve slopes
```{r}
positive_ramps <- subset(spatial_firing, lm_group_b == "Positive")
negative_ramps <- subset(spatial_firing, lm_group_b == "Negative")
unclassified_ramps <- subset(spatial_firing, lm_group_b == "Unclassified")
```

3. find apsolute ramp scores for -ve and +ve slopes
```{r}
positive_ramps_scores <- tibble(ramp_score = as.numeric(positive_ramps$start_ramp_score))
negative_ramps_scores <- tibble(ramp_score = as.numeric(negative_ramps$start_ramp_score))
unclassied_ramps_scores <- tibble(ramp_score = as.numeric(unclassified_ramps$start_ramp_score))
```

4. Find mean of the ramp scores
```{r}
abs_pos_ramp <- positive_ramps_scores  %>%
  summarise(mean(ramp_score, na.rm = TRUE), sd(ramp_score, na.rm = TRUE), min(ramp_score, na.rm = TRUE), max(ramp_score, na.rm = TRUE), range(ramp_score, na.rm = TRUE))
abs_neg_ramp <- negative_ramps_scores  %>%
  summarise(mean(ramp_score, na.rm = TRUE), sd(ramp_score, na.rm = TRUE), min(ramp_score, na.rm = TRUE), max(ramp_score, na.rm = TRUE), range(ramp_score, na.rm = TRUE))
abs_un_ramp <- unclassied_ramps_scores  %>%
  summarise(mean(ramp_score, na.rm = TRUE), sd(ramp_score, na.rm = TRUE), min(ramp_score, na.rm = TRUE), max(ramp_score, na.rm = TRUE), range(ramp_score, na.rm = TRUE)) 
```


## perform ttest on ramp scores 
_compare to unclassified cells_

1. negative ramps 
```{r}
# paired t-test
t.test(positive_ramps_scores,unclassied_ramps_scores) # where y1 & y2 are numeric

```
2. positive ramps
```{r}
# paired t-test
t.test(negative_ramps_scores,unclassied_ramps_scores) # where y1 & y2 are numeric

```

_compare to shuffled data_
```{r}
positive_ramps_scores <- as.numeric(positive_ramps$start_ramp_score)
shuff_positive_ramps_scores <- as.numeric(shuff_positive_ramps$start_ramp_score)

# paired t-test for positive ramps
t.test(positive_ramps_scores,shuff_positive_ramps_scores) # where y1 & y2 are numeric


negative_ramps_scores <- as.numeric(negative_ramps$start_ramp_score)
shuff_negative_ramps_scores <- as.numeric(shuff_negative_ramps$start_ramp_score)

# paired t-test for negative ramps
t.test(negative_ramps_scores,shuff_negative_ramps_scores) # where y1 & y2 are numeric
```



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

_Figure 2 stats_

### How many neurons passed criteria in the linear mixed effect ? 

1. Split data by linear model result
```{r}
pos <-subset(spatial_firing, lm_group_b == "Positive")
neg <-subset(spatial_firing, lm_group_b == "Negative")
none <-subset(spatial_firing, lm_group_b == "None")
```

2. Then for each linear model classification (Positive/Negative) find the number of cells that encode position, speed etc (P,A,PS etc)
Do fist for positive slope neurons
```{r}
pos_P_num <- nrow(subset(pos, final_model_o_b == "P"))
pos_PA_num <- nrow(subset(pos, final_model_o_b == "PA"))
pos_PS_num <- nrow(subset(pos, final_model_o_b == "PS"))
pos_PSA_num <- nrow(subset(pos, final_model_o_b == "PSA"))
pos_SA_num <- nrow(subset(pos, final_model_o_b == "SA"))
pos_S_num <- nrow(subset(pos, final_model_o_b == "S"))
pos_A_num <- nrow(subset(pos, final_model_o_b == "A"))
```

3. same as above for negative slope neurons
```{r}
neg_P_num <- nrow(subset(neg, final_model_o_b == "P"))
neg_PA_num <- nrow(subset(neg, final_model_o_b == "PA"))
neg_PS_num <- nrow(subset(neg, final_model_o_b == "PS"))
neg_PSA_num <- nrow(subset(neg, final_model_o_b == "PSA"))
neg_SA_num <- nrow(subset(neg, final_model_o_b == "SA"))
neg_S_num <- nrow(subset(neg, final_model_o_b == "S"))
neg_A_num <- nrow(subset(neg, final_model_o_b == "A"))
```


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

_Figure 4B-E stats_
position encoding ramp neurons were classified with a similar slope on the non-beaconed and probe trials as on beaconed trials ....
1. extract position cells
```{r}

position_cells <- subset(spatial_firing, final_model_o_b == "P" | final_model_o_b == "PS" | final_model_o_b == "PA" | final_model_o_b == "PSA")

position_ramp_cells <- subset(position_cells, lm_group_b == "Negative" | lm_group_b == "Positive")

```

2. from above, extract position neurons that are ramp like in beaconed trials
```{r}
pos_position_cells <- subset(position_ramp_cells, lm_group_b == "Positive")
neg_position_cells <- subset(position_ramp_cells, lm_group_b == "Negative")
```

3. find neurons that maintain ramp polarity in the non beaconed trial
```{r}
posnb_position_cells <- nrow(subset(pos_position_cells, lm_group_nb == "Positive"))
negnb_position_cells <- nrow(subset(neg_position_cells, lm_group_nb == "Negative"))

```

4. calculate numbers
```{r}
total <- nrow(pos_position_cells) + nrow(neg_position_cells)
neurons_same <- (posnb_position_cells + negnb_position_cells)
```

5. the same as above but for neurons that switch ramp polarity in the non beaconed trial
```{r}
posnb_position_cells <- nrow(subset(pos_position_cells, lm_group_nb == "Unclassified"))
negnb_position_cells <- nrow(subset(neg_position_cells, lm_group_nb == "Unclassified"))
```

```{r}
total <- nrow(pos_position_cells) + nrow(neg_position_cells)
neurons_diff <- (posnb_position_cells + negnb_position_cells)
```

```{r}
un_position_cells <- subset(spatial_firing, lm_group_b == "Unclassified")
unnb_position_cells <- nrow(subset(un_position_cells, lm_group_nb == "Positive" | lm_group_nb == "Negative"))

```




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

_Figure 4C stats_
compare coefficients for position, speed and acceleration from lmer model
_focus on all position neurons (Figure 3) that have positive or negative slope (Figure 2)_

```{r}
position_cells <- subset(spatial_firing, final_model_o_b == "P" | final_model_o_b == "PS" | final_model_o_b == "PA" | final_model_o_b == "PSA" )

position_ramp_cells <- subset(position_cells, lm_group_b == "Negative" | lm_group_b == "Positive")

```

2. take only columns we need for the rest of the analysis
```{r}
spatial_firing_save <- select(position_ramp_cells, 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()

```

1. put data into a tibble for plotting
```{r}
data_coef <- tibble(session_id = rep(unlist(spatial_firing_save$unique_id), times=3), 
                    cluster_id = rep(spatial_firing_save$cluster_id, times=3), 
                    coef = c(as.character(spatial_firing_save$o_b_mod_coefs_pos),
                             as.character(spatial_firing_save$o_b_mod_coefs_speed),
                             as.character(spatial_firing_save$o_b_mod_coefs_accel)),
                    group = c(as.character(spatial_firing_save$final_model_o_b),
                              as.character(spatial_firing_save$final_model_o_b),
                              as.character(spatial_firing_save$final_model_o_b)), 
                    coef_type = rep(c("P", "S", "A"),
                                    each = nrow(spatial_firing_save)),
                    lm_result = rep(as.character(spatial_firing_save$lm_group_b), times=3))
```

2. run anova
```{r}
one.way <- aov(coef ~ coef_type, data = data_coef)

summary(one.way)
```

3. Post hoc Tukey's test
```{r}
tukey.test <- TukeyHSD(one.way)
tukey.test$coef_type
```


same as above but split by positive or negative slope
_focus on all position neurons (Figure 3) that have positive or negative slopes (Figure2)_

```{r}
position_cells <- subset(spatial_firing, final_model_o_b == "P" | final_model_o_b == "PS" | final_model_o_b == "PA" | final_model_o_b == "PSA")
```

```{r}
pos_position_cells <- subset(spatial_firing, lm_group_b == "Positive")
neg_position_cells <- subset(spatial_firing, lm_group_b == "Negative")

```

2. take only columns we need for the rest of the analysis
```{r}
spatial_firing_save <- select(neg_position_cells, 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()

```

1. put data into a tibble for plotting
```{r}
data_coef <- tibble(session_id = rep(unlist(spatial_firing_save$unique_id), times=3), 
                    cluster_id = rep(spatial_firing_save$cluster_id, times=3), 
                    coef = c(as.character(spatial_firing_save$o_b_mod_coefs_pos),
                             as.character(spatial_firing_save$o_b_mod_coefs_speed),
                             as.character(spatial_firing_save$o_b_mod_coefs_accel)),
                    group = c(as.character(spatial_firing_save$final_model_o_b),
                              as.character(spatial_firing_save$final_model_o_b),
                              as.character(spatial_firing_save$final_model_o_b)), 
                    coef_type = rep(c("P", "S", "A"),
                                    each = nrow(spatial_firing_save)),
                    lm_result = rep(as.character(spatial_firing_save$lm_group_b), times=3))
```

2. run anova
```{r}
one.way <- aov(coef ~ coef_type, data = data_coef)

summary(one.way)
```

3. Post hoc Tukey's test
```{r}
tukey.test <- TukeyHSD(one.way)
tukey.test$coef_type
```


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

## Calculate number of neurons that areset/switch/continuous across outbound and homebound

_Figure 5_

```{r}
position_cells <- subset(spatial_firing, final_model_o_b == "P" | final_model_o_b == "PS" | final_model_o_b == "PA" | final_model_o_b == "PSA")
```

```{r}
pi_data <- subset(position_cells, cue_group_o == "Positive-PI" | cue_group_o == "Negative-PI")
```


```{r}
speed_cells <- subset(spatial_firing, final_model_o_b == "S" | final_model_o_b == "SA" | final_model_o_b == "A")
```


```{r}
pi_data <- subset(speed_cells, lm_group_b == "Positive" | lm_group_b == "Negative")
```


```{r}

pospos <-subset(pi_data, lm_group_b == "Positive" & lm_group_b_h == "Positive")
negneg <-subset(pi_data, lm_group_b == "Negative" & lm_group_b_h == "Negative")

posneg <-subset(pi_data, lm_group_b == "Positive" & lm_group_b_h == "Negative")
negpos <-subset(pi_data, lm_group_b == "Negative" & lm_group_b_h == "Positive")

posnon <-subset(pi_data, lm_group_b == "Positive" & lm_group_b_h == "Unclassified")
negnon <-subset(pi_data, lm_group_b == "Negative" & lm_group_b_h == "Unclassified")


pospos_r <-subset(pi_data, lm_group_b == "Positive" & lm_group_b_h == "Positive" & reset_group == "Reset")
pospos_c <-subset(pi_data, lm_group_b == "Positive" & lm_group_b_h == "Positive" & reset_group == "Continuous")
negneg_r <-subset(pi_data, lm_group_b == "Negative" & lm_group_b_h == "Negative" & reset_group == "Reset")
negneg_c <-subset(pi_data, lm_group_b == "Negative" & lm_group_b_h == "Negative" & reset_group == "Continuous")

```

_Figure 5_

```{r}
switch_slope <- nrow(posneg) + nrow(negpos)
only_outbound_slope <- nrow(posnon) + nrow(negnon)
same_slope <- nrow(pospos) + nrow(negneg)
same_slope_reset <- nrow(pospos_r) + nrow(negneg_r)
same_slope_continuous <- nrow(pospos_c) + nrow(negneg_c)

```


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

## Figure 7 stats

_Figure 7_

```{r}
ramps <- subset(spatial_firing, lm_group_b == "Negative" | lm_group_b == "Positive")
ramps_with_enough_trials <- subset(ramps, run_trials_number >= 3 & try_trials_number >= 3)
```

```{r}
sig_ramps <- nrow(subset(ramps, reward_interaction_id == "Significant"))/nrow(ramps_with_enough_trials)*100
nonsig_ramps <- nrow(subset(ramps, reward_interaction_id == "Not-Significant"))/nrow(ramps_with_enough_trials)*100
```
```{r}
sig_ramps <- nrow(subset(ramps_with_enough_trials, reward_interaction_id == "Significant"))
nonsig_ramps <- nrow(subset(ramps_with_enough_trials, reward_interaction_id == "Not-Significant"))
```



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

## calculate total number of trials for each trial outcome 
## Figure 7 stats

```{r}
trial_adder <- function(trials) {
  x<-length(trials)
  #print(x)
  return(x)
}

```


```{r}
ramps <- spatial_firing %>%
  select(run_through_trialid, try_trialid,rewarded_locations) %>%
  mutate(run_trials_number = map(run_through_trialid,trial_adder),
         try_trials_number = map(try_trialid,trial_adder),
         hit_trials_number = map(rewarded_locations,trial_adder))
```

```{r}
ramps_sum <- ramps %>%
  #unnest(c(run_trials_number,try_trials_number,hit_trials_number))
  summarise(total_run = sum(unlist(run_trials_number)), total_try = sum(unlist(try_trials_number)), total_hit = sum(unlist(hit_trials_number)))


```