11  推算値法

11.1 MLEやEAPの課題

前の章で見たように， IRTの能力推定法としてよく利用されるMLEやEAPには， 母集団の能力推定に利用すると推定値が偏るという欠点があります。 なぜ偏りが生じるか，簡単に説明しておきましょう。

無限に項目を用意できれば別ですが，普通はテストに出題できる項目数には限りがあります。 そのためテストの結果から推定された能力値には，必ず受検者の真の能力値（\(\theta\)）とのズレ（＝測定誤差）があります。 MLEやWLEによって推定された個人の能力値\(\hat{\theta}_i\)は， 真の能力値\(\theta_i\)に測定誤差\(\epsilon_i\)を加えたもの （\(\hat{\theta}_i = \theta_i + \epsilon_i\)）になっていると考えられます¹。 このとき，母集団における能力の分散は， 個々人の能力の分散に測定誤差の分散を加えたものとして推定されます。 つまり，\(Var(\hat{\theta}_i) = Var(\theta_i) + Var(\epsilon_i)\)なので， MLEやWLEで推定された能力値の分散は，実際の分散より大きくなってしまうのです。

EAPはどうかというと， 前の章で少し触れたように，EAPは受検者の能力分布に正規分布などの分布を仮定します。 事前に能力分布を仮定しているので分散が過大推定されることはないのですが， 今度は事前に想定した分布に推定された能力値が引き寄せられます。 一般的には，測定誤差が大きいほどEAPは母集団の分散を過小推定するそうです^[1]。

以上のような現象は測定誤差の問題なので，項目数を増やしていけば， いずれの推定法も同じような推定値に近づいていきます。 先行研究によると，だいたい180項目ほど出題すれば，MLE・WLE・EAPの差は無視できるそうです^[1]。 ただ180項目というのは相当な量で，あまり現実的ではありません。

要するに，MLEやEAP（あるいはWLE）といった能力推定法では， 母集団の特性を適切に推定することは難しいということです。 ここで登場する手法が，推算値（Plausible Values: PVs）です。

PVsは，あらかじめ受検者の能力値に正規分布などの分布を仮定します。 その意味では，EAPとよく似た手法と言えます。 PVsとEAPの違いは， EAPが唯一の受検者の能力を推定するのに対し， PVsは推定された受検者の能力分布から無作為に受検者の「あり得る能力」を抽出するという点です。 PISAやTIMSSでは，個々人の受検者に対して5個（PISAは2015以降は10個）の能力値が推定されています。

欠測値の処理に詳しい方の中には，多重代入法（Multiple Imputation: MI）^[2]という手法をご存じの方もいるでしょう。 MIは，データに含まれる欠測に伴う推定の偏りを補正するため， 欠測データの分布から無作為に複数の値を取り出し欠測値を置き換えるという手法です。 ここで重複テスト分冊法を思い出してほしいのですが， LSAでは受検者は出題された一部の項目にしか答えていません。 つまり全体としてみれば，個々の受検者の回答データは多数の欠測を含んだ不完全なものになっています。 ここで受検者の能力\(\theta\)を推定しようという試みは， 多数の欠測を含んだデータから適切な推定値を得ようとする試みと同義です。 受検者の能力分布から複数の「あり得る能力」を抽出するというPVsの考え方は， 欠測データを処理するためのMIの枠組みを教育測定に応用したものなのです。

なお，PVsが有効なのは重複テスト分冊法のように，受検者のデータが多数の欠測を含む場合に限りません。 テストの項目数が少ない場合や，テストの難易度が受検者の能力値に対して低すぎる（あるいは高すぎる） ような場合であっても，PVsを使うことで母集団の特性を適切に推定できることが明らかになっています^[3]。

11.2 Rによる確認

