---
title: "Fitting the meta-d' model"
output: rmarkdown::html_vignette
bibliography: citations.bib
link-citations: true
vignette: >
  %\VignetteIndexEntry{Fitting the meta-d' model}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  dpi = 150, fig.width = 8, fig.height = 6,
  fig.align = "center", out.width = "75%"
)
```


## Introduction
This vignette demonstrates how to use the `hmetad` package to fit the meta-d' model [@maniscalco2012] to a canonical metacognition experiment which requires a binary decision together with a confidence rating on each trial.



## Data preparation

To get a better idea of what kind of datasets the `hmetad` package is designed for, we can start by simulating one (see `help('sim_metad')` for a description of the data simulation function):
```{r setup, message=FALSE, warning=FALSE}
library(tidyverse)
library(tidybayes)
library(hmetad)

d <- sim_metad(
  N_trials = 1000, dprime = .75, c = -.5, log_M = -1,
  c2_0 = c(.25, .75, 1), c2_1 = c(.5, 1, 1.25)
)
```

```{r data, echo=FALSE}
d <- d |> select(trial, stimulus, response, confidence)
d
```
As you can see, our dataset has a column for the `trial` number, the
presented `stimulus` on each trial (`0` or `1`), the participant's type 1 response (`0` or `1`), and the corresponding type 2 response (confidence; `1:K`). The trials in this dataset are sorted by `stimulus`, `response`, and `confidence` because this data set is simulated, but otherwise this should look very similar to the kind of data that you would get from running your own experiment.

### Type 1, type 2, and joint responses
One wrinkle is that some paradigms do not collect a separate decision (i.e., type 1 response) and confidence rating (i.e., type 2 response)---rather, they collect a single rating reflecting both the primary decision and confidence. For example, instead of a binary type 1 response and a type 2 response ranging from `1` to `K` (where `K` is the maximum confidence level), sometimes participants are asked to make a rating on a scale from `1` to `2*K`, where `1` represents a confidence `"0"` response, `K` represents an uncertain `"0"` response, `K+1` represents an uncertain `"1"` response, and `2*K` represents a confident `"1"` response. We will refer to this as a *joint response*, as it is a combination of the type 1 response and the type 2 response.

If you would like to convert joint response data into separate type 1 and type 2 responses, you can use the corresponding functions `type1_response` and `type2_response`. For example, if instead we had a dataset that looked like this:

```{r joint_response, echo=FALSE}
d.joint_response <- d |>
  ungroup() |>
  mutate(joint_response = joint_response(response, confidence, max(confidence))) |>
  select(trial, joint_response)
d.joint_response
```

Then we could convert our joint response like so:
```{r convert_joint_response}
d.joint_response |>
  mutate(
    response = type1_response(joint_response, K = 4),
    confidence = type2_response(joint_response, K = 4)
  )
```

Similarly, you can also convert the separate responses into a joint response:
```{r convert_separate_responses}
d |>
  mutate(joint_response = joint_response(response, confidence, K = 4))
