[Basic of IRT using R] 3. Estimating Item Parameters

3.1 Introduction

여태까지는 모수(난이도, 변별도, 추측도)를 안다는 가정하에 그래프를 그렸지만 실제로 모수는 추정되어야한다.

이 모수추정방법을 배워보자.

 

3.2 Maximum Likelihood Estimation of Item Parameters

N명의 사람이 J개의 문항에 응답한다고 해보자

이 때 사람을 능력에 따라 G 그룹으로 나눌 수 있다.

각 그룹은 같은 능력을 갖고 있고, \(G = \{g_1, g_2, \dots, G\}\)

능력 \(\theta_g\)을 가진 사람은 정답률이 \(p(\theta_g) = \frac{r_g}{f_g}\)이다.

\(r_g\) 올바르게 답한 횟수이며, \(f_g\)는 총 문제수이다.

 

처음으로 해야할 일은 문항 특성 곡선을 찾는 것이다. 

MLE를 통해 이를 할 수 있는데 아쉽게 이 책에서는 MLE라고 이름을 써놓고 그 방법에 대해서는 말하지 않고 있다.

반대로 RMSE를 통한 방법으로 찾고 있는데 이 방법을 한 번 봐보자.

 

 

우리의 목적은 위처럼 자료에 대해 가장 잘 들어맞는 문항특성 곡선을 찾는 것이다.

이를 RMSE를 통해 해보자.

 

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

theta=np.arange(-3,3.1875,0.1875)
f=[21 for _ in range(len(theta))]
wb=-0.39
wa=1.27

가중치를 변별도는 1.27, 난이도는 -0.39로 하고 그에 맞는 데이터를 만들 것이다.

 

for g in range(1,len(theta)):
  p=1/(1+np.exp(-wa*(theta-wb)))
  r= np.random.binomial(f, p)

array([ 1, 1, 0, 1, 3, 3, 4, 5, 3, 7, 3, 7, 8, 9, 12, 12, 11, 13, 16, 12, 16, 20, 17, 18, 21, 18, 20, 21, 21, 20, 21, 20, 21])

 

np.random.binomial은 이항 분포로부터 표본을 생성하는 NumPy 함수입니다. 이항 분포는 이진 결과를 갖는 실험에서 성공의 횟수를 나타내는 확률 분포입니다. 이 함수는 주어진 확률(p)로 성공하는 베르누이 시행을 n번 반복하여 성공의 횟수를 반환합니다. [챗지피티]

 

한 그룹에는 21명이 있고, 그룹은 33개가 있다.

그리고 각 사람이 정답을 맞춘 비율을 확률로 만든다.

P=r/f
P

array([0.04761905, 0.04761905, 0.        , 0.04761905, 0.14285714,
       0.14285714, 0.19047619, 0.23809524, 0.14285714, 0.33333333,
       0.14285714, 0.33333333, 0.38095238, 0.42857143, 0.57142857,
       0.57142857, 0.52380952, 0.61904762, 0.76190476, 0.57142857,
       0.76190476, 0.95238095, 0.80952381, 0.85714286, 1.        ,
       0.85714286, 0.95238095, 1.        , 1.        , 0.95238095,
       1.        , 0.95238095, 1.        ])

전체 21개의 문항을 만들었고 그 중에 다 맞으면 1, 다 틀리면 0이 된다.

 

여기서 방법이 특이한데, MSE를 사용하지 않는다.

로지스틱회귀에서 경사하강법을 사용하는데,

\(\frac{1}{1+e^{-(\beta_0 + \beta_1 \times x})}\)을 갖는다. 

이 회귀식 \(\beta_0 + \beta_1 \times x\)를 미분하면 경사하강법을 사용할 수 있다.

그 방법은 아래 블로그에 회귀계수 추정에 있다.

https://eunjin123123.tistory.com/605

 

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

class LogisticRegression():
  def __init__(self,lr=0.01,n_iters=1000):
    self.lr=lr
    self.n_iters=n_iters
    self.weight=None
    self.bias=None

  def fit(self,X,y):
    n_samples =X.shape[0]
    self.weight = 0
    self.bias = 0

    for _ in range(self.n_iters):
      linear_pred = np.dot(X,self.weight) + self.bias
      predictions=sigmoid(linear_pred)
      dw = (1/n_samples) * np.dot(X,predictions-y) 
      db = (1/n_samples) *np.sum(predictions-y) 

      self.weight = self.weight - self.lr* dw
      self.bias = self.bias - self.lr*db
  def predict(self,X,y):
    linear_pred = np.dot(X,self.weight) + self.bias
    predictions=sigmoid(linear_pred)

  def show(self):
    return self.weight, self.bias

간단히 로지스틱회귀를 돌려서 그 결과를 보자.

 

A=LogisticRegression(n_iters=100000)
A.fit(theta,P)
beta_1,beta_0=A.show()
b=- beta_0/beta_1
a = beta_1
a,b

(1.2241365711060053, -0.36016386636452385)

실제값과 꽤 근사하는 것 같다.

a는 그대로 \(\beta_1\)의 값을 쓴다.

