先日の記事はおかげさまで好評をいただいたんですが、勉強中の身で教科書を確認せずに書いたこともあり多数ツッコミをいただきました。ツッコミをいただけるというのはもちろん大変良い勉強の機会になるということで*1、今回もひとしきり勉強してみました。
ということで、自戒も込めて備忘録的に勉強したことをまとめておこうと思います。今回はあまり広く読んでもらう内容ではないので、不親切かもしれませんがごめんなさい。ただし、あまりにも理論的側面ばかり色々書いても何なので、インターネット広告業界の言葉で喩えて言うなら「クリック数*2をモデリングしたい場合」と「コンバージョン数*3をモデリングしたい場合」とに分けた、と理解してもらえたら良いかなと思ってます。
今回も参考文献は久保本です。一般化線形モデルまわりではこの本より分かりやすい本は依然としてないと思います。
データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)
- 作者: 久保拓弥
- 出版社/メーカー: 岩波書店
- 発売日: 2012/05/19
- メディア: 単行本
- 購入: 16人 クリック: 163回
- この商品を含むブログ (17件) を見る
前回記事を書いた時には手元になかったので、今回は脇に置いて調べながら書くという慎重さ(笑)。と言っても、結構肝心要の大事なところが抜けていたりして困る局面もありましたが。。。
早見表・改とモデル選択手順を作ってみた(前回のバージョンとは違います!)
前回の早見表は久保先生の資料p.17に依拠して作ったんですが、調べていくうちに実は結構問題があると思ったのでその「改」バージョンを作ってみました。
確率分布 | 特徴 |
---|---|
ポアソン分布 | データが正の離散値、平均値30ぐらいまで、標本平均=標本分散 |
負の二項分布 | データが正の離散値、平均値30ぐらいまで、標本平均<標本分散 |
二項分布 | データが離散値、ゼロ以上でしかも有限 (0, 1, 2, ... N) |
正規分布 | データが連続値もしくは離散値でも平均値が十分大*4 (-∞~∞) |
対数正規分布 | 同上、ただし正の値、範囲 (0~∞) |
ガンマ分布 | データが連続値、範囲 (0~∞) |
このように変えた理由はこの後ゆっくり説明していきます。要は「杓子定規なルールでは必ずしも決められない」ということです。これに基づくモデル選択手順は以下の通りです。
- 早見表・改で当たりをつける
- 目的変数の分布をヒストグラムから確認する
- 目的変数の平均値とサンプルサイズのデータを元に候補となる分布のシミュレーションを行ってみて、目的変数のヒストグラムと比較して、改めて最適な分布を定める
- 残差平方和、AIC、そしてover-dispersion検定を行って、最適と定めた分布が実際に最適か否かを判定する
この手順については基本的にはid:uncorrelatedさんの別の記事(銀座で働くデータサイエンティストのモデル選択について)を参考にしました。この記事でご指摘を頂いてから他の資料も当たって勉強し直してみましたが、おそらくこの方法論で概ね問題ないだろうと思ってます。
モデル選択手順について
では、具体的にその手順の各ステップについて見ていきましょう。結構長いですよー。
まず目的変数の分布を見て、最適な分布の当たりをつける
久保先生のQ&Aにもあるように、目的変数の分布を何かしらに仮定し、その上でその分布のパラメータを推定するというのが(一般化)線形モデルの基本ということなんですね。なので、何はともあれまずは目的変数のヒストグラムをプロットしてみて、どの分布に従うかを大雑把に予想する必要があります。
ちなみにid:uncorrelatedさんからのツッコミにもありましたが、リンク関数としてlogが指定されているからといってポアソン回帰と対数正規分布回帰とは同じ結果になりません。そもそものモデルの前提として、「何を推定するつもりか」が違うからです。この辺についてはid:uncorrelatedさんの記事(ポアソン回帰で推定しているモノはλの式)を読んでみて下さい。
ということで、とにかく目的分布のヒストグラムを眺めてみて、どの分布を当てはめたらベストになるかをまずは考えてみるべきだ、ということのようです。一番手っ取り早いのが、同じサンプルサイズの、同じ平均値の候補となる分布(正規分布・ポアソン分布・負の二項分布・二項分布・ガンマ分布etc.)をRなどでシミュレーションしてみて、どれくらい分布形状が一致するかor異なるかを見るという方法ですね。
サンプルサイズが小さければこれで十分にいけますが、サンプルサイズが何百万とか大きい場合は難しいかもしれません。その場合はサンプリングで対処するべきなんでしょうが、その抽出サンプルサイズは例えば{pwr}などで決定するべきだと思います。
一応残差平方和を見る
次に、モデルの適合度の目安として残差平方和を計算するというのもありだと思います。例えばですが、あるデータdに対して一般化線形モデル(正規分布線形回帰を含む)を推定したモデルのオブジェクトをd.glmとし、目的変数ベクトルをyとすれば、
> sum((y-fitted.values(d.glm))^2)
で残差平方和を求めることができます。端的に言えば、残差平方和が小さいモデルの方が当てはまりそのものとしては良いということになります。
ただし、残差平方和が小さくてもオーバーフィッティングだと意味がないので、本丸はここではなくその次のAICです。
AICを出してみる
AICは皆さんもよくご存知の情報量基準であり、AICが小さければ小さいほど予測モデルとしての当てはまりは良くなる、というのも色々な教科書に書かれている通りです。
上記の通り既にd.glmが推定済みであれば、そのAICを求めるのは
> AIC(d.glm)
とすればOKです。なので、当てはめの良し悪しを比較したいモデルのAICをズラリと並べて、その中で最もAICの小さいものを選べば良いということになります。ちなみに、同じ分布に基づくモデル群内での説明変数の組み合わせを決めたいだけならstep(){stats}やstepAIC(){MASS}でいけます。
ただし!久保先生のQ&Aページにもありますが、このAICの比較は連続分布vs.連続分布もしくは離散分布vs.離散分布という場合には問題ないんですが、連続分布vs.離散分布(例えば正規分布線形回帰モデルvs.ポアソン回帰モデル)という場合には使えません。
そういう場合は杓子定規にAICの大小で決めるだけでなく、残差平方和を比較したり、そもそも目的変数の分布形状によって定性的に評価するといった総合的な判断が必要になるようです。
二項分布・ポアソン分布・負の二項分布などの指数分布族を仮定する場合はover-dispersion検定をする
原理的な説明はweb上で参照できる日本語文献だとあまり多くないんですが、例えばこの論文に詳しく書いてあります。要は、平均が決まれば分散も自動的に定まる指数分布族を仮定する場合には、推定された以外に定数項が存在しない(切片項が0でなければいけない)という仮定を満たすかどうかが重要なわけです。
もし仮定を満たしていれば、無事ポアソン分布で表現される範囲に観測値は収まるはずです。そうでなければ、ポアソン分布では表現し切れない何かしらの過大なばらつき(過分散:over-dispersion)がある、ということになるわけです。
id:uncorrelatedさんがRで手入力でやる方法を紹介なさってますが、ちろっと調べた範囲だと生データに対して直接検定しに行く関数と、推定済みモデルに対して検定する関数とが、それぞれあるようです。qcc.overdispersion.test(){qcc}は前者のタイプで、dispersiontest(){AER}とかodTest(){pscl}*5は後者のようです。
まとめ:最終的には「どう線形モデルを一般化してあてはめたいのか」が大事
まず、扱うデータ(特に目的変数)がどんな状況のものでどんな分布に従いそうかを考えなきゃいけない、というのが今回の最大の反省のひとつです。何を当たり前のことを言ってるんだと馬鹿にされそうでアレなんですが、前回の記事の早見表を鵜呑みにしてはいけないということですね。
ツッコミにもあったように、ポアソン分布だろうが負の二項分布だろうが、平均値がでかくなれば基本的には正規分布に近い分布形状になっていきます。理由は単純で、ポアソン分布などの離散分布であっても個々の観測値は独立試行の結果なので、これの合計値は大きくなればなるほど中心極限定理に従って正規分布に漸近していくということです。なので、総click数とか総PVとか総impression数とか大きな数字を相手にする場合のように、平均値がかなり大きいようであればあっさり正規分布に従うと仮定してlm(){base}でやってしまって多分問題ないだろう、ということになると思います。
ところが、総CV数のような場合によっては(time binによっては)1ケタとかものすごく小さい値になるような事例では、ポアソン回帰や負の二項分布回帰が有効になるのでしょう。実際、マーケティングリサーチのパネル調査データなどでサンプルサイズが5000ぐらいまでとかだと、そういうことはザラにあるので。
確認のため、適当にrnorm(), rpois()でシミュレーションしてみるとこんな感じになります。まず、平均値が3の場合。
次に平均値が30の場合。
最後に平均値が5000の場合。
確かに、平均値が大きくなるにつれて正規分布とポアソン分布の形状が互いに似てくることが分かります。「連続値か離散値か」という分類も、ここまで来るとあまり大きな効果はなさそうです。
ちなみに平均値もサンプルサイズも大きくて、正規分布で大丈夫だろうと思われる場合は、初めからKolmogorov-Smirnov(コロモゴロフ・スミノフ)検定で目的変数が正規分布に従っているかどうか検定してしまうという方法もアリだと思います。
> ks.test(x,"pnorm",mean=mean(x),sd=sd(x))
ということで、ある意味当たり前なんですが「目的変数の平均値が4ケタ5ケタと大きければ正規分布線形回帰で問題ないが、1ケタ2ケタと小さい時はポアソン回帰や負の二項分布回帰も疑ってかかるべき」という結論になります。
ところで、この結論が実は問題なんですね。これって、久保先生が提唱している早見表とは必ずしも合致しないのです(冒頭の早見表・改を参照のこと)。元々の久保先生の早見表だと「カウントデータならばポアソン分布・負の二項分布・二項分布」という説明になっていますが、カウントデータであっても、平均値とサンプルサイズの双方が大きくなれば、基本的には中心極限定理に従って正規分布に漸近していきますし、実用上はそう仮定してもほとんど差し支えありません。
なので、カウントデータ(離散値)であっても平均値が大きいのであれば正規分布を選択した方が良いし、そうでない典型的なポアソン分布に従うようなデータであればそのように扱うべきだ、ということになります。
改めて手順を書くと、
- 早見表・改で当たりをつける
- 目的変数の分布をヒストグラムから確認する
- 目的変数の平均値とサンプルサイズのデータを元に候補となる分布のシミュレーションを行ってみて、目的変数のヒストグラムと比較して、改めて最適な分布を定める
- 残差平方和、AIC、そしてover-dispersion検定を行って、最適と定めた分布が実際に最適か否かを判定する
そうそう、正規分布線形回帰をやる際の留意点。Rで計算すると、lm(){base}はOLSを用いるんですが、glm(){base}では最尤法を用いることになるようです(正規分布ならfamily=gaussianとすることでいけます)。計算負荷の都合を考えると、純粋に正規分布仮定が成立する時はlm()の方がコストが低くて適している模様です。ただし、厳密にはOLSには色々な制約条件があるため*6、どうしても気になる場合はglm()で最尤法を用いた方が無難なケースもあるかもです。たかだか数百ぐらいのサンプルサイズならどちらでやっても大差ないと思いますが、数万ぐらいになるとさすがに最尤法だと遅いかもしれません。。。
それでもダメだったら
要するにどのように分布をいじってもover-dispersionを解決できないくらい個体差がでかすぎるということなので、混合効果モデル(GLMM)を使うべきだということですね。これはもう仕方ないので、{glmmML}とか{R2WinBUGS}とか{MCMCpack}の出番になります。
もっともそこまで行くと今回の記事の趣旨からだいぶ逸脱してしまいますので、この辺の話題は今回は一旦置いときます。また実務で使う機会が増えて勉強の必要に迫られたら記事にしようかなー、と。最近データ分析業界でも個体差の大きいデータを扱うところが増えていて、ちょくちょくGLMM、MCMCそして階層ベイズの話題が持ち上がることがあります。いずれみっちりやることになるんだろうなぁと思ってます。。。
サンプルデータで実践してみる
ということで、改めてやってみようと思います。Rでやるだけですが、一応{MASS}, {qcc}, {AER}, {pscl}をインストールして展開しておいて下さい。
正規分布が当てはまりそうなケース
サンプルデータ1をGitHubに上げてあります。これはクリック数をモデリングしたいケースを想定しています。落としてきて、dat1という名前でRに読み込んでおきましょう。
まず、平均と分散・標準偏差を見てみます。
> mean(dat1$click) [1] 49977.2 > var(dat1$click) [1] 43304.16 > sd(dat1$click) [1] 208.0965
何とも言いようがない感じですが、平均が50000ぐらいと非常に大きく、サンプルサイズも70あるので多分正規分布線形モデルで大丈夫でしょう。それではやってみます。
> dat1.lm<-lm(click~.,dat1) # OLSでやってみる > dat1.lmg<-glm(click~.,dat1,family=gaussian) # 念のため最尤法でもやっておく > summary(dat1.lm) Call: lm(formula = click ~ ., data = dat1) Residuals: Min 1Q Median 3Q Max -337.52 -89.91 -0.11 106.76 300.47 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 49848.856 96.564 516.224 < 2e-16 *** d1 230.964 48.901 4.723 1.56e-05 *** d2 -80.567 50.037 -1.610 0.1129 d3 41.770 41.550 1.005 0.3190 d4 31.553 38.441 0.821 0.4152 d5 -35.350 44.419 -0.796 0.4294 d6 72.892 50.037 1.457 0.1507 d7 4.004 40.195 0.100 0.9210 d8 -22.608 42.271 -0.535 0.5948 d9 97.991 50.498 1.941 0.0573 . d10 -43.447 44.629 -0.974 0.3344 d11 -30.645 45.615 -0.672 0.5044 d12 25.072 40.450 0.620 0.5378 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 146.1 on 57 degrees of freedom Multiple R-squared: 0.593, Adjusted R-squared: 0.5073 F-statistic: 6.921 on 12 and 57 DF, p-value: 1.588e-07 > summary(dat1.lmg) Call: glm(formula = click ~ ., family = gaussian, data = dat1) Deviance Residuals: Min 1Q Median 3Q Max -337.52 -89.91 -0.11 106.76 300.47 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 49848.856 96.564 516.224 < 2e-16 *** d1 230.964 48.901 4.723 1.56e-05 *** d2 -80.567 50.037 -1.610 0.1129 d3 41.770 41.550 1.005 0.3190 d4 31.553 38.441 0.821 0.4152 d5 -35.350 44.419 -0.796 0.4294 d6 72.892 50.037 1.457 0.1507 d7 4.004 40.195 0.100 0.9210 d8 -22.608 42.271 -0.535 0.5948 d9 97.991 50.498 1.941 0.0573 . d10 -43.447 44.629 -0.974 0.3344 d11 -30.645 45.615 -0.672 0.5044 d12 25.072 40.450 0.620 0.5378 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for gaussian family taken to be 21335.18) Null deviance: 2987987 on 69 degrees of freedom Residual deviance: 1216105 on 57 degrees of freedom AIC: 910.04 Number of Fisher Scoring iterations: 2 # 最尤法でも特に推定結果は変わらない
そんなに問題はなさそうです。なお、step()で説明変数を絞り込むとこんな感じになります。
> step(dat1.lm) # 中略 Step: AIC=699.22 click ~ d1 + d2 + d6 + d9 Df Sum of Sq RSS AIC <none> 1321709 699.22 - d6 1 48206 1369915 699.72 - d2 1 55684 1377393 700.10 - d9 1 141383 1463092 704.33 - d1 1 467889 1789598 718.43 Call: lm(formula = click ~ d1 + d2 + d6 + d9, data = dat1) Coefficients: (Intercept) d1 d2 d6 d9 49824.84 210.31 -74.21 61.44 109.61 > dat1.lm2<-lm(formula = click ~ d1 + d2 + d6 + d9, data = dat1) > summary(dat1.lm2) Call: lm(formula = click ~ d1 + d2 + d6 + d9, data = dat1) Residuals: Min 1Q Median 3Q Max -376.94 -93.50 -14.26 99.95 303.01 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 49824.84 57.59 865.100 < 2e-16 *** d1 210.31 43.84 4.797 9.79e-06 *** d2 -74.21 44.84 -1.655 0.1028 d6 61.44 39.90 1.540 0.1285 d9 109.61 41.57 2.637 0.0105 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 142.6 on 65 degrees of freedom Multiple R-squared: 0.5577, Adjusted R-squared: 0.5304 F-statistic: 20.49 on 4 and 65 DF, p-value: 5.871e-11
ちなみに残差平方和をチェックすると、
> sum((y-fitted.values(dat1.lm))^2) [1] 174791188668 > sum((y-fitted.values(dat1.lm2))^2) [1] 174791086420
となって、説明変数を絞り込んだ方がAICで見ても残差平方和で見ても、予測モデルとしての当てはまりは良いという結論に達します。変数重要度としてはd1とd9が重要らしいな、ということで大体OKでしょう。
以下はおまけ。決定木とランダムフォレストで回帰をかけた結果です。
> dat1.rp<-rpart(click~.,dat1) > plot(dat1.rp,uniform=T,margin=0.2) > text(dat1.rp,uniform=T,use.n=F,all=F) > print(dat1.rp) n= 70 node), split, n, deviance, yval * denotes terminal node 1) root 70 2987987.00 49977.20 2) d1< 0.5 28 532000.70 49809.21 4) d8>=0.5 22 395160.40 49790.73 8) d9< 0.5 18 321428.90 49772.06 16) d12< 0.5 10 134969.60 49748.20 * 17) d12>=0.5 8 173654.90 49801.88 34) d6>=0.5 6 109034.80 49750.17 68) d4>=0.5 3 50078.00 49669.00 * 69) d4< 0.5 3 19428.67 49831.33 * 35) d6< 0.5 2 450.00 49957.00 * 9) d9>=0.5 4 39216.75 49874.75 * 5) d8< 0.5 6 101752.00 49877.00 10) d9>=0.5 3 56114.00 49804.00 * 11) d9< 0.5 3 13664.00 49950.00 * 3) d1>=0.5 42 1139090.00 50089.19 6) d9< 0.5 11 301856.70 49946.45 12) d4< 0.5 4 83355.00 49814.50 * 13) d4>=0.5 7 109054.90 50021.86 26) d3< 0.5 3 13712.67 49917.33 * 27) d3>=0.5 4 37984.75 50100.25 * 7) d9>=0.5 31 533602.20 50139.84 14) d11>=0.5 27 385837.90 50122.93 28) d10>=0.5 23 208232.90 50100.30 56) d8< 0.5 14 70951.43 50070.57 * 57) d8>=0.5 9 105652.20 50146.56 114) d4>=0.5 6 32666.83 50098.83 * 115) d4< 0.5 3 31992.00 50242.00 * 29) d10< 0.5 4 98158.00 50253.00 * 15) d11< 0.5 4 87910.00 50254.00 * > plotcp(dat1.rp) > dat1.rp2<-rpart(click~.,dat1,cp=0.21) > plot(dat1.rp2,uniform=T,margin=0.2) > text(dat1.rp2,uniform=T,use.n=F,all=F) > print(dat1.rp2) n= 70 node), split, n, deviance, yval * denotes terminal node 1) root 70 2987987.0 49977.20 2) d1< 0.5 28 532000.7 49809.21 * 3) d1>=0.5 42 1139090.0 50089.19 * > tuneRF(dat1[,-13],dat[,13],doBest=T) mtry = 4 OOB error = 187657.3 Searching left ... mtry = 2 OOB error = 191852 -0.02235303 0.05 Searching right ... mtry = 8 OOB error = 209201.1 -0.1148043 0.05 Call: randomForest(x = x, y = y, mtry = res[which.min(res[, 2]), 1]) Type of random forest: regression Number of trees: 500 No. of variables tried at each split: 4 Mean of squared residuals: 182209 % Var explained: 40.29 > dat1.rf<-randomForest(click~.,dat1,mtry=4) > print(dat1.rf) Call: randomForest(formula = click ~ ., data = dat1, mtry = 4) Type of random forest: regression Number of trees: 500 No. of variables tried at each split: 4 Mean of squared residuals: 25730.91 % Var explained: 39.72 > importance(dat1.rf) IncNodePurity d1 699095.07 d2 299057.38 d3 108777.38 d4 111309.79 d5 107866.45 d6 103819.53 d7 94657.48 d8 113031.15 d9 462585.78 d10 118602.27 d11 121147.19 d12 84900.28
機械学習で見ても、d1とd9が重要らしいということが分かります。最後に、コルモゴロフ・スミノフ検定で目的変数clickが正規分布に綺麗に従うかどうかのダメ押しをしておきます。
> ks.test(dat1$click,"pnorm",mean=mean(dat1$click),sd=sd(dat1$click)) One-sample Kolmogorov-Smirnov test data: dat1$click D = 0.09, p-value = 0.6224 # 正規分布だよ、という結果 alternative hypothesis: two-sided 警告メッセージ: In ks.test(dat1$click, "pnorm", mean = mean(dat1$click), sd = sd(dat1$click)) : コルモゴロフ・スミノフ検定において、タイは現れるべきではありません >hist(dat1$click,breaks=20)
ということで、dat1に関して言えば(たぶん)正規分布線形回帰モデルで問題なかったのではないか、ということになりました。
ポアソン分布が当てはまりそうなケース
サンプルデータ2をGitHubに上げてあります。これはコンバージョン数をモデリングしたいケースを想定しています。落としてきて、dat2という名前でRに読み込んでおきましょう。
> mean(dat2$cv) [1] 7.285714 > var(dat2$cv) [1] 8.207039 > sd(dat2$cv) [1] 2.864793
目的変数の平均が30以下で、平均と分散がほぼ同じ。ということはポアソン分布が当てはまるのではないかと考えられます。ヒストグラムを描いてもっと細かく見てみましょう。
> hist(dat2$cv,breaks=20) # 生データのヒストグラム > hist(rnorm(70,mean=mean(dat2$cv)),breaks=20) # 正規分布でシミュレーション > hist(rpois(70,lambda=mean(dat2$cv)),breaks=20) # ポアソン分布でシミュレーション
これはポアソン分布が妥当っぽいですね。ということで、まずはポアソン回帰モデルを当てはめてみることにします。ただし比較のための当て馬として正規分布線形回帰モデルも当てはめておきます。
> dat2.lm<-lm(cv~.,dat2) # 正規分布線形回帰モデル(当て馬) > dat2.glm<-glm(cv~.,dat2,family=poisson) # ポアソン回帰モデル > summary(dat2.lm) Call: lm(formula = cv ~ ., data = dat2) Residuals: Min 1Q Median 3Q Max -4.2030 -1.3532 -0.3425 0.9910 5.1315 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 5.85072 1.38622 4.221 8.84e-05 *** d1 3.28149 0.70200 4.675 1.85e-05 *** d2 -1.03809 0.71830 -1.445 0.1539 d3 0.56057 0.59646 0.940 0.3513 d4 0.46381 0.55184 0.840 0.4042 d5 -0.66893 0.63765 -1.049 0.2986 d6 1.26401 0.71830 1.760 0.0838 . d7 -0.01514 0.57701 -0.026 0.9792 d8 -0.47852 0.60682 -0.789 0.4336 d9 1.12945 0.72491 1.558 0.1248 d10 -0.77967 0.64066 -1.217 0.2286 d11 -0.42796 0.65482 -0.654 0.5160 d12 0.36019 0.58067 0.620 0.5375 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 2.097 on 57 degrees of freedom Multiple R-squared: 0.5574, Adjusted R-squared: 0.4643 F-statistic: 5.983 on 12 and 57 DF, p-value: 1.293e-06 > summary(dat2.glm) Call: glm(formula = cv ~ ., family = poisson, data = dat2) Deviance Residuals: Min 1Q Median 3Q Max -1.8245 -0.4835 -0.1034 0.5031 1.4941 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 1.667603 0.250712 6.651 2.9e-11 *** d1 0.487817 0.130344 3.743 0.000182 *** d2 -0.104083 0.120958 -0.860 0.389523 d3 0.092941 0.111606 0.833 0.404977 d4 0.060384 0.096983 0.623 0.533532 d5 -0.067516 0.113376 -0.596 0.551507 d6 0.152919 0.121863 1.255 0.209534 d7 -0.002984 0.104275 -0.029 0.977170 d8 -0.060244 0.103727 -0.581 0.561381 d9 0.174512 0.133993 1.302 0.192783 d10 -0.104137 0.114455 -0.910 0.362903 d11 -0.071120 0.118223 -0.602 0.547458 d12 0.050842 0.099898 0.509 0.610791 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 78.643 on 69 degrees of freedom Residual deviance: 34.416 on 57 degrees of freedom AIC: 324.06 Number of Fisher Scoring iterations: 4
当て馬として正規分布線形回帰モデルをぶつけてみましたが、残差平方和を比べてみると
> sum((y-fitted.values(dat2.lm))^2) [1] 250.6107 # 正規分布線形回帰モデル > sum((y-fitted.values(dat2.glm))^2) [1] 236.3826 # ポアソン回帰モデル
ということで、ポアソン回帰モデルの方が分布形状から大ざっぱに推測した通り、当てはまりが良いようです。念のため負の二項分布回帰モデルでも当てはめてみます。
> dat2.glmnb<-glm.nb(cv~.,dat2) 警告メッセージ: 1: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached 2: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached > summary(dat2.glmnb) Call: glm.nb(formula = cv ~ ., data = dat2, init.theta = 244630.2469, link = log) Deviance Residuals: Min 1Q Median 3Q Max -1.8245 -0.4835 -0.1034 0.5031 1.4941 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 1.667603 0.250716 6.651 2.9e-11 *** d1 0.487816 0.130346 3.742 0.000182 *** d2 -0.104083 0.120960 -0.860 0.389528 d3 0.092941 0.111607 0.833 0.404984 d4 0.060384 0.096984 0.623 0.533540 d5 -0.067516 0.113378 -0.595 0.551516 d6 0.152918 0.121865 1.255 0.209545 d7 -0.002984 0.104276 -0.029 0.977173 d8 -0.060243 0.103729 -0.581 0.561390 d9 0.174511 0.133995 1.302 0.192791 d10 -0.104136 0.114457 -0.910 0.362917 d11 -0.071119 0.118225 -0.602 0.547469 d12 0.050842 0.099899 0.509 0.610796 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for Negative Binomial(244630.2) family taken to be 1) Null deviance: 78.641 on 69 degrees of freedom Residual deviance: 34.415 on 57 degrees of freedom AIC: 326.06 Number of Fisher Scoring iterations: 1 Theta: 244630 Std. Err.: 3981615 Warning while fitting theta: iteration limit reached 2 x log-likelihood: -298.058 > sum((y-fitted.values(dat2.glmnb))^2) [1] 236.3828 # 負の二項分布回帰モデルの残差平方和
当てはまりが良過ぎて警告が出ました(笑)。ただし、残差平方和とAICで見るとほんの僅かですがポアソン回帰モデルの方が良いようです。最後に、over-dispersion検定をかけてみましょう。
> odTest(dat2.glmnb) Likelihood ratio test of H0: Poisson, as restricted NB model: n.b., the distribution of the test-statistic under H0 is non-standard e.g., see help(odTest) for details/references Critical value of test statistic at the alpha= 0.05 level: 2.7055 Chi-Square Test Statistic = -0.0011 p-value = 0.5 # 負の二項分布回帰モデルでover-dispersionなし > dispersiontest(dat2.glm) Overdispersion test data: dat2.glm z = -6.7175, p-value = 1 # over-dispersionなし alternative hypothesis: true dispersion is greater than 1 sample estimates: dispersion 0.4793123 > dispersiontest(dat2.glm,alternative="two.sided") Dispersion test data: dat2.glm z = -6.7175, p-value = 1.849e-11 # 過小分散だよという検定結果 alternative hypothesis: true dispersion is not equal to 1 sample estimates: dispersion 0.4793123 > qcc.overdispersion.test(dat2$cv) Overdispersion test Obs.Var/Theor.Var Statistic p-value poisson data 1.126456 77.72549 0.22071
過小分散の疑いがあるんですが、久保本にも解決策は載ってないので今回は無視します、すみません。最後に、説明変数を絞り込みたいのでstepAIC(){MASS}にかけておきます。
> stepAIC(dat2.glm) # 中略 Step: AIC=309.94 cv ~ d1 + d9 Df Deviance AIC <none> 40.294 309.94 - d9 1 45.744 313.39 - d1 1 55.883 323.52 Call: glm(formula = cv ~ d1 + d9, family = poisson, data = dat2) Coefficients: (Intercept) d1 d9 1.5579 0.4297 0.2404 Degrees of Freedom: 69 Total (i.e. Null); 67 Residual Null Deviance: 78.64 Residual Deviance: 40.29 AIC: 309.9 > dat2.glm2<-glm(formula = cv ~ d1 + d9, family = poisson, data = dat2) > summary(dat2.glm2) Call: glm(formula = cv ~ d1 + d9, family = poisson, data = dat2) Deviance Residuals: Min 1Q Median 3Q Max -1.8061 -0.4361 -0.1019 0.5505 1.9973 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 1.55788 0.08945 17.417 < 2e-16 *** d1 0.42970 0.11088 3.875 0.000107 *** d9 0.24043 0.10400 2.312 0.020786 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for poisson family taken to be 1) Null deviance: 78.643 on 69 degrees of freedom Residual deviance: 40.294 on 67 degrees of freedom AIC: 309.94 Number of Fisher Scoring iterations: 4
d1とd9が重要らしい、という同じ結論になりました。最後におまけでランダムフォレストでもチェックしてみます。
> dat2.rf<-randomForest(cv~.,dat2) > print(dat2.rf) Call: randomForest(formula = cv ~ ., data = dat2) Type of random forest: regression Number of trees: 500 No. of variables tried at each split: 4 Mean of squared residuals: 5.280869 % Var explained: 34.72 > importance(dat2.rf) IncNodePurity d1 131.98383 d2 54.74000 d3 21.09103 d4 19.17511 d5 22.73686 d6 26.62820 d7 15.26824 d8 26.49701 d9 76.65858 d10 28.83372 d11 19.89582 d12 17.11731
ランダムフォレストだと分散の35%ぐらいしか説明できていないので、普通にポアソン回帰モデルの方が良さそうですね。
ともあれ、dat2についてはポアソン回帰モデルdat2.glmがベストらしい、という結論になりました。
二項分布(ロジスティック回帰)が当てはまるケース
そもそも目的変数が二値データやカウントデータかつ上限値があるような生存確率系データであれば、基本的にはロジスティック回帰の他に方法を選びようがないので、サンプルデータでの演習はここでは割愛します。
ただし、二項分布もポアソン分布や負の二項分布と同様に平均が決まれば分散も決定してしまう指数分布族なので、同様にover-dispersionの問題が生じ得る*7点に注意が必要です。その場合は、前述のようにGLMMなどへと発展していくことになります。
最後に種明かし
実はdat1$clickもdat2$cvもポアソン分布乱数から発生させたデータをソートしただけの代物です*8。つまり今回の結論はものすごく当然だったりします。
なのですが、平均値の大きい前者は正規分布線形回帰モデルで問題なく、一方で平均値の小さい後者がポアソン回帰モデルが妥当というように、結果がそれなりに分かれたのは興味深かったです。
ちなみに、前回の記事で用いたサンプルデータのclickは実は線形関数に正規乱数を足し算して平均を数千のオーダーにしたものという、玉虫色のデータだったのでした。そりゃあどれもうまく当てはまらないよなぁ、と。。。
最後に、id:uncorrelatedさんのツッコミは大変勉強になりました。有難うございました。来たるGLMM編・階層ベイズ編でも突っ込んで頂ければ幸いです*9。
おまけ
先に予言しておきますが、多分このエントリのはてブは伸びないと思います(笑)。
そして、陽にモデル式を出す必要がなければ、普通に機械学習系の手法使いたいんですけどね。。。世の中「Excelでも実装可能な陽に表せる回帰モデル式」の納品を求められるケースが結構多いので、SVMとかRVMとかランダムフォレストとかでは納品できないという辛さ。
(※この後の部分は妥当ではないという指摘があり、自分自身で検討してみても確信が持てなかったので削除しました。ごめんなさい)
追記
別の角度からのツッコミをいただきました。
統計解析は分野ごとに風習があって、GLMが良く使われるとGLMを中心に話が進むのだと思うのだが、一旦GLMから離れて推定モデルを見た方がより一般的な議論になると思える。あまり考えないでガツガツと分析結果を出す方が生産的だけど(´▽’)アッハン
ということで、単に僕が以前からなじみ深いGLM*10をメインに使っているだけ&あまり考えないでガツガツと分析結果を出さなきゃいけない状況に来ているということで、この辺の議論はまたの機会にします。。。MCMCとか扱う段階が来たら振り返りますので。
*1:裏先生の回 http://tjo.hatenablog.com/entry/2013/08/15/001338も参照のこと
*3:データセットごとだとほんの数件ぐらいというオーダーに留まりやすい
*5:ただしglm.nb(){MASS}で推定した負の二項分布回帰モデル限定
*6:例えばこの辺など http://aoki2.si.gunma-u.ac.jp/lecture/mb-arc/arc013/168.html の192
*7:久保本p.148で触れられています
*8:もっと言ってしまうと説明変数はどちらも完全に同じ(笑)
*9:他力本願モード(笑)
*10:研究者だった頃の常套手法:ただし全部ツールに丸投げorz