---
title: "Bayesian regression models"
author: "Bruno Nicenboim / Shravan Vasishth"
date: "`r Sys.Date()`"
output:
  bookdown::beamer_presentation2:
    theme: "metropolis"
    keep_tex: yes 
    latex_engine: xelatex
    slide_level: 2
    incremental: no
    number_sections: true
    includes: 
      in_header: top-matter.tex
fontsize: 12pt
classoption: aspectratio=169
bibliography: ["BayesCogSci.bib", "packages.bib"]
---




# A first linear model: Does attentional load affect pupil size? {#sec:pupil}
<!-- https://bookdown.org/yihui/rmarkdown/beamer-presentation.html -->

```{r setup, include=FALSE, cache=FALSE}
knitr::opts_chunk$set(tidy = "styler",
                      cache=TRUE,
                      size = "small"
                      )

## #Hack to avoid compatibility issues with tikz
## knitr::knit_hooks$set(document = function(x) {
##     sub('\\usepackage{color}', '\\usepackage[table]{xcolor}', x, fixed = TRUE)
## })


## Reduces the size of the font in code
## https://stackoverflow.com/questions/25646333/code-chunk-font-size-in-rmarkdown-with-knitr-and-latex
def.chunk.hook  <- knitr::knit_hooks$get("chunk")
knitr::knit_hooks$set(chunk = function(x, options) {
    x <- def.chunk.hook(x, options)
    ifelse(options$size != "normalsize", paste0("\n \\", options$size,"\n\n", x, "\n\n \\normalsize"), x)
})

options(
    htmltools.dir.version = FALSE,
    formatR.indent = 2,
    width = 75,
    digits = 2,
    signif =2,
    warnPartialMatchAttr = FALSE,
    warnPartialMatchDollar = FALSE,
    # Don't use scientific notation:
    scipen=10000,
    # tibbles:
    tibble.width = Inf,
    tibble.print_max = 5,
    tibble.print_min = 5
)
library(papaja)
library(bookdown)
ggplot2::theme_set(ggplot2::theme_light())
library(partitions)
```




```{r load-internal, cache =FALSE, message=FALSE, echo = FALSE}
set.seed(42)
library(MASS)
##be careful to load dplyr after MASS
library(dplyr)
library(tidyr)
library(purrr)
library(readr)
library(extraDistr)
library(ggplot2)
library(brms)
library(rstan)
## Save compiled models:
rstan_options(auto_write = TRUE)
## Parallelize the chains using all the cores:
options(mc.cores = parallel::detectCores())
library(bayesplot)

```



---



<!-- We'll look at the effect of cognitive processing on human pupil size to illustrate the use of Bayesian linear regression models. Although pupil size is mostly related to the amount of light that reaches the retina or the distance to a perceived object, pupil sizes are also systematically influenced by cognitive processing: It has been found that increased cognitive load leads to an increase in the pupil size [for a review, see @mathotPupillometryPsychologyPhysiology2018]. -->

### Data:

One participant's pupil size of the control experiment of @wahnPupilSizesScale2016 averaged by trial

### Task:

A participant covertly tracked between zero and five objects among several randomly moving objects on a computer screen;  multiple object tracking--MOT-- [@pylyshynTrackingMultipleIndependent1988] task
<!-- several objects appear in the screen, and a subset of them -->
<!-- are indicated as "targets" at the beginning. Then, the objects start moving randomly across the screen and become indistinguishable. After several seconds, the objects stop moving and the participant needs to indicate which objects were the targets. See also Figure \@ref(fig:mot). -->


### Research question: 

How does the number of moving objects being tracked (attentional load) affect pupil size?

---

(ref:mot) Flow of events in a trial where two objects needs to be tracked. Adapted from @Blumberg2015; licensed under CC BY 4.0.

```{r mot, fig.cap = "(ref:mot)", out.width = "80%", echo = FALSE, fig.align = "center"}
knitr::include_graphics("cc_figure/MOT.png", dpi =1000)
```



----

<!-- We will model pupil size as normally distributed, because we are not expecting a skew, and we have no further information available about the distribution of pupil sizes. (Notice that pupil sizes cannot be of size zero or negative, so we know for sure that this choice is not exactly right.) For simplicity, we are also going to assume a linear relationship between load and the pupil size.  -->

### Assumptions:

1.  There is some average pupil size represented by $\alpha$.
2.  The increase of attentional load  has a linear relationship with pupil size, determined by $\beta$.
3.  There is some noise in this process, that is, variability around the true pupil size i.e., a scale, $\sigma$.
4.  The noise is normally distributed.

## Formal model

### Likelihood for each observation $n$:

\begin{equation}
p\_size_n \sim Normal(\alpha + c\_load_n \cdot \beta,\sigma)
\end{equation}

where $n$ indicates the observation number with $n = 1 \ldots N$


<!-- This means that the formula that we'll use in `brms` will be  `p_size ~ 1 + c_load`, where `1` represents the intercept, $\alpha$, which doesn't depend on a covariate or predictor, and `c_load` is our covariate  that is multiplied by $\beta$.  We will generally indicate with the prefix `c_`, that a covariate (in this case load) is centered (i.e., we subtract from each  value the mean of all values). If load is centered,  the intercept represents the pupil size at the average load in the experiment (because at the average load, the centered load is zero, and then $\alpha + 0 \cdot \beta$). Alternatively, if the load would not have been centered (i.e., starts with no load, then one, two, etc), then the intercept would represent the pupil size when there is no load. Although this formula would be enough to fit a frequentist model with `lm(p_size ~ 1 + c_load, dataset)`, when we fit a Bayesian model, we have to specify  priors for each of the parameters. -->


### How do we decide on priors?

## Priors