実際にPVsが有効かどうか，シミュレーションで確認してみましょう。 最初に仮の受検者の回答データとして，respを作成します。 母集団における能力値の分布は，実際の学力調査を意識して，以下のような複雑な状況を設定します。 まず，受検者のサンプルサイズは4000とし，うち男子2000名・女子2000名とします。 さらに能力値の平均と分散は男女で違い， 男子は平均-0.2・分散1に対し，女子は平均0.2・分散1になっています。 また能力値以外にSESという変数があると仮定し， 能力値とSESの相関係数は，男子0.5，女子0.4としています。 状況設定が複雑なので，コードは閉じています。

なお，コード実行のためにはMIRTやTAMが必要です。 インストールしていない場合は，項目反応理論の章を参考に設定を済ませてください。

set.seed(123)

# sample size
size <- 4000
# Item parameters
item_n <- 20
b <- seq(-3, 3, length.out = item_n)
a <- 0.9
# Population parameters
sig1 <- matrix(c(1.0, 0.4, 0.4, 1.0), nrow = 2)
sig2 <- matrix(c(1.0, 0.5, 0.5, 1.0), nrow = 2)
sex <- c(rep(1, size / 2), rep(-1, size / 2))
sim1 <- mvtnorm::rmvnorm(n = size / 2, mean = c(0.2, 0), sigma = sig1)
sim2 <- mvtnorm::rmvnorm(n = size / 2, mean = c(-.2, 0), sigma = sig2)
sim <- rbind(sim1, sim2)
theta <- sim[, 1]
ses <- sim[, 2]

# resp pattern
resp <- mirt::simdata(
  a = rep(1.7 * a, item_n),
  d = -b * 1.7 * a,
  Theta = theta,
  itemtype = "2PL"
)

それでは生成したrespに対して，TAMで能力値を推定してみましょう。 Rのコードは，次のようになります。 なお，通常はPVsは複数の値を生成するのですが， 処理が重くなるのでここでは1つの値だけを生成（nplausible = 1）しています。

# estimation
mod <- TAM::tam.mml.2pl(resp, verbose = FALSE)

# ability estimation
mle <- TAM::tam.wle(mod, progress = FALSE)$theta
eap <- mod$person$EAP
pv <- TAM::tam.pv(mod, nplausible = 1, verbose = FALSE)$pv$PV1.Dim1

# 能力推定
ability <- list(
  "theta" = theta,
  "MLE" = mle,
  "EAP" = eap,
  "PV" = pv
)

続いて，TAMを使って推定した能力値が， 適切に母集団の値を復元しているかどうか検討してみましょう。 最初に，母集団の平均と分散を見てみます。 なお以下の出力では， thetaが真の値，MLE・EAP・PVはそれぞれの推定法による推定値を示します。

# 平均値
round(sapply(ability, mean), 2)

theta   MLE   EAP    PV 
 0.01  0.00  0.00  0.00 

# 分散
round(sapply(ability, var), 2)

theta   MLE   EAP    PV 
 1.01  1.21  0.82  1.00 

先の章と同じく，MLEやEAPの推定値は，平均は適切に推定していますが， 分散を過大推定しています。 これに対してPVの推定値は，やや過小推定（0.96）ではあるものの， MLEやEAPに比べればthetaに近い分散を推定しています。

続いて，男女別の平均と分散，そして能力値とSESの相関を検討してみましょう。

# 男女別のスコア，SESとの相関を計算する関数
calc_stats <- function(ability, stat) {
  if (identical(stat, cor)) {
    sapply(ability, function(x) {
      cor(x, ses)
    })
  } else {
    sapply(ability, function(x) {
      tapply(x, sex, stat)
    })
  }
}

round(calc_stats(ability, mean), 2)

   theta   MLE   EAP    PV
-1 -0.20 -0.21 -0.17 -0.16
1   0.21  0.21  0.17  0.16

round(calc_stats(ability, var), 2)

   theta  MLE  EAP   PV
-1  0.97 1.16 0.78 0.96
1   0.96 1.19 0.80 1.00

round(calc_stats(ability, cor), 2)

theta   MLE   EAP    PV 
 0.42  0.38  0.38  0.36 

男女別の能力値の平均を見てみると，MLEはthetaとほぼ同じ（男子-0.21・女子0.21）値になっています。 他方，EAPとPVはいずれもthetaと値がズレています。

