Ma.08　XxYyの悩みへの対応（演習版に修正前）

データセットを結合できるのことは、データサイエンス的な一つのゴールと言って良いと思っています。
必要な変数を自分で考えて、データセットを構築すること、に意味があると思います。

そのうえで、悩ましい問題に対応していかなければなりません。
XxYy問題です。（僕が勝手に呼んでいるだけです）

Ma01で作成したデータセットは、naimanという個票データ（＝個人の事項）と、「県の合計特殊出生率」「県の自殺率」という集団レベルのデータの2つのレベルの事柄を含みます。

このように、データセットを作るとっていると、別のレベルのものが入り込むことがしばしば起こります。
異なるレベルの変数が混在するデータセットは、「個人問題（レベル1）」と「集団問題（レベル2）」という、悩ましい課題を提示してきます。

個人の結果 (y): 例えば、個人のK6とか、子育ての不満度
集団の特性 (X): 例えば、県の合計特殊出生率とか、自殺率

「集団の特性（X）」が「個人の結果（y）」に影響しているだろう。と、分析しようとする時、注意するべきことがあります。

生態学的誤謬へ（集計データへの逃げ）

「集団の特性（X）」と「個人の結果（y）」が混在する「XxYy問題」に直面したとき、最も単純で手っ取り早い解決策は、分析の単位を「集団」に統一してしまうことです。
つまり、N=10,633ある個人データ（y）をすべて集団レベル（Y）に丸めてしまう（平均値などをとる）アプローチです。

なぜこのアプローチを取りがちなのか？

集団に丸める方法の魅力は、「単純さ」です。
例えば…
①「個人のK6（y）」を47都道府県ごとに平均し、「県のK6平均（Y）」という新しい変数を作ります。
②　加えて、「県の自殺率（X）」という既存の集団レベル変数とを組み合わせます。
　　⇛　N=10600の複雑な階層データは、N=47の非常にシンプルな「県レベル」のデータセットに変わります。
③　あとは、このN=47のデータで相関分析や回帰分析（lm()）を実行するだけです。

なぜこれに注意が必要なのか？

このアプローチは「生態学的誤謬（Ecological Fallacy）」という、問題に当たります。
この分析で分かるのは、あくまで「集団レベルの相関」であり、「個人レベルの因果」ではありません。

• 県の分析（ X → Y）: 「自殺率が高い県は、K6平均も高い県だった」という発見!?
• 個人の現実（x → y）: その県の中でK6が高いのは、自殺率とは無関係な「特定の職業x1」や「特定の経済状況x2」の人だけかもしれません。

集団レベルの分析（X → Y）は、マクロなトレンドや「何が起こっているか」の示唆を与えてくれますが、
それが「なぜ起こっているのか」を個人レベルのメカニズム（X → y やx → y）で語ることを考えると、慎重な扱いが求められます。

つまり、仮説を見つけたり、EDAの範囲として「診る・見る」分にはOKですが、結論を導く際の手法としては避けるべきだと考えています。

Nの肥大化（個票データへの無理な結合）

逆に「個人のデータ（N=10,600）」に「県の自殺率（X）」をそのまま結合し、回帰分析にかけるのはどうでしょうか？

# 危険なコード（やりがち）
model_wrong <- glm(k6_h ~ suicide_rate + age + sex, data = data)

Rの勘違いを助長する

Rは、「N=10,600人分の、バラバラな自殺率情報がある」と勘違いします。
しかし、自殺率の情報は実質N=47種類しかなく、同じ県に住む1,000人は全員が同じ値を持っています。彼らは「独立」ではありません。
この「独立性の仮定」を違反すると、Rはサンプルサイズ（情報量）を過大評価します。
その結果、本当は偶然のノイズに過ぎない差を「統計的に極めて有意（p < .001）」と誤判定します。
（大きすぎるデータでは、ノイズを有意判定するのは、これまでに経験していると思います）

レベルが違う変数をどう扱うべきか？

まず取り組むべきは、扱っている「y（個人）」のばらつきは、どれくらい「集団（X）」によって説明できるか？をチェックすることです。
指標はICC（級内相関係数）を活用します。
この授業では、「K6（個人の精神的な苦痛）」の違いのうち、何%が「住んでる都道府県」の違いで説明できるのか？ってことで考えました。

たったの数行で、重要な示唆を与えてくれます。

library(lme4)
library(performance)
model_null <- glmer(k6_h ~ (1 | tiiki),
                    data = naiman,
                    family = binomial)

summary(model_null)
icc(model_null)

(1 | tiiki) の簡単な説明

これは、Rのlme4パッケージ（glmerやlmer）で使われるおまじないで、「ランダム切片モデル」を意味します。
一言でいえば、「地域（tiiki）ごとに、平均的な水準（ベースライン）が違うことを許可します」という意味です。

• 1: これは「切片 (Intercept)」を意味します。
• |: これは「～ごとに (given)」と読み、「～を条件とする」という意味です。
• tiiki: これは「グループ（クラスター）」となる変数です。

