Poisson adaptation of IRT model: how to interpret parameters?

Question

For my work I made an assessment in which users have to stack blocks in a certain configuration in as few steps as possible. If the user does this in the minimally required number of steps, they get the perfect score of 0. With every extra step they need, their score increases by 1.

To measure the skill of the participants, we use an IRT (item response theory) model. Since the test isn't dichotomous, scores are discrete, scores on every item can range from 0 to $\infty$, and scores seem roughly normally distributed, we adapted the IRT model to use a poisson distribution.

That is, we use the normal MaxLogLikelihood estimation, but instead of using the 2-parameter IRT fomula(1) $$ P(X|a,b,\theta) = {e^{a(\theta - b)} \over 1+ e^{a(\theta - b)}} $$

we calculate P according to the regular poisson formula (2) $$ P(X=n) = {\lambda^n e^{-\lambda} \over n!} $$

where $\lambda$ is ... well, that's the question. Based on a textbook I have been using (3) $$ \lambda = e^{(b - a\theta)} $$

My questions are about the relationship between $b$, $\lambda$ and $\theta$:
- does the formula at (3) make sense, or should it be something else?
- in the 2plm IRT model (1) the roles of the parameters are very clear to me: $\theta$ is assumed/defined to be N(0,1) so ranges roughly from -3 to 3. $b$ is defined such that if $b = \theta$, the chance of getting a correct answer is exactly 50%. And $a$ defines the slope. Makes perfect sense. In the poisson model however, I lack these definitions. If $b$ still represents the difficulty of an item, and suppose $b = \theta$, what $\lambda$ should I expect? What does $a$ represent? I'm finding it hard to wrap my head around the conceptual meaning.

In short: keeping $\theta$ defined as N(0,1), what would be the best/most sensible way to define the other parameters of a poisson-based IRT model?

p.s. first post here, feel free to give other feedback.

Answer 1

I recently came across a paper by Daniel Bauer and Andrea Hussong (2010, available for free on PubMed Central) that appears to tackle this topic. As mentioned in comments, I think that what you want to do can be described as a non-linear confirmatory factor analysis. The paper treated 5 variables as indicators of alcohol involvement, and the variables are all of different types (one is binary, one is Poisson, two are unordered categorical).

Their equation 3 shows how you'd model item responses in general terms (altered notation to be consistent with original question):

$g_i(\lambda_{ij}) = v_i + d_i\theta_j$

Where i indexes items/questions and j indexes persons

$\lambda{_ij}$ is the mean of item i for the j-th person

$g_i$ is a link function like you might find in a GLM (e.g. log link for Poisson, logit for binary items)

$v_i$ is the intercept for the i-th item, which I believe is a parallel concept to difficulty or severity (note @Heteroskedastic_Jim's comment that in the OP's context, because higher count = lower ability, this should probably be interpreted as easiness rather than difficulty)

$d_i$ is the factor loading for the i-th item, which I believe is parallel to discrimination (note: the paper used $\lambda_i$ for this term, but I'm reserving this for the Poisson parameter; also note Jim's comment that this would be on the log scale given the context)

And $\theta_j$ is the j-th person's trait level, whatever trait it is you're interested in.

---

For Poisson items, I would suggest the equation doesn't really need to be rearranged other than specifying the log link. In regular Poisson regression, our model is:

$ln(E(Y|X)) = ln(\lambda|X) = \alpha + X\beta$

In the original poster's context, our model is now:

$ln(\lambda_{ij}) = v_i + d_i\theta_j$

So, if I am correct, the original question gave the wrong formula for $\lambda_i$, but the correct probability mass function.

---

For Bernoulli items, like you'd find in a basic 1- or 2-parameter logistic (1- and 2-PL) IRT model, the link is the logit link, i.e.

$ln(\frac{\lambda_{ij}}{1 - \lambda_{ij}}) = v_i + d_i\theta_j$

Rearrange terms and you should get the IRT parameterization that many people are familiar with. For Gaussian items, then $g_i$ is most likely the identity link, and the equation becomes something you'd see in linear SEM. In principle, any family and link function supported by your software could be used to model the items, e.g. you could model count items as negative binomial, or you could model Gaussian item + log link.

I don't know any specific R packages, but they have to exist. The authors used SAS PROC NLMIXED to fit their model. Stata's gsem command will also fit this model, as demonstrated at the link.