Rのglmとsurveyは、何をしているのか？ 同じデータ・同じ式で比べてみた

禁酒で生活は劇的に良くなるのか。そこに数字で迫ってみたくて、アメリカの公開調査データBRFSSで、飲酒の頻度・量と生活満足度の関係を調べたのが記事1でした。今回はその裏側、技術寄りの話です。

やってみたかったのは、同じデータ・同じ式で、普通のロジスティック回帰（glm）と、調査データ向けのRパッケージsurveyの回帰（svyglm）を比べたら、何が変わるのかということでした。調査の重み・層・PSUを入れると、割合も信頼区間も、どれくらい動くんだろう。

記事1のおさらい

• データ：BRFSS 2023（CDCのBRFSS 2023年版のページと質問票の版ごとのデータのページから取得）。生活満足度の質問を実施した39州とワシントンDCが対象です。
• 対象：すべての分析で共通の255,168人（過去30日に飲酒なし123,182人、飲酒者131,986人）
• 評価項目：生活に「不満」「非常に不満」と答えた人の割合
• モデル：カテゴリ変数を1つだけ入れたロジスティック回帰が2つ。背景は調整しておらず、基準のカテゴリは「過去30日に飲酒なし」です。
  • 頻度モデル（life_dissatisfied ~ frequency_group）：飲酒なし、1〜3日、3日超〜7日、7日超〜15日、15日超〜30日
  • 量モデル（life_dissatisfied ~ quantity_group）：飲酒なし、飲酒日に平均1、2、3〜4、5以上のドリンク

比べる前に、コードをレビューする

記事1と同じく、コードはAIに手伝ってもらいながら書き、中身は自分でレビューしています。2つの方法を比べるなら、比べる土台がそろっていないと話にならないので、レビューではまずそこを見ました。

同じ人たちで比べているか

glmとsurveyで対象者が違っていたら、結果の差が手法のせいなのか、対象のせいなのか分からなくなってしまいます。なので、glmとsurveyが同じ255,168人を使っているかを確認しました。回答者IDも全変数も一致しています。

surveyでは、いったん全303,632件で調査の設計（重み・層・抽出単位）を作ってから、subset()で255,168人に絞っています。分析しない人を先に外すと、調査全体の情報まで落ちてしまう。なので、調査の設計を残したまま、対象者を絞るわけです。

予測は、集計とぴったり一致しているか

カテゴリ変数1つだけのロジスティック回帰なら、予測割合はカテゴリごとの集計と同じになるはずです。glmなら「不満と答えた人数÷人数」、surveyなら「不満と答えた人の重みの合計÷重みの合計」。突き合わせてみると、両モデル・両手法で、ほぼ誤差なく一致していました。

surveyの結果は、重み付きの集計や、surveyの関数を使わずに別に計算した分散などと突き合わせるチェックも、58項目走らせています。チェックのコードもAIに手伝ってもらったものですが、少なくとも数字が設計どおりに計算できていることは、これで確認できました。

surveyは何をしているのか

BRFSSの大きなデータより、小さな架空の例で考えるほうが、重みの役割は見えやすそうです。

ここで使う数字は説明のための架空の例で、BRFSSの値ではありません。

人口が同じA・Bの2つの地域があって、Aから80人、Bから20人の回答が集まったとします。不満ありは、Aで8人、Bで6人。

全員を同じ1人として数えると、こうです。

(8 + 6) / (80 + 20) = 14%

ただ、2つの地域の人口は同じでした。ということは、Bの回答者1人は、Aの回答者1人の4倍の人数を代表していることになります。そこでA=1、B=4という重みを付けて、人口の構成に合わせて数え直すと……

(8×1 + 6×4) / (80×1 + 20×4) = 32 / 160 = 20%

14%が20%になりました。重みって、信頼区間だけの話じゃなくて、割合そのものを動かすんですね。

重みは「ぶれ」にも効いてくる

とはいえ、Bの20人に重み4を付けたからといって、Bで80人に聞いたことにはなりません。重みを付けると、回答が1人変わったときの割合の動き方まで変わります。Aの1人が変われば1/160で0.625パーセントポイント（pp。割合同士の引き算の単位）、Bの1人なら4/160で2.5pp。1人の回答の影響がこれだけ違うので、推定の不確かさの計算にも重みを反映させます。

