Up

Three Aspects of Reliability in Item Response Theory

Precision, Consistency, Predictive Power

Yasuharu Okamoto, 2018

 

Three aspects of reliability, i.e., precision, consistency, and predictive power, in item response theory (IRT) can be estimated by simulation (Okamoto, 2017). Estimations by simulation of the three aspects are based on the following notation for IRT.

 

:   Test item i, which may have more than two response categories.

:   Used items for person p, i.e.,  denotes the j-th item from the test items s, which person p takes. In the case of a fixed set of test items,  is the constant set of test items, i.e., the set of all test items. In the case of an adaptive test (e.g., Okamoto, 2007),  depends on the performance of person p, so is considered to be determined stochastically depending on person pfs performance.

:    The true value of person pfs ability, which may be represented by a vector when a multifactor model is employed. In this study, a one-factor model is considered, because it is sufficient to use a simple model to show this studyfs essential point.

:     Estimate of . Popular methods of estimation are the maximum likelihood and the Bayesian. When the personfs responses are all the same extreme ones,  cannot be estimated by the maximum likelihood method without ad hoc assumptions. Hence, this study employs the Bayesian method.

 

The modelling above is represented by the following diagram.

Using diagram (1), three types of coefficients of reliability, precision, stability, and predictive power, are given as follows.

 

Precision.  Precision is the popular mathematical definition of coefficient of reliability. In IRT, this coefficient is given for the maximum likelihood (ML) method (Toyoda, 1989). For the Bayesian method, the coefficient of reliability is proposed as given by

where  and  are variances of post and prior distributions of  (Okamoto, 2015).  and  correspond to variances of observed scores and error terms in the case of CTT, respectively.

Taking expectation of , we have

A Python script for estimation of precision is given in the section below.

 

Stability.  Stability can be represented by correlation coefficient between two parallel tests. A sequence of a parallel test for sequence (1) can be represented by diagram (1f)

Person p takes two parallel tests, and results are denoted by  and , respectively. Corresponding to  and , two estimates  and  are given. Test stability is given by the correlation coefficient between   and .

Samejima (1994) discusses correlation (4) in the framework of the ML method.

A Python script for estimation of stability (or consistency) is given in the section below.

 

Predictive Power.  Prediction of the true value  from the estimate  can be represented as regression model

The predictive power of model (5) is given by the coefficient of determination  of model (5), which equals the square of the correlation coefficient between  and .

Model (5) represents how well the true value  can be inferred from the estimate . On the other hand, the following model (7) represents how strongly the estimate  is constrained by the true value .

The coefficient of determination of model (7) is the same as that of model (5). The coefficient  denotes strength of relation between the true value  and the estimate , which is a common value for models (5) and (7).

A Python script for estimation of predictive power is given in the section below.

 

A Basic Model of Item Response Theory for the Simulation

Consider  test items, item 1 to item , each having  response categories, 1 to . Denote a response to item  by . When the number of response categories is , . Set the following graded response model (Samejima, 1969) or Thurstonian model.

where

and

fs denote shifts of difficulty corresponding to categories 1 to  and satisfy the following condition.

 

For model (8), a Stan script EstTheta.stan to estimate the posterior distribution of  for responses  was developed as shown in listing 1.

 

Listing 1.  The Stan script EstTheta.stan for estimation of the posterior distribution of .

data {

    int K;

    int Nitm;

    real C[K - 1];

    real<lower = 0.0> a[Nitm];

    real b[Nitm];

    int<lower = 1, upper = K> X[Nitm];

}

parameters {

    real theta;

}

transformed parameters {

    vector[K] p[Nitm];

    for (i in 1:Nitm){

        p[i][K] = inv_logit(1.7 * a[i] * (theta - b[i] - C[K - 1]));

        p[i][1] = 1.0 - inv_logit(1.7 * a[i] * (theta - b[i] - C[1]));

        if (K > 2){

            for (k in 2:(K - 1)) {

                p[i][k] = inv_logit(1.7 * a[i] * (theta - b[i] - C[k - 1]))

                            - inv_logit(1.7 * a[i] * (theta - b[i] - C[k]));

            }

        }

    }

}

model {

    theta ~ normal(0.0, 1.0);

    for (i in 1:Nitm)

        X[i] ~ categorical(p[i]);

}

 

Using the Stan script EstTheta.stan, reliability coefficients for M items, to which model (8) are applied, are estimated by simulation as shown follows.

In the following examples, the following parameter values for Model (8) were used.