分散の推定値については，母集団の推定と同じく， MLEは過大推定，EAPは過小推定の傾向が見られます。 PVsは0.93とやや過小推定ですが，それでも比較的thetaに近い値が得られています。

相関係数については，thetaの値が0.42に対し， MLEとEAPは0.38でやや過小推定になっています。 PVsは0.34で更に過小推定であり，母集団の特性を適切に推定しているとは言えません。

11.3 推算値の条件付

先ほどのシミュレーションから言えることは， PVsは母集団の特性を適切に推定するわけではないということです。 分散はわずかに過小推定するにとどまっていますが， 男女の平均や相関係数の偏りなど，MLEやEAPと比べてPVsが優れているとは言えません。

それではなぜLSAではPVsが有効だと言われるのでしょうか。 実は先ほどのシミュレーションで推定したPVsが偏った推定値を返したのは， 母集団の能力値の分布について誤った仮定をおいていたからです。 すでに述べたように，EAPやPVsでは事前に母集団の能力分布について， 正規分布などの仮定を起きます。 しかし先のシミュレーションのように，男女で能力値の平均が異なるといった状況はありえます。 この場合，母集団の能力分布は，異なる2つの正規分布が重なった分布になっており， 単純な正規分布とは異なる形状になっています。

LSAでPVsを利用する場合，性別やSESのような母集団の能力分布に影響を与えるであろう要因を 事前に考慮して能力分布を設定し，項目パラメータや能力値の推定を行わなければなりません。 こうした母集団の要因を考慮することを，一般に「条件づけ（Conditioning）」と呼びます。 なお，LSAのようにサンプルサイズが大きいデータに対して，条件付けを行った推定を行うと， 長い時間がかかったり，推定がうまくいかない場合があります。 そこで一般的なLSAでは，いったん条件付けを行わずに項目パラメータを推定し， そこで得られた推定値を利用しつつ条件付けを行った能力推定を行うという 2段階の推定法が採用されています^[4]。 2段階の推定は，RのTAMパッケージに実装されており，tam.latreg関数を使うことで実行可能です。

以下では，条件付けを行って推定したPVs（PV2とします）を使い， 母集団の推定を行ってみましょう。 最初に，母集団の平均と分散を確認します。

likeli <- CDM::IRT.likelihood(mod)
mod2 <- TAM::tam.latreg(likeli, Y = cbind(sex, ses), verbose = FALSE)
# 後で使うのでPVを5つ生成し，そのうち1つだけを利用
pvs <- TAM::tam.pv(mod2, nplausible = 5, verbose = FALSE)
pv2 <- pvs$pv$PV1.Dim1

ability2 <- c(ability, list("PV2" = pv2))

# 平均値
round(sapply(ability2, mean), 2)

theta   MLE   EAP    PV   PV2 
 0.01  0.00  0.00  0.00 -0.01 

# 分散
round(sapply(ability2, var), 2)

theta   MLE   EAP    PV   PV2 
 1.01  1.21  0.82  1.00  0.99 

母集団の平均・分散ともに，PV2がthetaに非常に近い値を返していることがわかります。 同様に，男女別の平均と分散，SESとの相関係数も検討してみましょう。

round(calc_stats(ability2, mean), 2)

   theta   MLE   EAP    PV   PV2
-1 -0.20 -0.21 -0.17 -0.16 -0.21
1   0.21  0.21  0.17  0.16  0.20

round(calc_stats(ability2, var), 2)

   theta  MLE  EAP   PV  PV2
-1  0.97 1.16 0.78 0.96 0.94
1   0.96 1.19 0.80 1.00 0.96

round(calc_stats(ability2, cor), 2)

theta   MLE   EAP    PV   PV2 
 0.42  0.38  0.38  0.36  0.42 

いずれの値も，PV2がもっともthetaと近い値を返しています。

11.4 推算値の使い方