層は、調査のときに分けて抽出した区分のことです。Aから80人、Bから20人を別々に抽出したなら、その分け方も計算に入れます。「もう一度同じ調査をしたら、答えはどれくらいぶれるだろう」と考えるときも、AとBを分けた調査として計算するわけです。この区分は調査を設計する側が決めたもので、飲酒量のカテゴリや、後から作った年代区分とは役割が違います。

surveyがやっていること

今回surveyでやっていることを、大きく分けるとこうなります。

1. 重みで、各回答をどれくらい数えるかを変える（→ 割合が変わる）
2. 重み・層・抽出単位（PSU）から、推定値の不確かさを計算する（→ 信頼区間が変わる）

今回は信頼区間が広がりました。調査設計によって幅の変わり方は異なり、層化によって推定のぶれが小さくなる場合もあります。

それと、最初に勘違いしそうになったのがfamily = quasibinomial()です。これは、重み付きの0/1データでbinomialを使うと出る警告を避けるための指定でした。重み・層・PSUを渡しているのは、designのほうです。

BRFSSでの指定

BRFSSには、重み・層・PSUの情報があらかじめ入っています。今回は次のように指定しました。

• 重み：CDCが公開している、主ファイルと質問票の版ごとの重み。生活満足度を複数の版で聞いた州では、版ごとの人数の比で重みを調整しています。
• 層：州と、州の中の層（_STSTR）の組み合わせ。
• PSU：州と_PSUの組み合わせ。このデータでは、PSUは回答者ごとに1つずつでした。

library(survey)
options(survey.lonely.psu = "adjust", survey.adjust.domain.lonely = TRUE)

# d は生活満足度モジュールの全303,632件。ここから設計を作る
# （stratum_id は州＋_STSTR、psu_id は州＋_PSU、weight は調整済みの重み）
full <- svydesign(ids = ~psu_id, strata = ~stratum_id, weights = ~weight,
                  data = d, nest = TRUE)

# glm と同じ255,168人に絞る
design <- subset(full, analysis_eligible)

# 同じ式を glm と svyglm に渡す（S は glm 用の255,168人のデータ）
fit_glm <- glm(life_dissatisfied ~ quantity_group, data = S, family = binomial())
fit_svy <- svyglm(life_dissatisfied ~ quantity_group, design = design,
                  family = quasibinomial())

最初のoptions()の1行は、層の中に抽出単位が1つしかないときの扱いを決める設定です。

割合と、その「幅」

この記事では、カテゴリごとの不満割合と、飲まない人のカテゴリとの差に、それぞれ95%信頼区間を付けています。ざっくり言うと、「本当の値はだいたいこの範囲に入っていそう」という幅です。

この幅は、AIに書いてもらったコードで出しています。重み付きの集計や、別のやり方で求めた分散と突き合わせて、入力したものに対して何が出てくるかを確かめました。計算式そのものの導出は、ここでは追っていません。

glmの結果もsvyglmの結果も、区間を計算する同じ関数にそのまま渡せます。使い方は同じでも、中の「ばらつきの見積もり」（Rでいうとvcov()の中身）が違う。glmは全員を同じ1人として数え、svyglmは調査の重みや層を反映して見積もっています。

差のほうは、その幅に0が入っているかどうかを見ています。0が入っていなければ「飲まない人とは差がありそう」、入っていれば「差があるとは言い切れない」と読んでいます。

glmとsurveyを並べてみる

ここからが、glmとsurveyを並べた結果です。図の左が各カテゴリの不満割合、右が「過去30日に飲酒なし」との差。点が推定値、横線が95%信頼区間で、青の●が通常のglm、橙の▲がsurveyです。右の図の縦の点線が差0の位置です。

頻度

