統計勉強会
～GLM（一般化線形モデル）～
2010.3.8
中野貴瑛（理論生態学研究室）
目次
• はじめに
• 統計モデルとは
• モデル選択
• 一般化線形モデルとは
• どんなGLMを使おうか
• GLMM（一般化線形混合モデル）とは
• モデルの検査
• モデル選択の悩みどころ
• 調べるとよさそうなもの
• 架空データを用いた解析例
2
はじめに
• このスライドの目的
– 一般化線形（混合）モデルGLM（M）の概念をつかむ
– GLM（M）を使いたいときに、あまり困らないで済むようにする
• 注意事項
– 不確かな情報がある可能性あり
– あくまで参考程度にとどめること
• 参考になりそうなもの
– 北大・久保さんのページ http://hosho.ees.hokudai.ac.jp/~kubo/ce/FrontPage.html
3
統計モデルとは
•
観測データを解析するということは、統計モデリングをすることである
•
統計モデル：観測データのパターンをうまく説明できるようなモデル
•
たとえば、ちょっと歪んだサイコロを50回振ったとする。このとき・・・
– どの目が出る確率も、同じ 1/6 だ
– 目の値が小さい方が出る確率が高い
0.5
出目の確率分布（n=50）
0.4
などが統計モデル。
0.2
0.1
2の目が出る確率が有意に高いだろう、
などと言って検定をすれば、それは
「どの目が出る確率も、同じ 1/6だ」という
モデルを帰無仮説とした検定をしている
ことになる。
0.0
–
0.3
実は検定でも 統計モデリングをしている
頻度
•
1
2
3
4
5
6
目の値
4
モデル選択
検定
モデルを作ったら
モデルAとモデルBが有意に異なってる？
どっち？
モデル選択
最適なモデルはどれかな？
モデル選択とは
考えうるモデルの中から、簡潔かつデータによく当てはまっているモデルを選ぶこと
•
そうして得られた最適なモデルの挙動や予測性をみたり、
他のモデルとの比較を行ったりして考察する
•
よく用いられる基準に、AIC（赤池情報量基準）がある
•
具体的な方法は後述
5
モデル選択
検定のあやしさ
• 例：ある魚の体長が、一腹卵数に関係していそう。
• こんなときに、こんなことしていない？
80
回帰直線を引く
40
20
2
0
4
有意差！
ワーイ＼ (≧∇≦)／
15
20
25
30
35
40
45
体長（cm）
0
√(一腹卵数)
6
8
相関係数の検定
一腹卵数（個）
60
10
平方根変換
15
20
25
30
35
40
45
体長（cm）
体長が小さくなったら
卵数がマイナスになっちゃうよ？？
関係があるかないかも確かに重要だけど、
どんな風に関係しているかが知りたいのでは？
6
モデル選択
検定のあやしさ
• パラメトリック検定では
– 正規分布に従い等分散性のあるデータだけしか扱えない
• パラメトリック検定が無理なら、ひとくくりにノンパラメトリック検定
– 検出力の低下・等分散性のなさは気にしない
• さらに、データをあれこれいじることも（区分け・割り算・変換etc...）
– そうやって出てきた有意差って意味あるの？
• そもそも、有意水準って自分で決めるもの
– 本当に客観的なの？
• ある説明変数が有意に効いてるってことも大事だろうけど・・・
– それがどのような振る舞いをしているのかが、生態学にとって重要なのでは？
7
モデル選択
検定とモデル選択
•
•
「検定はダメダメでモデル選択は万能」というわけではない
使い分けが重要
•
仮説をしぼった実験（調査）から、
新しい要因や関係の存在を確実に示したいとき
検定
– 例）新薬が効いているのかを知りたい
•
予測精度が高いモデルを作りたいときや、
仮説を対等に見て、どちらかを選びたいとき
•
生態学、とくに野外での非実験的な調査・観察から得たデータは
検定には不向きである場合が多い
モデル選択
GLM(一般化線形モデル）を知って、
モデリング＆モデル選択をしよう！
8
一般化線形モデルとは
GLM（Generalized Linear Model：一般化線形モデル）
25
• 根本的な概念は直線回帰と同じ
5
10
右の例におけるパラメータは
切片と傾きの２つ
15
正規分布を仮定した最尤推定は
特に最小二乗法と呼ばれる
応答変数
20
– 最尤推定法でパラメータを推定
•
右の例は、 yi 1 2 xi というモデル
5
10
線形予測子
•
15
20
25
説明変数
しかし、さらに幅広いことができる
– 応答変数が正規分布に従わなさそうな場合にも対応可能
9
一般化線形モデルとは
ポアソン回帰モデルとも呼ぶ
log e yi 1 2 xi
60
応答変数
なぜexp?
→カウントデータは非負だから
20
yi exp 1 2 xi
0
モデル：
40
線形予測子
x = xi のときの平均値
80
• たとえば応答変数がカウントデータの時：
と表現できるので、link関数はlog
10
15
20
25
30
35
40
45
説明変数
ロジスティック回帰モデルとも呼ぶ
log e
•
1 exp 1 2 xi
yi
logit ( yi ) 1 2 xi
1 yi
0.6
0.4
yi
線形予測子
0.2
モデル：
1
0.0
x = xi のときの成功確率
応答変数
0.8
1.0
• たとえば応答変数が成功・失敗データの時：
0
と表現できるので、link関数はlogit
1
2
3
4
説明変数
いずれも、link関数というものを用いることで、線形の土俵に持ち込むことができる
10
一般化線形モデルとは
さらに・・・
• 複数の説明変数があってもOK
– いわゆる“重回帰“的なことも可能
• 例）ある植物の結実数を、日当たり・土中のリン濃度・植物の高さ、で説明
• 説明変数がカテゴリカルでもOK
– いわゆる“分散分析”的なことも可能
• 例）処理A～Cによる違いをみる
• 一口にGLMと言っても、非常に幅広い解析ができる
• もちろん、複数の説明変数にカテゴリカル変数が混ざっていてもよい
11
どんなGLMを使おうか
よく用いられる確率分布
• 正規分布
– 応答変数が連続データ
– 負の値を許す
• ポアソン分布
– 応答変数が離散データ
– 0以上で上限とくになし
– 平均と分散がほぼ等しい
• 二項分布
– 応答変数が離散データ
– 0以上で有限
応答変数が従っていそうな確率分布を考える
最終的に示したいモノ
魚Aの体長、動物Bの体重、
鳥Cの採餌にかかる時間 etc…
両生類Dの卵数、アブの訪花個体数、
イラガの繭数、トカゲの尾振り回数 etc…
生存率、成功率、
n個のうち何個が○○だった、みたいな
データ（要するに確率が知りたい） etc...
• とりあえずはこの3つで十分かも
• 確率分布はモデルの骨格になるので、しっかり勉強しよう・・・
12
どんなGLMを使おうか
すべては応答変数の確率分布を
仮定するところから始まる・・・
応答変数が・・・
正規分布： １
-∞ ~ +∞
値の範囲が
連続値？
ガンマ分布： ５
0 ~ +∞
0~1
離散値？
ベータ分布： ０
No
0, 1, 2, ・・・, N
値の範囲が
二項分布： ２
0 , 1, 2, ・・・, +∞
過（大）分散？
over-dispersion
ポアソン分布： ３
過（大）分散？
over-dispersion
Yes
負の二項分布： ４
or
GLMM： Z
No
ゼロがやたらと｛多く・少なく｝
No しかも思い当たる理由あり
ZIP model： α
Yes
ベータ二項分布： ０
or
GLMM： Z
hurdle model： β
13
どんなGLMを使おうか
前ページからの続き
Rの、パラメータ推定関数
０ ：ベータ分布・ベータ二項分布
不明
１ ：正規分布
glm(family = gaussian), lm()
２ ：二項分布
glm(family = binomial)
３ ：ポアソン分布
glm(family = poisson)
４ ：負の二項分布
glm.nb() @library(MASS)
５ ：ガンマ分布
glm(family = gamma)
α：ZIP(zero inflated poisson) model
zeroinfl(dist = "poisson") @library(pscl)
β：hurdle model
hurdle(dist = "poisson") @library(pscl)
Z ：GLMM(Generaized Linear Mixed Model)
glmmML(family = binomial or poisson) @library(glmmML),
lmer(family = gaussian) @library(lme4)
特殊なポアソン回帰モデル
GLMMで手に負えないデータは、
階層ベイズモデルで・・・
14
GLMM（一般化線形混合モデル）とは
Generalized Linear Mixed Model
• GLMでは考慮していなかったランダム効果を組み込んだモデル
–
ランダム効果：調査地間のばらつきとか、定量化できていない個体差とかの、
興味はないけど存在はするばらつき
• ランダム効果が一つなら、glmmML() @library(glmmML)がよさげ
• lmer() @ library(lme4)は、複数の・あるいはネストされたランダム効果も
組み込み可
– でも挙動不審な場合がある？
• ランダム効果が複数あるなら、ベイズ推定がオススメ（R + WinBUGS）
15
モデルの検査
•
モデルをたくさん作ってモデル選択をして、ベストモデルが分かればそれでよし？
•
モデル選択は相対的なもの
ベストモデルが本当にしっかりとしたものなのかを検査する必要がある
•
過大分散と過小分散
–
そのモデルから予想されるよりも、実際のデータのばらつきが大きい（過大分散）
•
•
–
そのモデルから予想されるよりも、実際のデータのばらつきが小さい（過小分散）
•
•
重要な説明変数を見落としている可能性
ランダム効果を考慮する必要があるかも
サンプルの独立性が保たれていない可能性
–
–
実際によくあるのは過大分散
単に過分散という場合には過大分散を指すことが多い
–
過分散について議論できるのは、平均が決まると自動的にばらつきも決まってしまう確率分布（ポ
アソン分布、二項分布）を用いたモデルだけ
他にもあるはず・・・
–
モデルの検査はないがしろにされているのが現状？
16
モデル選択の悩みどころ
結論がビシッときまらない
•
検定ならば・・・
–
•
「有意差があったので、これこれは効いてます！」と言える
モデル選択だと・・・
–
「ベストモデルに説明変数A が含まれていたということは、A を含んでいないモデルよりもベストモデルが現状をよく
説明出来ていたわけで、つまりそれは、A が重要であることを意味しているんですよ・・・」
–
このスタンスが納得できないと、あるいは他の良い結論方法がわからないと、
モデル選択は難しいかもしれない（気持ち的に）
17
調べるとよさそうなもの
•
説明変数、応答変数
•
固定効果（fixed effect）、ランダム効果（random effect）
•
AIC（赤池情報量基準）
•
尤度
•
最尤推定法
•
確率分布（特に、正規分布・二項分布・ポアソン分布）
•
信頼区間
•
モデルの検査方法
•
リンク関数、線形予測子
18
架空データを用いた解析例
解析の流れ
A~Eの一部あるいは全部を、２つの例で紹介
仮説を立てる
この流れ・方法が絶対というわけではない。
あくまで参考として。
データ収集
A
解析しやすいように
データをまとめる
B
色々と作図をして
傾向をチェック
・仮説→絞ったデータ収集→解析
・仮説→たくさんデータ収集→解析
・たくさんデータ収集→解析→言えそうなことを探す
どれがいいのだろうか・・・
No
C
Rに放り込んで
モデルを量産
D
モデル選択と
その結果の確認
OK？
Yes
考察・まとめ
E
モデルの検査
19
架空データを用いた解析例（ケース１）
ケース１
ある魚Aを川で計100匹捕獲し、一腹卵数を数えた。
一腹卵数のモデルを作る。
収集したデータ：①卵数、②親魚の体長、③捕獲した場所（上流・中流・下流）
A. 解析しやすいようにデータをまとめる
1.
エクセルなどで入力
データフレームを意識
2.
項目名
テキスト（タブ区切り）形式で保存
・項目名はアルファベット推奨
・カテゴリカルデータでは、文
字列を使う（Rが数値と勘違い
しないように）
・場合によってはカテゴリカル
データはダミー変数化したほう
がいいかも（アブ長談）
20
架空データを用いた解析例（ケース１）
A. 解析しやすいようにデータをまとめる
3.
Rに取り込む
これで準備は完了
21
架空データを用いた解析例（ケース１）
B. 色々と作図をして傾向をチェック
pairs(data1, panel=panel.smooth)
300
500
1.0
1.5
2.0
2.5
3.0
60 80 100
100
plot()で一つずつ作図してもいいけど、
500
0
20 40
No.
300
pairs()で対散布図が描ける。
eggs
10 15 20 25 30
100
panel=panel.smooth で平滑線が引ける
関係がありそう！
2.0
2.5
3.0
5
size
1.5
site
もし説明変数間に相関が見られた場合、同じような影
響の二重評価などにより、正確な予測ができなくなる。
1.0
その場合、説明変数を減らすなどの対策が必要。
0
20
40 60
80 100
5
10 15 20 25 30
22
架空データを用いた解析例（ケース１）
B. 色々と作図をして傾向をチェック
attach(data1)
plot(size, eggs, type = "n")
points(size[site=="upper"], eggs[site=="upper"], col = 2)
points(size[site=="middle"], eggs[site=="middle"], col = 3, pch = 20)
500
points(size[site=="lower"], eggs[site=="lower"], col = 4, pch = 18)
400
カテゴリカル変数を色分けで対処してみたり・・・
100
200
eggs
300
他にも、思いつく限り作図をしてみることが大切
5
10
15
20
size
25
30
23
架空データを用いた解析例（ケース１）
C. Rに放り込んでモデルを量産
今回の場合、
おそらく５通りのモデルが考えられる
線形予測子が切片のみの場合、“１“と記述
扱うデータフレームを指定する
応答変数を最初に書き、チルダ
以降に説明変数を記述
仮定する応答変数の確率分布
model1 <- glm(eggs~1
, data = data1, family = poisson)
model2 <- glm(eggs~size
, data = data1, family = poisson)
model3 <- glm(eggs~site
, data = data1, family = poisson)
model4 <- glm(eggs~size + site
, data = data1, family = poisson)
model5 <- glm(eggs~size + site + size:site, data = data1, family = poisson)
交互作用項
全部まとめて、size*siteと書くことも可能
24
架空データを用いた解析例（ケース１）
D. モデル選択とその結果の確認
＊ 全てのモデルのAICを確認→選択
＊ stepAIC() @library(MASS) で楽をする
はじめはこれがオススメ
＊ update() を使う
summary(モデル)で、結果を確認
・詳細情報は個別にアクセス可能
・項目名はnames(モデル)で分かる
25
架空データを用いた解析例（ケース１）
D. モデル選択とその結果の確認
model
model1
model2
model3
model4
model5
AIC
6128.3
1545.3
6116.2
1532.6
1533.4
model2 <- glm(eggs~size
, ...)
model4 <- glm(eggs~size + site
, ...)
model5 <- glm(eggs~size + site + size:site, ...)
・AICが低いモデルほど、簡潔かつデータによく当てはまっているモデル
・一番AICが低いモデルだけを見るのではなく、それに近いAICのものもチェック
・それらのモデルが持つ項やパラメータ等を比較し、考察へつなげる
26
架空データを用いた解析例（ケース１）
D. モデル選択とその結果の確認
解析結果の見方
カテゴリカル変数は、一つのカテゴリーが基準（つまり係
数は0）として扱われる。そして、それ以外のカテゴリーで
基準に対する係数がそれぞれ推定される。
係数
推定値
切片
標準誤差
p値
p値は、その推定値が0から有意に離れているかどうか
の目安。でも、かなりいい加減らしい。
量的な変数
カテゴリカル変数
この結果では、siteがmiddleの場合の係数が0に近い
ように見える。そこで、middleとlower(基準)を合わせ
て一つのカテゴリーとしたモデルを作り、AICを比較し
てみる、という方法が考えられる。
27
架空データを用いた解析例（ケース１）
D. モデル選択とその結果の確認
カテゴリカル変数（site）のカテゴリーを1つ減らした新たなモデル。
AICが小さくなったため、このモデルを最適なモデルとして、次に検査を行う。
lower と middle を合わせてmidlow という名前にした。
それが基準となり、それに対する upper の係数が推定されている。
※カテゴリーの変更に伴い、site→site2, data1→data2 と名前を変えた。
推定された係数から、この魚Aの卵数の予測は
体長（size）を用いて以下のように表せる
中下流： exp 3.10031 0.09702 * size
上流： exp 3.10031 0.09702 * size 0.06839
28
架空データを用いた解析例（ケース１）
D. モデル選択とその結果の確認
ベストモデルの予測を図示
推定された係数から、この魚Aの卵数の予測は
体長（size）を用いて以下のように表せる
赤が中・下流、青が上流の個体
eggs
300
400
中下流： exp 3.10031 0.09702 * size
上流： exp 3.10031 0.09702 * size 0.06839
100
200
中下流の魚の卵数と上流のそれとでは
大差がないようにも見える。
10
15
20
25
30
size(cm)
29
架空データを用いた解析例（ケース１）
E. モデルの検査
とりあえず今回は・・・
過分散のチェック
> ovdis(model6, 0.99)
confidence interval level: 99 %
outers: 41 in 100 samples
500
予測値に対する応答変数の実測値及びモデルの信頼区間
300
99%信頼区間から外れているプロットが明らかに
全体の1%以上存在している
200
過分散
GLMMにする、負の二項分布でのGLMにする、など
の対策が必要か
100
応答変数の実測値
400
「予測値が大きいほど実測値が予測値から外れる」
などの傾向は見られない
満足する？ Bへ戻る？
100
200
300
predictive value(予測値)
400
ケース１ 終
30
架空データを用いた解析例（ケース２）
ケース２
ある植物Bは、花を咲かせた後に1つの実をつけるが、つけない場合もあることがわかっている。
そこで、植物Bを100本観察し、花の数を数えた。
その数日後、つけていた実の数を数えた。
結実確率のモデルを作る。
収集したデータ：①花の数、②実の数、③植物の高さ（cm）
A. 解析しやすいようにデータをまとめる
1.
エクセルなどで入力
データフレームを意識
2.
テキスト（タブ区切り）形式で保存
知りたいのは結実確率だが、
fluits / flowers のように勝手に割り算しない
分母がサンプルごとに異なっていても平気
31
架空データを用いた解析例（ケース２）
Bは飛ばして
C. Rに放り込んでモデルを量産
今回の場合、
おそらく２通りのモデルが考えられる
応答変数部分を、
cbind(結実した数, 結実しなかった数)
となるように記述
確率を扱いたいので、二項分布を仮定する
flowers ではなく、「結実しなかった数」を
データ項目としておいてもよい
model1 <- glm(cbind(fluits, flowers-fluits)~1
, family = binomial, data = case2)
model2 <- glm(cbind(fluits, flowers-fluits)~height, family = binomial, data = case2)
～した
有り
成功
・
・
・
～しなかった
無し
失敗
・
・
・
D. モデル選択とその結果の確認
model1
model2
AIC
1006.6
149.4
・AICが低いモデルほど、簡潔かつデータによく当てはまっているモデル
・今回の場合、heightを説明変数にもつmodel2がベストモデルと言えよう
32
架空データを用いた解析例（ケース２）
D. モデル選択とその結果の確認
推定された係数から、この植物Bの結実
確率の予測は植物の高さ（height）を用い
て以下のように表せる
1
1 exp 14.68065 0.29346 * height
33
架空データを用いた解析例（ケース２）
D. モデル選択とその結果の確認
0.6
推定された係数から、この植物Bの結実
確率の予測は植物の高さ（height）を用い
て以下のように表せる
0.2
0.4
1
1 exp 14.68065 0.29346 * height
0.0
結実確率の予測
0.8
1.0
ベストモデルの予測を図示
30
40
50
60
height(cm)
70
80
34
架空データを用いた解析例（ケース２）
E. モデルの検査
今回も・・・
過分散のチェック
> ovdis(model2, 0.99)
confidence interval level: 99 %
outers: 0 in 100 samples
15
予測値に対する応答変数の実測値及びモデルの信頼区間
5
10
過分散ではない
0
応答変数の実測値
99%信頼区間から外れているプロットはひとつもない
0
5
10
predictive value(予測値)
15
満足する？ Bへ戻る？
ケース２ 終
35
おまけ
過分散とかを調べるところで使ったovdis()
#過分散（overdispersion）を調べる関数。lmを継承したモデリング関数なら代入して平気なはず。
#model：過分散を調べたいモデル, level：信頼区間
#family = poisson(link = log), binomial(link = logit)) にのみ対応
ovdis <- function(model, level = 0.95){
if((level < 0) || (level > 1)){
cat("level: error\n")
}else{
if(model$family$family == "poisson"){
Y <- model$y
pred <- model$fitted
area.over <- qpois(1 - (1-level)/2, lambda = pred)
area.under <- qpois(0 + (1-level)/2, lambda = pred)
}else if(model$family$family == "binomial"){
Y <- model$model[1][,1][,1]
N <- model$model[1][,1][,2]
YN <- Y+N
pred <- model$fitted*YN
area.over <- qbinom(1 - (1-level)/2, size = YN, prob = pred/YN)
area.under <- qbinom(0 + (1-level)/2, size = YN, prob = pred/YN)
}
info <- data.frame(Y, pred, area.over, area.under)
info <- info[order(info$pred),]
with(info,{
plot(pred, Y, xlab = "predictive value(予測値)",
main = "過分散チェック", ylab = "応答変数の実測値＆モデルの信頼区間")
points(pred, pred, type = "l")
lines(pred, area.over, lty = 2, col = 2)
lines(pred, area.under, lty = 2, col = 4)
})
samples <- length(Y)
outers <- length(Y[Y < area.under | area.over < Y])
cat(c("confidence interval level:", level*100, "%\n"))
cat(c("outers:", outers, "in", samples, "samples\n"))
}
}
36
© Copyright 2024 ExpyDoc
