I’m trying to generalize a standard two-parameter logistic IRT model (with one-dimensional latent ability) to a multi-dimensional version. However, I’ve been stuck for days now with the slicing and shaping of the involved arrays, and would greatly appreciate some help.
Here’s the problem. The classic 2PL is
where D=-1.7, \alpha_j and \delta_j are the discrimination and difficulty parameters, respectively, for item j = 1, 2, \ldots, n_{items}, while \theta_{i} is the latent ability for examinee i = 1, 2, \ldots, n_{persons}. The PyMC code is to implement this is:
import numpy as np
import pymc as pm
n_persons = 100 # Number of persons
n_items = 25 # Number of items
# Random data for testing
responses = np.random.randint(0, 2, size=(n_persons, n_items))
with pm.Model() as IRT_2PL:
# Parameters priors
discrim = pm.LogNormal('discrim', mu=0, sigma=1, shape=n_items)
diff = pm.Normal('diff', mu=0, sigma=1, shape=n_items)
# Abilities prior
theta = pm.Normal('theta', mu=0, sigma=1, shape=n_persons)
# 2PL Model
p = pm.math.invlogit(1.7 * discrim * (theta[:, None] - diff))
# Likelihood
obs = pm.Bernoulli('obs', p=p, observed=responses)
trace = pm.sample(2000)
The generalization to a multidimensional item response (MIRT) model is straightforward in linear algebra
where \vec{\alpha}_{j} = \begin{bmatrix} \alpha_{j1} & \alpha_{j2} \end{bmatrix}^{\mathsf{T}}, \vec{\theta}_{i} = \begin{bmatrix} \theta_{i1} & \theta_{i2} \end{bmatrix}^{\mathsf{T}}, and \vec{\gamma}_{j} = - \left( \alpha_{j1} \delta_{j1} + \alpha_{j2} \delta_{j2} \right) (more details in de Ayala, 2022, chapter 10).
I’m struggling with the implementation of this in an elegant vectorized way. Because a crude version would obviously be:
with pm.Model() as crude_MIRT_2PL:
# Parameters priors
discrim = pm.LogNormal('discrim', mu=0, sigma=1, shape=n_items)
discrim2 = pm.LogNormal('discrim2', mu=0, sigma=1, shape=n_items)
diff = pm.Normal('diff', mu=0, sigma=1, shape=n_items)
diff2 = pm.Normal('diff2', mu=0, sigma=1, shape=n_items)
# Abilities prior
theta = pm.Normal('theta', mu=0, sigma=1, shape=n_persons)
theta2 = pm.Normal('theta2', mu=0, sigma=1, shape=n_persons)
# 2PL Model
p = pm.math.invlogit(1.7 * discrim * (theta[:, None] - diff) + 1.7 * discrim2 * (theta2[:, None] - diff2))
# Likelihood
obs = pm.Bernoulli('obs', p=p, observed=responses)
But this non-vectorized version just makes it so much harder to later implement structural zeros among the discrimination parameters to exclude some of the thetas from certain items. So ideally one would like to work with 2D-arrays (matrices) for the parameters and abilities, i.e.
# ...
# Priors for item parameters
discrim = pm.Normal('discrim', mu=0, sigma=1, shape=(n_items, 2))
diff = pm.Normal('diff', mu=0, sigma=1, shape=(n_items, 2))
# Priors for person abilities
theta = pm.Normal('theta', mu=0, sigma=1, shape=(n_persons, 2))
# ...
But for the life of me I can’t figure out how to slice and dice these arrays in the right shape so that the p=... equation works. Any help would be appreciated.
