---
title: "Spacek et al., 2021, Figure 5-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}

# Read data
tib = get_data("../csv/fig5S2.csv")
```

```{r tidy, include = FALSE}

# Extract data for control condition
tb_control <- tib %>% filter(opto == FALSE)

# Extract data with suppression
tb_supp <- tib %>% filter(opto == TRUE)

```

# Figure 5-Supplement 2d
## Reliability - V1 control

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

tb2d = tb_control %>% select(mid, sid, eid, uid, mseu, tqrel, bqrel)

# Turn into long format, and create binary predictor for top vs bottom quartile
tb2d = tb2d %>% pivot_longer(cols = c(tqrel, bqrel), names_to = "quartiles", values_to = "rel")
tb2d = tb2d %>% mutate(bottom_quartile = ifelse(quartiles == "bqrel", 1, 0))
```

```{r fit_model_5_S2d}

# Random intercept for single neurons, 
# random intercept for experiments, nested within series
lmer.5_S2d = lmer(rel ~ bottom_quartile + (1 | uid) + (1 | sid/eid),
               data = tb2d %>% drop_na(rel))

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

```{r get_predicted_average_effect_5_S2d, include=F}

mTop = fixef(lmer.5_S2d)[1]
diffMeans = fixef(lmer.5_S2d)[2]
mBottom  = fixef(lmer.5_S2d)[1] + diffMeans
```

Top quartile: reliability of `r format(mTop, digits=2, nsmall=2)` \newline
Bottom quartile: reliability of `r format(mBottom, digits=2, nsmall=2)` \newline
n = `r nrow(tb2d %>% drop_na(rel) %>% count(uid))` neurons from `r nrow(tb2d %>% drop_na(rel) %>% count(mid))` mice

\newpage
# Figure 5-Supplement 2e
## Signal-to-noise ratio - V1 control

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

tb2e = tb_control %>% select(mid, sid, eid, uid, mseu, tqsnr, bqsnr)

# Turn into long format, and create binary predictor for top vs bottom quartile
tb2e = tb2e %>% pivot_longer(cols = c(tqsnr, bqsnr), names_to = "quartiles", values_to = "snr")
tb2e = tb2e %>% mutate(bottom_quartile = ifelse(quartiles == "bqsnr", 1, 0))
```

```{r fit_model_5_S2e}

# Random intercept for single neurons, 
# random intercept for experiments
lmer.5_S2e = lmer(snr ~ bottom_quartile + (1 | uid) + (1 | eid),
               data = tb2e %>% drop_na(snr))

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

```{r get_predicted_average_effect_5_S2e, include=F}

mTop = fixef(lmer.5_S2e)[1]
diffMeans = fixef(lmer.5_S2e)[2]
mBottom  = fixef(lmer.5_S2e)[1] + diffMeans
```

Top quartile: SNR of `r format(mTop, digits=2, nsmall=2)` \newline
Bottom quartile: SNR of `r format(mBottom, digits=2, nsmall=2)` \newline
n = `r nrow(tb2e %>% drop_na(snr) %>% count(uid))` neurons from `r nrow(tb2e %>% drop_na(snr) %>% count(mid))` mice

\newpage
# Figure 5-Supplement 2f
## Reliability - V1 suppressed

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

tb2f = tb_supp %>% select(mid, sid, eid, uid, mseu, tqrel, bqrel)

# Turn into long format, and create binary predictor for top vs bottom quartile
tb2f = tb2f %>% pivot_longer(cols = c(tqrel, bqrel), names_to = "quartiles", values_to = "rel")
tb2f = tb2f %>% mutate(bottom_quartile = ifelse(quartiles == "bqrel", 1, 0))
```

```{r fit_model_5_S2f}

# Random intercept for single neurons, 
# random intercept for experiments, nested within series
lmer.5_S2f = lmer(rel ~ bottom_quartile + (1 | uid) + (1 | sid/eid),
               data = tb2f %>% drop_na(rel))

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

```{r get_predicted_average_effect_5_S2f, include=F}

mTop = fixef(lmer.5_S2f)[1]
diffMeans = fixef(lmer.5_S2f)[2]
mBottom  = fixef(lmer.5_S2f)[1] + diffMeans
```

Top quartile: reliability of `r format(mTop, digits=2, nsmall=2)` \newline
Bottom quartile: reliability of `r format(mBottom, digits=2, nsmall=2)` \newline
n = `r nrow(tb2f %>% drop_na(rel) %>% count(uid))` neurons from `r nrow(tb2f %>% drop_na(rel) %>% count(mid))` mice

\newpage
# Figure 5-Supplement 2g
## Signal-to-noise ratio - V1 suppressed

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

tb2g = tb_supp %>% select(mid, sid, eid, uid, mseu, tqsnr, bqsnr)

# Turn into long format, and create binary predictor for top vs bottom quartile
tb2g = tb2g %>% pivot_longer(cols = c(tqsnr, bqsnr), names_to = "quartiles", values_to = "snr")
tb2g = tb2g %>% mutate(bottom_quartile = ifelse(quartiles == "bqsnr", 1, 0))
```

```{r fit_model_5_S2g}

# Random intercept for single neurons, 
# random intercept for experiments, nested within series
lmer.5_S2g = lmer(snr ~ bottom_quartile + (1 | uid) + (1 | sid/eid),
               data = tb2g %>% drop_na(snr))

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

```{r get_predicted_average_effect_5_S2g, include=F}

mTop = fixef(lmer.5_S2g)[1]
diffMeans = fixef(lmer.5_S2g)[2]
mBottom  = fixef(lmer.5_S2g)[1] + diffMeans
```

Top quartile: SNR of `r format(mTop, digits=2, nsmall=2)` \newline
Bottom quartile: SNR of `r format(mBottom, digits=2, nsmall=2)` \newline
n = `r nrow(tb2g %>% drop_na(snr) %>% count(uid))` neurons from `r nrow(tb2g %>% drop_na(snr) %>% count(mid))` mice