natfeedback_code_eLife/R/figure_6_S1.Rmd

---
title: "Spacek et al., 2021, Figure 6-Supplement 1"
output: pdf_document
---

```{r setup, include=FALSE}

knitr::opts_chunk$set(echo = TRUE)

library(arm)
library(lmerTest)
library(tidyverse)
library(data.table)

source('get_data.R')
```

```{r read_data, include=FALSE}

tib = get_data("../csv/fig6.csv")
tib <- tib %>% filter(stimtype == "grt") %>% select(mid, sid, eid, uid, mseu, mi, meanrate, meanburstratio)
```

```{r tidy_for_a1, include = FALSE}

# Select relevant columns
tb <- tib %>% filter(mi == "suppressionrmi" | mi == "feedbackrmi")

# Turn long format into wide format, using data.tables and dcast:
foo = dcast(setDT(tb), mid + sid + eid + uid + mseu ~ mi, value.var = "meanrate")
tbw = as_tibble(foo)
```

# Figure 6-S1a$_1$
## Comparing RMI during suppression against 0

```{r, fit_model_a1_1}

# Fixed effect intercept only, 
# random intercept for neurons,
# random intercept for experiments, nested within series
lmer.a1.1 = lmer(suppressionrmi ~ 1 + (1 | uid) + (1 | sid/eid), 
                  data = tbw %>% drop_na(suppressionrmi))

display(lmer.a1.1)
```

## Mean firing rate RMI
Suppression: RMI = `r format(fixef(lmer.a1.1)[1], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.a1.1)[1], digits=2, nsmall=2)` \newline
n = `r nrow(tbw %>% drop_na(suppressionrmi) %>% count(uid))` neurons from `r nrow(tbw %>% drop_na(suppressionrmi) %>% count(mid))` mice

\newpage
# Figure 6-S1a$_2$
## Comparing RMI during suppression against 0

```{r tidy_for_a2, include = FALSE}

# Turn long format into wide format, using data.tables and dcast:
foo = dcast(setDT(tb), mid + sid + eid + uid + mseu ~ mi, value.var = "meanburstratio")
tbw = as_tibble(foo)
```

```{r fit_model_a2_1}

# LMM gives singular fits, even with a single random intercept for neurons. 
# We do ordinary regression instead:
lm.a2.1 = lm(suppressionrmi ~ 1, data = tbw %>% drop_na(suppressionrmi))

display(lm.a2.1)
```

## Mean burst ratio RMI
Suppression: RMI = `r format(summary(lm.a2.1)$coefficients[1], digits=2, nsmall=2)` $\pm$ `r format(2 * summary(lm.a2.1)$coefficients[2], digits=1, nsmall=1)` \newline
n = `r nrow(tbw %>% drop_na(suppressionrmi) %>% count(uid))` neurons from `r nrow(tbw %>% drop_na(suppressionrmi) %>% count(mid))` mice

\newpage
# Figure 6-S1b$_1$
## Slope of regression line

```{r tidy_for_b1_1, include = FALSE}

tb <- tib %>% filter(mi == "sitfmi" | mi == "runfmi")

# Turn long format into wide format, using data.tables and dcast:
foo = dcast(setDT(tb), mid + sid + eid + uid + mseu ~ mi, value.var = "meanrate")
tbw = as_tibble(foo)
```

```{r fit_model_b1_1}

# Random intercept for neurons,
# random intercept for experiments
lmer.b1.1 = lmer(runfmi ~ sitfmi + (1 | uid) + (1 | eid), 
                  data = tbw %>% drop_na(runfmi, sitfmi))

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

``` {r store_coefficients_b1_1, include=FALSE}

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

Slope of `r format(fixef(lmer.b1.1)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.b1.1)[2], digits=2, nsmall=2)` \newline
n = `r nrow(tbw %>% drop_na(runfmi, sitfmi) %>% count(uid))` neurons from `r nrow(tbw %>% drop_na(runfmi, sitfmi) %>% count(mid))` mice

\newpage
# Figure 6-S1b$_2$
## Slope of regression line

```{r tidy_for_b2_1, include = FALSE}

tb <- tib %>% filter(mi == "sitfmi" | mi == "runfmi")

# Turn long format into wide format, using data.tables and dcast:
foo = dcast(setDT(tb), mid + sid + eid + uid + mseu ~ mi, value.var = "meanburstratio")
tbw = as_tibble(foo)

# Remove NaNs
tbw <- tbw %>% filter(sitfmi != "Na" & runfmi != "Na")
```

```{r fit_model_b2_1}

# Random intercept for neurons,
# random intercept for experiments, nested in series
lmer.b2.1 = lmer(runfmi ~ sitfmi + (1 | uid) + (1 | sid/eid), 
                  data = tbw %>% drop_na(runfmi, sitfmi))

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

```{r store_coefficients_b2_1, include=FALSE}

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

Slope of `r format(fixef(lmer.b2.1)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.b2.1)[2], digits=2, nsmall=2)` \newline
n = `r nrow(tbw %>% drop_na(runfmi, sitfmi) %>% count(uid))` neurons from `r nrow(tbw %>% drop_na(runfmi, sitfmi) %>% count(mid))` mice


\newpage
# Figure 6-S1c$_1$
## Slope of regression line

```{r tidy_for_c1_1, include = FALSE}

tb <- tib %>% filter(mi == "rmi" | mi == "fmi")

# Turn long format into wide format, using data.tables and dcast:
foo = dcast(setDT(tb), mid + sid + eid + uid + mseu ~ mi, value.var = "meanrate")
tbw = as_tibble(foo)

```

```{r fit_model_c1_1}

# Random intercept for neurons,
# random intercept for experiments
lmer.c1.1 = lmer(fmi ~ rmi + (1 | uid) + (1 | eid), 
                data = tbw %>% drop_na(fmi, rmi))

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

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

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

Slope of `r format(fixef(lmer.c1.1)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.c1.1)[2], digits=2, nsmall=2)` \newline
n = `r nrow(tbw %>% drop_na(fmi, rmi) %>% count(uid))` neurons from `r nrow(tbw %>% drop_na(fmi, rmi) %>% count(mid))` mice

\newpage
# Figure 6-S1c$_2$
## Slope of regression line

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

tb <- tib %>% filter(mi == "rmi" | mi == "fmi")

# Turn long format into wide format, using data.tables and dcast:
foo = dcast(setDT(tb), mid + sid + eid + uid + mseu ~ mi, value.var = "meanburstratio")
tbw = as_tibble(foo)
```

```{r fit_model_c2_1}

# Random intercept for neurons, 
# random intercept for experiments, nested in series
lmer.c2.1 = lmer(fmi ~ rmi + (1 | uid) + (1 | sid/eid), 
                data = tbw %>% drop_na(fmi, rmi))

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

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

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

Slope of `r format(fixef(lmer.c2.1)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.c2.1)[2], digits=2, nsmall=2)` \newline
n = `r nrow(tbw %>% drop_na(fmi, rmi) %>% count(uid))` neurons from `r nrow(tbw %>% drop_na(fmi, rmi) %>% count(mid))` mice