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

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

```{r tidy_for_1_S3, include = FALSE}

# Remove observations with NAs
tb <- tib %>% filter(sbc != "Na")

# Recode suppressed-by-contrast boolean into numeric, binary predictor
tb$sbc = ifelse(tb$sbc == TRUE, 1, 0)
```

```{r fit_model_1_S3a}

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

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

```{r get_predicted_average_change_1_S3a, include=F}

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

```{r store_coefficients_1_S3a, include=FALSE}

pred_means_df = data.frame("non_sbc" = m_sbc0, "sbc" = m_sbc1, row.names = "")
write_csv(pred_means_df, "_stats/figure_1_S3a_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(mvi_meanrate) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(mvi_meanrate) %>% count(mid))` mice

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

```{r fit_model_1_S3b}

# Random intercept for series
lmer.1_S3b = lmer(mvi_meanrate ~ depth  + (1 | sid), 
                  data = tib %>% drop_na(mvi_meanrate, depth))

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

```{r store_coefficients, include=FALSE}

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

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

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

``` {r, fit_model_1_S3c}

# Random intercept for series
lmer.1_S3c = lmer(mvi_meanrate ~ dsi + (1 | sid), 
                  data = tib %>% drop_na(mvi_meanrate, dsi))

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

```{r store_coefficients_1_S3c, include=FALSE}

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

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

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

``` {r, fit_model_1_S3d}

# Random intercept for series nested within mice
lmer.1_S3d = lmer(mvi_meanrate ~ rfdist + (1 | mid/sid), 
                  data = tib %>% drop_na(mvi_meanrate, rfdist))
display(lmer.1_S3d)
anova(lmer.1_S3d)
```

```{r store_coefficients_1_S3d, include=FALSE}

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

```

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

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

``` {r, fit_model_1_S3e}

# Random intercept for series nested within mice
lmer.1_S3e = lmer(mvi_meanrate ~ mvi_meanrate_raw + (1 | mid/sid), 
                  data = tib %>% drop_na(mvi_meanrate, mvi_meanrate_raw))

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

```{r store_coefficients_1_S3e, include=FALSE}

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

```

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

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

```{r fit_model_1_S3f}

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

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

```{r get_predicted_average_effect_1_S3f, include=F}

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

```{r store_coefficients_1_S3f, include=FALSE}

pred_means_df = data.frame("non_sbc" = m_sbc0, "sbc" = m_sbc1, row.names = "")
write_csv(pred_means_df, "_stats/figure_1_S3f_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(mvi_meanburstratio) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(mvi_meanburstratio) %>% count(mid))` mice

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

```{r fit_model_1_S3g}

# Random intercept for series nested within mice
lmer.1_S3g = lmer(mvi_meanburstratio ~ depth  + (1 | mid/sid), 
                  data = tib %>% drop_na(mvi_meanburstratio, depth))

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

```{r store_coefficients_1_S3g, include=FALSE}

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

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

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

``` {r, fit_model_1_S3h}

# Random intercept for series nested within mice
lmer.1_S3h = lmer(mvi_meanburstratio ~ dsi + (1 | mid/sid), 
                  data = tib %>% drop_na(mvi_meanburstratio, dsi))

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

```{r store_coefficients_1_S3h, include=FALSE}

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

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

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

``` {r, fit_model_1_S3i}

# Random intercept for series nested within mice
lmer.1_S3i = lmer(mvi_meanburstratio ~ rfdist + (1 | mid/sid), 
                  data = tib %>% drop_na(mvi_meanburstratio, rfdist))

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

```{r store_coefficients_1_S3i, include=FALSE}

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

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

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

``` {r, fit_model_1_S3j}

# Random intercept for series
lmer.1_S3j = lmer(mvi_meanburstratio ~ mvi_meanburstratio_raw + (1 | sid), 
                  data = tib %>% drop_na(mvi_meanburstratio, mvi_meanburstratio_raw))

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

```{r store_coefficients_1_S3j, include=FALSE}

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

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