生態学のデータ解析 - FAQ 一般化線形モデル

• ここでは R の glm() を使って解析した場合の説明をしてみます

  

• 参照: FAQ 系ペイジ一覧, GLM 参照, summary(glm()) の星
• この説明もしっかり読みましょう！信頼区間って難しい…

[項目]

• 研究発表で「GLM を使った」と説明するときにはどうしたらよいでしょうか?
• GLM で得られた結果を発表・説明するときにはどうしたらよいでしょうか?
• 説明変数，応答変数って何ですか?
• family で指定する確率分布は「誤差の分布」ですか?
• family 指定はどうすればよいのでしょうか?
• 応答変数のばらつきが family 指定ではうまく表現できないときはどうすればいいのでしょうか?
• (一般化) 線形モデルは必ず交互作用項を含んでいなければならないのですか?
• glm() とかで Y ~ X1 + X1:X2 というふうに X2 は使わないときに交互作用項 X1:X2 だけを使ってもよいのでしょうか?
• glm() なんかの formula 定義で使う * 演算子は何をあらわしてるのですか?
• summary(glm(...)) の出力にある (Dispersion parameter for XXX family taken to be XXX) とは?
• summary(glm(...)) の係数の推定値の表にあらわれる z は何ですか?
• ロジスティック回帰がうまくできない状況は?

研究発表で「GLM を使った」と説明するときにはどうしたらよいでしょうか? 

• よくあるマズい発表の例: 「GLM で解析しました (オワリ)」「GLM の尤度比検定しました (意味不明ぎみ)」
• 一般化線形モデル (GLM) を使った解析では以下のことを示すと良いだろうと思います
  • GLM を使った (あるいはロジスティック回帰，ポアソン回帰など)
  • 応答変数はこれこれ ……
  • 応答変数は二項分布 (あるいはポアソン分布などなど) に従うと仮定した
    • R の glm() の family で何を指定したか，ということ
  • link 関数は logit (あるいは log などなど) を使った
    • R の glm() の family(link=...) で何を指定したか，ということ
  • 説明変数はこれこれ ……
    • そしてその説明変数がどういったタイプなのか? つまり R のタイプで言えば因子型 (factor class) なのか数値型 (numeric class) なのか，といったことなど

GLM で得られた結果を発表・説明するときにはどうしたらよいでしょうか? 

• よくあるマズい発表の例: 「○○がキいてました (オワリ)」
• 必ず示すこと: 解析につかった標本数と推定値一覧
• 推定値 (estimate) に関して
  • 推定値そのものを示すのが最良
  • 発表媒体の余白の関係で省略した書きかたをしなければならないときでも，結果に影響している要因のプラス・マイナスどちらに影響しているのかはっきりさせる
  • 推定値のばらつき (すなわち標準誤差) なんかも示してよいかも
• モデル選択の場合
  • FAQ モデル選択 参照
• 検定の場合
  • 尤度比検定を使ったのか Wald 検定を使ったのか
  • 尤度比検定の場合は帰無仮説における尤度比分布をどのように算出したのか (χ自乗分布近似なのか， parametric bootstrap 法なのか……)
  • FAQ モデル選択 にも書きましたが，検定とモデル選択をまぜてはいけません

説明変数，応答変数って何ですか? 

• glm(y ~ x, ...) の x が説明変数で y が応答変数になります
• 回帰モデルとかについての昔の教科書では独立変数・従属変数とされていましたが，現在の教科書ではそれぞれ説明変数・応答変数となっています
  • 数学用語としては独立変数といった使いかたは必ずしも「まちがい」というわけではないようです (注意: この「独立」は確率論的な「独立」を意味しているわけではありません)
  • しかしながら，いつも混乱にみちている統計学ユーザーわーるどの一部では「複数の独立変数があるときには，それらは必ずたがいに確率論的に独立でなければならない」といったニセ情報がでまわったりしてましたので，安全のため 独立変数といういいかたは廃止したほうがよろしいかと
• 説明変数 (explanatory variable) に対応するのはふつー「応答変数 (response variable)」です
  • ということで，説明変数と従属変数をペアで使わないほうがよいでしょう
  • 応答変数ではなく目的変数 (objective variable?) といったいいかたもある (?) ようですが，これまた (optimization などで使われる) 目的関数とかとまぎらわしいような気もしますので，使わないほうが安全かも
  • 応答変数ではなく被説明変数といったいいかたもありましたが，こんにちでは (さいわいなことに) すたれつつあるようです

