natfeedback_code_eLife/R/figure_4_S1.Rmd

---
title: "Spacek et al., 2021, Figure 4-Supplement 1"
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_4_S1a, include=FALSE}

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

```{r tidy_for_4_S1, include = FALSE}

tb <- tib %>% filter(st8 == "none") %>% select(mid, sid, eid, uid, mseu, meanrate, meanburstratio)
```

# Figure 4-Supplement 1a
## Relation between firing rate FMI and burst ratio FMI

```{r fit_model_4_S1a}

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

display(lmer.4_S1a)
anova(lmer.4_S1a)
```

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

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

Slope of `r format(fixef(lmer.4_S1a)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.4_S1a)[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 4-Supplement 1b
## Relation between firing rates in response to movies vs gratings


```{r read_data_4_S1b, include=FALSE}

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

```{r tidy_for_4_S1b, include = FALSE}

# Turn stimulus type into binary predictor
tb <- tib %>% mutate(mvi = ifelse(stimtype == "mvi", 1, 0)) %>% select(mid, sid, uid, msu, mvi, meanrate)
```

```{r fit_model_4_S1b}

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

display(lmer.4_S1b)
anova(lmer.4_S1b)
```

```{r get_predicted_average_effect_4_S1b, include=F}

m_grt = fixef(lmer.4_S1b)[1]
diffMeans = fixef(lmer.4_S1b)[2]
m_mvi  = fixef(lmer.4_S1b)[1] + diffMeans
```

Movies: `r format(m_mvi, digits=2, nsmall=2)` spikes/s \newline
Gratings: `r format(m_grt, digits=2, nsmall=2)` spikes/s \newline
n = `r nrow(tb %>% drop_na(meanrate) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(meanrate) %>% count(mid))` mice

\newpage
# Figure 4-Supplement 1c,d
## Relation between firing rates during feedback vs suppression for separate epochs of the movie

```{r read_data_4_S1cd, include=FALSE}

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

```{r tidy_4_S1cd, include=FALSE}

# Get relevant variables, rename rate* columns, such that they have a common prefix, followed by a simple number
tmp1 = tib %>% select(mid, sid, eid, uid, mseu, trialis, opto, rate02s, rate35s) %>% drop_na(trialis, rate02s, rate35s) %>% rename(rate2 = rate02s, rate35 = rate35s)

# Turn rate* columns into long format
tmp2 = tmp1 %>% pivot_longer(cols = starts_with("rate"), names_to = "window", names_prefix = "rate", values_to = "trialrate")

# Turn booleans for 'optogenetic manipulation' and 'time window' into binary, numeric predictors
tmp2 = tmp2 %>% mutate(feedback = ifelse(opto == TRUE, 0, 1)) %>% mutate(earlywin = ifelse(window == '2', 1, 0))
```

```{r fit_model_4_S1cd}

lmer.4_S1cd = lmer(trialrate ~ feedback*earlywin + (1 + feedback*earlywin | uid) + (1 | mid/sid/eid),
                   data = tmp2 %>% drop_na(trialrate))

display(lmer.4_S1cd)
anova(lmer.4_S1cd)
```

n = `r nrow(tmp2 %>% drop_na(trialrate) %>% count(uid))` neurons from `r nrow(tmp2 %>% drop_na(trialrate) %>% count(mid))` mice

```{r post_hocs_4_S1cd, include=FALSE}

# Post-hoc: fit model to early window (0-2s)
foo = tmp2 %>% filter(earlywin == 1)

lmer.4_S1cd = lmer(trialrate ~ feedback + (1 + feedback | uid) + (1 | sid/eid),
                   data = foo %>% drop_na(trialrate))

display(lmer.4_S1cd)
anova(lmer.4_S1cd)

m_suppr = fixef(lmer.4_S1cd)[1]
diffMeans = fixef(lmer.4_S1cd)[2]
m_feedb  = fixef(lmer.4_S1cd)[1] + diffMeans

# Post-hoc: fit model to late window (3-5s)
foo = tmp2 %>% filter(earlywin == 0)

