项目反应理论参数的EM估计
写在前面:
本文主要描述了整个IRT使用EM算法参数的估计过程,其中涉及大量公式,如只是需要了解IRT相关基础知识,请转战wiki~~
预警: 大量公式来袭~~
IRT(项目反应理论)广泛应应用于心理测量学,相比CTT(传统测量理论)的主要优势在于不依赖于样本,能够为被试者提供具有一致性的测量结果。
IRT的基本假设:
- 单维性,即能力是单维度的。
- 局部独立性, 即项目间局部独立。
- 项目反应函数假设,即被试者对项目对反应符合项目反应方程。
IRT 模型:
常见的IRT模型包括2PL模型和3PL模型,其中2PL模型表达式如下:
其中 θ i \theta _ { i } θi表示被试者的能力, a j a _ { j } aj是项目的区分度, b j b _ { j } bj是项目的难度,D=1.7。
参数估计:
对于IRT参数估计有多种方法,如对于项目参数的估计有:边际极大似然估计(MMLE),极大边际后验法(MMAP),EM,MCMC等方法,对于能力参数的估计有极大似然(MLE),贝叶斯众数法(MAP),贝叶斯后验期望估计法(EAP)等。本文重点介绍对项目参数的EM估计。
建模:
- 被观测数据 Y \mathbf { Y } Y是N*J维的,N表示有N个被试者,J表示有J个题目, Y = ( y 1 , y 2 , … , y N ) \mathbf { Y } = \left( \mathbf { y } _ { 1 } , \mathbf { y } _ { 2 } , \dots , \mathbf { y } _ { N } \right) Y=(y1,y2,…,yN), y i = ( y i 1 , y i 2 , … , y i J ) \mathbf { y } _ { i }=\left( y _ { i 1 } , y _ { i 2 } , \ldots , y _ { i J } \right) yi=(yi1,yi2,…,yiJ),其中 y i j y _ { i j } yij表示第i个被试者回答第j道题目的答案;
- 不可观测变量为 θ = ( θ 1 , θ 2 , … , θ N ) \boldsymbol { \theta } = \left( \theta _ { 1 } , \theta _ { 2 } , \ldots , \theta _ { N } \right) θ=(θ1,θ2,…,θN) 其中 θ i \theta _ { i } θi是第i个被试者的能力值,模型中假定 θ i \theta _ { i } θi只能取离散的值 q k , k = 1 , … , K q _ { k } , k = 1 , \dots , K qk,k=1,…,K,对应取 q k q _ { k } qk的概率为 π k , k = 1 , … , K \pi _ { k } , k = 1 , \ldots , K πk,k=1,…,K,即 θ i \theta _ { i } θi为多项式分布,其概率为: π = ( π 1 , π 2 , … , π K ) \boldsymbol {\pi} = \left( \pi _ { 1 } , \pi _ { 2 } , \dots , \pi _ { K } \right) π=(π1,π2,…,πK);不可观测变量为 θ = ( θ 1 , θ 2 , … , θ N ) \boldsymbol { \theta } = \left( \theta _ { 1 } , \theta _ { 2 } , \ldots , \theta _ { N } \right) θ=(θ1,θ2,…,θN) 其中 θ i \theta _ { i } θi是第i个被试者的能力值,模型中假定 θ i \theta _ { i } θi只能取离散的值 q k , k = 1 , … , K q _ { k } , k = 1 , \dots , K qk,k=1,…,K,对应取 q k q _ { k } qk的概率为 π k , k = 1 , … , K \pi _ { k } , k = 1 , \ldots , K πk,k=1,…,K,即 θ i \theta _ { i } θi为多项式分布,其概率为: π = ( π 1 , π 2 , … , π K ) \boldsymbol {\pi} = \left( \pi _ { 1 } , \pi _ { 2 } , \dots , \pi _ { K } \right) π=(π1,π2,…,πK);
- 完全数据为 [ ( y 1 , θ 1 ) , ( y 2 , θ 2 ) , … , ( y N , θ N ) ] \left[ \left( \mathbf { y } _ { 1 } , \theta _ { 1 } \right) , \left( \mathbf { y } _ { 2 } , \theta _ { 2 } \right) , \ldots , \left( \mathbf { y } _ { N } , \theta _ { N } \right) \right] [(y1,θ1),(y2,θ2),…,(yN,θN)];完全数据为 [ ( y 1 , θ 1 ) , ( y 2 , θ 2 ) , … , ( y N , θ N ) ] \left[ \left( \mathbf { y } _ { 1 } , \theta _ { 1 } \right) , \left( \mathbf { y } _ { 2 } , \theta _ { 2 } \right) , \ldots , \left( \mathbf { y } _ { N } , \theta _ { N } \right) \right] [(y1,θ1),(y2,θ2),…,(yN,θN)];
- 待估计参数为项目参数 Δ = ( δ 1 , δ 2 , . . . δ j ) \boldsymbol {\Delta}=(\boldsymbol { \delta } _ { 1 },\boldsymbol { \delta } _ { 2 },...\boldsymbol { \delta } _ { j }) Δ=(δ1,δ2,...δj),其中 δ j \boldsymbol { \delta } _ { j } δj为第j道题的参数。待估计参数为项目参数 Δ = ( δ 1 , δ 2 , . . . δ j ) \boldsymbol {\Delta}=(\boldsymbol { \delta } _ { 1 },\boldsymbol { \delta } _ { 2 },...\boldsymbol { \delta } _ { j }) Δ=(δ1,δ2,...δj),其中 δ j \boldsymbol { \delta } _ { j } δj为第j道题的参数。
EM算法应用于IRT中是迭代估计
Δ
\boldsymbol {\Delta}
Δ 和
π
\boldsymbol {\pi}
π:
E step: 根据给定的缺失数据的分布,观察数据和参数初始值,求完全数据的对数似然函数的条件期望。
M step: 极大化E-step给出的完全数据的对数似然函数的条件期望,求参数的值。
似然函数
给定观测变量 y = ( y 1 , y 2 , … , y J ) \mathbf { y } = \left( y _ { 1 } , y _ { 2 } , \ldots , y _ { J } \right) y=(y1,y2,…,yJ)并假定能力参数分布为: π = ( π 1 , π 2 , … , π J ) \pi = \left( \pi _ { 1 } , \pi _ { 2 } , \dots , \pi _ { J } \right) π=(π1,π2,…,πJ),项目参数为 Δ = ( δ 1 , δ 2 , . . . δ j ) \boldsymbol {\Delta}=(\boldsymbol { \delta } _ { 1 },\boldsymbol { \delta } _ { 2 },...\boldsymbol { \delta } _ { j }) Δ=(δ1,δ2,...δj),则观测变量的条件概率分布为:
由IRT的局部独立性假设,由于给定能力变量下项目反应结果是相互独立的,可得:
其中 P ( q k ∣ δ j ) P \left( q _ { k } | \boldsymbol { \delta } _ { j } \right) P(qk∣δj)是题目j的项目特征曲线,描述了给定能力下项目j答对的概率。由上式可以推出,完全数据的似然函数为: L ( Y , θ ∣ Δ , π ) = ∏ i = 1 N ∏ j = 1 J P ( θ i ∣ δ j ) y i j [ 1 − P ( θ i ∣ δ j ) ] 1 − y i j f ( θ i ∣ π ) = ∏ j = 1 J ∏ i = 1 N P ( θ i ∣ δ j ) y i j [ 1 − P ( θ i ∣ δ j ) ] 1 − y i j f ( θ i ∣ π ) = ∏ j = 1 J ∏ k = 1 K P ( θ k ∣ δ j ) r j k [ 1 − P ( q k ∣ δ j ) ] n k − r j k π k n k \begin{aligned} L ( \mathbf { Y } , \boldsymbol { \theta } | \mathbf { \Delta } , \boldsymbol { \pi } ) & = \prod _ { i = 1 } ^ { N } \prod _ { j = 1 } ^ { J } P \left( \theta _ { i } | \boldsymbol { \delta } _ { j } \right) ^ { y _ { i j } } \left[ 1 - P \left( \theta _ { i } | \boldsymbol { \delta } _ { j } \right) \right] ^ { 1 - y _ { i j } } f \left( \theta _ { i } | \boldsymbol { \pi } \right) \\ & = \prod _ { j = 1 } ^ { J } \prod _ { i = 1 } ^ { N } P \left( \theta _ { i } | \boldsymbol { \delta } _ { j } \right) ^ { y _ { i j } } \left[ 1 - P \left( \theta _ { i } | \boldsymbol { \delta } _ { j } \right) \right] ^ { 1 - y _ { i j } } f \left( \theta _ { i } | \boldsymbol { \pi } \right) \\ & = \prod _ { j = 1 } ^ { J } \prod _ { k = 1 } ^ { K } P \left( \theta _ { k } | \boldsymbol { \delta } _ { j } \right) ^ { r _ { j k } } \left[ 1 - P \left( q _ { k } | \boldsymbol { \delta } _ { j } \right) \right] ^ { n _ { k } - r _ { j k } } \pi _ { k } ^ { n _ { k } } \end{aligned} L(Y,θ∣Δ,π)=i=1∏Nj=1∏JP(θi∣δj)yij[1−P(θi∣δj)]1−yijf(θi∣π)=j=1∏Ji=1∏NP(θi∣δj)yij[1−P(θi∣δj)]1−yijf(θi∣π)=j=1∏Jk=1∏KP(θk∣δj)rjk[1−P(qk∣δj)]nk−rjkπknk
其中, f ( θ i ∣ π ) = π k f \left( \theta _ { i } | \boldsymbol { \pi } \right) = \pi _ { k } f(θi∣π)=πk if θ i = q k \theta _ { i } = q _ { k } θi=qk, n k n _ { k } nk是N个被试者中能力为 q k q _ { k } qk的人数, r j k \boldsymbol { r }_{ j k} rjk则是 n k n _ { k } nk人中答题正确的人数。
为方便计算,将上式化成对数似然函数为:
E-Step:
E步估计隐变量
θ
\boldsymbol { \theta}
θ,在给定项目参数
Δ
(
s
)
\boldsymbol {\Delta}^{(s)}
Δ(s) 和能力分布
π
(
s
)
\boldsymbol {\pi}^{(s)}
π(s)(由M步迭代估计得到)的和观测变量下隐变量的条件分布为:
f
(
q
k
∣
y
i
,
Δ
(
s
)
,
π
(
s
)
)
=
f
(
y
i
∣
q
k
,
Δ
(
s
)
)
π
k
(
s
)
∑
k
′
=
1
K
f
(
y
i
∣
q
k
′
,
Δ
(
s
)
)
π
k
′
(
s
)
=
π
k
(
s
)
∏
j
=
1
J
P
(
q
k
∣
δ
j
(
s
)
)
y
i
j
[
1
−
P
(
q
k
∣
δ
j
(
s
)
)
]
1
−
y
i
j
∑
k
′
=
1
K
π
k
′
(
s
)
∏
j
=
1
J
P
(
q
k
′
∣
δ
j
(
s
)
)
y
i
j
[
1
−
P
(
q
k
′
∣
δ
j
(
s
)
)
]
1
−
y
i
j
\begin{aligned} f \left( q _ { k } | \mathbf { y } _ { i } , \Delta ^ { ( s ) } , \pi ^ { ( s ) } \right) & = \frac { f \left( \mathbf { y } _ { i } | q _ { k } , \mathbf { \Delta } ^ { ( s ) } \right) \pi _ { k } ^ { ( s ) } } { \sum _ { k ^ { \prime } = 1 } ^ { K } f \left( \mathbf { y } _ { i } | q _ { k ^ { \prime } } , \Delta ^ { ( s ) } \right) \pi _ { k ^ { \prime } } ^ { ( s ) } } \\ & = \frac { \pi _ { k } ^ { ( s ) } \prod _ { j = 1 } ^ { J } P \left( q _ { k } | \boldsymbol { \delta } _ { j } ^ { ( s ) } \right) ^ { y _ { i j } } \left[ 1 - P \left( q _ { k } | \boldsymbol { \delta } _ { j } ^ { ( s ) } \right) \right] ^ { 1 - y _ { i j } } } { \sum _ { k ^ { \prime } = 1 } ^ { K } \pi _ { k ^ { \prime } } ^ { ( s ) } \prod _ { j = 1 } ^ { J } P \left( q _ { k ^ { \prime } } | \boldsymbol { \delta } _ { j } ^ { ( s ) } \right) ^ { y _ { i j } } \left[ 1 - P \left( q _ { k ^ { \prime } } | \boldsymbol { \delta } _ { j } ^ { ( s ) } \right) \right] ^ { 1 - y _ { i j } } } \end{aligned}
f(qk∣yi,Δ(s),π(s))=∑k′=1Kf(yi∣qk′,Δ(s))πk′(s)f(yi∣qk,Δ(s))πk(s)=∑k′=1Kπk′(s)∏j=1JP(qk′∣δj(s))yij[1−P(qk′∣δj(s))]1−yijπk(s)∏j=1JP(qk∣δj(s))yij[1−P(qk∣δj(s))]1−yij
n
k
(
s
)
=
E
(
n
k
∣
Y
,
Δ
(
s
)
,
π
(
s
)
)
=
∑
i
=
1
N
f
(
q
k
∣
y
i
,
Δ
(
s
)
,
π
(
s
)
)
=
∑
i
=
1
N
f
(
y
i
∣
q
k
,
Δ
(
s
)
)
π
k
(
s
)
∑
k
′
=
1
K
f
(
y
i
∣
q
k
′
,
Δ
(
s
)
)
π
k
′
(
s
)
=
∑
i
=
1
N
π
k
(
s
)
∏
j
=
1
J
P
(
q
k
∣
δ
j
(
s
)
)
y
i
j
[
1
−
P
(
q
k
∣
δ
j
(
s
)
)
]
1
−
y
i
j
∑
k
′
=
1
K
π
k
′
(
s
)
∏
j
=
1
J
P
(
q
k
′
∣
δ
j
(
s
)
)
y
i
j
[
1
−
P
(
q
k
′
∣
δ
j
(
s
)
)
]
1
−
y
i
j
n _ { k } ^ { ( s ) } = E \left( n _ { k } | \mathbf { Y } , \Delta ^ { ( s ) } , \boldsymbol { \pi } ^ { ( s ) } \right)= \sum _ { i = 1 } ^ { N } f \left( q _ { k } | \mathbf { y } _ { i } , \mathbf { \Delta } ^ { ( s ) } , \boldsymbol { \pi } ^ { ( s ) } \right)= \sum _ { i = 1 } ^ { N } \frac { f \left( \mathbf { y } _ { i } | q _ { k } , \mathbf { \Delta } ^ { ( s ) } \right) \pi _ { k } ^ { ( s ) } } { \sum _ { k ^ { \prime } = 1 } ^ { K } f \left( \mathbf { y } _ { i } | q _ { k ^ { \prime } } , \boldsymbol { \Delta } ^ { ( s ) } \right) \pi _ { k ^ { \prime } } ^ { ( s ) } }= \sum _ { i = 1 } ^ { N } \frac { \pi _ { k } ^ { ( s ) } \prod _ { j = 1 } ^ { J } P \left( q _ { k } | \boldsymbol { \delta } _ { j } ^ { ( s ) } \right) ^ { y _ { i j } } \left[ 1 - P \left( q _ { k } | \boldsymbol { \delta } _ { j } ^ { ( s ) } \right) \right] ^ { 1 - y _ { i j } } } { \sum _ { k ^ { \prime } = 1 } ^ { K } \pi _ { k ^ { \prime } } ^ { ( s ) } \prod _ { j = 1 } ^ { J } P \left( q _ { k ^ { \prime } } | \boldsymbol { \delta } _ { j } ^ { ( s ) } \right) ^ { y _ { i j } } \left[ 1 - P \left( q _ { k ^ { \prime } } | \boldsymbol { \delta } _ { j } ^ { ( s ) } \right) \right] ^ { 1 - y _ { i j } } }
nk(s)=E(nk∣Y,Δ(s),π(s))=∑i=1Nf(qk∣yi,Δ(s),π(s))=∑i=1N∑k′=1Kf(yi∣qk′,Δ(s))πk′(s)f(yi∣qk,Δ(s))πk(s)=∑i=1N∑k′=1Kπk′(s)∏j=1JP(qk′∣δj(s))yij[1−P(qk′∣δj(s))]1−yijπk(s)∏j=1JP(qk∣δj(s))yij[1−P(qk∣δj(s))]1−yij
r
j
k
(
s
)
=
y
i
j
n
k
(
s
)
r _ { j k } ^ { ( s ) } = {y}_{ij}n _ { k } ^ { ( s ) }
rjk(s)=yijnk(s)
M-Step:
M步是使用E步计算的 n k ( s ) n _ { k } ^ { ( s ) } nk(s)和 r j k ( s ) r _ { j k } ^ { ( s ) } rjk(s)去迭代估计 Δ , π \Delta , \pi Δ,π。根据上述的对数似然函数可得:
其中
l ( π ) = ∑ k = 1 K n k ( s ) log [ π k ] l ( \boldsymbol { \pi } ) = \sum _ { k = 1 } ^ { K } n _ { k } ^ { ( s ) } \log \left[ \pi _ { k } \right] l(π)=∑k=1Knk(s)log[πk]
使用极大似然估计求 Δ , π \Delta , \pi Δ,π,由于 π \pi π是服从多项式分布的,所以很容易得到
为求 Δ \Delta Δ,对 l ( δ j ) l \left( \boldsymbol { \delta } _ { j } \right) l(δj)关于 δ j \boldsymbol { \delta } _ { j } δj求导:
∑ k = 1 K r j k ( s ) − n k ( s ) P ( q k ∣ δ j ) [ 1 − P ( q k ∣ δ j ) ] P ( q k ∣ δ j ) ∂ P ( q k ∣ δ j ) ∂ δ t j = 0 \sum _ { k = 1 } ^ { K } \frac { r _ { j k } ^ { ( s ) } - n _ { k } ^ { ( s ) } P \left( q _ { k } | \boldsymbol { \delta } _ { j } \right) } { \left[ 1 - P \left( q _ { k } | \boldsymbol { \delta } _ { j } \right) \right] P \left( q _ { k } | \boldsymbol { \delta } _ { j } \right) } \frac { \partial P \left( q _ { k } | \boldsymbol { \delta } _ { j } \right) } { \partial \delta _ { t j } } = 0 ∑k=1K[1−P(qk∣δj)]P(qk∣δj)rjk(s)−nk(s)P(qk∣δj)∂δtj∂P(qk∣δj)=0
结果:
使用12000名学生的约100万条做题记录进行IRT参数估计,其中涉及题目约1000道左右,首先使用EM算法对题目参数进行估计,然后使用EAP对学生能力进行估计,结果分别如下:
上图中,ability为学生能力对评估结果,acc为学生实际答题的准确率,length为做题记录的长度,从图中可以看出,估计出的学生能力与实际答题的准确率正相关,且估计出的学生能力分布呈正态。相同做题准确率下,不同做题长度对评估结果也有影响。
上图为估计的题目参数,a表示题目的区分度,b表示题目的难度,可以看出题目难度和题目的准确率严格正相关,区分度和准确率呈倒U型。
参考文献:
1. http://www.openirt.com/b-a-h/papers/note9801.pdf IRT Parameter Estimation using the EM Algorithm
2. https://zhuanlan.zhihu.com/p/29887184 Python与项目反应理论:基于EM和MCMC的参数估计算法
[3]. https://www.cnblogs.com/jiangxinyang/p/9621997.html IRT模型的参数估计方法(EM算法和MCMC算法)
[4]. https://www.john-uebersax.com/stat/irt2.htm How to Calculate Expected A Posteriori (EAP) Scores for Unidimensional and Multidimensional Latent Trait Models
[5]. https://en.wikipedia.org/wiki/Item_response_theory Item response theory