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

```{r tidy, include = FALSE}

# Transform opto's "false" and "true" into a binary, numeric predictor
tb <- tib %>% mutate(feedback = ifelse(opto == TRUE, 0, 1))
```

# Figure 3c
## Feedback effects on firing rate

```{r fit_model_3c}

# Random-intercept, random-slope for single neurons, 
# random intercept for experiments, nested in series
lmer.3c = lmer(rates ~ feedback + (1 + feedback | uid) + (1 | sid/eid), 
               data = tb %>% drop_na(rates))

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

```{r predicted_average_effect_3c, include=F}

mSuppr = fixef(lmer.3c)[1]
diffMeans = fixef(lmer.3c)[2]
mActive  = fixef(lmer.3c)[1] + diffMeans
```

Feedback: `r format(mActive, digits=1, nsmall=1)` spikes/s \newline
Suppression: `r format(mSuppr, digits=1, nsmall=1)` spikes/s \newline
n = `r nrow(tb %>% drop_na(rates) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(rates) %>% count(mid))` mice

\newpage
# Figure 3d
## Feedback effects on burst ratio

```{r fit_model_3d}

# Random-intercept, random-slope for single neurons, 
# random intercept for experiments, nested in series
lmer.3d = lmer(burstratios ~ feedback + (1 + feedback | uid) + (1 | sid/eid), 
               data = tb %>% drop_na(burstratios))

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

```{r predicted_average_effect_3d, include=F}

mSuppr = fixef(lmer.3d)[1]
diffMeans = fixef(lmer.3d)[2]
mActive  = fixef(lmer.3d)[1] + diffMeans
```

Feedback: `r format(mActive, digits=2, nsmall=2)` \newline
Suppression: `r format(mSuppr, digits=2, nsmall=2)` \newline
n = `r nrow(tb %>% drop_na(burstratios) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(burstratios) %>% count(mid))` mice

\newpage
# Figure 3e
## Feedback effects on orientation selectivity

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

# 'Orientation selectivity', and 'F1/F0 ratio' are not computed on a trial-by-trial basis. In the csv-file, 
# these two measures are therefore identical across trials, so we simply pull out the first trial of each neuron

tbeg = tb %>% select(mid, sid, eid, uid, mseu, feedback, osi, f1f0) %>% distinct(mseu, feedback, .keep_all = TRUE)
```

```{r fit_model_3e}

# Random-intercept for single neurons, 
# random intercept for experiments, nested in series
lmer.3e = lmer(osi ~ feedback + (1 | uid) + (1 | sid/eid), 
               data = tbeg %>% drop_na(osi))

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

```{r predicted_average_effect_3e, include=F}

mSuppr = fixef(lmer.3e)[1]
diffMeans = fixef(lmer.3e)[2]
mActive  = fixef(lmer.3e)[1] + diffMeans
```

Feedback: OSI = `r format(mActive, digits=2, nsmall=2)` \newline
Suppression: OSI = `r format(mSuppr, digits=2, nsmall=2)` \newline
n = `r nrow(tbeg %>% drop_na(osi) %>% count(uid))` neurons from `r nrow(tbeg %>% drop_na(osi) %>% count(mid))` mice

\newpage
# Figure 3g
## Feedback effects on F1/F0-ratio

```{r fit_model_3g}

# Random intercept, random slope for single neurons,
# random intercept for series
lmer.3g = lmer(f1f0 ~ feedback + (1 + feedback | uid) + (1 | sid), 
               data = tbeg %>% drop_na(f1f0))

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

```{r predicted_average_effect_3g, include=F}

mSuppr = fixef(lmer.3g)[1]
diffMeans = fixef(lmer.3g)[2]
mActive  = fixef(lmer.3g)[1] + diffMeans
```

Feedback: $F_1/F_0$-ratio = `r format(mActive, digits=2, nsmall=2)` \newline
Suppression: $F_1/F_0$-ratio = `r format(mSuppr, digits=2, nsmall=2)` \newline
n = `r nrow(tbeg %>% drop_na(f1f0) %>% count(uid))` neurons from `r nrow(tbeg %>% drop_na(f1f0) %>% count(mid))` mice


\newpage
# Figure 3i
## Feedback effects on distribution of cycle average phase differences

```{r tidy_3i, include=F}

# Remove redundant rows, select relevant variables
tb3i <- tib %>% filter(opto == "TRUE") %>% select(mid, sid, eid, uid, mseu, dbphi, dnbphi) %>% distinct(mseu, .keep_all = TRUE)

# Grab bursting phis
tb3i_db = tb3i %>% filter(dbphi != "NA")

# Grab non-bursting phis
tb3i_dnb = tb3i %>% filter(dnbphi != "NA")
```

```{r test_3i_bursting}

# Of all bursting phis, how many are phase-advanced?
b_adv = sum(tb3i_db$dbphi > 0)
resb = binom.test(b_adv, length(tb3i_db$dbphi), 0.5)
print(resb)
```

\vspace{0.5cm}

```{r test_3i_non_bursting}

# Of all non-bursting phis, how many are phase-advanced?
nb_adv = sum(tb3i_dnb$dnbphi > 0)
resnb = binom.test(nb_adv, length(tb3i_dnb$dnbphi), 0.5)
print(resnb)
```