Anatomy of PISA chapter3
Daiki Nakamura
2020/10/12
このドキュメントでは、袰岩・篠原・篠原（2019）『PISA調査の解剖』の第3章のコードを再現する。
なお、筆者の環境に合わせてコードを一部加筆修正している。
Microsoft R Open 4.0.2 を使用。
データの用意
2章で作成したCSVファイルを読み込む。
#read csv file
JPN2012 <- read.csv("JPN2012.csv", header = TRUE)スクリプト3-1 パッケージ「survey」のインストール
surveyパッケージをインストールする。
#Install packages
install.packages("survey")## package 'survey' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
## C:\Users\turid\AppData\Local\Temp\RtmpWekdZ7\downloaded_packages#Load library
library(survey)標本誤差、測定誤差、標準誤差
PISA調査の場合、生徒の得点が関係する推定値はすべてPVを伴っており、標本誤差と測定誤差を合わせたものを標準誤差と呼んでいる。それぞれの計算式は以下の通りである。 式2の第1項の分母でー１をしているのは不偏分散と同じ。 式3の第2項より、PVの数が増えるほど標準誤差が相対的に小さくなる。
\[
PVを用いた推定値の標本誤差=\sqrt{\frac{1}{PVの数}\times \sum{各PVの標本誤差}^2}
\\
PVの測定誤差=\sqrt{\frac{1}{PVの数-1}\times \sum{(各PVの推定値-推定値の平均)^2}}
\\
PVを用いた推定値の標準誤差=\sqrt{PVを用いた標本誤差^2+(1+\frac{1}{PVの数})\times 測定誤差^2}
\]
スクリプト3-2 単純無作為抽出を仮定した場合の標準誤差
2章では重み付けを行わずPV1のみを用いた標本誤差の計算を行ったが、ここでは重み付けを行い5つのPVを用いた計算を行う。本文のコードとは異なり、標本誤差と標準誤差を算出する関数を作成している。
数学的リテラシーの平均得点は536.41、標本誤差は1.20、標準誤差は1.21である。
書籍P.102に記載の標準誤差は1.22となっているが、正しくは1.21の誤りだと考えられる。
#Simple random sampling
sam_d <- svydesign(id = ~0, data = JPN2012, weights = ~W_FSTUWT)
#Calculate mean
p1 <- svymean(~PV1MATH, sam_d)
p2 <- svymean(~PV2MATH, sam_d)
p3 <- svymean(~PV3MATH, sam_d)
p4 <- svymean(~PV4MATH, sam_d)
p5 <- svymean(~PV5MATH, sam_d)
mean(c(p1, p2, p3, p4, p5))## [1] 536.4069#Calculate sampling error
samplerror <- function(sam_d){
p1 <- svymean(~PV1MATH, sam_d)
p2 <- svymean(~PV2MATH, sam_d)
p3 <- svymean(~PV3MATH, sam_d)
p4 <- svymean(~PV4MATH, sam_d)
p5 <- svymean(~PV5MATH, sam_d)
se <- sqrt(mean(c(SE(p1)^2, SE(p2)^2, SE(p3)^2, SE(p4)^2, SE(p5)^2)))
return(se)
}
samplerror(sam_d)## [1] 1.196722#Calculate standard error
standerror <- function(sam_d){
p1 <- svymean(~PV1MATH, sam_d)
p2 <- svymean(~PV2MATH, sam_d)
p3 <- svymean(~PV3MATH, sam_d)
p4 <- svymean(~PV4MATH, sam_d)
p5 <- svymean(~PV5MATH, sam_d)
samplerror <- mean(c(SE(p1)^2, SE(p2)^2, SE(p3)^2, SE(p4)^2, SE(p5)^2))
measurerror <- var(c(p1, p2, p3, p4, p5))
se <- sqrt(samplerror+(1+1/5)*measurerror)
return(se)
}
standerror(sam_d)## [1] 1.21181上記スクリプトの冒頭1行を書き換えることで、様々な標本抽出法に対応したデータと標準誤差を計算できる。 以下、データの読み込み方のみ紹介。
スクリプト3-3 二段抽出を仮定した場合の標準誤差
#two stage sampling
sam_d <- svydesign(id = ~SCHOOLID, data = JPN2012, weights = ~W_FSTUWT)
samplerror(sam_d)## [1] 5.008712standerror(sam_d)## [1] 5.012338スクリプト3-4 層化抽出を仮定した場合の標準誤差
#stratified sampling
sam_d <- svydesign(id = ~0, strata = ~STRATUM, data = JPN2012, weights = ~W_FSTUWT)
samplerror(sam_d)## [1] 1.171785standerror(sam_d)## [1] 1.18719スクリプト3-5 層化ニ段抽出を仮定した場合の標準誤差
#stratified two stage sampling
sam_d <- svydesign(id = ~SCHOOLID, strata = ~STRATUM, data = JPN2012, weights = ~W_FSTUWT)
samplerror(sam_d)## [1] 4.842416standerror(sam_d)## [1] 4.846166複雑な標本抽出法における標本誤差の計算方法
線形化法、BRR法、ブートストラップ法、ジャックナイフ法、などがある。
スクリプト3-6 ブートストラップ法による標準誤差の計算
\[ ブートストラップ法の標本誤差=\sqrt{\frac{\sum(複製の推定値-標本の推定値)^2}{複製数-1}} \]
#Bootstrap method
sam_d <- svydesign(id = ~SCHOOLID, strata = ~STRATUM, data = JPN2012, weights = ~W_FSTUWT)
sam_d <- as.svrepdesign(sam_d, type = "subbootstrap", replicates = 1000)
samplerror(sam_d)## [1] 4.939538standerror(sam_d)## [1] 4.943215スクリプト3-7 ジャックナイフ法による標準誤差の計算
\[ JKn法の標本誤差=\sqrt{\sum_{s=1}^{層の数}\left\{\frac{層内の要素数-1}{層内の要素数}\times \sum(複製の推定値-標本の推定値)^2 \right\}} \]
#JackKnife method
sam_d <- svydesign(id = ~SCHOOLID, strata = ~STRATUM, data = JPN2012, weights = ~W_FSTUWT)
sam_d <- as.svrepdesign(sam_d, type = "JKn")
samplerror(sam_d)## [1] 4.843792standerror(sam_d)## [1] 4.847542BRR法とFayの修正法の標本誤差
\[ BRR法の標本誤差=\sqrt{\frac{\sum(複製の推定値-標本の推定値)^2}{複製数}} \\ Fayによる修正法の標本誤差=\sqrt{\frac{1}{(1-修正値)^2}\times \frac{\sum(複製の推定値-標本の推定値)^2}{複製数}} \]
スクリプト3-8 疑似層の作成
JPN2012 <- JPN2012[order(JPN2012$STIDSTD),]
JPN2012$p_STR <- JPN2012$WVARSTRR
JPN2012$p_SCH <- as.character(JPN2012$SCHOOLID)
#Change school ID and pseudo stratum ID
JPN2012[JPN2012$SCHOOLID=="125",]$p_STR <- 81
JPN2012[JPN2012$SCHOOLID=="138",]$p_STR <- 81
JPN2012[JPN2012$SCHOOLID=="15",]$p_STR <- 82
JPN2012[JPN2012$SCHOOLID=="56",]$p_STR <- 82
JPN2012[1669:1685,]$p_SCH <- "192" #SCHOOLID 50
JPN2012[5899:5916,]$p_SCH <- "193" #SCHOOLID 178
JPN2012[JPN2012$SCHOOLID=="152",]$p_STR <-83
JPN2012[JPN2012$SCHOOLID=="57",]$p_STR <-83
JPN2012[JPN2012$SCHOOLID=="136",]$p_STR <-84
JPN2012[4468:4484,]$p_SCH <- "194" #SCHOOLID 136
JPN2012[JPN2012$SCHOOLID=="109",]$p_STR <- 85
JPN2012[JPN2012$SCHOOLID=="60",]$p_STR <- 85
JPN2012[1585:1601,]$p_SCH <- "195" #SCHOOLID 48
JPN2012[5666:5680,]$p_SCH <- "196" #SCHOOLID 171
JPN2012[JPN2012$SCHOOLID=="34",]$p_STR <- 86
JPN2012[JPN2012$SCHOOLID=="69",]$p_STR <- 86
JPN2012[2019:2035,]$p_SCH <- "197" #SCHOOLID 61
JPN2012[JPN2012$SCHOOLID=="67",]$p_STR <- 87
JPN2012[JPN2012$SCHOOLID=="37",]$p_STR <- 87
JPN2012[2054:2069,]$p_SCH <- "198" #SCHOOLID 62
JPN2012[JPN2012$SCHOOLID=="24",]$p_STR <- 88
JPN2012[JPN2012$SCHOOLID=="83",]$p_STR <- 88
JPN2012[203:220,]$p_SCH <- "199" #SCHOOLID 7
JPN2012[JPN2012$SCHOOLID=="92",]$p_STR <- 89
JPN2012[3045,]$p_SCH <- "200" #SCHOOLID 136
JPN2012[JPN2012$SCHOOLID=="122",]$p_STR <- 90
JPN2012[JPN2012$SCHOOLID=="135",]$p_STR <- 90
JPN2012[JPN2012$SCHOOLID=="77",]$p_STR <- 91
JPN2012[JPN2012$SCHOOLID=="74",]$p_STR <- 91
JPN2012[JPN2012$SCHOOLID=="68",]$p_STR <- 92
JPN2012[JPN2012$SCHOOLID=="31",]$p_STR <- 92
JPN2012[729:745,]$p_SCH <- "201" #SCHOOLID 22
JPN2012[JPN2012$SCHOOLID=="47",]$p_STR <- 93
JPN2012[JPN2012$SCHOOLID=="40",]$p_STR <- 93
JPN2012[JPN2012$SCHOOLID=="45",]$p_STR <- 94
JPN2012[JPN2012$SCHOOLID=="95",]$p_STR <- 94
JPN2012[4606:4622,]$p_SCH <- "202" #SCHOOLID 140
JPN2012[JPN2012$SCHOOLID=="161",]$p_STR <- 95
JPN2012[JPN2012$SCHOOLID=="93",]$p_STR <- 95
JPN2012[4081:4097,]$p_SCH <- "203" #SCHOOLID 124
JPN2012[JPN2012$SCHOOLID=="17",]$p_STR <- 96
JPN2012[JPN2012$SCHOOLID=="158",]$p_STR <- 96
JPN2012[JPN2012$SCHOOLID=="29",]$p_STR <- 97
JPN2012[JPN2012$SCHOOLID=="177",]$p_STR <- 97
JPN2012[JPN2012$SCHOOLID=="21",]$p_STR <- 98
JPN2012[JPN2012$SCHOOLID=="114",]$p_STR <- 98
JPN2012[JPN2012$SCHOOLID=="134",]$p_STR <- 99
JPN2012[JPN2012$SCHOOLID=="131",]$p_STR <- 99
JPN2012[JPN2012$SCHOOLID=="186",]$p_STR <- 100
JPN2012[JPN2012$SCHOOLID=="150",]$p_STR <- 100
JPN2012[JPN2012$SCHOOLID=="156",]$p_STR <- 101
JPN2012[5155:5171,]$p_SCH <- "204" #SCHOOLID 156スクリプト3-9 BRR法による標準誤差の計算
#BRR method
sam_d <- svydesign(id = ~p_SCH, strata = ~p_STR, data = JPN2012, weights = ~W_FSTUWT)
sam_d <- as.svrepdesign(sam_d, type = "BRR")
samplerror(sam_d)## [1] 3.629782standerror(sam_d)## [1] 3.634784スクリプト3-10 Fayによる修正法による標準誤差の計算
#Fay`s method
sam_d <- svydesign(id = ~p_SCH, strata = ~p_STR, data = JPN2012, weights = ~W_FSTUWT)
sam_d <- as.svrepdesign(sam_d, type = "Fay", fay.rho = 0.5)
samplerror(sam_d)## [1] 3.628948standerror(sam_d)## [1] 3.633951スクリプト3-11 架空の有限母集団(シミュレーションデータ)の作成
#Make a finit population
#set stratum and size
N_str4 <- 4 #4 strata
N_str100 <- 100 #100 strata
N_sch <- 7000 #7000 schools
N_stu <- 150 #150 students in each schools
ID <- 1:(N_sch*N_stu) #ID number
str4 <- rep(1:N_str4, each=N_sch*N_stu/N_str4) #stratum ID
str100 <- rep(1:N_str100, each=N_sch*N_stu/N_str100) #stratum ID
sch <- rep(1:N_sch, each=N_stu) #school ID
stu <- rep(1:N_stu, times=N_sch) #student ID
#0-100 percentile, delete 0 & 100 percentile
pe_sch <- seq (0, 1, length=N_sch+2)[2:(N_sch+1)]
#make school means
mean_sch <- rep(qnorm(pe_sch, mean=500, sd=100), each=N_stu)
#O-100 percentile, delete 0 & 100 percentile
pe_stu <- rep(seq(0, 1, length=N_stu+2)[2:(N_stu+1)], times=N_sch)
#make scores of all students
score <- qnorm(pe_stu, mean=mean_sch, sd=100)
#make a data
population <- data.frame(ID, str4, str100, sch, stu, score)
#calculate mean & sd
mean(population$score)## [1] 500sqrt(var(population$score)*(1050000-1)/1050000)## [1] 139.3332sd(population$score)## [1] 139.3333#Fit linear mixed model
library(lme4)
m_model1 <- lmer(score~(1|sch), data = population)
summary(m_model1)## Linear mixed model fit by REML ['lmerMod']
## Formula: score ~ (1 | sch)
## Data: population
##
## REML criterion at convergence: 12632121
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.5652 -0.6867 0.0000 0.6867 2.5652
##
## Random effects:
## Groups Name Variance Std.Dev.
## sch (Intercept) 9916 99.58
## Residual 9499 97.46
## Number of obs: 1050000, groups: sch, 7000
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 500.000 1.194 418.8#Histogram
hist(population$score, breaks = seq(-200, 1200, 1))
スクリプト3-12 多段抽出（復元あり）のシミュレーション
#Install packages
install.packages("foreach")## package 'foreach' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
## C:\Users\turid\AppData\Local\Temp\RtmpWekdZ7\downloaded_packagesinstall.packages("doParallel")## package 'doParallel' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
## C:\Users\turid\AppData\Local\Temp\RtmpWekdZ7\downloaded_packagesinstall.packages("sampling")## package 'sampling' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
## C:\Users\turid\AppData\Local\Temp\RtmpWekdZ7\downloaded_packages#Load library
library(foreach)
library(doParallel)
library(sampling)
library(survey)
#sampling simulation
ii <- 100 #number of simulation. original textbook is 10000.
n_psu <- 200 #number of psu
s_psu <- 35 #target cluster size
#two level sampling with replacement
str_brr <- rep(1:(n_psu/2), each=s_psu*2) #pseudo stratum for BRR
weight <- rep(1, times=n_psu*s_psu) #weight for mean calculation
q_sch <- rep(1:n_psu, each=s_psu) #school ID for replacement
mean_sample <- NULL #reset sample mean
cl <- makeCluster(detectCores())
registerDoParallel(cl)
mean_sample <- foreach(i = 1:ii, .combine = "rbind",
.packages = c("sampling", "survey")) %dopar% {
sample_score <- NULL #reset sample score
#school sampling
sample_sch <- sample(1:N_sch, n_psu, replace = TRUE)
#student sampling within each sampled schools
for(r in sample_sch){
sample_score <- c(sample_score,
sample(population[population$sch==r,]
$score, s_psu, replace = TRUE))
}
#make a sample data
sample_brr <- data.frame(q_sch, sample_score, str_brr, weight)
#BRR method
design0_brr <- svydesign(id=~q_sch, weights = ~weight,
strata = ~str_brr, data = sample_brr)
design1_brr <- as.svrepdesign(design0_brr, type = "BRR")
x <- svymean(~sample_score, design1_brr) #compute a sample mean
c(x, SE(x)) #mean and SE
}
stopCluster(cl)
write.table(mean_sample, "simulation01.txt", quote = F, col.names = F) #output
mean(mean_sample[,2])## [1] 7.104747スクリプト3-13 多段抽出（復元なし）のシミュレーション
#two level sampling without replacement
weight <- rep(1, times=n_psu*s_psu) #weight for mean calculation
mean_sample <- NULL #reset sample mean
stage_st <- list("cluster", "cluster") #set a sampling design
stage_na <- list("sch", "stu") #set stratum and cluster
size_s <- list(size1=n_psu, size2=rep(s_psu, times=n_psu)) #set sample size
method_s <- list("srswor", "srswor") #set sampling method
cl <- makeCluster(detectCores())
registerDoParallel(cl)
mean_sample <- foreach(i = 1:ii, .combine = "rbind",
.packages = c("sampling", "survey")) %dopar% {
#eliminate the effect of systematic sampling
str_brr <- rep(1:(n_psu/2), each=2)
str_brr <- rep(sample(str_brr, n_psu, replace = FALSE), each=s_psu)
#select number
ID_sample <- mstage(population, stage = stage_st, varnames = stage_na,
size = size_s, method = method_s)
#make a sample data
sample <- getdata(population, ID_sample)[[2]]
#BRR method
sample_brr <- cbind(sample, str_brr, weight)
design0_brr <- svydesign(id=~sch, weights = ~weight,
strata = ~str_brr, data = sample_brr)
design1_brr <- as.svrepdesign(design0_brr, type = "BRR")
x <- svymean(~score, design1_brr) #compute a sample mean
c(x, SE(x)) #mean and SE
}
stopCluster(cl)
write.table(mean_sample, "simulation02.txt", quote = F, col.names = F) #output
mean(mean_sample[,2])## [1] 7.146815スクリプト3-14 層化二段抽出（4層）のシミュレーション
#stratified two level sampling with 4 strata
weight <- rep(1, times=n_psu*s_psu) #weight for mean calculation
mean_sample <- NULL #reset sample mean
stage_st <- list("stratified", "cluster", "cluster") #set a sampling design
stage_na <- list("str4", "sch", "stu") #set stratum and cluster
size_s <- list(size1=table(population$str4),
size2=rep(n_psu/N_str4, times=N_str4),
size3=rep(s_psu, times=n_psu))#set sample size
method_s <- list("", "srswor", "srswor") #set sampling method
cl <- makeCluster(detectCores())
registerDoParallel(cl)
mean_sample <- foreach(i = 1:ii, .combine = "rbind",
.packages = c("sampling", "survey")) %dopar% {
#eliminate the effect of systematic sampling
str1_brr <- rep(1:(n_psu/N_str4/2), each=2)
str1_brr <- sample(str1_brr, n_psu/N_str4, replace = FALSE)
str2_brr <- rep((n_psu/N_str4/2+1):(n_psu/N_str4/2*2), each=2)
str2_brr <- sample(str2_brr, n_psu/N_str4, replace = FALSE)
str3_brr <- rep((n_psu/N_str4/2*2+1):(n_psu/N_str4/2*3), each=2)
str3_brr <- sample(str3_brr, n_psu/N_str4, replace = FALSE)
str4_brr <- rep((n_psu/N_str4/2*3+1):(n_psu/N_str4/2*4), each=2)
str4_brr <- sample(str4_brr, n_psu/N_str4, replace = FALSE)
str_brr <- rep(c(str1_brr, str2_brr, str3_brr, str4_brr), each=s_psu)
#select number
ID_sample <- mstage(population, stage = stage_st, varnames = stage_na,
size = size_s, method = method_s)
#make a sample data
sample <- getdata(population, ID_sample)[[3]]
#BRR method
sample_brr <- cbind(sample, str_brr, weight)
design0_brr <- svydesign(id=~sch, weights = ~weight,
strata = ~str_brr, data = sample_brr)
design1_brr <- as.svrepdesign(design0_brr, type = "BRR")
x <- svymean(~score, design1_brr) #compute a sample mean
c(x, SE(x)) #mean and SE
}
stopCluster(cl)
write.table(mean_sample, "simulation03.txt", quote = F, col.names = F) #output
mean(mean_sample[,2])## [1] 2.800764スクリプト3-15 層化二段抽出（100層）のシミュレーション
#stratified two level sampling with 2 psu strata
str_brr <- rep(1:(n_psu/2), each=s_psu*2) #pseudo stratum for BRR
weight <- rep(1, times=n_psu*s_psu) #weight for mean calculation
mean_sample <- NULL #reset sample mean
stage_st <- list("stratified", "cluster", "cluster") #set a sampling design
stage_na <- list("str100", "sch", "ID") #set stratum and cluster
size_s <- list(size1=table(population$str100),
size2=rep(n_psu/N_str100, times=N_str100),
size3=rep(s_psu, times=n_psu)) #set sample size
method_s <- list("", "srswor", "srswor") #set sampling method
prob <- list("", rep(table(population$sch), times=table(population$sch))) #set probability
cl <- makeCluster(detectCores())
registerDoParallel(cl)
mean_sample <- foreach(i = 1:ii, .combine = "rbind",
.packages = c("sampling", "survey")) %dopar% {
#select number
ID_sample <- mstage(population, stage = stage_st, varnames = stage_na,
size = size_s, method = method_s, pik = prob, description = FALSE)
#make a sample data
sample <- getdata(population, ID_sample)[[3]]
#BRR method
sample_brr <- cbind(sample, str_brr, weight)
design0_brr <- svydesign(id=~sch, weights = ~weight,
strata = ~str_brr, data = sample_brr)
design1_brr <- as.svrepdesign(design0_brr, type = "BRR")
x <- svymean(~score, design1_brr) #compute a sample mean
c(x, SE(x)) #mean and SE
}
stopCluster(cl)
write.table(mean_sample, "simulation04.txt", quote = F, col.names = F) #output
mean(mean_sample[,2])## [1] 1.061727スクリプト3-16 系統抽出のシミュレーション
#stratified two level sampling with systematic sampling
str_brr <- rep(1:(n_psu/2), each=s_psu*2) #pseudo stratum for BRR
weight <- rep(1, times=n_psu*s_psu) #weight for mean calculation
mean_sample <- NULL #reset sample mean
stage_st <- list("stratified", "cluster", "cluster") #set a sampling design
stage_na <- list("str4", "sch", "ID") #set stratum and cluster
size_s <- list(size1=table(population$str4),
size2=rep(N_stu*n_psu/N_str4, times=N_str4),
size3=rep(s_psu, times=n_psu)) #set sample size
method_s <- list("", "systematic", "srswor") #set sampling method
prob <- list("", rep(table(population$sch), times=table(population$sch))) #set probability
cl <- makeCluster(detectCores())
registerDoParallel(cl)
mean_sample <- foreach(i = 1:ii, .combine = "rbind",
.packages = c("sampling", "survey")) %dopar% {
#select number
ID_sample <- mstage(population, stage = stage_st, varnames = stage_na,
size = size_s, method = method_s, pik = prob, description = FALSE)
#make a sample data
sample <- getdata(population, ID_sample)[[3]]
#BRR method
sample_brr <- cbind(sample, str_brr, weight)
design0_brr <- svydesign(id=~sch, weights = ~weight,
strata = ~str_brr, data = sample_brr)
design1_brr <- as.svrepdesign(design0_brr, type = "BRR")
x <- svymean(~score, design1_brr) #compute a sample mean
c(x, SE(x)) #mean and SE
}
stopCluster(cl)
write.table(mean_sample, "simulation05.txt", quote = F, col.names = F) #output
mean(mean_sample[,2])## [1] 1.097657