lmer.4_S1cd = lmer(trialrate ~ feedback + (1 + feedback | uid) + (1 | sid/eid),
                   data = foo %>% drop_na(trialrate))

display(lmer.4_S1cd)
anova(lmer.4_S1cd)

m_suppr = fixef(lmer.4_S1cd)[1]
diffMeans = fixef(lmer.4_S1cd)[2]
m_feedb  = fixef(lmer.4_S1cd)[1] + diffMeans
```

\newpage
# Figure 4-Supplement 1e
## Burst ratio FMIs separately for grating versus movie presentations

```{r read_data_4_S1ef, include=FALSE}

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

```{r tidy_for_4_S1e, include=FALSE}

t4S1e <- tib %>% filter(blank == "False" & meanburstratio != "NA")

# Turn stimulus type into binary predictor
t4S1e = t4S1e %>% mutate(mvi = ifelse(stimtype == "mvi", 1, 0))
```

```{r fit_model_4_S1e}

# Random intercept for neurons,
# random intercept for mice
lmer.4_S1e = lmer(meanburstratio ~ mvi + (1 | uid) + (1 | mid), 
                  data = t4S1e %>% drop_na(meanburstratio))

display(lmer.4_S1e)
anova(lmer.4_S1e)
```

```{r predicted_average_effect_4_S1e, include=F}

m_grt = fixef(lmer.4_S1e)[1]
diffMeans = fixef(lmer.4_S1e)[2]
m_mvi  = fixef(lmer.4_S1e)[1] + diffMeans
```

```{r save_coefficients_4_S1e, include=FALSE}

pred_means_df = data.frame("mvi" = m_mvi, "grt" = m_grt, row.names = "")
write_csv(pred_means_df, "_stats/figure_4_S1e_pred_means.csv")
```

FMI movies: `r format(m_mvi, digits=2, nsmall=2)` \newline
FMI gratings: `r format(m_grt, digits=2, nsmall=2)` \newline
n = `r nrow(t4S1e %>% drop_na(meanburstratio) %>% count(uid))` neurons from `r nrow(t4S1e %>% drop_na(meanburstratio) %>% count(mid))` mice 

\newpage
# Figure 4-Supplement 1f
## Burst ratio FMIs across blank screen conditions

```{r tidy_for_4_S1f, include=FALSE}

# Code the 'blank condition' as integers 1, 2, 3
t3f = tib %>% mutate(blank_condition = case_when(
  stimtype == 'mvi' & blank == 'prestim' ~ "mvi",
  stimtype == 'grt' & blank == 'prestim' ~ "grt",
  stimtype == 'grt' & blank == 'cond' ~ "grt0c"
))

# Remove NAs
t3f = t3f %>% filter(meanburstratio != "Na" & blank_condition != "Na")

# And turn them into a factor for dummy coding
t3f$blank_condition = as.factor(t3f$blank_condition)
```


```{r fit_model_4_S1f}

# Random intercept for neurons,
# random intercept for series
lmer.4_S1f = lmer(meanburstratio ~ blank_condition + (1 | uid) + (1 | sid), 
                  data = t3f %>% drop_na(meanburstratio))

display(lmer.4_S1f)
anova(lmer.4_S1f)
```

```{r get_predicted_average_effect_4_S1f, include=FALSE}

m_grt = fixef(lmer.4_S1f)[1]
m_grt0c = fixef(lmer.4_S1f)[1] + fixef(lmer.4_S1f)[2]
m_mvi = fixef(lmer.4_S1f)[1]+ fixef(lmer.4_S1f)[3]
```

```{r save_coefficients_4_S1f, include=FALSE}

pred_means_df = data.frame("grt" = m_grt, "grt0c" = m_grt0c, "mvi" = m_mvi, row.names = "")
write_csv(pred_means_df, "_stats/figure_4_S1f_pred_means.csv")
```

FMI pre-movie:`r format(m_mvi, digits=2, nsmall=2)` \newline
FMI pre-grating: `r format(m_grt, digits=2, nsmall=2)` \newline
FMI blank grating: `r format(m_grt0c, digits=2, nsmall=2)` \newline
n = `r nrow(t3f %>% drop_na(meanburstratio) %>% count(uid))` neurons from `r nrow(t3f %>% drop_na(meanburstratio) %>% count(mid))` mice 

