Modeling performance on a multiple choice exam

This notebook includes the code for the Bayesian Workflow book Chapter 4 Introduction to workflow: Modeling performance on a multiple choice exam.

1 Introduction

We demonstrate with an example from Section 4.15 of Gelman and Vehtari (2024), assessing the grading of a multiple-choice test. We analyze data from a 24-question final exam from a class of 32 students, where our applied goal is to check that the individual test questions are doing a good job at discriminating between poorly- and well-performing students. No external data are available on the students, so we assess their abilities using their total score on the exam.

Each item is a multiple choice question with 4 possible answers and is scored either 1 (correct) or 0 (incorrect). The total scores across the 24 question exam range from 12 to 21, with an average of 16. Across the questions the hardest question has 4 of the 32 students answering it correctly, while all students manage to answer the easiest question correctly.

Load packages

library(rprojroot)
root <- has_file(".Bayesian-Workflow-root")$make_fix_file()
library(cmdstanr)
options(mc.cores = 4)
library(posterior)
library(arm)
set.seed(123)
# Bring in plotting functions from a separate file
source(root("multiple_choice", "plot_functions.R"))

print_stan_code <- function(code) {
  if (isTRUE(getOption("knitr.in.progress")) &
        identical(knitr::opts_current$get("results"), "asis")) {
    # In render: emit as-is so Pandoc/Quarto does syntax highlighting
    block <- paste0("```stan", "\n", paste(code, collapse = "\n"), "\n", "```")
    knitr::asis_output(block)
  } else {
    writeLines(code)
  }
}

2 Stan models

Compile all Stan programs to use throughout the file

logit_0 <- cmdstan_model(root("multiple_choice", "logit_0.stan"))
logit_prior <- cmdstan_model(root("multiple_choice", "logit_prior.stan"))
logit_guessing <- cmdstan_model(root("multiple_choice", "logit_guessing.stan"))
logit_guessing_uncentered <- cmdstan_model(root("multiple_choice", "logit_guessing_uncentered.stan"))
logit_guessing_multilevel <- cmdstan_model(root("multiple_choice", "logit_guessing_multilevel.stan"))
logit_guessing_uncentered_multilevel <- cmdstan_model(root("multiple_choice", "logit_guessing_uncentered_multilevel.stan"))
logit_guessing_multilevel_bivariate <- cmdstan_model(root("multiple_choice", "logit_guessing_multilevel_bivariate.stan"))
logit_guessing_multilevel_bivariate_cholesky <- cmdstan_model(root("multiple_choice", "logit_guessing_multilevel_bivariate_cholesky.stan"))
irt_guessing <- cmdstan_model(root("multiple_choice", "irt_guessing.stan"))
irt_guessing_discrimination <- cmdstan_model(root("multiple_choice", "irt_guessing_discrimination.stan"))

3 Data

Read in data and construct score for each student

responses <- read.csv(root("multiple_choice", "data", "final_exam_responses.csv"))
answers <- read.csv(root("multiple_choice", "data", "final_exam_answers.csv"))
J <- nrow(responses)  # number of students
K <- ncol(responses)  # number of items
correct <- array(NA, c(J,K))
for (k in 1:K){
  correct[,k] <- ifelse(responses[,k] == as.character(answers[k]), 1, 0)
}
score <- rowSums(correct)
item <- colSums(correct)
summary(score)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  12.00   14.75   16.00   16.09   17.00   21.00 

summary(item)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   4.00   15.25   24.00   21.46   27.25   32.00 

score_jitt <- score + jitter(rep(0, J), amount = 0.3)
score_adj <- (score - mean(score)) / sd(score)
score_adj_jitt <- (score_jitt - mean(score)) / sd(score)
data <- list(J = J, x = score, y = correct)

item_id_0 <- LETTERS[1:J]  # Only works here because J is no more than 26!

4 Simple models

4.1 Base model logit_0

print_stan_code(logit_0$code())

data {
  int J;
  array[J] int<lower=0, upper=1> y;
  vector[J] x;
}
parameters {
  real a, b;
}
model {
  y ~ bernoulli_logit(a + b*x);
}