*  pupil sizes range between 2 and 5 millimeters,
*  but the Eyelink-II eyetracker measures the pupils in arbitrary units [@hayesMappingCorrectingInfluence2016]
* we either need estimates from a previous  analysis or look at some measures of pupil sizes

---

### Pilot data:
Some measurements of the same participant with no attentional load for the first 100ms, each 10 ms,  in `pupil_pilot.csv`: 

```{r, message = FALSE}
df_pupil_pilot <- read_csv("./data/pupil_pilot.csv")
df_pupil_pilot$p_size %>% summary()
```

## Prior for $\alpha$

<!-- With this information we can set a regularizing prior for $\alpha$. We center the prior around 1000 to be in the right order of magnitude.^[The average pupil size will probably be higher than 800, since this measurement was with no load, but, in any case, the exact number won't matter, any mean between 500-1500 would be fine if the standard deviation is large.] Since we don't know how much pupil sizes are going to vary by load yet, we include a rather wide prior by defining it as a normally distribution and setting its standard deviation as $500$. -->


\begin{equation}
\alpha \sim Normal(1000, 500) 
\end{equation}


<!-- Given that our  covariate load is centered, with the prior for $\alpha$, we are saying that we suspect that the  -->

### Meaning:

We expect that the average pupil size  for the average load in the experiment would  be in a 95% central interval limited by approximately $1000 \pm 2 \cdot 500 = [20, 2000]$ units:


```{r}
c(qnorm(.025, 1000,500), qnorm(.975, 1000, 500))
```


## Prior for $\sigma$

<!-- We know that the measurements of the pilot data are strongly correlated because they were taken together just some milliseconds apart. For this reason, they won't tell us how much the pupil size can vary. We set up a quite weak prior for $\sigma$ that encodes our lack of precise information: $\sigma$  is surely larger than zero and has to be in the order of magnitude of the pupil size with no load. -->

\begin{equation}
\sigma \sim Normal_+(0, 1000)
\end{equation}

### Meaning:

We expect that the standard deviation of the pupil sizes should be in the following 95% interval. 
```{r}
c(qtnorm(.025, 0, 1000, a = 0),
  qtnorm(.975, 70,1000, a = 0))
```


## Prior for $\beta$

<!-- We still need to set a prior for $\beta$, the change in pupil size produced by the attentional load. Given that pupil size changes are not easily perceptible (we don't see them in our day-to-day life), we expect them to be much smaller than the pupil size,  so we use the following prior: -->

\begin{equation}
\beta \sim Normal(0, 100)
\end{equation}


### Meaning:

We don't really know if the attentional load will increase or even decrease the pupil size, but we are only saying that one unit of load  will potentially change the pupil size consistently with the following 95% interval:

```{r}
c(qnorm(.025, 0,100), qnorm(.975, 0,100))
```
<!-- That is, we don't expect changes in size that increase or decrease the pupil size in more than 200 units. -->


## Fitting the model

```{r, message = FALSE, size = "scriptsize"}
df_pupil_data <- read_csv("data/pupil.csv")
df_pupil_data <- df_pupil_data %>%
    mutate(c_load = load - mean(load))
df_pupil_data
```

## Specifying the model in brms

```{r fitpupil, message = FALSE, results = "hide", size = "scriptsize"}
fit_pupil <- brm(p_size ~ 1 + c_load,
                 data = df_pupil_data,
                 family = gaussian(),
                 prior = c(
                     prior(normal(1000, 500), class = Intercept),
                     prior(normal(0, 1000), class = sigma),
                     prior(normal(0, 100), class = b, coef = c_load)
                 )) 
```

<!-- The only difference from our previous models is that we now have a predictor in the formula and in the priors. Priors for predictors are indicated with `class = b`, and the specific predictor with `coef = c_load`. If we want to set the same priors to different predictors we can omit the argument `coef`. We can remove the `1` of the formula, and `brm()` will fit the exact same model as when we specify `1` explicitly. If we really want to remove the intercept we indicate this with `0 +...` or `-1 +...`. See also the box \@ref(thm:intercept) for more details about the treatment of the intercepts by `brms`. -->

<!-- We can inspect the output of our model now: -->

---

```{r}
plot(fit_pupil)
```

----

```{r, size = "scriptsize"}
fit_pupil
```



##  How to communicate the results?

```{r, echo = FALSE}
mean_load <- posterior_summary(fit_pupil)["b_c_load","Estimate"] %>%
    round(2)
load_l <- posterior_summary(fit_pupil)["b_c_load","Q2.5"]%>%
    round(2)
load_h <- posterior_summary(fit_pupil)["b_c_load","Q97.5"]%>%
    round(2)
```


### Research question: 
"What is the effect of attentional load on the participant’s pupil size?" 

We'll need to examine what happens with $\beta$ (`c_load`):


##  How to communicate the results?

* The most likely values of $\beta$ will be around the mean of the posterior, `r mean_load`, and we can be 95% certain that the true value of $\beta$ *given the model and the data* lies between `r load_l` and `r load_h`.
* We see that as the attentional load increases, the pupil size of the participant becomes larger. 

---


### How likely it is that the pupil size increased rather than decreased?
<!-- we can examine the proportion of samples above zero. (Notice that the intercept and the slopes, are always preceded by `b_` in `brms`.  One can see all the names of parameters being estimated with `parnames()`.) -->

```{r}
mean(posterior_samples(fit_pupil)$b_c_load > 0)
```

\color{red}
Take into account that this probability ignores the possibility of the participant not being affected at all by the manipulation, this is because $P(\beta=0)=0$.


## Descriptive adequacy {#sec:pupiladq}

<!-- Our model converged and we obtained a posterior distribution, there is, however, no guarantee that our model was adequate to represent our data. We can use posterior predictive checks to verify this.  -->

<!-- Sometimes it's useful, build our own posterior predictive check to visualize the fit of our model, as opposed to use the `pp_check` functions as we did before in section \@ref(sec:ppd). For example, here we use `posterior_predict()` to generate 1000 posterior predictive distributions, and we convert them from an array to a long data frame. -->

```{r, size = "scriptsize"}
# we start from an array of 1000 samples by 41 observations
df_pupil_pred <- posterior_predict(fit_pupil, nsamples = 1000) %>%
    # we convert it to a list of length 1000, with 41 observations in each element:
    array_branch(margin = 1) %>%
    # We iterate over the elements (the predicted distributions)
    # and we convert them into a long data frame similar to the data,
    # but with an extra column `iter` indicating from which iteration
    # the sample is coming from.
    map_dfr( function(yrep_iter) {
        df_pupil_data %>%
            mutate(p_size = yrep_iter)
    }, .id = "iter") %>%
    mutate(iter = as.numeric(iter))
```

----

<!-- Then we plot 100 of the densities of the predicted distributions in blue, and the distribution of our data in black for the five levels of load in Figure \@ref(fig:postpreddens). We don't have enough data to derive a strong conclusion: Notice that both the predictive distributions and our data look very wide, and it hard to tell if the distribution of the observations could have been generated by our model. For now we can say that it doesn't look too bad. -->

\vspace{.1in}

(ref:postpreddens) The plot shows 100 predicted distributions in blue density plots, the distribution of pupil size data in black density plots, and the observed pupil sizes in black dots for the five levels of attentional load.

```{r postpreddens, fig.cap ="(ref:postpreddens)" , message= FALSE, fig.height =3.5, size = "scriptsize", tidy=FALSE }
df_pupil_pred %>% filter(iter < 100) %>%
    ggplot(aes(p_size, group=iter)) + 
  geom_line(alpha = .05, stat="density", color = "blue") +
    geom_density(data=df_pupil_data, aes(p_size), inherit.aes = FALSE, size =1)+
    geom_point(data=df_pupil_data, aes(x=p_size, y = -0.001), alpha =.5, inherit.aes = FALSE) +
    coord_cartesian(ylim=c(-0.002, .01))+ facet_grid(load ~ .) 
```

## Distribution of statistics

<!-- We can instead look  at the  distribution of a  statistic, such as mean pupil size by load: -->

```{r, size = "scriptsize"}
# predicted means:
df_pupil_pred_summary <- df_pupil_pred %>%
    group_by(iter, load) %>%
    summarize(av_p_size = mean(p_size))
# observed means:
(df_pupil_summary <- df_pupil_data %>%
    group_by(load) %>%
    summarize(av_p_size = mean(p_size)))
```

---


(ref:postpredmean) Distribution of posterior predicted means in gray and observed pupil size means in black lines by load.

```{r postpredmean, fig.cap ="(ref:postpredmean)", message= FALSE, fig.height =4, size = "scriptsize"}
ggplot(df_pupil_pred_summary, aes(av_p_size)) +
    geom_histogram(alpha=.5)+
    geom_vline(aes(xintercept= av_p_size),data= df_pupil_summary)+
    facet_grid(load ~ .)
```

---

* the observed means for no load and for a load of two are falling in the tails of the distributions. 
*  the data might be indicating that the relevant difference is between (i) no load, (ii) a load between two and three, and then (iii) a load of four, and (iv) of five.
* but beware of overinterpreting noise.


## Value of posterior predictive distributions

* If we look hard enough, we'll find failures of descriptive adequacy.^[all models are wrong]
* Posterior predictive accuracy can be used to generate new hypotheses and to compare different models.

----

### \color{blue} Exercises

4.6.1.1 Our priors for this experiment were quite arbitrary. How do the prior predictive distributions look like? Do they make sense?

4.6.1.2 Is our posterior distribution sensitive to the priors that we selected? Perform a sensitivity analysis to find out whether the posterior is affected by our choice of prior for the $\sigma$.

4.6.1.3 Our dataset includes also a column that indicates the trial number. Could it be that trial has also an effect on the pupil size? As in `lm`, we indicate another main effect with a `+` sign. How would you communicate the new results?

# Log-normal model: Does trial affect reaction times? {#sec:trial}


---


We revisit the small experiment, where a participant repeatedly pressed the space bar as fast as possible, without paying attention to the stimuli.  

### New research question: 

Does the participant tend to speedup (practice effect) or slowdown (fatigue effect)? 

## Formal model

### Likelihood:

\begin{equation}
rt_n \sim LogNormal(\alpha + c\_trial_n \cdot \beta,\sigma)
(\#eq:rtloglik)
\end{equation}

### Priors 

\begin{equation}
\begin{aligned}
\alpha &\sim Normal(6, 1.5) \\
\sigma &\sim Normal_+(0, 1)\\
\beta &\sim \ldots
\end{aligned}
\end{equation}

## Prior for $\beta$

<!-- We still need a prior for $\beta$, but notice that effects are multiplicative rather than additive when we assume a log-normal likelihood  and that means that we need to take into account $\alpha$ in order to interpret $\beta$; see Box \@ref(thm:lognormal). We are going to try to understand how all our priors interact together generating some prior predictive distributions. We start with the following  prior centered in zero, a prior agnostic regarding the direction of the effect, which allows for both a slowdowns ($\beta>0$) or a speedups ($\beta<0$): -->

\begin{equation}
\beta \sim Normal(0, 1)
\end{equation}

----

We edit our `normal_predictive_distribution_fast` from section and make it log-normal and dependent on trial:

```{r, size="scriptsize", tidy= FALSE}
lognormal_model_pred <- function(alpha_samples,
                                 beta_samples,
                                 sigma_samples,
                                 N_obs) {
    # pmap extends map2 (and map) for a list of lists:
    pmap_dfr(list(alpha_samples, beta_samples, sigma_samples),
             function(alpha, beta, sigma) {
                 tibble(
                     trialn = seq_len(N_obs),
                     # we center trial:
                     c_trial = trialn - mean(trialn),
                     # we change the likelihood: 
                     # Notice rlnorm and the use of alpha and beta
                     rt_pred = rlnorm(N_obs, alpha + c_trial * beta, sigma))
             }, .id = "iter") %>%
    # .id is always a string and needs to be converted to a number
        mutate(iter = as.numeric(iter))}
```

-----

This is our first attempt for a prior predictive distribution:

```{r, size ="scriptsize"}
N_obs <- 361
N <- 800
alpha_samples <- rnorm(N, 6, 1.5)
sigma_samples <- rtnorm(N, 0, 1, a =0)
beta_samples <- rnorm(N, 0, 1)
prior_pred <- lognormal_model_pred(
    alpha_samples = alpha_samples,
    beta_samples = beta_samples, 
    sigma_samples = sigma_samples,
    N_obs = N_obs)
```

---------

```{r, size ="scriptsize"}
(median_effect <-
     prior_pred %>%
     group_by(iter) %>%
     mutate(diff = rt_pred - lag(rt_pred)) %>%
     summarize(
         median_rt = median(diff, na.rm = TRUE)
 ))
```

<!-- We plot it in Figure \@ref(fig:priorbeta), and as expected is center in zero (as our prior), but we see that the distribution of possible medians for the effect  is too spread and includes values that are too extreme. -->

---

\vspace{.1in}

(ref:priorbeta) Prior predictive distribution of the median effect  of the log-normal model with $\beta \sim Normal(0, 1)$.

```{r priorbeta,fig.cap="(ref:priorbeta)", message = FALSE, fig.height =4.5, size = "scriptsize" }
median_effect %>%
    ggplot(aes(median_rt)) +
    geom_histogram()
```


## Another prior for $\beta$

\begin{equation}
\beta \sim Normal(0, .01)
\end{equation}

<!-- We repeat the same procedure with $\beta \sim Normal(0,.01)$, and we plot it in Figure \@ref(fig:priorbeta2). The prior predictive distribution shows us that the prior is still quite vague, it is, howeverm at least in the right order of magnitude. Notice that we are using a distribution of medians because they are less affected by the variance in the posterior predicted distribution; distributions of means will have much more spread. If we want to make the distribution of means more realistic, we would also need to find a more accurate prior for the scale, $\sigma$. -->


```{r, echo = FALSE, size = "scriptsize"}
beta_samples2 <- rnorm(800, 0, .01)
prior_pred2 <- lognormal_model_pred(
    alpha_samples = alpha_samples,
    beta_samples = beta_samples2, 
    sigma_samples = sigma_samples,
    N_obs = N_obs)
```


(ref:priorbeta2) Prior predictive distribution of the median effect  of the log-normal model with $\beta \sim Normal(0, .01)$.

```{r priorbeta2,fig.cap="(ref:priorbeta2)", message = FALSE, echo = FALSE, fig.height = 3.8, size = "scriptsize" }
prior_pred2 %>%
    group_by(iter) %>%
    mutate(diff = rt_pred - lag(rt_pred)) %>%
    summarize(
        median_rt = median(diff, na.rm = TRUE)
    ) %>%
    ggplot(aes(median_rt))+
    geom_histogram()
```

## Prior selection

Prior selection might look daunting and a lot of work. However...

* priors can be informed by the estimates from previous experiments;
* this work is usually done only the first time we encounter an experimental paradigm;
* we will generally use very similar (or identical priors) for analyses dealing with the same type of task;
* when in doubt, do a sensitivity analysis.



## Fitting the  model

<!-- We are now relatively satisfied with the priors for our model, and we can fit the data with `brms`. Notice that we need to specify that the family is `lognormal()`. -->



```{r,  message = FALSE, size= "scriptsize"}
df_noreading_data <- read_csv("./data/button_press.csv")
df_noreading_data <- df_noreading_data %>%
    mutate(c_trial = trialn - mean(trialn))
fit_press_trial <- brm(rt ~ 1 + c_trial,
  data = df_noreading_data,
  family = lognormal(),
  prior = c(
    prior(normal(6, 1.5), class = Intercept),
    prior(normal(0, 1), class = sigma),
    prior(normal(0, .01), class = b, coef = c_trial)
  )
)
```

---

<!-- Instead of printing out the complete output from the model, look at the estimates from the posteriors for the parameters $\alpha$, $\beta$, and $\sigma$. Notice that these parameters are on the log scale: -->

```{r,size="scriptsize"}
posterior_summary(fit_press_trial)[,c("Estimate","Q2.5","Q97.5")]
```


---

```{r}
plot(fit_press_trial)
```

<!-- Next, we turn to the question of what we can report as our results, and what we can conclude from the data. -->

##  How to communicate the results?

```{r, echo=FALSE, results="hide"}
alpha_samples<- posterior_samples(fit_press_trial)$b_Intercept
beta_samples<- posterior_samples(fit_press_trial)$b_c_trial

beta_ms<- exp(alpha_samples) - exp(alpha_samples-beta_samples)

beta_msmean <- round(mean(beta_ms),5)
beta_mslow <- round(quantile(beta_ms,prob=0.025),5)
beta_mshigh <- round(quantile(beta_ms,prob=0.975),5)

beta_mean <- round(mean(beta_samples),5) %>% format() 
beta_low <- round(quantile(beta_samples,prob=0.025),5) %>% format() 
beta_high <- round(quantile(beta_samples,prob=0.975),5) %>% format() 
```

<!-- As shown above, the first step is to summarize the posteriors in a table or graphically (or both).  If the research relates to the effect estimated by the model,  -->

### We focus on the effect of trial:

* $\hat\beta = `r beta_mean`$, 95% CrI = $[ `r beta_low` , `r beta_high` ]$.

* But in most cases, the effect is easier to interpret in milliseconds. 

<!-- We can transform the estimates back to the millisecond scale from the log scale, but we need to take into account that the scale is not linear, and that the effect between two button presses will differ depending on where we are in the experiment. -->

---

We calculate an estimate if we consider the difference between reaction times in a trial at the middle of the experiment (when the centered trial number is zero) and the previous one (when the centered trial number is minus one).

```{r, size="scriptsize"}
alpha_samples<- posterior_samples(fit_press_trial)$b_Intercept
beta_samples<- posterior_samples(fit_press_trial)$b_c_trial
effect_middle_ms <- exp(alpha_samples) - 
  exp(alpha_samples - 1* beta_samples)
## ms effect in the middle of the expt (mean trial vs. mean trial - 1 ) 
c(mean = mean(effect_middle_ms), quantile(effect_middle_ms, c(.025,.975)))
```

---

Alternatively we consider the difference between the second trial and the first one:

```{r, size="scriptsize"}
first_trial <- min(df_noreading_data$c_trial)
second_trial <- min(df_noreading_data$c_trial) +1
effect_beginning_ms <- exp(alpha_samples+  second_trial * beta_samples) -
    exp(alpha_samples+  first_trial * beta_samples)
## ms effect from first to second trial:
c(mean = mean(effect_beginning_ms), quantile(effect_beginning_ms, c(.025,.975)))
```

There is a slowdown in both cases.

## Reporting results

### We can

* present the posterior mean and the 95% credible interval;
* assess if the observed estimates are consistent with the prediction from our theory;
* assess the practical relevance of the effect for the research question; (only after 100 button presses we see a slowdown of `r round(mean(effect_middle_ms),2) * 100` ms on average ($`r mean(effect_middle_ms)` \cdot 100$), with a 95% credible interval ranging from `r quantile(effect_middle_ms, .025)*100` to `r quantile(effect_middle_ms, .975)*100`);
* establish the presence or absence of an effect (Bayes factor)

<!-- Sometimes, researchers are only interested in establishing that there is an effect; the magnitude and uncertainty of the estimate is of secondary interest. Here, the goal is to argue that there is  **evidence** of a slowdown. The word evidence has a special meaning in statistics [@Royall], and in null hypothesis significance testing, a likelihood ratio test is the standard way to argue that one has evidence for an effect. In the Bayesian data analysis context, a Bayes factor hypothesis test must be carried out. We’ll come back to this issue in the model comparison section \@ref(sec:?). -->

----

### \color{blue} Exercises

4.6.2.1 Estimate the slowdown in milliseconds between the last two times the subject pressed the space bar in the experiment.

4.6.2.2 How would you change your model (keeping the log-normal likelihood) so that it includes centered log-transformed trial numbers or square-root-transformed trial numbers (instead of centered trial numbers)? Does the effect in milliseconds change?

# Logistic regression: Does set size affect free recall? {#sec:logistic}

---

We'll look at the capacity limit of working memory to illustrate one special case of GLMs, logistic regression.

Subset of the data of @oberauerWorkingMemoryCapacity2019:

**Data**
One participants recall success (1 success, 0 failure)

**Task:**
word lists of varying lengths (2, 4, 6, and 8 elements** were presented, and the participant was asked to recall a word given its position on the list

**Research question:**
How does the number of items to be held in working memory affects recall accuracy? 

<!-- It is well established that as the number of items to be held in working memory increases, performance, that is accuracy, decreases [among others @oberauerkliegel2001]. We will investigate whether we can establish this finding with data from only one subject. -->

---

(ref:oberauer) Flow of events in a trial with memory set size 4 and free recall. Adapted from @oberauerWorkingMemoryCapacity2019; licensed under CC BY 4.0.

```{r oberauer, fig.cap = "(ref:oberauer)", fig.height =5, echo = FALSE}
knitr::include_graphics("cc_figure/fig1_oberauer_2019_modified.png")
```

----

```{r,  message = FALSE, warning = FALSE, size = "scriptsize"}
df_recall_data <- read_csv("./data/PairsRSS1_all.csv") %>%
    # We ignore the type of incorrect responses (the focus of the paper)
    mutate(correct = if_else(response_category ==1, 1, 0)) %>%
    # and we only use the data from the free recall task:
    # (when there was no list of possible responses)
    filter(response_size_list + response_size_new_words == 0) %>%
    # We select one subject
    filter(subject == 10) %>%
    mutate(c_set_size = set_size - mean(set_size)) %>%
    select(subject, set_size, c_set_size, correct, trial)
```

----

```{r,  message = FALSE, warning = FALSE, size = "scriptsize"}
# Set sizes in the dataset:
df_recall_data$set_size %>%
    unique
# Trials by set size
df_recall_data %>%
    group_by(set_size) %>%
    count()
```

-----

<!-- The data look like this: the column `correct` records the 0 (incorrect) or 1 (correct) responses, and the column `c_set_size` records the centered memory set size; these latter scores have continuous values -3, -1, 1, and 3. These continuous values are centered versions of 2, 4, 6, and 8. -->

```{r, size = "scriptsize"}
df_recall_data
```

<!-- We want to model the trial by trial accuracy and examine whether the probability of recalling a word is related to the number of words in the set that the subject needs to remember. -->

## The likelihood for the logistic regression model

Recall that the Bernoulli likelihood generates a 0 or 1 response with a particular probability $\theta$ (here N = 10 trials with 50% chances of getting a one):

```{r}
# We use as.numeric to get zeros and ones rather than FALSE and TRUE
rbernoulli(n = 10, p = 0.5) %>% as.numeric()
```


## Formal model

**The likelihood for each observation $n$:**

\begin{equation}
correct_n \sim Bernoulli(\theta_n)
(\#eq:bernoullilik)
\end{equation}


- $\theta_n$ is bounded to be between 0 and 1 

*How do we fit a regression model?*

<!-- (it is a probability), we cannot just fit a regression model using the normal or lognormal likelihood as we did in the preceding examples. Such a model would be inappropriate because it would assume that the data range from $-\infty$ to $+\infty$,  rather than from 0 to 1.  -->

## The generalized linear modeling framework

- A **link function** $g(\cdot)$  connects the linear model (real numbers ranging from $(-\infty,+\infty)$) to the quantity to be estimated (here, the probabilities $\theta_n$ in $[0,1]$). 
- A (common) link function in this case is the **logit link**: 

\begin{equation}
\eta_n = g(\theta_n) = \log\left(\frac{\theta_n}{1-\theta_n}\right)
\end{equation}


<!-- The term $\frac{\theta_n}{1-\theta_n}$ is called the **odds**.^[Odds are defined to be the ratio of the probability of success to the probability of failure. For example, the odds of obtaining a one in a fair six-sided die are $\frac{1/6}{1-1/6}=1/5$. The odds of obtaining a heads in a fair coin are 1/1.]  The logit link function is therefore a log-odds; it maps real numbers ranging from $(-\infty,+\infty)$ to probability values ranging from $[0,1]$. Figure \@ref(fig:logisticfun) shows the logit link function, $\eta = g(\theta)$, and the inverse logit, $\theta = g^{-1}(\eta)$, which is called the **logistic function**; the relevance of this logistic function will become clear in a moment. -->

----

```{r logisticfun, fig.cap = "The logit and inverse logit (logistic) function." ,echo = FALSE,warning=FALSE, fig.height=4.5}
x<-seq(0.001,0.999,by=0.001)
y<-log(x/(1-x))
logistic_dat<-data.frame(theta=x,eta=y)

p1<-qplot(logistic_dat$theta,logistic_dat$eta,geom="line")+xlab(expression(theta))+ylab(expression(eta))+ggtitle("The logit link")+
  annotate('text', 
           x = 0.3, y = 4, 
        label = expression(paste(eta,"=",g(theta))),parse = TRUE,size=8)
  

p2<-qplot(logistic_dat$eta,logistic_dat$theta,geom="line")+xlab(expression(eta))+ylab(expression(theta))+ggtitle("The inverse logit link (logistic)")+annotate('text', 
           x = -3.5, y = 0.80, 
        label = expression(paste(theta,"=",g^-1, eta)),parse = TRUE,size=8)

gridExtra::grid.arrange(p1,p2,ncol=2)

```


## Formal model
<!-- In summary, the generalized linear model with the logit link fits the following Bernoulli likelihood: -->

<!-- The model is fit on the log-odds scale, $\eta_n = \alpha + c\_set\_size_n \cdot \beta$.  -->
<!-- Once $\eta_n$ has been estimated, the inverse logit or the logistic function is used to compute the probability estimates  -->
<!-- $\theta_n =  \log(\frac{\exp(\eta_n)}{1+\exp(\eta_n)})$`.  An example of the calculations will be shown in the next section. -->

**The likelihood for each observation $n$:**

\begin{equation}
\eta_n = \log\left(\frac{\theta_n}{1-\theta_n}\right) = \alpha + \beta \cdot c\_set\_size
\end{equation}

\begin{equation}
\theta_n = g^{-1}(\eta_n) =  \log\left(\frac{\exp(\eta_n)}{1+\exp(\eta_n)}\right)
\end{equation}

\begin{equation}
correct_n \sim Bernoulli(\theta_n)
(\#eq:bernoullilogislik)
\end{equation}


## Priors for logistic regression

<!-- - $\alpha$ and $\beta$ probabilities or proportions, but *log-odds* -->

<!-- In order to decide on priors for $\alpha$ and $\beta$ we need to take into account that these parameter do not represent probabilities or proportions, but *log-odds*, the x-axis in Figure \@ref(fig:logisticfun) (right-hand side figure).  As shown in  the figure, the relationship between log-odds and probabilities is not linear.  -->

<!-- There are two functions in R that implement the logit and inverse logit functions: `qlogis(p)` for the logit function and `plogis(x)` for the inverse logit or logistic function. -->

<!-- Now we need to set priors for $\alpha$ and $\beta$.  -->

- $\alpha$ represents the *log-odds* of correctly recalling one word in a random position for the average set size of five (because we centered the predictor and since $5 = \frac{2+4+6+8}{4}$). (It's telling us  how difficult the task is. Let's assume (a 50/50 chance) with a great deal of uncertainty:
 <!-- which, incidentally, was not presented in the experiment. -->
 <!-- This is  one case where the intercept doesn't have a clear interpretation if we leave the prediction uncentered: With non-centered set size, the intercept will be the probability of recalling one word in a set of *zero* words.  -->

<!-- The prior for $\alpha$ will depend on how difficult the recall task is. If we are not sure, we could assume that the probability of recalling a word for an average set size, $\alpha$, is centered in .5 (a 50/50 chance) with a great deal of uncertainty. The `R` command `plogis(.5)` tells us that .5 corresponds to zero in log-odds. How do we include a great deal of uncertainty? We could look at Figure \@ref(fig:logisticfun), and decide on  a standard deviation of 4 in a normal distribution centered in zero: -->

We use `qlogis(p)` for the inverse logit or logistic function:

```{r}
qlogis(.5)
```

## Prior for $\alpha$

\begin{equation}
\alpha \sim Normal(0, 4) 
\end{equation}


<!-- Let's plot this prior in log-odds and in probability scale by drawing random samples. -->


(ref:logoddspriorsf) Prior for $\alpha \sim Normal(0, 4)$ in log-odds and in probability space.

```{r logoddspriorsf, fig.cap= "(ref:logoddspriorsf)",fig.show='hold', out.width = "45%", fig.width =3, fig.height =2, size = "scriptsize", tidy = FALSE}
samples_logodds <- tibble(alpha = rnorm(100000, 0, 4))
samples_prob <- tibble(p = plogis(rnorm(100000, 0, 4)))
ggplot(samples_logodds, aes(alpha)) + geom_density()
ggplot(samples_prob, aes(p)) + geom_density()
```

<!-- Figure \@ref(fig:logoddspriorsf) shows that our prior assigns more probability mass to extreme probabilities of recall than to intermediate values. Clearly, this is not what we intended. -->


## Prior for $\alpha$

<!-- We could try several values for standard deviation of the prior, until we find a prior that make sense for us. Reducing the standard deviation to 1.5 seems to make sense as shown in Figure \@ref(fig:logoddspriorsf2). -->

We try with:

\begin{equation}
\alpha \sim Normal(0, 1.5)
\end{equation}

  
(ref:logoddspriorsf2) Prior for $\alpha \sim Normal(0, 1.5)$ in log-odds and in probability space.

```{r logoddspriorsf2, fig.cap= "(ref:logoddspriorsf2)",fig.show='hold', echo=FALSE, out.width = "45%", fig.width =3, fig.height =2}
samples_logodds <- tibble(alpha = rnorm(100000, 0, 1.5))
samples_prob <- tibble(p = plogis(rnorm(100000, 0, 1.5)))
ggplot(samples_logodds, aes(alpha)) +
    geom_density()
ggplot(samples_prob, aes(p)) +
    geom_density()
```

## Prior for $\beta$

* $\beta$ represents the effect in log-odds of increasing the set size. 

(a) $\beta \sim Normal(0, 1)$
(b) $\beta \sim Normal(0, .5)$
(c) $\beta \sim Normal(0, .1)$
(d) $\beta \sim Normal(0, .01)$
(e) $\beta \sim Normal(0, .001)$

----

Edited version of the earlier  `normal_predictive_distribution_fast`:

```{r, size = "scriptsize", tidy = FALSE}
logistic_model_pred <- function(alpha_samples,
                                beta_samples,
                                set_size,
                                N_obs) {
    map2_dfr(alpha_samples, beta_samples,
             function(alpha, beta) {
                 tibble(set_size = set_size,
                                        # we center size:
                        c_set_size = set_size - mean(set_size),
                                        # change the likelihood: 
                                        # Notice the use of a link function for alpha and beta
                        theta = plogis(alpha + c_set_size * beta),
                        correct_pred = rbernoulli(N_obs,  p = theta))
             }, .id = "iter") %>%
    # .id is always a string and needs to be converted to a number
        mutate(iter = as.numeric(iter))
}
```

----

Let's assume 800 observations with 200 observation of each set size:
```{r, size = "scriptsize"}
N_obs <- 800
set_size <- rep(c(2,4,6,8),200)
```
We iterate over the four possible standard deviations of $\beta$:

```{r, size = "scriptsize", tidy = FALSE}
alpha_samples <- rnorm(1000, 0, 1.5)
sds_beta <- c(1, 0.5, 0.1,0.01, 0.001) 
prior_pred <- map_dfr(sds_beta, function(sd) {
    beta_samples <- rnorm(1000, 0, sd)
    logistic_model_pred(alpha_samples = alpha_samples,
                        beta_samples = beta_samples,
                        set_size = set_size,
                        N_obs = N_obs) %>%
        mutate(prior_beta_sd = sd)
})
```

----

And we calculate the accuracy  for each one of the priors we want to examine, for each iteration, and for each set size.

```{r, size = "scriptsize", tidy = FALSE}
(mean_accuracy <- prior_pred %>%
   group_by(prior_beta_sd, iter, set_size) %>% 
   summarize(accuracy = mean(correct_pred)) %>%
   mutate(prior = paste0("Normal(0, ",prior_beta_sd,")")))
```

-----

<!-- We plot it in Figure \@ref(fig:priors4beta), and as expected the priors are centered at zero. We see that the distribution of possible accuracies for the prior that has a standard deviation of one is problematic: There is too much probability mass concentrated near zero and one for set sizes of 2 and 8. -->


```{r priors4beta, message = FALSE, size = "scriptsize", fig.height =3.8}
mean_accuracy %>%
    ggplot(aes(accuracy)) +
    geom_histogram() +
    facet_grid(set_size~prior)
```

## Prior predicted differences in accuracy

<!-- It's hard to tell the differences between the other priors, and it might be more useful to look at the predicted differences in accuracy between set sizes. We calculate them as follows: -->

```{r, size = "scriptsize", tidy = FALSE}
(diff_accuracy <- mean_accuracy %>%
    arrange(set_size) %>%
    group_by(iter, prior_beta_sd) %>%
    mutate(diffaccuracy = accuracy - lag(accuracy) ) %>%
    mutate(diffsize = paste(set_size,"-",  lag(set_size))) %>%
    filter(set_size >2))
```

-----


```{r priors4beta2, message = FALSE, size = "scriptsize", fig.height =3.8}
diff_accuracy %>%
    ggplot(aes(diffaccuracy)) +
    geom_histogram() +
    facet_grid(diffsize~prior)
```


<!-- We plot them in Figure \@ref(fig:priors4beta2). If we are not sure whether the increase of set size could produce something between a null effect and a  relatively large effect, we can choose the prior with a standard deviation of $0.5$. Thus we settle on the following priors:  -->

-----

\begin{equation}
\begin{aligned}
\alpha &\sim Normal(0, 1.5) \\
\beta &\sim Normal(0, 0.1) 
\end{aligned}
\end{equation}


## The `brms` model

<!-- Having decided on the likelihood, the link function, and the priors, the model can now be fit using `brms`. Notice that we need to specify that the family is `bernoulli()`, and the link is `logit`.  -->


```{r,  message = FALSE }
fit_recall <- brm(correct ~ 1 + c_set_size,
  data = df_recall_data,
  family = bernoulli(link = logit),
  prior = c(
    prior(normal(0, 1.5), class = Intercept),
    prior(normal(0, .1), class = b, coef = c_set_size)
  )
)
```

------


<!-- Next, look at the summary of the posteriors of each of the parameters. Keep in mind that the parameters are in log-odds space: -->

```{r}
posterior_summary(fit_recall, pars = c("b_Intercept", "b_c_set_size"))
```

-----

<!-- Plot the posteriors as well: -->

```{r}
plot(fit_recall)
```

<!-- Next, we turn to the question of what we can report as our results, and what we can conclude from the data. -->

##  How to communicate the results?

```{r, echo=FALSE, results="hide"}
alpha_samples<- posterior_samples(fit_recall)$b_Intercept
beta_samples<- posterior_samples(fit_recall)$b_c_set_size
beta_mean <- round(mean(beta_samples),5)
beta_low <- round(quantile(beta_samples,prob=0.025),5)
beta_high <- round(quantile(beta_samples,prob=0.975),5)
```

If we want to talk about the effect estimated by the model in log-odds space, we summarize the posterior of $\beta$ in the following way: 

- $\hat\beta = `r mean(beta_samples)`$, 95% CrI = $[ `r beta_low` , `r beta_high` ]$.

## Effect in proportions

<!-- However, the effect might be easier to understand in proportions rather than in log-odds.   -->

Average accuracy for the task:

```{r, size = "scriptsize"}
alpha_samples<- posterior_samples(fit_recall)$b_Intercept
av_accuracy <- plogis(alpha_samples)
c(mean = mean(av_accuracy), quantile(av_accuracy, c(.025,.975)))
```

## Effect in proportions
 
Effect of our manipulation

  <!-- - the scale is not linear -->
  <!-- - the effect of increasing the set size depends on the average accuracy, and the set size that we start from. -->

```{r, size = "scriptsize"}
beta_samples<- posterior_samples(fit_recall)$b_c_set_size
effect_av_set_size <- plogis(alpha_samples) - plogis(alpha_samples - beta_samples)
c(mean = mean(effect_av_set_size), quantile(effect_av_set_size, c(.025,.975)))
```

Notice the interpretation here: 

*if we increase the set size from the average set size minus one to the average set size (5), we get a reduction in the accuracy of recall of $`r mean(effect_av_set_size)`$, 95% CrI = $[ `r quantile(effect_av_set_size, .025)` , `r quantile(effect_av_set_size, .975)` ]$.* 

## Effect in proportions

Recall that the average set size, 5, was not presented to the subject! 

Decrease in accuracy from a set size of 2 to 4:


```{r, size = "scriptsize"}
set4 <- 4 - mean(df_recall_data$set_size)
set2 <- 2 - mean(df_recall_data$set_size)
effect_4m2 <- plogis(alpha_samples + set4 * beta_samples) -
              plogis(alpha_samples + set2  * beta_samples)
c(mean = mean(effect_4m2), quantile(effect_4m2, c(.025,.975)))
```

We see that increasing the set size does have a detrimental effect in recall, as we suspected.

## Descriptive adequacy

We could also make predictions for other conditions not presented in the actual experiment, such as set sizes that weren't tested:

<!-- We could  then verify if our model was right with another experiment. To make predictions for other set sizes,  -->

- We extend our dataset adding rows with set sizes of 3, 5, and 7: we add 23 trials of each new set size

- Notice is that **we need to center our predictor based on the original mean set size**

<!-- To be consistent with the data of the other set sizes  in the experiment, we add 23 trials of each new set size (this is the number of trial by set sizes in the dataset). Something important to notice is that **we need to center our predictor based on the original mean set size**. This is because we want to maintain our interpretation of the intercept. We extend the data as follows, and we summarize the data and plot it  in Figure \@ref(fig:postpredsum2). -->

----

```{r, size = "scriptsize"}
df_recall_data_ext <- df_recall_data %>%
    bind_rows(tibble(set_size = rep(c(3,5,7),23),
                     c_set_size = set_size - mean(df_recall_data$set_size)))
df_recall_pred_ext <- posterior_predict(fit_recall,
                                 newdata = df_recall_data_ext,
                                 nsamples = 1000) %>%
    array_branch(margin = 1) %>%
    map_dfr( function(yrep_iter) {
        df_recall_data_ext %>%
            mutate(correct = yrep_iter)
    }, .id = "iter") %>%
    mutate(iter = as.numeric(iter))
```

----



```{r,  message = FALSE, size = "scriptsize", fig.height = 4}
(df_recall_pred_ext_summary <- df_recall_pred_ext %>%
    group_by(iter, set_size) %>%
    summarize(accuracy = mean(correct)))
```

----


```{r,  message = FALSE, size = "scriptsize", fig.height = 4}
# observed means:
(df_recall_summary<- df_recall_data %>%
    group_by(set_size) %>%
    summarize(accuracy = mean(correct)))
```

----

```{r postpredsum2,  message = FALSE, size = "scriptsize", fig.height = 4}
ggplot(df_recall_pred_ext_summary, aes(accuracy)) +
    geom_histogram(alpha=.5)+
    geom_vline(aes(xintercept= accuracy),data= df_recall_summary)+
    facet_grid(set_size ~ .)
```



<!-- ### \color{blue} Exercises -->

<!-- As with our log-normal dataset, our current dataset includes also a column that indicates the trial number. Could it be that trial has also an effect on the recall accuracy? How would you communicate the new results? How does recall changes from the first to the second trial? How does it change in the last two trials?  -->





## References
