今さら人に聞けない「重回帰分析の各手法の使い分け」 - 渋谷駅前で働くデータサイエンティストのブログ

（※※※続編記事書きました→「使い分け」ではなく「妥当かどうか」が大事：重回帰分析＆一般化線形モデル選択まわりの再まとめ）

今ちょうどadtech tokyo 2013の会期中で、職場からも近い＆会社から行ってこいという指示が出たということで僕も色々セッションを聞いたり企業ブースのお話を聞いたりしてる*1ところです。

ところで、いくつかのセッションの中でキーワードとして「重回帰分析」という言葉が出てきてました。ま、それ自体はこのブログでもRによるデータ分析絡みで頻出だし、ぶっちゃけありふれた手法と言って良いでしょう。やりようによっては普通にExcelでもできますし、それだけ人口に膾炙していると言って良いのかもですね。

ただし。意外にも内部のパラメータというか細かい手法の分岐というか、それこそ普通の線形モデルvs.一般化線形モデル（バリエーション多数）があることを無視して漫然と重回帰分析をやると、結果が変わってしまったり間違った答えにつながってしまったりするんですね。。。ということで、今回は軽めのネタとして「今さら人に聞けない『重回帰分析の各手法の使い分け』」についてRによる実践例を交えながら書いてみようと思います。要は先日のポアソン回帰モデルの記事の続きです（笑）。

参考図書はこちら。前回の記事でも参考にした北大の久保先生の名著と、定番の古典『統計学入門』（通称赤本）です。

データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)

作者: 久保拓弥
出版社/メーカー: 岩波書店
発売日: 2012/05/19
メディア: 単行本
購入: 16人クリック: 163回
この商品を含むブログ (17件) を見る

統計学入門 (基礎統計学)

作者: 東京大学教養学部統計学教室
出版社/メーカー: 東京大学出版会
発売日: 1991/07/09
メディア: 単行本
購入: 158人クリック: 3,604回
この商品を含むブログ (77件) を見る

はっきり言って今回の記事なんか読まなくても、この久保先生の本と赤本を読めば事足りるんじゃないかと思うんですが*2、もっと手っ取り早く！という人のために適当にまとめてみました。

（※長年のとんでもない勘違いを含めて一部修正してあります。多分この後もさらに修正が入ります。。。僕もまだまだ勉強しないと。。。）

（※※いただいたツッコミによると、今回のサンプルデータに対するモデル選択には色々とややこしい点があるようです。そして最適モデル選択のプロセスにもかなり抜けている部分があると分かりました。いずれ問題点をまとめて改めてエントリを立てて論じたいと思いますので、悪しからずご了承を。。。）

必要なRパッケージ

{MASS}をインストールして展開して下さい（興味のある人は{betareg}もインストールして試してみても良いかも）。他は今回は特に要りません。

とりあえずこの早見表に従おう

以前の記事でも参照したページ（GLM入門編1 - Murakami's Memorandum）＆久保先生の資料のp.17に倣うと、こんな感じです。

確率分布	特徴
ポアソン分布	データが離散値、ゼロ以上、上限なし、標本平均＝標本分散
負の二項分布	データが離散値、ゼロ以上、上限なし、標本平均＜標本分散
二項分布	データが離散値、ゼロ以上で有限 (0, 1, 2, ... N)
正規分布	データが連続値、範囲 (-∞～∞)
対数正規分布	データが連続値、範囲 (0～∞)
ガンマ分布	データが連続値、範囲 (0～∞)
ベータ分布	データが連続値、範囲 (0～1)

ぶっちゃけこれでおしまいです（笑）。この早見表に従って、Rで言えばlm(){base}を使うか、glm(){base}の引数familyを指定するか、さもなくばglm.nb(){MASS}を使うかを決めれば良いだけです。それぞれの分布の詳細についてはいちいち挙げていくと面倒なので、『統計学入門』あたりを参照してもらえれば良いかと。

glm(){base}の引数familyについてCRANのヘルプを参照すると、以下のように書かれています。

Usage

family(object, ...)

binomial(link = "logit")
gaussian(link = "identity")
Gamma(link = "inverse")
inverse.gaussian(link = "1/mu^2")
poisson(link = "log")

今回はquasi-については割愛しました。負の二項分布の場合はglm.nb(){MASS}でいけます。またベータ分布の場合はbetareg.fit(){betareg}を使えばOK。ともあれ、こういう感じでどの分布に基づいて一般化線形モデルの計算をするかをRでは指定することができます。

試しにやってみると？

さらに簡単のため、普通の線形回帰と3つの一般化線形モデル（ポアソン分布・負の二項分布・二項分布＝ロジスティック回帰）の4パターンだけを今回は取り上げます。

