---
title: "IRT Model Re-Analysis"
output: html_notebook
---

```{r setup}
library(tidyverse)
library(dplyr)
library(brms)
library(cmdstanr)
library(tidybayes)
library(corrr)
library(ggplot2)
library(patchwork)
library(tidyr)
library(distributional)
library(ggdist)
library(broom)
library(ggbeeswarm)



theme_set(theme_ggdist())

options(mc.cores = 4)
```

# Fit models

```{r}
fit_calvi    <- readRDS("models/calviModel.rds")
fit_critique <- readRDS("models/critiqueReanalysisModel.rds")
fit_drawing  <- readRDS("models/drawingModel.rds")
fit_minivlat <- readRDS("models/minivlatModel.rds")
```

# Extract participant ability estimates

```{r extract-abilities}
draws_calvi    <- spread_draws(fit_calvi,    r_participant_id__z[participant_id, Intercept]) %>%
  rename(ability_calvi    = r_participant_id__z)

draws_critique <- spread_draws(fit_critique, r_participant_id[participant_id, Intercept]) %>%
  rename(ability_critique = r_participant_id)

draws_drawing  <- spread_draws(fit_drawing,  r_participant_id[participant_id, Intercept]) %>%
  rename(ability_drawing  = r_participant_id)

draws_minivlat <- spread_draws(fit_minivlat, r_participant_id__z[participant_id, Intercept]) %>%
  rename(ability_minivlat = r_participant_id__z)

# Join all four by participant and draw
draws_combined <- draws_calvi %>%
  inner_join(draws_critique, by = c("participant_id", ".draw")) %>%
  inner_join(draws_drawing,  by = c("participant_id", ".draw")) %>%
  inner_join(draws_minivlat, by = c("participant_id", ".draw"))

# Summarize to posterior means per participant
ability_summary <- draws_combined %>%
  group_by(participant_id) %>%
  summarize(
    ability_calvi    = mean(ability_calvi),
    ability_critique = mean(ability_critique),
    ability_drawing  = mean(ability_drawing),
    ability_minivlat = mean(ability_minivlat),
    .groups = "drop"
  )

ability_summary
```

# Correlations between question sets

```{r correlation-table}
cor_table <- ability_summary %>%
  select(starts_with("ability")) %>%
  rename(
    Calvi    = ability_calvi,
    Critique = ability_critique,
    Drawing  = ability_drawing,
    Minivlat = ability_minivlat
  ) %>%
  correlate(method = "pearson") %>%
  shave() %>%
  fashion()

cor_table
```

```{r correlation-plot}
ability_summary %>%
  select(starts_with("ability")) %>%
  rename(
    Calvi    = ability_calvi,
    Critique = ability_critique,
    Drawing  = ability_drawing,
    Minivlat = ability_minivlat
  ) %>%
  correlate(method = "pearson") %>%
  rplot(print_cor = TRUE) +
  theme_classic() +
  labs(title = "Correlation of participant ability estimates across question sets") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
```

# Parallel coordinates plot

```{r parallel-coordinates}
# Join group from one of the source datasets

df_calvi    <- read_csv("irt_calvi.csv")
df_critique <- read_csv("irt_critique.csv")
df_drawing  <- read_csv("irt_drawing.csv")
df_minivlat <- read_csv("irt_minivlat.csv")

groups_lookup <- bind_rows(
  df_calvi    %>% select(participant_id, group),
  df_critique %>% select(participant_id, group),
  df_drawing  %>% select(participant_id, group),
  df_minivlat %>% select(participant_id, group)
) %>%
  distinct(participant_id, group)

ability_long <- ability_summary %>%
  left_join(groups_lookup, by = "participant_id") %>%
  pivot_longer(
    cols = starts_with("ability"),
    names_to  = "test",
    values_to = "ability"
  ) %>%
  mutate(test = recode(test,
    ability_drawing  = "Drawing",
    ability_critique = "Critique",
    ability_calvi    = "Calvi",
    ability_minivlat = "Minivlat"
  )) %>%
  mutate(test = factor(test, levels = c("Drawing", "Critique", "Calvi", "Minivlat")))

ability_long

p <- ggplot(ability_long, aes(x = test, y = ability, group = participant_id, color = group)) +
  geom_vline(xintercept = 1:4, color = "gray80", linewidth = 0.5) +
  geom_line(alpha = 0.4, linewidth = 0.7) +
  scale_x_discrete(expand = c(0.05, 0.05)) +
  scale_color_brewer(palette = "Set2") +
  labs(
    title = "Participant latent ability across question sets",
    x     = NULL,
    y     = "Latent ability",
    color = "Group"
  ) +
  theme_classic() +
  theme(
    panel.grid.major.y = element_line(color = "gray90")
  )
p
```

