4
$\begingroup$

I'm trying to better understand Item response Theory (IRT) from a Bayesian perspective. Hypothetically, suppose I want to use a 1PL model and my data is a binary matrix

data = np.array([[1,1,1,1],
                 [1,1,1,0],
                 [1,1,0,0],
                 [1,0,0,0],
                 [0,0,0,0]])

There are five children and four questions. Children are depicted by rows, whereas questions are depicted by columns; if a 1 is present, the child correctly answered the question (else 0.) The first row corresponds to the smartest child and the last column corresponds to the hardest question.

My understanding of IRT is that we are simultaneously evaluating latent features, namely, namely child ability and question difficulty. The response variable, y, is {0,1} correct or incorrect and as such we use the 1PL model where the characteristic curve is described by $$ p(\theta|b_i) = \frac{\exp(\theta-b_i)}{1+\exp(\theta-b_i)} $$

For priors, I have arbitrarily selected normal distribution(s) with mean 1 and sigma 1, encouraging my latent feature variables to take on positive values. For a likelihood function, much like logistic regression, I've selected Bernoulli.

To explore the concepts, I've written a Metropolis sampler:

import numpy as np
import random

def PL1(ability, difficulty):
  return ability - difficulty

def sigmoid(z):
  return 1/(1 + np.exp(-z))

def normal(x,mu,sigma):
  num = np.exp(-1/2*((x-mu)/sigma)**2)
  den = np.sqrt(2*np.pi)*sigma
  return num/den

def bernoulli(y,p):
  return p**y*(1-p)**(1-y)

def cum_log_lik(A,D,Y):
  log_lik = 0
  for idx_a in range(len(A)):
    for idx_d in range(len(D)):
      z = sigmoid(PL1(A[idx_a],D[idx_d]))
      log_lik += np.log(bernoulli(y=Y[idx_a,idx_d],p=z))
  return log_lik

def cum_log_prior(A,D):
  log_prior = 0
  for a in A:
    log_prior += np.log(normal(x=a,mu=1,sigma=1))
  for d in D:
    log_prior += np.log(normal(x=d,mu=1,sigma=1))

  return log_prior

def MCMC(data,hops=10_000):
  u_dim = data.shape[0]
  v_dim = data.shape[1] 
  U = np.random.uniform(low=0,high=1,size=u_dim)
  V = np.random.uniform(low=0,high=1,size=v_dim)
  Y = data
  curr_log_lik = cum_log_lik(U,V,Y)
  curr_log_prior = cum_log_prior(U,V)
  current = curr_log_lik + curr_log_prior
  U_arr = []
  V_arr = []

  for epoch in range(hops):
    U_arr.append(U)
    V_arr.append(V)

    if epoch%2==0: #update U
      mov_U = U + np.random.uniform(low=-0.25,high=0.25,size=u_dim)
      mov_V = V
    
    else: #update V
      mov_U = U
      mov_V = V + np.random.uniform(low=-0.25,high=0.25,size=v_dim)
    
    mov_log_lik = cum_log_lik(mov_U,mov_V,Y)
    mov_log_prior = cum_log_prior(mov_U,mov_V)
    movement = mov_log_lik + mov_log_prior
    ratio = np.exp(movement - current)

    event = random.uniform(0,1)
    if event <= ratio:
      U = mov_U
      V = mov_V
      current = movement

  return np.array(U_arr), np.array(V_arr)
    
A, D = MCMC(data,hops = 50_000)

Now, to evaluate my sampler's performance:

def get_estimate(arr,idx):
  vec = [arr[i][idx] for i in range(len(arr))]
  return sum(vec)/len(vec)

for a in range(5):
  print(get_estimate(A,a))

>>>
2.356836411120115
1.4854360638445205
0.8823022398184828
0.40257074505614127
-0.14228691392908904

for a in range(4):
  print(get_estimate(D,a))

>>>
0.28806026673506735
0.7268234141444485
1.215012903954542
1.8960656959448172

My code does work. It accurately evaluates child ability and question difficulty. The problem I'm running into is that, I've been told by multiple sources that 1PL only evaluates one parameter, namely, question difficulty.

How can this be? Do we treat child ability as a constant, model it as a variable but not include it in our results, something else...?

Bonus points if you can update the MH sampler above to reflect the correct design.

$\endgroup$
  • 2
    $\begingroup$ Please post formulas as text, not as a photograph or screenshot (see here). Also, I switched the [mcmc] tag for [psychometrics], as I think that will help bring your Q to the attention of people who can best answer it. If you disapprove, switch it back w/ my apologies. $\endgroup$ – gung - Reinstate Monica Sep 28 at 13:02
  • $\begingroup$ Of possible interest: py-irt. $\endgroup$ – chl 2 days ago
1
$\begingroup$

In his original formulation, Rasch treated ability as fixed --- and such model are nowadays fitted using a conditional approach, which is what pure Rasch modelers prefer, for theoretical reason --- hence the name 1-PL (for item difficulty), but there are other approaches like the joint maximum likelihood technique (poorly recommended) or mixed-effect models, which are more flexible and allow to fit a wider range of models. In the latter case, individuals are treated as random effects and ability estimates are drawn from the posterior distribution (expected a posteriori, EAP). Bayesian models were developed in the 80s but newer models became available after 2000, especially for multidimensional and multilevel IRT models. Generally, they also rely on EAP estimation using MCMC. See Fox, J.-P., Bayesian item response modeling: theory and applications, Springer, (2010), for a review.

More generally the k in k-PL refers to the number of item characteristics that are estimated from the sample: 1-PL for difficulty alone (intercept of the item characteristic curve), 2-PL for difficulty and discrimination (intercept and slope), and 3-PL for the 2-PL with a guessing parameter (intercept, slope and lower asymptote). There's even a 1-PL model with fixed but pre-defined discrimination parameter (usually it is fixed to 1 by most software) --- I forgot the name but it was developed by the Cito team in the Netherlands.

| cite | improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.