Copy new R folder from main natfeedback_code repo

R/README.md

@@ -0,0 +1,5 @@
+The R markdown files (`.Rmd`) contain code for all statistical analyses in R, separately for each figure in the manuscript. The output created by fitting linear mixed-effects models is saved in the corresponding `.pdf` files.
+
+The data exported from Python for use as input to the R analyses are stored in `.csv` files in the parent folder's `csv` folder. The `.csv` input files are parsed by a helper function, `get_data.R`, which creates unique IDs for mice, series, experiments, and neurons, required to fit the models.
+
+For easy access, coefficients of the linear models are also saved in *.csv files in the _stats directory.

R/_stats/datasets_figure_1_S6.csv

@@ -0,0 +1,18 @@
+selected_datasets
+PVCre_2017_0006_s03_e03
+PVCre_2017_0006_s03_e04
+PVCre_2017_0006_s03_e05
+PVCre_2017_0008_s09_e04
+PVCre_2017_0015_s03_e05
+PVCre_2017_0015_s03_e07
+PVCre_2018_0001_s05_e03
+PVCre_2018_0001_s05_e04
+PVCre_2018_0001_s05_e05
+PVCre_2018_0003_s02_e03
+PVCre_2018_0003_s03_e03
+PVCre_2017_0015_s03_e04
+PVCre_2018_0001_s02_e02
+PVCre_2018_0003_s02_e02
+PVCre_2018_0003_s03_e02
+PVCre_2019_0002_s08_e03
+PVCre_2019_0002_s08_e07

R/_stats/figure_1_S2c_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+-0.2457854386539932,-0.17986111875873842

R/_stats/figure_1_S2d_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+-0.06223747841067918,-0.6173501501806536

R/_stats/figure_1_S2e_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+-0.06269487224311047,-0.01821533586388691

R/_stats/figure_1_S2f_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+-0.045695622194387574,-0.17727295791250702

R/_stats/figure_1_S2g_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+0.01915848342619288,0.1882790439473608

R/_stats/figure_1_S2i_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+0.05132551908770216,0.8305517796810173

R/_stats/figure_1_S3a_pred_means.csv

@@ -0,0 +1,2 @@
+non_sbc,sbc
+0.20161435670076305,0.0624194164739951

R/_stats/figure_1_S3b_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+0.2348176913104044,-3.147750019901174e-4

R/_stats/figure_1_S3c_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+0.1609107525696011,-0.03413824250103585

R/_stats/figure_1_S3d_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+0.20605334338705786,-0.003460023299540402

R/_stats/figure_1_S3e_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+0.1644391041812277,5.211564992533738e-4

R/_stats/figure_1_S3f_pred_means.csv

@@ -0,0 +1,2 @@
+non_sbc,sbc
+-0.3634075947334936,-0.40321204001616967

R/_stats/figure_1_S3g_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+-0.1684340138455133,-6.747549334719865e-4

R/_stats/figure_1_S3h_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+-0.37049667040931566,-0.05736952846352956

R/_stats/figure_1_S3i_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+-0.18442864711091392,-0.00814730269666998

R/_stats/figure_1_S3j_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+-0.374407649177406,1.0999815229647711

R/_stats/figure_1_S6g_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+0.07925863585442632,1.6471789123627942

R/_stats/figure_1_S6h_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+0.04644824623122267,0.29631085406903324

R/_stats/figure_1_S6i_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+-0.1887857793345154,-3.2649145501015013

R/_stats/figure_1_S6j_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+-0.3058486607963653,-1.2708625585749849

R/_stats/figure_2f_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+0.7002145479673938,-1.287741727070784

R/_stats/figure_3S1j_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+-0.12899792701643514,-1.452777469167875

R/_stats/figure_3_S1a_pred_means.csv

@@ -0,0 +1,2 @@
+non_sbc,sbc
+0.0509409277650741,0.07963562553369906

R/_stats/figure_3_S1b_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+0.05429713915145572,-2.8149815801408974e-6

R/_stats/figure_3_S1c_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+0.036812103372930755,0.11058852822150413

R/_stats/figure_3_S1d_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+0.047165571440243324,-3.731607453157406e-4

R/_stats/figure_3_S1e_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+0.044212121102171216,9.405837704745487e-4

R/_stats/figure_3_S1f_pred_means.csv

@@ -0,0 +1,2 @@
+non_sbc,sbc
+-0.2348290360561482,-0.4869505492936642

R/_stats/figure_3_S1g_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+-0.3701503384003855,4.283080267277868e-4

R/_stats/figure_3_S1h_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+-0.13573493974832399,-0.18021415793983098

R/_stats/figure_3_S1i_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+0.05963345795663311,-0.013030073412489903

R/_stats/figure_4_S1a_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+-0.13684596765847823,0.02857561066748279

R/_stats/figure_4_S1e_pred_means.csv

@@ -0,0 +1,2 @@
+mvi,grt
+-0.33523648349498647,-0.5046560877238382

R/_stats/figure_4_S1f_pred_means.csv

@@ -0,0 +1,2 @@
+grt,grt0c,mvi
+-0.6822892767332962,-0.5842384764304593,-0.6733001529210123

R/_stats/figure_4a_pred_means.csv

@@ -0,0 +1,2 @@
+grt,mvi
+0.052692014000876995,0.14621024628170237

R/_stats/figure_4b_pred_means.csv

@@ -0,0 +1,2 @@
+grt,grt0c,mvi
+0.30222494039645886,0.3576470233298834,0.27223803394804535

R/_stats/figure_5_S1c_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+-0.34162543882112567,0.40563920834914047

R/_stats/figure_5_S1d_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+-0.08612971600787364,-0.11036131665953618

R/_stats/figure_5_S1e_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+-0.14262999310888036,0.5938165913303274

R/_stats/figure_5_S1f_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+-0.13107617771839958,0.5464073341682911

R/_stats/figure_5_S1g_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+0.03173914746079965,0.12340650228886234

R/_stats/figure_5_S1i_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+0.050125039084015055,-0.4614203079725272

R/_stats/figure_6_S1b1_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+0.035808842752692964,0.5158535164912845

R/_stats/figure_6_S1b2_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+-0.2790437623666967,0.518588916898765

R/_stats/figure_6_S1c1_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+0.005919488081091819,0.1771943698966603

R/_stats/figure_6_S1c2_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+-0.3803348746278902,0.2539043923256366

R/_stats/figure_6a1_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+0.13594572766966995,0.513053959410822

R/_stats/figure_6a2_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+-0.1985638572774878,0.37522818202607183

R/_stats/figure_6a3_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+-0.060655497652279296,0.44291431093680106

R/_stats/figure_6a4_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+0.023204972477371223,0.49986997109767084

R/_stats/figure_6b1_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+0.08970968183131361,0.7222281004248075

R/_stats/figure_6b2_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+-0.20775415610743406,0.3443468318242497

R/_stats/figure_6b3_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+-0.005570036539061711,0.8531136959545401

R/_stats/figure_6b4_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+0.014957605204392746,0.4264300616042192

R/_stats/figure_6c1_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+0.1542138178913914,0.053650320172864085

R/_stats/figure_6c2_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+-0.2739246637014166,-0.0988712061764629

R/_stats/figure_6c3_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+-0.15422619151954292,0.004968785469243109

R/_stats/figure_6c4_coefs.csv

@@ -0,0 +1,2 @@
+intercept,slope
+-0.043157942154458166,-0.09511573187769935

R/figure_1.Rmd

@@ -0,0 +1,147 @@
+---
+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
+

R/figure_1_S2.Rmd

