TIMSS報告書の標準誤差（PVsを使った場合）は、正規の式による計算ではなく、ショートカット計算になっています。
TIMSSを使った分析をしていて、報告書と数値が合わない・・・となった方の参考になるかもしれません。
TIMSS報告書の標準誤差の計算
PVsを使うとき、最終的な分散は、下記の式で計算される。
最終的な分散は、σ(μ)とσ(test)の和である。ここでσ(μ)とσ(test)は、それぞれ次のようになる。
まず、σ(test)は、PV1からPV5の推定値の不偏分散である。
Rで計算すると以下のようになる。ここでは、TIMSS2011の日本のデータを用い、小4の科学得点（ASSSCI）の平均値を計算する。
library(foreign)
jp11 <- read.spss("asgjpnm5.sav", use.value.labels=F, to.data.frame=T, use.missings=F)
pv <- paste("ASSSCI0", 1:5, sep="")
se <- sapply(pv, function(i){weighted.mean(jp11[[i]], jp11$TOTWGT)})
var(se) # σ(test)
だいたい、1.35程度の数値になる。続いて、σ(μ)を考える。
σ(μ)は、PV1からPV5の分散の平均である。
Rで計算すると以下のようになる。
# Replicate Weightを作成し、75回分の平均を計算する関数を作成
jk <- function(x){sapply(1:75, function(i){
weighted.mean(x,
ifelse(jp11[["JKZONE"]]==i,
2 * jp11[["TOTWGT"]] * jp11[["JKREP"]],
jp11[["TOTWGT"]])
)
}
)
}
rw <- sapply(pv, function(i){jk(jp11[[i]])}) # 75×5回分の計算
v <- mean(sapply(1:5, function(i){sum((rw[,i] - se[i]) ^ 2)})) # σ(μ)
sqrt(v + 1.2 * var(se)) # 正式な値
ところが、この計算では、TIMSSの報告書と値が合わない。
TIMSS報告書では、ショートカットが使用されており、σ(μ)の代用として、PV1のみが使用されているのである。
Rのコードを書くと、次のようになっている。
sum((rw[,1] - se[1]) ^ 2) # これがσ(μ)の代用にされている sqrt(sum((rw[,1] - se[1]) ^ 2) + 1.2 * var(se)) # 最終的な標準誤差
このショートカット計算は、PISA Data Analysis Manualのp.129にも記載されている。
おそらく、TIMSSが開始された1995年はコンピュータの性能が低く、Replicate Weightsを使ったり、PVsを使ったりすると、計算に時間がかかりすぎるためにショートカットが利用されていたのだろう。
大抵の場合、ショートカットで計算した値は、PV1からPV5を使用したときの標準誤差と、ほとんど変わらないため、大きな問題にはならない。ただ、報告書の値を再現しようとしたときに、このことに気づかないと延々と悩むことになる。
ちなみに、intsvy packageを使った場合は、ショートカット計算の値が適用されるので、報告書の値と一致するはずである。
