natfeedback_code_eLife/R/figure_2.Rmd

---
title: "Spacek et al., 2021, Figure 2"
output: pdf_document
---

```{r setup, include=FALSE}

knitr::opts_chunk$set(echo = TRUE)

library(arm)
library(lmerTest)
library(tidyverse)

source('get_data.R')
```

```{r read_data, include=FALSE}

tib = get_data("../csv/fig2.csv")
```

# Figure 2e

``` {r tidy_2e, include=FALSE}

tb_all = tib %>% filter(fmode == "all" )
```

## (I) Fit intercept-only model to slope (all spikes)

``` {r fit_model_2e_slope}

# Random intercept for neurons,
# random intercept for experiments, nested in series
lmer.2e_slope = lmer(slope ~ 1 + (1 | uid) + (1 | sid/eid), 
                     data = tb_all %>% drop_na(slope))

display(lmer.2e_slope)
```

```{r, include=FALSE}

print(paste("Intercept: ", format(fixef(lmer.2e_slope), digits=2), " [", format(fixef(lmer.2e_slope) - 2 * se.fixef(lmer.2e_slope), digits=2), ", ", format(fixef(lmer.2e_slope) + 2 * se.fixef(lmer.2e_slope), digits=2), "]", sep=""))
```
## 95-% confidence interval on slope

`r format(fixef(lmer.2e_slope)[1], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.2e_slope)[1], digits=2, nsmall=2)` \newline
n = `r nrow(tb_all %>% drop_na(slope) %>% count(uid))` neurons from `r nrow(tb_all %>% drop_na(slope) %>% count(mid))` mice

\newpage
## (II) Fit intercept-only model to threshold (all spikes)

```{r, fit_model_2e_threshold}

# Random intercept for neurons,
# random intercept for experiments, nested in mice
lmer.2e_thresh = lmer(thresh ~ 1 + (1 | uid) + (1 | mid/eid), 
                      data = tb_all %>% drop_na(thresh))

display(lmer.2e_thresh)
```

```{r, include=FALSE}

print(paste("Intercept: ", format(fixef(lmer.2e_thresh), digits=2), " [", format(fixef(lmer.2e_thresh) - 2 * se.fixef(lmer.2e_thresh) , digits=2), ", ", format(fixef(lmer.2e_thresh) + 2 * se.fixef(lmer.2e_thresh), digits=2), "]", sep=""))
```
## 95-% confidence interval on threshold

`r format(fixef(lmer.2e_thresh)[1], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.2e_thresh)[1], digits=2, nsmall=2)` \newline
n = `r nrow(tb_all %>% drop_na(thresh) %>% count(uid))` neurons from `r nrow(tb_all %>% drop_na(thresh) %>% count(mid))` mice

\newpage
# Figure 2f
## Goodness-of-fit vs burst ratio during suppression

``` {r, fit_model_2f}

# Random intercept for neurons,
# random intercept for experiments, nested in series
lmer.2f = lmer(rsq ~ suppression_meanburstratio + (1 | uid) + (1 | sid/eid), 
               data = tb_all %>% drop_na(rsq, suppression_meanburstratio))

display(lmer.2f)
anova(lmer.2f)
```

## 95-% confidence interval on slope

`r format(fixef(lmer.2f)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.2f)[2], digits=2, nsmall=2)` \newline
n = `r nrow(tb_all %>% drop_na(rsq, suppression_meanburstratio) %>% count(uid))` neurons from `r nrow(tb_all %>% drop_na(rsq, suppression_meanburstratio) %>% count(mid))` mice

```{r, store_coeffs_2f, include=F}

# store results in file
coef_df = data.frame("intercept" = fixef(lmer.2f)[1], "slope" = fixef(lmer.2f)[2], row.names = "")
write_csv(coef_df, "_stats/figure_2f_coefs.csv")
```

\newpage
# Figure 2g
## Goodness-of-fit with and without removal of burst spikes

``` {r tidy_2g, include = F}

# Extract two firing modes: all spikes, and non-burst spikes, and turn them into numeric binary predictors
tb <- tib %>% filter(fmode == "all" | fmode == 'nonburst') %>% mutate(allSpikes = ifelse(fmode == "all", 1, 0))
```

``` {r fit_model_2g}

# Random intercept for neurons,
# random intercept for experiments, nested in series
lmer.2g = lmer(rsq ~ allSpikes + (1 | uid) + (1 | sid/eid), 
               data = tb %>% drop_na(rsq, allSpikes))

display(lmer.2g)
anova(lmer.2g)
```


```{r get_predicted_average_effect_2g, include=F}

mNonBurst = fixef(lmer.2g)[1]
diffMeans = fixef(lmer.2g)[2]
mAll = fixef(lmer.2g)[1] + diffMeans
```

# Predicted average effect

All spikes: $R^2$ = `r format(mAll, digits=2, nsmall=2)` \newline
non-burst spikes: $R^2$ = `r format(mNonBurst, digits=2, nsmall=2)` \newline
n = `r nrow(tb %>% drop_na(rsq, allSpikes) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(rsq, allSpikes) %>% count(mid))` mice