@@ -0,0 +1,281 @@
+---
+title: "Spacek et al., 2021, Figure 1-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}
+
+tib = get_data("../csv/fig1.csv")
+```
+
+```{r tidy_for_1_S2ab, include = FALSE}
+
+# Turn booleans for 'optogenetic manipulation' into a binary predictor
+tb <- tib %>% mutate(feedback = ifelse(opto == TRUE, 0, 1))
+
+# Signal-to-noise ratio and mean peak-width are not computed on a trial-by-trial basis. In the csv-file, 
+# these measures are therefore identical across trials, so we simply pull out the first trial of each neuron
+tb = tb %>% select(mid, sid, eid, uid, mseu, feedback, snr, meanpkw) %>% distinct(mseu, feedback, .keep_all = TRUE)
+```
+
+# Figure 1-Supplement 2a
+## Average effect of feedback on signal-to-noise ratio (SNR)
+
+```{r, fit_model_1_S2a}
+
+# We fit a random-intercept model with two random effects: 
+# (1) Neurons (uid) can have different baseline firing rates
+# (2) Mean firing rates are allowed to differ across recording sessions (sid)
+# More complex models with random slopes for neurons (or with experiments nested in
+# sessions, nested in mice) give singular fits.
+lmer.1_S2a = lmer(snr ~ feedback + (1 | uid) + (1 | sid), 
+                  data = tb %>% drop_na(snr))
+
+display(lmer.1_S2a)
+anova(lmer.1_S2a)
+```
+
+```{r get_predicted_average_effect1_S2a, include=F}
+
+mSuppr = fixef(lmer.1_S2a)[1]
+diffMeans = fixef(lmer.1_S2a)[2]
+mActive  = fixef(lmer.1_S2a)[1] + diffMeans
+```
+
+Feedback SNR: `r format(mActive, digits=2, nsmall=2)` \newline
+Suppression SNR: `r format(mSuppr, digits=2, nsmall=2)` \newline
+n = `r nrow(tb %>% drop_na(snr) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(snr) %>% count(mid))` mice
+
+\newpage
+# Figure 1-Supplement 2b
+## Average effect of feedback on PSTH mean peak width
+
+```{r, fit_model_1_S2b}
+
+# Random-intercept for single neurons, 
+# random intercept for experiments, nested in series
+lmer.1_S2b = lmer(meanpkw ~ feedback + (1 | uid) + (1 | sid/eid), 
+                  data = tb %>% drop_na(meanpkw))
+
+display(lmer.1_S2b)
+anova(lmer.1_S2b)
+```
+
+```{r get_predicted_average_effect_1_S2b, include=F}
+
+mSuppr = fixef(lmer.1_S2b)[1]
+diffMeans = fixef(lmer.1_S2b)[2]
+mActive  = fixef(lmer.1_S2b)[1] + diffMeans
+```
+
+Feedback mean peak width: `r format(mActive, digits=2, nsmall=2)` \newline
+Suppression mean peak width: `r format(mSuppr, digits=2, nsmall=2)` \newline
+n = `r nrow(tb %>% drop_na(meanpkw) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(meanpkw) %>% count(mid))` mice
+
+\newpage
+# Figure 1-Supplement 2c
+## Relation between firing rate FMI and burst ratio FMI
+
+```{r read_data_1_S2c-g, include=FALSE}
+
+tib = get_data("../csv/mviFMI.csv")
+
+# filter based on 'state'
+tb <- tib %>% filter(st8 == 'none')
+```
+
+```{r fit_model_1_S2c}
+
+# Random-intercept for single neurons, 
+# random intercept for experiments, nested in series
+lmer.1_S2c = lmer(meanburstratio ~ meanrate + (1 | uid) + (1 | sid/eid), 
+                  data = tb %>% drop_na(meanburstratio, meanrate))
+
+display(lmer.1_S2c)
+anova(lmer.1_S2c)
+```
+
+```{r store_coefficients_1_S2c, include=FALSE}
+
+coef_df = data.frame("intercept" = fixef(lmer.1_S2c)[1], "slope" = fixef(lmer.1_S2c)[2], row.names = "")
+write_csv(coef_df, "_stats/figure_1_S2c_coefs.csv")
+```
+
+Slope of `r format(fixef(lmer.1_S2c)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.1_S2c)[2], digits=2, nsmall=2)` (95\%-confidence interval) \newline
+n = `r nrow(tb %>% drop_na(meanburstratio, meanrate) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(meanburstratio, meanrate) %>% count(mid))` mice
+
+\newpage
+# Figure 1-Supplement 2d
+## Relation between firing rate FMI and sparseness FMI
+
+```{r fit_model_1_S2d}
+
+# Random-intercept for single neurons, 
+# random intercept for experiments, nested in series
+lmer.1_S2d = lmer(spars ~ meanrate + (1 | uid) + (1 | sid/eid), 
+                  data = tb %>% drop_na(spars, meanrate))
+
+display(lmer.1_S2d)
+anova(lmer.1_S2d)
+```
+
+```{r store_coefficients_1_S2d, include=FALSE}
+
+coef_df = data.frame("intercept" = fixef(lmer.1_S2d)[1], "slope" = fixef(lmer.1_S2d)[2], row.names = "")
+write_csv(coef_df, "_stats/figure_1_S2d_coefs.csv")
+```
+
+Slope of `r format(fixef(lmer.1_S2d)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.1_S2d)[2], digits=2, nsmall=2)` (95\%-confidence interval) \newline
+n = `r nrow(tb %>% drop_na(spars, meanrate) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(spars, meanrate) %>% count(mid))` mice
+
+\newpage
+# Figure 1-Supplement 2e
+## Relation between firing rate FMI and reliability FMI
+
+```{r fit_model_1_S2e}
+
+# Random-intercept for single neurons, 
+# random intercept for series, nested in mice
+lmer.1_S2e = lmer(rel ~ meanrate + (1 | uid) + (1 | mid/sid), 
+                  data = tb %>% drop_na(rel, meanrate))
+
+display(lmer.1_S2e)
+anova(lmer.1_S2e)
+```
+
+```{r store_coefficients_1_S2e, include=FALSE}
+
+coef_df = data.frame("intercept" = fixef(lmer.1_S2e)[1], "slope" = fixef(lmer.1_S2e)[2], row.names = "")
+write_csv(coef_df, "_stats/figure_1_S2e_coefs.csv")
+```
+
+Slope of `r format(fixef(lmer.1_S2e)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.1_S2e)[2], digits=2, nsmall=2)` (95\%-confidence interval) \newline
+n = `r nrow(tb %>% drop_na(rel, meanrate) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(rel, meanrate) %>% count(mid))` mice
+
+\newpage
+# Figure 1-Supplement 2f
+## Relation between firing rate FMI and SNR FMI
+
+```{r fit_model_1_S2f}
+
+# Random intercept for neurons, 
+# random intercept for series
+lmer.1_S2f = lmer(snr ~ meanrate + (1 | uid) + (1 | sid), 
+                  data = tb %>% drop_na(snr, meanrate))
+
+display(lmer.1_S2f)
+anova(lmer.1_S2f)
+```
+
+```{r store_coefficients_1_S2f, include=FALSE}
+
+coef_df = data.frame("intercept" = fixef(lmer.1_S2f)[1], "slope" = fixef(lmer.1_S2f)[2], row.names = "")
+write_csv(coef_df, "_stats/figure_1_S2f_coefs.csv")
+```
+
+Slope of `r format(fixef(lmer.1_S2f)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.1_S2f)[2], digits=2, nsmall=2)` (95\%-confidence interval) \newline
+n = `r nrow(tb %>% drop_na(snr, meanrate) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(snr, meanrate) %>% count(mid))` mice
+
+\newpage
+# Figure 1-Supplement 2g
+## Relation between firing rate FMI and mean peak witdth FMI
+
+```{r fit_model_1_S2g}
+
+# Random intercept for neurons,
+# random intercept for experiments, nested in series
+lmer.1_S2g = lmer(meanpkw ~ meanrate + (1 | uid) + (1 | sid/eid), 
+                  data = tb %>% drop_na(meanpkw, meanrate))
+
+display(lmer.1_S2g)
+anova(lmer.1_S2g)
+```
+
+```{r store_coefficients_1_S2g, include=FALSE}
+
+coef_df = data.frame("intercept" = fixef(lmer.1_S2g)[1], "slope" = fixef(lmer.1_S2g)[2], row.names = "")
+write_csv(coef_df, "_stats/figure_1_S2g_coefs.csv")
+```
+
+Slope of `r format(fixef(lmer.1_S2g)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.1_S2g)[2], digits=2, nsmall=2)` (95\%-confidence interval) \newline
+n = `r nrow(tb %>% drop_na(meanpkw, meanrate) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(meanpkw, meanrate) %>% count(mid))` mice
+
+\newpage
+# Figure 1-Supplement 2h
+## Effect of feedback on eye position variability
+
+```{r read_data_1_S2h, include=FALSE}
+
+tib = get_data("../csv/ipos_opto.csv")
+
+# Turn feedback into a binary variable
+tb = tib %>% mutate(feedback = ifelse(opto == "FALSE", 1, 0))
+
+# The standard deviation of eye position is not computed on a trial-by-trial basis. In the csv-file, 
+# std is therefore identical across trials, so we simply pull out the first trial of each neuron
+tb = tb %>% select(mid, sid, eid, mse, feedback, std_xpos_cross) %>% distinct(mse, feedback, .keep_all = TRUE)
+
+```
+
+```{r fit_model_1_S2h}
+
+# Random intercept for experiments, nested in series, nested in mice
+lmer.8 = lmer(std_xpos_cross ~ feedback + (1 | mid/sid/eid), 
+              data = tb %>% drop_na(std_xpos_cross))
+
+display(lmer.8)
+anova(lmer.8)
+
+```
+
+```{r get_predicted_average_effect_1_S2h, include=F}
+mSuppr = fixef(lmer.8)[1]
+diffMeans = fixef(lmer.8)[2]
+mFeedback  = fixef(lmer.8)[1] + diffMeans
+
+```
+
+Mean eye position standard deviation with feedback: `r format(mFeedback, digits=2, nsmall=2)`$^{\circ}$ \newline
+Mean eye position standard deviation with suppression: `r format(mSuppr, digits=2, nsmall=2)`$^{\circ}$  \newline
+n = `r nrow(tb %>% drop_na(std_xpos_cross) %>% count(eid))` experiments from `r nrow(tb %>% drop_na(std_xpos_cross) %>% count(mid))` mice
+
+\newpage
+# Figure 1-Supplement 2i
+## Relation between feedback effects on eye position and feedback effects on reliability
+
+```{r read_data_1_S2i, include=FALSE}
+
+tib = get_data("../csv/iposmi.csv")
+```
+
+```{r fit_model_1S2i}
+
+# Random intercept for neurons,
+# random intercept for experiments nested in series
+lmer.1_S2i = lmer(relfmi ~ iposfmi + (1 | uid) + (1 | sid/eid), 
+                  data = tib %>% drop_na(relfmi, iposfmi))
+
+display(lmer.1_S2i)
+anova(lmer.1_S2i)
+```
+
+```{r store_coefficients, include=FALSE}
+
+coef_df = data.frame("intercept" = fixef(lmer.1_S2i)[1], "slope" = fixef(lmer.1_S2i)[2], row.names = "")
+write_csv(coef_df, "_stats/figure_1_S2i_coefs.csv")
+```
+
+Slope of `r format(fixef(lmer.1_S2i)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.1_S2i)[2], digits=2, nsmall=2)` (95\%-confidence interval) \newline
+n = `r nrow(tib %>% drop_na(relfmi, iposfmi) %>% count(uid))` neurons from `r nrow(tib %>% drop_na(relfmi, iposfmi) %>% count(mid))` mice
+

R/figure_1_S3.Rmd

@@ -0,0 +1,277 @@
+---
+title: "Spacek et al., 2021, Figure 1-Supplement 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/fig1S33S1.csv")
+```
+
+# Figure 1-Supplement 3a
+## Firing rate FMIs separated by cell types
+
+```{r tidy_for_1_S3, include = FALSE}
+
+# Remove observations with NAs
+tb <- tib %>% filter(sbc != "Na")
+
+# Recode suppressed-by-contrast boolean into numeric, binary predictor
+tb$sbc = ifelse(tb$sbc == TRUE, 1, 0)
+```
+
+```{r fit_model_1_S3a}
+
+# Random intercept for series
+lmer.1_S3a = lmer(mvi_meanrate ~ sbc  + (1 | sid), 
+                  data = tb %>% drop_na(mvi_meanrate))
+
+display(lmer.1_S3a)
+anova(lmer.1_S3a)
+```
+
+```{r get_predicted_average_change_1_S3a, include=F}
+
+m_sbc0 = fixef(lmer.1_S3a)[1]
+diffMeans = fixef(lmer.1_S3a)[2]
+m_sbc1  = fixef(lmer.1_S3a)[1] + diffMeans
+```
+
+```{r store_coefficients_1_S3a, include=FALSE}
+
+pred_means_df = data.frame("non_sbc" = m_sbc0, "sbc" = m_sbc1, row.names = "")
+write_csv(pred_means_df, "_stats/figure_1_S3a_pred_means.csv")
+```
+
+FMI sbc: `r format(m_sbc1, digits=3, nsmall=3)` \newline
+FMI non-sbc: `r format(m_sbc0, digits=3, nsmall=3)` \newline
+n = `r nrow(tb %>% drop_na(mvi_meanrate) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(mvi_meanrate) %>% count(mid))` mice
+
+\newpage
+# Figure 1-Supplement 3b
+## Relation between firing rate FMI and recording depth
+
+```{r fit_model_1_S3b}
+
+# Random intercept for series
+lmer.1_S3b = lmer(mvi_meanrate ~ depth  + (1 | sid), 
+                  data = tib %>% drop_na(mvi_meanrate, depth))
+
+display(lmer.1_S3b)
+anova(lmer.1_S3b)
+```
+
+```{r store_coefficients, include=FALSE}
+
+coef_df = data.frame("intercept" = fixef(lmer.1_S3b)[1], "slope" = fixef(lmer.1_S3b)[2], row.names = "")
+write_csv(coef_df, "_stats/figure_1_S3b_coefs.csv")
+```
+
+Slope of `r format(fixef(lmer.1_S3b)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.1_S3b)[2], digits=2, nsmall=2)` (95\%-confidence interval) \newline
+n = `r nrow(tib %>% drop_na(mvi_meanrate, depth) %>% count(uid))` neurons from `r nrow(tib %>% drop_na(mvi_meanrate, depth) %>% count(mid))` mice
+
+\newpage
+# Figure 1-Supplement 3c
+## Relation between firing rate FMI and direction selectivity
+
+``` {r, fit_model_1_S3c}
+
+# Random intercept for series
+lmer.1_S3c = lmer(mvi_meanrate ~ dsi + (1 | sid), 
+                  data = tib %>% drop_na(mvi_meanrate, dsi))
+
+display(lmer.1_S3c)
+anova(lmer.1_S3c)
+```
+
+```{r store_coefficients_1_S3c, include=FALSE}
+
+coef_df = data.frame("intercept" = fixef(lmer.1_S3c)[1], "slope" = fixef(lmer.1_S3c)[2], row.names = "")
+write_csv(coef_df, "_stats/figure_1_S3c_coefs.csv")
+```
+
+Slope of `r format(fixef(lmer.1_S3c)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.1_S3c)[2], digits=2, nsmall=2)` (95\%-confidence interval) \newline
+n = `r nrow(tib %>% drop_na(mvi_meanrate, dsi) %>% count(uid))` neurons from `r nrow(tib %>% drop_na(mvi_meanrate, dsi) %>% count(mid))` mice
+
+\newpage
+# Figure 1-Supplement 3d
+## Relation between firing rate FMI and receptive field location
+
+``` {r, fit_model_1_S3d}
+
+# Random intercept for series nested within mice
+lmer.1_S3d = lmer(mvi_meanrate ~ rfdist + (1 | mid/sid), 
+                  data = tib %>% drop_na(mvi_meanrate, rfdist))
+display(lmer.1_S3d)
+anova(lmer.1_S3d)
+```
+
+```{r store_coefficients_1_S3d, include=FALSE}
+
+# store results in file
+coef_df = data.frame("intercept" = fixef(lmer.1_S3d)[1], "slope" = fixef(lmer.1_S3d)[2], row.names = "")
+write_csv(coef_df, "_stats/figure_1_S3d_coefs.csv")
+
+```
+
+Slope of `r format(fixef(lmer.1_S3d)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.1_S3d)[2], digits=2, nsmall=2)` (95\%-confidence interval) \newline
+n = `r nrow(tib %>% drop_na(mvi_meanrate, rfdist) %>% count(uid))` neurons from `r nrow(tib %>% drop_na(mvi_meanrate, rfdist) %>% count(mid))` mice
+
+\newpage
+# Figure 1-Supplement 3e
+## Relation between firing rate FMI and firing rate
+
+``` {r, fit_model_1_S3e}
+
+# Random intercept for series nested within mice
+lmer.1_S3e = lmer(mvi_meanrate ~ mvi_meanrate_raw + (1 | mid/sid), 
+                  data = tib %>% drop_na(mvi_meanrate, mvi_meanrate_raw))
+
+display(lmer.1_S3e)
+anova(lmer.1_S3e)
+```
+
+```{r store_coefficients_1_S3e, include=FALSE}
+
+coef_df = data.frame("intercept" = fixef(lmer.1_S3e)[1], "slope" = fixef(lmer.1_S3e)[2], row.names = "")
+write_csv(coef_df, "_stats/figure_1_S3e_coefs.csv")
+
+```
+
+Slope of `r format(fixef(lmer.1_S3e)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.1_S3e)[2], digits=2, nsmall=2)` (95\%-confidence interval) \newline
+n = `r nrow(tib %>% drop_na(mvi_meanrate, mvi_meanrate_raw) %>% count(uid))` neurons from `r nrow(tib %>% drop_na(mvi_meanrate, mvi_meanrate_raw) %>% count(mid))` mice
+
+\newpage
+# Figure 1-Supplement 3f
+## Burst ratio FMIs separated by cell types
+
+```{r fit_model_1_S3f}
+
+# Random intercept for series
+lmer.1_S3f = lmer(mvi_meanburstratio ~ sbc  + (1 | sid), 
+                  data = tb %>% drop_na(mvi_meanburstratio))
+
+display(lmer.1_S3f)
+anova(lmer.1_S3f)
+```
+
+```{r get_predicted_average_effect_1_S3f, include=F}
+
+m_sbc0 = fixef(lmer.1_S3f)[1]
+diffMeans = fixef(lmer.1_S3f)[2]
+m_sbc1  = fixef(lmer.1_S3f)[1] + diffMeans
+```
+
+```{r store_coefficients_1_S3f, include=FALSE}
+
+pred_means_df = data.frame("non_sbc" = m_sbc0, "sbc" = m_sbc1, row.names = "")
+write_csv(pred_means_df, "_stats/figure_1_S3f_pred_means.csv")
+```
+
+FMI sbc: `r format(m_sbc1, digits=3, nsmall=3)` \newline
+FMI non-sbc `r format(m_sbc0, digits=3, nsmall=3)` \newline
+n = `r nrow(tb %>% drop_na(mvi_meanburstratio) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(mvi_meanburstratio) %>% count(mid))` mice
+
+\newpage
+# Figure 1-Supplement 3g
+## Relation between burst ratio FMI and recording depth
+
+```{r fit_model_1_S3g}
+
+# Random intercept for series nested within mice
+lmer.1_S3g = lmer(mvi_meanburstratio ~ depth  + (1 | mid/sid), 
+                  data = tib %>% drop_na(mvi_meanburstratio, depth))
+
+display(lmer.1_S3g)
+anova(lmer.1_S3g)
+```
+
+```{r store_coefficients_1_S3g, include=FALSE}
+
+coef_df = data.frame("intercept" = fixef(lmer.1_S3g)[1], "slope" = fixef(lmer.1_S3g)[2], row.names = "")
+write_csv(coef_df, "_stats/figure_1_S3g_coefs.csv")
+```
+
+Slope of `r format(fixef(lmer.1_S3g)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.1_S3g)[2], digits=2, nsmall=2)` (95\%-confidence interval) \newline
+n = `r nrow(tib %>% drop_na(mvi_meanburstratio, depth) %>% count(uid))` neurons from `r nrow(tib %>% drop_na(mvi_meanburstratio, depth) %>% count(mid))` mice
+
+\newpage
+# Figure 1-Supplement 3h
+## Relation between burst ratio FMI and direction selectivity
+
+``` {r, fit_model_1_S3h}
+
+# Random intercept for series nested within mice
+lmer.1_S3h = lmer(mvi_meanburstratio ~ dsi + (1 | mid/sid), 
+                  data = tib %>% drop_na(mvi_meanburstratio, dsi))
+
+display(lmer.1_S3h)
+anova(lmer.1_S3h)
+```
+
+```{r store_coefficients_1_S3h, include=FALSE}
+
+coef_df = data.frame("intercept" = fixef(lmer.1_S3h)[1], "slope" = fixef(lmer.1_S3h)[2], row.names = "")
+write_csv(coef_df, "_stats/figure_1_S3h_coefs.csv")
+```
+
+Slope of `r format(fixef(lmer.1_S3h)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.1_S3h)[2], digits=2, nsmall=2)` (95\%-confidence interval) \newline
+n = `r nrow(tib %>% drop_na(mvi_meanburstratio, dsi) %>% count(uid))` neurons from `r nrow(tib %>% drop_na(mvi_meanburstratio, dsi) %>% count(mid))` mice
+
+\newpage
+# Figure 1-Supplement 3i
+## Relation between burst ratio FMI and receptive field location
+
+``` {r, fit_model_1_S3i}
+
+# Random intercept for series nested within mice
+lmer.1_S3i = lmer(mvi_meanburstratio ~ rfdist + (1 | mid/sid), 
+                  data = tib %>% drop_na(mvi_meanburstratio, rfdist))
+
+display(lmer.1_S3i)
+anova(lmer.1_S3i)
+```
+
+```{r store_coefficients_1_S3i, include=FALSE}
+
+coef_df = data.frame("intercept" = fixef(lmer.1_S3i)[1], "slope" = fixef(lmer.1_S3i)[2], row.names = "")
+write_csv(coef_df, "_stats/figure_1_S3i_coefs.csv")
+```
+
+Slope of `r format(fixef(lmer.1_S3i)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.1_S3i)[2], digits=2, nsmall=2)` (95\%-confidence interval) \newline
+n = `r nrow(tib %>% drop_na(mvi_meanburstratio, rfdist) %>% count(uid))` neurons from `r nrow(tib %>% drop_na(mvi_meanburstratio, rfdist) %>% count(mid))` mice
+
+\newpage
+# Figure 1-Supplement 3j
+## Relation between burst ratio FMI and burst ratio
+
+``` {r, fit_model_1_S3j}
+
+# Random intercept for series
+lmer.1_S3j = lmer(mvi_meanburstratio ~ mvi_meanburstratio_raw + (1 | sid), 
+                  data = tib %>% drop_na(mvi_meanburstratio, mvi_meanburstratio_raw))
+
+display(lmer.1_S3j)
+anova(lmer.1_S3j)
+```
+
+```{r store_coefficients_1_S3j, include=FALSE}
+
+coef_df = data.frame("intercept" = fixef(lmer.1_S3j)[1], "slope" = fixef(lmer.1_S3j)[2], row.names = "")
+write_csv(coef_df, "_stats/figure_1_S3j_coefs.csv")
+```
+
+Slope of `r format(fixef(lmer.1_S3j)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.1_S3j)[2], digits=2, nsmall=2)` (95\%-confidence interval) \newline
+n = `r nrow(tib %>% drop_na(mvi_meanburstratio, mvi_meanburstratio_raw) %>% count(uid))` neurons from `r nrow(tib %>% drop_na(mvi_meanburstratio, mvi_meanburstratio_raw) %>% count(mid))` mice
+

R/figure_1_S4.Rmd

