前回に引き続き、PISA Data Analysis Manual SPSS 2nd Editionの分析を、Rで再現してみます。
今回は、PVs（Plausible Values: 推算値法）に関わる第8章を扱います。
使用するデータ
今回もData Analysis Manualと同じく、PISA2006のベルギー（BEL）のデータを使う。データは、Database （PISA 2006）からダウンロードできる。
以下では、ベルギーのPISA2006の結果が2006BEL.csvというファイルに保存されているという前提で話をすすめる。
PVs(1)
PISAの得点は、ラッシュモデル（Rasch Model）およびPVs（Plausible Values）という技法によって算出されている。その詳細については、Data Analysis Manualの5章・6章に詳しい。
ここでは理屈にはあまり触れず、とりあえず、RでPISAの得点を計算する方法について述べることにする。
PISAでは、個々のデータにPV1からPV5までの5つの得点が含まれている。たとえば科学リテラシーであれば、PV1SCIEからPV5SCIEまで5つの変数がある。
PISAで得点を計算する場合、これら5つの変数それぞれについて、推定値を計算した上で、さらに前回述べたReplicate Weightによる標準誤差の推定を行う必要がある。最終的な推定値は、これら5つの推定値と5つの標準誤差を合算することで求められる。
具体的な方法は以下の通り。ここでは、PISA2006のベルギーの科学リテラシーの平均点を計算することを考える。
# P.119 Table 8.2
bel2006 <- read.csv("2006BEL.csv")
# まずPV1
weighted.mean(bel2006$PV1SCIE, bel2006$W_FSTUWT)
weighted.mean(bel2006$PV1SCIE, bel2006$W_FSTR1)
weighted.mean(bel2006$PV1SCIE, bel2006$W_FSTR2)
# 以下
# 延々と・・・
# 80回
# 繰り返していく
weighted.mean(bel2006$PV1SCIE, bel2006$W_FSTR80)
# 続いてPV2
weighted.mean(bel2006$PV1SCIE, bel2006$W_FSTUWT)
# 以下
# 延々と・・・
# 80回
# 繰り返していく
# PV3からPV5についても繰り返す。
weighted.mean(bel2006$PV5SCIE, bel2006$W_FSTR80)
これもあまりに面倒である。intsvy packageを使って、もう少し簡単にしてみる。
library(intsvy)
pv1 <- pisa.mean("PV1SCIE", data=bel2006)
pv2 <- pisa.mean("PV2SCIE", data=bel2006)
pv3 <- pisa.mean("PV3SCIE", data=bel2006)
pv4 <- pisa.mean("PV4SCIE", data=bel2006)
pv5 <- pisa.mean("PV5SCIE", data=bel2006)
# 最終的な推定値は、PV1からPV5の推定値の平均になる。
m <- (pv1+pv2+pv3+pv4+pv5)[1,2] / 5
# 標準誤差の計算。まず、PV1からPV5の標準誤差の2乗の平均を計算する。
u <- (pv1[1,3]^2 + pv2[1,3]^2 + pv3[1,3]^2 + pv4[1,3]^2 + pv5[1,3]^2) / 5
# imputation varianceの計算。
t <- ((pv1[1,2]-m)^2 + (pv2[1,2]-m)^2 + (pv3[1,2]-m)^2 + (pv4[1,2]-m)^2 + (pv5[1,2]-m)^2) / 4
# 最終的な分散は、u + 1.2 * t になる。最終的な標準誤差は、sqrt(e)
e <- u + 1.2 * t
# 推定値の表示
c(m, sqrt(e))
PVsの概念については、PISA Data Analysis Manualを読んでもよいし、あるいはMultiple Imputation（多重代入法）の文献を読んでもよい。基本的な発想は同じである。
ここまでの作業は、intsvy package内にある、pisa.mean.pvを使えば簡単にできる。この関数は、おまけに標準偏差まで出力してくれる。
#P.121 Table 8.3 pisa.mean.pv(pvlabel="SCIE", data=bel2006)
男女別（ST04Q01）に成績を出力する場合は、次の通り。
#P.121 Table 8.4 pisa.mean.pv(pvlabel="SCIE", by="ST04Q01", data=bel2006)
PVs(2)
続いて回帰分析を扱う。考え方は、先ほどとほとんど同じである。
今回は、従属変数を理科の点数、独立変数をジェンダーおよびHISEI（生徒の社会経済的指標）として分析しよう。
# 変数の下準備
library(car)
bel2006$gender <- recode(bel2006$ST04Q01, "2=0")
bel2006$HISEI[bel2006$HISEI > 90] <- NA
# PV1からPV5の推定
reg1 <- pisa.reg(y="PV1SCIE", x=c("HISEI", "gender"), data=bel2006)
reg2 <- pisa.reg(y="PV2SCIE", x=c("HISEI", "gender"), data=bel2006)
reg3 <- pisa.reg(y="PV3SCIE", x=c("HISEI", "gender"), data=bel2006)
reg4 <- pisa.reg(y="PV4SCIE", x=c("HISEI", "gender"), data=bel2006)
reg5 <- pisa.reg(y="PV5SCIE", x=c("HISEI", "gender"), data=bel2006)
# 最終的な推定値はPV1からPV5の平均
m2 <- (reg1 + reg2 + reg3 + reg4 +reg5)[,1] / 5
# 標準誤差の計算。まず、PV1からPV5の標準誤差の2乗の平均を計算する。
u2 <- (reg1[,2] ^ 2 + reg2[,2] ^ 2 + reg3[,2] ^ 2 + reg4[,2] ^ 2 + reg5[,2] ^ 2) / 5
# imputation varianceの計算。
t2 <- ((reg1[,1] - m2)^2 + (reg2[,1] - m2)^2 + (reg3[,1] - m2)^2 + (reg4[,1] - m2)^2 + (reg5[,1] - m2)^2) / 4
# 最終的な分散は、u2 + 1.2 * t2 になる。最終的な標準誤差は、sqrt(e2)
e2 <- u2 + 1.2 * t2
cbind(m2, sqrt(e2))
こちらもpisa.reg.pvを使えば、簡単にできる。
#P.123 Table 8.7
pisa.reg.pv(pvlabel="SCIE", x=c("HISEI", "gender"), data=bel2006)
いずれもSPSSやSASによる推定値と一致している。
