[Basic of IRT using R] 5. Estimating an Examinee's Ability

5.1 Introduction

문항반응이론에서 시험 응시자를 척도 위에 올려놓는 것이 가장 중요한 목적이다. 

응시자의 능력을 추정할 수 있어야 선별이나 순위 등을 매길 수 있다.

 

5.2 Ability Estimation Procedures

MLE를 통해서 응시자의 능력을 추정해보자.

 

$$ \hat{\theta}_{s+1} = \hat{\theta_s} - \displaystyle{\frac{\sum_{j=1}^J a_j[u_j - P_j(\hat{\theta_s})]}{-\sum_{j=1}^J a_j^2 P_j(\hat{\theta_s})Q_j(\hat{\theta_s})}}$$

 

\(\hat{\theta}_s\): 추정된 능력, s번 반복

\(a_j\): 변별도 모수

\(u_j\): 응시자가 문항 j에 한 응답

\(u_j=1 \) 정답일 경우

\(u_j=0\) 오답일 경우

\(P_j(\hat{\theta}_s)\): 문항j의 정답률

\(Q_j(\hat{\theta_s}) = 1- P_j(\hat{\theta}_s)\) 

 

초기값은 임의의 숫자로 한다(ex=1).

 

$$\hat{\theta}_{s+1} = \hat{\theta}_s - \Delta \hat{\theta}_s $$

 

위의 식을 반복하면 차이가 줄어든다. 많이 줄어들 때 까지 반복을 돌린다고 한다.

 

세개의 문항을 만들어서 테스트를 해보자.

\(b_1 = -1.0, b_2 =0.0, b_3 = 1.0\)

\(a_1 = 1.0, a_2 = 1.2, a_3 =0.8\)

 

응시자의 반응은

\(u_1 =1, u_2=0, u_3=1\)

\(\hat{\theta}_s = 1.0\)으로 설정한 후에 이를 반복문을 돌려보자.

\(P_j(\hat{\theta}_s) = \frac{1}{1+e^{-a_j(\hat{\theta}_s - b_j)}}\)

 

이를 표로 만들어보자.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

b=[-1.0,0.0,1.0]
a=[1.0,1.2,0.8]
u=[1,0,1]
theta=[1.0,1.0,1.0]
p=[]
ja=[]
mo=[]
q=[]

for i in range(len(theta)):
  P=1/(1+np.exp(-a[i]*(theta[i]-b[i])))
  bunja=a[i]*(u[i]-P)
  Q=1-P
  bunmo=(a[i]**2) * P * Q
  p.append(P)
  q.append(Q)
  ja.append(bunja)
  mo.append(bunmo)

pd.DataFrame({'u':u,'p':p,'q':q,'a[u-p]':ja,'a^2Pq':mo})

 

u	p	q	a[u-p]	a^2Pq012
1	0.880797	0.119203	0.119203	0.104994
0	0.768525	0.231475	-0.922230	0.256168
1	0.500000	0.500000	0.400000	0.160000

 

 

한번 더 루프를 돌려보자.

theta = theta - (sum(ja)/-sum(mo))
theta

p=[]
ja=[]
mo=[]
q=[]
for i in range(len(b)):
  P=1/(1+np.exp(-a[i]*(theta-b[i])))
  bunja=a[i]*(u[i]-P)
  Q=1-P
  bunmo=(a[i]**2) * P * Q
  p.append(P)
  q.append(Q)
  ja.append(bunja)
  mo.append(bunmo)
pd.DataFrame({'u':u,'p':p,'q':q,'a[u-p]':ja,'a^2Pq':mo})

u	p	q	a[u-p]	a^2Pq
0	1	0.773236	0.226764	0.226764	0.175342
1	0	0.567587	0.432413	-0.681104	0.353422
2	1	0.350086	0.649914	0.519931	0.145617

 

 

이제 표준오차를 구해보자.

$$ SE(\hat{\theta} = \frac{1}{\displaystyle{\sqrt{\displaystyle{\sum_{j=1}^J} a_j^2 P_j(\hat{\theta})Q_j(\hat{\theta})}}}$$

 

표준오차는 추정된 값이 좋은지 안좋은지를 알아본다. 

이 값이 작아야 좋은 것이다.

위의 식을 계산해보면 

\(SE=\hat{\theta} = 0.32) = 1/\sqrt{0.6616}=1.22944624\)

표준오차가 굉장히 높다. 문항이 3개뿐이라 정확한 추측이 어렵다.

 

 

5.3 Item Invariance of an Examinee's Ability Estimate

문항반응이론은 문항 불변성을 갖고 있다.

고전이론에서는 쉬운 시험에서 높은 점수를 받고 어려운 시험에서 낮은 점수를 받다. 그러면 이것을 가지고 능력을 알아내기는 어렵다. 그러나 문항반응이론은 난이도 등의 모수 등을 활용해 개인을 척도 위에 위치시킬 수 있다.

 

 

