日々の記録簿

日々のできごとの記録

PISA

RによるPISA Data Analysis Manualの再現(第7章: 標準誤差の推定)

投稿日:

学力調査の基本を勉強する上で大変参考になる、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はファイルサイズがかなり大きいので、これは深刻な問題である。この点をどう解決するかについては、また別の機会に触れることにする。

スポンサーリンク

スポンサーリンク

-PISA

執筆者:


comment

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

関連記事

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

三重大の奥村晴彦氏のホームページにて、RでPISAデータを読み込む方法(PISAデータを読む)が公開されていました。 こちらは、Rubyで読み込んでいるようです(ちなみに私はRubyはサッパリです・・ …

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

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

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

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

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

RによるPISAの分析(データの準備)で公開したpythonコードですが、その後、PISA2000のspsファイルに対応していなかったとか、すべての国のデータをまとめて変換したほうが便利だとか、ファイ …

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

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