infer.netの例題シリーズ6 項目反応理論(IRT)＋α

今回は「Difficulty versus ability」をBUGSで実装します。この例題のもとはDARE modelというもので以下の論文になります。 "How To Grade a Test Without Knowing the Answers --- A Bayesian Graphical Model for Adaptive Crowdsourcing and Aptitude Testing" by Bachrach et al (ICML 2012)

内容は項目応答理論（IRT; Item Response Theory）の亜種です。項目応答理論については、ここやwikipediaの記事を読んでください。ざっくりいいますと、たとえば40人が100問ある4択問題を回答したデータがあって、そのデータから人の能力や問題の難しさを推定するものです。今回は正解が与えられていないとして、正解も予測します。

今回の実装もほぼ項目応答理論の基本部分を踏襲しています。すなわち、以下の部分です。

• 人による能力の違い（ability[i]）
• 問題による難しさの違い（difficulty[j]）
• 上記の論文では ability[i] - difficulty[j] を advantage と呼んでいます。
• 問題による識別の良さ（discrimination[j];　値が高いとスパッと正答者群と誤答者群が分かれる）
• advantageにdiscriminationをかけてものが正答するかどうかの確率変数のlogit（もしくはprobit）になる。

DARE modelが項目応答理論と若干異なるところは、正答できなかった時に回答をランダム（=すべてが等確率のカテゴリカル分布; ≒サイコロ）で決めるところです。空欄（NA）もたまにあります。しかも、回答の分布から正解を予測します。ただし、能力が高い人たちの多数派が正答するという前提です。そのため、優秀な人ほどひっかかる問題とか、ほとんどの人が同じひっかけ選択肢を選んでしまうような場合には使えません。

DARE modelのグラフィカルモデルは以下になります。上記論文のFigure.1からの引用です。

Infer.netの方はデータが陽に与えられていないため、今回は自分でデータを作成して解きました（記事の最後にscriptを載せました、自由に使ってください。）。40人、100問、4択問題のデータです。経験からdiscriminationは対数正規分布で作ってみました。データの中身は以下のようになっています。行が人に対応し、列が問題に対応します。1-4の数字は回答です。まれに空欄（NA）を入れてあります。

Item.1	Item.2	…	Item.99	Item.100
2	4	…	2	3
2	4	…	4	3
3	4	…	2	3
…	…	…	…	…
4	3	…	4	3
2	2	…	2	3
1	1	…	2	3

BUGSコードは以下になりました。

model {
   for (i in 1:N) {
      ability[i] ~ dnorm(0, 1)
      for (j in 1:N.question) {
         Y[i,j] ~ dcat(q[i,j,])
         for (c in 1:N.choice) {
            q[i,j,c] <- correct[i,j] * equals(correct.answer[j], c) + (1-correct[i,j]) * Theta[c]
         }
         correct[i,j] ~ dbern(p[i,j])
         p[i,j] <- phi(discrimination[j] * (ability[i] - difficulty[j]))
      }
   }
   for (j in 1:N.question){
      correct.answer[j] ~ dcat(Theta[])
      difficulty[j] ~ dnorm(mu.difficulty, tau[1])
      log(discrimination[j]) <- log.discrimination[j]
      log.discrimination[j] ~ dnorm(mu.log.discrimination, tau[2])
   }
   mu.difficulty ~ dnorm(0, 1.0E-4)
   mu.log.discrimination ~ dnorm(0, 1.0E-4)
   for (k in 1:N.tau){
      tau[k] <- pow(s[k], -2)
      s[k] ~ dunif(0, 1000)
   }
}

• 9行目: 正答する（correct=1）かどうかが決まります。
• 7行目: 正答・誤答を両方加味してカテゴリカル分布のパラメータを決めています。
• 5行目: このパラメータにしたがって実際選んだ選択肢（Y[i,j]）が決まります。

このBUGSコードはカテゴリカル分布のパラメータをいじっているので、WinBUGSで実行すると簡単にTrapが出ます。修正も容易ではありません。このような場合はJAGSの出番です。JAGSはカテゴリカル分布のパラメータの和が1になるように自動で補正します。

キックするRコードは以下になります。

library(rjags)

