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

library(tidyverse)
theme_set(theme_light())

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

Let's pick up from where we left off last week.

In week 3, you learnt how to use ggplot2 to create bar charts and how to create different panels. Bar charts are great for categorical data.

This week you will learn how to make density and violin plots, which work very well with numeric data in combination of categorical data (for example, when you want to show values of a numeric variable across different groups).

Then, you will use the `brm()` function from the [brms](https://paul-buerkner.github.io/brms/) package to model continuous variables.

But first, let's read some data and transform it so that we can use it in the tutorial.

## Mutating data

### The `alb_vot` data

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

Let's play around with the Albanian VOT data we used in class.

The data comes from the paper by Coretta and colleagues, 2022. *Northern Tosk Albanian*. DOI: [10.1017/S0025100322000044](https://doi.org/10.1017/S0025100322000044). This is an IPA Illustration of Northern Tosk Albanian. It contains data from voiceless and voiced consonants of Nothern Tosk Albanian, as uttered in a set of words by 5 speakers.

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

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

I leave to you creating 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.

Now you can create a new code chunk and read the `alb_vot.csv` file in R. Remember how to do this?

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

alb_vot <- ...
```

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

It's one of the `read_*()` functions.
:::

```{r}
#| label: alb-vot
#| echo: false
#| message: false

alb_vot <- read_csv("data/alb_vot.csv")
```

Great!

```{r}
#| label: alb-vot-view

alb_vot
```

The data frame has the following columns:

-   `speaker`: the speaker ID.
-   `file`: the sound file the measurements are take from.
-   `label`: the label in the interval the measurements are taken from.
-   `release`: the time of the consonant's closure release from the beginning of the sound file in seconds.
-   `voi_onset`: the onset of voicing of the vowel following the target consonant from the beginning of the sound file in seconds.
-   `consonant`: the IPA transcription of the target consonant.

Now we need to create a new column with the Voice Onset Time (`vot`) values. Note that VOT is simply the difference between the onset of voicing (`voi_onset`) and the time of release (`release`). We will do this in the next section.

### Create a new column

To change or create columns, we can use the `mutate()` function from the [dplyr](https://dplyr.tidyverse.org) package. This is another tidyverse package, and it's attached when you attach the tidyverse. The `summarise()` and `count()` functions you used in Week 2 are from this same package.

Let's calculate VOT now.

```{r}
#| label: vot

alb_vot <- alb_vot %>%
  mutate(
    vot = voi_onset - release
  )
```

Since we are mutating `alb_vot`, we need to assign (`<-`) the output of `mutate()` back into `alb_vot`. Check the `alb_vot` data frame now to view the newly created column `vot` (to view the data frame, click on `alb_vot` in the Environment tab, in the top-right panel of Rstudio).

Well done!

::: callout-warning
Note that the following code, without the pipe, gives you the same result.

``` r
alb_vot <- mutate(alb_vot, vot = voi_onset - release)
```
:::

Normally, VOT is measured in milliseconds (ms), but the `vot` column is now in seconds (because `release` and `voi_onset` are in seconds).

Let's change this. We can overwrite the `vot` column like so:

```{r}
#| label: vot-2

alb_vot <- alb_vot %>%
  mutate(
    # Multiply by 1000 to get ms from s
    vot = (voi_onset - release) * 1000
  )
```

Check `alb_vot` now to see that it worked.

::: callout-tip
#### Mutate

The function `mutate()` allows you to create new columns or to change existing ones.
:::

## Density plots

::: callout-tip
#### Density plots

**Density plots** show the distribution (i.e. the probability density) of the values of a continuous variable**.**

They are created with `geom_density()`.
:::

VOT is a numeric continuous variable so density plots are appropriate.

To plot the probability density of a continuous variable, you can use the `density` geometry. Remember, all geometry functions start with `geom_`.

Fill in the ... in the following code to create a density plot of VOT values in `alb_vot`.

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

alb_vot %>%
  ggplot(aes(x = vot)) +
  ...
```

Note that to create a density plot, you only need to specify the `x`-axis. The `y`-axis is the probability density, which is automatically calculated (a bit like counts in bar charts, remember?).

### Make things cosy with a rug

The density line shows you a smoothed representation of the data distribution over the VOT values, but you might also want to see the raw data.

You can do so by adding the `rug` geometry. Go ahead and add a rug...

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

alb_vot %>%
  ggplot(aes(vot)) +
  geom_density() +
  ...
```

You should get the following:

```{r}
#| label: vot-dens-rug
#| echo: false

alb_vot %>%
  ggplot(aes(vot)) +
  geom_density() +
  geom_rug()
```

Nice huh?

::: callout-tip
#### Rug

Raw data can be shown with a **rug**, i.e. ticks on the axes that mark where the data is.

You can add a rug with `geom_rug()`.
:::

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

What can you notice about the distribution of VOT values?

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

Are there multiple peaks in the distribution?
:::
:::

## Filtering data

We are getting closer and closer to start modelling our VOT data to learn more about them.

Normally, you would approach data modelling by using the entire data set and including as many variables as it makes sense. Since you are just starting to learn how to model data with [brms](https://paul-buerkner.github.io/brms/), we will first take a pedagogical approach and start small.

You will have noticed in the density plot that there seems to be two groups of data: and you would be right! The data contains VOT values for both voiceless and voiced plosives.

For now, we will filter the data so that we include *only voiceless plosives that don't have a VOT longer than 20 ms*. (Again, we are doing this for pedagogical reasons, not because this is what you would do with this data in real-life research!).

Before we can move on onto filtering data, we first need to learn about *logical operators*.

### Logical operators

There are four main logical operators:

-   `x == y`: `x` equals `y`.

-   `x != y`: `x` is not equal to `y`.

-   `x > y`: `x` is greater than `y`.

-   `x < y`: `x` is smaller than `y`.

Logical operators return `TRUE` or `FALSE` depending on whether the statement they convey is true or false. Remember, `TRUE` and `FALSE` are logical values.

Try these out in the Console:

```{r}
#| label: logic

# This will return FALSE
1 == 2

# FALSE
"apples" == "oranges"

# TRUE
10 > 5

# FALSE
10 > 15

# TRUE
3 < 4

```

::: callout-tip
#### Logical operators

**Logical operators** are symbols that compare two objects and return either `TRUE` or `FALSE`.

The most common logical operators are `==`, `!=`, `>`, and `<`.
:::

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

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

opts_2a <- c(
  "`3 > 1`",
  answer = '`"a" = "a"`',
  '`"b" != "b"`',
  "`19 < 2`"
)

cat("a. Which of the following is not a logical operator?", longmcq(opts_2a))

opts_2b <- c(
  "`3 > c(1, 5)`",
  answer = '`c("a", "b") != c("a")`',
  '`"apple" != "apple"`',
  '`c("A", ";") == c("B", ";")`'
)

cat("b. Which of the following returns `c(FALSE, TRUE)`?", longmcq(opts_2b))
```

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

**2a.**

Check for errors in the logical operators.

**2b.**

I tricked you! The first `;` is the Roman semicolon, and the second `;` is the Greek question mark! (Have you copy-pasted from the tutorial? If not, do that!)
:::
:::

Now let's see how these work with `filter()`!

### The `filter()` function

Filtering in R with the tidyverse is straightforward. You can use the `filter()` function.

`filter()` takes one or more statements with logical operators.

Let's try this out. The following code filters the `consonant` column so that only the consonant `p` is included.

```{r}
#| label: vot-filt-1

alb_vot %>%
  filter(consonant == "p")
```

Neat! What if we want to include *all consonants except* `p`? Easy, we use the non-equal operator `!=`.

```{r}
#| label: vot-filt-2

alb_vot %>%
  filter(consonant != "p")
```

And if we want only `k` from speaker `s02`? We can include multiple statements separated by a comma!

```{r}
#| label: vot-filt-3

alb_vot %>%
  filter(consonant == "k", speaker == "s02")
```

Combining statements like this will give you only those rows where both conditions apply. You can add as many statements as you need.

Now try to filter the data so that you include only speaker `s04` and VOT values that are shorter than 20 ms.

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

alb_vot %>%
  filter(...)
```

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

Do you remember the logical operator for "smaller than"?
:::

This is all great, but what if we want to include multiple consonants?

To do that we need another operator: `%in%`.

### The `%in%` operator

::: callout-tip
#### `%in%`

The `%in%` operator is a special logical operator that returns `TRUE` if the value to the left of the operator is one of the values in the vector to its right, and `FALSE` if not.
:::

Try these in the Console:

```{r}
#| label: in

# TRUE
5 %in% c(1, 2, 5, 7)

# FALSE
"apples" %in% c("oranges", "bananas")
```

But `%in%` is even more powerful because the value on the left does not have to be a single value, but it can also be a vector! We say `%in%` is *vectorised* because it can work with vectors (most functions and operators in R are vectorised).

```{r}
#| label: in-vct

# TRUE, TRUE
c(1, 5) %in% c(4, 1, 7, 5, 8)

stocked <- c("durian", "bananas", "grapes")
needed <- c("durian", "apples")

# TRUE, FALSE
needed %in% stocked

```

Try to understand what is going on in the code above before moving on.

### Now filter the data

Now we can filter `alb_vot` to include only the voiceless plosives `p, t, k` and VOT values smaller than 20 ms.

```{r}
#| label: vot-filt-vl

alb_vot_vl <- alb_vot %>%
  filter(
    consonant %in% c("p", "t", "k"),
    vot < 20
  )
```

This should not look too alien! The first statement, `consonant %in% c("p", "t", "k")` looks at the `consonant` column and, for each row, it returns `TRUE` if the current row value is in `c("p", "t", "k")`, and `FALSE` if not.

::: {.callout-warning collapse="true"}
## Explanation

If it's not very clear what is going on, try the following code.

``` r
alb_vless <- alb_vot %>%
  mutate(voiceless = consonant %in% c("p", "t", "k"))
alb_vless

alb_vless %>%
  filter(voiceless == TRUE)
```
:::

The second statement, `vot < 20`, filters in VOT values that are less than 20 ms.

You might have noticed that we need to assign the output of `filter()`, which is a data frame, to a new variable `alb_vot_vl` ("vl" for "voiceless"). This will be the data frame you will use in the next section.

## Modelling with `brm()`

In the lecture, we have seen that we can assume voiceless VOT values to be distributed according to a Gaussian distribution (you may also have encountered Gaussian distributions under the name "normal distribution").

$$vot \sim Gaussian(\mu, \sigma)$$

We can calculate the sample mean and standard deviation easily, using the `summarise()` function, as you learnt in Week 2. Go ahead and complete the code.

```{r}
#| label: vot-msd
#| eval: false

alb_vot_vl %>%
  summarise(
    ...
  )
```

But now we want to make inference from our sample (N = 24) to the entire population of VOT values of the Albanian voiceless plosives /p/, /t/, /k/.

We want to estimate the following probability distributions and parameters:

$$
\begin{align}
vot & \sim Gaussian(\mu, \sigma) \\
\mu & \sim Gaussian(\mu_1, \sigma_1) \\
\sigma & \sim TruncGaussian(\mu_2, \sigma_2)
\end{align}
$$

In other words, we are trying to estimate the following hyperparameters: $\mu_1$, $\sigma_1$, $\mu_2$, $\sigma_2$.

We can achieve this by modelling the data using the `brm()` function from the brms package (BRM stands for Bayesian Regression Model). We will go into more details about what this means in the weeks to follow. For now, note that regression is a synonym of linear model, a type of modelling approach that you are learning in this course.

Add `library(brms)` in the `setup` chunk and run that chunk again so that the brms package is attached. Then, run the following code.

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

alb_vot_bm <- brm(
  vot ~ 1,
  family = gaussian(),
  data = alb_vot_vl,
  backend = "cmdstanr"
)
```

```{r}
#| label: vot-bm-2
#| echo: false
#| message: false

alb_vot_bm <- brm(
  vot ~ 1,
  family = gaussian(),
  data = alb_vot_vl,
  backend = "cmdstanr",
  file = "cache/vot_bm.rds"
)
```

You will see the message `Compiling Stan program...` followed by `Start sampling`. Something like in the figure below.

![](/img/mcmc-messages.png)

We will get back to what all this means in the coming weeks. For now know that [Stan](https://mc-stan.org) is a programming language that can run Bayesian models. `brm()` is an R interface to Stan, so that you can use R and R syntax to run Bayesian models without having to learn Stan! (If you feel adventurous, nobody stops you from learning Stan too, although it is not required for this course).

The "sampling" part is about the Markov-Chain Monte Carlo (or MCMC) sampling algorithm. For a quick intro on MCMC, check the Extra box below.

::: {.callout-important collapse="true"}
#### Extra: Markov-Chain Monte Carlo

::: {.callout-tip icon="false"}
#### MCMC

Markov-Chain Monte Carlo refers to a set of algorithms used to estimate parameters of unknown distributions, based on repeated draws from parameter values and observed data.
:::

You can read a more extended explanation by Elizabeth Pankratz [here](https://elizabethpankratz.github.io/bayes_stat/day1/mcmc.html).

You don't need to understand the details of this. For now, just remember that several draws are made (i.e. the algorithm is run for several iterations), and that the model runs 4 chains of these iterations.

We will learn more about chains and iterations in the coming weeks and you will know how to use information about these to diagnose the robustness of the estimated parameters.
:::

### Model summary

Now we can inspect the model output using the `summary()` function (don't confuse this with the `summarise()` function, used to get statistical summaries of data).

::: callout-warning
The summary *might look* something like this, but don't worry if your numbers are slightly different from these!

That happens due to the randomness inherent in the MCMC sampling process.
:::

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

summary(alb_vot_bm)

```

We have seen the summary of this model during the lecture. Do you remember how to interpret it?

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

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

cat("a. The 95% CrI of $\\mu$ (`Intercept`) includes 12.9.", torf(FALSE), "\n\n")

cat("b. The standard deviation $\\sigma_2$ of $\\sigma$ (`sigma`) is 0.44", torf(TRUE))
```
:::

You can use the `tidy()` function from the [broom.mixed](https://github.com/bbolker/broom.mixed) package to extract the information on the estimated parameters as a table (remember to attach the package in the `setup` chunk; you might also need to install the package).

```{r}
#| label: vot-bm-tidy

tidy(alb_vot_bm)
```

The `term` column gives us the two main parameters we wanted to estimate:

-   `(Intercept)` is actually $\mu$.

-   `sd__Observation` is $\sigma$.

Note that `sd__Observation` in the tidy table corresponds to `sigma` in the model summary above.

Now, for each of the two parameters, you get `estimate` and `std.error`. There are the hyperparameters of the probability distributions for $\mu$ and $\sigma$.

-   `(Intercept)` ($\mu$): `estimate` is $\mu_1$: 11.62 ms, `std.error` is $\sigma_1$: 0.57 ms.

-   `sd__Observation` ($\sigma$): `estimate` is $\mu_2$: 2.86 ms, `std.error` is $\sigma_2$: 0.44 ms.

In notation:

$$
\begin{align}
vot & \sim Gaussian(\mu, \sigma) \\
\mu & \sim Gaussian(11.62, 0.57) \\
\sigma & \sim TruncGaussian(2.86, 0.44)
\end{align}
$$

Based on the model summary and on the tidy summary, we also know the **95% Credible Intervals** (CrIs) of the probability distributions of $\mu$ and $sigma$. They are in the columns of the tidy summary called `conf.low` and `conf.high`.

::: {.callout-tip icon="false"}
#### Credible Interval

A **95% (Bayesian) Credible Interval** indicates that we can be **95% confident** that the true parameter value is within that interval.

You can also say that there is a **95% probability** that the true parameter value is within that interval.
:::

-   At 95% confidence, the mean ($\mu$) VOT value for Albanian voiceless stops is between 10.48 and 12.78 ms.

-   As for the standard deviation ($\sigma$), we can be 95% confident that it is between 2.18 and 3.88 ms.

You can also say:

-   There is a 95% probability that the mean ($\mu$) VOT value for Albanian voiceless stops is between 10.48 and 12.78 ms.

-   As for the standard deviation ($\sigma$), there is a 95% probabilty that it is between 2.18 and 3.88 ms.

::: {.callout-important collapse="true"}
## Extra: Frequentist Confidence Intervals

If you are familiar with or have heard of Confidence Intervals, these are not the same as Credible Intervals (even though we talk about *confidence* when talking about Credible Intervals.

Credible Intervals are *Bayesian*, while Confidence Intervals are *frequentist* (remember we mentioned the frequentist and Bayesian views of probabilities in class?).

If you want to know more about frequentist Confidence Intervals, check this web app: \<<https://rpsychologist.com/d3/ci/>\>.

Funnily enough, even experienced researchers don't really understand frequentist Confidence Intervals and tend to interpret them as Bayesian Credible Intervals. Check this paper for more info: [Robust misinterpretation of confidence intervals](https://link.springer.com/article/10.3758/s13423-013-0572-3).
:::

### Reporting

How would you report all that?

Reporting statistical models is quite formulaic, so you can use the following as a template and replace values as needed when you need to report other models.

> We fitted a Bayesian model of the VOT of Albanian voiceless plosive, using a Gaussian distribution as the distribution family of the model. According to the model, there is a 95% probability that the VOT mean is between 10.48 and 12.78 ms ($\beta$ = 11.6, SD = 0.59), and there is a 95% probability that the VOT SD is between 2.18 and 3.88 ($\beta$ = 2.87, SD = 0.43).

It is common to report both the CrIs and the estimate and standard deviation of each parameter.

It is also helpful to visualise probability distributions by plotting probability densities. You will learn how to do this next week!

In the meantime, read on to find out how to create violin plots and customise your plots!

## Violin plots

An efficient way of showing the distribution of continuous variables depending on discrete groups (like `consonant`) are the so-called **violin plots**.

Go ahead and run the following code.

```{r}
#| label: vot-consonant

alb_vot %>%
  ggplot(aes(consonant, vot)) +
  geom_violin()
```

They are called violin plots because they look like violins! The geometry for violin plots is... `geom_violin()`.

Usually, you need the following `aes`thetics:

-   `x`-axis: a categorical variable, like `consonant`.

-   `y`-axis: a continuous numeric variable, like `vot`.

### Lay over the raw data

You can lay over raw data on top of the violins. To do so, you can use the `jitter` geometry.

`geom_jitter()` creates so-called strip charts. Let's first see how these look like without violins.

```{r}
#| label: vot-jitter
alb_vot %>%
  ggplot(aes(consonant, vot)) +
  geom_jitter()
```

That doesn't look great right? The points are spread too wide.

We can fix that by specifying the argument `width` in the jitter geometry to be something smaller than `0.5`. Try a few values until you think it looks better.

```{r}
#| label: vot-jitter-1-ex
#| eval: false

alb_vot %>%
  ggplot(aes(consonant, vot)) +
  geom_jitter(width = ...)
```

When you are satisfied, create a violin plot with overlaid strip charts.

```{r}
#| label: vot-vio-jit-ex
#| eval: false

alb_vot %>%
  ggplot(aes(consonant, vot)) +
  geom_violin() +
  geom_jitter(width = ...)
```

Note that the order between the violin and jitter geometries is important. If you add the jitter geometry before the violin geometry, the violins will just cover the jittered points (see why it is called "layered" grammar of graphics?).

Now, try to create a violin and strip chart plot so that each `speaker` is plotted separately on the *y*-axis instead of `consonant`.

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

Normally, the violins are vertical (meaning, that `speaker` would be on the *x*-axis), but here I ask to create a plot with `speaker` on the *y*-axis!

You can achieve this by inverting the axes aesthetics (i.e. `x` and `y`) in `aes()`.
:::

### If-else

It's clear from this last plot where each speaker is plotted separately that there are two subgroups of consonants: there are voiced and voiceless consonants (and you have also noticed this in the density plot above).

So, what if we want to plot voiceless and voiced consonants separately?

We can use faceting, as you learnt in Week 3, or we can use the aesthetics `fill` and/or `colour`. However, first we need to create a new column that indicates if a consonant is voiceless or voiced!

We will use the `mutate()` function to create a new `voicing` column. But how do we get the values for this new column (without having to manually code each row as voiceless or voiced)?

We can use the `ifelse()` function!

The `ifelse()` function works a bit like the English sentence "If X is true, do Y, otherwise do Z". The "X" has to be a logical statement, i.e. an expression using the logical operators you learnt above.

Look at the following code.

```{r}
#| label: ifelse
#| eval: false

ifelse(10 > 9, "correct!", "wrong...")
```

You can read that as: If `10 > 9` is `TRUE`, then return `correct!`, if it's `FALSE`, then return `wrong...`.

What will that code return? Think about this and then try the code to see that you got it right.

`ifelse()` is vectorised, so you can even test multiple values, like so:

```{r}
#| label: ifelse-2

ifelse(c(1, 3, 10) > 9, "correct!", "wrong...")
```

We get `wrong...` for `1 > 9` and `3 > 9`, and `correct!` for `10 > 9`.

You can even use `%in%`!

```{r}
#| label: ifelse-3

ifelse(c("apple", "banana", "cherry") %in% c("apple", "durian", "passion fruit"), "got it!", "need to buy")
```

We will use `ifelse()` in conjunction with `mutate()` to create a new column `voicing`. But first, we need to get the list of consonants in the data set!

(You could find these manually, but why do so when you can automate things with code?)

### Finding unique values

So, the data contains the following consonants:

```{r}
#| label: cons
cons_ipa <- alb_vot$consonant

unique(cons_ipa)
```

Wait, what's that `$` in `alb_vot$consonant` there? That is a base R way of extracting a single column (in this case `consonant`) from a data frame (`alb_vot`).

::: callout-tip
#### The dollar sign \`\$\`

You can use the dollar sign `$` to extract a single column from a data frame as a vector.
:::

If you check `cons_ipa` in the Environment tab, it will say that this is `chr` vector, i.e. a character vector. This makes sense, because `consonant` is a character column in `alb_vot`.

::: callout-tip
#### Find unique values

The `unique()` function shows the unique values of a vector (it works with any type of vector, numeric, character, logical, ...).
:::

This is great, we now know all the consonants in the `consonant` column, but maybe it will be easier to type them later if we use the custom transcription in `label` rather than IPA.

```{r}
#| label: cons-2
cons_lab <- alb_vot$label

unique(cons_lab)
```

Great! Now, to make things easier, let me tell you which consonants are voiceless and which are voiced:

-   voiceless: `c("ts", "k", "p", "t", "tS", "c")`.

-   voiced: `c("b", "g", "d", "dz", "J", "dZ")`.

Now with this knowledge we can create a new column with `voiceless` and `voiced` values for each row in the data frame. Like so:

```{r}
#| label: voi

alb_vot <- alb_vot %>%
  mutate(
    voicing = ifelse(
      # IF this is TRUE
      label %in% c("ts", "k", "p", "t", "tS", "c"),
      # Write
      "voiceless",
      # OR ELSE write
      "voiced"
    )
  )
```

The `ifelse()` function is doing the following:

-   If the value in the column `label` is among (`%in%`) the values in the vector `c("ts", "k", "p", "t", "tS", "c")`.

-   Write `voiceless` in the `voicing` column.

-   Otherwise write `voiced`.

Check out the `voicing` column in `alb_vot`.

Now that we have the info we need to plot voiceless and voiced consonants separately, let's do that!

### Plotting voiceless and voiced VOT

Go ahead and fill in the code to create a plot with violins and overlayed strip chart with `voicing` on the *x*-axis.

```{r}
#| label: vot-vio-jitt-voi-ex
#| eval: false

alb_vot %>%
  ggplot(...) +
  ...
```

You should get the following:

```{r}
#| label: vot-vio-jitt-voi
#| echo: false

alb_vot %>%
  ggplot(aes(voicing, vot)) +
  geom_violin() +
  geom_jitter(width = 0.05)
```

### Let's add some colour

Let's **fill** the violins with colour. For now, just colour the violins depending on `voicing`.

```{r}
#| label: vot-vio-jit-2-ex
#| eval: false

alb_vot %>%
  ggplot(aes(voicing, vot, ...)) +
  geom_violin() +
  geom_jitter(width = 0.05)
```

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

Use the `fill` aesthetic in `aes()`.
:::

You should get this:

```{r}
#| label: vot-vio-jit-2
#| echo: false

alb_vot %>%
  ggplot(aes(voicing, vot, fill = voicing)) +
  geom_violin() +
  geom_jitter(width = 0.05)
```

Not bad! But maybe we can improve by making the fill a bit transparent?

You can achieve that with the `alpha` argument in `geom_violin()`, which takes a number between 0 (completely transparent) and 1 (completely opaque).

Add it to the code above to make the fill transparent. Try different values until you are satisfied with the output. You could also add `alpha` to `geom_jitter()` to make the points a bit transparent.

```{r}
#| label: vot-vio-jit-3-ex
#| eval: false

alb_vot %>%
  ggplot(aes(voicing, vot, fill = voicing)) +
  geom_violin(...) +
  geom_jitter(width = 0.05)
```

Wonderful! Now, what can you tell about VOT in voiceless and voiced consonants in Albanian based on this plot?

Next week we will learn how to model variables depending on different groups!

But for now, let's talk about legends and labels.

## Legends and labels in plots

### Plot titles and labels

The `labs()` function allows you to set the plot title, subtitle, axis labels and more.

Since these are strings (i.e. text), they must be quoted with `""`. Replace the `...` below with an appropriate title, like `"VOT of Albanian stops"`.

```{r}
#| label: vot-title
#| eval: false

alb_vot %>%
  ggplot(aes(voicing, vot, fill = voicing)) +
  geom_violin() +
  geom_jitter(width = 0.05, alpha = 0.2) +
  labs(
    title = ...
  )
```

### Axis labels

To change the the axis labels, you can specify a string for `x` and `y` in the `labs()` function.

I have changed the `x` label here below. Go ahead and change `y` to `"VOT (ms)"`. (Remember to also add the title).

```{r}
#| label: vot-labs
#| eval: false

alb_vot %>%
  ggplot(aes(voicing, vot, fill = voicing)) +
  geom_violin() +
  geom_jitter(width = 0.05, alpha = 0.2) +
  labs(
    title = ...,
    x = "Consonant voicing",
    ...
  )
```

### Positioning the legend

The position of the legend can be set with the `legend.position` argument in the `theme()` function.

```{r}
#| label: vot-legend

alb_vot %>%
  ggplot(aes(voicing, vot, fill = voicing)) +
  geom_violin() +
  geom_jitter(width = 0.05, alpha = 0.2) +
  theme(legend.position = "bottom")
```

Sometimes you don't want the legend at all. You can achieve that with `"none"` as the value to the `legend.position` argument.

```{r}
#| label: vot-legend-none
#| eval: false

alb_vot %>%
  ggplot(aes(voicing, vot, fill = voicing)) +
  geom_violin() +
  geom_jitter(width = 0.05, alpha = 0.2) +
  theme(legend.position = ...)

```

## Summary

::: {.callout-note appearance="minimal"}
**Transforming data**

-   `mutate()` creates new columns or mutates existing ones.

-   `filter()` can filter data based on given values.

-   With `ifelse()`, you can output values depending on a logical statement.

-   Use logical operators to create logical statements: `==`, `!=`, `>`, `<`.

-   `unique()` returns the unique values of a vector.

-   The `%in%` operator returns `TRUE` if the value on the operator's left-hand side is one of the values in the vector on the right-hand side, and `FALSE` if not.

**Plotting**

-   Use `geom_density()` to create density plots.

-   Use `geom_violin()` to create violin plots.

-   Use `geom_jitter()` to create strip charts.

-   `labs()` allows you to specify title and axis labels.

-   `theme(legend.position = …)` lets you specify the position of the legend.

**Modelling**

-   Use `brm()` with `family = gaussian()` to model the mean and standard deviation of Gaussian data.

-   `summary()` returns the summary of a `brm` model.

-   With `tidy()` from broom.mixed you can obtain a tidy table with the model summary.
:::