@@ -0,0 +1,239 @@
+---
+title: "Spacek et al., 2021, Figure 1-Supplement 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_1_S4_e-h, include=FALSE}
+
+tib = get_data("../csv/fig1S4mvi.csv")
+```
+
+```{r tidy_for_1_S4ef, include = FALSE}
+
+# Turn feedback into a binary variable
+tb <- tib %>% mutate(feedback = ifelse(opto == TRUE, 0, 1))
+```
+
+# Figure 1-Supplement 4e
+## Effect of suppression on firing rate - movies 
+
+```{r fit_model_1_S4e}
+
+# Random intercept, random slope for neurons,
+# random intercept for experiments
+lmer.1_S4e = lmer(rates ~ feedback + (1 + feedback | uid) + (1 | eid), 
+                  data = tb %>% drop_na(rates))
+
+display(lmer.1_S4e)
+anova(lmer.1_S4e)
+```
+
+```{r get_predicted_average_effect_1_S4e, include=F}
+
+mSuppr = fixef(lmer.1_S4e)[1]
+diffRate = fixef(lmer.1_S4e)[2]
+mActive  = fixef(lmer.1_S4e)[1] + diffRate
+```
+
+Feedback: `r format(mActive, digits=2, nsmall=2)` spikes/s \newline
+Suppression: `r format(mSuppr, digits=2, nsmall=2)` 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 1-Supplement 4f
+## Effect of suppression on burst ratio - movies
+
+```{r fit_model_1_S4f}
+
+# Random intercept, random slope for neurons,
+# random intercept for experiments, nested in series
+lmer.1_S4f = lmer(burstratios ~ feedback + (1 + feedback | uid) + (1 | sid/eid), 
+                  data = tb %>% drop_na(burstratios))
+
+display(lmer.1_S4f)
+anova(lmer.1_S4f)
+```
+
+```{r get_predicted_average_effect_1_S4f, include=F}
+
+mSuppr = fixef(lmer.1_S4f)[1]
+diffRate = fixef(lmer.1_S4f)[2]
+mActive  = fixef(lmer.1_S4f)[1] + diffRate
+```
+
+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 1-Supplement 4g
+## Effect of suppression on sparseness - movies
+
+```{r tidy_for_1_S4gh, 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
+
+tbgh = tb %>% select(mid, sid, eid, uid, mseu, feedback, spars, rel) %>% distinct(mseu, feedback, .keep_all = TRUE)
+
+```
+
+```{r fit_model_1_S4g}
+
+# Random intercept, random slope for single neurons,
+# random intercept for series
+lmer.1_S4g = lmer(spars ~ feedback + (1 + feedback | uid) + (1 | sid), 
+                  data = tbgh %>% drop_na(spars))
+
+display(lmer.1_S4g)
+anova(lmer.1_S4g)
+```
+
+```{r get_predicted_average_change_1_S4g, include=F}
+
+mSuppr = fixef(lmer.1_S4g)[1]
+mDiff = fixef(lmer.1_S4g)[2]
+mFeedback  = fixef(lmer.1_S4g)[1] + mDiff
+```
+
+Feedback: `r format(mFeedback, digits=2, nsmall=2)` \newline
+Suppression: `r format(mSuppr, digits=2, nsmall=2)` \newline
+n = `r nrow(tbgh %>% drop_na(spars) %>% count(uid))` neurons from `r nrow(tbgh %>% drop_na(spars) %>% count(mid))` mice
+
+\newpage
+# Figure 1-Supplement 4h
+## Effect of suppression on reliability - movies
+
+
+```{r fit_model_1_S4h}
+
+# Random intercept, random slope for neurons,
+# random intercept for experiments, nested in mice
+lmer.1_S4h = lmer(rel ~ feedback + (1 +feedback | uid) + (1 | mid/eid), 
+                  data = tbgh %>% drop_na(rel))
+
+display(lmer.1_S4h)
+anova(lmer.1_S4h)
+```
+
+```{r get_predicted_average_effect_1_S4h, include=F}
+
+mSuppr = fixef(lmer.1_S4h)[1]
+mDiff = fixef(lmer.1_S4h)[2]
+mFeedback  = fixef(lmer.1_S4h)[1] + mDiff
+```
+
+Feedback: `r format(mFeedback, digits=2, nsmall=2)` \newline
+Suppression: `r format(mSuppr, digits=2, nsmall=2)` \newline
+n = `r nrow(tbgh %>% drop_na(rel) %>% count(uid))` neurons from `r nrow(tbgh %>% drop_na(rel) %>% count(mid))` mice
+
+
+```{r read_data_1_S4_l-o, include=FALSE}
+
+tib = get_data("../csv/fig1S4grt.csv")
+```
+
+```{r tidy_for_1_S4_l-o, include = FALSE}
+
+# Turn feedback into a binary variable
+tb <- tib %>% mutate(feedback = ifelse(opto == TRUE, 0, 1))
+```
+
+\newpage
+# Figure 1-Supplement 4l
+## Effect of suppression on firing rate - gratings
+
+```{r fit_model_1_S4l}
+
+# Random intercept, random slope for neurons,
+# random intercept for experiments nested in series
+lmer.1_S4l = lmer(rates ~ feedback + (1 + feedback | uid) + (1 | sid/eid), 
+                  data = tb %>% drop_na(rates))
+
+display(lmer.1_S4l)
+anova(lmer.1_S4l)
+```
+
+```{r get_predicted_average_effect_1_S4l, include=F}
+
+mSuppr = fixef(lmer.1_S4l)[1]
+diffRate = fixef(lmer.1_S4l)[2]
+mActive  = fixef(lmer.1_S4l)[1] + diffRate
+```
+
+Feedback: `r format(mActive, digits=2, nsmall=2)` spikes/s \newline
+Suppression: `r format(mSuppr, digits=2, nsmall=2)` 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 1-Supplement 4m
+## Effect of suppression on burst ratio - gratings
+
+```{r fit_model_1_S4m}
+
+# Random intercept, random slope for neurons,
+# random intercept for experiments
+lmer.1_S4m = lmer(burstratios ~ feedback + (1 + feedback | uid) + (1 | eid), 
+                  data = tb %>% drop_na(burstratios))
+
+display(lmer.1_S4m)
+anova(lmer.1_S4m)
+```
+
+```{r predicted_average_effect_1_S4m, include=F}
+
+mSuppr = fixef(lmer.1_S4m)[1]
+diffRate = fixef(lmer.1_S4m)[2]
+mActive  = fixef(lmer.1_S4m)[1] + diffRate
+```
+
+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 1-Supplement 4n
+## Effect of suppression on F1/F0 - gratings
+
+```{r tidy_for_1_S4no, include=FALSE}
+
+# 'F1/F0 ratio' is not computed on a trial-by-trial basis. In the csv-file, 
+# this measure is therefore identical across trials, so we simply pull out the first trial of each neuron
+
+tbn = tb %>% select(mid, sid, eid, uid, mseu, feedback, f1f0) %>% distinct(mseu, feedback, .keep_all = TRUE)
+
+```
+
+```{r fit_model_1_S4n}
+
+# Random intercept for neurons,
+# random intercept for experiments, nested in series
+lmer.1_S4n = lmer(f1f0 ~ feedback + (1 | uid) + (1 | sid/eid), 
+                  data = tbn %>% drop_na(f1f0))
+
+display(lmer.1_S4n)
+anova(lmer.1_S4n)
+
+```
+
+```{r get_predicted_average_effect_1_S4n, include=F}
+
+mSuppr = fixef(lmer.1_S4n)[1]
+mDiff = fixef(lmer.1_S4n)[2]
+mFeedback  = fixef(lmer.1_S4n)[1] + mDiff
+```
+
+Feedback: `r format(mFeedback, digits=2, nsmall=2)` \newline
+Suppression: `r format(mSuppr, digits=2, nsmall=2)` \newline
+n = `r nrow(tbn %>% drop_na(f1f0) %>% count(uid))` neurons from `r nrow(tbn %>% drop_na(f1f0) %>% count(mid))` mice

R/figure_1_S5.Rmd

@@ -0,0 +1,137 @@
+---
+title: "Spacek et al., 2021, Figure 1-Supplement 5"
+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_1_S5_c-f, include=FALSE}
+
+tib = get_data("../csv/fig1S5mvi.csv")
+```
+
+```{r tidy_for_1_S5cd, include = FALSE}
+
+# Turn feedback into a binary variable
+tb <- tib %>% mutate(feedback = ifelse(opto == TRUE, 0, 1))
+```
+
+# Figure 1-Supplement 5c
+## Effect of suppression on firing rate - movies
+
+```{r fit_model_1_S5c}
+
+# Random intercept for neurons - including random slope gives singular fit
+lmer.1_S5c = lmer(rates ~ feedback + (1 | uid), 
+                  data = tb %>% drop_na(rates))
+
+display(lmer.1_S5c)
+anova(lmer.1_S5c)
+```
+
+```{r get_predicted_average_effect_1_S5c, include=F}
+
+mSuppr = fixef(lmer.1_S5c)[1]
+diffRate = fixef(lmer.1_S5c)[2]
+mActive  = fixef(lmer.1_S5c)[1] + diffRate
+```
+
+Feedback: `r format(mActive, digits=2, nsmall=2)` spikes/s \newline
+Suppression: `r format(mSuppr, digits=2, nsmall=2)` spikes/s  \newline
+n = `r nrow(tb %>% drop_na(rates) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(rates) %>% count(mid))` mouse
+
+\newpage
+# Figure 1-Supplement 5d
+## Effect of suppression on burst ratio - movies
+
+```{r fit_model_1_S5d}
+
+# Random intercept for neurons, including random slope gives singular fits
+lmer.1_S5d = lmer(burstratios ~ feedback + (1 | uid), 
+                  data = tb %>% drop_na(burstratios))
+
+display(lmer.1_S5d)
+anova(lmer.1_S5d)
+```
+
+```{r get_predicted_average_effect_1_S5d, include=F}
+
+mSuppr = fixef(lmer.1_S5d)[1]
+diffRate = fixef(lmer.1_S5d)[2]
+mActive  = fixef(lmer.1_S5d)[1] + diffRate
+```
+
+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))` mouse
+
+
+```{r read_data_1_S5_h-i, include=FALSE}
+
+tib = get_data("../csv/fig1S5grt.csv")
+```
+
+```{r tidy_for_1_S5_h-i, include = FALSE}
+
+# Turn feedback into a binary variable
+tb <- tib %>% mutate(feedback = ifelse(opto == TRUE, 0, 1))
+```
+
+\newpage
+# Figure 1-Supplement 5h
+## Effect of suppression on firing rates - gratings
+
+```{r fit_model_1_S5h}
+
+# Random intercept for neurons
+lmer.1_S5h = lmer(rates ~ feedback + (1 | uid), 
+                  data = tb %>% drop_na(rates))
+
+display(lmer.1_S5h)
+anova(lmer.1_S5h)
+```
+
+```{r get_predicted_average_effect_1_S5h, include=F}
+
+mSuppr = fixef(lmer.1_S5h)[1]
+diffRate = fixef(lmer.1_S5h)[2]
+mActive  = fixef(lmer.1_S5h)[1] + diffRate
+```
+
+Feedback: `r format(mActive, digits=2, nsmall=2)` spikes/s \newline
+Suppression: `r format(mSuppr, digits=2, nsmall=2)` spikes/s \newline
+n = `r nrow(tb %>% drop_na(rates) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(rates) %>% count(mid))` mouse
+
+\newpage
+# Figure 1-Supplement 5i
+## Effect of suppression on burst ratio - gratings
+
+```{r fit_model_1_S5i}
+
+# Random intercept for neurons
+lmer.1_S5i = lmer(burstratios ~ feedback + (1 | uid), 
+                  data = tb %>% drop_na(burstratios))
+
+display(lmer.1_S5i)
+anova(lmer.1_S5i)
+```
+
+```{r predicted_average_effect_1_S5i, include=F}
+
+mSuppr = fixef(lmer.1_S5i)[1]
+diffRate = fixef(lmer.1_S5i)[2]
+mActive  = fixef(lmer.1_S5i)[1] + diffRate
+```
+
+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))` mouse

R/figure_1_S6.Rmd

