集団と個人の階層データに対するマルチレベルモデリングの中心化は厄介
マルチレベル分析を必要とする研究仮説を扱おうとして，間違った解析結果を載せそうになった話．RQと推定したい指標は対応してるか？は基本の問だと思うんだけど，なかなか難しいなあ．
あまり細かく解説していないので，「おや？？？」と思った人は本や論文を読むことをおすすめします．気が向いたら，ちゃんとモデル式を書いた記事を作るかも．
結論
マルチレベルモデルは，自分が何を推定したいのか，モデルが何を意味しているのかを理解した上で使わないと，結構危ないものだと思います．自分も確認し直したからよかったけど，本当に間違えかけた．勉強は大事ですね．
まずはレベルごとの回帰式を書きましょう（マルチレベルモデリングの基本）．その上で，勉強しましょう．
僕が参考にしたものを提示します．
このブログを見るのも良いかもしれませんが，論文と本を読みましょう（暴論）
参考になるもの
- Enders & Tofighi, 2007は必読
-
Raudenbush, S. W., & Bryk, A. S. (2002). Hierarchical linear models: Applications and data analysis methods (Vol. 1). sage.
分厚い本ですけど，中心化の説明はとてもわかりやすいと感じました．
直後に紹介する，村山さんの資料はこの本の内容を汲んでいます． - 日本語の書籍なら，この本．第6章がとても素晴らしいと感じた．
日本語でググると出てくるもの（これらもかなりわかりやすいと思います）
僕が作ったshinyappも参考にどうぞ
論文や本だけだと，イマイチ飲み込めない感覚があると思います．具体的に図を作りながら，分析結果と合わせて確認できるwebアプリを作りました．ぜひご参考に．
扱うRQ
集団の文脈効果 (contextual effect) を求める研究仮説を扱います．
文脈効果とは？
Roux, 2002によると，contextual effectは次のように説明されています．
Term generally used to refer to the effects of variables defined at a higher level (usually at the group level) on outcomes defined at a lower level (usually at the individual level) after controlling for relevant individual level (lower level) confounders
個人と集団という階層データを考えると，個人レベルに依存しない集団レベルの影響と言えます．簡単な例だと，地域に文脈効果がある状況で，ある人が地域Aから地域Bに引っ越した場合の影響と考えれば良さそうです．
今回考える研究の設定
例として，職場と労働者という階層性があるデータを用いて，労働者の残業時間と職務満足度との関連を考えましょう．
個人単位で物事を考えると，残業時間が伸びると，職務満足度が下がりそうだなと感じます．一方，好きなお仕事を長く楽しくできて，満足する人もいるかもしれません．
職場単位で考えてみると，みんなが残業しまくる残業体質な職場だと，不満が溜まりそうですね．
このように，「労働者の残業時間と職務満足度」という中に，個人レベルの影響と職場レベルの影響の2つが存在することになります．そこでマルチレベル分析を採用するわけです（本当にざっくり説明です．観測値の独立性の問題とかも押さえようね）．
分析で扱う残業時間は次の処理を行うとします．
- 個人レベルで月間残業時間を評価．
- 集団レベルの残業時間の評価では，職場ごとに残業時間の平均値を使用．
つまり，職場Aの労働者の平均残業時間が職場Aの指標になります．
調べたいことは，「残業体質な職場で働く労働者は，個人の残業時間に依存せず，職務満足度が下がるか？」とします．よって，職場レベルの残業時間の文脈効果を推定していきます．
中心化
今回の研究トピックのように，次の状況では，"変数 - 変数の平均値"という中心化の処理が施されることが多いと感じます．
- 連続値で評価された個人レベルの説明変数
- 集団レベルの説明変数として，個人レベルを集計した値を使用
- マルチレベルモデリング
ここで，マルチレベルモデリングの場合は「平均」を複数計算できることに注意します．今回のように，個人と集団というレベルが2つ存在する場合は，集団ごとの平均値と，標本全体の平均値の2通りが存在します．
そこで，中心化に関して3種類の分析を実行して，結果を比較します．
- rawの値を投入する
- 個人レベルの残業時間: そのまま投入
- 職場レベルの残業時間: 各職場の個人の平均値をそのまま投入
- 集団平均中心化 (Centering Within Cluster, CWC)
- 個人レベル: 個人レベルの残業時間 - 職場の平均残業時間
- 職場レベル: 職場レベルの平均残業時間 - 全体の平均残業時間
- 全体平均中心化 (Centering at Grand Mean, CGM)
- 個人レベル: 個人レベルの残業時間 - 全体の平均残業時間
- 職場レベル: 職場レベルの平均残業時間 - 全体の平均残業時間
結論としては，contextual effectは次の値となります．
- rawの場合: 職場レベルの回帰係数
- CWCの場合: 職場レベルの回帰係数 - 個人レベルの回帰係数
- CGMの場合: 職場レベルの回帰係数
中心化によって，推定値が異なる！？では実際に，確認しましょう．
Shinyによる確認
こちらでお試しできます．
Rの知識は特に必要としていませんが，正しく理解しようとするなら，{lme4}パッケージとlmer関数の出力の理解が必要です．
データの見方
左側のパネルで，職場の数と，各職場の人数を変更でき，結果のシナリオを変更することができます．
まずは，データのプロットが表示されるはずです．ここでは横軸に残業時間を，縦軸に職務満足度を取っています．そして，各色が職場を表し，回帰直線は職場ごとに引かれています．
例えば，最初のシナリオである「1) 個人の残業時間が長いと満足度が上がるが，職場平均残業時間が長いと個人の満足度は下がる」の場合，各回帰直線は右肩上がりですが，残業時間が長い職場だと，満足度が下がっていくことがわかります．個人と集団レベルでは異なる関連が認められるケースですね．
解析結果
結論としては，Fixed effectを見てください．
上から順に切片，個人レベルの変数，職場レベルの変数，と並んでいます．
確認してほしいこと
まず，文脈効果 (contextual effect)は次の値となります．
- rawの場合: 職場レベルの回帰係数
- CWCの場合: 職場レベルの回帰係数 - 個人レベルの回帰係数
- CGMの場合: 職場レベルの回帰係数
その上で，以下を確認してください．
お伝えしたいこと
文脈効果を求めていないのに，「文脈効果を求めた」と書いている論文があります．僕もそう書きそうになりました．とても怖い話です．
自分の求めたestimandは何を表しているのか，理論に基づいて考え直すことが重要だなあと改めて実感しました．マルチレベルモデルは，自分が何を推定したいのか，モデルが何を意味しているのかを理解した上で使うべきものだと思います．
参考になるもの
- Enders & Tofighi, 2007は必読
-
Raudenbush, S. W., & Bryk, A. S. (2002). Hierarchical linear models: Applications and data analysis methods (Vol. 1). sage.
分厚い本ですけど，中心化の説明はとてもわかりやすいと感じました．
直後に紹介する，村山さんの資料はこの本の内容を汲んでいます． - 日本語の書籍なら，この本．第6章がとても素晴らしいと感じた．
日本語でググると出てくるもの（これらもかなりわかりやすいと思います）
「良い質問」は，良い研究・論文を考えることから始まると思う
所属している専攻の修士論文発表会では，質疑応答の時間に学生が学生に対して質問することが推奨されています．そして，教員が良い質問をした学生を評価する質問賞が設けられています．「あなた良い質問したね！」を褒めてもらえるわけです．良い試みですよね，楽しかったです．
私は質問することは結構大事だと思います．研究などの発表に質問をすることで，発表のさらなる改善を助けたり，自分の理解を深めることができると感じるからです．また，良い質問とは，その効果が大きいことを指すのでしょう．では，「良い質問」をするためには，何をどう考えれば良いのでしょう？今回は，個人的な意見を整理します．
良い悪いの前に質問自体が素晴らしい
良い質問/悪い質問の前に，基本的に他人の研究に対して質問をすること自体が素晴らしいことだと思います．もちろん攻撃的な質問（たとえば，「この研究に何の意味があるんですか？まったく体をなしていないと思います」）は，相手を傷つけるだけで何も生まないので，このような質問をするべきではありません．こうした攻撃的・否定的な質問を除けば，質問するということ自体が良いことだと感じます．
したがって，質問の良し悪しを考える前に，「〇〇がわからなかったので，教えて下さい」を言うだけでそもそも素晴らしい，ということが自分の意見です．1番シンプルな質問で全然okだと思います．悪い質問とか考えず，まずは素直に質問しようよ，と思います．
とはいえ，ゼミや発表会に参加していると，やっぱり良い質問ってありますよね．建設的な意見とか，ズバッと大事なことを質問している，とか．もう少し考えてみます．
"理想の研究・論文と目の前の発表内容におけるギャップ"を見つける
自分が理想と考える研究や論文の情報と，今聞いた発表との差を見つけることが，良い質問のスタート地点だと思います．理想と提示された研究発表とのギャップに基づいて「ここの情報 or 論拠 or 考察が足りなかったように思うのですが，どうお考えですか？」と聞けば，批判的に切り込んだ質問．ギャップを埋めるために「現在の発表だと，〇〇は足りなかったと感じました．□□のような内容を含めるともっと良くなると思いますが，ご意見いただけますか？」みたいに提示すると建設的．いずれもスタート地点は同じだと思います．
理想の研究って？
これは本当に個人の感想・考えです．領域や信念に従って自分なりの理想を考えればいいと思います．簡潔にまとめるなら，①新規性・有用性・尤もらしさの主張，②RQと手法の一致，③introductionとdiscussionが対応していること，の3点だと思っています．
Introduction
- big pictureを描けている
- 研究トピックの新規性・有用性・尤もらしさを提示できる
- 新規性: すでにわかりきったことを研究として発表されても困る．先行研究で何がわかっていて，何がわかっていないかを提示する．
- 有用性: 役に立たないであろう話をされても困る．これがわかるとどう嬉しいか，提示する．
- 尤もらしさ: その研究トピックがなぜ考えられるのか？突拍子もない話をされても困る．どうしてその仮説が尤もらしいのか提示する．可能なら理論に基づいて．
（descriptiveに提示するような研究であれば，ここは不要なことも．疾患の有病率を提示するような研究で「尤もらしさ」はない）
- big pictureからRQに向かって，話が一直線に絞られていく．
- 研究は全体の構成として砂時計型 (hourglass format) であるべき．
introductionでは，big pictureからRQに向かって，話題が単調に絞り込まれていくことになる． - 途中で「ところで」とか「あ，全然関係ないんやけど」といった割り込みはやめてほしい．
- 研究は全体の構成として砂時計型 (hourglass format) であるべき．
Method
- RQと手法が一致する（個人的に，methodはこの観点でしか聞いていない）
- 研究参加者の対象
- 変数の評価方法
- 解析手法
Results
（ここはあんまり文句言うことないかなあ．あるがままに記述して，きれいにtableとかfigureを作ればokだと思う）
Discussion
- methodとresultsに基づいて慎重に主張する．言えないことは言わない．
- introductionの伏線をちゃんと回収している．
- 結果の解釈
- 今回の結果と先行研究との一致・不一致．可能であれば効果・関連の大小．
「これまでに〇〇がわかっています」という新規性を主張する部分で引用した論文が登場するイメージ． - どうしてそういう結果になったのか，考察．
「この研究トピックは△△ということから考えられます」と仮説の尤もらしさをintroductionで提示しているので，結果を踏まえてさらに考察する．
- 今回の結果と先行研究との一致・不一致．可能であれば効果・関連の大小．
- 有用性の主張
- 結果がわかった上で，どう役立ちそうか主張する．
「このRQがわかると役に立つ」と有用性をintroductionで主張しているのだから，結果がわかった上でどうなのかを提示． - 将来の研究にとって役立つこと，実践にとって役立つこと，等領域ごとに書く内容は変わりうることに注意．
- 結果がわかった上で，どう役立ちそうか主張する．
- 結果の解釈
- 限界を整理する．可能ならjustifyする．
- 正当化していると面白い．
「XXの限界がある．とはいえ，XXによってこう結果が歪みうるので，今回の結論には問題がないだろう」みたいな論文大好き．
- 正当化していると面白い．
※自分の理想を整理すると，「気に食わねえなあ，この研究」と思ったときのモヤモヤ感を言語化しやすいのでオススメ
質問をまろやかにする
"自分の理想 - 目の前にある研究" の引き算を実行して，質問したいことができたとします．自分なりに理想と現状のギャップを埋めるためのアイデアがあれば，それを伝えましょう．建設的な質問・議論になって，とても有意義だと思います．
ギャップを埋めるためのアイデアは無いが，重要であろう質問を思いついたとします．これは少し危険です．そのまま質問をぶつけると，めちゃくちゃ怖い質問になりがちです．火の玉ストレートというやつです．
表現をマイルドにするテンプレートをいくつか用意しておきましょう．自分の中の一例を紹介します．
質問の導入
- とりあえず褒める
- 怖そうな態度を取らない
- 「わからなかったので教えていただきたい」
- 「自分の理解が追いついているか曖昧なので，確認させてください」
- 「興味のある分野でぜひ勉強したいから，教えて下さい」
質問の内容
基本具体的にしてあげると，急にマイルドになると思っています．（具体例を自分で作らなくちゃいけないのでやや面倒）
- NG: 〇〇バイアスがあると限界で言っていますが，どういうメカニズムでそのバイアスが起こりますか？
- → 〇〇バイアスについて，例えばXXXがあると思いますけど，その他に考えているものがあれば教えて下さい
- NG②: この手法で評価しているものは，正しく測りたいものを評価できているのでしょうか？
- ここで評価しているYYYは，やや限界が残ると思います．例えばαとβの問題があるのではないかと感じました．その点についてご意見伺ってもよろしいでしょうか？
その他，自分の中でありがちな質問をいくつか．なお，全然NGと実際の質問が一致していない気がするのですが，ふんわり質問から入ると，本題に切り込んだ場合の恐怖感を減らせているんじゃないかと感じます．
- 尤もらしさの質問のNG例: なぜこの仮説が導かれるのですか？
- この研究のアイデアを思いついた経緯を教えてください
- 有用性の質問のNG例: この結果は，どう役に立つんですか？
- この結果を踏まえて，どう活かしたいですか？
- この結果や考察を，誰に1番届けたいですか？どういうアクションを期待しますか？
- 限界に対する質問のNG例: 〇〇を限界としていましたが，どういう意味で限界なんですか？これによって，結果がどう歪みうるのですか？
- もしデータセットや資金・時間などが理想的な状態で，同じ研究仮説を調べられるとすると，どうしたいですか？
- RQと手法の不一致のNG例: 仮説として〇〇を提示しているのに，手法がうまく対応していないので，この仮説を正しく調べられていないと思います．なぜこの手法を選んだのですか？
- NG例②のように，できるだけ具体的に．
まとめ
- 「良い質問」とは，理想の研究と現状とのギャップを考えることから始まると思う．
- できるだけ答えやすくしてあげよう
- 具体的に指摘
- ふんわり質問を用意しておく
普段から「良い質問だ！」と思ったら蓄積する，という習慣をつけていると上手くなれるかもしれませんね．
稀でない2値イベントの回帰モデリング
前回，「"稀な疾病"とはどの程度稀なのか？」という記事を書きました．
今回は稀でないイベントに対する回帰モデルの簡単な紹介をします．2値アウトカムに対する回帰モデリングではロジスティック回帰がよく使われていますね．イベント発生割合が10%以下であれば，ロジスティック回帰を適用することに基本的に問題はないと思います．しかし稀でないイベントでオッズ比を求めてしまうとリスク比との乖離が大きく，効果の指標として不適切かもしれない，という話を前回やりました．
今回は，イベント発生割合が10%を超え，rare disease assumptionを仮定できない場合の回帰モデルとして，修正ポアソン回帰 (modified poisson regression) を紹介します．あくまで導入的な内容であることにご注意を．
目次
データ: タイタニック号乗客者の生存状況
有名なデータセットである，Titanicを使います．欠損値の話は本題ではないので，今回はcomplete case analysisとします．
調べる仮説としては「タイタニック号では，女性に比べて男性で死亡というイベントが起こりやすかった」としましょうか．次の変数を用いて，この仮説を調べます．
- アウトカム: 生存/死亡．死亡のリスクを評価する．
- 説明変数: 性別
- 共変量: 子ども (age < 20) か，大人か (age ≥ 20)
Rだと{carData}をはじめ，様々なパッケージにTitanicSurvivalデータが入っています．今回はそちらを利用します．
記述統計
Table1を提示します．
今回のデータはイベント (dead) が60%近く，稀でないイベントとなっていますね．
また年齢で欠損が263件あります．今回，欠測値の扱いが本題ではないので，リストワイズ除去（欠測値をすべて除外する）によるcomplete case analysisをしてしまいます（本当は良くないですよ）．
2×2表
リストワイズ除去を行ったn=1,046のデータに対し，年齢が20歳未満を「子ども」20歳以上を「大人」とカテゴリ化して，2×2表を作ります．
(RR: 女性に対する男性のリスク比; OR: 女性に対する男性のオッズ比)
子ども: RR=2.39, OR=6.01
| 死亡 (人) | 生存 (人) | 合計 (人) | |
|---|---|---|---|
| 男性 | 88 | 34 | 122 |
| 女性 | 31 | 72 | 103 |
大人: RR=3.56, OR=14.58
| 死亡 (人) | 生存 (人) | 合計 (人) | |
|---|---|---|---|
| 男性 | 435 | 101 | 536 |
| 女性 | 65 | 220 | 285 |
Overall:
- RR=3.21 [95%CI: 2.69–3.84], OR=11.78 [CI: 8.74–15.88] (crude)
- adjusted-RR=3.23 [CI: 2.70–3.87], adjusted-OR=11.41 [CI: 8.47–15.38] (Mantel-Haenszel)
| 死亡 (人) | 生存 (人) | 合計 (人) | |
|---|---|---|---|
| 男性 | 523 | 135 | 658 |
| 女性 | 96 | 292 | 388 |
稀でないイベントを評価しているので，リスク比とオッズ比の乖離がとんでもないことになっています…！こうしたデータで「ロジスティック回帰を用いてオッズ比を...」とするのは要注意！，という話ですね．
回帰分析へ
先程の2×2表による分析では，Mantel-Haenszel法による層別解析を行いましたが，層別解析ではなく回帰分析を行いたい場面は多々あると思います．例えば連続的な共変量を調整したいときなどです．今回はRRを導く回帰分析として，修正ポアソン回帰を紹介します．
修正ポアソン回帰
概要
一般に，ポアソン回帰は0以上の整数値を取るデータに対して適用される回帰モデルで，特に「稀なイベント」に使われています．例えば，交差点での交通事故発生件数や製品製造ラインでの故障品などが代表的でしょうか．
なお，故障品のモデリングなど「全体でいくつ作っているか」も重要になる場合は，オフセット項を導入することで，割合もモデリングすることが可能になります．
[FYI: J-STAGEの解説記事]
しかし，今回の状況では，アウトカムが「0以上の整数値」ではなく「0/1の2値」です．少しイメージしにくいので，Petersonら (2008) の表現を参考にしてみましょう．
It is well known that when the prevalence is low and the sample size is large, probabilities from the Poisson distribution can often be used to approximate probabilities from the binomial distribution. Similarly, one can think of an existing sample of binomial data (0 or 1) as being approximately Poisson, where the probability of a value of 2 or greater is low enough that no values greater than 1 occurred in the obtained sample.
bmcmedresmethodol.biomedcentral.com
ポアソン分布は0以上の整数値を取る分布ですが，「2以上の整数値を取る確率が非常に小さい」と考えることでbinary outcomeのモデリングを可能にする，と捉えていけば良さそうです．しかし，単にポアソン分布でモデルを組むだけでは精度があまり良くないので，"修正"ポアソン回帰としていく，という考え方になります（信頼区間を広く評価してしまう）．
修正ポアソン回帰の手法はZou (2004) にまとめられています．
では，実際にやってみましょう！
"修正していない"ポアソン回帰の実践
修正ポアソン回帰の前に，"修正していない"ポアソン回帰を実行します．
poisson <- glm(death ~ sex + child, data = df, family = poisson(link = "log"))
summary(poisson)##
## Call:
## glm(formula = death ~ sex + child, family = poisson(link = "log"),
## data = dat)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.2649 -0.7067 0.2151 0.2151 1.1524
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.42303 0.12714 -11.193 <2e-16 ***
## sex 1.16420 0.11132 10.458 <2e-16 ***
## childadult 0.03573 0.10226 0.349 0.727
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 649.48 on 1045 degrees of freedom
## Residual deviance: 508.22 on 1043 degrees of freedom
## AIC: 1752.2
##
## Number of Fisher Scoring iterations: 5
リンク関数にlogを使っているため，計算された回帰係数の指数を取ります．
exp(coef(poisson)) # risk ratio## (Intercept) sex childadult
## 0.2409832 3.2033631 1.0363789性別 (女性に対する男性) のリスク比は3.20倍と，2×2表と概ね相違ない結果が得られました！
信頼区間も見てみましょう！
exp(confint(poisson)) # 95%CI, basic poisson regression## 2.5 % 97.5 %
## (Intercept) 0.1865359 0.307146
## sex 2.5887985 4.007126
## childadult 0.8514672 1.271813むむむ，95%信頼区間は [2.59–4.01] と，2×2表で得られた [CI: 2.70–3.87] より広くなってしまいました．これがZouが指摘しているこの内容ですね．
On the other hand, use of Poisson regression tends to provide conservative results
そこで，分散の推定方法を改善して，"修正"ポアソン回帰に進みます．
修正ポアソン回帰の実践
分散をsandwich estimatorによる方法で求めます．ここが"modified"たる所以ですね！
library(lmtest)
library(sandwich)
modified_poiss <- coeftest(poisson, vcov = sandwich)
modified_poiss##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.423028 0.102390 -13.8981 <2e-16 ***
## sex 1.164201 0.090693 12.8367 <2e-16 ***
## childadult 0.035733 0.061727 0.5789 0.5627
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
信頼区間の計算はパッケージでやってくれないようです．STATAはロバスト分散を指定するだけで計算できたはず．Rも僕の知識不足なだけで，簡単にできるんじゃないかな．
modpoiss_results <- round(cbind(exp(cbind(
RR = modified_poiss[,1],
LCI = modified_poiss[,1] + qnorm(0.05/2)*modified_poiss[,2],
UCI = modified_poiss[,1] - qnorm(0.05/2)*modified_poiss[,2])),
P = modified_poiss[,4]),4)
# RR: リスク比， LCI: Lower CI, UCI: Upper CI. これらはexpする
# P: p-value．これはexpしない
modpoiss_results## RR LCI UCI P
## (Intercept) 0.2410 0.1972 0.2945 0.0000
## sex 3.2034 2.6817 3.8265 0.0000
## childadult 1.0364 0.9183 1.1697 0.5627信頼区間を確認します．
修正ポアソン回帰によるsexのリスク比は 3.20 [95%CI: 2.68–3.83] と，2×2表のRR=3.23 [CI: 2.70–3.87] とそれほど差がない結果になりました．これでokですね！
その他の回帰分析は？
対数二項回帰 (log binomial regression) はおさえるべき分析方法だと思います．これはかなりベーシックなアイデアに基づく回帰モデルです．しかし，この手法はパラメータ空間が制限される（ は
から
を取り，実数全体を取るわけではない）ため，推定アルゴリズムが停止してしまい，エラーを吐いてしまう場合があります．
何を言っているかわからない人は「とりあえず，なんか困る場合があるんだな」で良いんじゃないかな．
ということで，ひとまず第一歩目としては，修正ポアソン回帰を勉強して，対数二項回帰もあるんだな，という理解が良いのではないでしょうか．その先の学習は，このブログの範疇を超える気がするので，自分で頑張っていってください（投げやり）．今回の記事はあくまで導入用ということで．
論文を読む人は...
まとめ的論文はNaimiらの論文．
Zouは修正ポアソン回帰を使うときにひたすら引用されている論文っぽい．
- Ashley I Naimi, Brian W Whitcomb, Estimating Risk Ratios and Risk Differences Using Regression, American Journal of Epidemiology, Volume 189, Issue 6, June 2020, Pages 508–510.
- Zou, G. (2004). A modified poisson regression approach to prospective studies with binary data. American journal of epidemiology, 159(7), 702-706.
