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

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

```{r tidy, include = FALSE}

# Turn booleans for 'optogenetic manipulation' into a binary predictor
tb <- tib %>% mutate(feedback = ifelse(opto == TRUE, 0, 1))
```

# Figure 1f
## Feedback effects on firing rate

```{r, fit_model_1f}

# We fit a random-intercept, random-slope model with two random effects: 
# (1) Neurons (uid) can have different baseline firing rates, 
# and the effect of feedback can vary across neurons. 
# (2) Mean firing rates are allowed to differ across experiments (eid), 
# which are nested within recording sessions (sid), 
# which are nested within animals (mid).

lmer.1f = lmer(rates ~ feedback + (1 + feedback | uid) + (1 | mid/sid/eid), 
               data = tb %>% drop_na(rates))

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

```{r get_predicted_average_effect_1f, include=F}

mSuppr = fixef(lmer.1f)[1]
diffMeans = fixef(lmer.1f)[2]
mFeedbk  = fixef(lmer.1f)[1] + diffMeans
```

Feedback: mean firing rate of `r format(mFeedbk, digits=2, nsmall=1)` spikes/s \newline
Suppression: mean firing rate of `r format(mSuppr, digits=2, 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 1g
## Feedback effects on burst ratio

```{r fit_model_1g}

# Random-intercept, random-slope for single neurons, 
# random intercept for experiments.
# Variability across series and mice is close to zero, 
# including random intercepts for those gives singular fits.

lmer.1g = lmer(burstratios ~ feedback + (1 + feedback | uid) + (1 | eid), 
               data = tb %>% drop_na(burstratios))

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

```{r get_predicted_average_effect_1g, include=F}

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

Feedback: mean burst ratio of `r format(mFeedbk, digits=1, nsmall=1)` \newline
Suppression: mean burst ratio of `r format(mSuppr, digits=1, nsmall=1)` \newline
n = `r nrow(tb %>% drop_na(burstratios) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(burstratios) %>% count(mid))` mice

\newpage
# Figure 1h
## Feedback effects on sparseness

```{r tidy_for_1hi, include=FALSE}

# 'Sparseness', and 'reliability' 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

tbhi = tb %>% select(mid, sid, eid, uid, mseu, feedback, spars, rel) %>% distinct(mseu, feedback, .keep_all = TRUE)
```

```{r fit_model_1h}

# Random-intercept, random-slope for single neurons, 
# random intercept for experiments, nested within series
lmer.1h = lmer(spars ~ feedback + (1 + feedback | uid) + (1 | sid/eid), 
               data = tbhi %>% drop_na(spars))

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

```{r get_predicted_average_effect_1h, include=F}

mSuppr = fixef(lmer.1h)[1]
diffMeans = fixef(lmer.1h)[2]
mFeedbk  = fixef(lmer.1h)[1] + diffMeans
```

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

\newpage
# Figure 1i
## Feedback effects on reliability

```{r fit_model_1i}

# Random-intercept, random-slope for single neurons, 
# random intercept for experiments, nested within series
lmer.1i = lmer(rel ~ feedback + (1 + feedback | uid) + (1 | sid/eid), 
               data = tbhi %>% drop_na(rel))

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

```{r get_predicted_average_effect_1i, include=F}
mSuppr = fixef(lmer.1i)[1]
diffMeans = fixef(lmer.1i)[2]
mFeedbk  = fixef(lmer.1i)[1] + diffMeans

```

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