---
title: "Bayesian hierarchical 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"]
header-includes:
  \setbeamerfont{caption}{size=\scriptsize}
---


# Bayesian hierarchical models  (also known as multilevel or mixed-effects models)

<!-- 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 = 55,
    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)

```




## The N400 effect (hierarchical normal likelihood)

In the EEG literature, it has been shown that words with low-predictability are accompanied by an *N400 effect* in comparison with high-predictable words, this is a relative negativity that peaks around 300-500 after word onset over central parietal scalp sites [first noticed in  @kutasReadingSenselessSentences1980;  for semantic anomalies and in @kutasBrainPotentialsReading1984 for low predictable word; for a review: @kutasThirtyYearsCounting2011].


1.  Example from @delongProbabilisticWordPreactivation2005 
    a. The day was breezy so the boy went outside to fly a kite.  
    b. The day was breezy so the boy went outside to fly an airplane. 

-----
 
```{r mot, fig.cap = "Typical ERP for the grand average across the N400 spatial window (central parietal electrodes: Cz, CP1, CP2, P3, Pz, P4, POz) for high and low predictability nouns (specifically from the constraining context of the experiment reported in Nicenboim, Vasishth, and Rösler 2020). The x-axis indicates time in seconds and the y-axis indicates voltage in microvolts (note that unlike many EEG/ERP plots, the negative polarity is plotted downwards). ", out.width = "80%",out.height="70%", echo = FALSE, fig.align = "center"}
knitr::include_graphics("cc_figure/N400noun-1.pdf", dpi =1000)
```

-----

* We simplify the high-dimensional EEG data by focusing on the average amplitude of the EEG signal at the typical spatio-temporal window of the N400.
* We  focus on the N400 effect for nouns from a subset of the data from @nieuwlandLargescaleReplicationStudy2018. (To speed-up computation,  we'll restrict the dataset to the participants from the Edinburgh lab)

-----

```{r, message = FALSE, size = "scriptsize"}
df_eeg_data <- read_tsv("data/public_noun_data.txt") %>%
    filter(lab=="edin") %>%
    mutate(c_cloze= cloze/100 - mean(cloze/100) )
df_eeg_data$c_cloze %>% summary() 
```

-----

One nice aspect of this dataset is that the dependent variable is roughly normally distributed:


```{r, fig.cap="Histogram of the N400 averages for every trial in gray; density plot of a normal distribution in red.", message=FALSE, echo = FALSE, fig.height = 4.5}
df_eeg_data %>% ggplot(aes(n400)) +
    geom_histogram(binwidth = 4, colour="gray", alpha = .5, aes(y = ..density..)) +
    stat_function(fun = dnorm, color ="red", args = list(
                                                 mean = mean(df_eeg_data$n400),
                                                 sd = sd(df_eeg_data$n400))) +
    xlab("Average voltage in microvolts for the N400 spatiotemporal window")
```

## A complete pooling model 

We'll start from the simplest model which is basically a linear regression. **Note that this model is incorrect for these data due to point 2 below.**

* Model $M_{cp}$ assumptions:

1. EEG averages for the N400 spatiotemporal window are normally distributed.
2. Observations are *independent*.
3. There is a linear relationship between cloze and the EEG average for the trial.

-----

* Likelihood:

 \begin{equation}
   signal_n \sim Normal(\alpha + c\_cloze_n \cdot \beta,\sigma)
   \end{equation}
 

<!-- Predictive priors here? -->

* Priors:

 \begin{equation}
 \begin{aligned}
 \alpha &\sim Normal(0,10)\\
 \beta  &\sim Normal(0,10)\\
 \sigma  &\sim Normal_{+}(0,50)
 \end{aligned}
 \end{equation}


## Fitting the model

```{r, message = FALSE, results = "hide", tidy = FALSE}
fit_N400_cp <- brm(n400 ~ c_cloze,
  prior =
    c(prior(normal(0, 10), class = Intercept),
      prior(normal(0, 10), class = b),
      prior(normal(0, 50), class = sigma)),
  data = df_eeg_data
)
```


-------

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

----


```{r, fig.height = 6 }
plot(fit_N400_cp)
```

## No pooling model 

* Model $M_{np}$ assumptions:

1. EEG averages for the N400 spatio-temporal window are normally distributed.
2. Observations depend *completely* on the participant. (Participants have nothing in common.)
3. There is a linear relationship between cloze and the EEG average for the trial.

-----


* Likelihood:

 \begin{equation}
 signal_n \sim Normal( \alpha_{i[n]} + c\_cloze_n \cdot \beta_{i[n]},\sigma)
 \end{equation}
 

* Priors:

 \begin{equation}
 \begin{aligned}
 \alpha_i &\sim Normal(0,10)\\
 \beta_i  &\sim Normal(0,10)\\
 \sigma  &\sim Normal_+(0,50)
 \end{aligned}
 \end{equation}

----


We fit it in brms by removing the common intercept with `0 +` and thus having an intercept and effect for each level of `subject`:

```{r, message = FALSE, results = "hide", tidy = FALSE, size ="scriptsize"}
fit_N400_np <- brm(n400 ~ 0 +
                     factor(subject) + c_cloze:factor(subject),
                 prior =
                     c(prior(normal(0, 10), class = b),
                       prior(normal(0, 50), class = sigma)),
                 data = df_eeg_data)
