Correlated dimensions in a multidimensional IRT model

Modeling
1

I’m trying to fit an IRT model for an assessment that measures 3 dimensions total, but each item measures only 1 dimension. The mean and variance of each dimension are fixed to 0 and 1 respectively. However, I would like to estimate the correlations between the three dimensions. Here is the current iteration of my Stan code:

data {
  int<lower=0> I;                       // number of items
  int<lower=0> J;                       // number of students
  int<lower=0> D;                       // number of dimensions
  int<lower=0> N;                       // number of observations
  int<lower=0,upper=I> ii;              // item for observation n
  int<lower=0,upper=J> jj;              // student for observation n
  int<lower=0,upper=D> dd;              // dimension for observation n
  int<lower=0,upper=1> y;               // score for observation n
}
transformed data {
  vector[D] theta_mean;
  vector[D] theta_scale;
  
  for (d in 1:D) {
    theta_mean[d] = 0;
    theta_scale[d] = 1;
  }
}
parameters {
  matrix[J, D] theta;                   // student ability
  corr_matrix[D] theta_corr;            // theta correlation matrix
  
  vector<lower=0>[I] alpha;             // discriminations
  vector[I] beta;                       // difficulty
}
model {
  vector[N] eta;
  
  // Priors
  theta_corr ~ lkj_corr(1);
  theta ~ multi_normal(theta_mean, quad_form_diag(theta_corr, theta_scale));
  alpha ~ lognormal(0.5, 1);
  beta ~ normal(0, 10);
  
  // Likelihood
  for (n in 1:N)
    eta[n] = alpha[ii[n]] * (theta[jj[n], dd[n]] - beta[ii[n]]);
  y ~ bernoulli_logit(eta);
}

However, I get this error when I try to run the model:


No matches for: 

  matrix ~ multi_normal(vector, matrix)

Available argument signatures for multi_normal:

  vector ~ multi_normal(vector, matrix)
  vector ~ multi_normal(row vector, matrix)
  row vector ~ multi_normal(vector, matrix)
  row vector ~ multi_normal(row vector, matrix)
  vector ~ multi_normal(vector[], matrix)
  vector ~ multi_normal(row vector[], matrix)
  row vector ~ multi_normal(vector[], matrix)
  row vector ~ multi_normal(row vector[], matrix)
  vector[] ~ multi_normal(vector, matrix)
  vector[] ~ multi_normal(row vector, matrix)
  row vector[] ~ multi_normal(vector, matrix)
  row vector[] ~ multi_normal(row vector, matrix)
  vector[] ~ multi_normal(vector[], matrix)
  vector[] ~ multi_normal(row vector[], matrix)
  row vector[] ~ multi_normal(vector[], matrix)
  row vector[] ~ multi_normal(row vector[], matrix)

require real scalar return type for probability function.
  error in 'modelf1103aa62d17_mirt' at line 32, column 77
  -------------------------------------------------
    30:   // Priors
    31:   theta_corr ~ lkj_corr(1);
    32:   theta ~ multi_normal(theta_mean, quad_form_diag(theta_corr, theta_scale));
                                                                                    ^
    33:   alpha ~ lognormal(0.5, 1);
  -------------------------------------------------

