SPM_IRT_analysis.Rmd

---
title: "Supplements to Analysing Standard Progressive Matrices (SPM-LS) with Bayesian Item Response Models"
author: "Paul Bürkner" 
output:
  html_document:
    theme: default
encoding: UTF-8
editor_options: 
  chunk_output_type: console
---

<style>
pre {
  overflow-x: auto;
}
pre code {
  word-wrap: normal;
  white-space: pre;
}
</style>

---

# {.tabset}

## Setup

```{r setup, include=FALSE, cache=FALSE}
knitr::opts_chunk$set(cache = TRUE, fig.width = 10)
```

```{r general_options, cache=FALSE}
library(dplyr)
library(tidyr)
library(tibble)
library(stringr)
library(ggplot2)
library(GGally)
library(mirt)
library(brms)
cores <- 4
iter <- 3000
warmup <- 1000
control <- list(adapt_delta = 0.95)
options(width = 160)
theme_set(theme_default())
```

```{r functions}
scale2 <- function(x, na.rm = TRUE) {
  (x - mean(x, na.rm = na.rm)) / sd(x, na.rm = na.rm)
}

center <- function(x, na.rm = TRUE) {
  x - mean(x, na.rm = na.rm)
}

r2z <- function(r) {
  log((1 + r) / (1 - r)) / 2
}

z2r <- function(z) {
  (exp(2 * z) - 1) / (exp(2 * z) + 1)
}

SW <- suppressWarnings

sum_coding <- function(x, lvls = levels(x)) {
  # codes the first category with -1
  nlvls <- length(lvls)
  stopifnot(nlvls > 1)
  cont <- diag(nlvls)[, -nlvls, drop = FALSE]
  cont[nlvls, ] <- -1
  cont <- cont[c(nlvls, 1:(nlvls - 1)), , drop = FALSE]
  colnames(cont) <- lvls[-1]
  x <- factor(x, levels = lvls)
  contrasts(x) <- cont
  x
}

dummy_coding <- function(x, lvls = levels(x)) {
  x <- factor(x, levels = lvls)
  contrasts(x) <- contr.treatment(levels(x))
  x
}

Sum <- function(df, vars, na.rm = TRUE) {
  # row-wise sum of variables vars stored in df
  vars <- enquo(vars)
  df <- select(df, !!vars)
  all_na <- apply(df, 1, function(x) all(is.na(x)))
  out <- rowSums(select(df, !!vars), na.rm = na.rm)
  ifelse(all_na, NA, out)
}

Mean <- function(df, vars, na.rm = TRUE) {
  # row-wise mean of variables vars stored in df
  vars <- enquo(vars)
  df <- select(df, !!vars)
  all_na <- apply(df, 1, function(x) all(is.na(x)))
  out <- rowMeans(select(df, !!vars), na.rm = na.rm)
  ifelse(all_na, NA, out)
}

ll <- function(y, p) {
  y * log(p) + (1 - y) * log(1 - p)
}

fit_statistic <- function(model, criterion, group, nsamples = NULL) {
  group <- enquo(group)
  subset <- NULL
  if (!is.null(nsamples)) {
    subset <- sample(seq_len(nsamples(model)), nsamples) 
  }
  ppe <- pp_expect(model, subset = subset) %>%
    t() %>%
    as.data.frame() %>%
    cbind(model$data) %>%
    gather("draw", "ppe", starts_with("V"))
  
  yrep <- posterior_predict(model, subset = subset) %>%
    t() %>%
    as.data.frame() %>%
    cbind(spm_long) %>%
    gather("draw", "yrep", starts_with("V"))
  
  ppe %>%
    mutate(yrep = yrep$yrep) %>%
    mutate(
      crit = criterion(response2, ppe),
      crit_rep = criterion(yrep, ppe)
    ) %>%
    group_by(!!group, draw) %>%
    summarise(
      crit = sum(crit), 
      crit_rep = sum(crit_rep),
      crit_diff = crit_rep - crit
    ) %>%
    mutate(draw = as.numeric(sub("^V", "", draw))) %>%
    arrange(!!group, draw) %>%
    identity()
}

theme_hist <- function(...) {
  bayesplot::theme_default() +
  theme(
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank(),
    axis.title.y = element_blank(),
    ...
  )
}
```

## Data Preparation

```{r}
answers <- c(7, 6, 8, 2, 1, 5, 1, 6, 3, 2, 4, 5)
names(answers) <- paste0("SPM", seq_along(answers))
```

```{r}
spm_long <- read.csv("data/dataset.csv")  %>%
  mutate(person = seq_len(n())) %>%
  gather("item", "response", SPM1:SPM12) %>%
  mutate(
    answer = answers[item],
    response2 = as.numeric(response == answer),
    item = as.factor(as.numeric(sub("^SPM", "", item)))
  )
```

```{r}
spm_wide <- spm_long %>% 
  select(person, item, response2) %>%
  spread("item", "response2") %>%
  select(-person)
```

## 1PL Models {.tabset}

### MCMC with hierarchical item priors

```{r}
formula_1pl <- bf(
  formula = response2 ~ 1 + (1 | item) + (1 | person),
  family = brmsfamily("bernoulli", link = "logit")
)

prior_1pl <- 
  prior("normal(0, 2)", class = "Intercept") +
  prior("normal(0, 3)", class = "sd", group = "person") + 
  prior("normal(0, 3)", class = "sd", group = "item")

fit_1pl <- brm(
  formula = formula_1pl,
  data = spm_long,
  prior = prior_1pl,
  chains = cores, iter = iter, warmup = warmup,
  cores = cores, control = control,
  file = "models/fit_1pl"
)

fit_1pl <- add_loo(fit_1pl)
```

```{r}
summary(fit_1pl)
plot(fit_1pl)
```

```{r}
person_pars_1pl <- 
  ranef(fit_1pl, summary = FALSE)$person[, , "Intercept"] 

person_sds_1pl <- apply(person_pars_1pl, 1, sd)

person_pars_1pl <- person_pars_1pl %>%
  sweep(1, person_sds_1pl, "/") %>%
  posterior_summary() %>%
	as_tibble() %>%
	rownames_to_column(var = "person") %>%
  mutate(person = as.numeric(person))

saveRDS(person_pars_1pl, "results/person_pars_1pl")
```

```{r}
person_pars_1pl %>%
  arrange(Estimate) %>%
  mutate(id2 = seq_len(n())) %>%
	ggplot(aes(id2, Estimate, ymin = Q2.5, ymax = Q97.5)) +
	geom_pointrange(alpha = 0.7) +
	coord_flip() +
	labs(x = "Person Number (sorted after Estimate)") +
  theme(
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  )
```

```{r}
person_fit_1pl <- fit_statistic(
  fit_1pl, criterion = ll, group = person,
  nsamples = 1000
)

person_diff_1pl <- person_fit_1pl %>%
  group_by(person) %>%
  summarise(bp = mean(crit_diff > 0))

person_max_1pl <- which.max(person_diff_1pl$bp)
person_fit_1pl %>% 
  filter(person == person_max_1pl) %>%
  ggplot(aes(crit_diff)) +
  geom_histogram() +
  theme_hist()
```

```{r, fig.height=4}
item_pars_1pl <- coef(fit_1pl, summary = FALSE)$item[, , "Intercept"] %>%
  posterior_summary() %>%
	as_tibble() %>%
	rownames_to_column() %>%
	rename(item = "rowname") %>%
	mutate(
	  item = as.numeric(item),
	  method = "MCMC-H"
	)

saveRDS(item_pars_1pl, "results/item_pars_1pl")
```

```{r}
item_fit_1pl <- fit_statistic(
  fit_1pl, criterion = ll, group = item,
  nsamples = 1000
)

item_fit_1pl %>% 
  ggplot(aes(crit_diff)) +
  geom_histogram() +
  facet_wrap("item", scales = "free") +
  theme_hist()
```

### MCMC with non-hierarchical item priors

```{r}
formula_1pl_nh <- bf(
  formula = response2 ~ 0 + item + (1 | person),
  family = brmsfamily("bernoulli", link = "logit")
)

prior_1pl_nh <- 
  prior("normal(0, 5)", class = "b") +
  prior("normal(0, 3)", class = "sd", group = "person")

fit_1pl_nh <- brm(
  formula = formula_1pl_nh,
  data = spm_long,
  prior = prior_1pl_nh,
  chains = cores, iter = iter, warmup = warmup,
  cores = cores, control = control,
  file = "models/fit_1pl_nh"
)

fit_1pl_nh <- add_loo(fit_1pl_nh)
```

```{r}
summary(fit_1pl_nh)
plot(fit_1pl_nh)
```

```{r}
person_pars_1pl_nh <- 
  ranef(fit_1pl_nh, summary = FALSE)$person[, , "Intercept"] 

person_sds_1pl_nh <- apply(person_pars_1pl_nh, 1, sd)

person_pars_1pl_nh <- person_pars_1pl_nh %>%
  sweep(1, person_sds_1pl_nh, "/") %>%
  posterior_summary() %>%
	as_tibble() %>%
	rownames_to_column(var = "person") %>%
  mutate(person = as.numeric(person))

saveRDS(person_pars_1pl_nh, "results/person_pars_1pl_nh")
```

```{r}
person_pars_1pl_nh %>%
  arrange(Estimate) %>%
  mutate(id2 = seq_len(n())) %>%
	ggplot(aes(id2, Estimate, ymin = Q2.5, ymax = Q97.5)) +
	geom_pointrange(alpha = 0.7) +
	coord_flip() +
	labs(x = "Person Number (sorted after Estimate)") +
  theme(
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  )
```

```{r}
person_fit_1pl_nh <- fit_statistic(
  fit_1pl_nh, criterion = ll, group = person,
  nsamples = 1000
)

person_diff_1pl_nh <- person_fit_1pl_nh %>%
  group_by(person) %>%
  summarise(bp = mean(crit_diff > 0))

person_max_1pl_nh <- which.max(person_diff_1pl_nh$bp)
person_fit_1pl_nh %>% 
  filter(person == person_max_1pl_nh) %>%
  ggplot(aes(crit_diff)) +
  geom_histogram() +
  theme_hist()
```

