---
title: "model_critique_reanalysis"
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_critique_reanalysis.csv")
df_model 
```

# Get prior

```{r}

formula = bf(
  score  ~ 1 + (1|i| task:rubric) + (1 | participant_id), 
  disc ~ 1 + (1 |i| task:rubric)
)

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

# Set prior and model formula

```{r}
prior_2pm <-
  prior("lkj(2)", class = "cor") +
  prior("normal(0, 1)",  class = "Intercept") +
  prior("normal(0, 1)",  class = "sd") +
  prior("constant(1)",  class = "sd", group = "participant_id" ) +
  prior("normal(0, 1)", class = "sd", group = "task:rubric") +
  prior("normal(0, 1)", class = "sd", dpar = "disc") +
  prior("normal(0, 1)", class = "sd", group = "task:rubric", dpar = "disc")
```

# Run model

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

```{r}
irt_brms_2pm_model
```

# Item easiness

```{r}
df_unique_items <- df_model %>% 
  select(task, rubric) %>% 
  unique()
df_unique_items

df_draws_easiness <- df_unique_items %>% 
  add_linpred_draws(irt_brms_2pm_model, re_formula = ~(1 |i| task:rubric))
df_draws_easiness
```

## Coef plot

```{r}
coef_plot_easiness <- df_draws_easiness %>%
  mutate(task_rubric = paste(task, rubric, sep = ":")) %>%
  ggplot(aes(x=.linpred, y=task_rubric)) +
  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"))) +
    coord_cartesian(xlim = c(-2, 4)) 
coef_plot_easiness
```

```{r}
df_draws_easiness %>%
  mutate(task_rubric = paste(task, rubric, sep = ":")) %>%
  group_by(task_rubric) %>%
  summarise(
    mean   = mean(.linpred),
    median = median(.linpred),
    q5     = quantile(.linpred, 0.05),
    q95    = quantile(.linpred, 0.95),
    .groups = "drop"
  ) %>%
  arrange(desc(median))
```

```         
```

```{r}

```

# Item discrimination

```{r}
df_draws_discr <- df_unique_items %>% 
  add_linpred_draws(irt_brms_2pm_model, dpar = "disc", re_formula = ~(1 |i| task:rubric))
df_draws_discr
```

## Coef plot

```{r}
coef_plot_discr <- df_draws_discr %>%
  mutate(task_rubric = paste(task, rubric, sep = ":")) %>%
  mutate(task_rubric = factor(task_rubric),
         discrimination = exp(disc)) %>%
  ggplot(aes(x=discrimination, y=task_rubric)) +
  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"))) +
  coord_cartesian(xlim = c(0, 3)) 

coef_plot_discr
```

```{r}
install.packages("svglite")

ggsave(
  "critique_easiness.svg",
  width  = 8,
  height = 5,
  device = svglite::svglite
)
```

```{r}
# Compute numeric summaries matching the plot
discr_summary <- df_draws_discr %>%
  mutate(
    task_rubric = paste(task, rubric, sep = ":"),
    discrimination = exp(disc)
  ) %>%
  group_by(task_rubric) %>%
  median_qi(discrimination, .width = 0.95) %>%  # median_qi from tidybayes
  arrange(task_rubric) %>%
  select(task_rubric, discrimination, .lower, .upper) %>%
  mutate(across(where(is.numeric), \(x) round(x, 2)))

# Add average row
avg_row <- discr_summary %>%
  summarise(
    task_rubric = "AVERAGE",
    discrimination = mean(discrimination),
    .lower = mean(.lower),
    .upper = mean(.upper)
  ) %>%
  mutate(across(where(is.numeric), \(x) round(x, 2)))

bind_rows(discr_summary, avg_row)
```

```{r}
# Compute numeric summaries matching the plot
discr_summary <- df_draws_discr %>%
  mutate(
    task_rubric = paste(task, rubric, sep = ":"),
    discrimination = exp(disc)
  ) %>%
  group_by(task_rubric) %>%
  median_qi(discrimination, .width = 0.95) %>%  # median_qi from tidybayes
  arrange(task_rubric) %>%
  select(task_rubric, discrimination, .lower, .upper) %>%
  mutate(across(where(is.numeric), \(x) round(x, 2)))

# Add average row
avg_row <- discr_summary %>%
  summarise(
    task_rubric = "AVERAGE",
    discrimination = mean(discrimination),
    .lower = mean(.lower),
    .upper = mean(.upper)
  ) %>%
  mutate(across(where(is.numeric), \(x) round(x, 2)))

bind_rows(discr_summary, avg_row)
```

```{r}
# Compute numeric summaries matching the plot
discr_summary <- df_draws_easiness %>%
  mutate(task_rubric = paste(task, rubric, sep = ":")) %>%
  group_by(task_rubric) %>%
  median_qi(easiness, .width = 0.95) %>%  # median_qi from tidybayes
  arrange(task_rubric) %>%
  select(task_rubric, easiness, .lower, .upper) %>%
  mutate(across(where(is.numeric), \(x) round(x, 2)))

# Add average row
avg_row <- discr_summary %>%
  summarise(
    task_rubric = "AVERAGE",
    easiness = mean(easiness),
    .lower = mean(.lower),
    .upper = mean(.upper)
  ) %>%
  mutate(across(where(is.numeric), \(x) round(x, 2)))

bind_rows(discr_summary, avg_row)
```