This is coded as following

Nitems = 5                                 #    Number of items

X = [1, 1, 1, 1, 1]                        #    Responses to be set in simulation

 

a = [1.0, 1.0, 1.0, 1.0, 1.0]              #    Discrimination parameters

b = [-0.5, -0.3, 0.0, 0.3, 0.5]            #    Difficulty parameters

K = 2                                      #    Number of response categories

C = [0.0]                                  #    Category boundaries

Of course, other parameter values can be used with the corresponding coding for variables, Nitems, X, a, b, K and C, in the Python script.

 

Estimation of Precision

Precision is estimated with diagram (1) and equations (2) and (3) by Python script in Listing 2.

Listing 2. Python script for estimation of precision

import pystan

from pystan import StanModel

import pickle

import random

import math

import matplotlib.pyplot as plt

 

def icc(t, a, b, c):

    return 1.0 / (1.0 + math.exp(-1.7 * a * (t - b - c)))

 

def CatResponse(t, a, b, c, K):

    res = 1

    v = random.random()

    while True:

        if v > icc(t, a, b, c[res - 1]):

            break

        res += 1

        if res >= K:

            break

    return res

 

def calcVar( L ):

    n = len(L)

    sum = 0.0

    for v in L:

        sum += v

    mean = sum / n

    ssum = 0.0

    for v in L:

        ssum += (mean - v) ** 2.0

    return ssum / (n - 1.0)

 

Nitems = 5                                 #    Number of items

X = [1, 1, 1, 1, 1]                        #    Responses to be set in simulation

 

a = [1.0, 1.0, 1.0, 1.0, 1.0]              #    Discrimination parameters

b = [-0.5, -0.3, 0.0, 0.3, 0.5]            #    Difficulty parameters

K = 2                                      #    Number of response categories

C = [0.0]                                  #    Category boundaries

#   K = 7

#   C = [-2.0, -1.0, -0.25, 0.25, 1.0, 2.0]

 

Data = {'K': K, 'Nitm': Nitems,  'C': C, 'a':a, 'b': b, 'X': X}

 

sm = StanModel(file = 'EstTheta.stan')

with open('model.pkl', 'wb') as f:

    pickle.dump(sm, f)

 

L_theta = []

L_precision = []

n_step = 1000

for step in range(n_step):

    if (step % 5 == 0):

        print("\nstep = {0}/{1}".format(step, n_step))

    t = random.gauss(0.0, 1.0)

    for i in range(Nitems):

        X[i] = CatResponse(t, a[i], b[i], C, K)      

       

    sm = pickle.load(open('model.pkl', 'rb'))

    fit = sm.sampling(data = Data,  n_jobs = 1)

    theta = fit['theta']

    var_theta = calcVar(theta)

 

    L_theta.append(t)

    L_precision.append(1.0 - var_theta)

 

sum = 0.0

for p in L_precision:

    sum += p

mean_precision = sum / len(L_precision)

 

 

plt.title("Mean precison = {0:>.3}".format(mean_precision))

plt.xlabel("theta")

plt.ylabel("precison")

plt.axis([-4.0, 4.0, 0.0, 1.0])

plt.plot(L_theta, L_precision, 'bo')

plt.show()

 

Run the script, then the results of estimation of precision will be shown graphically as in Figure 1.

Figure 1

 

Figure 1 is for items, which have the following parameter values for Model (8),

This is coded as following

Nitems = 5                                 #    Number of items

X = [1, 1, 1, 1, 1]                        #    Responses to be set in simulation

 

a = [1.0, 1.0, 1.0, 1.0, 1.0]              #    Discrimination parameters

b = [-0.5, -0.3, 0.0, 0.3, 0.5]            #    Difficulty parameters

K = 2                                      #    Number of response categories

C = [0.0]                                  #    Category boundaries

 

When the number of categories is increased from 2 to 7, in this case the above code concerning to K and C are changed to, e.g.,

K = 7

C = [-2.0, -1.0, -0.25, 0.25, 1.0, 2.0]

we will get the results shown in Figure 2.

Figure 2

 

The python script and the Stan script are archived in a zip file Precision.zip, which can be downloaded.

 

Estimation of Stability (or Consistency)

Stability (or Consistency) can be estimated with diagrams (1) and (1f), and equation (4) using a Python script shown in Listing 3.

 

Listing 3. Python script for estimation of consistency.

import pystan

from pystan import StanModel

import pickle

import random

import math

import matplotlib.pyplot as plt

 

