natfeedback_code_eLife/R/figure_6.Rmd

---
title: "Spacek et al., 2021, Figure 6"
output: pdf_document
---

```{r setup, include=FALSE}

knitr::opts_chunk$set(echo = TRUE)

library(arm)
library(lmerTest)
library(tidyverse)
library(data.table)

source('get_data.R')
```

```{r read_data, include=FALSE}

tib = get_data("../csv/fig6.csv")
tib <- tib %>% filter(stimtype == "mvi") %>% select(mid, sid, eid, uid, mseu, mi, meanrate, meanburstratio, spars, rel)
```

```{r tidy_for_6a_1, include = FALSE}

# Select relevant columns
tb <- tib %>% filter(mi == "suppressionrmi" | mi == "feedbackrmi")

# Turn long format into wide format, using data.tables and dcast:
foo = dcast(setDT(tb), mid + sid + eid + uid + mseu ~ mi, value.var = "meanrate")
tbw = as_tibble(foo)

```

# Figure 6a$_1$
## (1) Comparing RMI during suppression against 0

```{r, fit_model_6a1_1}

# Fixed effect intercept only, 
# random intercept for neurons,
# random intercept for experiments, nested within series
lmer.6a1.1 = lmer(suppressionrmi ~ 1 + (1 | uid) + (1 | sid/eid), 
                  data = tbw %>% drop_na(suppressionrmi))

display(lmer.6a1.1)

```

## Mean firing rate RMI
Suppression: RMI = `r format(fixef(lmer.6a1.1)[1], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.6a1.1)[1], digits=1, nsmall=1)` \newline
n = `r nrow(tbw %>% drop_na(suppressionrmi) %>% count(uid))` neurons from `r nrow(tbw %>% drop_na(suppressionrmi) %>% count(mid))` mice


\newpage
# Figure 6a$_1$
## (2) Slope of regression line

```{r fit_model_6a1_2}

# Random intercept for neurons,
# random intercept for experiments, nested in series
lmer.6a1.2 = lmer(feedbackrmi ~ suppressionrmi + (1 | uid) + (1 | sid/eid), 
                  data = tbw %>% drop_na(feedbackrmi, suppressionrmi))

display(lmer.6a1.2)
anova(lmer.6a1.2)
```

``` {r save_coefficients_6a1_2, include=FALSE}

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

Slope = `r format(fixef(lmer.6a1.2)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.6a1.2)[2], digits=2, nsmall=2)` \newline
n = `r nrow(tbw %>% drop_na(feedbackrmi, suppressionrmi) %>% count(uid))` neurons from `r nrow(tbw %>% drop_na(feedbackrmi, suppressionrmi) %>% count(mid))` mice

\newpage
# Figure 6a$_1$
## (3) Average effect of feedback on firing rate RMI

```{r, tidy_for_6a1_3, include=F}

# Turn wide format back into long format
foo = melt(setDT(tbw), measure=c("feedbackrmi", "suppressionrmi"), value.name="meanrate", variable.name = "feedback")
tbl = as_tibble(foo)

# Code feedback as binary factor
tbl$feedback = ifelse(tbl$feedback == "feedbackrmi", 1, 0)

```

```{r fit_model_6a1_3}

# Random intercept for neurons,
# random intercept for experiments, nested in series
lmer.6a1.3 = lmer(meanrate ~ feedback + (1 | uid) + (1 | sid/eid), 
                  data = tbl %>% drop_na(meanrate))

display(lmer.6a1.3)
anova(lmer.6a1.3)
```

```{r get_predicted_average_effect_6a1_3, include=F}

m_suppress = fixef(lmer.6a1.3)[1]
diffMeans = fixef(lmer.6a1.3)[2]
m_feedback  = fixef(lmer.6a1.3)[1] + diffMeans
```

## Predicted average effect on firing rate RMI
Feedback: RMI = `r format(m_feedback, digits=2, nsmall=2)` \newline
Suppression: RMI = `r format(m_suppress, digits=2, nsmall=2)` \newline
n = `r nrow(tbl %>% drop_na(meanrate) %>% count(uid))` neurons from `r nrow(tbl %>% drop_na(meanrate) %>% count(mid))` mice

\newpage
# Figure 6a$_2$

```{r tidy_for_6a2, include = FALSE}