```

-----

\vspace{.1in}
```{r, output.lines=12, size = "tiny"}
fit_N400_np
```

-----

We plot the estimates using `bayesplot`.

\bigskip

```{r, eval = FALSE, size ="scriptsize"}
# I first peek at the internal names of the parameters.
# parnames(fit_N400_np)
ind_effects_np <- paste0("b_factorsubject",
                         unique(df_eeg_data$subject), ":c_cloze")
mcmc_intervals(fit_N400_np, pars=ind_effects_np,
               prob = 0.8,
               prob_outer = 0.95,  
               point_est = "mean"
)
```

--------


```{r, echo = FALSE}
# I first peek at the internal names of the parameters.
# parnames(fit_N400_np)
ind_effects_np <- paste0("b_factorsubject",unique(df_eeg_data$subject), ":c_cloze")
mcmc_intervals(fit_N400_np, pars=ind_effects_np,
               prob = 0.8,
               prob_outer = 0.95,  
               point_est = "mean"
)
```

--------


We can then calculate the  average of the $\beta$'s, even though the model doesn't assume that there's one common $\beta$:

```{r, tidy = FALSE, size= "footnotesize"}
average_beta_across_subj <-
    posterior_samples(fit_N400_np,
                      pars = ind_effects_np) %>%
    rowMeans()
c(mean=mean(average_beta_across_subj),
  quantile(average_beta_across_subj,
           c(.025,.975)))
```

## Varying intercept and varying slopes model ($M_{v}$)


* Model $M_{v}$ assumptions:

1. EEG averages for the N400 spatio-temporal window are normally distributed.
3. Each subject deviates to some extent (this is made precise below) from the grand mean and from the mean effect of predictability.
4. There is a linear relationship between cloze and the EEG average for the trial.


----

* Likelihood:

\begin{equation}
 signal_n \sim Normal(\alpha + u_{0,i[n]} + c\_cloze_n \cdot (\beta+ u_{1,i[n]}),\sigma)
 \end{equation}
 
* Prior:

\begin{equation}
 \begin{aligned}
 \alpha &\sim Normal(0,10)\\
 \beta  &\sim Normal(0,10)\\
 u_0 &\sim Normal(0,\tau_{u_0})\\
 u_1 &\sim Normal(0,\tau_{u_1})\\
 \tau_{u_0} &\sim Normal_+(0,20) \\
 \tau_{u_1} &\sim Normal_+(0,20) \\
 \sigma  &\sim Normal_+(0,50)
 \end{aligned}
 \end{equation}


----------

Some important (and sometimes confusing) points:

* Why does $u$ have a mean of 0? 

Because we want $u$ to capture only differences between subjects, we could achieve the same by assuming that

\begin{equation}
\begin{aligned}
\mu_n &= \alpha_{i[n]} + \beta_{i[n]} \cdot c\_cloze_n \text{and} \\
\alpha_i &\sim Normal(\alpha,\tau_{u_0})\\
\alpha &\sim Normal(0,10)\\
\beta_i &\sim Normal(\beta,\tau_{u_1})\\
\beta &\sim Normal(0,10)
\end{aligned}
\end{equation}

And in fact, that's another common way to write the model.

-----

* Why do the adjustments $u$ have a normal distribution?

Mostly because of "convention", that's the way it's implemented in most frequentist mixed models. 

But also because if we don't know anything about the distribution besides its mean and variance, the normal distribution is the most conservative assumption [see also chapter 9 of @mcelreath2015statistical].
 

--------

Let's see how we need to set up the priors:

```{r, size = "scriptsize"}
get_prior(n400 ~ c_cloze + (c_cloze || subject), data=df_eeg_data)
```

------

```{r, message = FALSE, results = "hide", tidy=FALSE, size = "scriptsize"}
fit_N400_v <- brm(n400 ~ c_cloze + (c_cloze || subject),
               prior =
                 c(prior(normal(0, 10), class = Intercept),
                   prior(normal(0, 10), class = b, coef = c_cloze),
                   prior(normal(0, 50), class = sigma),
                   prior(normal(0, 20), class = sd, coef = Intercept,
                         group = subject),
                   prior(normal(0, 20), class = sd, coef = c_cloze,
                         group = subject)
                   ),
               data = df_eeg_data)