```{r, fig.height=4}
item_pars_1pl_nh <- brms::fixef(fit_1pl_nh, summary = FALSE) %>%
  posterior_summary() %>%
	as_tibble() %>%
	rownames_to_column() %>%
	rename(item = "rowname") %>%
	mutate(
	  item = as.numeric(item),
	  item = item - 0.1,
	  method = "MCMC-NH"
	)

saveRDS(item_pars_1pl_nh, "results/item_pars_1pl_nh")
```

```{r}
item_fit_1pl_nh <- fit_statistic(
  fit_1pl_nh, criterion = ll, group = item,
  nsamples = 1000
)

item_fit_1pl_nh %>% 
  ggplot(aes(crit_diff)) +
  geom_histogram() +
  facet_wrap("item", scales = "free") +
  theme_hist()
```


### MAP

```{r, results="hide"}
parprior_1pl_map <- list(c(seq(2, 46, 4), "norm", 0, 5))

fit_1pl_map <- mirt(
  spm_wide, model = 1, 
  itemtype = "Rasch", SE = TRUE,
  technical = list(NCYCLES = 10000), 
  parprior = parprior_1pl_map
)
```

```{r}
person_pars_1pl_map <- fscores(fit_1pl_map, full.scores.SE = TRUE) %>%
  as.data.frame() %>%
  rename(Estimate = "F1", Est.Error = "SE_F1") %>%
  mutate(
    # person parameters are not scaled appropriately by mirt for the 1PL model
    scale = sd(Estimate) / sd(person_pars_1pl$Estimate),
    Estimate = Estimate / scale,
    Est.Error = Est.Error / scale,
    Q2.5 = Estimate - 1.96 * Est.Error, 
    Q97.5 = Estimate + 1.96 * Est.Error,
    person = seq_len(n())
  ) %>%
  select(-scale)

saveRDS(person_pars_1pl_map, "results/person_pars_1pl_map")
```

```{r}
item_pars_1pl_map <- coef(fit_1pl_map, as.data.frame = TRUE) %>%
  as.data.frame() %>%
  rownames_to_column() %>%
  filter(str_detect(rowname, "\\.d$")) %>%
  mutate(
    item = str_match(rowname, "^X[[:digit:]]+") %>%
      str_remove("^X") %>%
      as.numeric(),
    item = item - 0.2,
    method = "MAP"
  ) %>%
  rename(Estimate = par, Q2.5 = CI_2.5, Q97.5 = CI_97.5)

saveRDS(item_pars_1pl_map, "results/item_pars_1pl_map")
```


### ML

```{r, results="hide"}
fit_1pl_ml <- mirt(spm_wide, model = 1, itemtype = "Rasch", SE = TRUE,
                   technical = list(NCYCLES = 10000))
```

```{r}
person_pars_1pl_ml <- fscores(fit_1pl_ml, full.scores.SE = TRUE) %>%
  as.data.frame() %>%
  rename(Estimate = "F1", Est.Error = "SE_F1") %>%
  mutate(
    # person parameters are not scaled appropriately by mirt for the 1PL model
    scale = sd(Estimate) / sd(person_pars_1pl$Estimate),
    Estimate = Estimate / scale,
    Est.Error = Est.Error / scale,
    Q2.5 = Estimate - 1.96 * Est.Error, 
    Q97.5 = Estimate + 1.96 * Est.Error,
    person = seq_len(n())
  ) %>%
  select(-scale)

saveRDS(person_pars_1pl_ml, "results/person_pars_1pl_ml")
```

```{r}
item_pars_1pl_ml <- coef(fit_1pl_ml, as.data.frame = TRUE) %>%
  as.data.frame() %>%
  rownames_to_column() %>%
  filter(str_detect(rowname, "\\.d$")) %>%
  mutate(
    item = str_match(rowname, "^X[[:digit:]]+") %>%
      str_remove("^X") %>%
      as.numeric(),
    item = item - 0.3,
    method = "ML"
  ) %>%
  rename(Estimate = par, Q2.5 = CI_2.5, Q97.5 = CI_97.5)

saveRDS(item_pars_1pl_ml, "results/item_pars_1pl_ml")
```


### Joint plots: Person parameters

```{r}
# plot comparison of person parameters
person_pars_1pl_all <- person_pars_1pl %>% 
  bind_cols(person_pars_1pl_ml) %>%
  mutate(
    UIW = Q97.5 - Q2.5,
    UIW1 = Q97.51 - Q2.51
  )

person_pars_1pl_all %>%
	ggplot(aes(Estimate, Estimate1, color = UIW)) +
	geom_point() +
	geom_abline() +
	scale_color_viridis_c() +
	lims(x = c(-3, 2), y = c(-3, 2)) +
	labs(
		x = "Estimate (MCMC-H)",
	  y = "Estimate (ML)",
		color = "UIW (MCMC-H)"
	)

person_pars_1pl_all %>%
	ggplot(aes(UIW, UIW1, color = Estimate)) +
	geom_point() +
	geom_abline() +
	scale_color_viridis_c() +
	lims(x = c(1, 2.4), y = c(1, 2.4)) +
	labs(
		x = "UIW (MCMC-H)",
	  y = "UIW (ML)",
		color = "Estimate (MCMC-H)"
	)
```

### Joint plots: Item parameters
  
```{r}
item_pars_1pl %>%
  bind_rows(
    item_pars_1pl_nh,
    item_pars_1pl_map, 
    item_pars_1pl_ml
  ) %>%
  mutate(
    method = factor(method, levels = c("MCMC-H", "MCMC-NH", "MAP", "ML"))
  ) %>%
	ggplot(aes(item, Estimate, ymin = Q2.5, ymax = Q97.5, color = method)) +
  scale_x_continuous(breaks = 1:12) +
	geom_pointrange() +
  scale_color_viridis_d(name = "Method") +
	coord_flip() +
	labs(x = "Item Number")
```


## 2PL Models {.tabset}

### MCMC with hierarchical item priors

```{r}
formula_2pl <- bf(
  response2 ~ beta + exp(logalpha) * theta,
  nl = TRUE,
  theta ~ 0 + (1 | person),
  beta ~ 1 + (1 |i| item),
  logalpha ~ 1 + (1 |i| item),
  family = brmsfamily("bernoulli", link = "logit")
)

prior_2pl <- 
  prior("normal(0, 2)", class = "b", nlpar = "beta") +
  prior("normal(0, 1)", class = "b", nlpar = "logalpha") +
  prior("normal(0, 1)", class = "sd", group = "person", nlpar = "theta") + 
  prior("normal(0, 3)", class = "sd", group = "item", nlpar = "beta") +
  prior("normal(0, 1)", class = "sd", group = "item", nlpar = "logalpha")

fit_2pl <- brm(
  formula = formula_2pl,
  data = spm_long,
  prior = prior_2pl,
  chains = cores, iter = iter, warmup = warmup,
  cores = cores, control = control, inits = 0,
  file = "models/fit_2pl"
)

fit_2pl <- add_loo(fit_2pl)
```

```{r}
summary(fit_2pl)
plot(fit_2pl)
```

```{r}
person_pars_2pl <- 
  ranef(fit_2pl, summary = FALSE)$person[, , "theta_Intercept"] 

person_sds_2pl <- apply(person_pars_2pl, 1, sd)

person_pars_2pl <- person_pars_2pl %>%
  sweep(1, person_sds_2pl, "/") %>%
  posterior_summary() %>%
	as_tibble() %>%
	rownames_to_column(var = "person") %>%
  mutate(person = as.numeric(person))

saveRDS(person_pars_2pl, "results/person_pars_2pl")
```

```{r}
person_pars_2pl %>%
  arrange(Estimate) %>%
  mutate(id2 = seq_len(n())) %>%
	ggplot(aes(id2, Estimate, ymin = Q2.5, ymax = Q97.5)) +
	geom_pointrange(alpha = 0.7) +
	coord_flip() +
	labs(x = "Person Number (sorted after Estimate)") +
  theme(
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  )
```

```{r}
person_fit_2pl <- fit_statistic(
  fit_2pl, criterion = ll, group = person,
  nsamples = 1000
)

person_diff_2pl <- person_fit_2pl %>%
  group_by(person) %>%
  summarise(bp = mean(crit_diff > 0))

person_max_2pl <- which.max(person_diff_2pl$bp)
person_fit_2pl %>% 
  filter(person == person_max_2pl) %>%
  ggplot(aes(crit_diff)) +
  geom_histogram() +
  theme_hist()
```

```{r, fig.height=4}
item_pars_2pl <- coef(fit_2pl, summary = FALSE)$item

# locations
beta <- item_pars_2pl[, , "beta_Intercept"] %>%
  posterior_summary() %>%
	as_tibble() %>%
	rownames_to_column()

# slopes
alpha <- item_pars_2pl[, , "logalpha_Intercept"] %>%
  exp() %>%
  sweep(1, person_sds_2pl, "*") %>%
  posterior_summary() %>%
	as_tibble() %>%
	rownames_to_column()

item_pars_2pl <- bind_rows(beta, alpha, .id = "nlpar") %>%
	rename(item = "rowname") %>%
	mutate(item = as.numeric(item)) %>%
	mutate(
	  nlpar = factor(nlpar, labels = c("Location", "Slope")),
	  method = "MCMC-H"  
	)

saveRDS(item_pars_2pl, "results/item_pars_2pl")
```

```{r}
item_fit_2pl <- fit_statistic(
  fit_2pl, criterion = ll, group = item,
  nsamples = 1000
)

item_fit_2pl %>% 
  ggplot(aes(crit_diff)) +
  geom_histogram() +
  facet_wrap("item", scales = "free") +
  theme_hist()
```


### MCMC with non-hierarchical item priors

```{r}
formula_2pl_nh <- bf(
  response2 ~ beta + exp(logalpha) * theta,
  nl = TRUE,
  theta ~ 0 + (1 | person),
  beta ~ 0 + item,
  logalpha ~ 0 + item,
  family = brmsfamily("bernoulli", link = "logit")
)

prior_2pl_nh <- 
  prior("normal(0, 5)", class = "b", nlpar = "beta") +
  prior("normal(0, 1)", class = "b", nlpar = "logalpha") +
  prior("normal(0, 1)", class = "sd", group = "person", nlpar = "theta")

fit_2pl_nh <- brm(
  formula = formula_2pl_nh,
  data = spm_long,
  prior = prior_2pl_nh,
  chains = cores, iter = iter, warmup = warmup,
  cores = cores, control = control, inits = 0,
  file = "models/fit_2pl_nh"
)

fit_2pl_nh <- add_loo(fit_2pl_nh)
```

