natfeedback_code_eLife/R/figure_1_S2.Rmd

---
title: "Spacek et al., 2021, Figure 1-Supplement 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/fig1.csv")
```

```{r tidy_for_1_S2ab, include = FALSE}

# Turn booleans for 'optogenetic manipulation' into a binary predictor
tb <- tib %>% mutate(feedback = ifelse(opto == TRUE, 0, 1))

# Signal-to-noise ratio and mean peak-width are not computed on a trial-by-trial basis. In the csv-file, 
# these measures are therefore identical across trials, so we simply pull out the first trial of each neuron
tb = tb %>% select(mid, sid, eid, uid, mseu, feedback, snr, meanpkw) %>% distinct(mseu, feedback, .keep_all = TRUE)
```

# Figure 1-Supplement 2a
## Average effect of feedback on signal-to-noise ratio (SNR)

```{r, fit_model_1_S2a}

# We fit a random-intercept model with two random effects: 
# (1) Neurons (uid) can have different baseline firing rates
# (2) Mean firing rates are allowed to differ across recording sessions (sid)
# More complex models with random slopes for neurons (or with experiments nested in
# sessions, nested in mice) give singular fits.
lmer.1_S2a = lmer(snr ~ feedback + (1 | uid) + (1 | sid), 
                  data = tb %>% drop_na(snr))

display(lmer.1_S2a)
anova(lmer.1_S2a)
```

```{r get_predicted_average_effect1_S2a, include=F}

mSuppr = fixef(lmer.1_S2a)[1]
diffMeans = fixef(lmer.1_S2a)[2]
mActive  = fixef(lmer.1_S2a)[1] + diffMeans
```

Feedback SNR: `r format(mActive, digits=2, nsmall=2)` \newline
Suppression SNR: `r format(mSuppr, digits=2, nsmall=2)` \newline
n = `r nrow(tb %>% drop_na(snr) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(snr) %>% count(mid))` mice

\newpage
# Figure 1-Supplement 2b
## Average effect of feedback on PSTH mean peak width

```{r, fit_model_1_S2b}

# Random-intercept for single neurons, 
# random intercept for experiments, nested in series
lmer.1_S2b = lmer(meanpkw ~ feedback + (1 | uid) + (1 | sid/eid), 
                  data = tb %>% drop_na(meanpkw))

display(lmer.1_S2b)
anova(lmer.1_S2b)
```

```{r get_predicted_average_effect_1_S2b, include=F}

mSuppr = fixef(lmer.1_S2b)[1]
diffMeans = fixef(lmer.1_S2b)[2]
mActive  = fixef(lmer.1_S2b)[1] + diffMeans
```

Feedback mean peak width: `r format(mActive, digits=2, nsmall=2)` \newline
Suppression mean peak width: `r format(mSuppr, digits=2, nsmall=2)` \newline
n = `r nrow(tb %>% drop_na(meanpkw) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(meanpkw) %>% count(mid))` mice

\newpage
# Figure 1-Supplement 2c
## Relation between firing rate FMI and burst ratio FMI

```{r read_data_1_S2c-g, include=FALSE}

tib = get_data("../csv/mviFMI.csv")

# filter based on 'state'
tb <- tib %>% filter(st8 == 'none')
```

```{r fit_model_1_S2c}

# Random-intercept for single neurons, 
# random intercept for experiments, nested in series
lmer.1_S2c = lmer(meanburstratio ~ meanrate + (1 | uid) + (1 | sid/eid), 
                  data = tb %>% drop_na(meanburstratio, meanrate))

display(lmer.1_S2c)
anova(lmer.1_S2c)
```

```{r store_coefficients_1_S2c, include=FALSE}

coef_df = data.frame("intercept" = fixef(lmer.1_S2c)[1], "slope" = fixef(lmer.1_S2c)[2], row.names = "")
write_csv(coef_df, "_stats/figure_1_S2c_coefs.csv")
```

Slope of `r format(fixef(lmer.1_S2c)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.1_S2c)[2], digits=2, nsmall=2)` (95\%-confidence interval) \newline
n = `r nrow(tb %>% drop_na(meanburstratio, meanrate) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(meanburstratio, meanrate) %>% count(mid))` mice

\newpage
# Figure 1-Supplement 2d
## Relation between firing rate FMI and sparseness FMI

```{r fit_model_1_S2d}

# Random-intercept for single neurons, 
# random intercept for experiments, nested in series
lmer.1_S2d = lmer(spars ~ meanrate + (1 | uid) + (1 | sid/eid), 
                  data = tb %>% drop_na(spars, meanrate))

display(lmer.1_S2d)
anova(lmer.1_S2d)
```

```{r store_coefficients_1_S2d, include=FALSE}

coef_df = data.frame("intercept" = fixef(lmer.1_S2d)[1], "slope" = fixef(lmer.1_S2d)[2], row.names = "")
write_csv(coef_df, "_stats/figure_1_S2d_coefs.csv")
```

Slope of `r format(fixef(lmer.1_S2d)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.1_S2d)[2], digits=2, nsmall=2)` (95\%-confidence interval) \newline
n = `r nrow(tb %>% drop_na(spars, meanrate) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(spars, meanrate) %>% count(mid))` mice

\newpage
# Figure 1-Supplement 2e
## Relation between firing rate FMI and reliability FMI

```{r fit_model_1_S2e}

