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

# Figure 4a
## Modulation of firing rate by feedback during movie or grating presentations

```{r tidy_4a, include=F}

# Grab conditions, in which a stimulus is shown
t4a = tib %>% filter(blank == "False" & meanrate != "NA")

# Turn stimtype into factor for dummy coding
t4a$stimtype = as.factor(t4a$stimtype)

```

```{r fit_model_4a}

# Random intercept for neurons
lmer.4a = lmer(meanrate ~ stimtype + (1 | uid), 
               data = t4a %>% drop_na(meanrate))

display(lmer.4a)
anova(lmer.4a)
```

```{r get_predicted_average_effect_4a, include=F}

m_grt = fixef(lmer.4a)[1]
m_mvi = fixef(lmer.4a)[1] + fixef(lmer.4a)[2]

# Store coefficients in file
pred_means_df = data.frame("grt" = m_grt, "mvi" = m_mvi, row.names = "")
write_csv(pred_means_df, "_stats/figure_4a_pred_means.csv")
```

Grating: FMI = `r format(m_grt, digits=2, nsmall=2)` \newline
Movie: FMI = `r format(m_mvi, digits=2, nsmall=2)` \newline
n = `r nrow(t4a %>% drop_na(meanrate) %>% count(uid))` neurons from `r nrow(t4a %>% drop_na(meanrate) %>% count(mid))` mice

\newpage
# Figure 4b
## Modulation of firing rate by feedback during gray screen conditions

```{r tidy_4b, include=FALSE}

# We code the 'blank condition' as integers 1, 2, 3
t4b = tib %>% mutate(blank_condition = case_when(
  stimtype == 'mvi' & blank == 'prestim' ~ 'mvi',
  stimtype == 'grt' & blank == 'prestim' ~ 'grt',
  stimtype == 'grt' & blank == 'cond' ~ 'grt0c'
))

# Remove NAs
t4b = t4b %>% filter(meanrate != "Na" & blank_condition != "Na")

# And turn blank_condition into a factor for dummy coding
t4b$blank_condition = as.factor(t4b$blank_condition)

```

```{r fit_model_4b}

# Random intercept for neurons
lmer.4b = lmer(meanrate ~ blank_condition + (1 | uid), 
               data = t4b %>% drop_na(meanrate))

display(lmer.4b)
anova(lmer.4b)
```

```{r get_predicted_average_effect_4b, include=F}

# Get predicted average burst ratios
m_grt = fixef(lmer.4b)[1]
m_grt0c = fixef(lmer.4b)[1] + fixef(lmer.4b)[2]
m_mvi = fixef(lmer.4b)[1]+ fixef(lmer.4b)[3]

# store results in file
pred_means_df = data.frame("grt" = m_grt, "grt0c" = m_grt0c, "mvi" = m_mvi, row.names = "")
write_csv(pred_means_df, "_stats/figure_4b_pred_means.csv")
```

Blank before movie: FMI = `r format(m_mvi, digits=2, nsmall=2)` \newline
Blank before grating: FMI = `r format(m_grt, digits=2, nsmall=2)` \newline
Zero contrast grating: FMI = `r format(m_grt0c, digits=2, nsmall=2)` \newline
n = `r nrow(t4b %>% drop_na(meanrate) %>% count(uid))` neurons from `r nrow(t4b %>% drop_na(meanrate) %>% count(mid))` mice


```{r tidy_for_I, include = F}

t4_p = tib %>% filter(meanrate != "Na")

t4_p = t4_p %>% mutate(stim_condition = case_when(
  stimtype == 'mvi' & blank == 'False' ~ 'mvi',
  stimtype == 'grt' & blank == 'False' ~ 'grt',
  (stimtype == 'mvi' | stimtype == 'grt') & blank != 'False' ~ 'blank'
))

t4_p$stim_condition = as.factor(t4_p$stim_condition)

# Is mean rate FMI for blanks different from the one for movies?
t4_p1 = t4_p %>% filter(stim_condition != 'grt')

```

\newpage
## (I) Pooling across blank conditions, and comparing against the FMI for movies


```{r fit_model_I}

# Random intercept for neurons
lmer.I = lmer(meanrate ~ stim_condition  + (1 | uid), 
               data = t4_p1 %>% drop_na(meanrate))

display(lmer.I)
anova(lmer.I)

```

n = `r nrow(t4_p1 %>% drop_na(meanrate) %>% count(uid))` neurons from `r nrow(t4_p1 %>% drop_na(meanrate) %>% count(mid))` mice

\newpage
## (II) Pooling across blank conditions, and comparing against the FMI for gratings

```{r tidy_for_II, include=F}

# Is mean rate FMI for blanks different from the one for gratings?
t4_p2 = t4_p %>% filter(stim_condition != 'mvi')
```

```{r fit_model_II}

# Random intercept for neurons
lmer.II = lmer(meanrate ~ stim_condition  + (1 | uid), t4_p2)
display(lmer.II)
anova(lmer.II)
```

n = `r nrow(t4_p2 %>% drop_na(meanrate) %>% count(uid))` neurons from `r nrow(t4_p2 %>% drop_na(meanrate) %>% count(mid))` mice