全部つなげると、「切片(1)は、地域(tiiki)ごとに異なってよい」というモデルを指示しています。

なお、lmer（連続変数）でもglmer（二値変数）でも、NULLモデルの書き方は (1 | tiiki) だけでOKです。

演習では、 K6の分析で、「A県はもともとメンタルが良い（ベースK6が低い）」「B県はもともとメンタルが悪い（ベースK6が高い）」といった、県ごとの「個性」をモデルに組み込むことを想定していました。
（本当は、「雪国」という変数で分析する予定だったのですが、、、悔やまれます。）

ICCの判定

繰り返しになりますが、ICCは、分析対象「個人の結果（y）」の総分散（全体のばらつき）のうち、どれだけの割合（%）が、所属する「集団（グループ）」の違いによって説明できるのかを定量的に把握する指標です。
ICCは、0から1までの値をとり、ICC>0.05くらい(総分散の5%)で所属する集団がyに効いてくると考えます。

ICCを見て、5%以上あれば、集団の構造を無視した通常の回帰分析（glmやlm）では、観測値の独立性の仮定が崩れ、標準誤差の過小評価という問題が生じます。
したがって、階層性をモデルに組み込むことが可能なマルチレベル分析（混合効果モデル（glmerやlmer））が適切です。

逆に、ICCが5%未満の場合は（演習の場合はほぼゼロでした）、集団の影響をマルチレベル分析を用いることはあまり意義がありません。
とはいえ、集団からの影響をどうしても無視できないと判断したときには、私はANCOVAに補正をかけて解釈を試みます。（後述）

マルチレベル分析の一例

※例として提示しているモデルはICC =0.0005のマルチレベル分析するべきでないものです。以下は手順とコードを示すのが目的です。

model2 <- glmer(k6_h ~ 
                     q1_1 + age_num + sex 
                     + (1 | tiiki),
                     data = naiman,
                     family = binomial)
summary(model2)

上記のように、Null_modelに個人要素（q1_1 + q21c + age_num + sex）の最低限のものをいれたものをモデル2に設定します。

その後、仮説に基づいて、個人要素に関連する説明変数を入れていき、相対的に比較していきます。(model3→model4)
最終モデルにたどり着くまでは、summaryだけで見る形でも問題有りません。

ただし、最終モデルでは、各変数のOR、信頼区間を計算する必要があります。
EstimateからORを計算するのは面倒なので、Rをにやってもらうために、
tidy(model2, exponentiate = TRUE, conf.int = TRUE)　# exponentiate = TRUE estimateをオッズ比に変換してもらいます。
また、performance::r2(model2)でモデルの評価も同時に行います。

library(performance)
library(broom.mixed)
tidy(model4, exponentiate = TRUE, conf.int = TRUE)　# exponentiate = TRUE オッズ比に変換
performance::r2(model4)

すると、以下のようなモデル評価指標がでてきます。
  Conditional R2: 0.349
     Marginal R2: 0.348

Conditional R^2（条件付きR^2）: 固定効果（Fixed effects）とランダム効果（Random effects）の両方で、どれだけK6の分散を説明できたか？
Marginal R^2: 固定効果（Fixed effects）だけで、どれだけK6の分散を説明できたか？

今回の演習のモデルではRandom effect（tiiki）は何もしていないということですね。。。。


ここまできて、地域の効果が残っている場合は、ようやく地域Xの因子を投入できます。
演習の例では、jisatu（県の自殺率）です。（さすがに、すでに地域の影響はないとなっているので苦しいですが。。。）

model5 <- glmer(k6_h ~ 
                  q1_1 + q21c + age_num + sex
                + jisatu
                + (1 | tiiki),
                data = naiman_ketsugo,
                family = binomial)
summary(model5)

ICCがほぼ0だが、集団の影響を無視できない場合

シンプルな回帰モデル（＝ANCOVA） glmやlmを使って、モデルを組みます。
その際に、Nが肥大化することを承知しながら　jisatu （県の自殺率）などのXの要素も入れ込みます。（オススメしません）
この場合は、最後に調整を要します。（オススメしません）

library(lmtest)
library(sandwich)

model_ancova <- glm(k6_h ~ 
                  　q1_1 + q21c + age_num + sex
                  　+ jisatu, 
                  　data = naiman_ketsugo,   family = binomial)

# coeftestでp値を補正する
coeftest(model_ancova, 
         vcov = vcovCL(model_ancova, cluster = ~tiiki))

このようにすることで、単純なANCOVAに異質のXを投入することができました。
最初からMLAを無視して、ANCOVAに投入したら良いのではないか、という議論はあると思います。
いくら微々たる調整をしたところで、「そもそも個人レベルと集団レベルを混ぜるべきでない」という正当な批判に耐えることは難しいです。

少なくとも演習の具材としては、適切では有りませんでしたが、今後の分析で「人種」のような
影響力が強いグループを変数に入れる場合には必要な処置と思いますので必要な方に見ていただければと思います。