```{r}
summary(fit_2pl_nh)
plot(fit_2pl_nh)
```

```{r}
person_pars_2pl_nh <- 
  ranef(fit_2pl_nh, summary = FALSE)$person[, , "theta_Intercept"] 

person_sds_2pl_nh <- apply(person_pars_2pl_nh, 1, sd)

person_pars_2pl_nh <- person_pars_2pl_nh %>%
  sweep(1, person_sds_2pl_nh, "/") %>%
  posterior_summary() %>%
	as_tibble() %>%
	rownames_to_column(var = "person") %>%
  mutate(person = as.numeric(person))

saveRDS(person_pars_2pl_nh, "results/person_pars_2pl_nh")
```

```{r}
person_pars_2pl_nh %>%
  arrange(Estimate) %>%
  mutate(id2 = seq_len(n())) %>%
	ggplot(aes(id2, Estimate, ymin = Q2.5, ymax = Q97.5)) +
	geom_pointrange(alpha = 0.7) +
	coord_flip() +
	labs(x = "Person Number (sorted after Estimate)") +
  theme(
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  )
```

```{r}
person_fit_2pl_nh <- fit_statistic(
  fit_2pl_nh, criterion = ll, group = person,
  nsamples = 1000
)

person_diff_2pl_nh <- person_fit_2pl_nh %>%
  group_by(person) %>%
  summarise(bp = mean(crit_diff > 0))

person_max_2pl_nh <- which.max(person_diff_2pl_nh$bp)
person_fit_2pl_nh %>% 
  filter(person == person_max_2pl_nh) %>%
  ggplot(aes(crit_diff)) +
  geom_histogram() +
  theme_hist()
```

```{r}
# extract item parameters
item_pars_2pl_nh <- brms::fixef(fit_2pl_nh, summary = FALSE)

# plot item parameters
# locations
sel_beta <- grepl("^beta", colnames(item_pars_2pl_nh))
beta <- item_pars_2pl_nh[, sel_beta] %>%
  posterior_summary() %>%
	as_tibble() %>%
	rownames_to_column()

# slopes
sel_alpha <- grepl("^logalpha", colnames(item_pars_2pl_nh))
alpha <- item_pars_2pl_nh[, sel_alpha] %>%
  exp() %>%
  sweep(1, person_sds_2pl_nh, "*") %>%
  posterior_summary() %>%
	as_tibble() %>%
	rownames_to_column()

item_pars_2pl_nh <- bind_rows(beta, alpha, .id = "nlpar") %>%
	rename(item = "rowname") %>%
	mutate(
	  item = as.numeric(item),
	  item = item - 0.1,
	  nlpar = factor(nlpar, labels = c("Location", "Slope")),
	  method = "MCMC-NH"  
	)

saveRDS(item_pars_2pl_nh, "results/item_pars_2pl_nh")
```

```{r}
item_fit_2pl_nh <- fit_statistic(
  fit_2pl_nh, criterion = ll, group = item,
  nsamples = 1000
)

item_fit_2pl_nh %>% 
  ggplot(aes(crit_diff)) +
  geom_histogram() +
  facet_wrap("item", scales = "free") +
  theme_hist()
```


### MAP

```{r, results="hide"}
parprior_2pl_map <- list(
  c(seq(1, 45, 4), "lnorm", 0, 2),
  c(seq(2, 46, 4), "norm", 0, 5)
)

fit_2pl_map <- mirt(
  spm_wide, model = 1, 
  itemtype = "2PL", SE = TRUE,
  technical = list(NCYCLES = 10000), 
  parprior = parprior_2pl_map
)
```

```{r}
person_pars_2pl_map <- fscores(fit_2pl_map, full.scores.SE = TRUE) %>%
  as.data.frame() %>%
  rename(Estimate = "F1", Est.Error = "SE_F1") %>%
  mutate(
    # person parameters are not scaled appropriately by mirt for the 2pl model
    scale = sd(Estimate) / sd(person_pars_2pl$Estimate),
    Estimate = Estimate / scale,
    Est.Error = Est.Error / scale,
    Q2.5 = Estimate - 1.96 * Est.Error, 
    Q97.5 = Estimate + 1.96 * Est.Error,
    person = seq_len(n())
  ) %>%
  select(-scale)

saveRDS(person_pars_2pl_map, "results/person_pars_2pl_map")
```

```{r, results="hide"}
item_pars_2pl_map <- coef(fit_2pl_map, as.data.frame = TRUE) %>%
  as.data.frame() %>%
  rownames_to_column() %>%
  filter(str_detect(rowname, "\\.(a1|d)$")) %>%
  mutate(
    item = str_match(rowname, "^X[[:digit:]]+") %>%
      str_remove("^X") %>%
      as.numeric(),
    item = item - 0.2,
    nlpar = str_match(rowname, "\\.[^\\.]+$") %>%
      str_remove("^\\.") %>%
      factor(
        levels = c("d", "a1"),
        labels = c("Location", "Slope")
      ),
    method = "MAP"
  ) %>%
  rename(Estimate = par, Q2.5 = CI_2.5, Q97.5 = CI_97.5) 

saveRDS(item_pars_2pl_map, "results/item_pars_2pl_map")
```

### ML

```{r, results="hide"}
fit_2pl_ml <- mirt(spm_wide, model = 1, itemtype = "2PL", SE = TRUE,
                   technical = list(NCYCLES = 10000))
```

```{r}
person_pars_2pl_ml <- fscores(fit_2pl_ml, full.scores.SE = TRUE) %>%
  as.data.frame() %>%
  rename(Estimate = "F1", Est.Error = "SE_F1") %>%
  mutate(
    # person parameters are not scaled appropriately by mirt for the 2pl model
    scale = sd(Estimate) / sd(person_pars_2pl$Estimate),
    Estimate = Estimate / scale,
    Est.Error = Est.Error / scale,
    Q2.5 = Estimate - 1.96 * Est.Error, 
    Q97.5 = Estimate + 1.96 * Est.Error,
    person = seq_len(n())
  ) %>%
  select(-scale)

saveRDS(person_pars_2pl_ml, "results/person_pars_2pl_ml")
```

```{r, results="hide"}
item_pars_2pl_ml <- coef(fit_2pl_ml, as.data.frame = TRUE) %>%
  as.data.frame() %>%
  rownames_to_column() %>%
  filter(str_detect(rowname, "\\.(a1|d)$")) %>%
  mutate(
    item = str_match(rowname, "^X[[:digit:]]+") %>%
      str_remove("^X") %>%
      as.numeric(),
    item = item - 0.3,
    nlpar = str_match(rowname, "\\.[^\\.]+$") %>%
      str_remove("^\\.") %>%
      factor(
        levels = c("d", "a1"),
        labels = c("Location", "Slope")
      ),
    method = "ML"
  ) %>%
  rename(Estimate = par, Q2.5 = CI_2.5, Q97.5 = CI_97.5) 

saveRDS(item_pars_2pl_ml, "results/item_pars_2pl_ml")
```


### Joint plots: Person parameters

```{r}
# plot comparison of person parameters
person_pars_2pl_all <- person_pars_2pl %>% 
  bind_cols(person_pars_2pl_ml) %>%
  mutate(
    UIW = Q97.5 - Q2.5,
    UIW1 = Q97.51 - Q2.51
  )

person_pars_2pl_all %>%
	ggplot(aes(Estimate, Estimate1, color = UIW)) +
	geom_point() +
	geom_abline() +
	scale_color_viridis_c() +
	lims(x = c(-3, 2), y = c(-3, 2)) +
	labs(
		x = "Estimate (MCMC-H)",
	  y = "Estimate (ML)",
		color = "UIW (MCMC-H)"
	)

person_pars_2pl_all %>%
	ggplot(aes(UIW, UIW1, color = Estimate)) +
	geom_point() +
	geom_abline() +
	scale_color_viridis_c() +
	lims(x = c(1, 2.4), y = c(1, 2.4)) +
	labs(
		x = "UIW (MCMC-H)",
	  y = "UIW (ML)",
		color = "Estimate (MCMC-H)"
	)
```

### Joint plots: Item parameters
  
```{r}
item_pars_2pl %>%
  bind_rows(
    item_pars_2pl_nh,
    item_pars_2pl_map, 
    item_pars_2pl_ml
  ) %>%
  mutate(
    method = factor(method, levels = c("MCMC-H", "MCMC-NH", "MAP", "ML"))
  ) %>%
	ggplot(aes(item, Estimate, ymin = Q2.5, ymax = Q97.5, color = method)) +
  scale_x_continuous(breaks = 1:12) +
	geom_pointrange() +
  facet_wrap("nlpar", scales = "free_x") +
  scale_color_viridis_d(name = "Method") +
	coord_flip() +
	labs(x = "Item Number")
```


## 3PL models with fixed guessing probability {.tabset}

### MCMC with hierarchical item priors

```{r}
formula_3pl_fixed <- bf(
  response2 ~ 0.125 + 0.875 * inv_logit(beta + exp(logalpha) * theta),
  nl = TRUE,
  theta ~ 0 + (1 | person),
  beta ~ 1 + (1 |i| item),
  logalpha ~ 1 + (1 |i| item),
  family = brmsfamily("bernoulli", link = "logit")
)

prior_3pl_fixed <- 
  prior("normal(0, 2)", class = "b", nlpar = "beta") +
  prior("normal(0, 1)", class = "b", nlpar = "logalpha") +
  prior("normal(0, 1)", class = "sd", group = "person", nlpar = "theta") + 
  prior("normal(0, 3)", class = "sd", group = "item", nlpar = "beta") +
  prior("normal(0, 1)", class = "sd", group = "item", nlpar = "logalpha")

fit_3pl_fixed <- brm(
  formula = formula_3pl_fixed,
  data = spm_long,
  prior = prior_3pl_fixed,
  chains = cores, iter = iter, warmup = warmup,
  cores = cores, control = control, inits = 0,
  file = "models/fit_3pl_fixed"
)

fit_3pl_fixed <- add_loo(fit_3pl_fixed)
```