b는 \(-\frac{\beta_0}{\beta_1}\)로 구한다.

 

이 구한 곡선이 실제로 잘 들어맞는지 알려면 카이스퀘어 적합도검정(chi-square goodness-of-fit)을 해야 한다.

우리가 맨날 보던 p값 등을 단순 로지스틱 회귀로는 구할 수 없으니 객관적 지표를 봐야한다.

 

$$ \chi^2 = \displaystyle{\sum_{g=1}^G}f_g \frac{[p(\theta_g) - P(\theta_g)]^2}{P(\theta_g)Q(\theta_g)}$$

 

G=능력의 그룹의 수이다.

\(theta_g\)= g그룹의 능력

\(f_g\) = 능력 \(theta_g\)를 시험자의 수

\(p(\theta_g)\)= 그룹 g에서 올바른 정답이 관측된 비율

\(P(\theta_g)\) = 모수 추정을 ㅌ오해 얻은 정답을 맞출 확률

\(Q(\theta_g) = 1-P(\theta_g)\)

 

위의 식을 통해 얻은 값이 너무 크다면 이 모델이 데이터에 잘 적합하지 못하는 것이다.

이런 경우는 두 가지 경우에 발생하는데

 

1) 적절치못한 문항 곡선 특성

2) 값들이 너무 넓게 퍼져있어서 어떤 모델을 써도 적합할 수 없는 경우

 

위의 데이터는 적합도가 29.97이 나왔다고 한다.

 

 

3.3 The Group Invariance of Item Parameters

문항의 모수는 응시자의 능력에 독립적이다. 

문항의 모수들은 그래서 집단 불변성(group invariance)을 갖고 있다.

 

 

그룹1의 데이터이다. 능력 -1부터 -3까지의 사람들이 모여있다.

이 곡선은 이 데이터에서도 적합해야한다.

 

 

그룹2는 능력1부터 능력3까지 모여있는데, 여기서도 이 값이 적합해야 한다.

이 둘을 이었을 때도 데이터가 선에 적합해야 한다.

 

 

3.4 Computer Session

 

theta = np.arange(-3,3.1875,0.1875)
f = np.repeat(21,len(theta))
wb = round(np.random.uniform(-3,3),2)
wa = round(np.random.uniform(0.2,2),2)
wc = round(np.random.uniform(0,0.35),2)
mdl=2
if(mdl==1 or mdl==2):
  wc=0
if(mdl ==1):
  wa=1
for g in range(1,len(theta)):
  P=wc+(1-wc)/(1+np.exp(-wa*(theta-wb)))
p = np.random.binomial(f,P)/f

plt.plot(theta,p)
plt.xlabel('Ability')
plt.ylabel('Probability of Correct response')
plt.xlim(-3,3)
plt.ylim(0,1)

 

 

그림을 그려보았는데 이게 맞나

 

cs=0
for g in range(1,len(theta)):
  v=f[g]*(p[g]-P[g])**2/ (P[g]*(1-P[g]))
  cs=cs+v
cs=np.round(cs,2)
if(mdl==1):
  maintext = f'Chi-square={cs} \n b={wb}'
if(mdl==2):
  maintext = f'Chi-square={cs} \n a={wa} b={wb}'
if(mdl==3):
  maintext = f'Chi-square={cs} \n a={wa}, b={wb}, c={wc}'

plt.plot(theta,P)
plt.xlabel('Ability')
plt.ylabel('Probability of Correct response')
plt.title(maintext)
plt.xlim(-3,3)
plt.ylim(0,1)

이 코드는 카이스퀘어 검정을 한 건데,

 

 

$$ \chi^2 = \displaystyle{\sum_{g=1}^G}f_g \frac{[p(\theta_g) - P(\theta_g)]^2}{P(\theta_g)Q(\theta_g)}$$

 

이건 위의 식을 그대로한 것이고 책에는 좀 다르게 나와있다.

cs=0
for g in range(1,len(theta)):
  v=f[g]*(p[g]-P[g])**2/ (P[g]-P[g]**2)
  cs=cs+v
cs=np.round(cs,2)
if(mdl==1):
  maintext = f'Chi-square={cs} \n b={wb}'
if(mdl==2):
  maintext = f'Chi-square={cs} \n a={wa} b={wb}'
if(mdl==3):
  maintext = f'Chi-square={cs} \n a={wa}, b={wb}, c={wc}'

plt.plot(theta,P)
plt.xlabel('Ability')
plt.ylabel('Probability of Correct response')
plt.title(maintext)
plt.xlim(-3,3)
plt.ylim(0,1)

둘다 값은 똑같이 나온다.

 

별로 값이 적합하지는 않는 것 같다.

 

딱봐도 망해보인다

3.4.2 집단 불변성 테스트

부분으로 나누어서 잘되었는지 보자.(당연히 잘 안되었을 것이다.)

tll =-3
tlu=-1
lowerg1=0
for g in range(1,len(theta)):
  if(theta[g]<tll):
    lowerg1 = lowerg1+1