d <- read.delim("input/data.txt", as.is=T, header=T, sep="\t")
N <- nrow(d)
N.question <- ncol(d)
N.choice <- 4
N.tau <- 2

major.answer <- as.vector(apply(d, 2, function(x){ which.max(table(x)) }))

data <- list(
   N = N,
   N.question = N.question,
   N.choice = N.choice,
   N.tau = N.tau,
   Theta = rep(1/4, 4),
   Y = as.matrix(d)
)
inits <- list(
   correct.answer = major.answer,
   ability = rep(0, N),
   difficulty = rep(0, N.question),
   log.discrimination = rep(0, N.question),
   mu.difficulty = 0,
   mu.log.discrimination = 0,
   s = rep(1, N.tau)
)

m <- jags.model("model/model.bugs", data, inits, n.chains=3, n.adapt=1000)
update(m, 3000)
post.list <- coda.samples(m, c("correct.answer", "ability","difficulty","log.discrimination","mu.difficulty", "mu.log.discrimination", "s"), n.iter=100000, thin=100)
save.image("output/result.RData")

• 9行目, 20行目: 多数派の解答を正解として初期値を与えています。

---

結果は以下になります。

Fig. 上から順に ability, difficulty, log(discrimination)の真の値と推定値の比較 黒点: MCMCサンプルの中央値 黒実線: MCMCサンプルの95%信頼区間 カラー点: 真の値

difficulty, log(discrimination)の推定において、高い or 低い値がズレていますね。これらの標準偏差であるs[1],s[2]が低く見積もられているのが原因のようです。ぐぬぬ。どうしたものか。いつかStanでやってみるか。

parameter	true	median	mean	sd	2.5%	25%	50%	75%	97.5%	Rhat	n.eff
mu.difficulty	-0.3	-0.70	-0.71	0.45	-1.66	-0.99	-0.70	-0.41	0.13	1.00	3000
mu.log.discrimination	-1.4	-1.17	-1.18	0.21	-1.60	-1.32	-1.17	-1.04	-0.78	1.00	890
s[1]	4.5	3.32	3.40	0.72	2.26	2.89	3.32	3.80	5.01	1.00	1000
s[2]	1.3	0.92	0.94	0.23	0.53	0.78	0.92	1.08	1.44	1.00	3000

最後に推定された正解（の最頻値）と真の正解を混合表で比べました。最頻値を得るRコードはこんな感じです。

df <- data.frame(post.mcmc[,sprintf("correct.answer[%d]", 1:N.question)])
answer.mode <- apply(df, 2, function(x){ names(which.max(table(x))) })

	1	2	3	4
1	25	0	1	0
2	1	26	0	0
3	0	0	18	2
4	0	1	1	25

左が推定された正解, 上が真の正解です。こちらは思ったより良好です。

---

【参考】 データを生成したscript

library(Rlab)

set.seed(123)
N <- 40
N.question <- 100
N.choice <- 4
answer <- sample(N.choice, N.question, replace=T)
ability <- rnorm(N, 0, 1)
difficulty <- rnorm(N.question, -0.3, 4.5)
discrimination <- exp(rnorm(N.question, -1.4, 1.3))

ability.m <- matrix(rep(ability, N.question), N, N.question)
difficulty.m <- matrix(rep(difficulty, each=N), N, N.question)
discrimination.m <- matrix(rep(discrimination, each=N), N, N.question)

correct.p <- pnorm(discrimination.m * (ability.m - difficulty.m), 0, 1)
correct <- matrix(rbern(length(correct.p), correct.p), N, N.question)

data <- matrix(0, N, N.question)
index <- which(correct == 1, arr.ind=T)
data[index] <- answer[index[,2]]
data[correct == 0] <- sample(N.choice, sum(correct == 0), replace=T)
colnames(data) <- paste("Item.", 1:N.question, sep='')
rownames(data) <- paste("Person.", 1:N, sep='')
write.table(data, file="input/data.txt", quote=F, sep="\t", row.name=F)
write.table(answer, file="input/correct_answer.txt", quote=F, sep="\t", row.name=F, col.name=F)

• 17行目: rbern()を使うために{Rlab}パッケージを使っていますが、sample関数で十分と気づきました。