Stanで多次元項目反応理論

 

この記事は、Stan advent calendar 2016の19日目の記事です。

今回は、多次元項目反応理論をStanでやってみたいと思います。

ちょっと今回は時間取れなくてやっつけです。すみません。

 

項目反応理論とは

項目反応理論(IRT)は、テストを評価するための数理モデルです。複数の項目について正解・不正解がデータで与えられたとき、そこから回答者の潜在的な学力などを推定する、そして問題の特徴も推定する方法です。

IRTの確率モデルは、以下のようになります。

mirt_stan2

ここで、yはある回答者i、項目jに対する回答で、正解なら1、不正解なら0が入力されたデータです。これが正答率πをパラメータに持つベルヌーイ分布から生成されたと考えます。

そして、正答率πは、ロジスティックモデルか正規累積モデルで表現されます。上では、ロジスティックモデルで書いてますが、1.7という係数は、正規累積モデルに近似するための係数です。IRTでは計算が簡単なこちらがよくつかわれます。

さて、πは項目についてのパラメータα、β、そして回答者についてのパラメータθによって構造化されます。αは識別力といい、項目がテスト全体に与える寄与で、βは閾値パラメータと呼ばれます。なお、普通のIRTでは、‐1.7α(θ-β)と構造化し、その場合のβは困難度パラメータと呼ばれますが、今回はStanに実装しやすいように(あと多次元IRTに対応するように)このようにしています。

最後に、θは標準正規分布に従うと仮定します。

これを多次元に拡張するには、いろいろごにょごにょあるんですが、ここは省略して、とりあえず潜在変数θが多変量正規分布から生成されると仮定するだけでOKとしておきます。ただし、探索的因子分析(山口大学小杉先生のStanアドカレ記事)と同様、多次元IRTでも識別力に制約を与えないと識別不定になります。よって、ここでも共分散行列のコレスキー分解した要素を使って、制約を与えます。

ただし、IRTの場合、データの生成メカニズムに正規分布を使っていないため、潜在変数を積分消去して多変量正規分布を確率モデルにする、という手が使えません。なのでここでは潜在変数も推定していきます。

 

Stanでの実装

今回は多次元IRTということで、θが多次元の場合のモデルを考えます。よって、確率モデルは複雑になりますが、基本は変わりません。あと、2値だけではなく、多値の順序反応尺度でも使えるように拡張してみました。順序尺度に対応した確率モデルには、Stanには順序ロジスティック分布というものが用意されています。詳細はマニュアルを見てください。これを使うととってもらくちんです。

Stanコードは以下です。

 

注意点は、θは多変量正規分布から生成していますが、平均が0ベクトル、共分散行列が単位行列にしてあるので、結局は全部が独立した正規分布から生成しているのと同じです。あと順序ロジスティック分布は、リンク関数を確率モデルに含んでいるので、ロジスティックリンクは省略できます。

Rコードは以下です。

 

Stanコードはcategorical_fa2.stanというものを読み込んでいます。

また、多次元IRTはおそろしく時間がかかるので(潜在変数も同時に推定しているため)、MCMCではなくて変分ベイズで推定しています。変分ベイズは95%信頼区間などはかなり怪しいですが、点推定値については、回帰モデルでない限りはいい感じに推定してくれる(という僕の主観)と思うので、こういう因子分析モデルでは使えそうです。

ただ、vbのデフォルト設定はかなり甘いので、ちょっと設定をいじってます。

 

識別力の回転

多次元IRTの場合、上に書いたように因子負荷量に制約が必要です。これだと単純構造解にならず、いまいち解釈がしづらいです。そこで、ここでは簡易的に事後分布の平均をとってきて、それをGPArotationの回転関数に入れて回転させる方法を書いてみます。

次のRコードで簡単にできます。

 識別力を回転させて、その回転行列をθにかけてやれば、θも回転できます。

 

ちゃんと推定できてるかチェック

今回はテストように多次元IRTのモデルから乱数を生成して、真値に一致するかをチェックしてみましょう。ただ、Rには順序カテゴリカルの乱数を出してくれるものがないので、Stanを使って乱数を出してみます。

乱数生成用のStanコードは以下です。

 

ここにαとβ、θをいれてやれば、それに合ったデータyが生成されます。

Rコードは以下です。関数の形にしてみました。

 これでサンプルサイズ500、25項目のデータが生成できます。

上のStanとRのコードで分析した結果は以下です。

 識別力は該当する因子に対して、それぞれ1を真値としました。それなりにうまく推定できているようです。

続いて、閾値。-2.5、-0.8、0.8、2.5が真値です。

 

これもうまく推定できているようです。

最後に、θです。これは元のθとの相関を出してみました。対応する因子はズレてますが、0.9近い相関が得られています。Nを増やせばもっと相関は上がると思います。

 

雑な記事になってすみませんが、以上です。

 

 

This entry was posted in 心理統計学. Bookmark the permalink.