\newpage
## Post-hoc analysis, comparing mean across blank conditions against stimulus condition

```{r tidy_for_post_hoc_4_S1ef, include = F}

t4S1ef = tib %>% filter(meanburstratio != "Na")

t4S1ef = t4S1ef %>% mutate(stim_condition = case_when(
  stimtype == 'mvi' & blank == 'False' ~ 1,
  stimtype == 'grt' & blank == 'False' ~ 2,
  (stimtype == 'mvi' | stimtype == 'grt') & blank != 'False' ~ 3
))

t4S1ef$stim_condition = as.factor(t4S1ef$stim_condition)
```

```{r fit_model_4_S1ef, echo=F}

# Random intercept for neuron,
# random intercept for series
lmer.4_S1ef = lmer(meanburstratio ~ stim_condition + (1 | uid) + (1 | sid), 
                   data = t4S1ef %>% drop_na(meanburstratio))

display(lmer.4_S1ef)
anova(lmer.4_S1ef)
```

```{r predicted_average_effect_4_S1ef, include=FALSE}

m_mvi = fixef(lmer.4_S1ef)[1]
m_grt = fixef(lmer.4_S1ef)[1] + fixef(lmer.4_S1ef)[2]
m_blank = fixef(lmer.4_S1ef)[1]+ fixef(lmer.4_S1ef)[3]
```


\newpage
# Figure 4-Supplement 1g
## Feedback effect on firing rate during blank periods preceding movies


```{r read_data_4_S1gj, include=FALSE}

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

```{r tidy_for_4_S1gj, include = FALSE}

# Turn stimulus type into binary predictor, get relevant columns
tb <- tib %>% mutate(feedback = ifelse(opto == "TRUE", 0, 1)) %>% select(mid, sid, eid, uid, mseu, feedback, everything(), -opto)
```

```{r fit_model_4_S1g}

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

display(lmer.4_S1g)
anova(lmer.4_S1g)
```

```{r predicted_average_effect_4_S1g, include=F}

m_suppress = fixef(lmer.4_S1g)[1]
diffMeans = fixef(lmer.4_S1g)[2]
m_feedback  = fixef(lmer.4_S1g)[1] + diffMeans
```

Feedback: `r format(m_feedback, digits=2, nsmall=2)` spikes/s. \newline
Suppression: `r format(m_suppress, digits=2, nsmall=2)` spikes/s \newline
n = `r nrow(tb %>% drop_na(blankrates) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(blankrates) %>% count(mid))` mice 

\newpage
# Figure 4-Supplement 1j
## Feedback effect on burst ratio during blank periods preceding movies

```{r fit_model_4_S1j}

# Random intercept, random slope for neurons,
# random intercept for experiments
lmer.4_S1j = lmer(blankburstratios ~ feedback + (1 + feedback | uid) + (1 | eid), 
                  data = tb %>% drop_na(blankburstratios))

display(lmer.4_S1j)
anova(lmer.4_S1j)
```

```{r predicted_average_effect_4_S1j, include=F}

m_suppress = fixef(lmer.4_S1j)[1]
diffMeans = fixef(lmer.4_S1j)[2]
m_feedback  = fixef(lmer.4_S1j)[1] + diffMeans
```

Feedback: burst ratio of `r format(m_feedback, digits=2, nsmall=2)` \newline
Suppression: burst ratio of `r format(m_suppress, digits=2, nsmall=2)` \newline
n = `r nrow(tb %>% drop_na(blankburstratios) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(blankburstratios) %>% count(mid))` mice 

\newpage

```{r read_data_4_S1hkil, include=FALSE}

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

```{r tidy_for_$_S1hkil, include = FALSE}

# Turn stimulus type into binary predictor, get relevant columns
tb <- tib %>% mutate(feedback = ifelse(opto == "TRUE", 0, 1)) %>% select(mid, sid, eid, uid, mseu, feedback, everything(), -opto)
```

\newpage
# Figure 4-Supplement 1h
## Feedback effect on firing rate during blank periods preceding gratings

