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.
[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$