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