```{r}
summary(fit_3pl_fixed)
plot(fit_3pl_fixed)
```

```{r}
person_pars_3pl_fixed <- 
  ranef(fit_3pl_fixed, summary = FALSE)$person[, , "theta_Intercept"] 

person_sds_3pl_fixed <- apply(person_pars_3pl_fixed, 1, sd)

person_pars_3pl_fixed <- person_pars_3pl_fixed %>%
  sweep(1, person_sds_3pl_fixed, "/") %>%
  posterior_summary() %>%
	as_tibble() %>%
	rownames_to_column(var = "person") %>%
  mutate(person = as.numeric(person))

saveRDS(person_pars_3pl_fixed, "results/person_pars_3pl_fixed") 
```

```{r}
person_pars_3pl_fixed %>%
  arrange(Estimate) %>%
  mutate(id2 = seq_len(n())) %>%
	ggplot(aes(id2, Estimate, ymin = Q2.5, ymax = Q97.5)) +
	geom_pointrange(alpha = 0.7) +
	coord_flip() +
	labs(x = "Person Number (sorted after Estimate)") +
  theme(
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  )
```

```{r}
person_fit_3pl_fixed <- fit_statistic(
  fit_3pl_fixed, criterion = ll, group = person,
  nsamples = 1000
)

person_diff_3pl_fixed <- person_fit_3pl_fixed %>%
  group_by(person) %>%
  summarise(bp = mean(crit_diff > 0))

person_max_3pl_fixed <- which.max(person_diff_3pl_fixed$bp)
person_fit_3pl_fixed %>% 
  filter(person == person_max_3pl_fixed) %>%
  ggplot(aes(crit_diff)) +
  geom_histogram() +
  theme_hist()
```

```{r, fig.height=4}
item_pars_3pl_fixed <- coef(fit_3pl_fixed, summary = FALSE)$item

# locations
beta <- item_pars_3pl_fixed[, , "beta_Intercept"] %>%
  posterior_summary() %>%
	as_tibble() %>%
	rownames_to_column() 

# slopes
alpha <- item_pars_3pl_fixed[, , "logalpha_Intercept"] %>%
  exp() %>%
  sweep(1, person_sds_3pl_fixed, "*") %>%
  posterior_summary() %>%
	as_tibble() %>%
	rownames_to_column()

item_pars_3pl_fixed <- bind_rows(beta, alpha, .id = "nlpar") %>%
	rename(item = "rowname") %>%
	mutate(item = as.numeric(item)) %>%
	mutate(
	  nlpar = factor(nlpar, labels = c("Location", "Slope")),
	  method = "MCMC-H"  
	)

saveRDS(item_pars_3pl_fixed, "results/item_pars_3pl_fixed")
```

```{r}
item_fit_3pl_fixed <- fit_statistic(
  fit_3pl_fixed, criterion = ll, group = item,
  nsamples = 1000
)

item_fit_3pl_fixed %>% 
  ggplot(aes(crit_diff)) +
  geom_histogram() +
  facet_wrap("item", scales = "free") +
  theme_hist()
```


### MCMC with non-hierarchical item priors

```{r}
formula_3pl_fixed_nh <- bf(
  response2 ~ 0.125 + 0.875 * inv_logit(beta + exp(logalpha) * theta),
  nl = TRUE,
  theta ~ 0 + (1 | person),
  beta ~ 0 + item,
  logalpha ~ 0 + item,
  family = brmsfamily("bernoulli", link = "logit")
)

prior_3pl_fixed_nh <- 
  prior("normal(0, 5)", class = "b", nlpar = "beta") +
  prior("normal(0, 1)", class = "b", nlpar = "logalpha") +
  prior("normal(0, 1)", class = "sd", group = "person", nlpar = "theta")

fit_3pl_fixed_nh <- brm(
  formula = formula_3pl_fixed_nh,
  data = spm_long,
  prior = prior_3pl_fixed_nh,
  chains = cores, iter = iter, warmup = warmup,
  cores = cores, control = control, inits = 0,
  file = "models/fit_3pl_fixed_nh"
)

fit_3pl_fixed_nh <- add_loo(fit_3pl_fixed_nh) 
```

```{r}
summary(fit_3pl_fixed_nh)
```

```{r}
person_pars_3pl_fixed_nh <- 
  ranef(fit_3pl_fixed_nh, summary = FALSE)$person[, , "theta_Intercept"] 

person_sds_3pl_fixed_nh <- apply(person_pars_3pl_fixed_nh, 1, sd)

person_pars_3pl_fixed_nh <- person_pars_3pl_fixed_nh %>%
  sweep(1, person_sds_3pl_fixed_nh, "/") %>%
  posterior_summary() %>%
	as_tibble() %>%
	rownames_to_column(var = "person") %>%
  mutate(person = as.numeric(person)) 

saveRDS(person_pars_3pl_fixed_nh, "results/person_pars_3pl_fixed_nh")
```

```{r}
person_pars_3pl_fixed_nh %>%
  arrange(Estimate) %>%
  mutate(id2 = seq_len(n())) %>%
	ggplot(aes(id2, Estimate, ymin = Q2.5, ymax = Q97.5)) +
	geom_pointrange(alpha = 0.7) +
	coord_flip() +
	labs(x = "Person Number (sorted after Estimate)") +
  theme(
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  ) 
```

```{r}
person_fit_3pl_fixed_nh <- fit_statistic(
  fit_3pl_fixed_nh, criterion = ll, group = person,
  nsamples = 1000
)

person_diff_3pl_fixed_nh <- person_fit_3pl_fixed_nh %>%
  group_by(person) %>%
  summarise(bp = mean(crit_diff > 0))

person_max_3pl_fixed_nh <- which.max(person_diff_3pl_fixed_nh$bp)
person_fit_3pl_fixed_nh %>% 
  filter(person == person_max_3pl_fixed_nh) %>%
  ggplot(aes(crit_diff)) +
  geom_histogram() +
  theme_hist() 
```

```{r}
# extract item parameters
item_pars_3pl_fixed_nh <- brms::fixef(fit_3pl_fixed_nh, summary = FALSE)

# plot item parameters
# locations
sel_beta <- grepl("^beta", colnames(item_pars_3pl_fixed_nh))
beta <- item_pars_3pl_fixed_nh[, sel_beta] %>%
  posterior_summary() %>%
	as_tibble() %>%
	rownames_to_column() 

# slopes
sel_alpha <- grepl("^logalpha", colnames(item_pars_3pl_fixed_nh))
alpha <- item_pars_3pl_fixed_nh[, sel_alpha] %>%
  exp() %>%
  sweep(1, person_sds_3pl_fixed_nh, "*") %>%
  posterior_summary() %>%
	as_tibble() %>%
	rownames_to_column()

item_pars_3pl_fixed_nh <- bind_rows(beta, alpha, .id = "nlpar") %>%
	rename(item = "rowname") %>%
	mutate(
	  item = as.numeric(item),
	  item = item - 0.1,
	  nlpar = factor(nlpar, labels = c("Location", "Slope")),
	  method = "MCMC-NH"  
	)

saveRDS(item_pars_3pl_fixed_nh, "results/item_pars_3pl_fixed_nh") 
```

```{r}
item_fit_3pl_fixed_nh <- fit_statistic(
  fit_3pl_fixed_nh, criterion = ll, group = item,
  nsamples = 1000
)
 
item_fit_3pl_fixed_nh %>% 
  ggplot(aes(crit_diff)) +
  geom_histogram() +
  facet_wrap("item", scales = "free") +
  theme_hist()
```


### MAP

```{r, results="hide"}
parprior_3pl_fixed_map <- list(
  c(seq(1, 45, 4), "lnorm", 0, 2),
  c(seq(2, 46, 4), "norm", 0, 5)
)

fit_3pl_fixed_map <- mirt(
  spm_wide, model = 1, 
  itemtype = "2PL", guess = 0.11, 
  SE = TRUE, technical = list(NCYCLES = 10000), 
  parprior = parprior_3pl_fixed_map
)
```

```{r}
person_pars_3pl_fixed_map <- fscores(fit_3pl_fixed_map, full.scores.SE = TRUE) %>%
  as.data.frame() %>%
  rename(Estimate = "F1", Est.Error = "SE_F1") %>%
  mutate(
    # person parameters are not scaled appropriately by mirt for the 3pl_fixed model
    scale = sd(Estimate) / sd(person_pars_3pl_fixed$Estimate),
    Estimate = Estimate / scale,
    Est.Error = Est.Error / scale,
    Q2.5 = Estimate - 1.96 * Est.Error, 
    Q97.5 = Estimate + 1.96 * Est.Error,
    person = seq_len(n())
  ) %>%
  select(-scale)

saveRDS(person_pars_3pl_fixed_map, "results/person_pars_3pl_fixed_map")
```

```{r, results="hide"}
item_pars_3pl_fixed_map <- coef(fit_3pl_fixed_map, as.data.frame = TRUE) %>%
  as.data.frame() %>%
  rownames_to_column() %>%
  filter(str_detect(rowname, "\\.(a1|d)$")) %>%
  mutate(
    item = str_match(rowname, "^X[[:digit:]]+") %>%
      str_remove("^X") %>%
      as.numeric(),
    item = item - 0.2,
    nlpar = str_match(rowname, "\\.[^\\.]+$") %>%
      str_remove("^\\.") %>%
      factor(
        levels = c("d", "a1"),
        labels = c("Location", "Slope")
      ),
    method = "MAP"
  ) %>%
  rename(Estimate = par, Q2.5 = CI_2.5, Q97.5 = CI_97.5) 

saveRDS(item_pars_3pl_fixed_map, "results/item_pars_3pl_fixed_map")
```

### ML

```{r, results="hide"}
fit_3pl_fixed_ml <- mirt(
  spm_wide, model = 1, itemtype = "2PL", guess = 0.11, 
  SE = TRUE, technical = list(NCYCLES = 10000))
```