@@ -0,0 +1,310 @@
+---
+title: "Spacek et al., 2021, Figure 1-Supplement 6"
+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_pupil_diam, include=FALSE}
+
+# Read data
+tib_pupil = get_data("../csv/ipos_opto.csv")
+```
+
+```{r tidy_pupil_diam, include = FALSE}
+
+# Turn booleans for 'optogenetic manipulation' into a binary predictor
+tb_pupil <- tib_pupil %>% mutate(light_on = ifelse(opto == TRUE, 1, 0)) %>% drop_na(area_trialmean)
+```
+
+```{r do_stats_pupil_diam, include=FALSE}
+
+# On each data set (i.e., experiment) we run a 2-sample Kolmogorov-Smirnov test
+allds = unique(tb_pupil$mse)
+
+allstatistics = rep(NA, length(allds)) # KS statistic D
+allps = rep(NA, length(allds)) # p-values
+i = 1
+for (ids in allds) {
+  
+  df = tb_pupil %>% filter(mse == ids)
+  res = ks.test(df$area_trialmean[df$light_on == 0], df$area_trialmean[df$light_on == 1])
+  if (res$p.value >= 0.05) {
+    print(paste(ids, ': p =', format(res$p.value, digits=2, nsmall=2)), quote = FALSE)
+  }
+  allps[i] = res$p.value
+  allstatistics[i] = res$statistic
+  i = i + 1
+}
+
+# Put together an array of those experiments that have comparable pupil diameters (valid experiments)
+validExps = allds[allps >= 0.05]
+```
+
+For each experiment, we tested whether V1 suppression (i.e., turning on blue light) would affect pupil diameter.
+To do this, we performed a 2-sample Kolmogorov-Smirnov Test, comparing distributions of pupil diamter from trials with versus without V1 suppression.
+These are the experiments (n = `r length(validExps)` out of `r length(allds)`) with comparable distributions, across trials, of pupil diameter:
+
+```{r print_results, echo=FALSE}
+
+writeLines(sprintf("%s, D = %1.2f, p = %1.3f", validExps, allstatistics[allps >= 0.05], allps[allps >= 0.05]))
+
+# save to file
+validExps_df = data.frame("selected_datasets" = validExps)
+write_csv(validExps_df, "_stats/datasets_figure_1_S6.csv")
+
+```
+
+```{r read_neural_data, include=FALSE}
+
+# Read data
+tib = get_data("../csv/fig1.csv")
+```
+
+```{r tidy_neural_data, include = FALSE}
+
+# Turn booleans for 'optogenetic manipulation' into a binary predictor
+tb <- tib %>% mutate(feedback = ifelse(opto == TRUE, 0, 1))
+
+```
+
+\newpage
+# Figure 1-Supplement 6c
+## Feedback effects on firing rate
+
+```{r extract_valid_experiments, include=FALSE}
+  
+  tb_matched = filter(tb, grepl(paste(validExps, collapse="|"), mseu))
+```
+
+```{r fit_model_1_S6c}
+
+# We cannot simply repeat the identical analysis as for Figure 1f, 
+# because with this reduced data set, that model doesn't converge.
+#
+# Without the random intercept for mice, however, the model converges - so here we fit:
+# Random intercept, random slope for neurons,
+# random intercept for experiments, nested in series
+lmer.1_S6c = lmer(rates ~ feedback + (1 + feedback | uid) + (1 | sid/eid), 
+               data = tb_matched %>% drop_na(rates))
+
+
+display(lmer.1_S6c)
+anova(lmer.1_S6c)
+```
+
+```{r get_predicted_average_effect_1f, include=F}
+
+mSuppr = fixef(lmer.1_S6c)[1]
+diffMeans = fixef(lmer.1_S6c)[2]
+mFeedbk  = fixef(lmer.1_S6c)[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_matched %>% drop_na(rates) %>% count(uid))` neurons from `r nrow(tb_matched %>% drop_na(rates) %>% count(mid))` mice
+
+\newpage
+# Figure 1-Supplement 6d
+## Feedback effects on burst ratio
+
+```{r fit_model_1_S6d}
+
+# Random-intercept, random-slope for single neurons, 
+# random intercept for experiments
+lmer.1_S6d = lmer(burstratios ~ feedback + (1 + feedback | uid) + (1 | eid), 
+               data = tb_matched %>% drop_na(burstratios))
+
+display(lmer.1_S6d)
+anova(lmer.1_S6d)
+```
+
+```{r get_predicted_average_effect_1_S6d, include=F}
+
+mSuppr = fixef(lmer.1_S6d)[1]
+diffMeans = fixef(lmer.1_S6d)[2]
+mFeedbk  = fixef(lmer.1_S6d)[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_matched %>% drop_na(burstratios) %>% count(uid))` neurons from `r nrow(tb_matched %>% drop_na(burstratios) %>% count(mid))` mice
+
+\newpage
+# Figure 1-Supplement 6e
+## Feedback effects on sparseness
+
+```{r tidy_for_1_S6ef, 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
+
+tb_matched_ef = tb_matched %>% select(mid, sid, eid, uid, mseu, feedback, spars, rel) %>% distinct(mseu, feedback, .keep_all = TRUE)
+```
+
+
+```{r fit_model_1_S6e}
+
+# Random-intercept, random-slope for single neurons, 
+# random intercept for experiments, nested within series
+lmer.1_S6e = lmer(spars ~ feedback + (1 + feedback | uid) + (1 | sid/eid), 
+               data = tb_matched_ef %>% drop_na(spars))
+
+display(lmer.1_S6e)
+anova(lmer.1_S6e)
+```
+
+```{r get_predicted_average_effect_1_S6e, include=F}
+
+mSuppr = fixef(lmer.1_S6e)[1]
+diffMeans = fixef(lmer.1_S6e)[2]
+mFeedbk  = fixef(lmer.1_S6e)[1] + diffMeans
+```
+
+Feedback: `r format(mFeedbk, digits=2, nsmall=2)` \newline
+Suppression: `r format(mSuppr, digits=2, nsmall=2)` \newline
+n = `r nrow(tb_matched_ef %>% drop_na(spars) %>% count(uid))` neurons from `r nrow(tb_matched_ef %>% drop_na(spars) %>% count(mid))` mice
+
+\newpage
+# Figure 1-Supplement 6f
+## Feedback effects on reliability
+
+```{r fit_model_1S6_f}
+
+# Random-intercept, random-slope for single neurons, 
+# random intercept for experiments, nested within series
+lmer.1_S6f = lmer(rel ~ feedback + (1 + feedback | uid) + (1 | sid/eid), 
+               data = tb_matched_ef %>% drop_na(rel))
+
+display(lmer.1_S6f)
+anova(lmer.1_S6f)
+```
+
+```{r get_predicted_average_effect_1i, include=F}
+mSuppr = fixef(lmer.1_S6f)[1]
+diffMeans = fixef(lmer.1_S6f)[2]
+mFeedbk  = fixef(lmer.1_S6f)[1] + diffMeans
+
+```
+
+Feedback: `r format(mFeedbk, digits=2, nsmall=2)` \newline
+Suppression: `r format(mSuppr, digits=2, nsmall=2)` \newline
+n = `r nrow(tb_matched_ef %>% drop_na(rel) %>% count(uid))` neurons from `r nrow(tb_matched_ef %>% drop_na(rel) %>% count(mid))` mice
+
+\newpage
+# Figure 1-Supplement 6g
+## Relation between pupil area FMI and firing rate FMI
+
+```{r read_data_1_S6_gi, include=FALSE}
+
+# Read data
+tib = get_data("../csv/iposmi.csv")
+```
+
+```{r fit_model_1_S6g}
+
+# Random intercept for neurons,
+# random intercept for experiments
+lmer.1_S6g = lmer(meanratefmi ~ areafmi + (1 | uid) + (1 | eid), 
+                  data = tib %>% drop_na(meanratefmi, areafmi))
+
+display(lmer.1_S6g)
+anova(lmer.1_S6g)
+```
+
+```{r, store_coefficients_1_S6g, include=FALSE}
+
+coef_df = data.frame("intercept" = fixef(lmer.1_S6g)[1], "slope" = fixef(lmer.1_S6g)[2], row.names = "")
+write_csv(coef_df, "_stats/figure_1_S6g_coefs.csv")
+```
+
+Slope of `r format(fixef(lmer.1_S6g)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.1_S6g)[2], digits=2, nsmall=2)` \newline
+n = `r nrow(tib %>% drop_na(areafmi, meanratefmi) %>% count(uid))` neurons from `r nrow(tib %>% drop_na(areafmi, meanratefmi) %>% count(mid))` mice
+
+\newpage
+# Figure 1-Supplement 6i
+## Relation between pupil area FMI and burst ratio FMI
+
+```{r fit_model_1_S6i}
+
+# Random intercept for neurons,
+# random intercept for experiments, nested in series, nested in mice
+lmer.1_S6i = lmer(meanburstratiofmi ~ areafmi + (1 | uid) + (1 | mid/sid/eid), 
+                  data = tib %>% drop_na(meanburstratiofmi, areafmi))
+
+display(lmer.1_S6i)
+anova(lmer.1_S6i)
+```
+
+```{r, store_coefficients_1_S6i, include=FALSE}
+
+coef_df = data.frame("intercept" = fixef(lmer.1_S6i)[1], "slope" = fixef(lmer.1_S6i)[2], row.names = "")
+write_csv(coef_df, "_stats/figure_1_S6i_coefs.csv")
+```
+
+Slope of `r format(fixef(lmer.1_S6i)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.1_S6i)[2], digits=2, nsmall=2)` \newline
+n = `r nrow(tib %>% drop_na(areafmi, meanburstratiofmi) %>% count(uid))` neurons from `r nrow(tib %>% drop_na(areafmi, meanburstratiofmi) %>% count(mid))` mice
+
+\newpage
+# Figure 1-Supplement 6h
+## Relation between pupil area FMI and firing rate FMI (Ntsr1-Cre)
+
+```{r read_data_1_S6_hj, include=FALSE}
+
+# Read data
+tib = get_data("../csv/ntsrmvis/iposmi.csv")
+```
+
+```{r fit_model_1_S6h}
+
+# Random intercept for neurons,
+# random intercept for experiments
+lmer.1_S6h = lmer(meanratefmi ~ areafmi + (1 | uid), 
+                  data = tib %>% drop_na(meanratefmi, areafmi))
+
+display(lmer.1_S6h)
+anova(lmer.1_S6h)
+```
+
+```{r, store_coefficients_1_S6h, include=FALSE}
+
+coef_df = data.frame("intercept" = fixef(lmer.1_S6h)[1], "slope" = fixef(lmer.1_S6h)[2], row.names = "")
+write_csv(coef_df, "_stats/figure_1_S6h_coefs.csv")
+```
+
+Slope of `r format(fixef(lmer.1_S6h)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.1_S6h)[2], digits=2, nsmall=2)` \newline
+n = `r nrow(tib %>% drop_na(areafmi, meanratefmi) %>% count(uid))` neurons from `r nrow(tib %>% drop_na(areafmi, meanratefmi) %>% count(mid))` mice
+
+\newpage
+# Figure 1-Supplement 6j
+## Relation between pupil area FMI and burst ratio FMI (Ntsr1-Cre)
+
+```{r fit_model_1_S6j}
+
+# Random intercept for neurons,
+# random intercept for experiments, nested in series, nested in mice
+lmer.1_S6j = lmer(meanburstratiofmi ~ areafmi + (1 | uid) + (1 | sid/eid), 
+                  data = tib %>% drop_na(meanburstratiofmi, areafmi))
+
+display(lmer.1_S6j)
+anova(lmer.1_S6j)
+```
+
+```{r, store_coefficients_1_S6j, include=FALSE}
+
+coef_df = data.frame("intercept" = fixef(lmer.1_S6j)[1], "slope" = fixef(lmer.1_S6j)[2], row.names = "")
+write_csv(coef_df, "_stats/figure_1_S6j_coefs.csv")
+```
+
+Slope of `r format(fixef(lmer.1_S6j)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.1_S6j)[2], digits=2, nsmall=2)` \newline
+n = `r nrow(tib %>% drop_na(areafmi, meanburstratiofmi) %>% count(uid))` neurons from `r nrow(tib %>% drop_na(areafmi, meanburstratiofmi) %>% count(mid))` mice
+

R/figure_2.Rmd

@@ -0,0 +1,258 @@
+---
+title: "Spacek et al., 2021, Figure 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}
+
+tib = get_data("../csv/fig2.csv")
+```
+
+# Figure 2e
+
+``` {r tidy_2e, include=FALSE}
+
+tb_all = tib %>% filter(fmode == "all" )
+```
+
+## (I) Fit intercept-only model to slope (all spikes)
+
+``` {r fit_model_2e_slope}
+
+# Random intercept for neurons,
+# random intercept for experiments, nested in series
+lmer.2e_slope = lmer(slope ~ 1 + (1 | uid) + (1 | sid/eid), 
+                     data = tb_all %>% drop_na(slope))
+
+display(lmer.2e_slope)
+```
+
+```{r, include=FALSE}
+
+print(paste("Intercept: ", format(fixef(lmer.2e_slope), digits=2), " [", format(fixef(lmer.2e_slope) - 2 * se.fixef(lmer.2e_slope), digits=2), ", ", format(fixef(lmer.2e_slope) + 2 * se.fixef(lmer.2e_slope), digits=2), "]", sep=""))
+```
+## 95-% confidence interval on slope
+
+`r format(fixef(lmer.2e_slope)[1], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.2e_slope)[1], digits=2, nsmall=2)` \newline
+n = `r nrow(tb_all %>% drop_na(slope) %>% count(uid))` neurons from `r nrow(tb_all %>% drop_na(slope) %>% count(mid))` mice
+
+\newpage
+## (II) Fit intercept-only model to threshold (all spikes)
+
+```{r, fit_model_2e_threshold}
+
+# Random intercept for neurons,
+# random intercept for experiments, nested in mice
+lmer.2e_thresh = lmer(thresh ~ 1 + (1 | uid) + (1 | mid/eid), 
+                      data = tb_all %>% drop_na(thresh))
+
+display(lmer.2e_thresh)
+```
+
+```{r, include=FALSE}
+
+print(paste("Intercept: ", format(fixef(lmer.2e_thresh), digits=2), " [", format(fixef(lmer.2e_thresh) - 2 * se.fixef(lmer.2e_thresh) , digits=2), ", ", format(fixef(lmer.2e_thresh) + 2 * se.fixef(lmer.2e_thresh), digits=2), "]", sep=""))
+```
+## 95-% confidence interval on threshold
+
+`r format(fixef(lmer.2e_thresh)[1], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.2e_thresh)[1], digits=2, nsmall=2)` \newline
+n = `r nrow(tb_all %>% drop_na(thresh) %>% count(uid))` neurons from `r nrow(tb_all %>% drop_na(thresh) %>% count(mid))` mice
+
+\newpage
+# Figure 2f
+## Goodness-of-fit vs burst ratio during suppression
+
+``` {r, fit_model_2f}
+
+# Random intercept for neurons,
+# random intercept for experiments, nested in series
+lmer.2f = lmer(rsq ~ suppression_meanburstratio + (1 | uid) + (1 | sid/eid), 
+               data = tb_all %>% drop_na(rsq, suppression_meanburstratio))
+
+display(lmer.2f)
+anova(lmer.2f)
+```
+
+## 95-% confidence interval on slope
+
+`r format(fixef(lmer.2f)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.2f)[2], digits=2, nsmall=2)` \newline
+n = `r nrow(tb_all %>% drop_na(rsq, suppression_meanburstratio) %>% count(uid))` neurons from `r nrow(tb_all %>% drop_na(rsq, suppression_meanburstratio) %>% count(mid))` mice
+
+```{r, store_coeffs_2f, include=F}
+
+# store results in file
+coef_df = data.frame("intercept" = fixef(lmer.2f)[1], "slope" = fixef(lmer.2f)[2], row.names = "")
+write_csv(coef_df, "_stats/figure_2f_coefs.csv")
+```
+
+\newpage
+# Figure 2g
+## Goodness-of-fit with and without removal of burst spikes
+
+``` {r tidy_2g, include = F}
+
+# Extract two firing modes: all spikes, and non-burst spikes, and turn them into numeric binary predictors
+tb <- tib %>% filter(fmode == "all" | fmode == 'nonburst') %>% mutate(allSpikes = ifelse(fmode == "all", 1, 0))
+```
+
+``` {r fit_model_2g}
+
+# Random intercept for neurons,
+# random intercept for experiments, nested in series
+lmer.2g = lmer(rsq ~ allSpikes + (1 | uid) + (1 | sid/eid), 
+               data = tb %>% drop_na(rsq, allSpikes))
+
+display(lmer.2g)
+anova(lmer.2g)
+```
+
+
+```{r get_predicted_average_effect_2g, include=F}
+
+mNonBurst = fixef(lmer.2g)[1]
+diffMeans = fixef(lmer.2g)[2]
+mAll = fixef(lmer.2g)[1] + diffMeans
+```
+
+# Predicted average effect
+
+All spikes: $R^2$ = `r format(mAll, digits=2, nsmall=2)` \newline
+non-burst spikes: $R^2$ = `r format(mNonBurst, digits=2, nsmall=2)` \newline
+n = `r nrow(tb %>% drop_na(rsq, allSpikes) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(rsq, allSpikes) %>% count(mid))` mice
+
+
+\newpage
+# Figure 2h
+## Goodness of fit with and without removal of tonic spikes
+
+``` {r tidy_2h, include=F}
+
+# Extract two firing modes: all spikes, and non-random spikes
+tb <- tib %>% filter(fmode == "all" | fmode == 'nonrand') %>% mutate(allSpikes = ifelse(fmode == "all", 1, 0))
+```
+
+
+``` {r fit_model_2h}
+
+# Random intercept for neurons,
+# random intercept for experiments, nested in series
+lmer.2h = lmer(rsq ~ allSpikes + (1 | uid) + (1 | sid/eid), 
+               data = tb %>% drop_na(rsq, allSpikes))
+
+display(lmer.2h)
+anova(lmer.2h)
+```
+
+
+```{r get_predicted_average_effect_2h, include=F}
+
+mNonRand = fixef(lmer.2h)[1]
+diffMeans = fixef(lmer.2h)[2]
+mAllSpikes = fixef(lmer.2h)[1] + diffMeans
+```
+
+# Predicted average effect
+
+All spikes: $R^2$ = `r format(mAllSpikes, digits=2, nsmall=2)` \newline
+Randomly removed spikes: $R^2$ = `r format(mNonRand, digits=2, nsmall=2)` \newline
+n = `r nrow(tb %>% drop_na(rsq, allSpikes) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(rsq, allSpikes) %>% count(mid))` mice
+
+\newpage
+# Figure 2i
+
+``` {r tidy_2i, include=FALSE}
+
+tb_nonburst = tib %>% filter(fmode == "nonburst")
+```
+
+## (I) Fit intercept-only model to slope (non-burst spikes)
+
+```{r, fit_model_2i_slope}
+
+# Random intercept for neurons,
+# random intercepts for experiments, nested in series
+lmer.2i_slope = lmer(slope ~ 1 + (1 | uid) + (1 | sid/eid), 
+                     data = tb_nonburst %>% drop_na(slope))
+
+display(lmer.2i_slope)
+```
+
+## 95-% confidence interval
+
+```{r, include=FALSE}
+
+print(paste("Intercept: ", format(fixef(lmer.2i_slope), digits=2), " [", format(fixef(lmer.2i_slope) - 2 * se.fixef(lmer.2i_slope), digits=2), ", ", format(fixef(lmer.2i_slope) + 2 * se.fixef(lmer.2i_slope), digits=2), "]", sep=""))
+```
+
+Intercept of `r format(fixef(lmer.2i_slope)[1], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.2i_slope)[1], digits=1, nsmall=1)` \newline
+n = `r nrow(tb_nonburst %>% drop_na(slope) %>% count(uid))` neurons from `r nrow(tb_nonburst %>% drop_na(slope) %>% count(mid))` mice
+
+\newpage
+## (II) Fit intercept-only model to threshold (non-burst spikes)
+
+```{r, fit_model_2i_threshold}
+
+# Random intercept for neurons,
+# random intercept for experiments, nested in series, nested in mice
+lmer.2i_thresh = lmer(thresh ~ 1 + (1 | uid) + (1 | mid/sid/eid), 
+                      data = tb_nonburst %>% drop_na(thresh))
+
+display(lmer.2i_thresh)
+```
+
+## 95-% confidence interval
+
+```{r, include=FALSE}
+
+print(paste("Intercept: ", format(fixef(lmer.2i_thresh), digits=2), " [", format(fixef(lmer.2i_thresh) - 2 * se.fixef(lmer.2i_thresh), digits=2), ", ", format(fixef(lmer.2i_thresh) + 2 * se.fixef(lmer.2i_thresh), digits=2), "]", sep=""))
+```
+
+Intercept of `r format(fixef(lmer.2i_thresh)[1], digits=1, nsmall=1)` $\pm$ `r format(2 * se.fixef(lmer.2i_thresh)[1], digits=1, nsmall=1)`
+
+\newpage
+# Fit a single model to Figures 2e and 2i
+
+## (I) Predict slope based on firing mode (all spikes vs non-burst spikes)
+
+```{r, tidy_2ei, include=FALSE}
+
+# (1) Fit a single model to panels e and i to examine changes in slope with firing mode as a binary predictor. We need to turn the strings (= factors, = categorical predictors) into binary predictors (numeric) first:
+tb = tib %>% filter(fmode == "all" | fmode == "nonburst") %>% mutate(allSpikes = ifelse(fmode == "all", 1, 0))
+```
+
+```{r fit_model_2ei_slope}
+
+# Random intercept for neurons,
+# random intercept for experiments, nested in series, nested in mice
+lmer.2ei_slope = lmer(slope ~ allSpikes + (1 | uid) + (1 | mid/sid/eid), 
+                      data = tb %>% drop_na(slope, allSpikes))
+
+display(lmer.2ei_slope)
+anova(lmer.2ei_slope)
+```
+
+\newpage
+## (II) Predict threshold based on firing mode (all spikes vs non-burst spikes)
+
+```{r, fit_model_2ei_threshold}
+
+# Random intercept for neurons,
+# random intercept for experiments, nested in series, nested in mice
+lmer.2ei_thresh = lmer(thresh ~ allSpikes + (1 | uid) + (1 | mid/sid/eid), 
+                       data = tb %>% drop_na(thresh, allSpikes))
+
+display(lmer.2ei_thresh)
+anova(lmer.2ei_thresh)
+```
+