```

-----------

\vspace{.1in}
```{r,size="scriptsize"}
fit_N400_v
```

-----------

```{r, fig.height = 6 }
plot(fit_N400_v, N=6) 
```

## Individual effects


```{r, eval = FALSE, size ="scriptsize"}
#parnames(m_N400_v)
ind_effects_v <- paste0("r_subject[",unique(eeg_data$subject), ",c_cloze]")
mcmc_intervals(fit_N400_v, pars=ind_effects_v,
               prob = 0.8,
               prob_outer = 0.95,  
               point_est = "mean"
) 
```

-----

\vspace{.1in}
```{r, echo = FALSE, fig.height = 8}
#parnames(fit_N400_v)
ind_effects_v <- paste0("r_subject[",unique(df_eeg_data$subject), ",c_cloze]")
mcmc_intervals(fit_N400_v, pars=ind_effects_v,
               prob = 0.8,
               prob_outer = 0.95,  
               point_est = "mean"
) 
```



## Shrinkage

```{r comparison, message=F, fig.height=5,fig.width=6, echo = FALSE, out.width ="200%"}
# We'll need to make the plot "manually"
 
# No pooling model
par_np <- posterior_summary(fit_N400_np)[ind_effects_np,] %>%
    as_tibble() %>%
    mutate(model = "No pooling",
           subj = unique(df_eeg_data$subject))
# For the hierarchical model is more complicated, because we want the effect (beta) + adjustment:

par_h <- posterior_samples(fit_N400_v) %>%
    select(ind_effects_v)  %>%
    mutate_all( ~ . + posterior_samples(fit_N400_v)$b_c_cloze) %>%
    map_dfr(~ tibble(Estimate = mean(.),
                Q2.5 = quantile(.,.025),
                Q97.5 = quantile(., .975))) %>%
    mutate(model = "Hierarchical",
           subj = unique(df_eeg_data$subject))

by_subj_df <- bind_rows(par_h, par_np)


ggplot(by_subj_df, aes(ymin = Q2.5, ymax = Q97.5,x=subj, y = Estimate, color = model)) +
     geom_errorbar(position = position_dodge(1)) +
    geom_point(position = position_dodge(1)) +

# We'll also add the mean and 95% CrI of the overall difference to the plot:
    geom_hline(yintercept = posterior_summary(fit_N400_v)["b_c_cloze","Estimate"], linetype = "dotted") +
    geom_hline(yintercept = posterior_summary(fit_N400_v)["b_c_cloze","Q2.5"], linetype = "dotted",size = .5)+
    geom_hline(yintercept = posterior_summary(fit_N400_v)["b_c_cloze","Q97.5"], linetype = "dotted",size = .5) +
    xlab("N400 effect of predictability") +
    coord_flip()