サンプルデータはRではデフォルトで入っているairqualityと、前回の記事でも使ったサンプルデータ（dat1という名前で読み込む）と、ちょっと前の記事で使った別のサンプルデータ（dat2という名前で読み込む）を使いましょう。

まず普通の線形回帰の場合。airqualityに対してlm(){base}関数を使ってみましょう。

> data(airquality)
> data(airquality) # データ読み込み
> airq<-airquality[,1:4] # 月・日付のデータを外す
> airq.lm<-lm(Ozone~.,airq) # Ozone = a * Solar.R + b * Wind + c * Temp + dのモデルを推定する
> summary(airq.lm) # 結果を表示する

Call:
lm(formula = Ozone ~ ., data = airq)

Residuals:
    Min      1Q  Median      3Q     Max 
-40.485 -14.219  -3.551  10.097  95.619 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -64.34208   23.05472  -2.791  0.00623 ** 
Solar.R       0.05982    0.02319   2.580  0.01124 *  
Wind         -3.33359    0.65441  -5.094 1.52e-06 ***
Temp          1.65209    0.25353   6.516 2.42e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 21.18 on 107 degrees of freedom
  (42 observations deleted due to missingness)
Multiple R-squared:  0.6059,	Adjusted R-squared:  0.5948 
F-statistic: 54.83 on 3 and 107 DF,  p-value: < 2.2e-16

これは見たまんまですね。airqualityは典型的な線形回帰モデルが適用できるありがちな連続値データなので、普通の結果になっています。では、dat1に当てはめてみるとどうなるでしょうか？

> dat1.lm<-lm(click~.,dat1)
> summary(dat1.lm)

Call:
lm(formula = click ~ ., data = dat1)

Residuals:
    Min      1Q  Median      3Q     Max 
-676.00 -293.64  -35.29  203.50  971.09 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  729.253    250.522   2.911  0.00513 ** 
d1           668.344    126.867   5.268 2.19e-06 ***
d2          -212.361    129.814  -1.636  0.10737    
d3           177.564    107.794   1.647  0.10501    
d4             2.207     99.730   0.022  0.98242    
d5          -100.639    115.238  -0.873  0.38616    
d6           277.451    129.814   2.137  0.03688 *  
d7            49.890    104.279   0.478  0.63418    
d8           -18.690    109.666  -0.170  0.86528    
d9           216.890    131.009   1.656  0.10331    
d10         -160.735    115.783  -1.388  0.17046    
d11          -67.000    118.341  -0.566  0.57351    
d12          -98.109    104.941  -0.935  0.35379    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 378.9 on 57 degrees of freedom
Multiple R-squared:  0.6168,	Adjusted R-squared:  0.5361 
F-statistic: 7.646 on 12 and 57 DF,  p-value: 3.431e-08

パッと見では、うまく説明できているように見えます。けれども、このdat1というデータは元の記事でも触れられている通り、クリック数という「カウントデータ＝離散値データ」です。なので、基本的には普通の線形回帰だと正しく説明できない可能性があります。

そこで、まずポアソン回帰モデルにしてみましょう。

> dat1.glm<-glm(click~.,dat1,family=poisson)
> summary(dat1.glm)

Call:
glm(formula = click ~ ., family = poisson, data = dat1)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-23.149   -9.143   -2.092    5.004   32.048  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  6.411354   0.021168 302.884  < 2e-16 ***
d1           0.736294   0.011278  65.285  < 2e-16 ***
d2          -0.120134   0.009868 -12.174  < 2e-16 ***
d3           0.199693   0.009609  20.783  < 2e-16 ***
d4           0.004421   0.007941   0.557    0.578    
d5          -0.060423   0.009406  -6.424 1.33e-10 ***
d6           0.223852   0.010059  22.253  < 2e-16 ***
d7           0.049067   0.008805   5.572 2.51e-08 ***
d8          -0.016761   0.008551  -1.960    0.050 *  
d9           0.247857   0.011252  22.027  < 2e-16 ***
d10         -0.188730   0.009733 -19.391  < 2e-16 ***
d11         -0.079104   0.009869  -8.016 1.10e-15 ***
d12         -0.087549   0.008347 -10.489  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 21668.7  on 69  degrees of freedom
Residual deviance:  8720.3  on 57  degrees of freedom
AIC: 9350.6

Number of Fisher Scoring iterations: 5

だいぶ結果の見かけが変わりましたね。有意なβ値を示す説明変数がかなり増えています。ただし元の記事でも指摘したように、このデータは標本平均と標本分散の比が1よりかなり大きいので、これではちょっとまずいです。負の二項分布に変えてみましょう。