```

Note that in both cases we need to specify that our confidence scale has `K=4` levels (meaning that our joint type 1/type 2 scale has `8` levels).


### Signed and unsigned binary numbers
Often datasets will use `-1` and `1` instead of `0` and `1` to represent the two possible stimuli and type 1 responses. While the `hmetad` package is designed to use the *unsigned* (`0` or `1`) version, it provides helper functions to convert between the two:

```{r to_unsigned} 
to_unsigned(c(-1, 1))
```

```{r to_signed} 
to_signed(c(0, 1))
```


### Data aggregation
Finally, to ensure that the model runs efficiently, the `hmetad` package currently requires data to be aggregated. If it is easier, the `hmetad` package will aggregate your data for you when you fit your model. But if you would like to do so manually (e.g., for plotting or follow-up analyses), the `aggregate_metad` function can do this for you:

```{r aggregate_metad}
d.summary <- aggregate_metad(d)
```

```{r echo=FALSE}
d.summary
```

The resulting data frame has three columns: `N_0` is the number of trials with `stimulus==0`, `N_1` is the number of trials with `stimulus==1`, and `N` is a matrix containing the number of joint responses for each of the two possible stimuli (with column names indicating the `stimulus` and `joint_response`).

If you would like to use variable name other than `N` for the counts, you can change the name with the `.name` argument:
```{r aggregate_metad_response}
aggregate_metad(d, .name = "y")
```


If you have other columns in your dataset (e.g., `participant` or `condition` columns) that you would like to be aggregated separately, you can simply add them to the function call:
```{r aggregate_condition, eval=FALSE}
aggregate_metad(d, participant, condition)
```


Finally, note that `aggregate_metad` automatically estimates the number of confidence levels based on the maximum value of the confidence or joint response column in your data. This usually works fine, but may fail in cases with missing data (e.g., no participant gives a confidence rating of `3` on a `4`-point scale). The number of confidence levels can be specified manually using the argument `K`:
```{r aggregate_condition_K, eval=FALSE}
aggregate_metad(d, participant, condition, K = 4)
```



## Model fitting
To fit the model, we can use the `fit_metad` function. This function is simply a wrapper around `brms::brm`, so users are **strongly** encouraged to become familiar with [the `brms` package](https://paulbuerkner.com/brms/) before model fitting. In particular, users are likely to run into convergence errors using the default (flat) priors for model parameters, so we recommend doing careful prior predictive checks to set weakly informed priors [see @schad2021toward for more information].

Since `aggregate_metad` will place our dataset has our trial counts into a column named `N` by default, we can use `N` as our response variable even if our data is not yet aggregated. To fit a model with fixed values for each parameter, then, we can use the formula `N ~ 1`:
```{r model_fitting, results=FALSE, message=FALSE, warning=FALSE}
m <- fit_metad(N ~ 1,
  data = d,
  prior = prior(normal(0, 1), class = Intercept) +
    prior(normal(0, 1), class = dprime) +
    prior(normal(0, 1), class = c) +
    prior(lognormal(0, 1), class = metac2zero1diff) +
    prior(lognormal(0, 1), class = metac2zero2diff) +
    prior(lognormal(0, 1), class = metac2one1diff) +
    prior(lognormal(0, 1), class = metac2one2diff)
)
```

```{r, echo=FALSE}
summary(m)
```

Note that here we have arbitrarily chosen to use standard normal priors for all parameters. To get a better idea of how to set informed priors, please refer to `help('set_prior', package='brms')`.

In this model, `Intercept` is our estimate of $\textrm{log}(M) = \textrm{log}\frac{\textrm{meta-}d'}{d'}$, `dprime` is our estimate of $d'$, `c` is our estimate of $c$, `metac2zero1diff` and `metac2zero2diff` are the distances between successive confidence thresholds for `"0"` responses, and `metac2one1diff` and `metac2one2diff` are the distances between successive confidence thresholds for `"1"` responses. For each parameter, `brms` shows you the posterior means (`Estimate`), posterior standard deviations (`Est. Error`), upper- and lower-95% posterior quantiles (`l-95% CI` and `u-95% CI`), as well as some convergence metrics (`Rhat`, `Bulk_ESS`, and `Tail_ESS`).

### Manual model fitting
Most users can use `fit_metad` as above to fit their models. But in some cases, it might be preferable to call `brms::brm` directly. In such cases, the `fit_metad` function is roughly analogous to the following code:

```{r eval=FALSE}
# calculate number of confidence levels
K <- n_distinct(d$confidence)

m <- brm(bf(...),
  data = aggregate_metad(d, ...),
  family = metad(K = K, ...),
  stanvars = stanvars_metad(K = K, ...),
  ...
)
```

## Extract model estimates
Once we have our fitted model, there are many estimates that we can extract from it. Although `brms` provides its own functions for extracting posterior estimates, the `hmetad` package is designed to interface well with the `tidybayes` package to make it easier to work with model posterior samples.


### Parameter estimates
First, it is often useful to extract the posterior draws of the model parameters, which we can do with `linpred_draws_metad` (which is a wrapper around `tidybayes::linpred_draws`):
```{r parameters}
draws.metad <- tibble(.row = 1) |>
  add_linpred_draws_metad(m)