def icc(t, a, b, c):

    return 1.0 / (1.0 + math.exp(-1.7 * a * (t - b - c)))

 

def CatResponse(t, a, b, c, K):

    res = 1

    v = random.random()

    while True:

        if v > icc(t, a, b, c[res - 1]):

            break

        res += 1

        if res >= K:

            break

    return res

 

def calcCor(data1, data2):

    sum1 = 0.0

    sum2 = 0.0

    n = len(data1)

    for i in range(n):

        sum1 += data1[i]

        sum2 += data2[i]

    mean1 = sum1 / n

    mean2 = sum2 / n

    ssum1 = 0.0

    ssum2 = 0.0

    xsum = 0.0

    for i in range(n):

        ssum1 += (data1[i] - mean1) ** 2

        ssum2 += (data2[i] - mean2) ** 2

        xsum += (data1[i] - mean1) * (data2[i])

    return xsum / ((ssum1 * ssum2) ** 0.5)

 

Nitems = 5                         #    Number of items

X = [1, 1, 1, 1, 1]                #    Responses to be set in simulation

 

a = [1.0, 1.0, 1.0, 1.0, 1.0]      #    Discrimination parameters

b = [-0.5, -0.3, 0.0, 0.3, 0.5]    #    Difficulty parameters

K = 2                              #    Number of response categories

C = [0.0]                          #    Category boundaries

#   K = 7

#   C = [-2.0, -1.0, -0.25, 0.25, 1.0, 2.0]

 

Data = {'K': K, 'Nitm': Nitems,  'C': C, 'a':a, 'b': b, 'X': X}

 

sm = StanModel(file = 'EstTheta.stan')

with open('model.pkl', 'wb') as f:

    pickle.dump(sm, f)

 

L_theta1 = []

L_theta2 = []

 

n_step = 1000

for step in range(n_step):

    if (step % 5 == 0):

        print('\nstep = {0}/{1}'.format(step, n_step))

    t = random.gauss(0.0, 1.0)

    for i in range(Nitems):

        X[i] = CatResponse(t, a[i], b[i], C, K)

    sm = pickle.load(open('model.pkl', 'rb'))

    fit = sm.sampling(data = Data,  n_jobs = 1)

    theta = fit['theta']

    theta.sort()

    L_theta1.append(theta[int(len(theta) / 2)])

 

    for i in range(Nitems):

        X[i] = CatResponse(t, a[i], b[i], C, K)

    sm = pickle.load(open('model.pkl', 'rb'))

    fit = sm.sampling(data = Data,  n_jobs = 1)

    theta = fit['theta']

    theta.sort()

    L_theta2.append(theta[int(len(theta) / 2)])

   

r = calcCor(L_theta1, L_theta2)

 

plt.title("Cor(theta1, theta2) = {0:>.3}".format(r))

plt.xlabel("theta1")             

plt.ylabel("theta2")              

plt.plot(L_theta1, L_theta2, 'bo')   

plt.show()

 

Run the above script, then the results is shown graphically as in Figure 3.

Figure 3

 

Figure 3 is for items, which have the following parameter values for Model (8),

This is coded as following

Nitems = 5                                 #    Number of items

X = [1, 1, 1, 1, 1]                        #    Responses to be set in simulation

 

a = [1.0, 1.0, 1.0, 1.0, 1.0]              #    Discrimination parameters

b = [-0.5, -0.3, 0.0, 0.3, 0.5]            #    Difficulty parameters

K = 2                                      #    Number of response categories

C = [0.0]                                  #    Category boundaries

 

When the number of categories is increased from 2 to 7, in this case the code concerning to the number of categories K and category boundaries C are changed to, e.g.,

K = 7

C = [-2.0, -1.0, -0.25, 0.25, 1.0, 2.0]

we get the results as shown in Figure 4.

Figure 4

 

The Python script and the Stan script are archived in a zip file Consistency.zip, which can be downloaded.

 

Estimation of Predictive Power

Predictive power can be estimated with diagram (1), and equations (6) and (7), using a Python script in Listing 4.

 

Listing 4. Python script for estimation of predictive power.

import pystan

from pystan import StanModel

import pickle

import random

import math

import matplotlib.pyplot as plt

 

def icc(t, a, b, c):

    return 1.0 / (1.0 + math.exp(-1.7 * a * (t - b - c)))

 

def CatResponse(t, a, b, c, K):

    res = 1

    v = random.random()

    while True:

        if v > icc(t, a, b, c[res - 1]):

            break

        res += 1

        if res >= K:

            break

    return res

 

