日々の記録簿

日々のできごとの記録

PISA

RによるPISA Data Analysis Manualの再現(第8章: PVsを使う)

投稿日:

前回に引き続き、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による推定値と一致している。

スポンサーリンク

スポンサーリンク

-PISA

執筆者:


comment

メールアドレスが公開されることはありません。 が付いている欄は必須項目です

関連記事

PISA2015の分析方法とintsvyの使い方

Rのintsvy packageを使ってPISA2015を分析する方法です。 データの準備 まず、OECDのホームページからPISA2015のデータをダウンロードする。PISA2015からSPSSのデ …

SPSSとSASでsample weightsの扱いが違う

PISA Data Analysis Manual を読んでいたら興味深い記述がありました。 However, even with an equi-probabilistic sample, stat …

PISAの分析方法(お知らせ)

PISAデータの分析方法に関する資料を、githubにまとめることにしました。ぼちぼち追記していきますので、興味のある方はどうぞ。

RによるPISAの分析(PISAデータをCSVに変換する)

PISAデータをcsvに変換する作業が面倒になってきたので、以前紹介していたpythonスクリプトを書き直し、全自動で行えるようにしました。 ここからpisacsv.zipをダウンロードし、展開してく …

RによるPISAの分析(データの準備)

PISAをRで分析する際に、最初の壁になるのが、データファイルが大きすぎるという点でしょう。 すでに述べたように、PISAのデータは、(1)データが固定長のtxtファイルになっている、(2)参加国すべ …