natfeedback_code_eLife/R/figure_3_S1.Rmd

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

# This is the same file we used for Figure 1-Supplement 3. It contains the dependent variables for both, movies (Fig 1-Suppl 3) and gratings (Fig 3-Suppl 1)
tib = get_data("../csv/fig1S33S1.csv")
```

```{r tidy_for_3_S1, include = FALSE}

# Remove observations with NAs
tb <- tib %>% drop_na(sbc)

# Turn cell type (suppressed-by-contrast) into binary variable
tb$sbc = ifelse(tb$sbc == TRUE, 1, 0)
```

# Figure 3-Supplement 1a
## Firing rate FMIs separated by cell types

```{r fit_model_3_S1a}

# A mixed-effects model with any random intercept gives singular fits,
# we therefore revert to fitting an ordinary linear model

lmer.3_S1a = lm(grt_meanrate ~ sbc,
                  data = tb %>% drop_na(grt_meanrate))

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

```{r get_predicted_average_effect_3_S1a, include=F}

m_sbc0 = coef(lmer.3_S1a)[1]
diffMeans = coef(lmer.3_S1a)[2]
m_sbc1  = coef(lmer.3_S1a)[1] + diffMeans
```

```{r save_coefficients_3_S1a, include=FALSE}

pred_means_df = data.frame("non_sbc" = m_sbc0, "sbc" = m_sbc1, row.names = "")
write_csv(pred_means_df, "_stats/figure_3_S1a_pred_means.csv")
```

FMI sbc: `r format(m_sbc1, digits=3, nsmall=3)` \newline
FMI non-sbc: `r format(m_sbc0, digits=3, nsmall=3)` \newline
n = `r nrow(tb %>% drop_na(grt_meanrate) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(grt_meanrate) %>% count(mid))` mice

\newpage
# Figure 3-Supplement 1b
## Relation between firing rate FMI and recording depth

```{r fit_model_3S1b}

# Random intercept for mice
lmer.3_S1b = lmer(grt_meanrate ~ depth  + (1 | mid), 
                  data = tib %>% drop_na(grt_meanrate, depth))

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

```{r save_coefficients_3_S1b, include=FALSE}

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

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

\newpage
# Figure 3-Supplement 1c
## Relation between firing rate FMI and direction selectivity

``` {r fit_model_3_S1c}

# Random intercept for series, nested in mice
lmer.3_S1c = lmer(grt_meanrate ~ dsi + (1 | sid/mid), 
                  data = tib %>% drop_na(grt_meanrate, dsi))

display(lmer.3_S1c)
anova(lmer.3_S1c)
```

```{r save_coefficients_3_S1c, include=FALSE}

coef_df = data.frame("intercept" = fixef(lmer.3_S1c)[1], "slope" = fixef(lmer.3_S1c)[2], row.names = "")
write_csv(coef_df, "_stats/figure_3_S1c_coefs.csv")

```

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

\newpage
# Figure 3-Supplement 1d
## Relation between firing rate FMI and receptive field location

``` {r fit_model_3_S1d}

# Random intercept for series, nested in mice
lmer.3_S1d = lmer(grt_meanrate ~ rfdist + (1 | mid/sid), 
                  data = tib %>% drop_na(grt_meanrate, rfdist))

display(lmer.3_S1d)
anova(lmer.3_S1d)
```

```{r save_coefficients_3_S1d,include=FALSE}

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

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

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

``` {r, fit_model_3_S1e}

# Random intercept for series, nested in mice
lmer.3_S1e = lmer(grt_meanrate ~ grt_meanrate_raw + (1 | mid/sid), 
                  data = tib %>% drop_na(grt_meanrate, grt_meanrate_raw))

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

```{r save_coefficients_3_S1e, include=F}

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

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

\newpage
# Figure 3-Supplement 1f
## Burst ratio FMIs separated by cell types

```{r fit_model_3_S1f}

# Random intercept for series
lmer.3_S1f = lmer(grt_meanburstratio ~ sbc  + (1 | sid), 
                  data = tb %>% drop_na(grt_meanburstratio))

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

```{r get_predicted_average_effect_3_S1f, include=F}

m_sbc0 = fixef(lmer.3_S1f)[1]
diffMeans = fixef(lmer.3_S1f)[2]
m_sbc1  = fixef(lmer.3_S1f)[1] + diffMeans
```

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

pred_means_df = data.frame("non_sbc" = m_sbc0, "sbc" = m_sbc1, row.names = "")
write_csv(pred_means_df, "_stats/figure_3_S1f_pred_means.csv")
```

FMI sbc: `r format(m_sbc1, digits=3, nsmall=3)` \newline
FMI non-sbc: `r format(m_sbc0, digits=3, nsmall=3)` \newline
n = `r nrow(tb %>% drop_na(grt_meanburstratio) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(grt_meanburstratio) %>% count(mid))` mice

\newpage
# Figure 3-Supplement 1g
## Relation between burst ratio FMI and recording depth


```{r fit_model_3_S1g}

# Random intercept for series, nested in mice
lmer.3_S1g = lmer(grt_meanburstratio ~ depth  + (1 | mid/sid), 
                  data = tib %>% drop_na(grt_meanburstratio, depth))

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

```{r save_coefficients_3S1g, include=FALSE}

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

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

\newpage
# Figure 3-Supplement 1h
## Relation between burst ratio FMI and direction selectivity

``` {r, fit_model_3_S1h}

# Random intercept for series, nested in mice
lmer.3_S1h = lmer(grt_meanburstratio ~ dsi + (1 | mid/sid), 
                  data = tib %>% drop_na(grt_meanburstratio, dsi))

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

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

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

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

\newpage
# Figure 3-Supplement 1i
## Relation between burst ratio FMI and receptive field location

``` {r, fit_model_3_S1i}

# Random intercept for series, nested in mice
lmer.3_S1i = lmer(grt_meanburstratio ~ rfdist + (1 | mid/sid), 
                  data = tib %>% drop_na(grt_meanburstratio, rfdist))

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

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

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

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


\newpage
# Figure 3-Supplement 1j
## Relation between burst ratio FMI and burst ratio

``` {r, fit_model_3_S1j}

# Random intercept for series, nested in mice
lmer.3_S1j = lmer(grt_meanburstratio ~ grt_meanburstratio_raw + (1 | mid/sid), 
                  data = tib %>% drop_na(grt_meanburstratio, grt_meanburstratio_raw))

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

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

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

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