学力調査の基本を勉強する上で大変参考になる、PISA Data Analysis Manual SPSS 2nd Editionの分析を、Rで再現してみます。
今回は、標準誤差の推定に関わる第7章を扱います。
使用するデータ
PISA Data Analysis Manualと同じく、PISA2003のDEU（ドイツ）の結果を使う。データは、Database （PISA 2003）からダウンロードできる。
ただし、PISAのデータをRで扱うのは少々面倒である。なぜなら、(1)データが固定長のtxtファイルになっている、(2)参加国すべてのデータが1つのファイルにまとまっておりファイルサイズが大きい、という問題があるからだ。read.fwf関数で固定長のtxtデータを取り込むことも可能だが、PISAはファイルサイズが大きすぎる。
PISAデータをRで扱う場合は、データをRで扱いやすい形に変換する必要がある。
もっともわかりやすいのは、PISAの分析ソフトとして推奨されているSPSSかSASでデータを取り込み、そこから国別のcsvファイルに変換してしまうことである。
とりあえず以下では、ドイツのPISA2003の結果が2003DEU.csvというファイルに保存されているという前提で話をすすめる。
標準誤差の推定(1)
104ページから始まる第7章の分析結果をRで再現する。
PISAで推定値および標準誤差を計算する場合は、PISAで採用されているReplicate法（PISAではBRRを使用している）を考慮する必要がある。
具体的な平均値の推定を考えよう。ここでは、ドイツ（DEU）のHISEI（家族の社会経済指標: family socio-economic index）を例にとる。
はじめに、Final Weight（W_FSTUWT）および80のReplicate Weight（W_FSTR1から80まで）を使用した平均値を計算する。
# P.103 Box 7.1
deu2003 <- read.csv("2003DEU.csv")
deu2003$HISEI[deu2003$HISEI==99] <- NA
weighted.mean(deu2003$HISEI, deu2003$W_FSTUWT, na.rm=T)
weighted.mean(deu2003$HISEI, deu2003$W_FSTR1, na.rm=T)
weighted.mean(deu2003$HISEI, deu2003$W_FSTR2, na.rm=T)
# 以下
# 延々と・・・
# 80回
# 繰り返していく
weighted.mean(deu2003$HISEI, deu2003$W_FSTR80, na.rm=T)
PISAでは、こうして計算した81の値を使って、最終的な推定値を求めていく。
とは言え、これはあまりに面倒である。せっかくRを使っているのだから、もう少しRっぽく作業してみよう。
# P.105 Table 7.1
# Final Weightによる平均値。この値が最終的な推定値になる。
f_weight <- weighted.mean(deu2003$HISEI, deu2003$W_FSTUWT, na.rm=T)
# Replicate Weightによる平均値を処理する関数
deu03_hisei <- function(x){
weighted.mean(deu2003$HISEI, x, na.rm=T)
}
# 80回のReplicate Weightによる平均値。こちらは標準誤差を計算するために使う。
apply(deu2003[320:399], 2, deu03_hisei)
# P.106 Table 7.2
# 標準誤差の推定。まず、Final WeightとReplicate Weightの差の2乗を計算する。round関数で丸めておく。
round((apply(deu2003[320:399], 2, deu03_hisei) - f_weight) ^ 2, 4)
# さらに、その合計値を計算する。
sum((apply(deu2003[320:399], 2, deu03_hisei) - f_weight) ^ 2)
# これを20で割った数値が標準誤差の2乗になる。
# なぜ20かと言えば、BRRのFay係数が関係している。この点については、第4章を参照。
sum((apply(deu2003[320:399], 2, deu03_hisei) - f_weight) ^ 2) / 20
# 標準誤差を計算する
sam_var <- sum((apply(deu2003[320:399], 2, deu03_hisei) - f_weight) ^ 2) / 20
se <- sqrt(sam_var)
# 最終的な推定値は、49.33(0.42）となる。
c(f_weight, se)
実際は、intsvy packageを使うことで、こうした計算を意識することなく分析できる。
library(intsvy) pisa.mean(variable="HISEI", data=deu2003)
なお、intsvyを使えば、男女別（ST03Q01）にHISEIの平均値を計算することもできる。
deu2003$ST03Q01[deu2003$ST03Q01==8] <- NA pisa.mean(variable="HISEI", by="ST03Q01", data=deu2003)
標準誤差の推定（2）
続いて、男女比を推定する場合を考えよう。これもintsvyを使えば容易である。
# P.110 Table 7.6に対応 pisa.table(variable="ST03Q01", data=deu2003)
背後にある基本的な発想は変わらない。念の為、Rで細かく計算してみよう。
# 男子=0、女子=1のダミー変数を作成する
library(car)
deu2003$gender <- recode(deu2003$ST03Q01, "2=0")
# Final Weightによる平均値。この値が最終的な推定値になる。
f_weight2 <- weighted.mean(deu2003$gender, deu2003$W_FSTUWT, na.rm=T) * 100
# Replicate Weightによる平均値を処理する関数
deu03_gender <- function(x){
weighted.mean(deu2003$gender, x, na.rm=T)
}
# P.111 Table 7.7
# 80回のReplicate Weightによる平均値
round(apply(deu2003[320:399], 2, deu03_gender) * 100, 2)
# Final WeightとReplicate Weightの差の2乗
round((apply(deu2003[320:399], 2, deu03_gender) * 100 - f_weight2) ^ 2, 2)
# 標準誤差の2乗を計算
sum((apply(deu2003[320:399], 2, deu03_gender) * 100 - f_weight2) ^ 2) / 20
# 最終的な推定値
sam_var2 <- sum((apply(deu2003[320:399], 2, deu03_gender) * 100 - f_weight2) ^ 2) / 20
se2 <- sqrt(sam_var2)
c(f_weight2, se2)
男女ごとに、学年（ST01Q01）を見てみよう。ドイツは留年が存在するため、同じ15歳といえど、異なる学年に在籍することが多い。
# P.112 Table 7.8 deu2003$ST01Q01[deu2003$ST01Q01 > 90] <- NA pisa.table(variable="ST01Q01", by="ST03Q01", data=deu2003)
標準誤差の推定（3）
最後に回帰分析の推定値を考えよう。基本的な発想は同じなので、intsvyによる出力のみ紹介する。
従属変数に、30歳時点で生徒が期待する職業の社会経済指標（BSMJ）をとり、独立変数に、家族の社会経済指標（HISEI）とジェンダー（男子=0、女子=1）の回帰分析を考える。
# P.113 Table 7.9
deu2003$BSMJ[deu2003$BSMJ > 95] <- NA
pisa.reg(y="BSMJ", x=c("HISEI", "gender"), data=deu2003)
いずれもSPSSやSASによる推定値と一致していることが確認できる。
なお、SASによる分析を行ったことはないのだが、PISAを分析する際は、SPSSよりRのほうが高速かつ安定している。
SPSSのマクロにはなにかバグがあるらしく、ときどき原因不明のエラーがでることがある。また、SPSSでPISAの推定を繰り返すと、シンタックスは正しいのに推定値がおかしくなることがあるため、分析に使うには怖いところがある。
Rの難点は、冒頭でも触れたように、あまり大きなファイルを扱うことができないという点だろう。
PISAはファイルサイズがかなり大きいので、これは深刻な問題である。この点をどう解決するかについては、また別の機会に触れることにする。