Error in stanc(file = file, model_code = model_code, model_name = model_name,  : 
  failed to parse Stan model 'mirt' due to the above error.

So the issue is that theta can’t be a matrix when defining the multi_normal prior. However, I do need to be able to access specific elements of theta in the likelihood (i.e., theta[jj[n], dd[n]]). I tried following the example model in section 1.13 of the Stan manual, and defining theta as:

parameters {
  vector[D] theta[J];
  ...
}

This resolves the problem with the multi_normal prior, but I am then unable to access specific elements of theta (or at least I haven’t been able to figure it out).

Is there a way to estimate these correlations between dimensions, but also access the elements individually?

1 Like
2

Here is some R code to simulate example data and estimate the model:

library(tidyverse)
library(mvtnorm)
library(glue)
library(rstan)


### Simulation parameters ------------------------------------------------------
nstu <- 1000
ndim <- 3
nitm <- ndim * 10


### Functions ------------------------------------------------------------------
logit <- function(x) {
  log(x / (1 - x))
}
inv_logit <- function(x) {
  1 / (1 + exp(-x))
}


### Make Q-matrix --------------------------------------------------------------
qmatrix <- map_dfc(seq_len(ndim), function(x, nitm, ndim) {
  att_name <- glue("dim_{x}")
  totalone <- nitm / ndim
  
  tibble(item_id = seq_len(nitm),
         section = rep(seq_len(ndim), each = totalone)) %>%
    mutate(!!att_name := case_when(section == x ~ 1L, TRUE ~ 0L)) %>%
    select(!!att_name)
},
                   nitm = nitm, ndim = ndim) %>%
  rowid_to_column(var = "item_id")


### Simulate MIRT data ---------------------------------------------------------
# People
theta_cor <- matrix(data = c(1.0, 0.5, 0.2,
                             0.5, 1.0, 0.6,
                             0.2, 0.6, 1.0),
                    nrow = ndim, byrow = TRUE)

theta <- rmvnorm(n = nstu, mean = rep(0, ndim), sigma = theta_cor) %>%
  as_tibble(.name_repair = ~ glue("dim_{seq_len(ndim)}")) %>%
  rownames_to_column(var = "stu_id")

# Items
a_mean <- 1.25
a_sd <- 0.5
rl_loc <- log(a_mean^2 / sqrt(a_sd^2 + a_mean^2))
rl_shp <- sqrt(log(1 + (a_sd^2 / a_mean^2)))

item_params <- tibble(item_id = seq_len(nitm),
                      a = rlnorm(nitm, meanlog = rl_loc, sdlog = rl_shp),
                      b = rnorm(nitm, mean = 0, sd = 0.7))

# Response data
mirt_response <- theta %>%
  nest(-stu_id) %>%
  mutate(resp_data = map(data, function(x, item_params, qmatrix) {
    qmatrix %>%
      gather(key = "dim", value = "value", -item_id) %>%
      filter(value == 1L) %>%
      left_join(gather(x, key = "dim", value = "theta"), by = "dim") %>%
      left_join(item_params, by = "item_id") %>%
      mutate(prob = inv_logit(a * (theta - b)),
             rand = runif(n = nrow(.), min = 0, max = 1),
             score = case_when(rand <= prob ~ 1L, TRUE ~ 0L)) %>%
      select(item_id, score)
  },
                         item_params = item_params, qmatrix = qmatrix)) %>%
  select(-data) %>%
  unnest() %>%
  left_join(
    qmatrix %>%
      gather(key = "dim", value = "value", -item_id) %>%
      filter(value == 1L) %>%
      separate(dim, into = c("junk", "dim"), convert = TRUE) %>%
      select(item_id, dim),
    by = "item_id"
  )

# Stan data
stan_mirt <- list(
  I = nitm,
  J = nstu,
  D = ndim,
  N = nrow(mirt_response),
  ii = mirt_response$item_id,
  jj = mirt_response$stu_id,
  dd = mirt_response$dim,
  y = mirt_response$score
)

# Stan model
mirt_model <- stan(file = "stan-models/mirt.stan", data = stan_mirt,
                   chains = 4, iter = 2000, warmup = 1000, cores = 4)
3

You would have to do a nested loop inside the loop that creates the elements of eta.

1 Like
4

Hi wjakethompson
I am a student who is just learning stan for item response model, and I want to learn how to using stan to fit a multidimensional IRT model. So I find your post here. I learned a lot from your posting code.
Could you please tell me the right code and I am full of gratitude to you for helping me.

5

Hey @Lucius,

The code I ended up using is here. There is also some code that I used to simulate example data and estimate the models in the same repo. But note that this code is only set up to estimate a MIRT model where each item measures only one dimension. If you have items that measure multiple dimensions, this will need some modifications.

Hope this is helpful!

1 Like
6

Hey @wjakethompson
These codes are very beautiful, and these are really helpful!. Thank you very much for your help.