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

投稿日：2017-05-19

Rのintsvy packageを使ってPISA2015を分析する方法です。

データの準備

まず、OECDのホームページからPISA2015のデータをダウンロードする。PISA2015からSPSSのデータファイルが公開されており、【SPSS (TM) Data Files】からzipファイルをダウンロードできる。ダウンロード後に展開すれば、【CY6_MS_CMB_STU_QQQ.sav】や【CY6_MS_CMB_STU_QQQ.sav】というファイルが得られるので、これを使えばよい。

なお、PISA2015のデータファイルは1GBを超える大きさなので、これを毎回読み込むのは大変である。そこで以下では、国別のcsvファイルに変換して保存し、必要な国のデータだけ利用することにする。

library(haven)

filename <- c("CY6_MS_CMB_SCH_QQQ.sav",
              "CY6_MS_CMB_STU_QQQ.sav")

dat <- lapply(filename, function(x) read_spss(x))

cnt <- names(table(dat[[1]]$CNT))
for (i in cnt){
  sch <- subset(dat[[1]], dat[[1]]$CNT == i)
  stu <- subset(dat[[2]], dat[[2]]$CNT == i)
  write.table(sch, file = paste0("2015sch", i, ".csv"),
              quote = F, sep = ",", row.names = FALSE)
  write.table(stu, file = paste0("2015stu", i, ".csv"),
              quote = F, sep = ",", row.names = FALSE)
}

これで【2015stuJPN.csv】や【2015schUSA.csv】といったcsvファイルが得られる。

intsvyによる分析

intsvy packageの初期設定では、PISA2015は分析できない。なぜならPISA2015では、（1）PVsが5から10に、（2）Replicate Weightの変数名がW_FSTRからW_FSTURWTに、といった変更が加えられたからである。

現時点（2017年5月19日）では、githubにこの変更を反映したintsvyが公開されている。しかし、CRANのintsvyには反映されていない。

とは言え、intsvy関数の引数（config）をいじればCRANのintsvyでも対応可能である。具体的な方法は以下。

library(intsvy)
jpn <- read.csv("2015stuJPN.csv")

# 科学リテラシーの習熟度。Technical Reportに記載がある
scie <- c(260.54, 334.94, 409.54, 484.14, 558.73, 633.33, 707.93)

pisa2015_conf <- list(variables = list(pvlabelpref = "PV",
                                       pvlabelsuff = "SCIE",
                                       weightFinal = "W_FSTUWT",
                                       weightBRR = "W_FSTURWT"),
                      parameters = list(cutoffs = scie,
                                        percentiles = c(5, 10, 25, 75, 90, 95),
                                        PVreps = 10,
                                        BRRreps = 80,
                                        weights = "BRR",
                                        replication_scheme = "pisa")
                      )

# 回帰分析
intsvy.reg.pv(pvlabel = "SCIE", x = "ESCS", data = jpn, config = pisa2015_conf)

# 度数分布（性別：女=1，男=2）
intsvy.table(variable = "ST004D01T", data = jpn, config = pisa2015_conf)

# 習熟度（男女別）
intsvy.ben.pv(pvlabel = "SCIE", by = "ST004D01T", data = jpn, config = pisa2015_conf)