# Turn long format into wide format, using data.tables and dcast:
foo = dcast(setDT(tb), mid + sid + eid + uid + mseu ~ mi, value.var = "meanburstratio")
tbw = as_tibble(foo)

# remove outliers
tbw_clean <- tbw %>% filter(suppressionrmi < 0.95 & feedbackrmi < 0.95)
```

## (1) Comparing RMI during suppression against 0

```{r fit_model_6a2_1}

# Fixed effect intercept only, 
# random intercept for neurons,
# random intercept for experiments,nested in series
lmer.6a2.1 = lmer(suppressionrmi ~ 1 + (1 | uid) + (1 | sid/eid), 
                  data = tbw_clean %>% drop_na(suppressionrmi))

display(lmer.6a2.1)
```

## Mean burst ratio RMI
Suppression: RMI = `r format(fixef(lmer.6a2.1)[1], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.6a2.1)[1], digits=1, nsmall=1)` \newline
n = `r nrow(tbw_clean %>% drop_na(suppressionrmi) %>% count(uid))` neurons from `r nrow(tbw_clean %>% drop_na(suppressionrmi) %>% count(mid))` mice

\newpage
# Figure 6a$_2$
## (2) Slope of regression line

```{r fit_model_6a2_2}

# Random intercept for neurons,
# random intercept for experiments, nested in series
lmer.6a2.2 = lmer(feedbackrmi ~ suppressionrmi + (1 | uid) + (1 | sid/eid), 
                  data = tbw_clean %>% drop_na(feedbackrmi, suppressionrmi))

display(lmer.6a2.2)
anova(lmer.6a2.2)
```

``` {r save_coefficients_6a2_2, include=F}

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

Slope = `r format(fixef(lmer.6a2.2)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.6a2.2)[2], digits=1, nsmall=1)` \newline
n = `r nrow(tbw_clean %>% drop_na(feedbackrmi, suppressionrmi) %>% count(uid))` neurons from 
`r nrow(tbw_clean %>% drop_na(feedbackrmi, suppressionrmi) %>% count(mid))` mice

\newpage
# Figure 6a$_2$
## (3) Average effect of feedback on burst ratio RMI

```{r, tidy_for_6a2_3, include=F}

# Turn wide format back into long format
foo = melt(setDT(tbw_clean), measure=c("feedbackrmi", "suppressionrmi"), value.name="meanburstratio", variable.name = "feedback")
tbl = as_tibble(foo)

# Code feedback as binary variable
tbl$feedback = ifelse(tbl$feedback == "feedbackrmi", 1, 0)

```


```{r fit_model_6a2_3}

# Random intercept for neurons,
# random intercept for experiments, nested in series, nested in mice
lmer.6a2.3 = lmer(meanburstratio ~ feedback + (1 | uid) + (1 | mid/sid/eid), 
                  data = tbl %>% drop_na(meanburstratio))

display(lmer.6a2.3)
anova(lmer.6a2.3)

```

```{r get_predicted_average_effect_6a2_3, include=F}
m_suppress = fixef(lmer.6a2.3)[1]
diffMeans = fixef(lmer.6a2.3)[2]
m_feedback  = fixef(lmer.6a2.3)[1] + diffMeans

```

## Predicted average effect on burst ratio RMI

Feedback: RMI = `r format(m_feedback, digits=2, nsmall=2)` \newline
Suppression: RMI = `r format(m_suppress, digits=2, nsmall=2)`\newline
n = `r nrow(tbl %>% drop_na(meanburstratio) %>% count(uid))` neurons from `r nrow(tbl %>% drop_na(meanburstratio) %>% count(mid))` mice

\newpage
# Figure 6a$_3$

```{r tidy_for_6a3, include = FALSE}

# Turn long format into wide format, using data.tables and dcast:
foo = dcast(setDT(tb), mid + sid + eid + uid + mseu ~ mi, value.var = "spars")
tbw = as_tibble(foo)
```

## (1) Comparing RMI during suppression against 0 

```{r, fit_model_6a3_1}

# Fixed effect intercept only,
# random intercept for neurons,
# random intercept for experiments, nested in series
lmer.6a3.1 = lmer(suppressionrmi ~ 1 + (1 | uid) + (1 | sid/eid), 
                  data = tbw %>% drop_na(suppressionrmi))

display(lmer.6a3.1)
```
 
