日々の記録簿

日々のできごとの記録

日記

survey package(R)の使い方

投稿日:





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の分析は本当に難しい。

スポンサーリンク

スポンサーリンク

-日記

執筆者:


comment

メールアドレスが公開されることはありません。 が付いている欄は必須項目です

関連記事

全国学力・学習状況調査の都道府県別平均正答率csvファイル

ごく稀に必要になることがあるのですが、どこにあったっけ・・・と探すのがめんどくさいので、アップロードしてしまうことにしました。 全国学力・学習状況調査の都道府県別平均正答率csvファイルです。ただし、 …

全国学力・学習状況調査を悉皆で実施することの問題点

全国学力・学習状況調査(通称、全国学力テスト)は、対象となる全員を調査する調査方法(いわゆる悉皆調査)で実施されていますが、その問題点について整理しておきます。 悉皆調査の問題点はいろいろありますが、 …

beamerでQRコードを表示する

タイトルのとおりだが、beamerでプレゼンテーションをするときにQRコードを埋め込みたくなった。できるだけ簡単な方法を探ったところ、このサイトがヒット。 要は{textpos}を使って、任意の位置に …

棒グラフとヒストグラムの違い

棒グラフとヒストグラムの違いがわからないというコメントを受けることが多いので、どう説明したものかと思っていたのですが、次のような図を書いたら、理解してくれる人が多いような気がしたのでご紹介。 まずデー …

x1 carbon (6th gen)にLinuxをインストールしてWWANを使う

Thinkpad X1 carbon (6th gen)のWWAN(Fibocom L850-GL)はLinuxでは動作しない。しかし、最近この問題がクリアされたようである。例によってソースはArch …