```{r}
ability_long_groups <- ability_summary %>%
  left_join(groups_lookup, by = "participant_id") %>%
  pivot_longer(
    cols = starts_with("ability"),
    names_to  = "model",
    values_to = "ability"
  ) %>%
  mutate(model = recode(model,
    ability_minivlat = "Mini-VLAT",
    ability_calvi    = "CALVI",
    ability_critique = "Critique",
    ability_drawing  = "Drawing"
  ))


# Compute prolific mean per model (to anchor the CIs)
prolific_means <- ability_long_groups %>%
  filter(group == "prolific") %>%
  group_by(model) %>%
  summarize(prolific_mean = mean(ability), .groups = "drop")

ci_comparisons <- ability_long_groups %>%
  group_by(model) %>%
  group_modify(~ {
    prolific_vals <- .x %>% filter(group == "prolific") %>% pull(ability)
    students_vals <- .x %>% filter(group == "students") %>% pull(ability)
    experts_vals  <- .x %>% filter(group == "experts")  %>% pull(ability)

    experts_mean  <- mean(experts_vals)
    prolific_mean <- mean(prolific_vals)

    t_exp_pro <- t.test(experts_vals, prolific_vals) %>% broom::tidy()
    t_stu_pro <- t.test(students_vals, prolific_vals) %>% broom::tidy()
    t_exp_stu <- t.test(experts_vals, students_vals)  %>% broom::tidy()

    tibble(
      group = c("experts", "students", "experts"),
      comparison = c("vs prolific", "vs prolific", "vs students"),
      mean_anchored = c(
        experts_mean,
        prolific_mean + t_stu_pro$estimate1 - t_stu_pro$estimate2,
        experts_mean
      ),
      ci_low_anchored = c(
        experts_mean  - (t_exp_pro$estimate1 - t_exp_pro$estimate2 - t_exp_pro$conf.low),
        prolific_mean + t_stu_pro$conf.low,
        experts_mean  - (t_exp_stu$estimate1 - t_exp_stu$estimate2 - t_exp_stu$conf.low)
      ),
      ci_high_anchored = c(
        experts_mean  + (t_exp_pro$conf.high - (t_exp_pro$estimate1 - t_exp_pro$estimate2)),
        prolific_mean + t_stu_pro$conf.high,
        experts_mean  + (t_exp_stu$conf.high - (t_exp_stu$estimate1 - t_exp_stu$estimate2))
      )
    )
  }) %>%
  ungroup() %>%
  mutate(
    group      = factor(group, levels = c("experts", "students", "prolific")),
    model      = factor(model, levels = c("Mini-VLAT", "CALVI", "Critique", "Drawing")),
    comparison = factor(comparison, levels = c("vs prolific", "vs students"))
  )


# Add to your existing student means data
student_means <- ability_long_groups %>%
  filter(group == "students") %>%
  group_by(model) %>%
  summarize(student_mean = mean(ability), .groups = "drop") %>%
  mutate(model = factor(model, levels = c("Mini-VLAT", "CALVI", "Critique", "Drawing")))


```