```{r}
person_pars_3pl_fixed_ml <- fscores(fit_3pl_fixed_ml, full.scores.SE = TRUE) %>%
  as.data.frame() %>%
  rename(Estimate = "F1", Est.Error = "SE_F1") %>%
  mutate(
    # person parameters are not scaled appropriately by mirt for the 3pl_fixed model
    scale = sd(Estimate) / sd(person_pars_3pl_fixed$Estimate),
    Estimate = Estimate / scale,
    Est.Error = Est.Error / scale,
    Q2.5 = Estimate - 1.96 * Est.Error, 
    Q97.5 = Estimate + 1.96 * Est.Error,
    person = seq_len(n())
  ) %>%
  select(-scale)

saveRDS(person_pars_3pl_fixed_ml, "results/person_pars_3pl_fixed_ml")
```

```{r, results="hide"}
item_pars_3pl_fixed_ml <- coef(fit_3pl_fixed_ml, as.data.frame = TRUE) %>%
  as.data.frame() %>%
  rownames_to_column() %>%
  filter(str_detect(rowname, "\\.(a1|d)$")) %>%
  mutate(
    item = str_match(rowname, "^X[[:digit:]]+") %>%
      str_remove("^X") %>%
      as.numeric(),
    item = item - 0.3,
    nlpar = str_match(rowname, "\\.[^\\.]+$") %>%
      str_remove("^\\.") %>%
      factor(
        levels = c("d", "a1"),
        labels = c("Location", "Slope")
      ),
    method = "ML"
  ) %>%
  rename(Estimate = par, Q2.5 = CI_2.5, Q97.5 = CI_97.5) 

saveRDS(item_pars_3pl_fixed_ml, "results/item_pars_3pl_fixed_ml")
```


### Joint plots: Person parameters

```{r}
# plot comparison of person parameters
person_pars_3pl_fixed_all <- person_pars_3pl_fixed %>% 
  bind_cols(person_pars_3pl_fixed_ml) %>%
  mutate(
    UIW = Q97.5 - Q2.5,
    UIW1 = Q97.51 - Q2.51
  )

person_pars_3pl_fixed_all %>%
	ggplot(aes(Estimate, Estimate1, color = UIW)) +
	geom_point() +
	geom_abline() +
	scale_color_viridis_c() +
	lims(x = c(-3, 2), y = c(-3, 2)) +
	labs(
		x = "Estimate (MCMC-H)",
	  y = "Estimate (ML)",
		color = "UIW (MCMC-H)"
	)

person_pars_3pl_fixed_all %>%
	ggplot(aes(UIW, UIW1, color = Estimate)) +
	geom_point() +
	geom_abline() +
	scale_color_viridis_c() +
	lims(x = c(1, 2.4), y = c(1, 2.4)) +
	labs(
		x = "UIW (MCMC-H)",
	  y = "UIW (ML)",
		color = "Estimate (MCMC-H)"
	)
```

### Joint plots: Item parameters
  
```{r}
item_pars_3pl_fixed %>%
  bind_rows(
    item_pars_3pl_fixed_nh,
    item_pars_3pl_fixed_map, 
    item_pars_3pl_fixed_ml
  ) %>%
  mutate(
    method = factor(method, levels = c("MCMC-H", "MCMC-NH", "MAP", "ML"))
  ) %>%
	ggplot(aes(item, Estimate, ymin = Q2.5, ymax = Q97.5, color = method)) +
  scale_x_continuous(breaks = 1:12) +
	geom_pointrange() +
  facet_wrap("nlpar", scales = "free_x") +
  scale_color_viridis_d(name = "Method") +
	coord_flip() +
	labs(x = "Item Number")
```


## 3PL Models {.tabset}

### MCMC with hierarchical item priors

```{r}
formula_3pl <- bf(
  response2 ~ gamma + (1 - gamma) * 
    inv_logit(beta + exp(logalpha) * theta),
  nl = TRUE,
  theta ~ 0 + (1 | person),
  beta ~ 1 + (1 |i| item),
  logalpha ~ 1 + (1 |i| item),
  logitgamma ~ 1 + (1 |i| item),
  nlf(gamma ~ inv_logit(logitgamma)),
  family = brmsfamily("bernoulli", link = "identity")
)

prior_3pl <- 
  prior("normal(0, 2)", class = "b", nlpar = "beta") +
  prior("normal(0, 1)", class = "b", nlpar = "logalpha") +
  prior("normal(-2, 0.5)", class = "b", nlpar = "logitgamma") +
  prior("normal(0, 1)", class = "sd", group = "person", nlpar = "theta") + 
  prior("normal(0, 3)", class = "sd", group = "item", nlpar = "beta") +
  prior("normal(0, 1)", class = "sd", group = "item", nlpar = "logalpha") +
  prior("normal(0, 1)", class = "sd", group = "item", nlpar = "logitgamma")

fit_3pl <- brm(
  formula = formula_3pl,
  data = spm_long,
  prior = prior_3pl,
  chains = cores, iter = iter, warmup = warmup,
  cores = cores, control = control, inits = 0,
  file = "models/fit_3pl"
)

fit_3pl <- add_loo(fit_3pl)
```

```{r}
summary(fit_3pl)
plot(fit_3pl)
```

```{r}
person_pars_3pl <- 
  ranef(fit_3pl, summary = FALSE)$person[, , "theta_Intercept"] 

person_sds_3pl <- apply(person_pars_3pl, 1, sd)

person_pars_3pl <- person_pars_3pl %>%
  sweep(1, person_sds_3pl, "/") %>%
  posterior_summary() %>%
	as_tibble() %>%
	rownames_to_column(var = "person") %>%
  mutate(person = as.numeric(person))

saveRDS(person_pars_3pl, "results/person_pars_3pl")
```

```{r}
person_pars_3pl %>%
  arrange(Estimate) %>%
  mutate(id2 = seq_len(n())) %>%
	ggplot(aes(id2, Estimate, ymin = Q2.5, ymax = Q97.5)) +
	geom_pointrange(alpha = 0.7) +
	coord_flip() +
	labs(x = "Person Number (sorted after Estimate)") +
  theme(
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  )
```

```{r}
person_fit_3pl <- fit_statistic(
  fit_3pl, criterion = ll, group = person,
  nsamples = 1000
)

person_diff_3pl <- person_fit_3pl %>%
  group_by(person) %>%
  summarise(bp = mean(crit_diff > 0))

person_max_3pl <- which.max(person_diff_3pl$bp)
person_fit_3pl %>% 
  filter(person == person_max_3pl) %>%
  ggplot(aes(crit_diff)) +
  geom_histogram() +
  theme_hist()
```

```{r}
# extract item parameters
item_pars_3pl <- coef(fit_3pl, summary = FALSE)$item

# plot item parameters
# locations
beta <- item_pars_3pl[, , "beta_Intercept"] %>%
  posterior_summary() %>%
	as_tibble() %>%
	rownames_to_column()

# slopes
alpha <- item_pars_3pl[, , "logalpha_Intercept"] %>%
	exp() %>%
  sweep(1, person_sds_3pl, "*") %>%
  posterior_summary() %>%
	as_tibble() %>%
	rownames_to_column()

# guessing parameters
gamma <- item_pars_3pl[, , "logitgamma_Intercept"] %>%
	inv_logit_scaled() %>%
  posterior_summary() %>%
	as_tibble() %>%
	rownames_to_column()

item_pars_3pl <- bind_rows(beta, alpha, gamma, .id = "nlpar") %>%
	rename(item = "rowname") %>%
	mutate(item = as.numeric(item)) %>%
	mutate(
	  nlpar = factor(nlpar, labels = c("Location", "Slope", "Guessing")),
	  method = "MCMC-H"  
	)

saveRDS(item_pars_3pl, "results/item_pars_3pl")
```

```{r}
item_fit_3pl <- fit_statistic(
  fit_3pl, criterion = ll, group = item,
  nsamples = 1000
)

item_fit_3pl %>% 
  ggplot(aes(crit_diff)) +
  geom_histogram() +
  facet_wrap("item", scales = "free") +
  theme_hist()
```


### MCMC with non-hierarchical item priors

```{r}
formula_3pl_nh <- bf(
  response2 ~ gamma + (1 - gamma) * 
    inv_logit(beta + exp(logalpha) * theta),
  nl = TRUE,
  theta ~ 0 + (1 | person),
  beta ~ 0 + item,
  logalpha ~ 0 + item,
  logitgamma ~ 0 + item,
  nlf(gamma ~ inv_logit(logitgamma)),
  family = brmsfamily("bernoulli", link = "identity")
)

prior_3pl_nh <- 
  prior("normal(0, 5)", class = "b", nlpar = "beta") +
  prior("normal(0, 1)", class = "b", nlpar = "logalpha") +
  prior("normal(-2, 1)", class = "b", nlpar = "logitgamma") +
  prior("normal(0, 1)", class = "sd", group = "person", nlpar = "theta")

fit_3pl_nh <- brm(
  formula = formula_3pl_nh,
  data = spm_long,
  prior = prior_3pl_nh,
  chains = cores, iter = iter, warmup = warmup,
  cores = cores, control = control, inits = 0,
  file = "models/fit_3pl_nh"
)

fit_3pl_nh <- add_loo(fit_3pl_nh)
```

```{r}
summary(fit_3pl_nh)
plot(fit_3pl_nh)
```

```{r}
person_pars_3pl_nh <- 
  ranef(fit_3pl_nh, summary = FALSE)$person[, , "theta_Intercept"] 

person_sds_3pl_nh <- apply(person_pars_3pl_nh, 1, sd)

person_pars_3pl_nh <- person_pars_3pl_nh %>%
  sweep(1, person_sds_3pl_nh, "/") %>%
  posterior_summary() %>%
	as_tibble() %>%
	rownames_to_column(var = "person") %>%
  mutate(person = as.numeric(person))

saveRDS(person_pars_3pl_nh, "results/person_pars_3pl_nh")
```

```{r}
person_pars_3pl_nh %>%
  arrange(Estimate) %>%
  mutate(id2 = seq_len(n())) %>%
	ggplot(aes(id2, Estimate, ymin = Q2.5, ymax = Q97.5)) +
	geom_pointrange(alpha = 0.7) +
	coord_flip() +
	labs(x = "Person Number (sorted after Estimate)") +
  theme(
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  )
```