カテゴリ	glm 不満割合 %	survey 不満割合 %	glm 飲酒なしとの差 pp	survey 飲酒なしとの差 pp
過去30日に飲酒なし	6.3（6.2〜6.4）	6.7（6.4〜7.0）	基準	基準
1〜3日	5.2（5.0〜5.4）	6.0（5.6〜6.5）	−1.1（−1.4〜−0.9）	−0.7（−1.3〜−0.2）
3日超〜7日	4.2（4.0〜4.4）	5.0（4.5〜5.6）	−2.1（−2.4〜−1.9）	−1.7（−2.3〜−1.1）
7日超〜15日	3.9（3.6〜4.1）	4.5（4.0〜5.0）	−2.4（−2.7〜−2.2）	−2.2（−2.8〜−1.6）
15日超〜30日	4.7（4.5〜5.0）	5.9（5.3〜6.7）	−1.6（−1.9〜−1.3）	−0.8（−1.6〜−0.0）

「15日超〜30日」のsurveyの差は、表だと幅の上の端が「−0.0」になっていて紛らわしいのですが、丸める前は−0.029ppで、ぎりぎり0より下です。なので、0は入っていません。

量

カテゴリ	glm 不満割合 %	survey 不満割合 %	glm 飲酒なしとの差 pp	survey 飲酒なしとの差 pp
過去30日に飲酒なし	6.3（6.2〜6.4）	6.7（6.4〜7.0）	基準	基準
飲酒日に平均1ドリンク	3.7（3.5〜3.8）	4.4（4.1〜4.8）	−2.6（−2.8〜−2.4）	−2.3（−2.8〜−1.9）
平均2ドリンク	4.1（3.9〜4.3）	5.0（4.5〜5.5）	−2.2（−2.4〜−2.0）	−1.8（−2.3〜−1.2）
平均3〜4ドリンク	6.1（5.8〜6.4）	6.7（6.1〜7.4）	−0.2（−0.6〜＋0.1）	＋0.0（−0.7〜＋0.7）
平均5ドリンク以上	8.6（8.1〜9.2）	9.0（7.9〜10.3）	＋2.3（＋1.7〜＋2.9）	＋2.3（＋1.0〜＋3.6）

何が変わって、何が変わらなかったか

いちばん分かりやすい変化は、割合そのものです。同じ回答者・同じ式なのに、surveyでは全カテゴリで不満割合が0.4〜1.2pp高くなっています。重みで数え方が変わった結果です。

信頼区間のほうはもっと大きく動いていて、不満割合の95%信頼区間の幅はglmの2.09〜2.83倍。人数は同じなのに、ここまで変わるとは思っていませんでした。重みのばらつきが効いているんだろうな、とは想像していますが、どの要因がどれだけ効いたのかまでは、今回は分解していません。

変わらなかったのは、飲まない人のカテゴリとの差の判定です。95%区間が0を含むかどうかは、頻度・量のどのカテゴリでも同じでした。

じゃあ差の結果は同じかというと、そうでもありません。たとえば平均5ドリンク以上。差は約＋2.3pp（surveyでは6.7%と9.0%の差）でほぼ動いていないのに、95%区間は＋1.7〜＋2.9ppから＋1.0〜＋3.6ppへ広がっています。「有意かどうか」の二択で見ていたら、どちらも「有意」で終わってしまって、この違いには気づけなかったと思います。

数え方が変わると、割合も区間も変わる

同じ回答者・同じ式でも、何をどう数えるかで、割合も信頼区間も変わる。記事1では「何を比べた数字なのかを見たい」と書きましたが、手法を比べるときも、気にするところは同じでした。

こういう解析のいいところは、数字の出どころをたどれることだと思います。割合ひとつ、区間ひとつでも、どういう計算でそうなったのかを追いかけられる。今回も、重みや調査設計を入れると何が変わるのかを、集計とコードを見比べながら確かめられました。

実務でこの規模の調査データを扱う機会は、多くないかもしれません。今回は日本で目的に合う申請不要の個票データが見つからず、アメリカのデータを使いました。使えるデータは限られますが、調査設計を反映した推定を手元で確かめられる題材でした。

コードとデータ

• 分析のRコードと図：GitHubのsasai-lab/qiita-article-assetsのalcohol-life-satisfaction-surveyフォルダに置いています（glmとsurveyの推定と検証のコードです）。
• データ：CDCのBRFSS 2023年版のページと、質問票の版ごとのデータのページから取得しました。
• 環境：R 4.6.1、survey 4.5
• Rの環境構築について：LLM時代のR環境構築：コマンドは覚えなくていい