> dat1.glmnb<-glm.nb(click~.,dat1)
> summary(dat1.glmnb)

Call:
glm.nb(formula = click ~ ., data = dat1, init.theta = 6.40747075, 
    link = log)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.4451  -0.7618  -0.1643   0.4616   2.5026  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept)  6.48958    0.26216  24.754  < 2e-16 ***
d1           0.71567    0.13279   5.389 7.07e-08 ***
d2          -0.15583    0.13574  -1.148   0.2510    
d3           0.23391    0.11288   2.072   0.0382 *  
d4          -0.05671    0.10433  -0.544   0.5868    
d5          -0.04942    0.12058  -0.410   0.6819    
d6           0.15940    0.13581   1.174   0.2405    
d7           0.11804    0.10915   1.081   0.2795    
d8          -0.08133    0.11470  -0.709   0.4783    
d9           0.19946    0.13711   1.455   0.1457    
d10         -0.18792    0.12117  -1.551   0.1209    
d11         -0.08563    0.12386  -0.691   0.4893    
d12         -0.07408    0.10978  -0.675   0.4998    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(6.4075) family taken to be 1)

    Null deviance: 158.029  on 69  degrees of freedom
Residual deviance:  71.921  on 57  degrees of freedom
AIC: 1052.4

Number of Fisher Scoring iterations: 1


              Theta:  6.41 
          Std. Err.:  1.07 

 2 x log-likelihood:  -1024.355

これまた結果が大きく変わりました。しかも、有意な説明変数の個数が減った一方でその組み合わせは普通の線形回帰の場合とは異なります。上の早見表に従えば、この負の二項分布回帰モデルの方が妥当なので、これが理論的には最適なモデルということになります。

では、dat2はどうでしょうか？　一旦cvの列をfactor型からnumeric型に直した上で、普通に線形回帰モデルを推定してみると。。。

> dat2.lm<-lm(cv~.,dat2)
> summary(dat2.lm)

Call:
lm(formula = cv ~ ., data = dat2)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.08346 -0.08003 -0.02191  0.08174  1.00823 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.187056   0.014339  82.783  < 2e-16 ***
a1           0.052277   0.008507   6.145 9.03e-10 ***
a2          -0.027412   0.008482  -3.232  0.00124 ** 
a3           0.001369   0.008384   0.163  0.87028    
a4          -0.167875   0.010358 -16.207  < 2e-16 ***
a5           0.086887   0.008907   9.755  < 2e-16 ***
a6           0.754514   0.010085  74.813  < 2e-16 ***
a7           0.002729   0.008379   0.326  0.74469    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2292 on 2992 degrees of freedom
Multiple R-squared:  0.7903,	Adjusted R-squared:  0.7899 
F-statistic:  1611 on 7 and 2992 DF,  p-value: < 2.2e-16

何となく良さそうに見えますね。ところが。。。このdat2の目的変数cvは元記事にもあるように本来はyes/noの二値で表されるカテゴリカルデータです。なので、その分布は正規分布に従うとは限りません。なのでこれに対して線形回帰モデルを適用するのは危ないです。

もちろん、元記事に従ってここでは二項分布を充てるロジスティック回帰モデルを用いるのが妥当でしょう。

> dat2.glm<-glm(cv~.,dat2,family=binomial)
> summary(dat2.glm)

Call:
glm(formula = cv ~ ., family = binomial, data = dat2)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.6404  -0.2242  -0.0358   0.2162   3.1418  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept) -1.37793    0.25979  -5.304 1.13e-07 ***
a1           1.05846    0.17344   6.103 1.04e-09 ***
a2          -0.54914    0.16752  -3.278  0.00105 ** 
a3           0.12035    0.16803   0.716  0.47386    
a4          -3.00110    0.21653 -13.860  < 2e-16 ***
a5           1.53098    0.17349   8.824  < 2e-16 ***
a6           5.33547    0.19191  27.802  < 2e-16 ***
a7           0.07811    0.16725   0.467  0.64048    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 4158.9  on 2999  degrees of freedom
Residual deviance: 1044.4  on 2992  degrees of freedom
AIC: 1060.4

Number of Fisher Scoring iterations: 7

線形回帰モデルの推定結果と比べてみると、有意な説明変数の組み合わせとβ値の符号は一致していますが、その比率や有意確率などは皆異なっています。これではまずいですね。

ということで、同じデータに対して本来使われるべき一般化線形モデルではなく普通の線形回帰モデルを適用すると、結果が変わってしまうということがこれで分かりました。

リンク関数が違うだけで、全ての一般化線形モデルはノーマルな線形回帰に帰着する