5.4 Computer Session

 

b=[-1.0,0.0,1.0]
a=[1.0,1.2,0.8]
u=[1,0,1]
th=1.0
J=len(b)
S=10
ccrit=0.001
for s in range(S):
  sumnum=0.0
  sumdem=0.0
  for j in range(J):
    phat = 1/(1+np.exp(-a[j]*(th-b[j])))
    sumnum=sumnum+a[i]*(u[j]-phat)
    sumdem=sumdem - a[j]**2 * phat * (1.0-phat)
    delta = sumnum/sumdem
    th=th-delta
    print(f'th={th}, \n')
    if(abs(delta)<ccrit or s==S):
      se=1/np.sqrt(-sumdem)
      print(f'se={se}\n')
      break

th=1.9082682265892903, 
th=-0.8934974040486872, 
th=-0.8161104872416066, 
th=0.6495116945159266, 
th=0.3179490908014487, 
th=0.7725343041409092, 
th=1.7084557919515735, 
th=-0.48937598398457105, 
th=-0.4346371467216754, 
th=0.8198860155394085, 
th=0.2534086621522046, 
th=0.5907215191067869, 
th=1.553744341968402, 
th=-0.25591890477208823, 
th=-0.19049205994885055, 
th=0.9655695535905658, 
th=0.20240844704925476, 
th=0.46216043238276716, 
th=1.4475485008219575, 
th=-0.12893936597503597, 
th=-0.04434567520915192, 
th=1.0633026726195913, 
th=0.1598179153255388, 
th=0.3742005935685893, 
th=1.376634032139736, 
th=-0.05722855597655174, 
th=0.04475892120173748, 
th=1.1261802132337766, 
th=0.12719793595871487, 
th=0.3148971563439983, 

 

 

함수화

def ability(mdl,u,b,a=np.repeat(1,J),c=np.repeat(0,J)):
  J=len(b)
  x=sum(u)
  if(x==0):
    th=-np.log(2*J)
  if(x==J):
    th=np.log(2*J)
  if(x==0 or x==J):
    sumdem =0.0
    for j in range(J):
      pstar= 1/(1+np.exp(-a[j]*(th-b[j])))
      phat = c[j] + (1.0 - c[j])*pstar
      sumdem = sumdem - a[j]**2 * phat * (1.0 -phat) * (pstar/phat)**2
    se=1/np.sqrt(-sumdem)
  if(x!=0 and x != J):
    th=np.log(x/(J-x))
    S=10
    ccrit=0.001
    for s in range(S):
      sumnum=0.0
      sumdem=0.0
      for j in range(J):
        pstar= 1/(1+np.exp(-a[j]*(th-b[j])))
        phat = c[j] + (1.0 - c[j])*pstar
        sumnum=sumnum+a[j]*(u[j]-phat)*(pstar/phat)
        sumdem=sumdem - a[j]**2 * phat * (1.0 -phat) * (pstar/phat)**2

      delta=sumnum/sumdem
      th=th-delta
      if(abs(delta)<ccrit or s==S):
        se=1/np.sqrt(-sumdem)
        break
    print(f'th={th} \n')
    print(f'th={th} \n')
    thse=[th,se]
    return thse

u=[1,0,1]
b=[-1.0,0.0,1.0]
a=[1.0,1.2,0.8]
ability(2,u,b,a,)

th=0.32484627572812413 

th=0.32484627572812413 

[0.32484627572812413, 1.2295613112833954]

mdl=2
theta=0.5
b=[-0.5,-0.25,0.0,0.25,0.5]
a=[1.0,1.5,0.7,0.6,1.8]
J=len(b)
if(mdl==1 or mdl==2):
  c=np.repeat(0,J)
if(mdl==1):
  a=np.repeat(1,J)
sudemet=0.0
for j in range(J):
  Pstar=1/(1+np.exp(-a[j]*(theta-b[j])))
  p= c[j]+ (1-c[j])*Pstar
  sudemet=sudemet - a[j]**2 * P * (1.0 - P)* (Pstar/P)**2
set1 = 1/np.sqrt(-sudemet)
R=10
thr=np.repeat(0,R)
ser=np.repeat(0,R)
for r in range(R):
  u=np.repeat(0,J)
  for j in range(J):
    P = c[j]+ (1-c[j])/(1+np.exp(-a[j]*(theta-b[j])))
    u[j] =np.random.binomial(1,P)
  thse =ability(mdl,u,b,a,c)
  thr[r]=thse[0]
  ser[r]=thse[1]

theta, set1, thr, np.mean(thr), np.std(thr),ser

(0.5,
 0.5872898811308709,
 array([1, 0, 0, 0, 0, 0, 0, 0, 0, 0]),
 0.1,
 0.30000000000000004,
 array([1, 0, 0, 0, 0, 0, 0, 0, 0, 0]))