R/figure_3.Rmd

@@ -0,0 +1,172 @@
+---
+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)
+```

R/figure_3_S1.Rmd

@@ -0,0 +1,283 @@
+---
+title: "Spacek et al., 2021, Figure 3-Supplement 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}
+
+# This is the same file we used for Figure 1-Supplement 3. It contains the dependent variables for both, movies (Fig 1-Suppl 3) and gratings (Fig 3-Suppl 1)
+tib = get_data("../csv/fig1S33S1.csv")
+```
+
+```{r tidy_for_3_S1, include = FALSE}
+
+# Remove observations with NAs
+tb <- tib %>% drop_na(sbc)
+
+# Turn cell type (suppressed-by-contrast) into binary variable
+tb$sbc = ifelse(tb$sbc == TRUE, 1, 0)
+```
+
+# Figure 3-Supplement 1a
+## Firing rate FMIs separated by cell types
+
+```{r fit_model_3_S1a}
+
+# A mixed-effects model with any random intercept gives singular fits,
+# we therefore revert to fitting an ordinary linear model
+
+lmer.3_S1a = lm(grt_meanrate ~ sbc,
+                  data = tb %>% drop_na(grt_meanrate))
+
+display(lmer.3_S1a)
+anova(lmer.3_S1a)
+```
+
+```{r get_predicted_average_effect_3_S1a, include=F}
+
+m_sbc0 = coef(lmer.3_S1a)[1]
+diffMeans = coef(lmer.3_S1a)[2]
+m_sbc1  = coef(lmer.3_S1a)[1] + diffMeans
+```
+
+```{r save_coefficients_3_S1a, include=FALSE}
+
+pred_means_df = data.frame("non_sbc" = m_sbc0, "sbc" = m_sbc1, row.names = "")
+write_csv(pred_means_df, "_stats/figure_3_S1a_pred_means.csv")
+```
+
+FMI sbc: `r format(m_sbc1, digits=3, nsmall=3)` \newline
+FMI non-sbc: `r format(m_sbc0, digits=3, nsmall=3)` \newline
+n = `r nrow(tb %>% drop_na(grt_meanrate) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(grt_meanrate) %>% count(mid))` mice
+
+\newpage
+# Figure 3-Supplement 1b
+## Relation between firing rate FMI and recording depth
+
+```{r fit_model_3S1b}
+
+# Random intercept for mice
+lmer.3_S1b = lmer(grt_meanrate ~ depth  + (1 | mid), 
+                  data = tib %>% drop_na(grt_meanrate, depth))
+
+display(lmer.3_S1b)
+anova(lmer.3_S1b)
+```
+
+```{r save_coefficients_3_S1b, include=FALSE}
+
+coef_df = data.frame("intercept" = fixef(lmer.3_S1b)[1], "slope" = fixef(lmer.3_S1b)[2], row.names = "")
+write_csv(coef_df, "_stats/figure_3_S1b_coefs.csv")
+```
+
+Slope of `r format(fixef(lmer.3_S1b)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.3_S1b)[2], digits=2, nsmall=2)` (95\%-confidence interval) \newline
+n = `r nrow(tib %>% drop_na(grt_meanrate, depth) %>% count(uid))` neurons from `r nrow(tib %>% drop_na(grt_meanrate, depth) %>% count(mid))` mice
+
+\newpage
+# Figure 3-Supplement 1c
+## Relation between firing rate FMI and direction selectivity
+
+``` {r fit_model_3_S1c}
+
+# Random intercept for series, nested in mice
+lmer.3_S1c = lmer(grt_meanrate ~ dsi + (1 | sid/mid), 
+                  data = tib %>% drop_na(grt_meanrate, dsi))
+
+display(lmer.3_S1c)
+anova(lmer.3_S1c)
+```
+
+```{r save_coefficients_3_S1c, include=FALSE}
+
+coef_df = data.frame("intercept" = fixef(lmer.3_S1c)[1], "slope" = fixef(lmer.3_S1c)[2], row.names = "")
+write_csv(coef_df, "_stats/figure_3_S1c_coefs.csv")
+
+```
+
+Slope of `r format(fixef(lmer.3_S1c)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.3_S1c)[2], digits=2, nsmall=2)` (95\%-confidence interval) \newline
+n = `r nrow(tib %>% drop_na(grt_meanrate, dsi) %>% count(uid))` neurons from `r nrow(tib %>% drop_na(grt_meanrate, dsi) %>% count(mid))` mice
+
+\newpage
+# Figure 3-Supplement 1d
+## Relation between firing rate FMI and receptive field location
+
+``` {r fit_model_3_S1d}
+
+# Random intercept for series, nested in mice
+lmer.3_S1d = lmer(grt_meanrate ~ rfdist + (1 | mid/sid), 
+                  data = tib %>% drop_na(grt_meanrate, rfdist))
+
+display(lmer.3_S1d)
+anova(lmer.3_S1d)
+```
+
+```{r save_coefficients_3_S1d,include=FALSE}
+
+# store results in file
+coef_df = data.frame("intercept" = fixef(lmer.3_S1d)[1], "slope" = fixef(lmer.3_S1d)[2], row.names = "")
+write_csv(coef_df, "_stats/figure_3_S1d_coefs.csv")
+```
+
+Slope of `r format(fixef(lmer.3_S1d)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.3_S1d)[2], digits=2, nsmall=2)` (95\%-confidence interval) \newline
+n = `r nrow(tib %>% drop_na(grt_meanrate, rfdist) %>% count(uid))` neurons from `r nrow(tib %>% drop_na(grt_meanrate, rfdist) %>% count(mid))` mice
+
+\newpage
+# Figure 3-Supplement 1e
+## Relation between firing rate FMI and firing rate
+
+``` {r, fit_model_3_S1e}
+
+# Random intercept for series, nested in mice
+lmer.3_S1e = lmer(grt_meanrate ~ grt_meanrate_raw + (1 | mid/sid), 
+                  data = tib %>% drop_na(grt_meanrate, grt_meanrate_raw))
+
+display(lmer.3_S1e)
+anova(lmer.3_S1e)
+```
+
+```{r save_coefficients_3_S1e, include=F}
+
+coef_df = data.frame("intercept" = fixef(lmer.3_S1e)[1], "slope" = fixef(lmer.3_S1e)[2], row.names = "")
+write_csv(coef_df, "_stats/figure_3_S1e_coefs.csv")
+```
+
+Slope of `r format(fixef(lmer.3_S1e)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.3_S1e)[2], digits=2, nsmall=2)` (95\%-confidence interval) \newline
+n = `r nrow(tib %>% drop_na(grt_meanrate, grt_meanrate_raw) %>% count(uid))` neurons from `r nrow(tib %>% drop_na(grt_meanrate, grt_meanrate_raw) %>% count(mid))` mice
+
+\newpage
+# Figure 3-Supplement 1f
+## Burst ratio FMIs separated by cell types
+
+```{r fit_model_3_S1f}
+
+# Random intercept for series
+lmer.3_S1f = lmer(grt_meanburstratio ~ sbc  + (1 | sid), 
+                  data = tb %>% drop_na(grt_meanburstratio))
+
+display(lmer.3_S1f)
+anova(lmer.3_S1f)
+```
+
+```{r get_predicted_average_effect_3_S1f, include=F}
+
+m_sbc0 = fixef(lmer.3_S1f)[1]
+diffMeans = fixef(lmer.3_S1f)[2]
+m_sbc1  = fixef(lmer.3_S1f)[1] + diffMeans
+```
+
+```{r, save_coefficients_3_S1f, include=F}
+
+pred_means_df = data.frame("non_sbc" = m_sbc0, "sbc" = m_sbc1, row.names = "")
+write_csv(pred_means_df, "_stats/figure_3_S1f_pred_means.csv")
+```
+
+FMI sbc: `r format(m_sbc1, digits=3, nsmall=3)` \newline
+FMI non-sbc: `r format(m_sbc0, digits=3, nsmall=3)` \newline
+n = `r nrow(tb %>% drop_na(grt_meanburstratio) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(grt_meanburstratio) %>% count(mid))` mice
+
+\newpage
+# Figure 3-Supplement 1g
+## Relation between burst ratio FMI and recording depth
+
+
+```{r fit_model_3_S1g}
+
+# Random intercept for series, nested in mice
+lmer.3_S1g = lmer(grt_meanburstratio ~ depth  + (1 | mid/sid), 
+                  data = tib %>% drop_na(grt_meanburstratio, depth))
+
+display(lmer.3_S1g)
+anova(lmer.3_S1g)
+```
+
+```{r save_coefficients_3S1g, include=FALSE}
+
+# store results in file
+coef_df = data.frame("intercept" = fixef(lmer.3_S1g)[1], "slope" = fixef(lmer.3_S1g)[2], row.names = "")
+write_csv(coef_df, "_stats/figure_3_S1g_coefs.csv")
+```
+
+Slope of `r format(fixef(lmer.3_S1g)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.3_S1g)[2], digits=2, nsmall=2)` (95\%-confidence interval) \newline
+n = `r nrow(tib %>% drop_na(grt_meanburstratio, depth) %>% count(uid))` neurons from `r nrow(tib %>% drop_na(grt_meanburstratio, depth) %>% count(mid))` mice
+
+\newpage
+# Figure 3-Supplement 1h
+## Relation between burst ratio FMI and direction selectivity
+
+``` {r, fit_model_3_S1h}
+
+# Random intercept for series, nested in mice
+lmer.3_S1h = lmer(grt_meanburstratio ~ dsi + (1 | mid/sid), 
+                  data = tib %>% drop_na(grt_meanburstratio, dsi))
+
+display(lmer.3_S1h)
+anova(lmer.3_S1h)
+```
+
+```{r, save_coefficients_3_S1h, include=FALSE}
+
+coef_df = data.frame("intercept" = fixef(lmer.3_S1h)[1], "slope" = fixef(lmer.3_S1h)[2], row.names = "")
+write_csv(coef_df, "_stats/figure_3_S1h_coefs.csv")
+```
+
+Slope of `r format(fixef(lmer.3_S1h)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.3_S1h)[2], digits=2, nsmall=2)` (95\%-confidence interval) \newline
+n = `r nrow(tib %>% drop_na(grt_meanburstratio, dsi) %>% count(uid))` neurons from `r nrow(tib %>% drop_na(grt_meanburstratio, dsi) %>% count(mid))` mice
+
+\newpage
+# Figure 3-Supplement 1i
+## Relation between burst ratio FMI and receptive field location
+
+``` {r, fit_model_3_S1i}
+
+# Random intercept for series, nested in mice
+lmer.3_S1i = lmer(grt_meanburstratio ~ rfdist + (1 | mid/sid), 
+                  data = tib %>% drop_na(grt_meanburstratio, rfdist))
+
+display(lmer.3_S1i)
+anova(lmer.3_S1i)
+```
+
+```{r, save_coefficients_3_S1i, include=FALSE}
+
+coef_df = data.frame("intercept" = fixef(lmer.3_S1i)[1], "slope" = fixef(lmer.3_S1i)[2], row.names = "")
+write_csv(coef_df, "_stats/figure_3_S1i_coefs.csv")
+```
+
+Slope of `r format(fixef(lmer.3_S1i)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.3_S1i)[2], digits=2, nsmall=2)` (95\%-confidence interval) \newline
+n = `r nrow(tib %>% drop_na(grt_meanburstratio, rfdist) %>% count(uid))` neurons from `r nrow(tib %>% drop_na(grt_meanburstratio, rfdist) %>% count(mid))` mice
+
+
+\newpage
+# Figure 3-Supplement 1j
+## Relation between burst ratio FMI and burst ratio
+
+``` {r, fit_model_3_S1j}
+
+# Random intercept for series, nested in mice
+lmer.3_S1j = lmer(grt_meanburstratio ~ grt_meanburstratio_raw + (1 | mid/sid), 
+                  data = tib %>% drop_na(grt_meanburstratio, grt_meanburstratio_raw))
+
+display(lmer.3_S1j)
+anova(lmer.3_S1j)
+```
+
+```{r, save_coefficients, include=FALSE}
+
+coef_df = data.frame("intercept" = fixef(lmer.3_S1j)[1], "slope" = fixef(lmer.3_S1j)[2], row.names = "")
+write_csv(coef_df, "_stats/figure_3S1j_coefs.csv")
+```
+
+Slope of `r format(fixef(lmer.3_S1j)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.3_S1j)[2], digits=2, nsmall=2)` (95\%-confidence interval) \newline
+n = `r nrow(tib %>% drop_na(grt_meanburstratio, grt_meanburstratio_raw) %>% count(uid))` neurons from `r nrow(tib %>% drop_na(grt_meanburstratio, grt_meanburstratio_raw) %>% count(mid))` mice
+

R/figure_4.Rmd

@@ -0,0 +1,161 @@
+---
+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
+
+

R/figure_4_S1.Rmd