```{r fit_model_4_S1h}

# Random intercept for neurons,
# random intercept for experiments, nested in series, nested in mice
lmer.4_S1h = lmer(blankrates ~ feedback + (1 + feedback | uid) + (1 | mid/sid/eid), 
                  data = tb %>% drop_na(blankrates))

display(lmer.4_S1h)
anova(lmer.4_S1h)
```

```{r predicted_average_effect_4_S1h, include=F}

m_suppress = fixef(lmer.4_S1h)[1]
diffMeans = fixef(lmer.4_S1h)[2]
m_feedback  = fixef(lmer.4_S1h)[1] + diffMeans
```

Feedback: `r format(m_feedback, digits=2, nsmall=2)` spikes/s \newline
Suppression: `r format(m_suppress, digits=2, nsmall=2)` spikes/s \newline
n = `r nrow(tb %>% drop_na(blankrates) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(blankrates) %>% count(mid))` mice 

\newpage
# Figure 4-Supplement 1k
## Feedback effect on burst ratio during blank periods preceding gratings

```{r fit_model_4_S1k}

# Random intercept, random slope for neurons,
# random intercept for series, nested in mice
lmer.4_S1k = lmer(blankburstratios ~ feedback + (1 + feedback | uid) + (1 | mid/sid), 
                  data = tb %>% drop_na(blankburstratios))

display(lmer.4_S1k)
anova(lmer.4_S1k)
```

```{r predicted_average_effect_4_S1k, include=F}

m_suppress = fixef(lmer.4_S1k)[1]
diffMeans = fixef(lmer.4_S1k)[2]
m_feedback  = fixef(lmer.4_S1k)[1] + diffMeans
```

Feedback: burst ratio of `r format(m_feedback, digits=2, nsmall=2)` \newline
Suppression: burst ratio of `r format(m_suppress, digits=2, nsmall=2)` \newline
n = `r nrow(tb %>% drop_na(blankburstratios) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(blankburstratios) %>% count(mid))` mice 

\newpage
# Figure 4-Supplement 1i
## Feedback effect on firing rate during zero-contrast gratings


```{r fit_model_4_S1i}

# Random intercept, random slope for neurons,
# random intercept for experiments, nested in mice
lmer.4_S1i = lmer(blankcondrates ~ feedback + (1 + feedback | uid) + (1 | sid/eid), 
                  data = tb %>% drop_na(blankcondrates))

display(lmer.4_S1i)
anova(lmer.4_S1i)
```

```{r predicted_average_effect_4_S1i, include=F}

m_suppress = fixef(lmer.4_S1i)[1]
diffMeans = fixef(lmer.4_S1i)[2]
m_feedback  = fixef(lmer.4_S1i)[1] + diffMeans
```

Feedback: `r format(m_feedback, digits=2, nsmall=2)` spikes/s \newline
Suppression: `r format(m_suppress, digits=2, nsmall=2)` spikes/s \newline
n = `r nrow(tb %>% drop_na(blankcondrates) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(blankcondrates) %>% count(mid))` mice 

\newpage
# Figure 4-Supplement 1l
## Feedback effect on burst ratio during zero-contrast gratings

```{r fit_model_4_S1l}

# Random intercept for neurons,
# random intercept for series, nested in mice
lmer.4_S1l = lmer(blankcondburstratios ~ feedback + (1 | uid) + (1 | mid/sid), 
                  data = tb %>% drop_na(blankcondburstratios))

display(lmer.4_S1l)
anova(lmer.4_S1l)
```

```{r predicted_average_effect_4_S1l, include=F}

m_suppress = fixef(lmer.4_S1l)[1]
diffMeans = fixef(lmer.4_S1l)[2]
m_feedback  = fixef(lmer.4_S1l)[1] + diffMeans
```

Feedback: burst ratio of `r format(m_feedback, digits=2, nsmall=2)` \newline
Suppression: burst ratio of `r format(m_suppress, digits=2, nsmall=2)` \newline
n = `r nrow(tb %>% drop_na(blankcondburstratios) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(blankcondburstratios) %>% count(mid))` mice