upperg1=0
for g in range(1,len(theta)):
  if(theta[g]<tlu):
    upperg1=upperg1+1

theta1=theta[lowerg1:upperg1]
theta1
p1=p[lowerg1:upperg1]
p1

plt.plot(theta1,p1)
plt.xlabel('Ability')
plt.ylabel('Probability of Correct response')
plt.title(maintext)
plt.xlim(-3,3)
plt.ylim(0,1)

원래 적합해야할거보다 꽤 위로 가 있다.

 

다른 구간도 봐보자(볼 필요가 있나?).

t2l =1
t2u=3
lowerg1=0
for g in range(1,len(theta)):
  if(theta[g]<t2l):
    lowerg1 = lowerg1+1
upperg1=0
for g in range(1,len(theta)):
  if(theta[g]<t2u):
    upperg1=upperg1+1

theta2=theta[lowerg1:upperg1]
p2=p[lowerg1:upperg1]

plt.plot(theta2,p2)
plt.plot(theta,P)
plt.xlabel('Ability')
plt.ylabel('Probability of Correct response')
plt.title(maintext)
plt.xlim(-3,3)
plt.ylim(0,1)

 

3.4.3 함수화

def iccfit(mdl):
  theta = np.arange(-3,3.1875,0.1875)
  f=np.repeat(21,len(theta))
  wb=round(np.random.uniform(-3,3),2)
  wa=round(np.random.uniform(0.2,2.8),2)
  wc=round(np.random.uniform(0,0.35),2)

  if(mdl ==1or mdl==2 ):
    wc=0
  if(mdl==1):
    wa=1
  for g in range(len(theta)):
    P= wc + (1-wc) / (1+np.exp(-wa*(theta-wb)))
  p=np.random.binomial(f,P)/f

  plt.plot(theta,p,label=f'a={wa}, b={wb},c={wc}')
  plt.legend()
  plt.xlabel('Ability')
  plt.ylabel('Probability of Correct response')
  plt.xlim(-3,3)
  plt.ylim(0,1)
  
  cs=0
  for g in range(0,len(theta)):
    v=f[g]*(p[g]-P[g])**2 / (P[g]-P[g]**2)
    cs=cs+v
  cs=np.round(cs,2)
  if(mdl==1):
    maintext=f'Chi-square= {cs} \n b={wb}'
  if(mdl==2):
    maintext=f'Chi-square= {cs} \n a={wa} b={wb}'
  if(mdl==3):
    maintext=f'Chi-square= {cs} \n a={wa} b={wb} c={wc}'
  
  plt.plot(theta,P,'-')
  plt.xlim(-3,3)
  plt.ylim(0,1)
  plt.title(maintext)

iccfit(1)

iccfit(2)

이건 그나마 좀 적합도가 괜찮게 나온 거 같다.

 

iccfit(3)

얘는 뭐

 

2모수 모형이 가장 데이터에 잘적합하는 것 같다.

 

 

3.4.4 집단불변성 함수화

def groupinv(md1,t1l=-3,t1u=-1,t2l=1,t2u=3):
  theta = np.arange(-3,3.1875,0.1875)
  f=np.repeat(21,len(theta))
  wb=round(np.random.uniform(-3,3),2)
  wa=round(np.random.uniform(0.2,2.8),2)
  wc=round(np.random.uniform(0,0.35),2)

  if(mdl ==1or mdl==2 ):
    wc=0
  if(mdl==1):
    wa=1
  for g in range(len(theta)):
    P= wc + (1-wc) / (1+np.exp(-wa*(theta-wb)))
  p=np.random.binomial(f,P)/f
  
  lowerg1=0
  for g in range(len(theta)):
    if(theta[g]<=t1l):
      lowerg1=lowerg1+1
  upperg1=0
  for g in range(len(theta)):
    if(theta[g]<=t1u):
      upperg1=upperg1+1
  theta1=theta[lowerg1:upperg1]
  p1=p[lowerg1:upperg1]

  lowerg2=0
  for g in range(len(theta)):
    if(theta[g]<=t2l):
      lowerg2=lowerg2+1
  upperg2=0
  for g in range(len(theta)):
    if(theta[g]<=t2u):
      upperg2=upperg2+1
  theta2=theta[lowerg2:upperg2]
  p2=p[lowerg2:upperg2]
 
  
  plt.plot(theta1,p1,label=f'theta=1')
  plt.plot(theta2,p2,label=f'theta=2')
  plt.legend()
  plt.xlabel('Ability')
  plt.ylabel('Probability of Correct response')
  plt.yticks([0,0.5,1])
  plt.xlim(-3,3)
  plt.ylim(0,1)

  if(mdl==1):
    maintext=f'Pooled groups \n b={wb}'
  if(mdl==2):
    maintext=f'Pooled groups \n a={wa} b={wb}'
  if(mdl==3):
    maintext=f'Pooled groups \n a={wa} b={wb} c={wc}'

  plt.plot(theta,P)
  plt.legend()
  plt.xlim(-3,3)
  plt.ylim(0,1)
  plt.title(maintext)

groupinv(1)