---
title: "Using alternative signal distributions with the meta-d' model"
output: rmarkdown::html_vignette
bibliography: citations.bib
csl: apa.csl
link-citations: true
vignette: >
  %\VignetteIndexEntry{Using alternative signal distributions with 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
The standard meta-d' model assumes that the evidence for making type 1 decisions follows an equal-variance normal distribution and that the evidence for making type 2 decisions follows a truncated equal-variance normal distribution. However, there has been recent interest in using other distributions in signal detection theory. While not all distributions will be identifiable with the meta-d' model (e.g., one cannot simultaneously estimate unequal variances and $\textrm{meta-}d'$), the `hmetad` package allows one to specify any distribution that takes a single parameter defining the location (i.e., mean, median, or mode) of the distribution.

To demonstrate this functionality, here we implement the meta-d' model with the Gumbel-min distribution, which has been shown to provide a parsimonious explanation of recognition memory data [@meyer2026extreme].

We begin by loading necessary packages in `R`:
```{r setup, message=FALSE, warning=FALSE}
library(tidyverse)
library(brms)
library(tidybayes)
library(hmetad)
```

## Implementing a distribution function for use with `hmetad`
In order for a distribution to be used with the `hmetad` package, one needs to implement its cumulative distribution functions in both `R` and `Stan`. Specifically, since the `hmetad` package computes the model likelihood on the logarithmic scale, one must define:

* `<distribution>_lcdf(x, mu)`: the log cumulative distribution function defining $\textrm{log } P(X \le x)$ for a random variable $X \sim \textrm{distribution}(\mu)$.
* `<distribution>_lccdf(x, mu)`: the log complementary cumulative distribution function defining $\textrm{log } P(X \ge x)$ for a random variable $X \sim \textrm{distribution}(\mu)$.

Please note that for use with the `hmetad` package, these two functions **must** use the above naming scheme (i.e., must be named `<distribution>_l(c)cdf`). 

For example, we can write the gumbel min distribution functions in `R` as follows:
```{r gumbel_min_r}
gumbel_min_lcdf <- function(x, g) {
  log1p(-exp(-exp(x - g)))
}
gumbel_min_lccdf <- function(x, g) {
  -exp(x - g)
}
```

One also needs to implement the two functions in `Stan` using `brms::stanvar` so that they are available to `brms` during model fitting. Fortunately, the functions themselves will usually look almost identical to their definitions in `R`, with only minor syntactic changes and/or use of more efficient helper functions. The code below implements the same two functions in `Stan`: 
```{r gumbel_min_stan, results=FALSE}
gumbel_min <- stanvar(
  scode = "
real gumbel_min_lcdf(real x, real g) {
  return log1m_exp(-exp(x - g));
}
real gumbel_min_lccdf(real x, real g) {
  return -exp(x - g);
}",
  block = "functions"
)
```

Again, note that the name of the functions in `Stan` must match their corresponding names in `R`. With the `lcdf` and `lccdf` functions implemented in both `R` and `Stan`, the new distribution is ready to use with the `hmetad` package!


## Data simulation
Before continuing on to model fitting, this section describes how to simulate data from the meta-d' model with your custom distribution. This is a necessary step for parameter recovery to ensure that the meta-d' model is well-defined with respect to your distribution.

To simulate data, you can call the `sim_metad` function supplying the optional arguments `lcdf` and `lccdf`:
```{r data}
d <- sim_metad(
  N_trials = 10000, dprime = 1.5, c = .1, log_M = -.5,
  c2_0_diff = c(.25, .5, .25), c2_1_diff = c(.1, .5, .25),
  lcdf = gumbel_min_lcdf, lccdf = gumbel_min_lccdf
)
```


## Model fitting
Once you have some data, fitting the model is exactly the same as with the equal-variance normal distribution, only now we also need to specify two additional arguments to `fit_metad`. 

  1. The `distribution` argument should be the name of the distribution as a string. This should be the part of the function names preceding `"_lcdf"` and `"_lccdf"` in both `R` and in `Stan`.
  2. The `stanvars` argument should be the `stanvar` object you have created above containing the `Stan` code for your cumulative distribution functions.

Otherwise, one can call `fit_metad` just as with the equal-variance normal distribution! Please note, however, that the scale of all parameters will vary from distribution to distribution, so set priors accordingly. The code below shows how to fit the meta-d' model with our new `gumbel_min` distribution:
```{r model, 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(-1, 1), class = metac2zero1diff) +
    prior(lognormal(-1, 1), class = metac2zero2diff) +
    prior(lognormal(-1, 1), class = metac2zero3diff) +
    prior(lognormal(-1, 1), class = metac2one1diff) +
    prior(lognormal(-1, 1), class = metac2one2diff) +
    prior(lognormal(-1, 1), class = metac2one3diff),
  distribution = "gumbel_min", stanvars = gumbel_min,
)
```

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

The model summary can be interpreted just as with any other model, however here you can see that the model family is `metad__4__gumbel_min__absolute__multinomial`, indicating that this model indeed uses the `gumbel_min` distribution with four confidence levels and $\textrm{meta-}c = c$. 


## Model estimates
Once the model is fit, it can be post-processed like any other model from the `hmetad` package. Because alternative distributions are often understood in terms of their effects on the ROC, here we will focus on plotting both pseudo-type 1 and type 2 ROCs.

Looking at the pseudo-type 1 ROC, we can see that the `gumbel_min` distribution exhibits an asymmetry:
```{r roc1}
# psuedo type-1 ROC
tibble(.row = 1) |>
  add_roc1_draws(m, bounds = TRUE) |>
  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_point() +
  geom_line() +
  coord_fixed(xlim = 0:1, ylim = 0:1, expand = FALSE) +
  xlab("P(False Alarm)") +
  ylab("P(Hit)") +
  theme_bw(18)
```

Likewise, the `gumbel_min` distribution also has asymmetric type 2 ROCs:
```{r roc2}
# type 2 ROC
roc2_draws(m, tibble(.row = 1), bounds = TRUE) |>
  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_point() +
  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