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になるようです。
この記事が参考になりましたら幸いです。