本文為閱讀Baker與Kim（2017）《The Basics of Item Response Theory using R》第五章部份之筆記，本文聚焦模擬能力參數（ability parameter），想藉此說明能力參數的模擬與估計值變動性。
這一段主要描述一些模擬設定，分別為模型為2PL、此人能力值為0.5、共五題之難度參數（b）與鑑別力參數（a）。通過迴圈計算理論之答題正確率（P），再計算理論之標準誤（theroatical standard error, set）。
mdl = 2
theta = 0.5
b = c(-0.5, -0.25, 0.0, 0.25, 0.5)
a = c(1, 1.5, 0.7, 0.6, 1.8)
J = length(b)
if(mdl == 1| mdl == 2){ c = rep(0,J)}
if(mdl == 1){a = rep(1,J)}
sumdemt = 0
for(j in 1:J){
Pstar = 1/(1+exp(-a[j]*(theta-b[j])))
P=c[j]+(1-c[j]) *Pstar
sumdemt = sumdemt-a[j]**2*P*(1.0 - P)*(Pstar/P)**2
}
set=1/sqrt(-sumdemt)#SE-theoretical
此段這定模擬情境，模擬次數（R）為10。宣告thr（theta hat from the replication）與ser（standard error from the replication），記錄每次模擬之能力估計值與其標準誤。u 向量為由0與1組成之反應向量（0表該題錯誤；1表該題正確）。
我們可以藉由一開始對題目參數之設定計算出在能力為0.5之下，每題的理論答題正確機率（P）。其中，U服從二項分配（binominal distribution），以各題之理論答題正確機率（P）作為母體參數產生實現值（0、1）模擬答題反應，形成 u 向量。
最後再使用 ability()，以最大概似估計（maximum likelihood estimation）與牛頓法（Newton's method）估計能力估計值。ability()會回傳兩個值，第一個為能力估計值，第二個為標準誤，再將此二值分別儲存於thr 與 ser。
ability() 見【試題反應理論】Ch5 能力值之估計（R code）。
R=10 # the number of replications
thr=rep(0,R) # theta hat from the replication
ser=rep(0,R) # standard error from the replication
for(r in 1:R){
u=rep(0,J) #response vector u is generated from binominal distribution
for(j in 1:J){
P=c[j]+(1-c[j])/(1+exp(-a[j]*(theta-b[j])))
u[j]=rbinom(1,1,P)
}
thse=ability(mdl,u,b,a,c)
thr[r] = thse[1]
ser[r] = thse[2]
}
最後呈現模擬結果。我們模擬能力值為0.5的人之答題反應，故真實能力值為 0.5，得到之理論標準誤（set）為0.78。
模擬10次之能力值參數估計值為 0.50, 1.87, 1.87, 0.90, 2.30, -0.57, -1.7 , 0.32, -0.10, 0.4。計算期這10次估計之平均值為 0.58，與這10次估計值之標準差 sd(thr) 為 1.23，與其標準誤平均數 mean(ser) 為 1.05。
theta # 0.5
mean(thr)
set
sd(thr)
mean(ser)
我們可以發現模擬之平均能力估計值0.58與理論能力值 0.5 有些差距，且理論標準誤為0.78，與模擬估計值之標準標準差（1.23）與估計值之平均標準誤（1.05）也有差距。這可能是源於題數不足（J = 5）、模擬次數不足（R=10）、試題參數（b、c）與ability()設定之迭代次數與收斂標準有關係。
以下將難度參數b改為接近0.5，提高鑑別力參數。
#b = c(-0.5, -0.25, 0.0, 0.25, 0.5) #old
#a = c(1, 1.5, 0.7, 0.6, 1.8) #old
b = c(0.5, 0.45, 0.4, 0.55, 0.6)
a = c(1.4, 1.5, 1.7, 1.6, 1.8)
將模擬次數提升至2000次（R=2000），修改ability中的最大迭代次數（S=1000）與收斂標準（ccrit=.0001）。
S = 1000 # the maximum number of iterations
ccrit=.0001 #the converge criterion of Newton's method
最後計算出之平均能力估計值為0.48。理論之標準誤為0.56，能力估計值之標準差0.76，與估計值之平均標準誤為0.68。估計值較前面更接近理論值，且標準誤皆有下降。
筆者進一步將題數增加至50題，最後能力估計值之平均為0.49；理論標準誤為0.18，參數估計值之標準差為 0.18，平均參數估計值之標準誤為0.18。
b = c(0.5, 0.45, 0.4, 0.55, 0.6)
a = c(1.4, 1.5, 1.7, 1.6, 1.8)
b = rep(b,10)
a = rep(a,10)
J = length(b)
完整程式碼：
set.seed(20250913)
mdl = 2
theta = 0.5
#b = c(-0.5, -0.25, 0.0, 0.25, 0.5)
#a = c(1, 1.5, 0.7, 0.6, 1.8)
b = c(0.5, 0.45, 0.4, 0.55, 0.6)
a = c(1.4, 1.5, 1.7, 1.6, 1.8)
J = length(b)
if(mdl == 1| mdl == 2){ c = rep(0,J)}
if(mdl == 1){a = rep(1,J)}
sumdemt = 0
for(j in 1:J){
Pstar = 1/(1+exp(-a[j]*(theta - b[j])))
P = c[j] + (1 - c[j]) *Pstar
sumdemt = sumdemt - a[j]**2 *P*(1.0 - P) * (Pstar/P)**2
}
set = 1/sqrt(-sumdemt)#SE-theoretical
R=2000 # the number of replications
thr=rep(0,R) # theta hat from the replication
ser=rep(0,R) # standard error from replication
for(r in 1:R){
u=rep(0,J) #response vector u is generated from binominal distribution
for(j in 1:J){
P=c[j]+(1-c[j])/(1+exp(-a[j]*(theta-b[j])))
u[j]=rbinom(1,1,P)
}
thse=ability(mdl,u,b,a,c)
thr[r] = thse[1]
ser[r] = thse[2]
}
theta
#thr
mean(thr)
set
sd(thr)
mean(ser)
參考資料：
Baker, F. B., & Kim, S. H. (2017). The Basics of Item Response Theory Using R (Vol. 10). New York: Springer.
