GLMの予測区間 ～直感的理解と図示方法～ var. 2010.3.24

目次 はじめに family = poisson 編 family = binomial 編（以下、未完成・・・） ケース1（説明変数がない場合） ケース2（説明変数が1つある場合） ケース3（説明変数が複数ある場合） family = binomial 編（以下、未完成・・・）

はじめに このスライドの目的 タイトルの通り、GLMの予測区間について理解を深め、 Rによる作図ができるようになること。 前提とするもの よくつかわれる確率分布に関する、ある程度の知識 GLM(一般化線形モデル)に関する、ある程度の知識 予測区間とは 母集団を仮定した上で、将来観察されるであろう標本値に対して「どの範囲にあると予測されるか」を示すものである。（wikipedia より） これに対し、信頼区間とは、母集団の母数（標本から測定できない）に対して 「どの範囲にあると推定できるか」を示すものである。混同しないように注意。（wikipedia より） その他 “GLMの”予測区間と言ってるけど、考え方は他のモデルでも同じかも。 ただし、ランダム効果入りのモデルは別。要注意。

family = poisson 編 GLM(一般化線形モデル)は、“一般化”とつくだけあって、 いろいろな確率分布を扱うことができる。 ここでは、その中でも理解しやすいと思われる、ポアソン分布に従ったデータを題材として解説していく。

ケース１（説明変数がない場合） まさお君は釣りが好きであった。今までに手賀沼に釣りに行った回数はちょうど50回である。まさお君は、よく釣れるブルーギルの数を記録していた。そのデータを利用して、手賀沼で釣れるブルーギルの数をモデリングすることにした。ちなみに、釣りをする時間は毎回一定であった。 釣れたブルーギル数のヒストグラム gill <- rpois(50, 5)　#ポアソン分布に従う乱数生成 hist(gill) #ヒストグラム作成 mean(gill) #平均 var(gill) #分散 平均：4.74 分散：4.52 平均≒分散となっていて、ポアソン分布に従っていそう！ 注）乱数を使っているので、実際に試した場合とは微妙に値が異なると思います。

ケース１（説明変数がない場合） 釣れたブルーギルの匹数を応答変数として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の係数）。 beta2には、切片及びpeopleの係数を格納。 モデルから得られたパラメータを使って、横軸（people）がとる値ごとに95%予測区間を調べていく。 model2 <- glm(bass ~ people, family = poisson, data = d) #GLMを作成 beta2 <- model1$coefficients #このモデルのパラメータを格納 #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 編と共通しているはず。