## Mean sparseness RMI
Suppression: RMI = `r format(fixef(lmer.6a3.1)[1], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.6a3.1)[1], digits=1, nsmall=1)` \newline
n = `r nrow(tbw %>% drop_na(suppressionrmi) %>% count(uid))` neurons from `r nrow(tbw %>% drop_na(suppressionrmi) %>% count(mid))` mice

\newpage
# Figure 6a$_3$
## (2) Slope of regression line

```{r fit_model_6a3_2}

# Random intercept for neurons,
# random intercept for experiments, nested in series
lmer.6a3.2 = lmer(feedbackrmi ~ suppressionrmi + (1 | uid) + (1 | sid/eid), 
                  data = tbw %>% drop_na(feedbackrmi, suppressionrmi))

display(lmer.6a3.2)
anova(lmer.6a3.2)
```

``` {r save_coefficients_6a3_2, include=F}

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

## Regression line parameters

Slope of `r format(fixef(lmer.6a3.2)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.6a3.2)[2], digits=2, nsmall=2)` \newline
n = `r nrow(tbw %>% drop_na(feedbackrmi, suppressionrmi) %>% count(uid))` neurons from `r nrow(tbw %>% drop_na(feedbackrmi, suppressionrmi) %>% count(mid))` mice

\newpage
# Figure 6a$_3$
## (3) Average effect of feedback on sparseness RMI

```{r, tidy_for_6a3_3, include=F}

# Turn wide format back into long format
foo = melt(setDT(tbw), measure=c("feedbackrmi", "suppressionrmi"), value.name="spars", variable.name = "feedback")
tbl = as_tibble(foo)

# Code feedback as binary variable
tbl$feedback = ifelse(tbl$feedback == "feedbackrmi", 1, 0)
```


```{r, fit_model_6a3_3}

# Random intercept for neurons,
# random intercept for experiments, nested in series
lmer.6a3.3 = lmer(spars ~ feedback + (1 | uid) + (1 | sid/eid), 
                  data = tbl %>% drop_na(spars))

display(lmer.6a3.3)
anova(lmer.6a3.3)

```

```{r get_predicted_average_effect_6a3_3, include=F}
m_suppress = fixef(lmer.6a3.3)[1]
diffMeans = fixef(lmer.6a3.3)[2]
m_feedback  = fixef(lmer.6a3.3)[1] + diffMeans

```

## Predicted average effect on sparseness RMI

Feedback: Sparseness = `r format(m_feedback, digits=2, nsmall=2)` \newline
Suppression: Sparseness = `r format(m_suppress, digits=2, nsmall=2)` \newline
n = `r nrow(tbl %>% drop_na(spars) %>% count(uid))` neurons from `r nrow(tbl %>% drop_na(spars) %>% count(mid))` mice

\newpage
# Figure 6a$_4$

```{r, tidy_for_6a4, include=F}

# Turn long format into wide format, using data.tables and dcast:
foo = dcast(setDT(tb), mid + sid + eid + uid + mseu ~ mi, value.var = "rel")
tbw = as_tibble(foo)

# remove outliers
tbw_clean <- tbw %>% filter(suppressionrmi > -0.99 & suppressionrmi < 0.99 & feedbackrmi > -0.99 & feedbackrmi < 0.99 )
```

## (1) Comparing RMI during suppression against 0

```{r, fit_model_6a4_1}

# Fixed effect intercept only,
# random intercept for neurons
# random intercept for experiments
lmer.6a4.1 = lmer(suppressionrmi ~ 1 + (1 | uid) + (1 | eid), 
                  data = tbw_clean %>% drop_na(suppressionrmi))

display(lmer.6a4.1)

```

## Mean reliability RMI
Suppression: RMI = `r format(fixef(lmer.6a4.1)[1], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.6a4.1)[1], digits=1, nsmall=1)` \newline
n = `r nrow(tbw_clean %>% drop_na(suppressionrmi) %>% count(uid))` neurons from `r nrow(tbw_clean %>% drop_na(suppressionrmi) %>% count(mid))` mice


\newpage
# Figure 6a$_4$
## (2) Slope of regression line

```{r fit_model_6a4_2}

# Random intercept for neurons,
# random intercept for experiments
lmer.6a4.2 = lmer(feedbackrmi ~ suppressionrmi + (1 | uid) + (1 | eid), 
                  data = tbw_clean %>% drop_na(feedbackrmi, suppressionrmi))

display(lmer.6a4.2)
anova(lmer.6a4.2)
```

``` {r save_coefficients_6a4_2, include=FALSE}

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