```{r}
person_fit_3pl_nh <- fit_statistic(
  fit_3pl_nh, criterion = ll, group = person,
  nsamples = 1000
)

person_diff_3pl_nh <- person_fit_3pl_nh %>%
  group_by(person) %>%
  summarise(bp = mean(crit_diff > 0))

person_max_3pl_nh <- which.max(person_diff_3pl_nh$bp)
person_fit_3pl_nh %>% 
  filter(person == person_max_3pl_nh) %>%
  ggplot(aes(crit_diff)) +
  geom_histogram() +
  theme_hist()
```

```{r}
# extract item parameters
item_pars_3pl_nh <- brms::fixef(fit_3pl_nh, summary = FALSE)

# plot item parameters
# locations
sel_beta <- grepl("^beta", colnames(item_pars_3pl_nh))
beta <- item_pars_3pl_nh[, sel_beta] %>%
  posterior_summary() %>%
	as_tibble() %>%
	rownames_to_column()

# slopes
sel_alpha <- grepl("^logalpha", colnames(item_pars_3pl_nh))
alpha <- item_pars_3pl_nh[, sel_alpha] %>%
	exp() %>%
  sweep(1, person_sds_3pl_nh, "*") %>%
  posterior_summary() %>%
	as_tibble() %>%
	rownames_to_column()

# guessing parameters
sel_gamma <- grepl("^logitgamma", colnames(item_pars_3pl_nh))
gamma <- item_pars_3pl_nh[, sel_gamma] %>%
	inv_logit_scaled() %>%
  posterior_summary() %>%
	as_tibble() %>%
	rownames_to_column()

item_pars_3pl_nh <- bind_rows(beta, alpha, gamma, .id = "nlpar") %>%
	rename(item = "rowname") %>%
	mutate(
	  item = as.numeric(item),
	  item = item - 0.1,
	  nlpar = factor(nlpar, labels = c("Location", "Slope", "Guessing")),
	  method = "MCMC-NH"  
	)

saveRDS(item_pars_3pl_nh, "results/item_pars_3pl_nh")
```

```{r}
item_fit_3pl_nh <- fit_statistic(
  fit_3pl_nh, criterion = ll, group = item,
  nsamples = 1000
)

item_fit_3pl_nh %>% 
  ggplot(aes(crit_diff)) +
  geom_histogram() +
  facet_wrap("item", scales = "free") +
  theme_hist()
```


### MAP

```{r, results="hide"}
parprior_3pl_map <- list(
  c(seq(1, 45, 4), "lnorm", 0, 2),
  c(seq(2, 46, 4), "norm", 0, 5),
  c(seq(3, 47, 4), "norm", -2, 1)
)

fit_3pl_map <- mirt(
  spm_wide, model = 1, 
  itemtype = "3PL", SE = TRUE,
  technical = list(NCYCLES = 10000), 
  parprior = parprior_3pl_map
)
```

```{r}
person_pars_3pl_map <- fscores(fit_3pl_map, full.scores.SE = TRUE) %>%
  as.data.frame() %>%
  rename(Estimate = "F1", Est.Error = "SE_F1") %>%
  mutate(
    # person parameters are not scaled appropriately by mirt for the 3pl model
    scale = sd(Estimate) / sd(person_pars_3pl$Estimate),
    Estimate = Estimate / scale,
    Est.Error = Est.Error / scale,
    Q2.5 = Estimate - 1.96 * Est.Error, 
    Q97.5 = Estimate + 1.96 * Est.Error,
    person = seq_len(n())
  ) %>%
  select(-scale)

saveRDS(person_pars_3pl_map, "results/person_pars_3pl_map")
```

```{r, results="hide"}
item_pars_3pl_map <- coef(fit_3pl_map, as.data.frame = TRUE) %>%
  as.data.frame() %>%
  rownames_to_column() %>%
  filter(str_detect(rowname, "\\.(a1|d|g)$")) %>%
  mutate(
    item = str_match(rowname, "^X[[:digit:]]+") %>%
      str_remove("^X") %>%
      as.numeric(),
    item = item - 0.2,
    nlpar = str_match(rowname, "\\.[^\\.]+$") %>%
      str_remove("^\\.") %>%
      factor(
        levels = c("d", "a1", "g"),
        labels = c("Location", "Slope", "Guessing")
      ),
    method = "MAP"
  ) %>%
  rename(Estimate = par, Q2.5 = CI_2.5, Q97.5 = CI_97.5) 

saveRDS(item_pars_3pl_map, "results/item_pars_3pl_map")
```

### ML

```{r, results="hide"}
fit_3pl_ml <- mirt(spm_wide, model = 1, itemtype = "3PL", SE = TRUE,
                   technical = list(NCYCLES = 10000))
```

```{r}
person_pars_3pl_ml <- fscores(fit_3pl_ml, full.scores.SE = TRUE) %>%
  as.data.frame() %>%
  rename(Estimate = "F1", Est.Error = "SE_F1") %>%
  mutate(
    # person parameters are not scaled appropriately by mirt for the 3pl model
    scale = sd(Estimate) / sd(person_pars_3pl$Estimate),
    Estimate = Estimate / scale,
    Est.Error = Est.Error / scale,
    Q2.5 = Estimate - 1.96 * Est.Error, 
    Q97.5 = Estimate + 1.96 * Est.Error,
    person = seq_len(n())
  ) %>%
  select(-scale)

saveRDS(person_pars_3pl_ml, "results/person_pars_3pl_ml")
```

```{r, results="hide"}
item_pars_3pl_ml <- coef(fit_3pl_ml, as.data.frame = TRUE) %>%
  as.data.frame() %>%
  rownames_to_column() %>%
  filter(str_detect(rowname, "\\.(a1|d|g)$")) %>%
  mutate(
    item = str_match(rowname, "^X[[:digit:]]+") %>%
      str_remove("^X") %>%
      as.numeric(),
    item = item - 0.3,
    nlpar = str_match(rowname, "\\.[^\\.]+$") %>%
      str_remove("^\\.") %>%
      factor(
        levels = c("d", "a1", "g"),
        labels = c("Location", "Slope", "Guessing")
      ),
    method = "ML"
  ) %>%
  rename(Estimate = par, Q2.5 = CI_2.5, Q97.5 = CI_97.5)

saveRDS(item_pars_3pl_ml, "results/item_pars_3pl_ml")
```


### Joint plots: Person parameters

```{r}
# plot comparison of person parameters
person_pars_3pl_all <- person_pars_3pl %>% 
  bind_cols(person_pars_3pl_ml) %>%
  mutate(
    UIW = Q97.5 - Q2.5,
    UIW1 = Q97.51 - Q2.51
  )

person_pars_3pl_all %>%
	ggplot(aes(Estimate, Estimate1, color = UIW)) +
	geom_point() +
	geom_abline() +
	scale_color_viridis_c() +
	lims(x = c(-3, 2), y = c(-3, 2)) +
	labs(
		x = "Estimate (MCMC-H)",
	  y = "Estimate (ML)",
		color = "UIW (MCMC-H)"
	)

person_pars_3pl_all %>%
	ggplot(aes(UIW, UIW1, color = Estimate)) +
	geom_point() +
	geom_abline() +
	scale_color_viridis_c() +
	lims(x = c(1, 2.4), y = c(1, 2.4)) +
	labs(
		x = "UIW (MCMC-H)",
	  y = "UIW (ML)",
		color = "Estimate (MCMC-H)"
	)
```

### Joint plots: Item parameters
  
```{r}
item_pars_3pl %>%
  bind_rows(
    item_pars_3pl_nh,
    item_pars_3pl_map, 
    item_pars_3pl_ml
  ) %>%
  mutate(
    method = factor(method, levels = c("MCMC-H", "MCMC-NH", "MAP", "ML"))
  ) %>%
	ggplot(aes(item, Estimate, ymin = Q2.5, ymax = Q97.5, color = method)) +
  scale_x_continuous(breaks = 1:12) +
	geom_pointrange() +
  facet_wrap("nlpar", scales = "free_x") +
  scale_color_viridis_d(name = "Method") +
	coord_flip() +
	labs(x = "Item Number")
```


## 4PL Models {.tabset}

### MCMC with hierarchical item priors

```{r}
formula_4pl <- bf(
  response2 ~ gamma + (1 - gamma - psi) * 
    inv_logit(beta + exp(logalpha) * theta),
  nl = TRUE,
  theta ~ 0 + (1 | person),
  beta ~ 1 + (1 |i| item),
  logalpha ~ 1 + (1 |i| item),
  logitgamma ~ 1 + (1 |i| item),
  nlf(gamma ~ inv_logit(logitgamma)),
  logitpsi ~ 1 + (1 |i| item),
  nlf(psi ~ inv_logit(logitpsi)),
  family = brmsfamily("bernoulli", link = "identity")
)

prior_4pl <- 
  prior("normal(0, 2)", class = "b", nlpar = "beta") +
  prior("normal(0, 1)", class = "b", nlpar = "logalpha") +
  prior("normal(-2, 0.5)", class = "b", nlpar = "logitgamma") +
  prior("normal(-2, 0.5)", class = "b", nlpar = "logitpsi") +
  prior("normal(0, 1)", class = "sd", group = "person", nlpar = "theta") + 
  prior("normal(0, 3)", class = "sd", group = "item", nlpar = "beta") +
  prior("normal(0, 1)", class = "sd", group = "item", nlpar = "logalpha") +
  prior("normal(0, 1)", class = "sd", group = "item", nlpar = "logitgamma") +
  prior("normal(0, 1)", class = "sd", group = "item", nlpar = "logitpsi")

fit_4pl <- brm(
  formula = formula_4pl,
  data = spm_long,
  prior = prior_4pl,
  chains = cores, iter = iter, warmup = warmup,
  cores = cores, control = control, inits = 0,
  file = "models/fit_4pl"
)

fit_4pl <- add_loo(fit_4pl) 
```

```{r}
summary(fit_4pl)
plot(fit_4pl)
```

```{r}
person_pars_4pl <- 
  ranef(fit_4pl, summary = FALSE)$person[, , "theta_Intercept"] 

person_sds_4pl <- apply(person_pars_4pl, 1, sd)

person_pars_4pl <- person_pars_4pl %>%
  sweep(1, person_sds_4pl, "/") %>%
  posterior_summary() %>%
	as_tibble() %>%
	rownames_to_column(var = "person") %>%
  mutate(person = as.numeric(person))

saveRDS(person_pars_4pl, "results/person_pars_4pl")
```

