SPSSのウェイトの扱いについて調べていたら、「SPSSとかSASとか！そういう統計ソフトの！ ウェイトとか重み付けとかっていうのは！ 調査でいうところの！ いわゆるウェイトのことではないです！」と書いているサイトを見つけた。「SPSSでいう「ケースの重み付け」は確率ウェイティングのことではない。」という話は知っていたが、自分では検証していなかったのでありがたい。
それはともかく、Rのsurvey packageなら正しく標準誤差を算出できるのだろうか。気になったので、このサイトの例題をもとに確認してみた。
library(survey)
d <- data.frame(
ans = c(1, 2, 2, 3, 3, 4, 4, 1, 2, 3),
gender = c(rep(1, 7), rep(2, 3)),
w = c(rep(5 / 7, 7), rep(5 / 3, 3))
)
d$prob <- 1 / d$w # 抽出確率を計算（なんとなくやってみただけ）
# des1でもdes2でも同じこと
des1 <- svydesign(ids = ~1, strata = ~gender, weights = ~w, data = d)
des2 <- svydesign(ids = ~1, strata = ~gender, prob = ~prob, data = d)
svymean(~ans, des1)
svyby(~ans, by = ~gender, des1, svymean)
svymean(~ans, des2)
svyby(~ans, by = ~gender, des2, svymean)
どうやら「正解」を導けたようだ。ついでなので、ジャックナイフ法による標準誤差の算出も試してみた。以下では、あえてsvyrepdesignに落とし込んで計算している。
jkw <- as.svrepdesign(des1, type = "JKn")
# replicate weightsは下記コマンドで取り出せる
# weights(jkw, type = "analysis")
# https://stackoverflow.com/questions/73627746/how-do-i-get-the-weights-from-a-survey-design-object-in-r
d2 <- cbind(d, data.frame(weights(jkw, type = "analysis")))
# rscalesについては下記資料（p.13, p.28）を参照
# https://faculty.washington.edu/tlumley/survey-jsm-nup.pdf
des3 <- svrepdesign(
weights = ~w, repweights = "X[0-9]+", type = "JKn",
data = d2, rscales = jkw$rscales
)
svymean(~ans, des3)
svyby(~ans, by = ~gender, des3, svymean)
# replicate weightsを取り出さなくても可
svymean(~ans, jkw)
svyby(~ans, by = ~gender, jkw, svymean)
当たり前だが計算結果は一致する。いろいろやって気づいたのだが、少し設定をミスっただけで標準誤差がおかしくなる。survey dataの分析は本当に難しい。