Slope of `r format(fixef(lmer.6a4.2)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.6a4.2)[2], digits=2, nsmall=2)` \newline
n = `r nrow(tbw_clean %>% drop_na(feedbackrmi, suppressionrmi) %>% count(uid))` neurons from 
`r nrow(tbw_clean %>% drop_na(feedbackrmi, suppressionrmi) %>% count(mid))` mice

\newpage
# Figure 6a$_4$
## (3) Average effect of feedback on reliability RMI

```{r, tidy_for_6a4_3, include=F}

# Turn wide format back into long format
foo = melt(setDT(tbw_clean), measure=c("feedbackrmi", "suppressionrmi"), value.name="rel", variable.name = "feedback")
tbl = as_tibble(foo)

# Code feedback as binary variable
tbl$feedback = ifelse(tbl$feedback == "feedbackrmi", 1, 0)
```

```{r fit_model_6a4_3}

# Random intercept for neurons,
# random intercept for experiments, nested in series
lmer.6a4.3 = lmer(rel ~ feedback + (1 | uid) + (1 | sid/eid), 
                  data = tbl %>% drop_na(rel))

display(lmer.6a4.3)
anova(lmer.6a4.3)
```

```{r get_predicted_average_effect_6a4_4, include=F}

m_suppress = fixef(lmer.6a4.3)[1]
diffMeans = fixef(lmer.6a4.3)[2]
m_feedback  = fixef(lmer.6a4.3)[1] + diffMeans
```

## Predicted average effect on reliability RMI

Feedback: RMI = `r format(m_feedback, digits=2, nsmall=2)` \newline
Suppression: RMI = `r format(m_suppress, digits=2, nsmall=2)` \newline
n = `r nrow(tbl %>% drop_na(rel) %>% count(uid))` neurons from `r nrow(tbl %>% drop_na(rel) %>% count(mid))` mice

\newpage
# Figure 6b$_1$

```{r tidy_for_6b1_1, include = FALSE}

tb <- tib %>% filter(mi == "sitfmi" | mi == "runfmi")

# Turn long format into wide format, using data.tables and dcast:
foo = dcast(setDT(tb), mid + sid + eid + uid + mseu ~ mi, value.var = "meanrate")
tbw = as_tibble(foo)
```

## (1) Slope of regression line

```{r fit_model_6b1_1}

# Random intercept for neurons,
# random intercept for experiments, nested in series
lmer.6b1_1 = lmer(runfmi ~ sitfmi + (1 | uid) + (1 | sid/eid), 
                  data = tbw %>% drop_na(runfmi, sitfmi))

display(lmer.6b1_1)
anova(lmer.6b1_1)
```

``` {r store_coefficients_6b1_1, include=FALSE}

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

Slope of `r format(fixef(lmer.6b1_1)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.6b1_1)[2], digits=2, nsmall=2)` \newline
n = `r nrow(tbw %>% drop_na(runfmi, sitfmi) %>% count(uid))` neurons from `r nrow(tbw %>% drop_na(runfmi, sitfmi) %>% count(mid))` mice

\newpage
# Figure 6b$_1$
## (2) Average effect of locomotion state on firing rate FMI

```{r, tidy_for_6b1_2, include=F}

# Turn wide format back into long format
foo = melt(setDT(tbw), measure=c("sitfmi", "runfmi"), value.name="meanrate", variable.name = "run")
tbl = as_tibble(foo)

# Code locomotion state as binary variable
tbl$run = ifelse(tbl$run == "runfmi", 1, 0)
```

```{r fit_model_6b1_2}

# Random intercept for neurons,
# random intercept for experiments
lmer.6b1_2 = lmer(meanrate ~ run + (1 | uid) + (1 | eid), 
                  data = tbl %>% drop_na(meanrate))

display(lmer.6b1_2)
anova(lmer.6b1_2)
```

```{r get_predicted_average_effect_6b1_2, include=F}

m_sit = fixef(lmer.6b1_2)[1]
diffMeans = fixef(lmer.6b1_2)[2]
m_run  = fixef(lmer.6b1_2)[1] + diffMeans
```

## Predicted average effect on firing rate FMI

Locomotion: `r format(m_run, digits=2, nsmall=2)` \newline
Quiescence: `r format(m_sit, digits=2, nsmall=2)` \newline
n = `r nrow(tbl %>% drop_na(meanrate) %>% count(uid))` neurons from `r nrow(tbl %>% drop_na(meanrate) %>% count(mid))` mice

\newpage
# Figure 6b$_2$

```{r tidy_for_6b2_1, include = FALSE}

