yolo some ROC plots from a brmsfit

rocs.Rmd

---
title: "Untitled"
output: html_document
---

Good luck trying to use this for anything!

```{r setup, include=FALSE}
library(brms)
library(patchwork)
library(tidybayes)
library(scales)
library(tidyverse)
```

```{r}
sdt_roc <- function(x, y, constant, ...) {
  
  tab <- table(x, y) # Frequency table
  tab <- tab + constant  # Add constant? (if it is zero no constant is added)
  tab <- tab / rowSums(tab) # Proportions
  # Cumulate proportions from high rating to low rating
  roc <- apply(tab, 1, function(x) rev(cumsum(rev(x))))
  roc
}
```


```{r}
data(roc6, package = "MPTinR")
dat <- tibble(roc6) %>% 
  filter(exp=="Koen_2011") %>% 
  select(-exp) %>% 
  pivot_longer(-c(id), values_to = "count") %>% 
  separate(name, c("stimulus", "rating")) %>% 
  mutate(
    stimulus = factor(str_to_title(stimulus)),
    decision = ifelse(str_detect(rating, "old"), 1, 0),
    rating = factor(
      rating,
      levels = c("3new", "2new", "1new", "1old", "2old", "3old"),
      ordered = TRUE,
      labels = 1:6
    ),
  ) %>% 
  separate(id, c("subject", "experiment"), sep = ":") %>%
  mutate(subject = as.integer(subject)) %>% 
  select(-experiment) %>% 
  filter(subject <= 9)

dat_roc <- dat %>% 
  mutate(rating = as.integer(rating)) %>% 
  select(subject, stimulus, rating, count) %>% 
  uncount(count) %>%
  group_by(subject) %>% 
  group_modify(
    ~data.frame(
      sdt_roc(
        .x$stimulus, 
        .x$rating, 
        min = 1, max = 6, n_bins = 6, constant = 0
      )
    )
  ) %>% 
  rename(far = New, hr = Old) %>% 
  mutate(zfar = qnorm(far), zhr = qnorm(hr))
```

fit models

```{r}
fit_evsdt <- brm(
  bf(
    rating | weights(count) ~ stimulus + (stimulus | subject)
  ),
  data = dat, 
  family = cumulative("probit"),
  cores = 4, 
  init = 0,
  control = list(adapt_delta = 0.95),
  file = "fit-evsdt"
)
# fit_mix <- brm(
#   bf(
#     rating | weights(n) ~ 1,
#     mu1 ~ 1 + item + (1 + item | id),
#     mu2 ~ 1,
#     theta1 ~ 1 + (1 | id),
#     family = mixture(cumulative("probit"), nmix = 2, order = "mu")
#   ), 
#   counts, 
#   cores = 4,
#   control = list(adapt_delta = 0.9),
#   file = "fit-mix-subjects"
# )
```

get rocs from models

```{r}
# Get implied ROC from models' posteriors using tidybayes
grid <- tibble(
  x = seq(-2.5, 2.5, length.out = 30),
  far = pnorm(x, lower = FALSE),
  zfar = qnorm(far)
)

evsdt_roc <- spread_draws(
  fit_evsdt, r_subject[subject, parameter], b_stimulusOld
) %>% 
  filter(parameter == "stimulusOld") %>% 
  mutate(r_stimulusOld = b_stimulusOld + r_subject) %>% 
  # right_join(grid_sub) %>%
  crossing(grid) %>% 
  mutate(hr = pnorm(x, r_stimulusOld, lower = FALSE)) %>% 
  select(-contains("b_"), -contains("r_")) %>% 
  group_by(subject) %>% 
  mutate(zhr = qnorm(hr), zfar = qnorm(far)) %>% 
  group_by(subject, far, zfar) %>% 
  mean_qi(hr, zhr)

dat_roc %>% 
  ggplot(aes(far, hr)) +
  scale_x_continuous("FAR", breaks = pretty_breaks(3)) +
  scale_y_continuous("HR", breaks = pretty_breaks(3)) +
  geom_abline(lty = 2, size = .25) +
  coord_cartesian(ylim = 0:1, xlim = 0:1) +
  geom_ribbon(
    data = evsdt_roc, col = NA, fill = "#FDE725FF",
    aes(x=far, y=hr, ymin = hr.lower, ymax = hr.upper),
    alpha = .33
  ) +
  geom_line(
    data = evsdt_roc, 
    aes(x=far, y=hr)
    ) +
  geom_point(shape = 21) +
  # geom_line(alpha = .5) +
  theme(aspect.ratio = 1) +
  facet_wrap("subject")
```

This might not work and I don't quite remember what it does or how. Try at your own risk :)