このドキュメントでは、袰岩・篠原・篠原（2019）『PISA調査の解剖』の第2章のコードを再現する。
なお、筆者の環境に合わせてコードを一部加筆修正している。
書籍同様、不要になったデータはrm()で消去しながら作業を進めている。
ワーキングディレクトリを指定する。
ここでは、マイドキュメントのRフォルダ内にPISAというフォルダを作成して実行
#Set Working directory and read data file
getwd()## [1] "C:/Users/turid/Documents/R/PISA"setwd("C:/Users/turid/Documents/R/PISA")必要なパッケージとデータをインストールする。
githubからのインストールは時間がかかるので注意。
#Install package
install.packages("devtools", repos = "http://cran.us.r-project.org")## Installing package into 'C:/Users/turid/Documents/R/win-library/4.0'
## (as 'lib' is unspecified)## package 'devtools' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
## C:\Users\turid\AppData\Local\Temp\RtmpYRNbMQ\downloaded_packages#Load library
library(devtools)## Loading required package: usethis#Install package
install_github("pbiecek/PISA2012lite") #時間がかかる## Skipping install of 'PISA2012lite' from a github remote, the SHA1 (3ae23dd6) has not changed since last install.
## Use `force = TRUE` to force installation#Load library
library(PISA2012lite)ダウンロードしたデータの整形。
#Load data
data("student2012") #生徒質問紙と得点のデータ
data("student2012weights") #生徒に対する重みの情報
data("school2012") #学校質問紙のデータ
#Make dataset of PISA2012
stu2012 <- cbind(student2012[,-552:-555],
student2012weights[,-1:-4]) #データと重みを結合
pisa2012 <- merge(stu2012, school2012,
by=c("CNT","SUBNATIO","STRATUM","OECD","NC","SCHOOLID"),
all.x=TRUE) #生徒質問紙と学校質問紙を結合
rm("stu2012", "student2012", "student2012weights", "school2012")
#Make dataset of Japan
JPN2012 <- pisa2012[pisa2012$CNT=="Japan",]
rm("pisa2012")データフレームに変換し、どのような変数が含まれているか確認する。
#Load data
data("student2012dict")
data("school2012dict")
student2012dict <- as.data.frame(student2012dict)
school2012dict <- as.data.frame(school2012dict)
#View variable list
View(student2012dict)
View(school2012dict)国際的な大規模調査の分析にはintsvyを使うと便利。
ここでは、PISA2012の日本のデータを使って、5つのPVの平均と標準誤差、属性ごとの回答割合を求める。
#Install packages
install.packages("intsvy", repos = "http://cran.us.r-project.org")## Installing package into 'C:/Users/turid/Documents/R/win-library/4.0'
## (as 'lib' is unspecified)## package 'intsvy' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
## C:\Users\turid\AppData\Local\Temp\RtmpYRNbMQ\downloaded_packages#Load library
library(intsvy)
#Culculate mean and SD
pisa.mean.pv(pvlabel = "MATH", data = JPN2012) #平均得点と標準偏差を計算 MATH/READ/SCIE## Freq Mean s.e. SD s.e
## 1 6351 536.41 3.59 93.52 2.19#Culculate percent
pisa.table(variable = "STRATUM", data = JPN2012)## STRATUM Freq Percentage Std.err.
## 688 JPN0101 3137 48.08 0.77
## 689 JPN0202 1355 20.64 0.52
## 690 JPN0203 1532 26.88 0.90
## 691 JPN0204 327 4.41 0.20pisa.table(variable = "SC03Q01", data = JPN2012)## SC03Q01 Freq Percentage Std.err.
## 2 Small Town 95 1.64 0.98
## 3 Town 1689 25.75 2.69
## 4 City 3144 49.42 3.40
## 5 Large City 1423 23.20 3.31復元単純無作為抽出の場合、標本誤差は以下の式で求められる。 \[ 復元単純無作為抽出での標本誤差 = \sqrt \frac{分散 v}{サンプルサイズ n} \] ここでは、数学的リテラシーのPV1を対象に、復元単純無作為抽出を行う。
all_var <- var(JPN2012$PV1MATH) #分散
stuNum <- length(JPN2012$PV1MATH) #サンプルサイズ
sqrt(all_var/stuNum) #平方根## [1] 1.181374非復元単純無作為抽出の場合、標本誤差は以下の式で求められる。
有限母集団修正（１－抽出率）を積算しているところがポイント。
抽出率が高くなれば標本誤差は小さくなる。
全数調査などで抽出率が1.0になる場合、復元単純無作為抽出の標本誤差と一致する。 \[
非復元単純無作為抽出での標本誤差 = \sqrt {\frac{分散 v}{サンプルサイズ n} \times {\left(1- \frac{サンプルサイズn}{母集団サイズN}\right)}}
\] ## スクリプト2-7 標本誤差の計算（二段抽出） 無限母集団における二段抽出では、標本誤差を以下の式で計算できる。
得点の分散を学校間分散と学校内分散に分解しているところがポイント。
TCS（各学校で抽出する生徒数ｎ）を増やすと、第２項のみが小さくなるのに対して、学校数kを増やすと2つの項がどちらも小さくなるので、効果的に標本誤差を小さくできる。
ただし、学校数を増やすことは調査コストの増加にもつながる。 \[
二段抽出での標本誤差 = \sqrt {\frac{学校間分散 v_b}{学校数 k} + \frac{学校内分散 v_w}{生徒数 n}}
\] では、関数lmer()を用いて、学校間分散と学校内分散を計算してみよう。
#Install packages
install.packages("lme4", repos = "http://cran.us.r-project.org")## Installing package into 'C:/Users/turid/Documents/R/win-library/4.0'
## (as 'lib' is unspecified)## package 'lme4' successfully unpacked and MD5 sums checked## Warning: cannot remove prior installation of package 'lme4'## Warning in file.copy(savedcopy, lib, recursive = TRUE): C:
## \Users\turid\Documents\R\win-library\4.0\00LOCK\lme4\libs\x64\lme4.dll の C:
## \Users\turid\Documents\R\win-library\4.0\lme4\libs\x64\lme4.dll へのコピーに問題
## があります: Permission denied## Warning: restored 'lme4'##
## The downloaded binary packages are in
## C:\Users\turid\AppData\Local\Temp\RtmpYRNbMQ\downloaded_packages#Load library
library(lme4)## Loading required package: Matrix## Registered S3 methods overwritten by 'lme4':
## method from
## cooks.distance.influence.merMod car
## influence.merMod car
## dfbeta.influence.merMod car
## dfbetas.influence.merMod car#Fit linear mixed model
m_model1 <- lmer(PV1MATH~(1|SCHOOLID), data = JPN2012)
summary(m_model1)## Linear mixed model fit by REML ['lmerMod']
## Formula: PV1MATH ~ (1 | SCHOOLID)
## Data: JPN2012
##
## REML criterion at convergence: 71588.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.1009 -0.6316 0.0161 0.6439 3.7507
##
## Random effects:
## Groups Name Variance Std.Dev.
## SCHOOLID (Intercept) 4908 70.06
## Residual 4122 64.20
## Number of obs: 6351, groups: SCHOOLID, 191
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 534.065 5.137 104PISA調査の平均得点について、抽出法別に標本後誤差を算出すると、層化に一定の効果があることが分かる。
書籍中の表２－９を以下に示す。
| 単純無作為抽出 | 二段抽出（層化なし） | 層化二段抽出（BRR法） | |
|---|---|---|---|
| 標本誤差 | 1.174 | 5.104 | 3.587 |
多段抽出を行わない場合、平均得点の標本誤差は以下の式で求まる。 \[
層化抽出での標本誤差 = \sqrt {\sum_{s=1}^{層の数} \left\{\left(\frac{各層の母集団のサイズN_s}{母集団サイズN}\right)^2 \times \frac{各層内での分散v_s}{各層のサンプルサイズn_s} \times {\left(1- \frac{各層のサンプルサイズn_s}{各層の母集団サイズN_s}\right)} \right\}}
\] ここで、上式の第3項は有限母集団修正であり、無限母集団の場合は無視できる。
また、各層のサンプルサイズはサンプルサイズと各層の比率の積で表現できる。 \[
層化抽出での標本誤差 =\sqrt {\sum_{s=1}^{層の数} \left\{\left(\frac{各層の母集団のサイズN_s}{母集団サイズN}\right)^2 \times \frac{各層内での分散v_s}{各層のサンプルサイズn_s}\right\}}
\\
=\sqrt {\sum_{s=1}^{層の数} \left\{\left(\frac{各層の母集団のサイズN_s}{母集団サイズN}\right)^2 \times \frac{各層内での分散v_s}{{サンプルサイズn} \times {\frac{各層の母集団サイズN_s}{母集団サイズN}}}\right\}}
\\
=\sqrt {\sum_{s=1}^{層の数} \left(\frac{各層の母集団のサイズN_s}{母集団サイズN} \times \frac{各層内での分散v_s}{サンプルサイズn}\right)}
\] さらに、第1項は各層内の分散を足し合わせる際の重みであると解釈できるが、重み付けられた分散の総和は級内分散を意味する。 \[
層化抽出での標本誤差=\sqrt {\sum_{s=1}^{層の数} \left(重み \times \frac{各層内での分散v_s}{サンプルサイズn}\right)}
\\
=\sqrt{\frac{級内分散v_w}{サンプルサイズn}} =\sqrt{\frac{分散v-級間分散v_b}{サンプルサイズn}}
\] 上式より、層の平均得点の分散（級間分散）が大きいほど標本誤差が小さくなることが分かる。
PISA調査で主層を分けるのに日本は設置者（公立／国立・私立）×学科の種類（普通科／専門学科）＝４層を用いている。
一方、諸外国では学校所在地域によって主層を分けている国も多い。
では、日本のデータの場合、どちらの分け方が適切だろうか。
それぞれの分け方で層化抽出の標本誤差を計算してみよう。
m_model2 <- lmer(PV1MATH~(1|STRATUM), data = JPN2012)
m_model3 <- lmer(PV1MATH~(1|SC03Q01), data = JPN2012)
summary(m_model2)## Linear mixed model fit by REML ['lmerMod']
## Formula: PV1MATH ~ (1 | STRATUM)
## Data: JPN2012
##
## REML criterion at convergence: 75452.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8784 -0.6482 0.0460 0.6882 3.3307
##
## Random effects:
## Groups Name Variance Std.Dev.
## STRATUM (Intercept) 674.1 25.96
## Residual 8441.0 91.88
## Number of obs: 6351, groups: STRATUM, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 527.91 13.08 40.37summary(m_model3)## Linear mixed model fit by REML ['lmerMod']
## Formula: PV1MATH ~ (1 | SC03Q01)
## Data: JPN2012
##
## REML criterion at convergence: 75461
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8312 -0.6648 0.0327 0.7043 3.3753
##
## Random effects:
## Groups Name Variance Std.Dev.
## SC03Q01 (Intercept) 3926 62.65
## Residual 8447 91.91
## Number of obs: 6351, groups: SC03Q01, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 507.82 31.43 16.16desi_eff2 <- m_model2２つのモデルを比較するための指標として、デザイン効果を算出する。
デザイン効果は、単純無作為抽出との比較指標であり、1に近いほど標本誤差が単純無作為抽出に近く、1を下回るほど標本誤差が単純無作為抽出より小さいことを意味する。
相対的にデザイン効果指標が小さい方が層化の効果が高いと判断できる。
復元抽出（無限母集団）の場合、デザイン効果指標は以下の式で計算できる。
級間分散が大きいほどデザイン効果指標も小さくなる（層化の効果が高い）ことが分かる。 \[
デザイン効果=\frac{層化抽出での標本誤差^2}{単純無作為抽出での標本誤差^2}
\\
=\frac{級内分散v_w}{サンプルサイズn} \div \frac{分散v}{サンプルサイズn}
\\
=\frac{級内分散v_w}{級間分散v_b+級内分散v_w}
\] PV1MATHを用いてデザイン効果を計算すると以下の通り。
var2 <- as.data.frame(VarCorr(m_model2))
var2$vcov[2]/(var2$vcov[1]+var2$vcov[2])## [1] 0.9260449var3 <- as.data.frame(VarCorr(m_model3))
var3$vcov[2]/(var3$vcov[1]+var3$vcov[2])## [1] 0.6827192層ごとに習熟度レベル別生徒割合を求める。
学校所在地の方が、層ごとのばらつきが大きいことが分かる。
pisa.ben.pv(pvlabel = "MATH", by = "STRATUM", data = JPN2012)## STRATUM Benchmarks Percentage Std. err.
## 1 JPN0101 <= 357.77 3.06 0.74
## 2 JPN0101 (357.77, 420.07] 6.31 0.74
## 3 JPN0101 (420.07, 482.38] 14.30 1.05
## 4 JPN0101 (482.38, 544.68] 23.10 1.44
## 5 JPN0101 (544.68, 606.99] 25.85 1.25
## 6 JPN0101 (606.99, 669.3] 18.31 1.02
## 7 JPN0101 > 669.3 9.05 1.10
## 8 JPN0202 <= 357.77 4.51 1.04
## 9 JPN0202 (357.77, 420.07] 11.36 2.06
## 10 JPN0202 (420.07, 482.38] 24.02 2.38
## 11 JPN0202 (482.38, 544.68] 30.18 1.68
## 12 JPN0202 (544.68, 606.99] 20.32 1.95
## 13 JPN0202 (606.99, 669.3] 7.95 1.51
## 14 JPN0202 > 669.3 1.65 0.72
## 15 JPN0203 <= 357.77 2.02 0.99
## 16 JPN0203 (357.77, 420.07] 7.10 1.42
## 17 JPN0203 (420.07, 482.38] 15.32 2.03
## 18 JPN0203 (482.38, 544.68] 22.90 2.32
## 19 JPN0203 (544.68, 606.99] 23.72 2.01
## 20 JPN0203 (606.99, 669.3] 19.29 2.69
## 21 JPN0203 > 669.3 9.66 2.29
## 22 JPN0204 <= 357.77 4.80 2.15
## 23 JPN0204 (357.77, 420.07] 14.05 4.30
## 24 JPN0204 (420.07, 482.38] 21.85 4.46
## 25 JPN0204 (482.38, 544.68] 26.54 6.32
## 26 JPN0204 (544.68, 606.99] 15.69 3.32
## 27 JPN0204 (606.99, 669.3] 9.31 3.04
## 28 JPN0204 > 669.3 7.76 5.22pisa.ben.pv(pvlabel = "MATH", by = "SC03Q01", data = JPN2012)## SC03Q01 Benchmarks Percentage Std. err.
## 1 Small Town <= 357.77 23.77 4.87
## 2 Small Town (357.77, 420.07] 29.74 5.48
## 3 Small Town (420.07, 482.38] 27.21 3.73
## 4 Small Town (482.38, 544.68] 11.02 6.25
## 5 Small Town (544.68, 606.99] 6.84 1.92
## 6 Small Town (606.99, 669.3] 1.42 1.55
## 7 Small Town > 669.3 0.00 0.00
## 8 Town <= 357.77 3.28 1.01
## 9 Town (357.77, 420.07] 9.05 1.55
## 10 Town (420.07, 482.38] 21.81 2.24
## 11 Town (482.38, 544.68] 27.68 2.03
## 12 Town (544.68, 606.99] 22.02 1.97
## 13 Town (606.99, 669.3] 12.18 2.11
## 14 Town > 669.3 3.98 1.21
## 15 City <= 357.77 2.97 0.81
## 16 City (357.77, 420.07] 7.33 1.11
## 17 City (420.07, 482.38] 16.10 1.36
## 18 City (482.38, 544.68] 24.32 1.37
## 19 City (544.68, 606.99] 25.02 1.56
## 20 City (606.99, 669.3] 16.48 1.41
## 21 City > 669.3 7.79 1.34
## 22 Large City <= 357.77 1.97 0.69
## 23 Large City (357.77, 420.07] 6.33 1.54
## 24 Large City (420.07, 482.38] 12.50 1.99
## 25 Large City (482.38, 544.68] 23.01 2.93
## 26 Large City (544.68, 606.99] 23.90 2.17
## 27 Large City (606.99, 669.3] 20.41 2.95
## 28 Large City > 669.3 11.89 2.75母集団から標本が選ばれる確率がすべての生徒で等しく、欠損値が発生していない場合以外は、データを重み付けして補正する必要がある。
重み付けてによって、１）生徒が抽出される確率が異なる際の推定値のゆがみ（バイアス）を補正し、２）不参加やデータ破損などによる欠損を補完することができる。
PISA調査は国や地域ごとに生徒の抽出確率が異なるので、抽出確率の逆数を使って補正を行う。
PISA調査の重み付けは次の式で計算される。なお、i は学科、j は生徒を表す。
詳細は、PISA2015 technical report part8（P.115~）を参照 \[
W_{ij}=w_{1i}\times t_{1i}\times w_{2ij}\times f_{1j}\times f_{2ij}\times t_{2ij}
\] \[
\begin{align*}
&学校基本ウェイトw_{1i}=\begin{cases}
I/MOS_i & (MOS_i < I) \\
1 & (otherwise)
\end{cases}\\
&学校ウェイト修正係数t_{1i}=\begin{cases}
w_{1i}=I/3*MAX(TCS, MOS_i) & (調査時点での生徒数>3 \times MOS) \\
1 & (otherwise)
\end{cases}\\
&生徒基本ウェイトw_{2ij}=\frac{生徒名簿での生徒数enr_i}{選ばれた生徒数sam_i}\\
&学校不参加補正f_{1j}=\frac{層内の全学校の重み付き生徒数}{層内の参加学校の重み付き生徒数}=\frac{\sum_{k{\in \Omega}}(W_{1k}\times enr_k)}{\sum_{k{\in \Gamma}}(W_{1k}\times enr_k)}\\
&生徒不参加補正f_{2ij}=\frac{グループ内の全生徒の重みの合計}{グループ内の参加した生徒の重みの合計}=\frac{\sum_{k{\in X}}(f_{1i}\times w_{1i} \times w_{2ik})}{\sum_{k{\in \Delta}}(f_{1i}\times w_{1i} \times w_{2ik})}\\
&最終ウェイト修正係数t_{2ij}
\end{align*}
\] 上記の重みを積算した各生徒の最終ウェイトを調査データに追加して使用する。
上記の重みを使うことでゆがみの少ない母集団の推定値を得ることができる。
これまでの作業結果をcsvファイルに書き出す。
#Write csv file
write.csv(JPN2012,"JPN2012.csv", quote = TRUE, row.names = FALSE)
#reaad csv file
JPN2012_2 <- read.csv("JPN2012.csv", header = TRUE)PISA2012の日本の数学的リテラシーのデータを例に、重み付けの効果を検証する。
具体的には、大見付け前後の平均得点と標準偏差を比較する。
重み付け前後の差は平均得点・標準偏差ともに0.5程度であり、それほど違いは見られない。
#Culculate unighted mean
mean(c(mean(JPN2012$PV1MATH),
mean(JPN2012$PV2MATH),
mean(JPN2012$PV3MATH),
mean(JPN2012$PV4MATH),
mean(JPN2012$PV5MATH)))## [1] 535.9252#Culculate unweighted standard deviation
mean(c(sd(JPN2012$PV1MATH),
sd(JPN2012$PV2MATH),
sd(JPN2012$PV3MATH),
sd(JPN2012$PV4MATH),
sd(JPN2012$PV5MATH)))## [1] 94.05592#Install packages
install.packages("Hmisc", repos = "http://cran.us.r-project.org")## Installing package into 'C:/Users/turid/Documents/R/win-library/4.0'
## (as 'lib' is unspecified)## package 'Hmisc' successfully unpacked and MD5 sums checked## Warning: cannot remove prior installation of package 'Hmisc'## Warning in file.copy(savedcopy, lib, recursive = TRUE): C:
## \Users\turid\Documents\R\win-library\4.0\00LOCK\Hmisc\libs\x64\Hmisc.dll の C:
## \Users\turid\Documents\R\win-library\4.0\Hmisc\libs\x64\Hmisc.dll へのコピーに問
## 題があります: Permission denied## Warning: restored 'Hmisc'##
## The downloaded binary packages are in
## C:\Users\turid\AppData\Local\Temp\RtmpYRNbMQ\downloaded_packageslibrary(Hmisc)## Loading required package: lattice## Loading required package: survival## Loading required package: Formula## Loading required package: ggplot2##
## Attaching package: 'Hmisc'## The following objects are masked from 'package:base':
##
## format.pval, units#Culculate weighted mean
mean(c(wtd.mean(JPN2012$PV1MATH, JPN2012$W_FSTUWT),
wtd.mean(JPN2012$PV2MATH, JPN2012$W_FSTUWT),
wtd.mean(JPN2012$PV3MATH, JPN2012$W_FSTUWT),
wtd.mean(JPN2012$PV4MATH, JPN2012$W_FSTUWT),
wtd.mean(JPN2012$PV5MATH, JPN2012$W_FSTUWT)))## [1] 536.4069#Culculate weighted standard deviation
mean(c(sqrt(wtd.var(JPN2012$PV1MATH, JPN2012$W_FSTUWT)),
sqrt(wtd.var(JPN2012$PV2MATH, JPN2012$W_FSTUWT)),
sqrt(wtd.var(JPN2012$PV3MATH, JPN2012$W_FSTUWT)),
sqrt(wtd.var(JPN2012$PV4MATH, JPN2012$W_FSTUWT)),
sqrt(wtd.var(JPN2012$PV5MATH, JPN2012$W_FSTUWT))))## [1] 93.52401質問紙の場合は重み付けの前後で差が見られるだろうか。
生徒質問紙ST04Q01は、性別を尋ねる質問項目である。
重み付け前後の差は0.1であり、違いは見られない。
#Culculate unweighted propotion of ST04Q01
prop.table(table(JPN2012$ST04Q01))##
## Female Male
## 0.4756731 0.5243269#Culculate weighted propotion of ST04Q01
prop.table(wtd.table(JPN2012$ST04Q01,JPN2012$W_FSTUWT,
type = "table"))## Female Male
## 0.4743144 0.5256856日本のデータは実施率が高く欠損の補正が小さいことや、母集団から生徒が選ばれる抽出確率に大きな違いが無いことから、重み付けの効果が小さくなっている。
生徒の抽出確率が異なっている場合（オーバーサンプリングやアンダーサンプリングが行われている）や，回収率が低い場合は，重み付けの効果が重要になる。