family で指定する確率分布は「誤差の分布」ですか? 

• うーむ，「誤差分布をポアソン分布とした」といった記述をよくみかけますけど，これっていったいどうなんでしょうね?
• より無難ないいかたとしては，「応答変数の (したがう) 確率分布」でしょうか?
  • 例文: The response variable Y follows the Poisson distribution of mean m. (応答変数 Y は平均 m のポアソン分布にしたがう)
• というのも，ポアソン分布はカウントデータの「かぞえまちがい (測定誤差，観測者のしくじり)」をあらわしているわけではないのです
  • つまり「ある区域にホントは 4 個体の動物がいたんだけど，まちがって 5 個体とカウントした」といったことをあらわしているわけではありません
  • ポアソン分布があらわしているのは，「この区画の動物の密度は 4.12 なので，個体数が 5 とカウントされる確率は 0.16 である」といったことです
  • したがって「ポアソン分布を誤差分布」と言われると「このヒト，ホントに意味をわかってんのかなぁ」と不安になったりします
  • なお統計学用語である error distribution の (統計学用語としての) error とは必ずしも「測定誤差」だけの意味ではありませんね
• いっぽうで「応答変数の分布」が正規分布である場合は……
  • これは観測者のしくじり，つまり測定誤差をあらわしている場合もあります: たとえば「ホントは 177.5 cm なんだけど 176.9 cm と測定してしまった」といった場合には，誤差分布といういいかたでも異和感が少ない …… かな?
  • しかし「測定誤差だけ」考えればよい状況って少ないような気もします
    • たとえば「1000 人の身長の分布が正規分布に……」という場合には，個体の差と測定誤差がごちゃまぜになっています
    • 「同じ物体の重量を 10 回くりかえして測定したら，測定値は正規分布になった」といった場合には，まあ「測定誤差だけ」をあらわしているんでしょうね

family 指定はどうすればよいのでしょうか? 

• データの構造といいますか応答変数の種類で判別してください
  • binomial: 応答変数が離散値 (カウントデータ) で，とりうる値が {0, 1} もしくは {0, 1, 2, 3, ... N} (上限が N)
    • 「応答した・しなかった」「N 個中 k 個が生残した」といったデータ
  • poisson: 応答変数が離散値 (カウントデータ) で，とりうる値が {0, 1, 2, 3, ...} (上限不明)
    • おなじ平均値をとるデータのばらつきにも注意してください
    • 平均値とばらつきがだいたい同じぐらいであればポアソン分布でよい
    • 平均値よりばらつきのほうがずっと大きければ，たとえば負の二項分布をつかった統計モデルを適用してみてください (glm.nb() in library(MASS))
  • Gamma: 応答変数が連続値で正の値
    • 応答変数にゼロが含まれてるとエラーとなる
  • gaussian: 応答変数が連続値で値は ∞
    • まあ，そうでない範囲のデータに適用されることがほとんどですが

応答変数のばらつきが family 指定ではうまく表現できないときはどうすればいいのでしょうか? 

• 一般化線形混合モデル化・ベイズモデル化してください
  • これによって glm() の family で指定できる確率分布を「まぜあわせて」いろいろな確率分布をあつかえるようになります
  • つまり基本は二項分布やポアソン分布といった確率分布を使ったモデルで，それに手を加えていくといったかんじになります

(一般化) 線形モデルは必ず交互作用項を含んでいなければならないのですか? 

• そんなことはありません
• 交互作用項とは (一般化もふくむ) 線形モデルの説明変数の要因 X1 と要因 X2 の「積」をあらわす項です (R ふうに書くと X1:X2)
• みなさんの大好きな交互作用項ですが……
  • 交互作用項がある GLM というのはある意味「線形モデルとしては失敗している」と言えるかもしれません
  • 交互作用項さえなければ「要因 X1 の影響はこれこれ，要因 X2 の影響はこれこれ」と X1 と X2 の効果を分離して説明できます
  • しかし交互作用項のあるモデルでは (おおくの場合) このような説明が難しくなります

glm() とかで Y ~ X1 + X1:X2 というふうに X2 は使わないときに交互作用項 X1:X2 だけを使ってもよいのでしょうか? 