ここまでの分析では，簡単のためにPVsは1つだけを生成していました。 シミュレーションの結果からもわかるように，条件づけたPVsは1つの値しか生成しなくても 母集団の推定には有効です。 ただ，より正確な推定のためには複数の値を生成し，それぞれで得られた推定値を合成する必要があります。

PVsを使った推定値を合成する際は，Rubinのルールと呼ばれる方法を使います。 今，PVsの数を\(M\)，\(i\)番目のPVを使った推定値を\(\theta_i\)，\(i\)番目のPVを使った推定値の分散を\(V_i\)とすると， PVsを使った推定値（\(\theta\)），及びその分散（\(V\)）は，以下の式になります。

\[\theta = \frac{1}{M}\sum_{i=1}^M{\theta_i}\]

\[V = U + (1 + \frac{1}{M})B_m\] \[U = \frac{1}{M}\sum_{i=1}^M{V_i} \qquad B_m = \frac{1}{M-1}\sum_{i=1}^M(\theta_i - \theta) ^ 2\]

これでは何のことかわからないと思いますが， まず最終的な推定値（\(\theta\)）は，単純に個々のPVで得られた推定値（\(\theta_i\)）の平均です。 推定値の分散（つまり標準誤差の2乗）は，個々のPVで得られた推定値の分散（\(V_i\)）の平均（これが\(U\)です）に， \((1+\frac{1}{M})B_m\)を足したものになります。 \(B_m\)は個々のPVごとの推定値のばらつきを示し， PVが受検者の能力分布から得られたランダムな値であることによる推定値のブレを考慮しています。

PISA2012を例に，Rubinのルールを使った推定の例を見ておきましょう。 ここでは読解リテラシー（PV1READからPV5READ）の平均値を考えます。

url <- "https://raw.githubusercontent.com/kawa5902/LSAdata/refs/heads/main/pisa2012stuJPN.csv"
jpn2012 <- read.csv(url)

pv1 <- intsvy::pisa.mean("PV1READ", data = jpn2012)
pv2 <- intsvy::pisa.mean("PV2READ", data = jpn2012)
pv3 <- intsvy::pisa.mean("PV3READ", data = jpn2012)
pv4 <- intsvy::pisa.mean("PV4READ", data = jpn2012)
pv5 <- intsvy::pisa.mean("PV5READ", data = jpn2012)

est <- c(pv1$Mean, pv2$Mean, pv3$Mean, pv4$Mean, pv5$Mean)
# s.e.なので注意。s.eは分散の標準誤差
se <- c(pv1$s.e., pv2$s.e., pv3$s.e., pv4$s.e., pv5$s.e.)

round(est, 3) # 5つの推定値

[1] 538.116 537.872 538.473 538.064 537.732

round(se^2, 3) # 推定値の分散

[1] 13.441 13.372 13.082 13.344 13.742

5つのPVの推定値は，estおよびseに格納されています。 これを使って最終的な推定値を計算します。 なお，結果の解釈に必要なのは推定値の分散ではなく標準誤差の方ですので， 最後にvの平方根を取って標準誤差を算出しています。

f_est <- mean(est)
u <- mean(se^2)
b <- var(se^2)
v <- u + (1 + 1 / 5) * b
round(c(mean = f_est, var = v), 2)

  mean    var 
538.05  13.46 

round(c(mean = f_est, se = sqrt(v)), 2)

  mean     se 
538.05   3.67 

得られた推定値（平均538.05，標準誤差3.67）が， PISAの報告書にも記載されている適切な推定値です。 ここでは平均値を例に紹介しましたが，分散や回帰分析の係数など， 他の推定値であっても同じようにRubinのルールを適用して推定値を得ることができます ²。

さて，こうして得られた最終的な推定値（日本のPISA2012の読解リテラシーの平均値）を PV1READから得られた値と比べてみると，実はほとんど違いがありません。 PVを5回も計算するのが面倒だという場合は， とりあえずPV1のみを使って分析を行っても，深刻な解釈ミスにつながることはないと考えられます³。

pv1 # 最終的な推定値とあまり変わらない

  Freq     Mean     s.e.      SD      s.e
