前回SPSSの「ケースの重み付け」を調べて満足していたのだが、Rのコードを見ていてとんでもないことに気づいた。
weighted.se <- function(x, w, na.rm = FALSE) {
if (na.rm) {
w <- w[!is.na(x)]
x <- x[!is.na(x)]
}
v <- sum(w * (x - weighted.mean(x, w))^2) / (sum(w) - 1)
sqrt(v) / sqrt(sum(w))
}
標準誤差を算出する際、最後にウェイトの合計の平方根sqrt（sum(w))で割っているのである。大規模調査ではウェイトの合計が母集団の人数になっている。たとえばPISA2018の日本のデータだとウェイト（W_FSTUWT）の合計は1078921だ。前回の記事で使った例題ではウェイトの合計は10でサンプルサイズと一致していたが、実際の大規模データではウェイトの合計が1000になっていてもおかしくない。この場合、SPSSで標準誤差を算出すると「「正解」よりも少し小さい」では済まない値になる。ここでは単純に前回のウェイト（w1）を100倍したウェイト（w2）を使って標準誤差を算出してみる。
library(survey)
weighted.se <- function(x, w, na.rm = FALSE) {
if (na.rm) {
w <- w[!is.na(x)]
x <- x[!is.na(x)]
}
v <- sum(w * (x - weighted.mean(x, w))^2) / (sum(w) - 1)
sqrt(v) / sqrt(sum(w))
}
d <- data.frame(
ans = c(1, 2, 2, 3, 3, 4, 4, 1, 2, 3),
gender = c(rep(1, 7), rep(2, 3)),
w1 = c(rep(5 / 7, 7), rep(5 / 3, 3)),
w2 = c(rep(5 / 7, 7), rep(5 / 3, 3)) * 100
)
weighted.se(d$ans, d$w1) # 0.3319131
weighted.se(d$ans, d$w2) # 0.03150379
# survey package
des1 <- svydesign(ids = ~1, strata = ~gender, weights = ~w1, data = d)
des2 <- svydesign(ids = ~1, strata = ~gender, weights = ~w2, data = d)
svymean(~ans, des1) # 0.3571
svymean(~ans, des2) # 0.3571
# 重み付けなしの標準誤差
std.err <- function(x) sd(x, na.rm = TRUE) / sqrt(length(x[!is.na(x)]))
mean(d$ans) # 2.5
std.err(d$ans) # 0.341565
「正解」は0.357に対し、w1は0.332、w2は0.032（！）である。念の為にSPSSで確認してみたが、Rと同じ結果になった。これでは判断を誤る。幸い（？）推定値は正しいので、それぞれの大規模調査で推奨されているreplication methodsを使えば「正解」を導くことは可能だと思う。あるいは思い切って「ケースの重み付け」をしない（平均は2.5、標準誤差は0.342になる）というのも一つの手なのかもしれない。ウェイトの扱いはくれぐれも慎重に。
