# Correlated dimensions in a multidimensional IRT model

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);
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

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)
``````

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

1 Like

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.

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.