```

## Correlated varying intercept varying slopes model ($M_{h}$) {#sec:mcvivs}

* In $M_h$, we  model the EEG data with the following assumptions:

1. EEG averages for the N400 spatio-temporal window are normally distributed.
2. Some aspects of the signal voltage and the effect of predictability on the signal depend on the participant, and these two might be correlated, i.e., we assume random intercept, slope and correlation by-subject.
3. There is a linear relationship between cloze and the EEG average for the trial.

--------

* Likelihood:

 \begin{equation}
  signal_n \sim Normal(\alpha + u_{i[n],0} + c\_cloze_n \cdot  (\beta + u_{i[n],1}),\sigma)
  \end{equation}
 

We need to have priors on the adjustments for intercept and slopes, $u_{,0-1}$. 

* Priors:
 \begin{equation}
 \begin{aligned}
   \alpha & \sim Normal(0,10) \\
   \beta  & \sim Normal(0,10) \\
    \sigma  &\sim Normal_+(0,50)\\
    {\begin{pmatrix}
    u_{i,0} \\
    u_{i,1}
    \end{pmatrix}}
   &\sim {\mathcal {N}}
    \left(
   {\begin{pmatrix} 
    0\\
    0
   \end{pmatrix}}
 ,\boldsymbol{\Sigma_u} \right)
 \end{aligned}
 \end{equation}

 
-----------

\begin{equation}
\boldsymbol{\Sigma_u} = 
{\begin{pmatrix} 
\tau_{u_0}^2 & \rho_u \tau_{u_0} \tau_{u_1} \\ 
\rho_u \tau_{u_0} \tau_{u_1} & \tau_{u_1}^2
\end{pmatrix}}
\end{equation}

<!-- In addition, it has a correlation matrix associated with it:^[Since we have two random variables, there is only one correlation in the matrix between $u_1$ and $u_2$. With 3 random variables (i.e., adding $u_3$), we end up with a $3 \times 3$ correlation matrix that includes three different correlations $\rho_{u_{1,2}}=\rho_{u_{2,1}}$, $\rho_{u_{1,3}}=\rho_{u_{3,1}}$, and $\rho_{u_{2,3}}=\rho_{u_{3,2}}$.] -->

<!-- \begin{equation} -->
<!-- {\begin{pmatrix}  -->
<!-- 1 & \rho_u  \\  -->
<!-- \rho_u  & 1 -->
<!-- \end{pmatrix}} -->
<!-- \end{equation} -->

<!-- ----- -->

<!-- We still need to define a prior for $\boldsymbol{\Sigma_u}$. We can decompose the covariance matrix  in the following -->
<!-- way: -->

<!-- \begin{equation} -->
<!-- \begin{aligned} -->
<!-- \boldsymbol{\Sigma_u} &= diag\_matrix(\tau_u) \cdot \boldsymbol{\rho_u} \cdot diag\_matrix(\tau_u)\\ -->
<!-- &= -->
<!-- {\begin{pmatrix}  -->
<!-- \tau_{u_1} & 0 \\  -->
<!-- 0  & \tau_{u_2} -->
<!-- \end{pmatrix}} -->
<!-- {\begin{pmatrix}  -->
<!-- 1 & \rho_u  \\  -->
<!-- \rho_u  & 1 -->
<!-- \end{pmatrix}} -->
<!-- {\begin{pmatrix}  -->
<!-- \tau_{u_1} & 0 \\  -->
<!-- 0  & \tau_{u_2} -->
<!-- \end{pmatrix}} -->
<!-- \end{aligned} -->
<!-- \end{equation} -->

<!-- ------------- -->

---

And now we need priors for the $\tau_u$s and for $\rho_u$:

\begin{equation}
\begin{aligned}
\tau_{u_0} &\sim Normal_+(0,20)\\
\tau_{u_1} &\sim Normal_+(0,20)\\
\rho_u &\sim LKJcorr(2) 
\end{aligned}
\end{equation}


------

\vspace{.1in}

(ref:lkjviz) Visualization of the LKJ prior with four different values of the $\eta$ parameter.

```{r lkjviz,echo=FALSE, fig.cap ="(ref:lkjviz)", message= FALSE,warning=FALSE,results="asis",fig.width =4, fig.height=2,fig.show='hold', out.width='.48\\linewidth'}

## https://github.com/rmcelreath/rethinking/blob/1def057174071beb212532d545bc2d8c559760a2/R/distributions.r
# onion method correlation matrix
dlkjcorr <- function( x , eta=1 , log=FALSE ) {
    ll <- det(x)^(eta-1)
    if ( log==FALSE ) ll <- exp(ll)
    return(ll)
}

dlkjcorr2 <- function(rho, eta = 1 ) {
    purrr::map_dbl(rho, ~ matrix(c(1, .x,.x,1),ncol=2) %>%  
                                 dlkjcorr(., eta))
}