\newpage
# Figure 2h
## Goodness of fit with and without removal of tonic spikes

``` {r tidy_2h, include=F}

# Extract two firing modes: all spikes, and non-random spikes
tb <- tib %>% filter(fmode == "all" | fmode == 'nonrand') %>% mutate(allSpikes = ifelse(fmode == "all", 1, 0))
```


``` {r fit_model_2h}

# Random intercept for neurons,
# random intercept for experiments, nested in series
lmer.2h = lmer(rsq ~ allSpikes + (1 | uid) + (1 | sid/eid), 
               data = tb %>% drop_na(rsq, allSpikes))

display(lmer.2h)
anova(lmer.2h)
```


```{r get_predicted_average_effect_2h, include=F}

mNonRand = fixef(lmer.2h)[1]
diffMeans = fixef(lmer.2h)[2]
mAllSpikes = fixef(lmer.2h)[1] + diffMeans
```

# Predicted average effect

All spikes: $R^2$ = `r format(mAllSpikes, digits=2, nsmall=2)` \newline
Randomly removed spikes: $R^2$ = `r format(mNonRand, digits=2, nsmall=2)` \newline
n = `r nrow(tb %>% drop_na(rsq, allSpikes) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(rsq, allSpikes) %>% count(mid))` mice

\newpage
# Figure 2i

``` {r tidy_2i, include=FALSE}

tb_nonburst = tib %>% filter(fmode == "nonburst")
```

## (I) Fit intercept-only model to slope (non-burst spikes)

```{r, fit_model_2i_slope}

# Random intercept for neurons,
# random intercepts for experiments, nested in series
lmer.2i_slope = lmer(slope ~ 1 + (1 | uid) + (1 | sid/eid), 
                     data = tb_nonburst %>% drop_na(slope))

display(lmer.2i_slope)
```

## 95-% confidence interval

```{r, include=FALSE}

print(paste("Intercept: ", format(fixef(lmer.2i_slope), digits=2), " [", format(fixef(lmer.2i_slope) - 2 * se.fixef(lmer.2i_slope), digits=2), ", ", format(fixef(lmer.2i_slope) + 2 * se.fixef(lmer.2i_slope), digits=2), "]", sep=""))
```

Intercept of `r format(fixef(lmer.2i_slope)[1], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.2i_slope)[1], digits=1, nsmall=1)` \newline
n = `r nrow(tb_nonburst %>% drop_na(slope) %>% count(uid))` neurons from `r nrow(tb_nonburst %>% drop_na(slope) %>% count(mid))` mice

\newpage
## (II) Fit intercept-only model to threshold (non-burst spikes)

```{r, fit_model_2i_threshold}

# Random intercept for neurons,
# random intercept for experiments, nested in series, nested in mice
lmer.2i_thresh = lmer(thresh ~ 1 + (1 | uid) + (1 | mid/sid/eid), 
                      data = tb_nonburst %>% drop_na(thresh))

display(lmer.2i_thresh)
```

## 95-% confidence interval

```{r, include=FALSE}

print(paste("Intercept: ", format(fixef(lmer.2i_thresh), digits=2), " [", format(fixef(lmer.2i_thresh) - 2 * se.fixef(lmer.2i_thresh), digits=2), ", ", format(fixef(lmer.2i_thresh) + 2 * se.fixef(lmer.2i_thresh), digits=2), "]", sep=""))
```

Intercept of `r format(fixef(lmer.2i_thresh)[1], digits=1, nsmall=1)` $\pm$ `r format(2 * se.fixef(lmer.2i_thresh)[1], digits=1, nsmall=1)`

\newpage
# Fit a single model to Figures 2e and 2i

## (I) Predict slope based on firing mode (all spikes vs non-burst spikes)

```{r, tidy_2ei, include=FALSE}

# (1) Fit a single model to panels e and i to examine changes in slope with firing mode as a binary predictor. We need to turn the strings (= factors, = categorical predictors) into binary predictors (numeric) first:
tb = tib %>% filter(fmode == "all" | fmode == "nonburst") %>% mutate(allSpikes = ifelse(fmode == "all", 1, 0))
```

```{r fit_model_2ei_slope}

# Random intercept for neurons,
# random intercept for experiments, nested in series, nested in mice
lmer.2ei_slope = lmer(slope ~ allSpikes + (1 | uid) + (1 | mid/sid/eid), 
                      data = tb %>% drop_na(slope, allSpikes))

display(lmer.2ei_slope)
anova(lmer.2ei_slope)
```

\newpage
## (II) Predict threshold based on firing mode (all spikes vs non-burst spikes)

```{r, fit_model_2ei_threshold}

# Random intercept for neurons,
# random intercept for experiments, nested in series, nested in mice
lmer.2ei_thresh = lmer(thresh ~ allSpikes + (1 | uid) + (1 | mid/sid/eid), 
                       data = tb %>% drop_na(thresh, allSpikes))

display(lmer.2ei_thresh)
anova(lmer.2ei_thresh)
```