本文為閱讀Baker與Kim（2017）《The Basics of Item Response Theory using R》第五章部份之筆記，本文聚焦於估計能力參數（ability parameter）之R語言程式嗎。
設定ability 作為函式，其中參數有mdl （模型，即Rasch、2PL或3PL）；u 為反應向量（由0與1構成的向量，0表該題錯誤，1表該題正確）；b、a 與 c 表試題參數向量。為避免能力值（theta）在試題全對（即 u 向量全為1）或全錯（即 u 向量全為0）情況下使能力值之估計值出現無窮大之情況，我們設定其上下限，其中 J 為總題數，
th= -log(2*J) #theta hat lower bound
th= log(2*J) #theta hat upper bound
通過最大概似估計法（MLE）估計能力值，所得之一階條件（f.o.c）利用牛頓法迭代求數值解。我們採用 log(x/(J-x)) 作為迭代的起始值，其中 J 為總題數，x 為 u 向量內元素之總和（答對題數）。我們設定最大的迭代次數（S）為10，收斂標準為前後兩次迭代數值之差小於.001。
th =log(x/(J-x)) #starting value
S = 10 # the maximum number of iterations
ccrit=.001 #the converge criterion of Newton's method
印出每次之迭代之能力值，
cat(paste("th=",th,"\n"));flush.console()
待達成迭代終止條件後，計算其標準誤（standard error, SE），
cat(paste("se=",se,"\n"));flush.console()
我們可以簡單執行一個五題二參數模型之的範例，
u=c(1,0,1,1,0)
b=c(-1,0,1,0,1)
a=c(1,1.2,0.8,0.7,0.2)
ability(mdl=2,u=u,b=b,a=a)
得到的能力值估計值為0.6079441，其標準誤為1.1638076。
附上課本之完整 R 程式碼。
ability = function(mdl,u,b,a,c){
J = length(b)
if(mdl==1|mdl==2|missing(c)){
c=rep(0,J)
}
if(mdl==1|missing(a)){a=rep(1,J)}
x=sum(u) #total score(a summation of u_j)
if(x==0){
th= -log(2*J)
}
if(x==J){
th= log(2*J)
}
if(x==0|x==J){
sumdem=0
for(j in 1:J){
pstar = 1/(1+exp(-a[j]*(th-b[j])))
phat = c[j]+(1-c[j])*pstar
sumdem = sumdem - a[j]**2 * phat*(1-phat)*(pstar/phat)**2 #3PL formula
}
se=1/sqrt(-sumdem)
}
if(x!=0 & x != J){
th =log(x/(J-x))
S = 10
ccrit=.001
for(s in 1:S){
sumnum = 0
sumdem = 0
for(j in 1:J){
pstar = 1/(1+exp(-a[j]*(th-b[j])))
phat = c[j]+(1-c[j])*pstar
sumnum = sumnum + a[j]*(u[j]-phat)*(pstar/phat)
sumdem = sumdem - a[j]**2*phat*(1-phat)*(pstar/phat)**2
}
delta = sumnum/sumdem
th=th-delta
if(abs(delta) < ccrit |s==S){
se = 1/sqrt(-sumdem)
break
}
cat(paste("th=",th,"\n"));flush.console()
}
}
cat(paste("th=",th,"\n"));flush.console()
cat(paste("se=",se,"\n"));flush.console()
thse=c(th,se)
return(thse)
}
參考資料：
Baker, F. B., & Kim, S. H. (2017). The Basics of Item Response Theory Using R (Vol. 10). New York: Springer.