# Random-intercept for single neurons, 
# random intercept for series, nested in mice
lmer.1_S2e = lmer(rel ~ meanrate + (1 | uid) + (1 | mid/sid), 
                  data = tb %>% drop_na(rel, meanrate))

display(lmer.1_S2e)
anova(lmer.1_S2e)
```

```{r store_coefficients_1_S2e, include=FALSE}

coef_df = data.frame("intercept" = fixef(lmer.1_S2e)[1], "slope" = fixef(lmer.1_S2e)[2], row.names = "")
write_csv(coef_df, "_stats/figure_1_S2e_coefs.csv")
```

Slope of `r format(fixef(lmer.1_S2e)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.1_S2e)[2], digits=2, nsmall=2)` (95\%-confidence interval) \newline
n = `r nrow(tb %>% drop_na(rel, meanrate) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(rel, meanrate) %>% count(mid))` mice

\newpage
# Figure 1-Supplement 2f
## Relation between firing rate FMI and SNR FMI

```{r fit_model_1_S2f}

# Random intercept for neurons, 
# random intercept for series
lmer.1_S2f = lmer(snr ~ meanrate + (1 | uid) + (1 | sid), 
                  data = tb %>% drop_na(snr, meanrate))

display(lmer.1_S2f)
anova(lmer.1_S2f)
```

```{r store_coefficients_1_S2f, include=FALSE}

coef_df = data.frame("intercept" = fixef(lmer.1_S2f)[1], "slope" = fixef(lmer.1_S2f)[2], row.names = "")
write_csv(coef_df, "_stats/figure_1_S2f_coefs.csv")
```

Slope of `r format(fixef(lmer.1_S2f)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.1_S2f)[2], digits=2, nsmall=2)` (95\%-confidence interval) \newline
n = `r nrow(tb %>% drop_na(snr, meanrate) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(snr, meanrate) %>% count(mid))` mice

\newpage
# Figure 1-Supplement 2g
## Relation between firing rate FMI and mean peak witdth FMI

```{r fit_model_1_S2g}

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

display(lmer.1_S2g)
anova(lmer.1_S2g)
```

```{r store_coefficients_1_S2g, include=FALSE}

coef_df = data.frame("intercept" = fixef(lmer.1_S2g)[1], "slope" = fixef(lmer.1_S2g)[2], row.names = "")
write_csv(coef_df, "_stats/figure_1_S2g_coefs.csv")
```

Slope of `r format(fixef(lmer.1_S2g)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.1_S2g)[2], digits=2, nsmall=2)` (95\%-confidence interval) \newline
n = `r nrow(tb %>% drop_na(meanpkw, meanrate) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(meanpkw, meanrate) %>% count(mid))` mice

\newpage
# Figure 1-Supplement 2h
## Effect of feedback on eye position variability

```{r read_data_1_S2h, include=FALSE}

tib = get_data("../csv/ipos_opto.csv")

# Turn feedback into a binary variable
tb = tib %>% mutate(feedback = ifelse(opto == "FALSE", 1, 0))

# The standard deviation of eye position is not computed on a trial-by-trial basis. In the csv-file, 
# std is therefore identical across trials, so we simply pull out the first trial of each neuron
tb = tb %>% select(mid, sid, eid, mse, feedback, std_xpos_cross) %>% distinct(mse, feedback, .keep_all = TRUE)

```

```{r fit_model_1_S2h}

# Random intercept for experiments, nested in series, nested in mice
lmer.8 = lmer(std_xpos_cross ~ feedback + (1 | mid/sid/eid), 
              data = tb %>% drop_na(std_xpos_cross))

display(lmer.8)
anova(lmer.8)

```

```{r get_predicted_average_effect_1_S2h, include=F}
mSuppr = fixef(lmer.8)[1]
diffMeans = fixef(lmer.8)[2]
mFeedback  = fixef(lmer.8)[1] + diffMeans

```

Mean eye position standard deviation with feedback: `r format(mFeedback, digits=2, nsmall=2)`$^{\circ}$ \newline
Mean eye position standard deviation with suppression: `r format(mSuppr, digits=2, nsmall=2)`$^{\circ}$  \newline
n = `r nrow(tb %>% drop_na(std_xpos_cross) %>% count(eid))` experiments from `r nrow(tb %>% drop_na(std_xpos_cross) %>% count(mid))` mice

\newpage
# Figure 1-Supplement 2i
## Relation between feedback effects on eye position and feedback effects on reliability

```{r read_data_1_S2i, include=FALSE}

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

```{r fit_model_1S2i}

# Random intercept for neurons,
# random intercept for experiments nested in series
lmer.1_S2i = lmer(relfmi ~ iposfmi + (1 | uid) + (1 | sid/eid), 
                  data = tib %>% drop_na(relfmi, iposfmi))

display(lmer.1_S2i)
anova(lmer.1_S2i)
```

```{r store_coefficients, include=FALSE}

coef_df = data.frame("intercept" = fixef(lmer.1_S2i)[1], "slope" = fixef(lmer.1_S2i)[2], row.names = "")
write_csv(coef_df, "_stats/figure_1_S2i_coefs.csv")
```

Slope of `r format(fixef(lmer.1_S2i)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.1_S2i)[2], digits=2, nsmall=2)` (95\%-confidence interval) \newline
n = `r nrow(tib %>% drop_na(relfmi, iposfmi) %>% count(uid))` neurons from `r nrow(tib %>% drop_na(relfmi, iposfmi) %>% count(mid))` mice