123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297 |
- ---
- 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
|