```{r}
person_pars_4pl %>%
  arrange(Estimate) %>%
  mutate(id2 = seq_len(n())) %>%
	ggplot(aes(id2, Estimate, ymin = Q2.5, ymax = Q97.5)) +
	geom_pointrange(alpha = 0.7) +
	coord_flip() +
	labs(x = "Person Number (sorted after Estimate)") +
  theme(
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  )
```

```{r}
person_fit_4pl <- fit_statistic(
  fit_4pl, criterion = ll, group = person,
  nsamples = 1000
)

person_diff_4pl <- person_fit_4pl %>%
  group_by(person) %>%
  summarise(bp = mean(crit_diff > 0))

person_max_4pl <- which.max(person_diff_4pl$bp)
person_fit_4pl %>% 
  filter(person == person_max_4pl) %>%
  ggplot(aes(crit_diff)) +
  geom_histogram() +
  theme_hist()
```

```{r}
# extract item parameters
item_pars_4pl <- coef(fit_4pl, summary = FALSE)$item

# locations
beta <- item_pars_4pl[, , "beta_Intercept"] %>%
  posterior_summary() %>%
	as_tibble() %>%
	rownames_to_column()

# slopes
alpha <- item_pars_4pl[, , "logalpha_Intercept"] %>%
	exp() %>%
  sweep(1, person_sds_4pl, "*") %>%
  posterior_summary() %>%
	as_tibble() %>%
	rownames_to_column()

# guessing parameters
gamma <- item_pars_4pl[, , "logitgamma_Intercept"] %>%
	inv_logit_scaled() %>%
  posterior_summary() %>%
	as_tibble() %>%
	rownames_to_column()

# lapse parameters
psi <- item_pars_4pl[, , "logitpsi_Intercept"] %>%
	inv_logit_scaled() %>%
  posterior_summary() %>%
	as_tibble() %>%
	rownames_to_column()

item_pars_4pl <- bind_rows(beta, alpha, gamma, psi, .id = "nlpar") %>%
	rename(item = "rowname") %>%
	mutate(item = as.numeric(item)) %>%
	mutate(
	  nlpar = nlpar %>% 
	    factor(labels = c("Location", "Slope", "Guessing", "Lapse")),
	  method = "MCMC-H"
	)

saveRDS(item_pars_4pl, "results/item_pars_4pl")
```

```{r}
item_fit_4pl <- fit_statistic(
  fit_4pl, criterion = ll, group = item,
  nsamples = 1000
)

item_fit_4pl %>% 
  ggplot(aes(crit_diff)) +
  geom_histogram() +
  facet_wrap("item", scales = "free") +
  theme_hist()
```


### MCMC with non-hierarchical item priors

```{r}
formula_4pl_nh <- bf(
  response2 ~ gamma + (1 - gamma - psi) * 
    inv_logit(beta + exp(logalpha) * theta),
  nl = TRUE,
  theta ~ 0 + (1 | person),
  beta ~ 0 + item,
  logalpha ~ 0 + item,
  logitgamma ~ 0 + item,
  nlf(gamma ~ inv_logit(logitgamma)),
  logitpsi ~ 0 + item,
  nlf(psi ~ inv_logit(logitpsi)),
  family = brmsfamily("bernoulli", link = "identity")
)

prior_4pl_nh <- 
  prior("normal(0, 5)", class = "b", nlpar = "beta") +
  prior("normal(0, 1)", class = "b", nlpar = "logalpha") +
  prior("normal(-2, 1)", class = "b", nlpar = "logitgamma") +
  prior("normal(-2, 1)", class = "b", nlpar = "logitpsi") +
  prior("normal(0, 1)", class = "sd", group = "person", nlpar = "theta")

fit_4pl_nh <- brm(
  formula = formula_4pl_nh,
  data = spm_long,
  prior = prior_4pl_nh,
  chains = cores, iter = iter, warmup = warmup,
  cores = cores, control = control, inits = 0,
  file = "models/fit_4pl_nh"
)

fit_4pl_nh <- add_loo(fit_4pl_nh) 
```

```{r}
summary(fit_4pl_nh)
plot(fit_4pl_nh)
```

```{r}
person_pars_4pl_nh <- 
  ranef(fit_4pl_nh, summary = FALSE)$person[, , "theta_Intercept"] 

person_sds_4pl_nh <- apply(person_pars_4pl_nh, 1, sd)

person_pars_4pl_nh <- person_pars_4pl_nh %>%
  sweep(1, person_sds_4pl_nh, "/") %>%
  posterior_summary() %>%
	as_tibble() %>%
	rownames_to_column(var = "person") %>%
  mutate(person = as.numeric(person))

saveRDS(person_pars_4pl_nh, "results/person_pars_4pl_nh")
```

```{r}
person_pars_4pl_nh %>%
  arrange(Estimate) %>%
  mutate(id2 = seq_len(n())) %>%
	ggplot(aes(id2, Estimate, ymin = Q2.5, ymax = Q97.5)) +
	geom_pointrange(alpha = 0.7) +
	coord_flip() +
	labs(x = "Person Number (sorted after Estimate)") +
  theme(
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  )
```

```{r}
person_fit_4pl_nh <- fit_statistic(
  fit_4pl_nh, criterion = ll, group = person,
  nsamples = 1000
)

person_diff_4pl_nh <- person_fit_4pl_nh %>%
  group_by(person) %>%
  summarise(bp = mean(crit_diff > 0))

person_max_4pl_nh <- which.max(person_diff_4pl_nh$bp)
person_fit_4pl_nh %>% 
  filter(person == person_max_4pl_nh) %>%
  ggplot(aes(crit_diff)) +
  geom_histogram() +
  theme_hist()
```

```{r}
# extract item parameters
item_pars_4pl_nh <- brms::fixef(fit_4pl_nh, summary = FALSE)

# plot item parameters
# locations
sel_beta <- grepl("^beta", colnames(item_pars_4pl_nh))
beta <- item_pars_4pl_nh[, sel_beta] %>%
  posterior_summary() %>%
	as_tibble() %>%
	rownames_to_column()

# slopes
sel_alpha <- grepl("^logalpha", colnames(item_pars_4pl_nh))
alpha <- item_pars_4pl_nh[, sel_alpha] %>%
	exp() %>%
  sweep(1, person_sds_4pl_nh, "*") %>%
  posterior_summary() %>%
	as_tibble() %>%
	rownames_to_column()

# guessing parameters
sel_gamma <- grepl("^logitgamma", colnames(item_pars_4pl_nh))
gamma <- item_pars_4pl_nh[, sel_gamma] %>%
	inv_logit_scaled() %>%
  posterior_summary() %>%
	as_tibble() %>%
	rownames_to_column()

# lapse parameters
sel_psi <- grepl("^logitpsi", colnames(item_pars_4pl_nh))
psi <- item_pars_4pl_nh[, sel_psi] %>%
	inv_logit_scaled() %>%
  posterior_summary() %>%
	as_tibble() %>%
	rownames_to_column()

item_pars_4pl_nh <- bind_rows(beta, alpha, gamma, psi, .id = "nlpar") %>%
	rename(item = "rowname") %>%
	mutate(
	  item = as.numeric(item),
	  item = item - 0.1,
	  nlpar = nlpar %>% 
	    factor(labels = c("Location", "Slope", "Guessing", "Lapse")),
	  method = "MCMC-NH"
	) 

saveRDS(item_pars_4pl_nh, "results/item_pars_4pl_nh")
```

```{r}
item_fit_4pl_nh <- fit_statistic(
  fit_4pl_nh, criterion = ll, group = item,
  nsamples = 1000
)

item_fit_4pl_nh %>% 
  ggplot(aes(crit_diff)) +
  geom_histogram() +
  facet_wrap("item", scales = "free") +
  theme_hist()
```


### MAP

```{r, results="hide"}
parprior_4pl_map <- list(
  c(seq(1, 45, 4), "lnorm", 0, 2),
  c(seq(2, 46, 4), "norm", 0, 5),
  c(seq(3, 47, 4), "norm", -2, 1),
  c(seq(4, 48, 4), "norm", 2, 1)
)

fit_4pl_map <- mirt(
  spm_wide, model = 1, 
  itemtype = "4PL", SE = TRUE,
  technical = list(NCYCLES = 10000), 
  parprior = parprior_4pl_map
)
```

```{r}
person_pars_4pl_map <- fscores(fit_4pl_map, full.scores.SE = TRUE) %>%
  as.data.frame() %>%
  rename(Estimate = "F1", Est.Error = "SE_F1") %>%
  mutate(
    # person parameters are not scaled appropriately by mirt for the 4pl model
    scale = sd(Estimate) / sd(person_pars_4pl$Estimate),
    Estimate = Estimate / scale,
    Est.Error = Est.Error / scale,
    Q2.5 = Estimate - 1.96 * Est.Error, 
    Q97.5 = Estimate + 1.96 * Est.Error,
    person = seq_len(n())
  ) %>%
  select(-scale)

saveRDS(person_pars_4pl_map, "results/person_pars_4pl_map")
```

```{r, results="hide"}
item_pars_4pl_map <- coef(fit_4pl_map, as.data.frame = TRUE) %>%
  as.data.frame() %>%
  rownames_to_column() %>%
  filter(str_detect(rowname, "\\.(a1|d|g|u)$")) %>%
  mutate(
    item = str_match(rowname, "^X[[:digit:]]+") %>%
      str_remove("^X") %>%
      as.numeric(),
    item = item - 0.2,
    nlpar = str_match(rowname, "\\.[^\\.]+$") %>%
      str_remove("^\\.") %>%
      factor(
        levels = c("d", "a1", "g", "u"),
        labels = c("Location", "Slope", "Guessing", "Lapse")
      ),
    # invert guessing parameter
    par = ifelse(nlpar == "Lapse", 1 - par, par),
    CI_2.5 = ifelse(nlpar == "Lapse", 1 - CI_2.5, CI_2.5),
    CI_97.5 = ifelse(nlpar == "Lapse", 1 - CI_97.5, CI_97.5),
    method = "MAP"
  ) %>%
  rename(Estimate = par, Q2.5 = CI_2.5, Q97.5 = CI_97.5)

saveRDS(item_pars_4pl_map, "results/item_pars_4pl_map")
```