```

```{r, echo=FALSE}
print(draws.metad)
```
This `tibble` has a separate row for every posterior sample and a separate column for every model parameter. This format is useful for some purposes, but it will often be useful to pivot it so that we have a separate row for each model parameter and posterior sample:
```{r}
draws.metad <- tibble(.row = 1) |>
  add_linpred_draws_metad(m, pivot_longer = TRUE)
```

```{r, echo=FALSE}
print(draws.metad)
```
Now that all of the posterior samples are stored in a single column `.value`, it is easy to get posterior summaries using e.g. `tidybayes::median_qi`:
```{r}
draws.metad |>
  median_qi()
```


### Posterior predictions
One way to evaluate model fit is to perform a *posterior predictive check*: to simulate data from the model's posterior and compare our simulated and actual data. We can do this using the function `predicted_draws_metad` (which is a wrapper around `tidybayes::predicted_draws`):
```{r post_pred1}
draws.predicted <- predicted_draws_metad(m, d.summary)
```

```{r echo=FALSE}
draws.predicted
```

In this data frame, we have all of the columns from our aggregated data `d.summary` as well as `stimulus`, `joint_response`, `response`, and `confidence` (indicating the simulated trial type), as well as `.prediction` (indicating the number of simulated trials per trial type). From here, we can plot the posterior predictions (points and error-bars) against the actual data (bars):
```{r}
draws.predicted |>
  group_by(.row, stimulus, joint_response, response, confidence) |>
  median_qi(.prediction) |>
  group_by(.row) |>
  mutate(N = t(d.summary$N[.row, ])) |>
  ggplot(aes(x = joint_response)) +
  geom_col(aes(y = N), fill = "grey80") +
  geom_pointrange(aes(y = .prediction, ymin = .lower, ymax = .upper)) +
  facet_wrap(~stimulus, labeller = label_both) +
  theme_classic(18)
```

### Posterior expectations
Usually it will be simpler to compare response probabilities rather than raw response counts. To do this, we can use the same workflow as above but using `epred_draws_metad` (which is a wrapper around `tidybayes::epred_draws`):

```{r epred_draws}
draws.epred <- epred_draws_metad(m, newdata = tibble(.row = 1))
```

```{r echo=FALSE}
draws.epred
```

```{r epred}
draws.epred |>
  group_by(.row, stimulus, joint_response, response, confidence) |>
  median_qi(.epred) |>
  group_by(.row) |>
  mutate(.true = t(response_probabilities(d.summary$N[.row, ]))) |>
  ggplot(aes(x = joint_response)) +
  geom_col(aes(y = .true), fill = "grey80") +
  geom_pointrange(aes(y = .epred, ymin = .lower, ymax = .upper)) +
  scale_alpha_discrete(range = c(.25, 1)) +
  facet_wrap(~stimulus, labeller = label_both) +
  theme_classic(18)
```


### Mean confidence
One can also compute implied values of mean confidence from the meta-d' model using `mean_confidence_draws`:
```{r mean_confidence}
tibble(.row = 1) |>
  add_mean_confidence_draws(m) |>
  median_qi(.epred) |>
  left_join(d |>
    group_by(stimulus, response) |>
    summarize(.true = mean(confidence)))
```
Here, `.epred` refers to the model-estimated mean confidence per stimulus and response, and `.true` is the empirical mean confidence.

In addition, we can compute mean confidence marginalizing over stimuli:
```{r mean_confidence2}
tibble(.row = 1) |>
  add_mean_confidence_draws(m, by_stimulus = FALSE) |>
  median_qi(.epred) |>
  left_join(d |>
    group_by(response) |>
    summarize(.true = mean(confidence)))
```
over responses:
```{r mean_confidence3}
tibble(.row = 1) |>
  add_mean_confidence_draws(m, by_response = FALSE) |>
  median_qi(.epred) |>
  left_join(d |>
    group_by(stimulus) |>
    summarize(.true = mean(confidence)))
```

or both over stimuli and responses:
```{r mean_confidence4}
tibble(.row = 1) |>
  add_mean_confidence_draws(m, by_stimulus = FALSE, by_response = FALSE) |>
  median_qi(.epred) |>
  bind_cols(d |>
    ungroup() |>
    summarize(.true = mean(confidence)))