plot_logit(
  "Fit to item A on exam",
  "Score on exam",
  logit_0,
  list(J = data$J, x = data$x, y = data$y[, 1]),
  score_jitt,
  guessprob = 0
)

 variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
     lp__ -14.30 -13.92 1.20 0.78 -16.55 -13.23 1.01      681      646
     a    -11.85 -11.26 5.19 4.78 -21.06  -4.51 1.01      344      370
     b      0.85   0.81 0.35 0.32   0.36   1.47 1.01      341      367

Figure 1

4.2 Add priors

print_stan_code(logit_prior$code())

data {
  int J;
  array[J] int<lower=0, upper=1> y;
  vector[J] x;
  real mu_a, mu_b;
  real<lower=0> sigma_a, sigma_b;
}
transformed data {
  vector[J] x_adj = (x - mean(x))/sd(x);
}
parameters {
  real a, b;
}
model {
  a ~ normal(mu_a, sigma_a);
  b ~ normal(mu_b, sigma_b);
  y ~ bernoulli_logit(a + b*x_adj);
}

plot_logit(
  "Fit to item A:  rescaled predictor and weakly informative prior",
  "Standardized exam score",
  logit_prior,
  list(J = data$J, x = data$x, y = data$y[, 1], 
       mu_a = 0, sigma_a = 5, mu_b = 0, sigma_b = 5),
  score_adj_jitt,
  guessprob = 0
)

 variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
     lp__ -14.35 -14.03 1.05 0.76 -16.39 -13.34 1.00     1449     1767
     a      1.76   1.72 0.62 0.62   0.83   2.88 1.00     1635     1542
     b      1.87   1.82 0.71 0.70   0.79   3.09 1.00     1587     1697

Figure 2

plot_logit_grid(
  "Rescaled predictor and weakly informative prior",
  "Standardized exam score",
  logit_prior,
  c(data, list(mu_a = 0, sigma_a = 5, mu_b = 0, sigma_b = 5)),
  score_adj_jitt,
  item_id_0,
  guessprob = 0
)

Figure 3

plot_logit(
  "Fit to item G: rescaled predictor and weakly informative prior",
  "Standardized exam score",
  logit_prior,
  list(J = data$J, x = data$x, y = data$y[, 7], 
       mu_a = 0, sigma_a = 5, mu_b = 0, sigma_b = 5),
  score_adj_jitt,
  guessprob = 0
)

 variable  mean median   sd  mad    q5   q95 rhat ess_bulk ess_tail
     lp__ -1.74  -1.42 1.01 0.73 -3.85 -0.77 1.00     1010     1691
     a     7.68   7.27 2.72 2.68  3.96 12.73 1.00     1152     1717
     b     0.08   0.02 2.09 1.96 -3.31  3.64 1.00      824     1142

Figure 4

4.3 Fix the data problems

answers[c(5, 14, 17)] <- c("d", "d", "c")
correct <- array(NA, c(J,K))
for (k in 1:K){
  correct[,k] <- ifelse(responses[,k] == as.character(answers[k]), 1, 0)
}
score <- rowSums(correct)
score_jitt <- score + jitter(rep(0, J), amount = 0.3)
score_adj <- (score - mean(score)) / sd(score)
score_adj_jitt <- (score_jitt - mean(score)) / sd(score)
data <- list(J = J, x = score, y = correct)
item_id <- rank(colSums(correct), ties = "first")

plot_logit_grid(
  "After fixing the data problem",
  "Standardized exam score",
  logit_prior,
  c(data, list(mu_a = 0, sigma_a = 5, mu_b = 0, sigma_b = 5)),
  score_adj_jitt,
  item_id,
  guessprob = 0
)

Figure 5

4.4 Allow for guessing

print_stan_code(logit_guessing$code())

data {
  int J;
  array[J] int<lower=0, upper=1> y;
  vector[J] x;
  real mu_a, mu_b;
  real<lower=0> sigma_a, sigma_b;
}
transformed data {
  vector[J] x_adj = (x - mean(x))/sd(x);
}
parameters {
  real a, b;
}
model {
  a ~ normal(mu_a, sigma_a);
  b ~ normal(mu_b, sigma_b);
  y ~ bernoulli(0.25 + 0.75*inv_logit(a + b*x_adj));
}

plot_logit_grid(
  "Probabilities constrained to range from 0.25 to 1",
  "Standardized exam score",
  logit_guessing,
  c(data, list(mu_a = 0, sigma_a = 5, mu_b = 0, sigma_b = 5)),
  score_adj_jitt,
  item_id,
  guessprob = 0.25
)