tb <- tib %>% filter(mi == "sitfmi" | mi == "runfmi")

# Turn long format into wide format, using data.tables and dcast:
foo = dcast(setDT(tb), mid + sid + eid + uid + mseu ~ mi, value.var = "meanburstratio")
tbw = as_tibble(foo)

# Remove NaNs
tbw <- tbw %>% filter(sitfmi != "Na" & runfmi != "Na")

# remove two outliers
tbw_clean <- tbw %>% filter(runfmi < 0.99 & sitfmi < 0.99)
```

## (1) Slope of regression line

```{r fit_model_6b2_1}

# Random intercept for neurons,
# random effect for experiments, nested in series
lmer.6b2_1 = lmer(runfmi ~ sitfmi + (1 | uid) + (1 | sid/eid), 
                  data = tbw_clean %>% drop_na(runfmi, sitfmi))

display(lmer.6b2_1)
anova(lmer.6b2_1)
```

```{r store_coefficients_6b2_1, include=FALSE}

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


\textbf{Note: the 2 outliers sitting in the top left and bottom right corner have been exluded before fitting the model!}

Slope of `r format(fixef(lmer.6b2_1)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.6b2_1)[2], digits=2, nsmall=2)` \newline
n = `r nrow(tbw_clean %>% drop_na(runfmi, sitfmi) %>% count(uid))` neurons from `r nrow(tbw_clean %>% drop_na(runfmi, sitfmi) %>% count(mid))` mice

\newpage
# Figure 6b$_2$
## (2) Average effect of locomotion state on burst ratio FMI

```{r, tidy_for_6b2_2, include=F}

# Turn wide format back into long format
foo = melt(setDT(tbw_clean), measure=c("sitfmi", "runfmi"), value.name="meanburstratio", variable.name = "run")
tbl_clean = as_tibble(foo)

# Code locomotion state as binary variable
tbl_clean$run = ifelse(tbl_clean$run == "runfmi", 1, 0)
```


```{r fit_model_6b2_2}

# Random intercept for neurons,
# random intercept for series
lmer.6b2_2 = lmer(meanburstratio ~ run + (1 | uid) + (1 | sid), 
                  data = tbl_clean %>% drop_na(meanburstratio))

display(lmer.6b2_2)
anova(lmer.6b2_2)
```

```{r get_predicted_average_effect_6b2_2, include=F}

m_sit = fixef(lmer.6b2_2)[1]
diffMeans = fixef(lmer.6b2_2)[2]
m_run  = fixef(lmer.6b2_2)[1] + diffMeans
```

## Predicted average effect on burst ratio FMI

\textbf{Note: the 2 outliers sitting in the top left and bottom right corner have been exluded before fitting the model!}

Locomotion: `r format(m_run, digits=2, nsmall=2)` \newline
Quiescence: `r format(m_sit, digits=2, nsmall=2)` \newline
n = `r nrow(tbl_clean %>% drop_na(meanburstratio) %>% count(uid))` neurons from `r nrow(tbl_clean %>% drop_na(meanburstratio) %>% count(mid))` mice

\newpage
# Figure 6b$_3$

```{r tidy_for_6b3_1, include=FALSE}

tb <- tib %>% filter(mi == "sitfmi" | mi == "runfmi")

# Turn long format into wide format, using data.tables and dcast:
foo = dcast(setDT(tb), mid + sid + eid + uid + mseu ~ mi, value.var = "spars")
tbw = as_tibble(foo)
```

## (1) Slope of regression line

```{r fit_model_6b3_1}

# Random intercept for neurons,
# random intercept for experiments, nested in series
lmer.6b3_1 = lmer(runfmi ~ sitfmi + (1 | uid) + (1 | sid/eid), 
                  data = tbw %>% drop_na(runfmi, sitfmi))

display(lmer.6b3_1)
anova(lmer.6b3_1)
```