1 6351 538.1164 3.666172 99.1964 2.294498

5回もPVを計算するのは面倒だという人も多いでしょうが， 実際は，intsvyを使うことでPVの処理も自動化できます。 intsvyを使った自動化については，後の章で扱います。

intsvy::pisa.mean.pv(paste0("PV", 1:5, "READ"), data = jpn2012)

  Freq   Mean s.e.    SD  s.e
1 6351 538.05 3.67 98.69 2.27

なお，PVsを使うときに複数のPVを個人レベルで平均してはいけません。 PVsを個人レベルで平均すると，その値はEAPに近いものになります。 以下では，生成した5つのPVsをを個人レベルで平均したもの（PV3）を使った例を示します。 推定結果を見ると，PV3はEAPと同じように母集団の分散を過小推定していることがわかります⁴。

# 5つのPVsをlistに変換
pvs_list <- TAM::tampv2datalist(pvs)
# 5つのPVsを個人レベルで平均
pv3 <- Reduce("+", lapply(1:5, function(i) pvs_list[[i]]$PV.Dim1)) / 5
ability3 <- c(ability2, list("PV3" = pv3))

# 平均値
round(sapply(ability3, mean), 2)

theta   MLE   EAP    PV   PV2   PV3 
 0.01  0.00  0.00  0.00 -0.01 -0.01 

# 分散
round(sapply(ability3, var), 2)

theta   MLE   EAP    PV   PV2   PV3 
 1.01  1.21  0.82  1.00  0.99  0.87 

# 男女別の平均値
round(calc_stats(ability3, mean), 2)

   theta   MLE   EAP    PV   PV2   PV3
-1 -0.20 -0.21 -0.17 -0.16 -0.21 -0.21
1   0.21  0.21  0.17  0.16  0.20  0.20

# 男女別の分散
round(calc_stats(ability3, var), 2)

   theta  MLE  EAP   PV  PV2  PV3
-1  0.97 1.16 0.78 0.96 0.94 0.83
1   0.96 1.19 0.80 1.00 0.96 0.82

# SESとの相関係数
round(calc_stats(ability3, cor), 2)

theta   MLE   EAP    PV   PV2   PV3 
 0.42  0.38  0.38  0.36  0.42  0.45 

1. 袰岩晶・篠原真子・篠原康正. (2019年). PISA調査の解剖ー能力評価・調査のモデル. 東信堂.

2. 高橋将宜・渡辺美智子. (2017年). 欠測データ処理. 共立出版.

3. Wu, M. (2005年). The role of plausible values in large-scale surveys. Studies in educational Evaluation, 31(2-3), 114–128.

4. Zheng, X. (2024年). On generating plausible values for multilevel modelling with large-scale-assessment data. British Journal of Mathematical and Statistical Psychology, 77(1), 212–236.

---

1. 推定された能力値なので\(\theta\)ではなく，\(\hat{\theta}\)としています。

2. Rubinのルールには，推定値の分布が正規分布であるという前提があるため， たとえば相関係数には適用できません^[2]。 ただLSAの場合はサンプルサイズが大きいこともあって，相関係数（を含むその他の分析）についても Rubinのルールを適用することが一般的です。

3. LSAの分析はReplication MethodとPVsを使うので，それなりに時間がかかります。 たとえばPISAの場合，Replication Methodだけで80回の計算が必要です。 個々のPVについて80回の計算を行うとなると，PVが5つの場合は全部で400回の計算になります。 PISA2015以降はPVが10個になっているので，最終的な推定値を得るまでに800回以上の計算が必要となるのです。 PVを1つだけ使えば，計算回数を大きく抑えることができます。

4. 通常，LSAのデータセットに含まれるPVsは条件付けを行った能力分布から抽出されていますので， その個人レベルの平均は，通常のEAPではなく条件づけたEAPになります。 これが，シミュレーションでEAP（＝条件づけていないEAP）とPV3（≒条件づけたEAP）の値が微妙に違う理由です。