GLMの予測区間
～直感的理解と図示方法～
var. 2010.3.24
目次
• はじめに
• family = poisson 編
– ケース1（説明変数がない場合）
– ケース2（説明変数が1つある場合）
– ケース3（説明変数が複数ある場合）
• family = binomial 編（以下、未完成・・・）
2
はじめに
• このスライドの目的
– タイトルの通り、GLMの予測区間について理解を深め、
Rによる作図ができるようになること。
•
前提とするもの
– よくつかわれる確率分布に関する、ある程度の知識
– GLM(一般化線形モデル)に関する、ある程度の知識
•
予測区間とは
– 母集団を仮定した上で、将来観察されるであろう標本値に対して「どの範囲にあると予
測されるか」を示すものである。（wikipedia より）
• これに対し、信頼区間とは、母集団の母数（標本から測定できない）に対して 「どの範囲にあ
ると推定できるか」を示すものである。混同しないように注意。（wikipedia より）
•
その他
– “GLMの”予測区間と言ってるけど、考え方は他のモデルでも同じかも。
– ただし、ランダム効果入りのモデルは別。要注意。
3
family = poisson 編
• GLM(一般化線形モデル)は、“一般化”とつくだけあって、
いろいろな確率分布を扱うことができる。
• ここでは、その中でも理解しやすいと思われる、ポアソン分
布に従ったデータを題材として解説していく。
ケース１（説明変数がない場合）
•
まさお君は釣りが好きであった。今までに手賀沼に釣りに行った回数はちょうど
50回である。まさお君は、よく釣れるブルーギルの数を記録していた。そのデータ
を利用して、手賀沼で釣れるブルーギルの数をモデリングすることにした。ちなみ
に、釣りをする時間は毎回一定であった。
釣れたブルーギル数のヒストグラム
gill <- rpois(50, 5) #ポアソン分布に従う乱数生成
hist(gill)
#ヒストグラム作成
mean(gill)
#平均
var(gill)
#分散
平均：4.74
分散：4.52
平均≒分散となっていて、ポアソン分
布に従っていそう！
注）乱数を使っているので、実際に試した場合とは微妙に値が異なると思います。
5
ケース１（説明変数がない場合）
•
•
釣れたブルーギルの匹数を応答変数としてGLMを作成。
説明変数は無し（平均値のパラメータ:beta のみ）。
model1 <- glm(gill ~ 1, family = poisson) #GLMを作成
beta1 <- model1$coefficients
#このモデルのbetaの値を格納
#このモデルの確率分布を調べる
pred1 <- c()
#確率を入れる変数を用意
for(i in 0:max(gill)){
pred1[i+1] <- dpois(i, exp(beta1))
}
#lambda = exp(beta)のとき、
#釣れる数がi匹である確率を求める
#ちなみにこれは、for文を用いずに書くことも可能↓
pred1 <- sapply(0:max(gill), function(i) dpois(i, exp(beta1)))
6
ケース１（説明変数がない場合）
•
そうして求まった確率に、サンプル数をかければ予測値となる。
#hist()だとx軸がいい加減になっちゃうので、度数分布データに直して改めて作図
plot(table(factor(gill, levels = 0:max(gill))))
#先ほど求めた確率にサンプル数をかけて出した予測値を重ね描き
lines(0:max(gill), pred1*50, type = "b")
けっこういい感じな予測？
そして、次からがいよいよ予測区間！
7
ケース１（説明変数がない場合）
•
まず、話を予測値からに確率分布に戻す（サンプル数をかけていない状態）
plot(0:max(gill), pred1, type = "b")
この確率分布を見て分かること
・4匹釣れる確率が一番高い
・0匹や10匹の確率は低い
→ほとんど起こらない両側2.5%ずつに入る
部分を調べて、95%予測区間を求めよう！
（95%というのは、よく使われているだけで
あって、目的に応じて何%でも構わない。）
8
ケース１（説明変数がない場合）
•
両側の2.5%の境目を調べるには qpois() を使う。
qpois(0.025, exp(beta1))
#qpois(分位数, lambda) という風に使う
qpois(0.975, exp(beta1))
#結果、2.5%の方は1、 97.5%の方は9であった。
つまり・・・
この２本の線より内側の部分が、95%予測区間。
はみ出ているのは10匹釣れた2サンプルのみ
→ 2/50 = 4% なので、このモデルの予測は悪くないといえる。
ちなみに、はみ出ているサンプルが5%より遥かに多いと、過分散の可能性が疑われる
9
ケース２（説明変数が１つある場合）
•
まさお君は釣りが好きであった。今までに霞ヶ浦に釣りに行った回数はちょうど50
回である。まさお君は、よく釣れるブラックバスの数を記録していた。そのデータを
利用して、霞ヶ浦で釣れるブラックバスの数をモデリングすることにした。ちなみに、
釣りをする時間は毎回一定であった。さらに今回は、近くで釣りをしていた人の数
も記録してあり、データとして使うことができる！
#下準備・・・
beta2 <- c(1.5, -0.03) #架空データのパラメータ
people <- rpois(50, 2) #近くで釣りをしていた人の数を適当に生成
#パラメータと人の数から、釣れたブラックバスの数を生成
bass <- rpois(50, exp(beta2[1] + people * beta2[2]))
d <- data.frame(people, bass) #なんとなくデータフレームにしてみる
・
・
・
10
ケース２（説明変数が１つある場合）
•
今回のデータを図示してみよう！
plot(d$people, d$bass) #横軸に人の数、縦軸に釣れたブラックバスの数
なんとなく、近くで釣りをしていた人
が多いほど、ブラックバスの釣れた
数が少ないように思える・・・
でもちょっと待って！？
この図とケース1のような図との
関係、わかってますか？
11
ケース２（説明変数が１つある場合）
•
見えているプロットの数を数えると32個。でも、サンプル数は50だったはず・・・
→いくつかのプロットは重なってしまっている
→それが、この図では見られない“度数”
people = 0 のときのbassのヒストグラム
people = 1 のときのbassのヒストグラム
people = 2 のときのbassのヒストグラム
people = 3 のときのbassのヒストグラム
この図は、peopleが0~6のときのbassのヒスト
グラムを並べて、上から見たような図である！
12
ケース２（説明変数が１つある場合）
•
ということは、peopleが0~6の場合ごとに
確率分布あるいは予測値を考えることができる。
people = 0 のときのbassの予測値（てきとー）
people = 1 のときのbassの予測値（てきとー）
people = 2 のときのbassの予測値（てきとー）
people = 3 のときのbassの予測値（てきとー）
ということは・・・・・・
13
ケース２（説明変数が１つある場合）
•
予測区間も同様に考えることができる！
people = 0 のときのbassの95%予測区間（てきとー）
people = 1 のときのbassの95%予測区間（てきとー）
people = 2 のときのbassの95%予測区間（てきとー）
people = 3 のときのbassの95%予測区間（てきとー）
イメージがつかめれば、もう
予測区間は描けたも同然！
14
ケース２（説明変数が１つある場合）
•
•
釣れたブラックバスの匹数を応答変数としてGLMを作成。
説明変数はpeople（パラメータは切片及びpeopleの係数）。
model2 <- glm(bass ~ people, family = poisson, data = d) #GLMを作成
beta2 <- model1$coefficients
#このモデルのパラメータを格納
•
beta2には、切片及びpeopleの係数を格納。
•
モデルから得られたパラメータを使って、横軸（people）がとる値ごとに95%予測区
間を調べていく。
#people = i のときの、bassのヒストグラムの左側2.5%の境目
pred2.1 <- sapply(0:max(d$people), function(i){
qpois(0.025, exp(beta2[1]+ i * beta2[2]))
})
#people = i のときの、bassのヒストグラムの右側2.5%の境目
pred2.2 <- sapply(0:max(d$people), function(i){
qpois(0.975, exp(beta2[1]+ i * beta2[2]))
})
15
ケース２（説明変数が１つある場合）
•
得られた予測区間の両側はこんな感じ
•
あまりパッとしないけど、とりあえず作図してみる。
→→→
plot(d$people, d$bass)
lines(0:max(d$people), pred2.1, type = "b")
lines(0:max(d$people), pred2.2, type = "b")
この２本の線より内側の部分が、95%予
測区間。
はみ出ているのは4サンプル（見えるプ
ロットは3つだけど、先に述べたようにこ
の図では“度数”を確認することができな
い）。
→ 2/50 = 4% なので、このモデルの予
測は悪くないといえる。
では、説明変数が2つ、3つと増えた場合にはどうす
ればよいのだろう？ →次のスライドへ
16
ケース３（説明変数が複数ある場合）
•
•
説明変数が２つであれば、３次元の図を描けないこともない。
しかしそれ以上となると、各軸に説明変数をとることは厳しい。
•
ところで、「ポアソン分布に従うデータのモデル」の予測区間を決めているのは
いったいなんだろうか？
答え･･･ ポアソン分布のパラメータ：λ
•
•
だったら、横軸にλ、縦軸に応答変数をとればいいじゃないか。
今さらながら、ポアソン分布のパラメータλは、その分布の平均と分散を意味する。
そして、線形予測子を z とするならば、λ＝exp(z) である（リンク関数がlogの場合）。
•
まぁ、難しい話はおいといて、実践・・・
17
ケース３（説明変数が複数ある場合）
•
一日に見つけた鳥の数（tori）が、平均気温（kion）、雲量（kumo）、雨の有無
（ame）で決まるだろう、とかいう適当すぎるモデルを考える。
kion <- round(rnorm(20, 20, 5), 1)
#平均気温（℃）
kumo <- round(rnorm(20, 40, 20), 1) #雲量、単位は謎
ame <- rbinom(20, 1, 0.4) #雨が降った:1、降らなかった:0
tori <- rpois(20, exp(0.1 + 0.04 * kion + 0.04 * kumo - 0.06 * ame))
d3 <- data.frame(tori, kion, kumo, ame)
こんな感じの架空データができる。
次に、モデルを作って各パラメータの値を得る。
model3 <- glm(tori ~ kion + kumo + ame,
family = poisson, data = d3) #GLMを作成
beta3 <- model3$coefficients #各パラメータの値を格納
18
ケース３（説明変数が複数ある場合）
•
•
•
そして、各説明変数の取りうる値すべての組み合わせで95%予測区間を・・・なん
てことはしなくてもいい。
なぜなら、各説明変数がどんな値を取ろうとも、最終的にはλという一つの値に落
ち着くから。
model$fitted でモデルが予測したλの値を取り出すことができる。
kion=21.2, kuom = 35.1, ame = 1 で、各パラメータが
beta3のときのλ・・・といった感じに、この場合20個のλを取
り出すことができる。
・そのλの最大値と最小値を調べて、その間を適当な間隔で刻む
（λs とする）。
・λs それぞれで95%予測区間を調べる。
・横軸にλ、縦軸にtoriをとってプロットする。
・さらに、横軸にλs、縦軸にそれらの予測区間をとり、描く。
次のスライドで、コード及び結果を示す。
19
ケース３（説明変数が複数ある場合）
•
Rコード
model3$fitted #モデルの予測（λ）
fit.min <- min(model3$fitted) #モデルの予測したλの最小値
fit.max <- max(model3$fitted) #モデルの予測したλの最大値
lambdas <- seq(fit.min, fit.max, 1) #λの最小値から最大値までを1間隔で刻む
pred3.1 <- qpois(0.025, lambdas)
#刻んだλごとに予測区間を調べる
pred3.2 <- qpois(0.975, lambdas)
plot(model3$fitted, tori)
lines(lambdas, pred3.1, type = "b", col = 4, pch = 20)
lines(lambdas, pred3.2, type = "b", col = 4, pch = 20)
20
ケース３（説明変数が複数ある場合）
•
結果の図
この２本の線より内側の部分が、95%予
測区間。
＜重要事項＞
・しっかりと理解していないと、途中で何
をしているのかわからなくなる。
・この予測区間は過分散のチェックにも
使える。
・もし間違ったことが書いてあったら、被
害者を増やさないためにもすぐに修正し
てほしい。
family = poisson 編 完
21
family = binomial 編
• GLM(一般化線形モデル)は、“一般化”とつくだけあって、
いろいろな確率分布を扱うことができる。
• ここでは、よく扱われる二項分布に従ったデータを題材として
解説していく・・・予定だけど、時間がなく未完成。
– 有志による付け足しを望む（あるいは暇ができたら追加する）
• 基本的な考え方は family = poisson 編と共通しているはず。
© Copyright 2024 ExpyDoc