```{r store_coefficients_6b3_1, include=FALSE}

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

Slope of `r format(fixef(lmer.6b3_1)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.6b3_1)[2], digits=2, nsmall=2)` \newline
n = `r nrow(tbw %>% drop_na(runfmi, sitfmi) %>% count(uid))` neurons from `r nrow(tbw %>% drop_na(runfmi, sitfmi) %>% count(mid))` mice

\newpage
# Figure 6b$_3$
## (2) Average effect of locomotion state on sparseness

```{r, tidy_for_6b3_2, include=F}

# Turn wide format back into long format
foo = melt(setDT(tbw), measure=c("sitfmi", "runfmi"), value.name="spars", variable.name = "run")
tbl = as_tibble(foo)

# Code locomotion state as binary variable
tbl$run = ifelse(tbl$run == "runfmi", 1, 0)
```


```{r fit_model_6b3_2}

# Random intercept for neurons,
# random intercept for experiments, nested in series
lmer.6b3_2 = lmer(spars ~ run + (1 | uid) + (1 | sid/eid), 
                  data = tbl %>% drop_na(spars))

display(lmer.6b3_2)
anova(lmer.6b3_2)

```

```{r get_predicted_average_effect_6b3_2, include=F}
m_sit = fixef(lmer.6b3_2)[1]
diffMeans = fixef(lmer.6b3_2)[2]
m_run  = fixef(lmer.6b3_2)[1] + diffMeans

```

## Predicted average effect sparseness FMI

Quiescence: `r format(m_sit, digits=2, nsmall=2)` \newline
Locomotion: `r format(m_run, digits=2, nsmall=2)` \newline
n = `r nrow(tbl %>% drop_na(spars) %>% count(uid))` neurons from `r nrow(tbl %>% drop_na(spars) %>% count(mid))` mice

\newpage
# Figure 6b$_4$

```{r tidy_for_6b4_1, include = FALSE}

tb <- tib %>% filter(mi == "sitfmi" | mi == "runfmi")

# Turn long format into wide format, using data.tables and dcast:
foo = dcast(setDT(tb), mid + sid + eid + uid + mseu ~ mi, value.var = "rel")
tbw = as_tibble(foo)

# remove outliers
tbw_clean <- tbw %>% filter(runfmi < 0.99 & runfmi > -0.99 & sitfmi < 0.95 & sitfmi > -0.99)
```

## (1) Slope of regression line

```{r fit_model_6b4_1}

# Random intercept for neurons,
# random intercept for experiments
lmer.6b4_1 = lmer(runfmi ~ sitfmi + (1 | uid) + (1 | eid), 
                  data = tbw_clean %>% drop_na(runfmi, sitfmi))

display(lmer.6b4_1)
anova(lmer.6b4_1)
```

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

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

\textbf{Note: Outliers (11 observations) been exluded before fitting the model!}

Slope of `r format(fixef(lmer.6b4_1)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.6b4_1)[2], digits=2, nsmall=2)` \newline
n = `r nrow(tbw_clean %>% drop_na(runfmi, sitfmi) %>% count(uid))` neurons from `r nrow(tbw_clean %>% drop_na(runfmi, sitfmi) %>% count(mid))` mice

\newpage
# Figure 6b$_4$
## (2) Average effect of locomotion state on reliability

```{r, tidy_for_6b4_2, include=F}

# Turn wide format back into long format
foo = melt(setDT(tbw_clean), measure=c("sitfmi", "runfmi"), value.name="rel", variable.name = "run")
tbl_clean = as_tibble(foo)

# Code locomotion state as binary variable
tbl_clean$run = ifelse(tbl_clean$run == "runfmi", 1, 0)
```

```{r fit_model_6b4_2}

# Random intercept for neurons,
# random intercept for experiments, nested in series, nested in mice
lmer.6b4_2 = lmer(rel ~ run + (1 | uid) + (1 | mid/sid/eid), 
                  data = tbl_clean %>% drop_na(rel))

display(lmer.6b4_2)
anova(lmer.6b4_2)
```

```{r get_predicted_average_effect_6b4_2, include=F}
m_sit = fixef(lmer.6b4_2)[1]
diffMeans = fixef(lmer.6b4_2)[2]
m_run  = fixef(lmer.6b4_2)[1] + diffMeans

```

## Predicted average effect reliability FMI

\textbf{Note: Outliers (11 observations) been exluded before fitting the model!}

Quiescence:  `r format(m_sit, digits=2, nsmall=2)` \newline
Locomotion: `r format(m_run, digits=2, nsmall=2)` \newline
n = `r nrow(tbl_clean %>% drop_na(rel) %>% count(uid))` neurons from `r nrow(tbl_clean %>% drop_na(rel) %>% count(mid))` mice

\newpage
# Figure 6c$_1$
## Slope of regression line

```{r tidy_for_6c1, include = FALSE}