ひとまず3つの一般化線形モデルを見てきたわけですが、それぞれのモデルの数式を最後におさらいしてみましょう。

ポアソン回帰モデル

$log(\lambda_i) = \alpha + \beta_1 x_{1i} + \beta_2 x_{2i} + \beta_3 x_{3i} + \cdots + \beta_k x_{ki}$

負の二項分布回帰モデル

$log(\frac{\rho}{1+\rho}) = log(\frac{\mu}{\mu+k}) = \alpha + \beta_1 x_{1i} + \beta_2 x_{2i} + \beta_3 x_{3i} + \cdots + \beta_k x_{ki}$

ロジスティック回帰モデル

$log(\frac{p}{1-p}) = \alpha + \beta_1 x_{1i} + \beta_2 x_{2i} + \beta_3 x_{3i} + \cdots + \beta_k x_{ki}$

左辺はリンク関数で変換した結果なんですが、右辺は実は下記の線形回帰モデルと全く同じなんですね。

線形回帰モデル

$y_i = \alpha + \beta_1 x_{1i} + \beta_2 x_{2i} + \beta_3 x_{3i} + \cdots + \beta_k x_{ki}$

以前の記事でも書いたように、線形回帰モデル自体は最小二乗法(OLS)もしくは最尤法(MLE)に基づいて機械的に極めて容易に推定することができます（glm(){base}関数は最尤法に基づいて計算しています*3）。なので、いかなる分布を用いようとも計算プロセス的には~~基本的にはリンク関数による変換を経た後では~~同じプロセスをたどっていることになります。

にもかかわらず、その分布＆リンク関数のチョイスを間違えると結果も変わってしまうわけです。一口に「重回帰分析」と言っても、きちんとデータの性質に基づいて様々な手法を使い分けるべきだという所以です。なんてことを、adtech tokyoの会場で言ったらうざがられそうですが。。。（笑）

追記1

読み返していたら色々問題がある*4ことに気がついたので、そのうち補足の記事を一本立ててちゃんと書こうかと思います。気長にお待ちを。。。

ということでツッコミ大歓迎ですm(_ _)m

追記2

2件まとまったツッコミを頂きました。勉強中の身には大変助かります。。。

ということで、個人的にはこの先何を勉強するかの指針が立って至れり尽くせり。どうも有難うございました。なのですが、1番目のツッコミの結論を見てちょっと詰まってしまいました。自家製のサンプルデータが良くなかったんでしょうか。。。

それから、今回問題になったポイントを理解する上で久保先生の以下のページ2つが参考になると思われるので、併せて列挙しておきます。

この次にまとめ記事を作る際には、ここで指摘されているポイントも織り込む予定です。

追記3

ピンポン状態になってますが、uncorrelatedさんの1番目のブログ記事に加筆がされていました。

追記(2013/09/20 09:38):問題に気付かれたようでブログ主が「自家製のサンプルデータが良くなかったんでしょうか。。。」と追記でつぶやいていた。dat1$clickの値が大きすぎるので、例えば以下のような値にすればポアソン回帰が良く当てはまる。
dat1$click <- c(4,3,9,3,9,5,6,5,7,17,6,7,6,3,9,4,16,1,11,5,4,1,8,5,2,6,18,10,3,6,9,3,4,2,5,3,8,23,10,7,8,18,13,5,16,6,3,15,17,20,25,11,9,6,6,17,5,21,24,6,23,20,8,15,14,9,4,11,13,9)
当てはまりすぎて、glm.nb関数の方に警告が出るようになるが。

平均値が9ぐらいと小さいカウントデータであれば、ポアソン回帰でも負の二項分布回帰でも良く当てはまる、ということですね。

個人的にはこれは結構重要なポイントだと思っていて、click総数やPV数やimpression数といった平均値のでかいデータなら正規分布線形回帰で良いわけですが、CV実測数のようなとんでもなく少なくなることもある*5データの場合はやはりポアソン回帰や負の二項分布回帰が有効になるということかと。

だいぶ考えがすっきりしてきたので、そのうちまとめます。

*1:キャンペーンガールのお姉さんたちを眺めに行ってるわけではありません。。。多分（笑）

*2:そもそもRでの実践例とサンプルコードまである

*3:参考資料： http://www.geocities.jp/ancientfishtree/R_nra.html および http://mjin.doshisha.ac.jp/R/16.html および http://www.agr.nagoya-u.ac.jp/~seitai/document/R2009/t090614glmIntro.pdf

*4:もっと言ってしまうとサンプルデータの説明変数に多重共線性がある

*5:ビジネス区分によっては1日に1件しかCVがないなんてことも