---
title: "Computational Bayesian data analysis"
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"]
---





<!-- 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)
```

* Deriving  the posterior distribution analytically  is possible for only a very limited number of cases. 
* The denominator, the marginal likelihood, requires us to compute a possibly intractable integral:


\begin{equation}
p(\Theta|y) = \cfrac{ p(y|\Theta) \cdot p(\Theta) }{\int_{\Theta} p(y|\Theta) \cdot p(\Theta) d\Theta }
\end{equation}

## Alternative: Deriving the posterior through sampling



```{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)

```

```{r, echo = FALSE, message=FALSE}
lst_cloze_data <- list(k = 80, N = 100)
# Fit the model with the default values of number of chains and iterations
fit_cloze <- stan(
  file = "stan_models/binomial_cloze.stan", data = lst_cloze_data,
  warmup = 1000,
  iter = 21000
)
```

```{r, echo = FALSE}
samples_theta <- rstan::extract(fit_cloze)$theta
some_samples <- toString(round(head(samples_theta, 20), 3))
df_fit_cloze <- data.frame(theta = samples_theta)
diff_means <- format(mean(samples_theta) - 84 / (84 + 24), digits =1)
diff_var <- format(sd(samples_theta)^2 - 84 * 24 / ((84 + 24)^2 * (84 + 24 + 1)), digits=1)
```

### We want to derive  the posterior distribution of the Cloze probability of *"umbrella"*, $\theta$:
* Data: a word (e.g., *"umbrella"*) was answered 80 out of 100 times,
* Likelihood: a  binomial distribution 
* Prior for $\theta$:  $Beta(a=4,b=4)$

---

### We sample from the posterior distribution of $\theta$:
* We use a probabilistic programming language,
* given enough samples we will have a good approximation of the real posterior distribution,
* say we got 20000 samples from the posterior distribution of the Cloze probability, $\theta$:

\scriptsize

`r some_samples`, ... 

---

The approximation of the posterior looks quite similar to the real posterior.^[The difference between the true and the approximated mean and variance are `r diff_means` and `r diff_var` respectively]


(ref:betapost) Histogram of the samples of $\theta$ from the posterior distribution calculated through sampling in gray;  density plot of the exact posterior in red.

```{r betapost, fig.cap="(ref:betapost)",echo = FALSE, fig.height =3.5}
ggplot(df_fit_cloze, aes(theta)) +
  geom_histogram(binwidth = .01, colour = "gray", alpha = .5, aes(y = ..density..)) +
  stat_function(fun = dbeta, color = "red", args = list(
    shape1 = 84,
    shape2 = 24
  ))
```



## Computational Bayesian data analysis:

### Why has Bayes suddenly gained in prominence since the 2000s?

* increase in computing power 
* appearance of probabilistic programming languages: WinBUGS [@lunn2000winbugs], JAGS [@plummer2016jags], and more recently pymc3 [@Salvatier2016], Stan [@carpenter2017stan], INLA [@blangiardo2015spatial]. 

<!-- These statistical languages allow the user to define models without having to deal (for the most part) with the complexities of the sampling process. However, they require learning a new language since they the user has to fully specify the statistical model using a particular syntax.^[The python package pymc3 is an exception since it is fully integrated into python.] Furthermore, some knowledge of the sampling process is needed to correctly parameterize the models and to avoid convergence issues (these topics will be covered in detail later in this book).  -->

### Easier alternatives based on Stan:

* `rstanarm` [@rstanarm] 
* `brms` [@R-brms] 

# Bayesian Regression Models using 'Stan': brms

## Load the following:

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

```

# Examples 1:  A single participant pressing a button repeatedly (A simple linear model)

----

We have data from a participant repeatedly pressing the space bar as fast as possible, without paying attention to any stimuli. 

#### Data: 
reaction times in milliseconds  in each trial

### Question:
How long does it take to press a key when there is no decision involved?

-----

### Assumptions:

1. There is a true underlying time, $\mu$, that the participant needs to press the space bar.
2. There is some noise in this process.
3. The noise is normally distributed (this assumption is questionable given that reaction times are generally  skewed; we fix this assumption later).

## Formal model:

###  Likelihood for each observation $n$:

\begin{equation}
\begin{aligned}
rt_n \sim Normal(\mu, \sigma)
\end{aligned}
(\#eq:rtlik)
\end{equation}

<!-- where $n =1 \ldots N$, and $rt$ is the dependent variable (reaction times in milliseconds). The variable $N$ indexes the total number of data points. The letter $\mu$ indicates the *location* of the normal distribution function, the location  parameter shifts the distribution left or right on the horizontal axis. For the normal distribution, the location is also the mean of the distribution. The letter $\sigma$ indicates the *scale* of the distribution, as the scale decreases, the distribution gets narrower. The compressing approaches a spike (all the probability mass in one point) as the scale parameter goes to zero. For the normal distribution, the scale is also its standard deviation. -->

<!-- For a frequentist model that will give us the maximum likelihood estimate (the sample mean) of the time it takes to press the space bar, this would be enough information to write the formula in `R`, `rt ~ 1`, and plug it into the function `lm()` together with the data: `lm(rt ~ 1, data)`. The meaning of the `1` here is that there is no predictor associated with this parameter, and `lm` will estimate the so-called intercept of the model, in our case $\mu$. -->



<!-- For a Bayesian model, we will also need to define priors for the two parameters of our model. Let's say that we know for sure that the time it takes to press a key will be positive and lower than a minute (60000ms), but we don't want to make a commitment regarding which values are more likely. We encode what we know about the noise in the task in $\sigma$: we know that this parameter  must be positive and we'll assume that any value below 2000ms is equally likely. These priors are in general strongly discouraged because even when we know very little, a flat (or very wide) prior will almost never be the best approximation of what we know. We'll use them in this section for pedagogical purposes; the next chapter will show more realistic uses of priors. -->

###  (Bad) priors:

\begin{equation}
\begin{aligned}
\mu &\sim Uniform(0, 60000) \\
\sigma &\sim Uniform(0, 2000) 
\end{aligned}
(\#eq:rtpriors)
\end{equation}



## Fitting the model


We'll first load the data from `data/button_press.csv`: 

```{r, reading_noreading, message = FALSE}
df_noreading_data <-
    read_csv("./data/button_press.csv")
df_noreading_data
```

------------

```{r m1visualize, fig.cap="Visualizing the data", fig.height =5}
ggplot(df_noreading_data, aes(rt)) +
  geom_density() +
  ggtitle("Button-press data")
```

## Specifying the model in `brms`


```{r, message = FALSE, cache = TRUE}
fit_press <- brm(rt ~ 1,
  data = df_noreading_data,
  family = gaussian(),
  prior = c(prior(uniform(0, 60000), class = Intercept),
      prior(uniform(0, 2000), class = sigma)),
  chains = 4,
  iter = 2000,
  warmup = 1000
)
```

## Sampling and convergence in a nutshell

<!-- We start four chains independent from each other. Each chain "searches" for samples of the posterior in a multidimensional space, where each parameter corresponds to a dimension, and the shape of this space is determined by the priors and the likelihood. The chains start in random locations and in each iteration they take one sample each. The samples at the beginning do not belong to the posterior distribution. Eventually, the chains end up in the vicinity of the posterior distribution, and from that point onwards the samples will belong to the posterior. That means that at the beginning the samples from the different chains will be far from each other, but that *at some point* they will converge. While there are no guarantees that we are running the chains for enough iterations, the default values of `brms` (and Stan) are in many cases enough to achieve that, and when they are not, we will receive warnings with recommendations. If the chains converged to the same distribution, by removing the "warmup" (also called burn-in) samples--by default half of a total of 2000 iterations--, we make sure that we do not get samples from the path to the posterior distribution; see figure \@ref(fig:warmup). Stan runs diagnostics with the information from the chains, and if there are no warnings after fitting the model, we can be reasonable sure that the model converged and our samples are from the true posterior distribution. However, we do need to run more than one chain (preferably four), with a couple of thousands of iterations (at least) so that the diagnostics will work. -->

\scriptsize

1. Chains start in random locations;
2. in each iteration they take one sample each;
3. samples at the beginning do not belong to the posterior distribution;
4. eventually, the chains end up in the vicinity of the posterior distribution;
5. from that point onwards the samples will belong to the posterior. 


(ref:warmup) Trace plot of the `brms` model

```{r warmup, fig.cap = "(ref:warmup)", echo = FALSE, fig.height =3.5 }
library(bayesplot)
plot_all <- function(fit, xlab=500, ylab =100){
  chains <- as.mcmc(fit, inc_warmup = TRUE) %>%
    map_dfr(~ as.data.frame(.x) %>%
              as_tibble() %>%
              mutate(iter = 1:n()), .id = "chain") %>%
    rename(mu = b_Intercept) %>%
    dplyr::select(-`lp__`) %>%
    tidyr::gather("parameter", "value", -iter, -chain)

  ggplot(chains, aes(x = iter, y = value, color = chain)) + geom_line() +
    facet_wrap(~parameter, ncol = 1) +
    geom_vline(xintercept = 1000, linetype = "dashed") +
    xlab("Iteration number") +
    ylab("Sample value")+
    annotate("text", x=500, y =100, color= "black", label="Warm-up", size =5.2)
}

plot_all(fit_press)
```

----
 
(ref:warmup2) Trace plot of a  model that **did not** converge.

```{r warmup2, fig.cap = "(ref:warmup2)", echo = FALSE, warning = FALSE, message = FALSE}
data_mm <- tibble(rt = rnorm(500, c(5,3000), c(5,5)))
fit_press_bad <- brm(rt ~ 1,
  data = data_mm,
  family = gaussian(),
  prior = c(
    prior(uniform(0, 60000), class = Intercept),
    prior(uniform(0, 2000), class = sigma)
  ),
  chains = 4,
  iter = 2000,
  warmup = 1000
  )

plot_all(fit_press_bad, ylab = 500)
```

## Output of `brms`

<!-- The model converged (i.e., if we didn't have any warning messages), the output of the sampling process shows the samples of the posterior distributions of each of the parameters:  -->

```{r}
posterior_samples(fit_press) %>% str() 
```

<!-- Notice that `b_Intercept` corresponds to our $\mu$ and that `lp` is not really part of the posterior, it's the density of the unnormalized posterior for each iteration. -->

## Output of `brms`
```{r, fig.height = 4}
plot(fit_press)
```


## Output of `brms`
```{r, results = "hold", size = "scriptsize"}
fit_press
# posterior_summary(fit_press) is also useful
```

## Output of `brms`
Notice that the `Estimate` is just the mean of the posterior sample, and
CI are the 95% quantiles:

```{r}
posterior_samples(fit_press)$b_Intercept %>%
                               mean()
posterior_samples(fit_press)$b_Intercept %>%
                               quantile(c(0.025, .975))
```

---



## Important questions
\large
<!-- We see that we can fit our model without problems, and we get some posterior distributions for our parameters. However, we should ask ourselves the following questions: -->

1. What information are the priors encoding? Do the priors make sense?
2. Does the likelihood assumed in the model make sense for the data?

<!-- We'll try to answer these questions by looking at the *Prior and posterior predictive distributions*, and by doing sensitivity analyses as described in the following sections. -->


#  Prior predictive distributions {#sec:priorpred}

We had defined the following priors for our linear model:

\begin{equation}
\begin{aligned}
\mu &\sim Uniform(0, 60000) \\
\sigma &\sim Uniform(0, 2000) 
\end{aligned}
(\#eq:rtpriorsrepeated)
\end{equation}

These priors encode assumptions about the kind of data we would expect to see in a future study. 

<!-- To understand these assumptions, we are going to generate data from the model; such data, which is generated entirely by the prior distributions, is  called the prior predictive distribution. Generating prior predictive distributions repeatedly helps us to check whether the priors make sense. What we want to know here is,  -->

Do the priors generate realistic-looking data? 


##  Prior predictive distributions 

We want to know the density $p(\cdot)$ of a data vector $y=\langle y_1,\dots,y_n\rangle$, given a vector of priors $\Theta$ (e.g., $\Theta=\langle\mu,\sigma \rangle$) 

The prior predictive density is:

\begin{equation}
p(y)= \int p(y|\Theta) p(\Theta) \, d\Theta 
\end{equation}

We avoid doing the integration analytically; we will use Monte Carlo integration. We repeat the following:

1.  Take one sample from each of the priors.
2.  Plug those samples in the likelihood and generate a dataset $y_{pred}=\langle y_{pred_1},\ldots,y_{pred_n}\rangle$. 


---

\vspace{.1in}

```{r,size = "scriptsize"}
normal_predictive_distribution <- function(mu_samples, sigma_samples, N_obs) {
  # empty data frame with headers:
  df_pred <- tibble(trialn = numeric(0),
                    rt_pred = numeric(0),
                    iter = numeric(0))
  # i iterates from 1 to the length of mu_samples,
  # which we assume is identical to 
  # the length of the sigma_samples:
  for (i in seq_along(mu_samples)) {
    mu <- mu_samples[i]
    sigma <- sigma_samples[i]
    df_pred <- bind_rows(
      df_pred,
      tibble(trialn = seq_len(N_obs), #1, 2,... N_obs
        rt_pred = rnorm(N_obs, mu, sigma),
        iter = i)
    )
  }
  df_pred
}
```

---

This approach works, but it's quite slow:

```{r, results = "hold",size = "scriptsize", tidy = FALSE}
tic()
N_samples <- 1000
N_obs <- nrow(df_noreading_data)
mu_samples <- runif(N_samples, 0, 60000)
sigma_samples <- runif(N_samples, 0, 2000)
normal_predictive_distribution(mu_samples = mu_samples,
                               sigma_samples = sigma_samples,
                                N_obs = N_obs)
toc()
```

---

A more efficient version:

```{r, results = "hold",size = "scriptsize", tidy = FALSE}
normal_predictive_distribution_fast <- function(mu_samples,
                                                sigma_samples,
                                                N_obs) {
  # map_dfr works similarly to lapply, it essentially runs
  # a for-loop, and builds a dataframe with the output.
  # We iterate over the values of mu_samples and sigma_samples
  # simultaneously, and in each iteration we bind a new
  # data frame with N_obs observations.
  map2_dfr(mu_samples, sigma_samples, function(mu, sigma) {
    tibble(
      trialn = seq_len(N_obs),
      rt_pred = rnorm(N_obs, mu, sigma)
    )}, .id = "iter") %>%
    # .id is always a string and needs to be converted to a number
    mutate(iter = as.numeric(iter))
}
```


----
```{r, results = "hold",size = "scriptsize", tidy = FALSE}
tic()
(prior_pred <- normal_predictive_distribution_fast(
  mu_samples = mu_samples,
  sigma_samples = sigma_samples,
  N_obs))
toc()
```

----

\vspace{.1in}
<!-- Figure \@ref(fig:priorpred-simple) shows the first 18 samples of the prior predictive distribution. These are 18 predicted datasets.  -->

(ref:priorpred-simple) Eighteen samples from the prior predictive distribution.

```{r priorpred-simple, fig.cap = "(ref:priorpred-simple)", message = FALSE, fig.height = 4, size ="scriptsize" }
prior_pred %>%
  filter(iter <= 12) %>%
  ggplot(aes(rt_pred)) +
  geom_histogram() +
  facet_wrap(~iter, ncol = 3)
```

<!-- The prior predictive distribution in Figure \@ref(fig:priorpred-simple)  shows prior datasets that are not realistic: Besides the fact that the datasets show that reaction times distributions are symmetrical--and we know that they are generally right-skewed--, some datasets present reaction times that are unrealistically long, and worst yet, dataset 18 presents negative press time values. -->


## Distribution of statistics
\vspace{.1in}
<!-- We can also look at the distribution of statistics here. Even if we don't know beforehand what the data should look like, it's very likely that we have some expectations for possible mean, minimum, or maximum values. -->

```{r, size = "scriptsize", tidy = FALSE}
(prior_stat <- prior_pred %>%
  group_by(iter) %>%
  summarize(min_rt = min(rt_pred),
            max_rt = max(rt_pred),
            average_rt = mean(rt_pred)) %>%
  # we convert the previous data frame to a long one,
  # where min_rt, max_rt, average_rt are possible values
  # of the columns "stat"
  pivot_longer(cols = ends_with("rt"),
               names_to = "stat",
               values_to = "rt"))
```

----


(ref:priorpred-stats) Prior predictive distribution of averages, maximum, and minimum values.

```{r priorpred-stats,fig.cap="(ref:priorpred-stats)", message = FALSE,size = "scriptsize", fig.height = 4 }
prior_stat %>%
  ggplot(aes(rt)) +
  geom_histogram(binwidth = 500) +
  facet_wrap(~stat, ncol = 1)
```

----

### Why are our distributions so bad?

We used much less prior information than what we really had: our priors are clearly not very realistic given what we know about reaction times for such a button pressing task. 

### What priors should we have chosen? 


# The influence of priors: sensitivity analysis {#sec:sensitivity}

## Types of priors

1. **Flat uninformative priors**:  priors as uninformative as possible. <!-- The idea behind this approach is to let the data "speak for itself" and to not bias the statistical inference with "subjective" priors. There are several issues with this approach: First, the prior is as subjective as the likelihood, and in fact, different choices of likelihood might have a much stronger impact on the posterior than different choices of priors. Second, uninformative priors are in general unrealistic because they give equal weight to any value, ignoring the fact that we do have some minimal information about our parameters of interest, at the very least, the order of magnitude (reaction times will be in milliseconds and not days, EEG signals some microvolts and not volts, etc). Finally, uninformative priors make the sampling slower and might lead to problems of convergence. Unless we have a large amount of data, it would be wise to avoid them. In the real data analyses we present in this book,  -->


2. **Regularizing priors**: <!-- If we don't have much prior information, and we have enough data (what enough means here will presently become clear when we look at specific examples), it is fine to use so-called *regularizing priors*.  --> priors that downweight extreme values (that is, they provide regularization), they are not very informative, and mostly let the likelihood dominate in determining the posteriors. <!-- These priors are theory-neutral; that is, they do not bias the parameters to values supported by any theory.  The idea behind this type of prior is to help to stabilize computation. For many applications, they perform well, but as we will see later, they tend to be problematic if we want to use Bayes factors. -->

3. **Principled priors**:  priors that encode all (or most of) the theory-neutral information that we do have. <!-- Since we generally know how our data do and do not look like, we can build priors that truly reflect the properties of potential datasets.  -->

4. **Informative priors**: There are cases where we have a lot of prior knowledge, and not much data. <!-- In general, unless we have *very* good reasons for having informative priors, we don’t want our priors to have too much influence on our posterior. An example where informative priors would be important is when we are investigating a language-impaired population from which we can't get many participants. -->

<!-- These four options constitute a continuum. The model from section \@ref(sec:simplenormal) falls between flat uninformative and regularizing priors. The priors we used were flat but they allowed for values with at least the right order of magnitude. In practical data analysis situations, we are mostly going to choose priors that fall between regularizing and principled. -->


## Revisiting the button-pressing example with different priors

What would happen if we use even wider priors for the model?

<!-- We could assume that every mean between $-10^{10}$ and $10^{10}$ ms is equally likely. Regarding the standard deviation, we could assume that any value between $0$ and $10^{10}$ is equally likely. We keep the likelihood as it is, and we encode the following priors: -->

\begin{equation}
\begin{aligned}
\mu &\sim Uniform(-10^{10}, 10^{10}) \\
\sigma &\sim Uniform(0,  10^{10}) 
\end{aligned}
(\#eq:rtpriorsflat)
\end{equation}

### In brms:

```{r, message = FALSE, cache = TRUE, tidy = FALSE}
fit_press_unif <- brm(rt ~ 1,
  data = df_noreading_data,
  family = gaussian(),
  prior = c(
      prior(uniform(-10^10, 10^10), class = Intercept),
    prior(uniform(0, 10^10), class = sigma))
)
```

---

The output of the model is virtually identical!

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

----

### What happens if we use very informative priors and they are off? 


\begin{equation}
\begin{aligned}
\mu &\sim Normal(400, 10) \\
\sigma &\sim Normal_+(100, 10) 
\end{aligned}
(\#eq:infrtpriors)
\end{equation}


```{r, message = FALSE, cache = TRUE}
fit_press_inf <- brm(rt ~ 1,
  data = df_noreading_data,
  family = gaussian(),
  prior = c(
    prior(normal(400, 10), class = Intercept),
    # brms knows that SD needs to be bounded by zero:
    prior(normal(100, 10), class = sigma)
  )
)
```

---

Even in this case, the new estimates are just a couple of milliseconds away from our previous estimates:

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

---

This doesn't mean that priors never matter:

* When there is enough data for *a certain parameter*, the likelihood will dominate
* If we are not sure  about the extent to which the posterior is influenced by our priors, we can do a *sensitivity analysis*  [for a published  example in psycholinguistics, see @vasishthProcessingChineseRelative2013]. 
* We can use prior predictive distributions to see if we are on the right order of magnitude for our priors


---


# Posterior predictive distributions {#sec:ppd} 

 <!-- After we have seen the data and obtained the posterior distributions of the parameters, we can now use the *posterior distributions* to  generate future data from the model. In other words, given the posterior distributions of the parameters of the model, the posterior predictive distribution shows how future data might look like. -->

----

Once we have the posterior distribution $p(\Theta\mid y)$, we can derive the predictions based on this distribution:

\begin{equation}
p(D_{pred}\mid y )=\int_\Theta p(D_{pred}\mid \Theta) p(\Theta\mid y)\, d\Theta
\end{equation}


<!-- Note that we are conditioning $D_{pred}$ only on $y$, we do not condition on what we don't know ($\Theta$); we integrate out the unknown parameters. This posterior predictive distribution is different from  the frequentist approach, which gives only a predictive distribution of $D_{pred}$ given our maximum likelihood estimate of $\theta$ (a point value).  -->

---

Here we can also avoid the integration, and we can even use the same function that we created before:

```{r, size = "scriptsize"}
N_obs <- nrow(df_noreading_data)
mu_samples <- posterior_samples(fit_press)$b_Intercept
sigma_samples <- posterior_samples(fit_press)$sigma
(normal_predictive_distribution_fast(
  mu_samples = mu_samples,
  sigma_samples = sigma_samples,
  N_obs))
```


## Descriptive adequacy/posterior predictive checks

<!-- We can use the posterior predictive distribution to examine the "descriptive adequacy" of our models [@Gelman14, Chapter 6; @shiffrinSurveyModelEvaluation2008]; these are called posterior predictive checks, and what we want to establish here is that the posterior predictive data look more or less similar to the observed data. Achieving descriptive adequacy means that the current data could have been generated by the model. While passing a test of descriptive adequacy is not strong evidence in favor of a model, a major failure in descriptive adequacy can be interpreted as strong evidence against a model [@shiffrinSurveyModelEvaluation2008]. Thus, posterior predictive checking is an important sanity check to assess whether the model behavior is reasonable. -->


### Could the  current data have been generated by our model?
<!-- * Yes $->$ maybe keep looking -->
<!-- * No $->$ evidence against the model -->

----

<!-- In many cases, we can simply use the plot functions from `brms` and `bayesplot` that take the model as an argument for the visualization posterior predictive checks. For example, we can use `pp_check`  to investigate how well the observed distribution of reaction times fit our model based on some number (11 and 100) samples of the posterior predictive distributions; see figures \@ref(fig:normalppc) and \@ref(fig:normalppc2). -->

(ref:normalppc) Eleven samples from the posterior predictive distribution of the model `fit_press`.
```{r normalppc, fig.cap = "(ref:normalppc)", message = FALSE,fig.height = 6}
pp_check(fit_press, nsamples = 11, type = "hist")
```

-----

(ref:normalppc2) Posterior predictive check that shows the fit of the model `fit_press` in comparison to datasets from the posterior predictive distribution.
```{r normalppc2, fig.cap = "(ref:normalppc2)" , message = FALSE, fig.height = 6}
pp_check(fit_press, nsamples = 100)
```


<!-- Notice that the real data is slightly skewed and has no values shorter than 100 ms, while the predictive distributions are centered and symmetrical; see figures \@ref(fig:normalppc) and \@ref(fig:normalppc2). This posterior predictive check shows a slight mismatch between the observed and predicted data. Can we build a better model?  We’ll come back to this issue in the next section. -->

#  Comparing different likelihoods: The log-normal likelihood {#sec:lnfirst}

<!-- Since we know that the reaction times shouldn't be normally distributed, we can choose a more realistic distribution for the likelihood. A good candidate is the log-normal distribution since a variable (such as time) that is log-normally distributed takes only positive real values and is right skewed. -->


----

If $y$ is log-normally distributed, this means that $\log(y)$ is normally distributed.^[In fact, $\log_e(y)$ or $\ln(y)$, but  we'll write it as just $log()$] 

<!-- Thus, when we model some data $y$ using the log-normal likelihood, the parameters $\mu$ and $\sigma$ are on a different scale than the data $y$. -->

<!-- We can create a log-normal distribution by exponentiating the samples of a normal distribution. See Figure \@ref(fig:logndemo). -->

\begin{equation}
\begin{aligned}
\log(y) &\sim Normal( \mu, \sigma)\\
y &\sim \exp(Normal( \mu, \sigma)) \\
y &\sim LogNormal( \mu, \sigma)
\end{aligned}
\end{equation}

\color{red}
The log-normal distribution is again defined using $\mu$ and $\sigma$, but these correspond to the mean and standard deviation of the normally distributed logarithm of the data $y$: $\log(y)$.
\color{black}



## Re-fitting a single participant pressing a button repeatedly with a log-normal likelihood  {#sec:lognormal}

<!-- If we assume that reaction times are log-normally distributed, we’ll need to change our likelihood function as follows: -->

### New likelihood:

\begin{equation}
rt_n \sim LogNormal(\mu,\sigma)
\end{equation}


### New scale for the priors: 

<!-- But now the scale of our priors needs to change!  We'll continue with the  uniform priors for ease of exposition, even though, as we mentioned earlier, these are not recommended. -->

\begin{equation}
\begin{aligned}
\mu &\sim Uniform(0, 8) \\
\sigma &\sim Uniform(0, 1) \\
\end{aligned}
(\#eq:logpriorsunif)
\end{equation}

---

Because the parameters are in a different scale than the dependent variable, their  interpretation changes:


-   *The location, $\mu$*: In our previous linear model, $\mu$ represented the grand mean (or the grand median, or grand mode, since in a normal distribution the three coincide). But now, the grand mean is $\exp(\mu +\sigma ^{2}/2)$ and the grand median is $\exp(\mu)$. <!-- We could assume that the grand median, $\exp(\mu)$,  represents the underlying time it takes to press the space bar if there would be no noise, that is, if $\sigma$ would be 0. This also means that the prior of $\mu$ is not in milliseconds, but in log(milliseconds). -->
-   *The scale, $\sigma$*: This is the standard deviation of the normal distribution of $\log(y)$. The standard deviation of a log-normal distribution with *location* $\mu$ and *scale* $\sigma$ will be  $\exp(\mu +\sigma ^{2}/2)\times \sqrt(\exp(\sigma^2)- 1)$. <!-- It's important to notice that, unlike the normal distribution, the spread of the log-normal distribution depends on both $\mu$ and $\sigma$. -->

## Prior predictive distributions

<!-- To understand the meaning of our priors in milliseconds scale, we need to take into account both the priors and the likelihood. We can do this by generating a prior predictive distribution. Notice that we can just exponentiate the samples produced by `normal_predictive_distribution_fast()` (or, alternatively, we could have edited the function and replaced `rnorm` for `rlnorm`).  -->


```{r}
N_samples <- 1000
N_obs <- nrow(df_noreading_data)
mu_samples <- runif(N_samples, 0, 8)
sigma_samples <- runif(N_samples, 0, 1)
prior_pred_ln <- exp(normal_predictive_distribution_fast(
  mu_samples = mu_samples,
  sigma_samples = sigma_samples,
  N_obs
))
```

## Distribution of statistics

```{r, size = "scriptsize"}
(prior_pred_stat_ln <- 
prior_pred_ln %>%
  group_by(iter) %>%
  summarize(
    min_rt = min(rt_pred),
    max_rt = max(rt_pred),
    average_rt = mean(rt_pred),
    median_rt = median(rt_pred)
  ) %>%
  pivot_longer(cols = ends_with("rt"), names_to = "stat", values_to = "rt"))
``` 

-----

<!-- And then we plot the distribution of some representative statistics: -->
\vspace{.1in}

(ref:priorpredlogunif) Prior predictive distribution of averages, maximum, and minimum value of the log-normal model; the x-axis is log-transformed.

```{r priorpredlogunif,fig.cap="(ref:priorpredlogunif)", message = FALSE, fig.height=3.8, size = "scriptsize", tidy = FALSE}
prior_pred_stat_ln %>%
  ggplot(aes(rt)) +
  scale_x_continuous("Reaction times in ms",
    trans = "log", breaks = c(0.001, 1, 100, 1000, 10000, 100000)) +
  geom_histogram() +
  facet_wrap(~stat, ncol = 1)
```

**We cannot not generate negative values anymore, since $\exp($any number$) > 0$**



## Better regularizing priors for the log-normal model

\begin{equation}
rt_n \sim LogNormal(\mu,\sigma)
\end{equation}


\begin{equation}
\begin{aligned}
\mu &\sim Normal(6, 1.5) \\
\sigma &\sim Normal_+(0, 1) \\
\end{aligned}
(\#eq:logpriorsnorm)
\end{equation}

<!-- Notice that while $\mu$ can be negative, the dependent variable won't, since the exponent of a negative value, $\exp($some negative value$)$, is always greater than $0$. Even before generating the prior predictive distributions, we can calculate the values within which we are 95% sure that the expected median of the observations will lie. We do this by looking at what happens at two standard deviations away from the mean of the *prior*, $\mu$, that is $6 - 2\times 1.5$ and $6 + 2\times 1.5$, and exponentiating these values: -->

-----

### Median effect for our new priors:

```{r}
c(
  lower = exp(6 - 2 * 1.5),
  higher = exp(6 + 2 * 1.5)
)
```

<!-- This means that our prior for $\mu$ is still not too informative (these are medians; the actual values generated by the distribution can be much more spread out). We plot the distribution of some representative statistics in Figure  \@ref(fig:priorpredlognorm). -->

## Prior predictive distributions

```{r, size ="scriptsize"}
N_samples <- 1000
N_obs <- nrow(df_noreading_data)
mu_samples <- rnorm(N_samples, 6, 1.5)
sigma_samples <- rtnorm(N_samples, 0, 1, a = 0)
(prior_pred_ln_better <- exp(normal_predictive_distribution_fast(
    mu_samples = mu_samples,
    sigma_samples = sigma_samples,
    N_obs
 )))
```

----

```{r, size ="scriptsize"}
(prior_pred_stat_better_ln <- prior_pred_ln_better %>%
    group_by(iter) %>%
    summarize(
        min_rt = min(rt_pred),
        max_rt = max(rt_pred),
        average_rt = mean(rt_pred),
        median_rt = median(rt_pred)
    ) %>%
 pivot_longer(cols = ends_with("rt"),
              names_to = "stat", values_to = "rt"))
```

----

\vspace{.1in}

(ref:priorpredlognorm) Prior predictive distribution of averages, maximum, and minimum value of the log-normal model  with better priors.

```{r priorpredlognorm,fig.cap="(ref:priorpredlognorm)", message = FALSE, tidy=FALSE, fig.height = 3.8, size = "scriptsize"}
prior_pred_stat_better_ln %>% ggplot(aes(rt)) +
    scale_x_continuous(trans = "log",
                       breaks = c(0.001, 1, 100, 1000, 10000, 100000)) +
    geom_histogram() +
    facet_wrap(~stat, ncol = 1) +
    coord_cartesian(xlim = c(0.001, 300000))
```

<!-- We see that the priors that we are using are still too uninformative. We could do more iterations of choosing priors and generating posterior predictive distributions until we have priors that generate realistic data. However, for most cases, priors that generate data that whose statistics (mean, median, min, max, etc.) lie roughly in the correct order of magnitude are going to be acceptable. -->

----

### brms model with reasonable priors:^[ Notice that we need to specify that the family is `lognormal()`]

```{r, message = FALSE, cache = TRUE}
fit_press_ln <- brm(rt ~ 1,
  data = df_noreading_data,
  family = lognormal(),
  prior = c(
    prior(normal(6, 1.5), class = Intercept),
    prior(normal(0, 1), class = sigma)
  )
)
```

-----

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

----

### How long does it take to press the space bar in milliseconds?

<!-- we need to transform the $\mu$ (or `Intercept` in the model) to milliseconds. Since we know that the median of the log-normal distribution is $exp(\mu)$, we do the following to calculate an estimate in milliseconds: -->
 
```{r, message = FALSE, size = "scriptsize"}
estimate_ms <- exp(posterior_samples(fit_press_ln)$b_Intercept)
c(mean = mean(estimate_ms), quantile(estimate_ms, probs = c(.025, .975)))
```

## Posterior predictive checks

(ref:logppc) Posterior predictive distribution of `fit_noreading_ln`
```{r lognppc, message=FALSE, fig.cap="(ref:logppc)", fig.height = 4}
pp_check(fit_press_ln, nsamples = 100)
```

----

### *Are the posterior predicted data now more similar to the real data, compared to the case where  we had a Normal likelihood?*

<!-- It seems so, but it’s not easy to tell. Another way to examine this would be to look at the distribution of summary statistics. We  compare the distribution of representative summary statistics for the datasets generated by different models and compare them to the observed statistics. We suspect that the normal distribution would generate reaction times that are too fast (since it's symmetrical) and that the log-normal distribution may capture the long tail better than the normal model. Based on our hunch, we compute the distribution of minimum and maximum values for the posterior predictive distributions, and we compare them with the  minimum and maximum value respectively in the data. -->
\vspace{.1in}

We suspect that the normal distribution would generate reaction times that are too fast (since it's symmetrical) and that the log-normal distribution may capture the long tail better than the normal model.


----


\vspace{.1in}
```{r ppcheckmin,fig.cap="Distribution of minimum values in a posterior predictive check. The minimum in the data is 110 ms.", fig.show="hold", message=FALSE, out.width = "45%", fig.width=4, fig.height=3.8, size ="scriptsize"}
pp_check(fit_press, type = "stat", stat = "min") +
    ggtitle("Normal model")
pp_check(fit_press_ln, type = "stat", stat = "min") +
    ggtitle("Log-normal model")
```

-----

\vspace{.1in}

```{r ppcheckmax,fig.cap="Distribution of maximum values in a posterior predictive check. The maximum in the data is 409 ms.", fig.show="hold", message=FALSE, out.width = "45%", fig.width=4, fig.height=3.8, size ="scriptsize"}
pp_check(fit_press, type = "stat", stat = "max") +
    ggtitle("Normal model")
pp_check(fit_press_ln, type = "stat", stat = "max") +
    ggtitle("Log-normal model")
```

<!-- Figure \@ref(fig:ppcheckmin) shows that the log-normal likelihood does a slightly better job since the minimum value is contained in the bulk of the log-normal distribution and in the tail of normal one. Figure \@ref(fig:ppcheckmax) shows that both models  are unable to capture the maximum value of the observed data. One explanation for this is that the log-normal-ish observations in our data are being generated by the task of pressing as fast as possible, while the observations with long reaction times are being generated by lapses of attention.  -->


---

## What did we do?


* fitted and interpreted a normal model
* looked at the effect of priors:
   - prior predictive distributions
   - sensitivity analysis
* looked at the fit of the posterior:
   - posterior predictive distribution (descriptive adequacy)
* fitted and interpreted a log-normal model
*  compared a normal model with a log-normal one

## References