tb <- tib %>% filter(mi == "rmi" | mi == "fmi")

# Turn long format into wide format, using data.tables and dcast:
foo = dcast(setDT(tb), mid + sid + eid + uid + mseu ~ mi, value.var = "meanrate")
tbw = as_tibble(foo)

```

```{r fit_model_6c1}

# Random intercept for neurons,
# random intercept for series
lmer.6c1 = lmer(fmi ~ rmi + (1 | uid) + (1 | sid), 
                data = tbw %>% drop_na(fmi, rmi))

display(lmer.6c1)
anova(lmer.6c1)
```

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

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

Slope of `r format(fixef(lmer.6c1)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.6c1)[2], digits=2, nsmall=2)` \newline
n = `r nrow(tbw %>% drop_na(fmi, rmi) %>% count(uid))` neurons from `r nrow(tbw %>% drop_na(fmi, rmi) %>% count(mid))` mice

\newpage
# Figure 6c$_2$
## Slope of regression line

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

tb <- tib %>% filter(mi == "rmi" | mi == "fmi")

# Turn long format into wide format, using data.tables and dcast:
foo = dcast(setDT(tb), mid + sid + eid + uid + mseu ~ mi, value.var = "meanburstratio")
tbw = as_tibble(foo)
```

```{r fit_model_6c2}

# Random intercept for neurons, 
# random intercept for experiments, nested in series
lmer.6c2 = lmer(fmi ~ rmi + (1 | uid) + (1 | sid/eid), 
                data = tbw %>% drop_na(fmi, rmi))

display(lmer.6c2)
anova(lmer.6c2)
```

```{r, store_coefficients_6c2}

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

Slope of `r format(fixef(lmer.6c2)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.6c2)[2], digits=2, nsmall=2)` \newline
n = `r nrow(tbw %>% drop_na(fmi, rmi) %>% count(uid))` neurons from `r nrow(tbw %>% drop_na(fmi, rmi) %>% count(mid))` mice


\newpage
# Figure 6c$_3$
## Slope of regression line

```{r tidy_for_6c3, include = FALSE}

tb <- tib %>% filter(mi == "rmi" | mi == "fmi")

# Turn long format into wide format, using data.tables and dcast:
foo = dcast(setDT(tb), mid + sid + eid + uid + mseu ~ mi, value.var = "spars")
tbw = as_tibble(foo)
```

```{r fit_model_6c3}

# Random intercept for neurons,
# random intercept for experiments, nested in series
lmer.6c3 = lmer(fmi ~ rmi + (1 | uid) + (1 | sid/eid), 
                data = tbw %>% drop_na(fmi, rmi))

display(lmer.6c3)
anova(lmer.6c3)
```

```{r store_coefficients_6c3, include=FALSE}

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

Slope of `r format(fixef(lmer.6c3)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.6c3)[2], digits=2, nsmall=2)` \newline
n = `r nrow(tbw %>% drop_na(fmi, rmi) %>% count(uid))` neurons from `r nrow(tbw %>% drop_na(fmi, rmi) %>% count(mid))` mice

\newpage
# Figure 6c$_4$
## Slope of regression line

```{r tidy_for_6c4, include = FALSE}

tb <- tib %>% filter(mi == "rmi" | mi == "fmi")

# Turn long format into wide format, using data.tables and dcast:
foo = dcast(setDT(tb), mid + sid + eid + uid + mseu ~ mi, value.var = "rel")
tbw = as_tibble(foo)
```

```{r fit_model_6c4}

# Random intercept for neurons,
# random intercept for series, nested in mice
lmer.6c4 = lmer(fmi ~ rmi + (1 | uid) + (1 | mid/sid), 
                data = tbw %>% drop_na(fmi, rmi))

display(lmer.6c4)
anova(lmer.6c4)
```

```{r store_coefficients_6c4, include=FALSE}

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

Slope of `r format(fixef(lmer.6c4)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.6c4)[2], digits=2, nsmall=2)` \newline
n = `r nrow(tbw %>% drop_na(fmi, rmi) %>% count(uid))` neurons from `r nrow(tbw %>% drop_na(fmi, rmi) %>% count(mid))` mice