### ML

```{r, results="hide"}
fit_4pl_ml <- mirt(
  spm_wide, model = 1, itemtype = "4PL", SE = TRUE,
  technical = list(NCYCLES = 10000)
)
```

```{r}
person_pars_4pl_ml <- fscores(fit_4pl_ml, full.scores.SE = TRUE) %>%
  as.data.frame() %>%
  rename(Estimate = "F1", Est.Error = "SE_F1") %>%
  mutate(
    # person parameters are not scaled appropriately by mirt for the 4pl model
    scale = sd(Estimate) / sd(person_pars_4pl$Estimate),
    Estimate = Estimate / scale,
    Est.Error = Est.Error / scale,
    Q2.5 = Estimate - 1.96 * Est.Error, 
    Q97.5 = Estimate + 1.96 * Est.Error,
    person = seq_len(n())
  ) %>%
  select(-scale)

saveRDS(person_pars_4pl_ml, "results/person_pars_4pl_ml")
```

```{r, results="hide"}
item_pars_4pl_ml <- coef(fit_4pl_ml, as.data.frame = TRUE) %>%
  as.data.frame() %>%
  rownames_to_column() %>%
  filter(str_detect(rowname, "\\.(a1|d|g|u)$")) %>%
  mutate(
    item = str_match(rowname, "^X[[:digit:]]+") %>%
      str_remove("^X") %>%
      as.numeric(),
    item = item - 0.3,
    nlpar = str_match(rowname, "\\.[^\\.]+$") %>%
      str_remove("^\\.") %>%
      factor(
        levels = c("d", "a1", "g", "u"),
        labels = c("Location", "Slope", "Guessing", "Lapse")
      ),
    # invert guessing parameter
    par = ifelse(nlpar == "Lapse", 1 - par, par),
    # inverting the information matrix fails
    CI_2.5 = par,
    CI_97.5 = par,
    method = "ML"
  ) %>%
  rename(Estimate = par, Q2.5 = CI_2.5, Q97.5 = CI_97.5)

saveRDS(item_pars_4pl_ml, "results/item_pars_4pl_ml")
```


### Joint plots: Person parameters

```{r}
# plot comparison of person parameters
person_pars_4pl_all <- person_pars_4pl %>% 
  bind_cols(person_pars_4pl_ml) %>%
  mutate(
    UIW = Q97.5 - Q2.5,
    UIW1 = Q97.51 - Q2.51
  )

person_pars_4pl_all %>%
	ggplot(aes(Estimate, Estimate1, color = UIW)) +
	geom_point() +
	geom_abline() +
	scale_color_viridis_c() +
	lims(x = c(-3, 2), y = c(-3, 2)) +
	labs(
		x = "Estimate (MCMC-H)",
	  y = "Estimate (ML)",
		color = "UIW (MCMC-H)"
	)

person_pars_4pl_all %>%
	ggplot(aes(UIW, UIW1, color = Estimate)) +
	geom_point() +
	geom_abline() +
	scale_color_viridis_c() +
	lims(x = c(1, 2.4), y = c(1, 2.4)) +
	labs(
		x = "UIW (MCMC-H)",
	  y = "UIW (ML)",
		color = "Estimate (MCMC-H)"
	)
```

### Joint plots: Item parameters
  
```{r}
item_pars_4pl %>%
  bind_rows(
    item_pars_4pl_nh,
    item_pars_4pl_map, 
    item_pars_4pl_ml
  ) %>%
  mutate(
    method = factor(method, levels = c("MCMC-H", "MCMC-NH", "MAP", "ML"))
  ) %>%
	ggplot(aes(item, Estimate, ymin = Q2.5, ymax = Q97.5, color = method)) +
  scale_x_continuous(breaks = 1:12) +
	geom_pointrange() +
  facet_wrap("nlpar", scales = "free_x") +
  scale_color_viridis_d(name = "Method") +
	coord_flip() +
	labs(x = "Item Number")
```



## Comparison of models

### Out-Of-Sample Predictive Performance

```{r}
loo_compare(fit_1pl, fit_2pl, fit_3pl_fixed, fit_3pl, fit_4pl)
```

```{r}
loo_compare(fit_1pl_nh, fit_2pl_nh, fit_3pl_fixed_nh, fit_3pl_nh, fit_4pl_nh) 
```

### Person Parameters

```{r}
person_pars_all <- bind_rows(
  "1PL: brms" = person_pars_1pl,
  "1PL: mirt" = person_pars_1pl_ml, 
  "2PL: brms" = person_pars_2pl, 
  "2PL: mirt" = person_pars_2pl_ml, 
  "3PL: brms" = person_pars_3pl,
  "3PL: mirt" = person_pars_3pl_ml,
  "4PL: brms" = person_pars_4pl, 
  "4PL: mirt" = person_pars_4pl_ml, 
  .id = "model"
) %>%
  mutate(model = factor(model, levels = unique(model)))
```

```{r, fig.height=12, fig.width=12}
person_pars_all %>%
  select(Estimate, model, person) %>% 
  spread(model, Estimate) %>%
  select(-person) %>%
  GGally::ggpairs(aes(alpha = 0.4))
```

```{r, fig.height=12, fig.width=12}
person_pars_all %>%
  select(Est.Error, model, person) %>% 
  spread(model, Est.Error) %>%
  select(-person) %>%
  GGally::ggpairs(aes(alpha = 0.4))
```


## 4PL models in conditional formulation {.tabset}

### MCMC with hierarchical item priors

```{r}
formula_4pl_cond <- bf(
  response2 ~ sigma * zeta + (1 - sigma) * inv_logit(beta + exp(logalpha) * theta),
  nl = TRUE,
  theta ~ 0 + (1 | person),
  beta ~ 1 + (1 |i| item),
  logalpha ~ 1 + (1 |i| item),
  logitsigma ~ 1 + (1 |i| item),
  nlf(sigma ~ inv_logit(logitsigma)),
  logitzeta ~ 1 + (1 |i| item),
  nlf(zeta ~ inv_logit(logitzeta)),
  family = brmsfamily("bernoulli", link = "identity")
)

prior_4pl_cond <- 
  prior("normal(0, 2)", class = "b", nlpar = "beta") +
  prior("normal(0, 1)", class = "b", nlpar = "logalpha") +
  prior("normal(-2, 0.5)", class = "b", nlpar = "logitsigma") +
  prior("normal(-2, 0.5)", class = "b", nlpar = "logitzeta") +
  prior("normal(0, 1)", class = "sd", group = "person", nlpar = "theta") + 
  prior("normal(0, 3)", class = "sd", group = "item", nlpar = "beta") +
  prior("normal(0, 1)", class = "sd", group = "item", nlpar = "logalpha") +
  prior("normal(0, 1)", class = "sd", group = "item", nlpar = "logitsigma") +
  prior("normal(0, 1)", class = "sd", group = "item", nlpar = "logitzeta")

fit_4pl_cond <- brm(
  formula = formula_4pl_cond,
  data = spm_long,
  prior = prior_4pl_cond,
  chains = cores, iter = iter, warmup = warmup,
  cores = cores, control = control, inits = 0,
  file = "models/fit_4pl_cond"
)

fit_4pl_cond <- add_loo(fit_4pl_cond)
```

```{r}
summary(fit_4pl_cond)
```

```{r}
# extract and scale person parameters
person_pars_4pl_cond <- 
  ranef(fit_4pl_cond, summary = FALSE)$person[, , "theta_Intercept"]

person_sds_4pl_cond <- apply(person_pars_4pl_cond, 1, sd)

person_pars_4pl_cond <- person_pars_4pl_cond %>%
  sweep(1, person_sds_4pl_cond, "/") %>%
  posterior_summary() %>%
	as_tibble() %>%
	rownames_to_column(var = "person") %>%
  mutate(person = as.numeric(person))

# plot person parameters
person_pars_4pl_cond %>%
  arrange(Estimate) %>%
  mutate(id2 = seq_len(n())) %>%
	ggplot(aes(id2, Estimate, ymin = Q2.5, ymax = Q97.5)) +
	geom_pointrange(alpha = 0.7) +
	coord_flip() +
	labs(x = "Person Number (sorted after Estimate)") +
  theme(
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  )
```

```{r, fig.height=8}
# extract item parameters
item_pars_4pl_cond <- coef(fit_4pl_cond, summary = FALSE)$item

# locations
beta <- item_pars_4pl_cond[, , "beta_Intercept"] %>%
  posterior_summary() %>%
	as_tibble() %>%
	rownames_to_column()

# slopes
alpha <- item_pars_4pl_cond[, , "logalpha_Intercept"] %>%
	exp() %>%
  sweep(1, person_sds_4pl_cond, "/") %>%
  posterior_summary() %>%
	as_tibble() %>%
	rownames_to_column()

# false answers probability
sigma <- item_pars_4pl_cond[, , "logitsigma_Intercept"] %>%
	inv_logit_scaled()

# conditional guessing probability
zeta <- item_pars_4pl_cond[, , "logitzeta_Intercept"] %>%
	inv_logit_scaled()

# guessing parameters
gamma <- (sigma * zeta) %>%
  posterior_summary() %>%
	as_tibble() %>%
	rownames_to_column()

# lapse parameters
psi <- (sigma * (1 - zeta)) %>%
  posterior_summary() %>%
	as_tibble() %>%
	rownames_to_column()

# plot locations and slope next to each other
bind_rows(beta, alpha, gamma, psi, .id = "nlpar") %>%
	rename(item = "rowname") %>%
	mutate(item = as.numeric(item)) %>%
	mutate(nlpar = factor(
	  nlpar, labels = c("Location", "Slope", "Guessing", "Lapse")
	)) %>%
	ggplot(aes(item, Estimate, ymin = Q2.5, ymax = Q97.5)) +
  scale_x_continuous(breaks = 1:12) +
	facet_wrap("nlpar", scales = "free_x") +
	geom_pointrange() +
	coord_flip() +
	labs(x = "Item Number")
```