ggplot(tibble(rho = c(-.99,.99)), aes(rho)) +
    stat_function(fun = dlkjcorr2,  geom = "line", args = list(eta = 1)) +
    ylab("density") +
    ggtitle("eta = 1")

ggplot(tibble(rho = c(-.99,.99)), aes(rho)) +
    stat_function(fun = dlkjcorr2,  geom = "line", args = list(eta = 2)) +
    ylab("density") +
    ggtitle("eta = 2")

ggplot(tibble(rho = c(-.99,.99)), aes(rho)) +
    stat_function(fun = dlkjcorr2,  geom = "line", args = list(eta = 4)) +
    ylab("density") +
    ggtitle("eta = 4")

ggplot(tibble(rho = c(-.99,.99)), aes(rho)) +
    stat_function(fun = dlkjcorr2,  geom = "line", args = list(eta = .9)) +
    ylab("density") +
    ggtitle("eta = .9")


```




-----

Let's see how we need to set up the priors:

```{r, size ="scriptsize"}
get_prior(n400 ~ c_cloze + (c_cloze | subject), data=df_eeg_data)
```


## Fitting the model

```{r, message = FALSE, results = "hide", tidy = FALSE, size = "scriptsize"}
fit_N400_h <- brm(n400 ~ c_cloze + (c_cloze | subject),
                  prior =
                      c(prior(normal(0, 10), class = Intercept),
                        prior(normal(0, 10), class = b, coef = c_cloze),
                        prior(normal(0, 50), class = sigma),
                        prior(normal(0, 20), class = sd, coef = Intercept,
                              group = subject),
                        prior(normal(0, 20), class = sd, coef = c_cloze,
                              group = subject),
                        prior(lkj(2), class = cor,
                              group= subject)),
              data = df_eeg_data)
```

----

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

----


```{r, fig.height = 7}
plot(fit_N400_h, N=6)
```



<!-- TODO -->

## Why should we take the trouble of fitting a Bayesian hierarchical model?

* We can better characterize the generative process by adding the relevant clusters in our data (participants, items, maybe labs, etc)
* The same approach we used here can be used to extend any parameter of any model:
   * (generalized) linear models
   * non-linear/cognitive models

## How much structure should we add to our statistical models?

### The level of complexity depends on 

1. the answers we are looking for 
2. the size of the data at hand  
3. our computing power
4. our domain and experimental knowledge. 

\vspace{.3in}

> "Simplification is essential, but it comes at a cost, and real understanding depends in part on understanding the effects of the simplification" @mcclellandPlaceModelingCognitive2009


<!-- ### Why should we take the trouble of fitting a Bayesian hierarchical model? -->

<!-- Carrying out Bayesian data analysis clearly requires much more effort than fitting a frequentist model: we have to define priors, verify that our model works, and decide how to interpret the results. By comparison, fitting a linear mixed model using `lme4` consists of only a single line of code; the model fit using `lmer` makes many assumptions, but they are hidden from the user. We want to emphasize that there are important motivations for fitting Bayesian hierarchical models. -->

<!--  and also non-linear models,
including highly complex cognitive models [For an accessible introduction of Bayesian methods for cognitive modeling see @LeeWagenmarks2014].
.
 -->

<!-- * The same approach we used here can be used to extend any parameter of any  -->
<!--   model. This includes popular uses, such as logistic and Poisson regressions, and also useful models that are relatively rarely used in cognitive science such as multi-logistic regression (e.g., accuracy in some task with more than two answers), ordered logistic (e.g., ratings), and  models with a shifted log-normal distribution [see @NicenboimEtAl2016Frontiersb; @Rouder2005]. We provide examples of this flexibility in the coming chapters. -->
 
<!-- * Complex cognitive models can be extended hierarchically in a straightforward way, see @Lee2011 and @LeeWagenmakers2014. This is because, as we have seen with distributional regression models in section \@ref({#sec:distrmodel}), any parameter can have a group-level effect strucure. Some examples of hierarchical computational cognitive models in psycholinguistics are @LogacevVasishth2015, @nicenboimModelsRetrievalSentence2018, @VasishthEtAl2017Modelling, and @VasishthEtAl2017Feature. -->


## References