```{r}
# Plot
beeswarm_plot <- ability_long_groups %>%
  mutate(
    group = factor(group, levels = c("experts", "students", "prolific")),
    model = factor(model, levels = c("Mini-VLAT", "CALVI", "Critique", "Drawing"))
  ) %>%
  ggplot(aes(x = ability, y = group, fill = group)) +
  geom_dots(side = "topright", layout = "swarm", overflow = "compress", linewidth = 0, binwidth = unit(c(2), "mm"), dotsize = 0.9, stackratio = 1 / 0.9) +
  scale_fill_manual(values = c(
    "prolific" = "#4C78A8",
    "students" = "#F58518",
    "experts"  = "#54A24B"
  )) +
  # prolific mean marker
  geom_point(
    data = prolific_means %>%
      mutate(model = factor(model, levels = c("Mini-VLAT", "CALVI", "Critique", "Drawing"))),
    aes(x = prolific_mean, y = "prolific"),
    color = "black", shape = 18, size = 3,
    inherit.aes = FALSE
  ) +
  # comparison CIs, dodged so students-vs-prolific and students-vs-experts don't overlap
  geom_pointrange(
    data = ci_comparisons,
    aes(x = mean_anchored, xmin = ci_low_anchored, xmax = ci_high_anchored,
        y = group, linetype = comparison),
    color = "black", linewidth = 0.7, size = 0.4,
    position = position_dodge(width = 0.5),
    inherit.aes = FALSE
  ) +
  geom_vline(
    data = prolific_means %>%
      mutate(model = factor(model, levels = c("Mini-VLAT", "CALVI", "Critique", "Drawing"))),
    aes(xintercept = prolific_mean),
    linetype = "dashed", color = "gray50", linewidth = 0.5
  ) +
  geom_vline(
    data = student_means,
    aes(xintercept = student_mean),
    linetype = "dashed", color = "gray50", linewidth = 0.5
  ) +
  facet_wrap(~ model, nrow = 1) +
  theme_classic() +
  labs(x = "Latent ability", y = NULL, linetype = "Comparison") +
  theme(
    strip.background = element_blank(),
    panel.grid.major.y = element_line(color = "gray90"),
    legend.position = "bottom"
  )

beeswarm_plot
```

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

ggsave(
  "dotplot.svg",
  width  = ,
  height = 10,
  device = svglite::svglite
)
```

```{r}
ci_comparisons %>%
  mutate(
    diff = mean_anchored - case_when(
      comparison == "vs prolific" ~ prolific_means$prolific_mean[match(model, prolific_means$model)],
      comparison == "vs students" ~ student_means$student_mean[match(model, student_means$model)]
    ),
    ci_low_diff  = ci_low_anchored  - case_when(
      comparison == "vs prolific" ~ prolific_means$prolific_mean[match(model, prolific_means$model)],
      comparison == "vs students" ~ student_means$student_mean[match(model, student_means$model)]
    ),
    ci_high_diff = ci_high_anchored - case_when(
      comparison == "vs prolific" ~ prolific_means$prolific_mean[match(model, prolific_means$model)],
      comparison == "vs students" ~ student_means$student_mean[match(model, student_means$model)]
    )
  ) %>%
  mutate(
    ci_string = paste0(round(diff, 3), " [", round(ci_low_diff, 3), ", ", round(ci_high_diff, 3), "]")
  ) %>%
  select(model, group, comparison, ci_string) %>%
  arrange(model, group, comparison)
```

```{r}
# Average posterior SD of latent ability estimates
# Overall and by participant group

sd_summary <- draws_combined %>%
  group_by(participant_id) %>%
  summarize(
    sd_calvi    = sd(ability_calvi),
    sd_critique = sd(ability_critique),
    sd_drawing  = sd(ability_drawing),
    sd_minivlat = sd(ability_minivlat),
    .groups = "drop"
  ) %>%
  left_join(groups_lookup, by = "participant_id")

# Overall average SD
sd_overall <- sd_summary %>%
  summarize(across(starts_with("sd_"), mean)) %>%
  mutate(group = "overall")

# Average SD by group
sd_by_group <- sd_summary %>%
  group_by(group) %>%
  summarize(across(starts_with("sd_"), mean), .groups = "drop")

# Combined table
bind_rows(sd_overall, sd_by_group) %>%
  relocate(group) %>%
  rename(
    CALVI    = sd_calvi,
    Critique = sd_critique,
    Drawing  = sd_drawing,
    `Mini-VLAT` = sd_minivlat
  )
```

```{r}
ability_with_groups <- ability_summary %>%
  left_join(groups_lookup, by = "participant_id") %>%
  mutate(group = factor(group, levels = c("prolific", "students", "experts")))

group_colors <- c(
  "prolific" = "#4C78A8",
  "students" = "#F58518",
  "experts"  = "#54A24B"
)

# Build each panel manually using patchwork
make_scatter <- function(df, x_var, y_var) {
  r <- cor(df[[x_var]], df[[y_var]], use = "complete.obs") %>% round(2)
  ggplot(df, aes(x = .data[[x_var]], y = .data[[y_var]], color = group)) +
    geom_point(alpha = 0.5, size = 1.5) +
    geom_smooth(method = "lm", se = FALSE, color = "gray30", linewidth = 0.6) +
    scale_color_manual(values = group_colors) +
    annotate("text", x = Inf, y = Inf, label = paste0("r = ", r),
             hjust = 1.1, vjust = 1.5, size = 3.5, color = "gray20") +
    labs(x = x_var, y = y_var) +
    theme_classic() +
    theme(legend.position = "none")
}

