---
title: "Model preregistration"
output: html_notebook
---

```{r}
library(tidyverse)
library(dplyr)
library(brms)
library(cmdstanr)
options(mc.cores = 4)
library(tidybayes)
```

# Read pilot data

```{r}
df_model <- read_delim("irt_calvi.csv")

```

# Get prior

```{r}

formula = bf(
  score ~ exp(disc) * z,
  z  ~ 1 + (1|i| item) + (1 | participant_id), 
  disc ~ 1 + (1 |i| item),
  nl = TRUE
)

get_prior(
  formula = formula,
  data = df_model,
  family = bernoulli(link ="logit")
)
```

# Set prior and model formula

```{r}
prior_2pm <-
  prior("lkj(2)", class = "cor") +
  prior("normal(0, 1)",  class = "b", nlpar="disc") +
  prior("normal(0, 1)",  class = "b", nlpar="z") +
  prior("normal(0, 1)",  class = "sd", nlpar="disc") +
  prior("normal(0, 1)",  class = "sd", nlpar="z") +
  prior("constant(1)",  class = "sd", group = "participant_id", nlpar="z") +
  prior("normal(0, 1)", class = "sd", group = "item", nlpar = "disc") +
  prior("normal(0, 1)", class = "sd", group = "item", nlpar = "z")
```

# Run model

```{r}
irt_brms_2pm_model = brm(
  formula,
  data    = df_model, 
  prior   = prior_2pm,
  thin    = 5,
  iter    = 20000, 
  family = bernoulli(link ="logit"),
  cores   = 4,
  control = list(adapt_delta = 0.95),
  backend = "cmdstanr",
  file_refit = "on_change",
  file = "models/calviModel"
)
```

```{r}
irt_brms_2pm_model
```

# Item easiness

```{r}
df_unique_items <- df_model %>% 
  select(item) %>% 
  unique()
df_unique_items

df_draws_easiness <- df_unique_items %>% 
  add_linpred_draws(irt_brms_2pm_model, re_formula = ~(0 + Intercept|i| item))

df_draws_easiness
```

## Coef plot

```{r}
coef_plot_easiness <- df_draws_easiness %>%
  ggplot(aes(x=.linpred, y=item)) +
  stat_pointinterval() +
  theme_classic() +
  labs(x = "Item Easiness")  +
  theme(legend.position = "none",
        axis.title.y = element_blank(),
        panel.grid.major.y = element_line(),
        axis.text.y = element_text(margin = margin(t=0, r=0, b=0, l=30, unit="pt"))) 
coef_plot_easiness
```

# Item discrimination

```{r}
df_draws_discr <- add_linpred_draws(
  object = irt_brms_2pm_model,
  newdata = df_unique_items,
  nlpar = "disc",
  re_formula = ~(1 |i| item)
)
df_draws_discr
```

## Coef plot

```{r}
coef_plot_discr <- df_draws_discr %>%
  ungroup() %>%  # Make sure no lingering grouping
  mutate(item = factor(item),
         discrimination = exp(.linpred)) %>%
  ggplot(aes(x=discrimination, y=item)) +
  stat_pointinterval() +
  theme_classic() +
  labs(x = "Item Discrimination")  +
  theme(legend.position = "none",
        axis.title.y = element_blank(),
        panel.grid.major.y = element_line(color = "lightgray"),
        axis.text.y = element_text(margin = margin(t=0, r=0, b=0, l=30, unit="pt"))) 
coef_plot_discr
```