```

### Metacognitive bias
While mean confidence is often empirically informative, it is not recommended as a measure of metacognitive bias because it is known to be confounded by type 1 response characteristics (i.e., $d'$ and $c$) and by metacognitive sensitivity [i.e., $\textrm{meta-}d'$, @sherman2018]. Instead, we recommend a new measure of metacognitive bias, $\textrm{meta-}\Delta$, which is the distance between the average of the confidence criteria and $\textrm{meta-}c$. 

$\textrm{meta-}\Delta$ can be interpreted as lying between two extremes: when $\textrm{meta-}\Delta = 0$, the observer only uses the highest confidence rating, and when $\textrm{meta-}\Delta = \infty$, the observer only uses the lowest confidence rating. 

To obtain estimates of $\textrm{meta-}\Delta$, one can use the function `metacognitive_bias_draws`:
```{r metacognitive_bias_draws}
tibble(.row = 1) |>
  add_metacognitive_bias_draws(m) |>
  median_qi()
```






### Pseudo Type 1 ROC
To obtain type 1 performance as a pseudo-type 1 ROC, we can use `add_roc1_draws`:
```{r roc1_draws}
draws.roc1 <- tibble(.row = 1) |>
  add_roc1_draws(m)
```
```{r echo=FALSE}
draws.roc1
```

Again, we have a tidy tibble with columns `.chain`, `.iteration`, and `.draw` identifying individual posterior samples, `joint_response`, `response`, and `confidence` identifying the different points on the ROC, and `.row` identifying different ROCs (since our data frame has only one row, here there is only one ROC). In addition, we also have `p_hit` and `p_fa`, which contain posterior estimates of type 1 hit rate (i.e., the probability of a `"1"` response with `confidence >= c` given `stimulus==1`) and type 1 false alarm rate (i.e., the probability of a `"1"` response with `confidence >= c` given `stimulus==0`).

For visualization, we can get posterior summaries of the ROC using `tidybayes::median_qi` and then simply plot as a line:
```{r roc1}
draws.roc1 |>
  median_qi(p_fa, p_hit) |>
  ggplot(aes(
    x = p_fa, xmin = p_fa.lower, xmax = p_fa.upper,
    y = p_hit, ymin = p_hit.lower, ymax = p_hit.upper
  )) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
  geom_errorbar(orientation = "y", width = .01) +
  geom_errorbar(orientation = "x", width = .01) +
  geom_line() +
  coord_fixed(xlim = 0:1, ylim = 0:1, expand = FALSE) +
  xlab("P(False Alarm)") +
  ylab("P(Hit)") +
  theme_bw(18)
```

### Type 2 ROC
Finally, to plot type 2 performance as a type 2 ROC, we can use `add_roc2_draws`:
```{r roc2}
draws.roc2 <- tibble(.row = 1) |>
  add_roc2_draws(m)
```

```{r echo=FALSE}
draws.roc2
```
This tibble looks the same as for `roc1_draws`, except now there are columns for `p_hit2` representing the type 2 hit rate (i.e., the probability of a correct response with `confidence >= c` given `response`) and the type 2 false alarm rate (i.e., the probability of an incorrect response with `confidence >= c` given `response`). Note that this is the response-specific type 2 ROC, so there are two separate curves for the two type 1 responses. 

We can also plot the type 2 ROC similarly:
```{r}
draws.roc2 |>
  median_qi(p_hit2, p_fa2) |>
  mutate(response = factor(response)) |>
  ggplot(aes(
    x = p_fa2, xmin = p_fa2.lower, xmax = p_fa2.upper,
    y = p_hit2, ymin = p_hit2.lower, ymax = p_hit2.upper,
    color = response
  )) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
  geom_errorbar(orientation = "y", width = .01) +
  geom_errorbar(orientation = "x", width = .01) +
  geom_line() +
  coord_fixed(xlim = 0:1, ylim = 0:1, expand = FALSE) +
  xlab("P(Type 2 False Alarm)") +
  ylab("P(Type 2 Hit)") +
  theme_bw(18)
```


## References