@@ -0,0 +1,461 @@
+---
+title: "Spacek et al., 2021, Figure 4-Supplement 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_4_S1a, include=FALSE}
+
+tib = get_data("../csv/grtFMI.csv")
+```
+
+```{r tidy_for_4_S1, include = FALSE}
+
+tb <- tib %>% filter(st8 == "none") %>% select(mid, sid, eid, uid, mseu, meanrate, meanburstratio)
+```
+
+# Figure 4-Supplement 1a
+## Relation between firing rate FMI and burst ratio FMI
+
+```{r fit_model_4_S1a}
+
+# Random intercept for neurons,
+# random intercept for experiments, nested in series, nested in mice
+lmer.4_S1a = lmer(meanburstratio ~ meanrate + (1 | uid) + (1 | mid/sid/eid), 
+                  data = tb %>% drop_na(meanburstratio, meanrate))
+
+display(lmer.4_S1a)
+anova(lmer.4_S1a)
+```
+
+```{r, save_coefficients_4_S1a, include=FALSE}
+
+coef_df = data.frame("intercept" = fixef(lmer.4_S1a)[1], "slope" = fixef(lmer.4_S1a)[2], row.names = "")
+write_csv(coef_df, "_stats/figure_4_S1a_coefs.csv")
+```
+
+Slope of `r format(fixef(lmer.4_S1a)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.4_S1a)[2], digits=2, nsmall=2)` (95\%-confidence interval) \newline
+n = `r nrow(tb %>% drop_na(meanburstratio, meanrate) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(meanburstratio, meanrate) %>% count(mid))` mice
+
+\newpage
+# Figure 4-Supplement 1b
+## Relation between firing rates in response to movies vs gratings
+
+
+```{r read_data_4_S1b, include=FALSE}
+
+tib = get_data("../csv/fig4S1b.csv")
+```
+
+```{r tidy_for_4_S1b, include = FALSE}
+
+# Turn stimulus type into binary predictor
+tb <- tib %>% mutate(mvi = ifelse(stimtype == "mvi", 1, 0)) %>% select(mid, sid, uid, msu, mvi, meanrate)
+```
+
+```{r fit_model_4_S1b}
+
+# Random intercept for neurons, 
+# random intercept for series, nested in mice
+lmer.4_S1b = lmer(meanrate ~ mvi + (1 | uid) + (1 | mid/sid), 
+                  data = tb %>% drop_na(meanrate))
+
+display(lmer.4_S1b)
+anova(lmer.4_S1b)
+```
+
+```{r get_predicted_average_effect_4_S1b, include=F}
+
+m_grt = fixef(lmer.4_S1b)[1]
+diffMeans = fixef(lmer.4_S1b)[2]
+m_mvi  = fixef(lmer.4_S1b)[1] + diffMeans
+```
+
+Movies: `r format(m_mvi, digits=2, nsmall=2)` spikes/s \newline
+Gratings: `r format(m_grt, digits=2, nsmall=2)` spikes/s \newline
+n = `r nrow(tb %>% drop_na(meanrate) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(meanrate) %>% count(mid))` mice
+
+\newpage
+# Figure 4-Supplement 1c,d
+## Relation between firing rates during feedback vs suppression for separate epochs of the movie
+
+```{r read_data_4_S1cd, include=FALSE}
+
+tib = get_data("../csv/fig1.csv")
+```
+
+```{r tidy_4_S1cd, include=FALSE}
+
+# Get relevant variables, rename rate* columns, such that they have a common prefix, followed by a simple number
+tmp1 = tib %>% select(mid, sid, eid, uid, mseu, trialis, opto, rate02s, rate35s) %>% drop_na(trialis, rate02s, rate35s) %>% rename(rate2 = rate02s, rate35 = rate35s)
+
+# Turn rate* columns into long format
+tmp2 = tmp1 %>% pivot_longer(cols = starts_with("rate"), names_to = "window", names_prefix = "rate", values_to = "trialrate")
+
+# Turn booleans for 'optogenetic manipulation' and 'time window' into binary, numeric predictors
+tmp2 = tmp2 %>% mutate(feedback = ifelse(opto == TRUE, 0, 1)) %>% mutate(earlywin = ifelse(window == '2', 1, 0))
+```
+
+```{r fit_model_4_S1cd}
+
+lmer.4_S1cd = lmer(trialrate ~ feedback*earlywin + (1 + feedback*earlywin | uid) + (1 | mid/sid/eid),
+                   data = tmp2 %>% drop_na(trialrate))
+
+display(lmer.4_S1cd)
+anova(lmer.4_S1cd)
+```
+
+n = `r nrow(tmp2 %>% drop_na(trialrate) %>% count(uid))` neurons from `r nrow(tmp2 %>% drop_na(trialrate) %>% count(mid))` mice
+
+```{r post_hocs_4_S1cd, include=FALSE}
+
+# Post-hoc: fit model to early window (0-2s)
+foo = tmp2 %>% filter(earlywin == 1)
+
+lmer.4_S1cd = lmer(trialrate ~ feedback + (1 + feedback | uid) + (1 | sid/eid),
+                   data = foo %>% drop_na(trialrate))
+
+display(lmer.4_S1cd)
+anova(lmer.4_S1cd)
+
+m_suppr = fixef(lmer.4_S1cd)[1]
+diffMeans = fixef(lmer.4_S1cd)[2]
+m_feedb  = fixef(lmer.4_S1cd)[1] + diffMeans
+
+# Post-hoc: fit model to late window (3-5s)
+foo = tmp2 %>% filter(earlywin == 0)
+
+lmer.4_S1cd = lmer(trialrate ~ feedback + (1 + feedback | uid) + (1 | sid/eid),
+                   data = foo %>% drop_na(trialrate))
+
+display(lmer.4_S1cd)
+anova(lmer.4_S1cd)
+
+m_suppr = fixef(lmer.4_S1cd)[1]
+diffMeans = fixef(lmer.4_S1cd)[2]
+m_feedb  = fixef(lmer.4_S1cd)[1] + diffMeans
+```
+
+\newpage
+# Figure 4-Supplement 1e
+## Burst ratio FMIs separately for grating versus movie presentations
+
+```{r read_data_4_S1ef, include=FALSE}
+
+tib = get_data("../csv/fig4.csv")
+```
+
+```{r tidy_for_4_S1e, include=FALSE}
+
+t4S1e <- tib %>% filter(blank == "False" & meanburstratio != "NA")
+
+# Turn stimulus type into binary predictor
+t4S1e = t4S1e %>% mutate(mvi = ifelse(stimtype == "mvi", 1, 0))
+```
+
+```{r fit_model_4_S1e}
+
+# Random intercept for neurons,
+# random intercept for mice
+lmer.4_S1e = lmer(meanburstratio ~ mvi + (1 | uid) + (1 | mid), 
+                  data = t4S1e %>% drop_na(meanburstratio))
+
+display(lmer.4_S1e)
+anova(lmer.4_S1e)
+```
+
+```{r predicted_average_effect_4_S1e, include=F}
+
+m_grt = fixef(lmer.4_S1e)[1]
+diffMeans = fixef(lmer.4_S1e)[2]
+m_mvi  = fixef(lmer.4_S1e)[1] + diffMeans
+```
+
+```{r save_coefficients_4_S1e, include=FALSE}
+
+pred_means_df = data.frame("mvi" = m_mvi, "grt" = m_grt, row.names = "")
+write_csv(pred_means_df, "_stats/figure_4_S1e_pred_means.csv")
+```
+
+FMI movies: `r format(m_mvi, digits=2, nsmall=2)` \newline
+FMI gratings: `r format(m_grt, digits=2, nsmall=2)` \newline
+n = `r nrow(t4S1e %>% drop_na(meanburstratio) %>% count(uid))` neurons from `r nrow(t4S1e %>% drop_na(meanburstratio) %>% count(mid))` mice 
+
+\newpage
+# Figure 4-Supplement 1f
+## Burst ratio FMIs across blank screen conditions
+
+```{r tidy_for_4_S1f, include=FALSE}
+
+# Code the 'blank condition' as integers 1, 2, 3
+t3f = tib %>% mutate(blank_condition = case_when(
+  stimtype == 'mvi' & blank == 'prestim' ~ "mvi",
+  stimtype == 'grt' & blank == 'prestim' ~ "grt",
+  stimtype == 'grt' & blank == 'cond' ~ "grt0c"
+))
+
+# Remove NAs
+t3f = t3f %>% filter(meanburstratio != "Na" & blank_condition != "Na")
+
+# And turn them into a factor for dummy coding
+t3f$blank_condition = as.factor(t3f$blank_condition)
+```
+
+
+```{r fit_model_4_S1f}
+
+# Random intercept for neurons,
+# random intercept for series
+lmer.4_S1f = lmer(meanburstratio ~ blank_condition + (1 | uid) + (1 | sid), 
+                  data = t3f %>% drop_na(meanburstratio))
+
+display(lmer.4_S1f)
+anova(lmer.4_S1f)
+```
+
+```{r get_predicted_average_effect_4_S1f, include=FALSE}
+
+m_grt = fixef(lmer.4_S1f)[1]
+m_grt0c = fixef(lmer.4_S1f)[1] + fixef(lmer.4_S1f)[2]
+m_mvi = fixef(lmer.4_S1f)[1]+ fixef(lmer.4_S1f)[3]
+```
+
+```{r save_coefficients_4_S1f, include=FALSE}
+
+pred_means_df = data.frame("grt" = m_grt, "grt0c" = m_grt0c, "mvi" = m_mvi, row.names = "")
+write_csv(pred_means_df, "_stats/figure_4_S1f_pred_means.csv")
+```
+
+FMI pre-movie:`r format(m_mvi, digits=2, nsmall=2)` \newline
+FMI pre-grating: `r format(m_grt, digits=2, nsmall=2)` \newline
+FMI blank grating: `r format(m_grt0c, digits=2, nsmall=2)` \newline
+n = `r nrow(t3f %>% drop_na(meanburstratio) %>% count(uid))` neurons from `r nrow(t3f %>% drop_na(meanburstratio) %>% count(mid))` mice 
+
+\newpage
+## Post-hoc analysis, comparing mean across blank conditions against stimulus condition
+
+```{r tidy_for_post_hoc_4_S1ef, include = F}
+
+t4S1ef = tib %>% filter(meanburstratio != "Na")
+
+t4S1ef = t4S1ef %>% mutate(stim_condition = case_when(
+  stimtype == 'mvi' & blank == 'False' ~ 1,
+  stimtype == 'grt' & blank == 'False' ~ 2,
+  (stimtype == 'mvi' | stimtype == 'grt') & blank != 'False' ~ 3
+))
+
+t4S1ef$stim_condition = as.factor(t4S1ef$stim_condition)
+```
+
+```{r fit_model_4_S1ef, echo=F}
+
+# Random intercept for neuron,
+# random intercept for series
+lmer.4_S1ef = lmer(meanburstratio ~ stim_condition + (1 | uid) + (1 | sid), 
+                   data = t4S1ef %>% drop_na(meanburstratio))
+
+display(lmer.4_S1ef)
+anova(lmer.4_S1ef)
+```
+
+```{r predicted_average_effect_4_S1ef, include=FALSE}
+
+m_mvi = fixef(lmer.4_S1ef)[1]
+m_grt = fixef(lmer.4_S1ef)[1] + fixef(lmer.4_S1ef)[2]
+m_blank = fixef(lmer.4_S1ef)[1]+ fixef(lmer.4_S1ef)[3]
+```
+
+
+\newpage
+# Figure 4-Supplement 1g
+## Feedback effect on firing rate during blank periods preceding movies
+
+
+```{r read_data_4_S1gj, include=FALSE}
+
+tib = get_data("../csv/fig1.csv")
+```
+
+```{r tidy_for_4_S1gj, include = FALSE}
+
+# Turn stimulus type into binary predictor, get relevant columns
+tb <- tib %>% mutate(feedback = ifelse(opto == "TRUE", 0, 1)) %>% select(mid, sid, eid, uid, mseu, feedback, everything(), -opto)
+```
+
+```{r fit_model_4_S1g}
+
+# Random intercept, random slope for neurons,
+# random intercept for experiments, nested in series, nested in mice
+lmer.4_S1g = lmer(blankrates ~ feedback + (1 + feedback | uid) + (1 | mid/sid/eid), 
+                  data = tb %>% drop_na(blankrates))
+
+display(lmer.4_S1g)
+anova(lmer.4_S1g)
+```
+
+```{r predicted_average_effect_4_S1g, include=F}
+
+m_suppress = fixef(lmer.4_S1g)[1]
+diffMeans = fixef(lmer.4_S1g)[2]
+m_feedback  = fixef(lmer.4_S1g)[1] + diffMeans
+```
+
+Feedback: `r format(m_feedback, digits=2, nsmall=2)` spikes/s. \newline
+Suppression: `r format(m_suppress, digits=2, nsmall=2)` spikes/s \newline
+n = `r nrow(tb %>% drop_na(blankrates) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(blankrates) %>% count(mid))` mice 
+
+\newpage
+# Figure 4-Supplement 1j
+## Feedback effect on burst ratio during blank periods preceding movies
+
+```{r fit_model_4_S1j}
+
+# Random intercept, random slope for neurons,
+# random intercept for experiments
+lmer.4_S1j = lmer(blankburstratios ~ feedback + (1 + feedback | uid) + (1 | eid), 
+                  data = tb %>% drop_na(blankburstratios))
+
+display(lmer.4_S1j)
+anova(lmer.4_S1j)
+```
+
+```{r predicted_average_effect_4_S1j, include=F}
+
+m_suppress = fixef(lmer.4_S1j)[1]
+diffMeans = fixef(lmer.4_S1j)[2]
+m_feedback  = fixef(lmer.4_S1j)[1] + diffMeans
+```
+
+Feedback: burst ratio of `r format(m_feedback, digits=2, nsmall=2)` \newline
+Suppression: burst ratio of `r format(m_suppress, digits=2, nsmall=2)` \newline
+n = `r nrow(tb %>% drop_na(blankburstratios) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(blankburstratios) %>% count(mid))` mice 
+
+\newpage
+
+```{r read_data_4_S1hkil, include=FALSE}
+
+tib = get_data("../csv/fig3.csv")
+```
+
+```{r tidy_for_$_S1hkil, include = FALSE}
+
+# Turn stimulus type into binary predictor, get relevant columns
+tb <- tib %>% mutate(feedback = ifelse(opto == "TRUE", 0, 1)) %>% select(mid, sid, eid, uid, mseu, feedback, everything(), -opto)
+```
+
+\newpage
+# Figure 4-Supplement 1h
+## Feedback effect on firing rate during blank periods preceding gratings
+
+```{r fit_model_4_S1h}
+
+# Random intercept for neurons,
+# random intercept for experiments, nested in series, nested in mice
+lmer.4_S1h = lmer(blankrates ~ feedback + (1 + feedback | uid) + (1 | mid/sid/eid), 
+                  data = tb %>% drop_na(blankrates))
+
+display(lmer.4_S1h)
+anova(lmer.4_S1h)
+```
+
+```{r predicted_average_effect_4_S1h, include=F}
+
+m_suppress = fixef(lmer.4_S1h)[1]
+diffMeans = fixef(lmer.4_S1h)[2]
+m_feedback  = fixef(lmer.4_S1h)[1] + diffMeans
+```
+
+Feedback: `r format(m_feedback, digits=2, nsmall=2)` spikes/s \newline
+Suppression: `r format(m_suppress, digits=2, nsmall=2)` spikes/s \newline
+n = `r nrow(tb %>% drop_na(blankrates) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(blankrates) %>% count(mid))` mice 
+
+\newpage
+# Figure 4-Supplement 1k
+## Feedback effect on burst ratio during blank periods preceding gratings
+
+```{r fit_model_4_S1k}
+
+# Random intercept, random slope for neurons,
+# random intercept for series, nested in mice
+lmer.4_S1k = lmer(blankburstratios ~ feedback + (1 + feedback | uid) + (1 | mid/sid), 
+                  data = tb %>% drop_na(blankburstratios))
+
+display(lmer.4_S1k)
+anova(lmer.4_S1k)
+```
+
+```{r predicted_average_effect_4_S1k, include=F}
+
+m_suppress = fixef(lmer.4_S1k)[1]
+diffMeans = fixef(lmer.4_S1k)[2]
+m_feedback  = fixef(lmer.4_S1k)[1] + diffMeans
+```
+
+Feedback: burst ratio of `r format(m_feedback, digits=2, nsmall=2)` \newline
+Suppression: burst ratio of `r format(m_suppress, digits=2, nsmall=2)` \newline
+n = `r nrow(tb %>% drop_na(blankburstratios) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(blankburstratios) %>% count(mid))` mice 
+
+\newpage
+# Figure 4-Supplement 1i
+## Feedback effect on firing rate during zero-contrast gratings
+
+
+```{r fit_model_4_S1i}
+
+# Random intercept, random slope for neurons,
+# random intercept for experiments, nested in mice
+lmer.4_S1i = lmer(blankcondrates ~ feedback + (1 + feedback | uid) + (1 | sid/eid), 
+                  data = tb %>% drop_na(blankcondrates))
+
+display(lmer.4_S1i)
+anova(lmer.4_S1i)
+```
+
+```{r predicted_average_effect_4_S1i, include=F}
+
+m_suppress = fixef(lmer.4_S1i)[1]
+diffMeans = fixef(lmer.4_S1i)[2]
+m_feedback  = fixef(lmer.4_S1i)[1] + diffMeans
+```
+
+Feedback: `r format(m_feedback, digits=2, nsmall=2)` spikes/s \newline
+Suppression: `r format(m_suppress, digits=2, nsmall=2)` spikes/s \newline
+n = `r nrow(tb %>% drop_na(blankcondrates) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(blankcondrates) %>% count(mid))` mice 
+
+\newpage
+# Figure 4-Supplement 1l
+## Feedback effect on burst ratio during zero-contrast gratings
+
+```{r fit_model_4_S1l}
+
+# Random intercept for neurons,
+# random intercept for series, nested in mice
+lmer.4_S1l = lmer(blankcondburstratios ~ feedback + (1 | uid) + (1 | mid/sid), 
+                  data = tb %>% drop_na(blankcondburstratios))
+
+display(lmer.4_S1l)
+anova(lmer.4_S1l)
+```
+
+```{r predicted_average_effect_4_S1l, include=F}
+
+m_suppress = fixef(lmer.4_S1l)[1]
+diffMeans = fixef(lmer.4_S1l)[2]
+m_feedback  = fixef(lmer.4_S1l)[1] + diffMeans
+```
+
+Feedback: burst ratio of `r format(m_feedback, digits=2, nsmall=2)` \newline
+Suppression: burst ratio of `r format(m_suppress, digits=2, nsmall=2)` \newline
+n = `r nrow(tb %>% drop_na(blankcondburstratios) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(blankcondburstratios) %>% count(mid))` mice 
+
+

R/figure_5.Rmd

@@ -0,0 +1,255 @@
+---
+title: "Spacek et al., 2021, Figure 5"
+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/fig5.csv")
+```
+
+```{r tidy, include = FALSE}
+
+# Transform columns st8 and opto into binary, numerical predictors
+tb = tib %>% mutate(run = ifelse(st8 == "run", 1, 0))
+tb = tb  %>% mutate(feedback = ifelse(opto == "TRUE", 0, 1))
+
+tbf = tb %>% filter(feedback == 1) # feedback, panels c-f
+tbs = tb %>% filter(feedback == 0) # suppression, panels i-l
+```
+
+# Figure 5c
+## Effect of locomotion state on firing rate
+
+```{r fit_model_5c}
+
+# Random intercept, random slope for neurons,
+# random intercept for experiments, nested in series
+lmer.5c = lmer(rates ~ run + (1 + run | uid) + (1 | sid/eid), 
+              data = tbf %>% drop_na(rates))
+
+display(lmer.5c)
+anova(lmer.5c)
+```
+
+```{r get_predicted_average_effect_5c, include=F}
+
+mSit = fixef(lmer.5c)[1]
+diffMeans = fixef(lmer.5c)[2]
+mRun  = fixef(lmer.5c)[1] + diffMeans
+```
+
+Running: `r format(mRun, digits=1, nsmall=1)` spikes/s \newline
+Sitting: `r format(mSit, digits=1, nsmall=1)` spikes/s \newline
+n = `r nrow(tbf %>% drop_na(rates) %>% count(uid))` neurons from `r nrow(tbf %>% drop_na(rates) %>% count(mid))` mice
+
+\newpage
+# Figure 5d
+## Effect of locomotion state on burst ratio
+
+```{r fit_model_5d}
+
+# Random intercept, random slope for neurons,
+# random intercept for experiments, nested in series
+lmer.5d = lmer(burstratios ~ run + (1 + run | uid) + (1 | sid/eid), 
+               data = tbf %>% drop_na(burstratios))
+
+display(lmer.5d)
+anova(lmer.5d)
+```
+
+```{r get_predicted_average_effect_5d, include=F}
+
+mSit = fixef(lmer.5d)[1]
+diffMeans = fixef(lmer.5d)[2]
+mRun  = fixef(lmer.5d)[1] + diffMeans
+```
+
+Running: `r format(mRun, digits=2, nsmall=2)` \newline
+Sitting: `r format(mSit, digits=2, nsmall=2)` \newline
+n = `r nrow(tbf %>% drop_na(burstratios) %>% count(uid))` neurons from `r nrow(tbf %>% drop_na(burstratios) %>% count(mid))` mice
+
+\newpage
+# Figure 5e
+## Effect of locomotion state on sparseness
+
+```{r tidy_for_5ef, 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
+
+tbfef = tbf %>% select(mid, sid, eid, uid, mseu, run, spars, rel) %>% distinct(mseu, run, .keep_all = TRUE)
+
+```
+
+```{r fit_model_5e}
+
+# Random intercept for neurons,
+# random intercept for experiments, nested in series
+lmer.5e = lmer(spars ~ run + (1 | uid) + (1 | sid/eid), 
+               data = tbfef %>% drop_na(spars))
+
+display(lmer.5e)
+anova(lmer.5e)
+```
+
+```{r get_predicted_average_effect_5e, include=F}
+
+mSit = fixef(lmer.5e)[1]
+diffMeans = fixef(lmer.5e)[2]
+mRun  = fixef(lmer.5e)[1] + diffMeans
+```
+
+Running: `r format(mRun, digits=2, nsmall=2)` \newline
+Sitting: `r format(mSit, digits=2, nsmall=2)` \newline
+n = `r nrow(tbfef %>% drop_na(spars) %>% count(uid))` neurons from `r nrow(tbfef %>% drop_na(spars) %>% count(mid))` mice
+
+\newpage
+# Figure 5f
+## Effect of locomotion state on reliability
+
+```{r fit_model_5f}
+
+# Random intercept for neurons,
+# random intercept for experiments, nested in series
+lmer.5f = lmer(rel ~ run + (1 | uid) + (1 | sid/eid), 
+               data = tbfef %>% drop_na(rel))
+
+display(lmer.5f)
+anova(lmer.5f)
+```
+
+```{r get_predicted_average_effect_5f, include=F}
+
+mSit = fixef(lmer.5f)[1]
+diffMeans = fixef(lmer.5f)[2]
+mRun  = fixef(lmer.5f)[1] + diffMeans
+```
+
+Running: `r format(mRun, digits=2, nsmall=2)` \newline
+Sitting: `r format(mSit, digits=2, nsmall=2)` \newline
+n = `r nrow(tbfef %>% drop_na(rel) %>% count(uid))` neurons from `r nrow(tbfef %>% drop_na(rel) %>% count(mid))` mice
+
+\newpage
+# Figure 5i
+## Effect of locomotion state on firing rate during V1 suppression
+
+```{r fit_model_5i}
+
+# Random intercept, random slope for neurons,
+# random intercept for experiments nested in series
+lmer.5i = lmer(rates ~ run + (1 + run | uid) + (1 | sid/eid), 
+               data = tbs %>% drop_na(rates))
+
+display(lmer.5i)
+anova(lmer.5i)
+```
+
+```{r get_predicted_average_effect_5i, include=F}
+
+mSit = fixef(lmer.5i)[1]
+diffMeans = fixef(lmer.5i)[2]
+mRun  = fixef(lmer.5i)[1] + diffMeans
+```
+
+Running: `r format(mRun, digits=1, nsmall=1)` spikes/s \newline
+Sitting: `r format(mSit, digits=1, nsmall=1)` spikes/s \newline
+n = `r nrow(tbs %>% drop_na(rates) %>% count(uid))` neurons from `r nrow(tbs %>% drop_na(rates) %>% count(mid))` mice
+
+\newpage
+# Figure 5j
+## Effect of locomotion state on burst ratio during V1 suppression
+
+```{r fit_model_5j}
+
+# Random intercept, random slope for neurons,
+# random intercept for experiments, nested in series
+lmer.5j = lmer(burstratios ~ run + (1 + run | uid) + (1 | sid/eid), 
+               data = tbs %>% drop_na(burstratios))
+
+display(lmer.5j)
+anova(lmer.5j)
+```
+
+```{r get_predicted_average_effect_5j, include=F}
+
+mSit = fixef(lmer.5j)[1]
+diffMeans = fixef(lmer.5j)[2]
+mRun  = fixef(lmer.5j)[1] + diffMeans
+```
+
+Running: `r format(mRun, digits=2, nsmall=2)` \newline
+Sitting: `r format(mSit, digits=2, nsmall=2)` \newline
+n = `r nrow(tbs %>% drop_na(burstratios) %>% count(uid))` neurons from `r nrow(tbs %>% drop_na(burstratios) %>% count(mid))` mice
+
+```{r tidy_for_5kl, 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
+
+tbskl = tbs %>% select(mid, sid, eid, uid, mseu, run, spars, rel) %>% distinct(mseu, run, .keep_all = TRUE)
+```
+
+\newpage
+# Figure 5k
+## Effect of locomotion state on sparseness during V1 suppression
+
+```{r fit_model_5k}
+
+# Random intercept for neurons,
+# random intercept for experiments, nested in series, nested in mice
+lmer.5k = lmer(spars ~ run + (1 | uid) + (1 | mid/sid/eid), 
+               data = tbskl %>% drop_na(spars))
+
+display(lmer.5k)
+anova(lmer.5k)
+```
+
+```{r get_predicted_average_change_5k, include=F}
+
+mSit = fixef(lmer.5k)[1]
+diffMeans = fixef(lmer.5k)[2]
+mRun  = fixef(lmer.5k)[1] + diffMeans
+```
+
+Running: `r format(mRun, digits=2, nsmall=2)` \newline
+Sitting: `r format(mSit, digits=2, nsmall=2)` \newline
+n = `r nrow(tbskl %>% drop_na(spars) %>% count(uid))` neurons from `r nrow(tbskl %>% drop_na(spars) %>% count(mid))` mice
+
+\newpage
+# Figure 5l
+## Effect of locomotion state on reliability during V1 suppression
+
+```{r fit_model_5l}
+
+# Random intercept for neurons,
+# random intercept for experiments, nested in series
+lmer.5l = lmer(rel ~ run + (1 | uid) + (1 | sid/eid), 
+               data = tbskl %>% drop_na(rel))
+
+display(lmer.5l)
+anova(lmer.5l)
+```
+
+```{r get_predicted_average_effect_5l, include=F}
+
+mSit = fixef(lmer.5l)[1]
+diffMeans = fixef(lmer.5l)[2]
+mRun  = fixef(lmer.5l)[1] + diffMeans
+```
+
+Running: `r format(mRun, digits=2, nsmall=2)` \newline
+Sitting: `r format(mSit, digits=2, nsmall=2)` \newline
+n = `r nrow(tbskl %>% drop_na(rel) %>% count(uid))` neurons from `r nrow(tbskl %>% drop_na(rel) %>% count(mid))` mice
+

