六本木で働くデータサイエンティストのブログ

元祖「銀座で働くデータサイエンティスト」です / 道玄坂→銀座→東京→六本木

「使い分け」ではなく「妥当かどうか」が大事:重回帰分析&一般化線形モデル選択まわりの再まとめ

先日の記事はおかげさまで好評をいただいたんですが、勉強中の身で教科書を確認せずに書いたこともあり多数ツッコミをいただきました。ツッコミをいただけるというのはもちろん大変良い勉強の機会になるということで*1、今回もひとしきり勉強してみました。


ということで、自戒も込めて備忘録的に勉強したことをまとめておこうと思います。今回はあまり広く読んでもらう内容ではないので、不親切かもしれませんがごめんなさい。ただし、あまりにも理論的側面ばかり色々書いても何なので、インターネット広告業界の言葉で喩えて言うなら「クリック数*2をモデリングしたい場合」と「コンバージョン数*3をモデリングしたい場合」とに分けた、と理解してもらえたら良いかなと思ってます。


今回も参考文献は久保本です。一般化線形モデルまわりではこの本より分かりやすい本は依然としてないと思います。



前回記事を書いた時には手元になかったので、今回は脇に置いて調べながら書くという慎重さ(笑)。と言っても、結構肝心要の大事なところが抜けていたりして困る局面もありましたが。。。


今回参考にした資料など


前回の記事の最後にも列挙しましたが、改めて。


ぶっちゃけこちらを読んでもらえれば今回の僕の記事なんざ不要だと思いますが(笑)、しばしお付き合いのほどを。


早見表・改とモデル選択手順を作ってみた(前回のバージョンとは違います!)


前回の早見表は久保先生の資料p.17に依拠して作ったんですが、調べていくうちに実は結構問題があると思ったのでその「改」バージョンを作ってみました。

確率分布 特徴
ポアソン分布 データが正の離散値、平均値30ぐらいまで、標本平均=標本分散
負の二項分布 データが正の離散値、平均値30ぐらいまで、標本平均<標本分散
二項分布 データが離散値、ゼロ以上でしかも有限 (0, 1, 2, ... N)
正規分布 データが連続値もしくは離散値でも平均値が十分大*4 (-∞~∞)
対数正規分布 同上、ただし正の値、範囲 (0~∞)
ガンマ分布 データが連続値、範囲 (0~∞)


このように変えた理由はこの後ゆっくり説明していきます。要は「杓子定規なルールでは必ずしも決められない」ということです。これに基づくモデル選択手順は以下の通りです。

  1. 早見表・改で当たりをつける
  2. 目的変数の分布をヒストグラムから確認する
  3. 目的変数の平均値とサンプルサイズのデータを元に候補となる分布のシミュレーションを行ってみて、目的変数のヒストグラムと比較して、改めて最適な分布を定める
  4. 残差平方和、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上で参照できる日本語文献だとあまり多くないんですが、例えばこの論文に詳しく書いてあります。要は、平均が決まれば分散も自動的に定まる指数分布族を仮定する場合には、推定された\beta_i以外に定数項が存在しない(切片項が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の場合。


f:id:TJO:20130920131114p:plain


次に平均値が30の場合。


f:id:TJO:20130920131136p:plain


最後に平均値が5000の場合。


f:id:TJO:20130920131154p:plain


確かに、平均値が大きくなるにつれて正規分布ポアソン分布の形状が互いに似てくることが分かります。「連続値か離散値か」という分類も、ここまで来るとあまり大きな効果はなさそうです。


ちなみに平均値もサンプルサイズも大きくて、正規分布で大丈夫だろうと思われる場合は、初めからKolmogorov-Smirnov(コロモゴロフ・スミノフ)検定で目的変数が正規分布に従っているかどうか検定してしまうという方法もアリだと思います。

> ks.test(x,"pnorm",mean=mean(x),sd=sd(x))


ということで、ある意味当たり前なんですが「目的変数の平均値が4ケタ5ケタと大きければ正規分布線形回帰で問題ないが、1ケタ2ケタと小さい時はポアソン回帰や負の二項分布回帰も疑ってかかるべき」という結論になります。


ところで、この結論が実は問題なんですね。これって、久保先生が提唱している早見表とは必ずしも合致しないのです(冒頭の早見表・改を参照のこと)。元々の久保先生の早見表だと「カウントデータならばポアソン分布・負の二項分布・二項分布」という説明になっていますが、カウントデータであっても、平均値とサンプルサイズの双方が大きくなれば、基本的には中心極限定理に従って正規分布に漸近していきますし、実用上はそう仮定してもほとんど差し支えありません。


なので、カウントデータ(離散値)であっても平均値が大きいのであれば正規分布を選択した方が良いし、そうでない典型的なポアソン分布に従うようなデータであればそのように扱うべきだ、ということになります。


改めて手順を書くと、

  1. 早見表・改で当たりをつける
  2. 目的変数の分布をヒストグラムから確認する
  3. 目的変数の平均値とサンプルサイズのデータを元に候補となる分布のシミュレーションを行ってみて、目的変数のヒストグラムと比較して、改めて最適な分布を定める
  4. 残差平方和、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

f:id:TJO:20130922231902p:plain

f:id:TJO:20130922231912p:plain

f:id:TJO:20130922231921p:plain


機械学習で見ても、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)

f:id:TJO:20130922232820p:plain


ということで、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)
# ポアソン分布でシミュレーション

f:id:TJO:20130922232831p:plain

f:id:TJO:20130922233531p:plain

f:id:TJO:20130922233540p:plain


これはポアソン分布が妥当っぽいですね。ということで、まずはポアソン回帰モデルを当てはめてみることにします。ただし比較のための当て馬として正規分布線形回帰モデルも当てはめておきます。

> 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も参照のこと

*2:データセットごとに何千何万というオーダーになりやすい

*3:データセットごとだとほんの数件ぐらいというオーダーに留まりやすい

*4:中心極限定理が効き始める程度に大きいという意味

*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