```{r}
#| label: setup
#| message: false
#| echo: false

library(tidyverse)
theme_set(theme_light())
library(posterior)

library(brms)
library(broom.mixed)
library(webexercises)
```

## *Takete* and *maluma*

::: callout-warning
#### Important

When working through these tutorials, always **make sure you are in the course RStudio Project** you just created.

You know you are in an RStudio Project because you can see the name of the Project in the top-right corner of RStudio, next to the light-blue cube icon.

If you see `Project (none)` in the top-right corner, that means **you are not** in an RStudio Project.

To make sure you are in the RStudio project, you can open the project by going to the project folder in File Explorer (Win) or Finder (macOS) and double click on the `.Rproj` file.
:::

We will model data from Koppensteiner et al, 2016. *Shaking Takete and Flowing Maluma. Non-Sense Words Are Associated with Motion Patterns*. DOI: [10.1371/journal.pone.0150610](https://doi.org/10.1371/journal.pone.0150610). The study looked at the relationship between sound iconicity and body motion. Read the paper abstract for an overview.

To download the file with the data right-click on the following link and download the file: [takete_maluma.txt](../data/takete_maluma.txt). (Note that tutorial files are also linked in the [Syllabus](../content.qmd)). Remember to save the file in `data/` in the course project folder.

## Read and process the data

Create a new `.Rmd` file first, save it in `code/` and name it `tutorial-w07` (the extension `.Rmd` is added automatically!).

I leave to you creating title headings in your file as you please. Remember to to add `knitr::opts_knit$set(root.dir = here::here())` in the `setup` chunk and to attach the tidyverse.

Go ahead and read the data. Since this file is a `.txt` file, you will have to use the `read_tsv()` function, rather than the `read_csv()` function. (`read_tsv()` is in the package readr, which comes with the tidyverse, so you don't need to attach any extra packages for this.) Assign the data table to a variable called `motion`.

```{r}
#| label: read-data-ex
#| eval: false

motion <- ...
```

```{r}
#| label: read-data
#| echo: false
#| message: false

motion <- read_tsv("data/takete_maluma.txt")
motion
```

The data has the following columns:

-   `Tak_Mal_Stim`: the stimulus (`Maluma` vs `Takete`).

-   `Answer`: accuracy (`CORRECT` vs `WRONG`).

-   `Corr_1_Wrong_0`: same as `Answer` but coded as `1` and `0` respectively.

-   `Rater`: study design info.

-   `Female_0`: participant's gender (alas, as binary gender).

The study specifically analysed the accuracy of the responses (`Answer`) but in this tutorial we will look instead at the response itself (whether they selected `takete` or `maluma`).

Alas, this piece of information is not coded in a column in the data, but we can create a new column based on the available info.

-   When the stimulus is `Takete` and the answer is `CORRECT` then the participant's response was `Takete`.

-   When the stimulus is `Takete` and the answer is `WRONG` then the participant's response was `Maluma`.

-   When the stimulus is `Maluma` and the answer is `CORRECT` then the participant's response was `Maluma`.

-   When the stimulus is `Maluma` and the answer is `WRONG` then the participant's response was `Takete`.

Now, go ahead and create a new column called `Response` using the `mutate()` and the `case_when()` function.

We have not encountered this function yet, so here's a challenge for you: check out [its documentation](https://dplyr.tidyverse.org/reference/case_when.html) to learn how it works. You will also need to use the AND operator `&`: this allows you to put two statements together, like `Tak_Mal_Stim == "Takete" & Answer == "CORRECT"` for "if stimulus is Takete AND answer is CORRECT".

(If you are following the documentation but `case_when()` is still giving you mysterious errors, make sure that your version of dplyr is the most current one. To do this, run `packageVersion("dplyr")` in the console. You want the output to be `1.1.0`. If it's not, you'll need to update tidyverse by re-installing it.)

```{r}
#| label: resp-ex
#| eval: false

motion <- motion %>%
  mutate(
    ...
  )
```

```{r}
#| label: resp
#| echo: false

motion <- motion %>%
  mutate(
    Response = case_when(
      Tak_Mal_Stim == "Takete" & Answer == "CORRECT" ~ "Takete",
      Tak_Mal_Stim == "Takete" & Answer == "WRONG" ~ "Maluma",
      Tak_Mal_Stim == "Maluma" & Answer == "CORRECT" ~ "Maluma",
      Tak_Mal_Stim == "Maluma" & Answer == "WRONG" ~ "Takete"
    )
  )

```

The column `Tak_Mal_Stim` has quite a long and redundant name. Let's change it to something shorter: `Stimulus`.

You can do so with the `rename()` function from dplyr. The new name goes on the left of the equals sign `=`, and the current name goes on the right.

```{r}
#| label: rename

motion <- motion %>%
  rename(Stimulus = Tak_Mal_Stim)

```

If you have done things correctly, you should have a `Response` column that looks like this (only showing relevant columns).

```{r}
#| label: show-resp

motion %>%
  # select() allows you to select specific columns
  select(Stimulus, Answer, Response)
```

Before moving on, plot the proportion of response ("takete" vs "maluma") by stimulus ("takete" vs "maluma"). Come up with the code yourself to reproduce the following plot.

```{r}
#| label: motion-plot
#| echo: false

motion %>%
  ggplot(aes(Stimulus, fill = Response)) +
  geom_bar(position = "fill") +
  labs(y = "Proportion")

```

You will learn how to create an alternative plot of this type of data below.

Now we can use `Response` as our outcome variable!

## Model `Response`

We want to model `Response` as a function of `Stimulus`. Since `Response` is binary, we have:

$$
\begin{align}
\text{Resp} & \sim Bernoulli(p) \\
logit(p) & = \beta_0 + \beta_1 \cdot \text{Stim} \\
\beta_0 & \sim Gaussian(\mu_0, \sigma_0) \\
\beta_1 & \sim Gaussian(\mu_1, \sigma_1)
\end{align}
$$

Go through the formulae above and try to understand what is going on with each.

Now, it's time to model the data with `brm()` (you need to attach the brms package):

-   Add the formula (referring to lecture materials if needed).

-   Specify the family and data.

-   Add `backend = "cmdstanr"`.

-   Remember to save the model to a file using the `file` argument (if you make changes to the model formula after running the model, remember to delete the model file first or the code will just read in the save file instead of running the model again).

Just a reminder about a few things that you want to look out for:

-   R orders the levels of categorical variables following the alphabetical order by default. Since we have not specified the order ourselves, the levels in the variables to be included in the models will be in alphabetical order.

-   Note that both the outcome `Response` and the predictor `Stimulus` have two levels `Maluma` and `Takete`.

-   However, note the difference:

    -   The model will estimate the probability of getting the **second** level in `Response` (i.e. `Takete`).

    -   For the predictor `Stimulus`, the `Intercept` will be the probability of getting `Takete` when stimulus = "maluma".

    -   To summarise: the estimates refer to the probability of getting the **second level** of the outcome (not the reference level) and the intercept refers to the probability of getting the second level of the outcome when the predictors are at their reference level!

```{r}
#| label: resp-bm-ex
#| eval: false

resp_bm <- brm(
  ...
)

```

```{r}
#| label: resp-bm
#| echo: false

resp_bm <- brm(
  Response ~ Stimulus,
  family = bernoulli(),
  data = motion,
  backend = "cmdstanr",
  file = "cache/resp_bm",
  cores = 4
)

```

Run the model and then check the summary:

```{r}
#| label: resp-bm-summ

summary(resp_bm)

```

### Random MCMC

If you compare the output of your model with the output shown here, you might notice discrepancies in values between the two.

This is because of the random component of the MCMC algorithm: every time you re-run a Bayesian model, the results will be slightly different.

One way to make your results **reproducible** (i.e. they return exactly the same values every time), is to set a **seed**, i.e. a number that is used for (pseudo-)random generation of numbers (which is used in the MCMC algorithm).

The quickest way to set a seed is to use the `seed` argument of the `brm()` function. You can set it to any number. Note that by saving the model output to a file with the `file` argument you are also ensuring reproducibility, as long as you don't delete the file.

The best (and safest) practice is to both set a seed and save the output to file.

## Results

Now, look at the population-level effects.

::: callout-note
#### Quiz 1

```{r}
#| label: quiz-1
#| results: asis
#| echo: false

opts_1 <- c(
  answer = "$\\beta_0$ and $\\beta_1$",
  "$\\mu_0$ and $\\mu_1$",
  "$\\sigma_0$ and $\\sigma_1$"
)

cat("Which parameters from the formulae above do `Intercept` and `StimulusTakete` correspond to?", longmcq(opts_1))
```

::: {.callout-tip collapse="true"}
#### Hint

The model summary of the population-level effects includes information on the beta coefficients.
:::
:::

### `Intercept`

Let's go through the `Intercept` together. I will leave working out the interpretation of `StimulusTakete` to you.

Barring small differences in the results due to the random component of the MCMC, the estimate and error for the `Intercept` are --0.81 and 0.15.

::: callout-note
#### Quiz 2

```{r}
#| label: quiz-2
#| results: asis
#| echo: false

opts_2 <- c(
  "Probabilities.",
  answer = "Log-odds.",
  "Odds."
)

cat("Which unit are these estimates in?", longmcq(opts_2))
```

::: {.callout-important collapse="true"}
#### Answer

They are in **log-odds**. This is because of the $logit(p)$ part in the formula above: to model $p$, the **logit** (mathematical) function is used, which outputs log-odds.
:::
:::

We can convert log-odds back into probabilities with the `plogis()` function. Mathematically, this R function uses the logistic function to transform log-odds into probabilities.

The logistic function is the inverse of the **logit** function (which in R is `qlogis()`. You will very rarely have to use this R function).

As an exercise, convert the estimate of the `Intercept` to a probability. What probability does --0.81 log-odds correspond to? Again, do this before moving on.

Since the `Intercept` corresponds to $\beta_0$ and the `Intercept`'s *estimate* is $\mu_0$ from the formula above, you can say that the mean probability of $\beta_0$ is the probability that corresponds to --0.81 log-odds (which you should have found is approximately 0.3 or 30% probability).

$\beta_0$ is the probability of choosing "takete" when the stimulus is "maluma". Why? To understand this, think about the order of the levels in `Response` and `Stimulus`, and how they are coded with `0`s and `1`s. Remember that things work differently in *outcomes* vs *predictors*.

> Based on the model and data, there is 30% probability of choosing the response "takete" when the stimulus is "maluma", at 95% confidence.

### `StimulusTakete`

We have calculate the probability of getting "takete" when the stimulus is "maluma".

Now, what about the coefficient `StimulusTakete`? What does this tell you? Work through this for yourself, using what you've learned so far. As a hint, remember that beta coefficients that are not the intercept tell you about the *difference* of other levels relative to the reference level (i.e. the intercept).

## Plotting

Let's plot the results.

First, let's plot the posterior probabilities of the `Intercept` and `StimulusTakete`.

### Posterior probability distributions

The first step is to obtain the draws with `as_draws_df()`.

```{r}
#| label: resp-bm-draws

resp_bm_draws <- as_draws_df(resp_bm)

resp_bm_draws
```

Go ahead, plot the density of `b_Intercept` (use `geom_density()`).

```{r}
#| label: int-dens-ex
#| eval: false

resp_bm_draws %>%
  ggplot(...) +
  ...

```

And separately, plot the density of `b_StimulusTakete`.

### Conditional posterior probability distributions

Now, let's plot the conditional posterior probability distributions (aka conditional probabilities) of the log-odds when the stimulus is "maluma" and when it is "takete".

There are a few data-manipulation steps involved before we can get our data into a format that ggplot can handle. First up, we'll use our draws to compute the conditional posterior probability distributions for each level of `Stimulus`.

Do you remember how to do this? Here's the equations:

-   $logit(p_{stim = maluma}) = \beta_0$

-   $logit(p_{stim = takete}) = \beta_0 + \beta_1$

Identify which parts of `resp_bm_draws` correspond to $\beta_0$ and $\beta_1$.

Now, mutate the data following these equations, so that you have two new columns: one called `Maluma` and one called `Takete`, with the conditional log-odds of getting a "takete" response when the stimulus is "maluma" and "takete" respectively.

```{r}
#| label: cond-ex
#| eval: false

resp_bm_draws <- resp_bm_draws %>%
  mutate(
    Maluma = ...,
    Takete = ...
  )

```

```{r}
#| label: cond
#| echo: false

resp_bm_draws <- resp_bm_draws %>%
  mutate(
    Maluma = b_Intercept,
    Takete = b_Intercept + b_StimulusTakete
  )

```

Now we've got the data we need for plotting the conditional posterior distributions, but it's not quite in the right format yet. So let's reformat the data so that it's easier to plot!

At this point, you should have two columns, `Maluma` and `Takete`, each containing log-odds values. Something like this (not necessarily containing these values exactly):

```         
Maluma     Takete
-0.59      0.83
-0.53      0.93
...
```

But plotting would be easier if we could have instead something like:

```         
Stimulus     logodds
Maluma       -0.59
Maluma       -0.53
Takete       0.83
Takete       0.93
...
```

so that you could use `logodds` as the *x*-axis and `Stimulus` as fill.

We can achieve that by using the so-called [pivot functions](https://tidyr.tidyverse.org/articles/pivot.html): they are very powerful functions, but we will only look at their simple use here.

More specifically, we need the `pivot_longer()` function: this is used in cases where you have columns named after some groups (e.g., our `Maluma` and `Takete` columns here), and instead you want the group names to be in a column of their own (e.g., one called `Stimulus`) while the values are in another column (e.g., one called `logodds`).

Here's how it works (never mind the warning):

```{r}
#| label: pivot

resp_bm_draws_l <- resp_bm_draws %>%
  # Let's select the columns from .chain to Takete
  select(.chain:Takete) %>%
  # Now let's pivot the Maluma and Takete columns
  pivot_longer(
    cols = Maluma:Takete,
    # Name of the column to store the original column names in
    names_to = "Stimulus",
    # Name of the column to store the original values in
    values_to = "logodds"
  )

resp_bm_draws_l
```

This looks more like what we want!

Now we can easily plot the conditional probability distributions of the response "takete" when the stimulus is "maluma" and "takete".

```{r}
#| label: cond-plot

resp_bm_draws_l %>%
  ggplot(aes(logodds, fill = Stimulus)) +
  geom_density(alpha = 0.5) +
  geom_rug()

```

Now, this is well and good, but often it's easier to think about our results in terms of probabilities, rather than log-odds. So what do we do if we want to plot the **probability** of choosing "takete"? Before looking ahead, try to guess for yourself what the solution will be.

It's as we've seen before! We convert log-odds into probabilities with the `plogis()` function, like so:

```{r}
#| label: cond-plot-2

resp_bm_draws_l %>%
  mutate(
    prob = plogis(logodds)
  ) %>%
  ggplot(aes(prob, fill = Stimulus)) +
  geom_density(alpha = 0.5) +
  labs(
    title = "Probability of takete response"
  )
```

The distributions themselves look pretty similar, but have a look at the numbers on the x axis. We're back in probability space, hurrah!

The take-away message from this plot: We can be quite confident that different stimuli result in different probabilities of choosing "takete", and that it is very much more probable that participants choose "takete" when the stimulus is in fact of the "takete" type!

## Credible intervals

In the lecture, we learnt about the empirical rule, quantiles and credible intervals.

Now let's look in more detail at how to obtain credible intervals (CrIs) for any probability distribution.

For sake of illustration, let's calculate CrIs of different levels (60%, 80%, 90%) for the conditional posterior probability distributions of the model above. (Note that you can calculate CrIs also for the (non-conditional) posterior probability distributions).

The following code uses the `quantile2()` function from the [posterior](https://mc-stan.org/posterior/) package to calculate the 60% CrIs.

The posterior package should already be installed, but if not make sure you install it (to know if a package is installed, go to the Packages tab in the bottom-right panel of RStudio and use the search box).

Then you also need to **attach the package**. Add the code for attaching packages in the `setup` chunk and make sure you run the chunk to attach the package!

For this function, we'll again need to use the longer-format version of the draws we obtained from our model. This data is in `resp_bm_draws_l`.

```{r}
#| label: 60-cri

resp_bm_draws_l %>%
  group_by(Stimulus) %>%
  summarise(
    # Calculate the lower and upper bounds of a 60% CrI
    q60_lo = round(quantile2(logodds, probs = 0.2), 2),
    q60_hi = round(quantile2(logodds, probs = 0.8), 2),
    # Transform into probabilities
    p_q60_lo = round(plogis(q60_lo), 2),
    p_q60_hi = round(plogis(q60_hi), 2)
  )
```

We have already seen `group_by()` and `summarise()` before, so I leave understanding this part to you.

The new piece of code is `quantile2(logodds, probs = 0.2)` and `quantile2(logodds, probs = 0.8)`:

-   The `quantile2()` function gives you the percentile based on the specified probability (the `probs` argument). So `probs = 0.2` returns the 20th percentile and `probs = 0.8` returns the 80th percentile. These percentiles give you the 60% CrI.
-   The `round()` function is used to round numbers to the nearest integer or decimal:
    -   The second argument of `round()` is the number of digits you want to round to, here `2`.
-   Finally, we convert log-odds to probabilities with the `plogis()` function (which mathematically is the logistic function).

The output of the code gives us the 60% CrI of the probability of obtaining a "takete" response when the stimulus is "maluma" and "takete" respectively, both in log-odds and probabilities.

What about the 60% CrIs of the difference in probability between "takete" and "maluma" words?

To calculate that, you need to do the following:

$$
p_{takete-maluma} = logistic(\beta_0 + \beta_1) - logistic(\beta_0)
$$

Note that it is not correct to just do $logistic(\beta_1)$. This is because to get probabilities from log-odds you always have to include the `Intercept` $\beta_0$.

In R:

```{r}
#| label: 60-cri-diff

as_draws_df(resp_bm) %>%
  mutate(
    TM_d = plogis(b_Intercept + b_StimulusTakete) - plogis(b_Intercept)
  ) %>%
  summarise(
    TM_p_q60_lo = round(quantile2(TM_d, probs = 0.2), 2),
    TM_p_q60_hi = round(quantile2(TM_d, probs = 0.8), 2),
  )
```

So at 60% probability, there is an increase of 38 to 45 percent points in the probability of a "takete" response when the word is a "takete" word relative to when it is a "maluma" word.

We could then report the following:

> We can be 60% confident that the probability of a "takete" response is between 28% and 33% ($\beta$ = -0.81, SD = 0.15) with the "maluma" stimuli, and it is between 70% and 75% with the "takete" stimuli. At 60% confidence, the difference in probability between the "takete" and "maluma" stimuli is between 38 and 45 percent points ($\beta$ = 1.78, SD = 0.21).

Nice!

I leave to you calculating the 80% and 90% CrIs. Can you tell which values to use with the `probs` argument for an 80% and 90% intervals?

::: callout-note
#### Quiz 3

```{r}
#| label: quiz-3
#| results: asis
#| echo: false

opts_3 <- c(
  "a. 10th and 90th percentiles.",
  "b. 0th and 90th percentiles.",
  answer = "c. 5th and 95th percentiles."
)

cat("Which percentiles correspond to a 90% interval?", longmcq(opts_3))
```

::: {.callout-important collapse="true"}
#### Answer

a.  These are the percentiles for an 80% CrI.
b.  These correspond to the 90% left tail of the distribution, but not a CrI.
c.  This is the correct answer. $95-5=90$.
:::
:::

## Reporting

Try and fill in the missing parts in the following report.

> We fitted a Bayesian model to ... with a ... distribution family. We included ... as the only predictor, with the levels ... and .... The predictor was coded using the default ..., with ... as the reference level.
>
> According to the model, there is ...-... probability of a "takete" response when the stimulus is a "maluma" word, at 80% confidence ($\beta$ = ..., SD = ..., 80% CrI \[..., ...\]). The probability of a "takete" response when the stimulus is a "takete" word is between ... and .... The difference in probability between the two conditions is ...-... percent points, at 80% probability ($\beta$ = ..., SD = ..., 80% CrI \[..., ...\]).

## Extra: Plotting proportions of individual participants

An alternative way to plot proportion data is to plot the proportion for individual participants (`Rater` in the `motion` data frame).

To do so we can use a combination of `summarise()`, `geom_jitter()` and `stat_summary()`, which you have seen already.

First we need to calculate the proportion of "takete" response (our "success" level, corresponding to `1`) for each participant for each stimulus.

To calculate this proportion we need to create a new column with a numeric version of the `Response` column. Since we want to calculate the proportion of "takete" responses, we use `1` for takete and `0` for maluma.

We can do that using the `ifelse()` function.

```{r}
#| label: motion-num

motion <- motion %>%
  mutate(
    Response_num = ifelse(Response == "Takete", 1, 0)
  )
motion
```

Now we can calculate the proportion of takete responses by participant (`Rater`) with the simple formula: number of takete responses / total number of responses.

Since we have coded takete as `1` and maluma as `0`, we can get the number of takete responses by simply taking the sum of `Response_num`.

We can use a special tidyverse function, `n()`, which returns the number of rows in the data frame. Since we will group the data frame by `Rater` and `Stimulus`, `n()` returns the number of rows for each participant and each stimulus!

```{r}
#| label: motion-prop

motion_prop <- motion %>%
  group_by(Rater, Stimulus) %>%
  summarise(
    takete_prop = sum(Response_num) / n(),
    # The following drops the grouping created by group_by() which we don't
    # need anymore.
    .groups = "drop"
  )
motion_prop
```

Now we can plot the proportions and add mean and confidence intervals using `geom_jitter()` and `stat_summary()`.

Before proceeding, you need to install the **Hmisc** package. There is no need to attach it (it is used by `stat_summary()` under the hood).

```{r}
#| label: prop-plot

ggplot() +
  # Proportion of each participant
  geom_jitter(
    data = motion_prop,
    aes(x = Stimulus, y = takete_prop),
    width = 0.1, alpha = 0.5
  ) +
  # Mean proportion by stimulus with CI
  stat_summary(
    data = motion,
    aes(x = Stimulus, y = Response_num, colour = Stimulus),
    fun.data = "mean_cl_boot", size = 1
  ) +
  labs(
    title = "Proportion of takete responses by participant and stimulus",
    caption = "Mean proportion is represented by coloured points with 95% bootstrapped Confidence Intervals.",
    x = "Stimulus",
    y = "Proportion"
  ) +
  # ylim(0, 1) +
  theme(legend.position = "none")
```

You will notice something new: we have specified the data inside `geom_jitter()` and `stat_summary()` instead of inside `ggplot()`. This is because the two functions need different data: `geom_jitter()` needs the data with the proportion we calculated for each participant and stimulus; `stat_summary()` needs to calculate the mean and CIs from the overall data, rather than from the proportion data we created.

We are also specifying the aesthetics within each geom/stat function, because while `x` is the same, the `y` differs!

In `stat_summary()`, the `fun.data` argument lets you specify the function you want to use for the summary statistics to be added. Here we are using the `mean_cl_boot` function, which returns the mean proportion of `Response_num` and the 95% Confidence Intervals (CIs, there are *frequentist* confidence intervals, rather than Bayesian Credible Intervals) of that mean. The CIs are calculated using a bootstrapping procedure (if you are interested in learning what that is, check the documentation of `smean.sd` from the Hmisc package).

## Summary

::: {.callout-note appearance="minimal"}
**Read and process data**

-   Use `read_tsv()` for tab separated files.

-   `case_when()` is a more flexible version of `ifelse()` that allows you to specify multiple values based on multiple conditions.

-   The pivot functions allow you to reshape the data. `pivot_longer()` makes the data table longer by moving the names of specified columns into a column of its own and by moving the values to a column of its own. Check out `pivot_wider()` for the opposite procedure.

**Modelling**

-   We can model binary outcomes with `brm(family = bernoulli())`.

-   To make your analysis reproducible, set a seed using the `seed` argument and save the model to a file using the `file` argument.

-   You can transform log-odds back to probabilities using the `plogis()` function.

**Credible intervals**

-   Calculate the lower and upper bounds of CrIs using the `quantile2()` function from the [posterior](https://github.com/stan-dev/posterior) package.
:::
