R言語による項目反応理論

Rで項目反応理論（2PLM：2パラメータ・ロジスティックモデル）のパラメータを推定するスクリプトについて、これまで参考にさせていただいていたHPが削除されたようなので、自分用のメモを兼ねてまとめておきます。

lbrary(ltm)
result=ltm(resp~z1) #2パラメータモデル 

とりあえず、ライブラリltmを呼び出してパラメータの推定を行います。
上のrespは0/1の回答データのオブジェクトです。以下の記事でIRTモデルに従うランダムな回答データを生成するコードを紹介しています。

項目パラメータの推定結果はresultというオブジェクトに格納しています。直接resultの中身を表示すると以下のようになります。

> result
Call:
ltm(formula = resp ~ z1)
Coefficients:
         Dffclt  Dscrmn
Item 1    0.783   0.843
Item 2    0.560   1.312
Item 3    1.451   0.740
Item 4   -0.313   0.289
Item 5   -1.269   0.321
Item 6   -1.069   1.543
Item 7   -0.813   1.032
Item 8   -0.879   0.760
Item 9   -0.346   1.412
Item 10   0.801   0.254

Dffcltがbパラメータとか困難度パラメータと呼ばれる値、Dscrmnがaパラメータ、識別力パラメータと呼ばれる値です。

resultの中には、他にもいくつかのデータが入っています。

> names(result)
 [1] "coefficients" "log.Lik"      "convergence"  "hessian"      "counts"      
 [6] "patterns"     "GH"           "max.sc"       "ltst"         "X"           
[11] "control"      "IRT.param"    "formula"      "call" 

このうち、項目パラメータはcoefficientsに格納されています。
（IRT.paramというデータもありますが、中身はnullのようです。もしかするとsummaryやplotで中身を見られるのかもしれません。）
項目パラメータの推定結果を利用する場合は、coefficientsの中身を利用すればいいのですが、一つ問題があります。

> result$coefficients
         (Intercept)          z1
Item 1  -0.660092129  0.84301049
Item 2  -0.735092413  1.31187704
Item 3  -1.073375578  0.73980573
Item 4   0.090410302  0.28882351
Item 5   0.407555088  0.32120565
Item 6   1.648851420  1.54310697
Item 7   0.838497257  1.03178179
Item 8   0.667796344  0.75984712
Item 9   0.488267142  1.41187472
Item 10 -0.203312944  0.25376898

coefficientsの中身を見ると、何かオプションで制御できる可能性がありますが、元の推定結果のオブジェクト（result）ではDffcltとDscrmnのパラメータであったのが、result$coefficientsでは（intercept）とz1という形に表示が変わっています。理由は不明ですがltmでは推定結果を直接表示するのと、coefficientsを取り出す場合で想定するモデルの形式が変わっているようです。

※result=ltm(resp~z1,IRT.param=FALSE)のようにオプションでIRT.param=FALSEを指定すると、（intercept）とz1のほうに結果が統一されます。デフォルトはIRT.param=TRUEです。

coefficientsからDffcltとDscrmnを計算するには、例えば次のようにします。

temp=result$coefficients
dscrmn=temp[,2]
dffclt=-temp[,1]/dscrmn

結果を比較すると

> cbind(dffclt,dscrmn)
              dffclt      dscrmn
Item 1   0.783017691  0.84301049
Item 2   0.560336367  1.31187704
Item 3   1.450888425  0.73980573
Item 4  -0.313029580  0.28882351
Item 5  -1.268829152  0.32120565
Item 6  -1.068526972  1.54310697
Item 7  -0.812669177  1.03178179
Item 8  -0.878856186  0.75984712
Item 9  -0.345828943  1.41187472
Item 10  0.801173365  0.25376898

 #resultの結果 
Coefficients:
         Dffclt  Dscrmn
Item 1    0.783   0.843
Item 2    0.560   1.312
Item 3    1.451   0.740
Item 4   -0.313   0.289
Item 5   -1.269   0.321
Item 6   -1.069   1.543
Item 7   -0.813   1.032
Item 8   -0.879   0.760
Item 9   -0.346   1.412
Item 10   0.801   0.254

上記のように一致していることがわかります。教科書によってはD=1.7としてモデルに定数を入れているものがありますが、ltmの推定結果にはこの項は含まれておらずD=1になっています。D=1.7にしたい場合にはaパラメータ（dscrmn）を1.7で割っておく必要があります。

なお、個人ごとの能力特性値の推定は

factor.scores(result,resp.patterns=resp)$score.dat$z1

のようにします。

> mean(score)
[1] -0.005572732
> sd(score)
[1] 0.9233741

スコアの平均は、概ね0になるようです。

この記事が参考になりましたら幸いです。