• そのように使ってもまちがいではありません
• しかし「GLM とか使うってことは現象が線形モデルで表現できると思ってるんだろ? だったら Y ~ X1 + X2 + X1:X2 としろ」と「統計モデルがきちんと線形モデルである」ことを重視するヒトはいます (例えば stepAIC() による自動的なモデル選択，とか)
• 「形式的な線形性なんてどうでもいいんだ」と最初から投げてしまってるアナタは Y ~ X1 + X1:X2 というモデルを使えばいいんじゃないでしょうか (この統計モデルが表現していることをよく理解できているならば)
  • ところで，おそるべきことに，ともうしますか，いわゆる nested ANOVA で使われてる線形モデルはどうもこの Y ~ X1 + X1:X2 であるようです

glm() なんかの formula 定義で使う * 演算子は何をあらわしてるのですか? 

• X1*X2 は要因 X1 と X2 と交互作用 X1:X2 ぜんぶ，という意味です
• つまり下の二つの formula はまったく同じことをあらわしています
  • Y ~ X1*X2
  • Y ~ X1 + X2 + X1:X2
• つまり * 演算子は入力の手間をへらすために存在する演算子なのです

summary(glm(...)) の出力にある (Dispersion parameter for XXX family taken to be XXX) とは? 

• family = gaussian の場合
  • 不偏分散の推定値が表示されています
• family = Gamma の場合
  • Gamma 分布の GLM や KuboLog2012-11-21 などを参照
  • R の help() で gamma.dispersion() や gamma.shape も調べてみてください
• family = poisson または family = binomial の場合
  • ユーザーが summay(glm(...), disp = ...) というふうに指定する dispersion parameter の値が表示されます
  • この disp 指定によって推定値の SE などが変化します

summary(glm(...)) の係数の推定値の表にあらわれる z は何ですか? 

• 係数 (coefficient) の推定値 (Estimate) を標準誤差 (Std. Error) でわった値です
  • 推定値がどれぐらいばらつきそうか，というめやすです
• family = gaussian の場合，z ではなく t になっていて，これは t 検定の検定統計量とだいたい同じ値になります
  • つまり係数の推定値がゼロからずれているかどうかを検定している (t 検定に近い方法で)，ということです (ホントに正しく P 値を計算してるのかどうか，いまいち不明)
  • 下の非 gaussian な場合と同じく，信頼区間にゼロが含まれるかどうかと解釈してもよいでしょう
• family = gaussian ではない summary(glm(...)) の係数の z 値は Wald 統計量であり，Wald 信頼区間 (CI) の推定に使います
  • 参照: Pawitan. 2001. In All likelihood. p.42 http://books.google.co.jp/books?id=M-3pSCVxV5oC&pg=PA42
  • Pr(> |z|) はゼロという値が (Wald) 信頼限界になる信頼水準です
  • つまり Pr(> |z|) が 0.05 より小さい場合，「95% Wald 信頼区間にゼロが含まれない」と解釈できます
• いろいろと参考になる Wikipedia の記述: http://en.wikipedia.org/wiki/Wald_test
• 推定結果として得られた Wald 信頼区間 (CI) については，とりあえず以下のように説明したらいいんじゃないでしょうか? (もとネタ: ぎょーむ日誌 2009-01-08)
  • この確率をネイマン=ピアソンわくぐみの検定のたぐい，と考えてはいけない……つまり，むやみやたらに「ゆーい」「ゆーい」と言わない
  • 「説明変数○○の係数 (パラメーター) の推定値はこれこれで，その 95% Wald 信頼区間 にはゼロが含まれていませんでした (あるいは含まれていました)」とだけ述べ，あとは読者や聴き手の統計学的見識，思慮ぶかさ，あるいは浅はかさや注意力欠如な思考，「わかってないのにわかったふりをする」態度に解釈をゆだねる
  • 注意ぶかい聴き手が「95% Wald CI にゼロが含まれない，とはどういう意味ですか?」と質問してきたら，上記のような教科書的な信頼区間の解釈を述べる，ネイマン=ピアソン的なわくぐみにおける「ゆーいゆーい」とは別モノだと述べる，さらに「あなたは推定値の信頼区間をどう解釈しているのですか?」と質問者の統計学理論の理解のどあいをさぐる

ロジスティック回帰がうまくできない状況は? 

• 例 1: 説明変数が factor であるときに，ある水準の応答変数が全てゼロ (たとえば，対照区はゼロばかり，など)