R/figure_5_S1.Rmd

@@ -0,0 +1,297 @@
+---
+title: "Spacek et al., 2021, Figure 5 Supplement-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_5_S1ab, include=FALSE}
+
+tib = get_data("../csv/fig5.csv")
+```
+
+```{r tidy_for_5_S1ab, include = FALSE}
+
+# Get relevant conditions, turn locomotion state into a binary predictor
+tb = tib %>% filter(opto == "FALSE") %>% filter(st8 == "run" | st8 == "sit") %>% mutate(run = ifelse(st8 == "run", 1, 0))
+
+# 'Signal-to-noise ratio', and 'mean peak width' 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
+
+tbab = tb %>% select(mid, sid, eid, uid, mseu, run, snr, meanpkw) %>% distinct(mseu, run, .keep_all = TRUE)
+```
+
+# Figure 5-Supplement 1a
+## Effect of locomotion state on signal-to-noise ratio (SNR)
+
+```{r fit_model_5_S1a}
+
+# Random intercept for neurons,
+# random intercept for experiments, nested in series
+lmer.5_S1a = lmer(snr ~ run + (1 | uid) + (1 | sid/eid), 
+                  data = tbab %>% drop_na(snr))
+
+display(lmer.5_S1a)
+anova(lmer.5_S1a)
+```
+
+```{r predicted_average_effect_5_S1a, include=F}
+
+mSit = fixef(lmer.5_S1a)[1]
+diffMeans = fixef(lmer.5_S1a)[2]
+mRun  = fixef(lmer.5_S1a)[1] + diffMeans
+```
+
+SNR locomotion: `r format(mRun, digits=2, nsmall=2)` \newline
+SNR sitting: `r format(mSit, digits=2, nsmall=2)` \newline
+n = `r nrow(tb %>% drop_na(snr) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(snr) %>% count(mid))` mice 
+
+\newpage
+# Figure 5-Supplement 1b
+## Effect of locomotion state on mean peak width
+
+```{r fit_model_5_S1b}
+
+# Random for neurons,
+# random intercept for series
+lmer.5_S1b = lmer(meanpkw ~ run + (1 | uid) + (1 | sid), 
+                  data = tbab %>% drop_na(meanpkw))
+
+display(lmer.5_S1b)
+anova(lmer.5_S1b)
+```
+
+```{r predicted_average_effect_5_S1b, include=F}
+
+mSit = fixef(lmer.5_S1b)[1]
+diffMeans = fixef(lmer.5_S1b)[2]
+mRun  = fixef(lmer.5_S1b)[1] + diffMeans
+```
+
+Mean peak width running: `r format(mRun, digits=2, nsmall=2)` \newline
+Mean peak width sitting: `r format(mSit, digits=2, nsmall=2)` \newline
+n = `r nrow(tb %>% drop_na(meanpkw) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(meanpkw) %>% count(mid))` mice 
+
+\newpage
+# Figure 5-Supplement 1c
+## Relation between firing rate RMI and burst ratio RMI
+
+```{r read_data_5_S1cdefg, include=FALSE}
+
+tib = get_data("../csv/mviRMI.csv")
+```
+
+```{r tidy_for_5_S1cdefg, include=FALSE}
+
+tb <- tib %>% filter(opto == 'FALSE')
+```
+
+```{r fit_model_5_S1c}
+
+# Remove outliers
+tb_clean <- tb %>% filter(meanburstratio < 0.99 & meanburstratio > -0.99)
+
+# Random intercept for neurons,
+# random intercept for experiments, nested in series
+lmer.5_S1_c = lmer(meanburstratio ~ meanrate + (1 | uid) + (1 | sid/eid), 
+                   data = tb_clean %>% drop_na(meanburstratio, meanrate))
+
+display(lmer.5_S1_c)
+anova(lmer.5_S1_c)
+```
+
+```{r save_coefficients_5_S1c, include=FALSE}
+
+coef_df = data.frame("intercept" = fixef(lmer.5_S1_c)[1], "slope" = fixef(lmer.5_S1_c)[2], row.names = "")
+write_csv(coef_df, "_stats/figure_5_S1c_coefs.csv")
+```
+
+Slope of `r format(fixef(lmer.5_S1_c)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.5_S1_c)[2], digits=2, nsmall=2)` (95-\% confidence interval) \newline
+n = `r nrow(tb_clean %>% drop_na(meanburstratio, meanrate) %>% count(uid))` neurons from `r nrow(tb_clean %>% drop_na(meanburstratio, meanrate) %>% count(mid))` mice 
+
+\newpage
+# Figure 5-Supplement 1d
+## Relation between firing rate RMI and sparseness RMI
+
+```{r fit_model_5_S1d}
+
+# Random intercept for neurons,
+# random intercept for experiments, nested in series
+lmer.5_S1_d = lmer(spars ~ meanrate + (1 | uid) + (1 | sid/eid), 
+                   data = tb %>% drop_na(spars, meanrate))
+
+display(lmer.5_S1_d)
+anova(lmer.5_S1_d)
+```
+
+```{r save_coefficients_5_S1d, include=FALSE}
+
+coef_df = data.frame("intercept" = fixef(lmer.5_S1_d)[1], "slope" = fixef(lmer.5_S1_d)[2], row.names = "")
+write_csv(coef_df, "_stats/figure_5_S1d_coefs.csv")
+```
+
+Slope of `r format(fixef(lmer.5_S1_d)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.5_S1_d)[2], digits=2, nsmall=2)` (95\%-confidence interval) \newline
+n = `r nrow(tb %>% drop_na(spars, meanrate) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(spars, meanrate) %>% count(mid))` mice 
+
+\newpage
+# Figure 5-Supplement 1e
+## Relation between firing rate RMI and reliability RMI
+
+```{r fit_model_5_S1e}
+
+# Random intercept for neurons,
+# random intercept for experiments, nested in series,nested in mice
+lmer.5_S1e = lmer(rel ~ meanrate + (1 | uid) + (1 | mid/sid/eid), 
+                  data = tb %>% drop_na(rel, meanrate))
+
+display(lmer.5_S1e)
+anova(lmer.5_S1e)
+```
+
+```{r store_coefficients_5_S1e, include=F}
+
+coef_df = data.frame("intercept" = fixef(lmer.5_S1e)[1], "slope" = fixef(lmer.5_S1e)[2], row.names = "")
+write_csv(coef_df, "_stats/figure_5_S1e_coefs.csv")
+```
+
+Slope of `r format(fixef(lmer.5_S1e)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.5_S1e)[2], digits=2, nsmall=2)` (95\%-confidence interval) \newline
+n = `r nrow(tb %>% drop_na(rel, meanrate) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(rel, meanrate) %>% count(mid))` mice 
+
+\newpage
+# Figure 5-Supplement 1f
+## Relation between firing rate RMI and SNR RMI
+
+```{r fit_model_5_S1f}
+
+# Random intercept for neurons,
+# random intercept for experiments, nested in series
+lmer.5_S1f = lmer(snr ~ meanrate + (1 | uid) + (1 | sid/eid), 
+                  data = tb %>% drop_na(snr, meanrate))
+
+display(lmer.5_S1f)
+anova(lmer.5_S1f)
+```
+
+```{r store_coefficents_5_S1f, include=F}
+
+coef_df = data.frame("intercept" = fixef(lmer.5_S1f)[1], "slope" = fixef(lmer.5_S1f)[2], row.names = "")
+write_csv(coef_df, "_stats/figure_5_S1f_coefs.csv")
+```
+
+Slope of `r format(fixef(lmer.5_S1f)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.5_S1f)[2], digits=2, nsmall=2)` (95\%-confidence interval) \newline
+n = `r nrow(tb %>% drop_na(snr, meanrate) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(snr, meanrate) %>% count(mid))` mice 
+
+\newpage
+# Figure 5-Supplement 1g
+## Relation between firing rate RMI and peak width RMI
+
+```{r fit_model_5_S1g}
+
+# Random intercept for mice
+lmer.5_S1g = lmer(meanpkw ~ meanrate + (1 | mid), 
+                  data = tb %>% drop_na(meanpkw, meanrate))
+
+display(lmer.5_S1g)
+anova(lmer.5_S1g)
+```
+
+```{r store_coefficients_5_S1g, include=FALSE}
+
+coef_df = data.frame("intercept" = fixef(lmer.5_S1g)[1], "slope" = fixef(lmer.5_S1g)[2], row.names = "")
+write_csv(coef_df, "_stats/figure_5_S1g_coefs.csv")
+```
+
+Slope of `r format(fixef(lmer.5_S1g)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.5_S1g)[2], digits=2, nsmall=2)` (95\%-confidence interval) \newline
+n = `r nrow(tb %>% drop_na(meanpkw, meanrate) %>% count(uid))` neurons from `r nrow(tb %>% drop_na(meanpkw, meanrate) %>% count(mid))` mice 
+
+\newpage
+# Figure 5-Supplement 1h
+# Distributions of eye position variability, separated by locomotion state
+
+```{r read_data_5_S1h, include=FALSE}
+
+tib = get_data("../csv/ipos_st8.csv")
+
+# Turn locomotion state into a binary predictor
+tb = tib %>% mutate(run = ifelse(st8 == "run", 1, 0))
+
+# 'Standard deviation of eye position is not computed on a trial-by-trial basis. In the csv-file, 
+# this measure is therefore identical across trials, so we simply pull out the first trial of each experiment
+
+tbh = tb %>% select(mid, sid, eid, mse, run, std_xpos_cross) %>% distinct(mse, run, .keep_all = TRUE)
+
+```
+
+```{r fit_model_5_S1h}
+
+# Random intercept for series
+lmer.5_S1h = lmer(std_xpos_cross ~ run + (1 | sid), 
+                  data = tbh %>% drop_na(std_xpos_cross))
+
+display(lmer.5_S1h)
+anova(lmer.5_S1h)
+```
+
+```{r predicted_average_effect_5_S1h, include=F}
+
+mSit = fixef(lmer.5_S1h)[1]
+diffMeans = fixef(lmer.5_S1h)[2]
+mRun  = fixef(lmer.5_S1h)[1] + diffMeans
+```
+
+Locomotion: mean eye position standard deviation of `r format(mRun, digits=2, nsmall=2)` (95\%-confidence interval) \newline
+Sitting:  mean eye position standard deviation of `r format(mSit, digits=2, nsmall=2)` \newline
+n = `r nrow(tbh %>% drop_na(std_xpos_cross) %>% count(eid))` experiments from `r nrow(tbh %>% drop_na(std_xpos_cross) %>% count(mid))` mice  
+
+\newpage
+# Figure 5-Supplement 1i
+# Relation between locomotion effects on eye position variability and firing rate variability
+
+```{r read_data_5_S1i, include=FALSE}
+
+tib = get_data("../csv/iposmi.csv")
+```
+
+```{r fit_model_5_S1i}
+
+# Random effect of experiment, units partially crossed
+lmer.5_S1i = lmer(relrmi ~ iposrmi + (1 | uid), 
+                  data = tib %>% drop_na(relrmi, iposrmi))
+
+display(lmer.5_S1i)
+anova(lmer.5_S1i)
+```
+
+```{r save_coefficients_5_S1i, include=FALSE}
+
+coef_df = data.frame("intercept" = fixef(lmer.5_S1i)[1], "slope" = fixef(lmer.5_S1i)[2], row.names = "")
+write_csv(coef_df, "_stats/figure_5_S1i_coefs.csv")
+
+```
+
+
+Slope of `r format(fixef(lmer.5_S1i)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.5_S1i)[2], digits=2, nsmall=2)` (95\%-confidence interval)
+
+```{r assess_coefficient, include = FALSE}
+
+foo = tib %>% filter(iposrmi != "NA")
+sd_ipos = sd(foo$iposrmi)
+
+# expected difference in rel RMI corresponding to a 1-sd difference in eye pos signm
+exp_diff = fixef(lmer.5_S1i)[2] * sd_ipos
+
+print(paste("The expected difference in relRMI corresponding to a 1 standard deviation difference in eyePosSigmaRMI is ", format(exp_diff, digits=2)))
+
+```
+Expected difference in reliability RMI corresponding to a 1 standard deviation difference in eye position $\sigma$ RMI is `r format(exp_diff, digits=2, nsmall=2)`, the standard deviation of the residuals is `r format(sigma(lmer.5_S1i), digits=2, nsmall=2)`. \newline \newline
+n = `r nrow(tib %>% drop_na(relrmi, iposrmi) %>% count(uid))` neurons from `r nrow(tib %>% drop_na(relrmi, iposrmi) %>% count(mid))` mice  
+

R/figure_5_S2.Rmd

@@ -0,0 +1,170 @@
+---
+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

R/figure_6.Rmd

@@ -0,0 +1,847 @@
+---
+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
+
+

R/figure_6_S1.Rmd

@@ -0,0 +1,214 @@
+---
+title: "Spacek et al., 2021, Figure 6-Supplement 1"
+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 == "grt") %>% select(mid, sid, eid, uid, mseu, mi, meanrate, meanburstratio)
+```
+
+```{r tidy_for_a1, 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 6-S1a$_1$
+## Comparing RMI during suppression against 0
+
+```{r, fit_model_a1_1}
+
+# Fixed effect intercept only, 
+# random intercept for neurons,
+# random intercept for experiments, nested within series
+lmer.a1.1 = lmer(suppressionrmi ~ 1 + (1 | uid) + (1 | sid/eid), 
+                  data = tbw %>% drop_na(suppressionrmi))
+
+display(lmer.a1.1)
+```
+
+## Mean firing rate RMI
+Suppression: RMI = `r format(fixef(lmer.a1.1)[1], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.a1.1)[1], digits=2, nsmall=2)` \newline
+n = `r nrow(tbw %>% drop_na(suppressionrmi) %>% count(uid))` neurons from `r nrow(tbw %>% drop_na(suppressionrmi) %>% count(mid))` mice
+
+\newpage
+# Figure 6-S1a$_2$
+## Comparing RMI during suppression against 0
+
+```{r tidy_for_a2, 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)
+```
+
+```{r fit_model_a2_1}
+
+# LMM gives singular fits, even with a single random intercept for neurons. 
+# We do ordinary regression instead:
+lm.a2.1 = lm(suppressionrmi ~ 1, data = tbw %>% drop_na(suppressionrmi))
+
+display(lm.a2.1)
+```
+
+## Mean burst ratio RMI
+Suppression: RMI = `r format(summary(lm.a2.1)$coefficients[1], digits=2, nsmall=2)` $\pm$ `r format(2 * summary(lm.a2.1)$coefficients[2], 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 6-S1b$_1$
+## Slope of regression line
+
+```{r tidy_for_b1_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)
+```
+
+```{r fit_model_b1_1}
+
+# Random intercept for neurons,
+# random intercept for experiments
+lmer.b1.1 = lmer(runfmi ~ sitfmi + (1 | uid) + (1 | eid), 
+                  data = tbw %>% drop_na(runfmi, sitfmi))
+
+display(lmer.b1.1)
+anova(lmer.b1.1)
+```
+
+``` {r store_coefficients_b1_1, include=FALSE}
+
+coef_df = data.frame("intercept" = fixef(lmer.b1.1)[1], "slope" = fixef(lmer.b1.1)[2], row.names = "")
+write_csv(coef_df, "_stats/figure_6_S1b1_coefs.csv")
+```
+
+Slope of `r format(fixef(lmer.b1.1)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.b1.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 6-S1b$_2$
+## Slope of regression line
+
+```{r tidy_for_b2_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")
+```
+
+```{r fit_model_b2_1}
+
+# Random intercept for neurons,
+# random intercept for experiments, nested in series
+lmer.b2.1 = lmer(runfmi ~ sitfmi + (1 | uid) + (1 | sid/eid), 
+                  data = tbw %>% drop_na(runfmi, sitfmi))
+
+display(lmer.b2.1)
+anova(lmer.b2.1)
+```
+
+```{r store_coefficients_b2_1, include=FALSE}
+
+# store results in file
+coef_df = data.frame("intercept" = fixef(lmer.b2.1)[1], "slope" = fixef(lmer.b2.1)[2], row.names = "")
+write_csv(coef_df, "_stats/figure_6_S1b2_coefs.csv")
+```
+
+Slope of `r format(fixef(lmer.b2.1)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.b2.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 6-S1c$_1$
+## Slope of regression line
+
+```{r tidy_for_c1_1, 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_c1_1}
+
+# Random intercept for neurons,
+# random intercept for experiments
+lmer.c1.1 = lmer(fmi ~ rmi + (1 | uid) + (1 | eid), 
+                data = tbw %>% drop_na(fmi, rmi))
+
+display(lmer.c1.1)
+anova(lmer.c1.1)
+```
+
+```{r, store_coefficients_c1_1, include=FALSE}
+
+# store results in file
+coef_df = data.frame("intercept" = fixef(lmer.c1.1)[1], "slope" = fixef(lmer.c1.1)[2], row.names = "")
+write_csv(coef_df, "_stats/figure_6_S1c1_coefs.csv")
+```
+
+Slope of `r format(fixef(lmer.c1.1)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.c1.1)[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 6-S1c$_2$
+## Slope of regression line
+
+```{r, tidy_for_c2_1, 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_c2_1}
+
+# Random intercept for neurons, 
+# random intercept for experiments, nested in series
+lmer.c2.1 = lmer(fmi ~ rmi + (1 | uid) + (1 | sid/eid), 
+                data = tbw %>% drop_na(fmi, rmi))
+
+display(lmer.c2.1)
+anova(lmer.c2.1)
+```
+
+```{r, store_coefficients_c2_1, include=FALSE}
+
+coef_df = data.frame("intercept" = fixef(lmer.c2.1)[1], "slope" = fixef(lmer.c2.1)[2], row.names = "")
+write_csv(coef_df, "_stats/figure_6_S1c2_coefs.csv")
+```
+
+Slope of `r format(fixef(lmer.c2.1)[2], digits=2, nsmall=2)` $\pm$ `r format(2 * se.fixef(lmer.c2.1)[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
+

R/get_data.R

@@ -0,0 +1,56 @@
+get_data <- function(csv_file_name) {
+  # Reads data from csv file into a tibble, and creates IDs for mice, series, experiments, and neurons
+  
+  tib = read_csv(csv_file_name)
+  
+  if (has_name(tib, 'mseu')) {
+  
+    # add mouse ID - this is based on m(ouse) only
+    tib <- tib %>% mutate(mid = as.numeric(as.factor(str_split_fixed(tib$mseu, '_s', 2)[,1])))
+    
+    # add series ID - this is based on m and s(eries)
+    tib <- tib %>% mutate(sid = as.numeric(as.factor(str_split_fixed(tib$mseu, '_e', 2)[,1])))
+    
+    # add exp ID - this is based on m, s, and e(xperiment)
+    tib <- tib %>% mutate(eid = as.numeric(as.factor(str_split_fixed(tib$mseu, '_u', 2)[,1])))
+    
+    # add unit ID - this is based on m, s, and u(nit), ignoring e
+    tib <- tib %>% mutate(uid = as.numeric(as.factor(paste(str_split_fixed(tib$mseu, '_e', 2)[,1], '_u' , str_split_fixed(tib$mseu, '_u', 2)[,2], sep=""))))
+    
+    # put iDs up front
+    tib <- tib %>% select(mid, sid, eid, uid, everything())
+  
+    # If we don't have 'experiments', we need to do a slightly different parsing
+  } else if (has_name(tib, 'msu')) {
+    
+    # add mouse ID - this is based on m only
+    tib <- tib %>% mutate(mid = as.numeric(as.factor(str_split_fixed(tib$msu, '_s', 2)[,1])))
+    
+    # add series ID - this is based on m and s
+    tib <- tib %>% mutate(sid = as.numeric(as.factor(str_split_fixed(tib$msu, '_u', 2)[,1])))
+    
+    # add unit ID - this is based on m, s, and u
+    tib <- tib %>% mutate(uid = as.numeric(as.factor(tib$msu)))
+    
+    # put IDs up front
+    tib <- tib %>% select(mid, sid, uid, everything()) # put id up front
+    
+    # We don't have units - this requires a slightly differnt parsing
+  } else if (has_name(tib, 'mse')) {
+    
+    # add mouse ID
+    tib <- tib %>% mutate(mid = as.numeric(as.factor(str_split_fixed(tib$mse, '_s', 2)[,1])))
+    
+    # add series ID
+    tib <- tib %>% mutate(sid = as.numeric(as.factor(str_split_fixed(tib$mse, '_e', 2)[,1])))
+    
+    # add experiment ID
+    tib <- tib %>% mutate(eid = as.numeric(as.factor(tib$mse)))
+    
+    # put IDs up front
+    tib <- tib %>% select(mid, sid, eid, everything()) 
+    
+  }
+
+  return(tib)
+}