vars <- c("ability_calvi", "ability_critique", "ability_drawing", "ability_minivlat")
labels <- c("CALVI", "Critique", "Drawing", "Mini-VLAT")
names(vars) <- labels

# Generate all lower-triangle pairs
pairs <- expand.grid(y = labels, x = labels, stringsAsFactors = FALSE) %>%
  filter(match(x, labels) < match(y, labels))

plots <- pmap(pairs, function(x, y) {
  make_scatter(ability_with_groups, vars[[x]], vars[[y]]) +
    labs(x = x, y = y)
})

# Arrange with patchwork
wrap_plots(plots, ncol = 3) +
  plot_layout(guides = "collect") &
  theme(legend.position = "bottom") &
  plot_annotation(title = "Scatterplot matrix of latent ability estimates")
A few things worth noting:

make_scatter is a helper function that builds one panel — it takes two variable names as strings, computes r, draws the points colored by group, adds an OLS line, and annotates with the correlation value.
.data[[x_var]] is an R tidyverse idiom for referencing a column by a string variable name inside ggplot(aes(...)) — the equivalent of df[x_var] in pandas.
pmap iterates over the rows of the pairs dataframe, calling make_scatter for each lower-triangle pair — similar to Array.map() over a list of argument objects in JS.
The result is 6 panels (all lower-triangle pairs from 4 instruments) assembled by patchwork's wrap_plots.
Claude Fable 5 is currently unavailable.Learn more(opens in new tab)
```

```{r}
# Full 4x4 scatterplot matrix matching reference layout

vars   <- c("ability_calvi", "ability_minivlat", "ability_critique", "ability_drawing")
labels <- c("CALVI", "Mini-VLAT", "Critique", "Drawing")
n      <- length(labels)

group_colors <- c(
  "prolific" = "#4C78A8",
  "students" = "#F58518",
  "experts"  = "#54A24B"
)

ability_with_groups <- ability_summary %>%
  left_join(groups_lookup, by = "participant_id") %>%
  mutate(group = factor(group, levels = c("prolific", "students", "experts")))

make_panel <- function(row_i, col_i) {
  x_var <- vars[col_i]
  y_var <- vars[row_i]

  if (row_i == col_i) {
    p <- ggplot(ability_with_groups, aes(x = .data[[x_var]], y = .data[[y_var]])) +
      geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "gray40", linewidth = 0.7) +
      scale_x_continuous(limits = c(0, 1)) +
      scale_y_continuous(limits = c(0, 1)) +
      coord_fixed(ratio = 1) +
      theme_bw() +
      theme(
        panel.grid.major = element_line(color = "gray88"),
        axis.title       = element_blank(),
        axis.text        = element_text(size = 7),
        legend.position  = "none"
      )
  } else {
    r <- cor(
      ability_with_groups[[x_var]],
      ability_with_groups[[y_var]],
      use = "complete.obs"
    )
    r_label <- paste0("r=", ifelse(r >= 0, "+", ""), sprintf("%.2f", r), "*")

    p <- ggplot(ability_with_groups, aes(x = .data[[x_var]], y = .data[[y_var]], color = group)) +
      geom_point(alpha = 0.6, size = 1.5) +
      scale_color_manual(values = group_colors) +
      scale_x_continuous(limits = c(0, 1)) +
      scale_y_continuous(limits = c(0, 1)) +
      coord_fixed(ratio = 1) +
      annotate("text", x = 0.02, y = 0.97, label = r_label,
               hjust = 0, vjust = 1, size = 3.2, color = "gray10", fontface = "bold") +
      theme_bw() +
      theme(
        panel.grid.major = element_line(color = "gray88"),
        axis.title       = element_blank(),
        axis.text        = element_text(size = 7),
        legend.position  = "none"
      )
  }
  p
}

# Build all 16 panels
panel_list <- vector("list", n * n)
for (row_i in 1:n) {
  for (col_i in 1:n) {
    panel_list[[(row_i - 1) * n + col_i]] <- make_panel(row_i, col_i)
  }
}

# Assemble with patchwork, add outer labels via plot_annotation + tag_levels
wrap_plots(panel_list, nrow = n, ncol = n) +
  plot_layout(guides = "collect") &
  theme(legend.position = "bottom")
```