def calcCoeff(x, y):

    sumx = 0.0

    sumy = 0.0

    n = len(x)

    for i in range(n):

        sumx += x[i]

        sumy += y[i]

    meanx = sumx / n

    meany = sumy / n

    sumxx = 0.0

    sumyy = 0.0

    sumxy = 0.0

    for i in range(n):

        sumxx += (x[i] - meanx) ** 2

        sumyy += (y[i] - meany) ** 2

        sumxy += (x[i] - meanx) * (y[i] - meany)

    Cxy = sumxy / n

    Sx = (sumxx / n) ** 0.5

    Sy = (sumyy / n) ** 0.5

    r = Cxy / (Sx * Sy)

    a = r * Sy / Sx

    b = -Cxy * meanx / (Sx ** 2) + meany

       

    return r, a, b

 

Nitems = 5                         #    Number of items

X = [1, 1, 1, 1, 1]                #    Responses to be set in simulation

 

a = [1.0, 1.0, 1.0, 1.0, 1.0]      #    Discrimination parameters

b = [-0.5, -0.3, 0.0, 0.3, 0.5]    #    Difficulty parameters

K = 2                              #    Number of response categories

C = [0.0]                          #    Category boundaries

#   K = 7

#   C = [-2.0, -1.0, -0.25, 0.25, 1.0, 2.0]

 

Data = {'K': K, 'Nitm': Nitems,  'C': C, 'a':a, 'b': b, 'X': X}

 

sm = StanModel(file = 'EstTheta.stan')

with open('model.pkl', 'wb') as f:

    pickle.dump(sm, f)

 

L_theta1 = []

L_thetaTrue = []

n_step = 1000

for step in range(n_step):

    if (step % 5 == 0):

        print('\nstep = {0}/{1}'.format(step, n_step))

    t = random.gauss(0.0, 1.0)

    L_thetaTrue.append(t)

    for i in range(Nitems):

        X[i] = CatResponse(t, a[i], b[i], C, K)

    sm = pickle.load(open('model.pkl', 'rb'))

    fit = sm.sampling(data = Data,  n_jobs = 1)

    theta = fit['theta']

    theta.sort()

    L_theta1.append(theta[int(len(theta) / 2)])

 

r, coeff_a, coeff_b = calcCoeff(L_thetaTrue, L_theta1)

 

min_t = L_thetaTrue[0]

max_t = min_t

n = len(L_thetaTrue)

for i in range(n):

    if min_t > L_thetaTrue[i]:

        min_t = L_thetaTrue[i]

    if max_t < L_thetaTrue[i]:

        max_t = L_thetaTrue[i]

 

regression_x = [min_t, max_t]

regression_y = [coeff_a * min_t + coeff_b, coeff_a * max_t + coeff_b]

plt.plot(regression_x, regression_y, 'b-')

 

plt.title("Predictive power = {0:>.3}\nr = {1:>.3}".format(r**2, r))

plt.xlabel("True theta")   

plt.ylabel("Est. theta")    

plt.plot(L_thetaTrue, L_theta1, 'bo')   

plt.show() 

 

Run the above script, then the results is shown graphically as shown in Figure 5.

Figure 5

 

Figure 5 is for items, which have the following parameter values for Model (8),

This is coded as following

Nitems = 5                                 #    Number of items

X = [1, 1, 1, 1, 1]                        #    Responses to be set in simulation

 

a = [1.0, 1.0, 1.0, 1.0, 1.0]              #    Discrimination parameters

b = [-0.5, -0.3, 0.0, 0.3, 0.5]            #    Difficulty parameters

K = 2                                      #    Number of response categories

C = [0.0]                                  #    Category boundaries

 

When the number of categories is increased from 2 to 7, in this case the code for the number of categories K and the category boundaries are changed to, e.g.,

K = 7

C = [-2.0, -1.0, -0.25, 0.25, 1.0, 2.0]

we get the following results as shown in Figure 6.

Figure 6

 

The files of the Python script and the Stan script are archived in a zip file Predictive.zip, which can be downloaded.

 

 

Reference

Okamoto, Y. (2017).  Three Aspects of Reliability in Item Response Theory. Proceedings of the 15th Annual Conference of the Japan Association for Research on TestingA152|155

Samejima, F. (1994). Estimation of reliability coefficients using the test information function and its modifications. Applied Psychological Measurement, 18, 229-244.

 

Up