ディスプレイ広告のクリック数改善のためのデザイン最適化に、ポアソン分布 / è² ã®äºŒé …åˆ†å¸ƒã«ã‚ˆã‚‹GLMを使ってみる

これはディスプレイ広告に限った話ではないと思うんですが、あるPC / スマホ上の何かしらのクリエイティブに対するクリック数がそのデザインの良し悪しによって左右されるということは、web業界ではよく知られているかと思います。

そういう場合「どんなデザインがクリック数を増やすのに有効か？」というのは、厳密にはきちんと条件統制をかけて実験計画法に基づいてデザインしたA/Bテストなどで調べるべきなんでしょうが、そこまで綿密にやっている余裕のない現場も結構多いはず。

そこで、今回は既に計測済みの各広告のクリック数（CTRでしか得られていないようであれば実クリック数に直すものと想定する）データが得られているものと仮定して、それを各広告のデザイン要素を表すインデックス（二値orカテゴリカルデータ）のデータと組み合わせて、「どんなデザインをすればクリック数が増えるか？」を推定するというケースを想定してRでどうすれば良いかを考えてみようと思います。

ちなみに、今回用いるポアソン回帰モデルを初めとするGLM周りの基礎知識についてはこちらの本が最も分かりやすいと思います。書評を途中まで書いて、恐れ多くて結局書くのをやめてしまったんですが、分かりやすさでは他の追随を許さない名著です。

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

また、同じ久保先生の手による解説ページ（生態学データ解析 - FAQ 一般化線形モデル）も大いに参考になるのではないかと。

必要なRパッケージ

{MASS}と{randomForest}をインストールして展開して下さい。あ、タイトルにわざと書かなかったのがこれでバレちゃいましたね（笑）。

ポアソン分布によるGLM（ポアソン回帰モデル）でやってみる

いつも通りサンプルをGitHubに上げておきましたので落としてきてください*1。名前は何でもいいんですが、今回は適当にdat1とかしておきます。

データのイメージとしては、ディスプレイ広告のデザイン要素を12種類に分けて、それぞれの「ありorなし」の組み合わせで個々の広告を表現しているものとします（d1-d12の各カラムで「あり」なら1＆「なし」なら0）。そして、それぞれの広告のクリック数の実数が最後に記録されているものとします（clickカラム）。いかにもありがちなデータですね。

d1	d2	d3	d4	d5	d6	d7	d8	d9	d10	d11	d12	click
0	1	1	1	0	1	0	1	0	1	0	1	175
0	1	1	0	0	0	0	1	0	0	0	0	350
1	1	1	0	1	0	1	1	0	1	1	0	675
0	1	0	1	1	0	1	1	0	0	1	0	250
0	1	0	1	1	1	1	1	0	1	1	1	275
...	...	...	...	...	...	...	...	...	...	...	...	...

まずは、「目的変数：カウントデータ / 説明変数：二値データ」ということでセオリー通りポアソン回帰モデルでやってみようと思います。その回帰式は、説明変数を、目的変数をとすると以下の通り。

要は、ポアソン分布のパラメータλを線形モデルで説明可能なものと仮定して、その線形モデルを推定する目的変数を対数変換したものに対して、通常の線形（重）回帰モデルを推定する*2ものです。この辺の理屈は例えばこちらのページとか、もしくはid:uncorrelatedさんのブログ記事を参照のこと（ご指摘有難うございました）。ちなみに線形回帰モデルの基礎については以前の記事でちろっと触れているので、そちらもお読みあれ。

Rでの計算自体は、glm(){stats}でfamily=poissonとすれば簡単に実行できます。

> 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

一応、step()で変数選択を行ってみます。久保先生のコメントに従い、一旦ここでは交互作用のことは忘れたことにします。

> dat1.glm2<-step(dat1.glm)
Start:  AIC=9350.63
click ~ d1 + d2 + d3 + d4 + d5 + d6 + d7 + d8 + d9 + d10 + d11 + 
    d12

       Df Deviance     AIC
- d4    1   8720.6  9348.9
<none>      8720.3  9350.6
- d8    1   8724.1  9352.5
- d7    1   8751.6  9379.9
- d5    1   8761.3  9389.6
- d11   1   8784.1  9412.5
- d12   1   8831.1  9459.5
- d2    1   8868.7  9497.0
- d10   1   9093.1  9721.5
- d3    1   9162.3  9790.6
- d6    1   9206.5  9834.8
- d9    1   9206.6  9835.0
- d1    1  13026.4 13654.7

Step:  AIC=9348.94
click ~ d1 + d2 + d3 + d5 + d6 + d7 + d8 + d9 + d10 + d11 + d12

       Df Deviance     AIC
<none>      8720.6  9348.9
- d8    1   8724.5  9350.8
- d7    1   8752.5  9378.9
- d5    1   8762.0  9388.3
- d11   1   8784.2  9410.5
- d12   1   8831.8  9458.2
- d2    1   8869.2  9495.5
- d10   1   9101.1  9727.5
- d3    1   9172.8  9799.2
- d9    1   9206.9  9833.2
- d6    1   9209.6  9835.9
- d1    1  13027.5 13653.9
> summary(dat1.glm2)

Call:
glm(formula = click ~ d1 + d2 + d3 + d5 + d6 + d7 + d8 + d9 + 
    d10 + d11 + d12, family = poisson, data = dat1)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-23.118   -9.118   -2.065    4.919   32.058  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)  6.413185   0.020914 306.643  < 2e-16 ***
d1           0.736299   0.011281  65.271  < 2e-16 ***
d2          -0.120217   0.009867 -12.183  < 2e-16 ***
d3           0.200392   0.009528  21.032  < 2e-16 ***
d5          -0.059299   0.009187  -6.455 1.08e-10 ***
d6           0.224234   0.010039  22.336  < 2e-16 ***
d7           0.049432   0.008782   5.629 1.81e-08 ***
d8          -0.016869   0.008549  -1.973   0.0485 *  
d9           0.247561   0.011240  22.026  < 2e-16 ***
d10         -0.189420   0.009656 -19.616  < 2e-16 ***
d11         -0.078933   0.009865  -8.002 1.23e-15 ***
d12         -0.087681   0.008344 -10.509  < 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.6  on 58  degrees of freedom
AIC: 9348.9

Number of Fisher Scoring iterations: 5

基本的には偏回帰係数であるEstimateがプラスの値であればクリック数増に貢献し、逆にマイナスの値であればクリック数減につながってしまいます。またEstimateの値の大きさ（絶対値）がクリック数への影響の大小を表します。

ということを踏まえて読み解いてみると、どの偏回帰係数も統計的には有意という判定になっている点が気になるものの、とりあえずd1のデザイン要素がクリック数増に大きく貢献していることが分かります。他にもクリック増に貢献しそうなのはd9, d6, d3、逆にクリック減につながってしまいそうなのはd10,d2ということも見て取れます。

負の二項分布によるGLMでやってみる

ところで、そもそもこのデータのclickカラムについて以下のように計算してみると。。。

> mean(dat1[,13])
[1] 1067.5
> var(dat1[,13])
[1] 309571.6

となるのを見れば分かる通り、平均値より分散の方がずっと大きく、ポアソン分布をただ当てはめるのは本来ならばあまり適切ではないようです。またポアソン回帰モデルによる結果を見ると"Residual deviance" / "degree of freedom"の値が非常に大きく*3、モデルの当てはまりとしては良くないということが分かります。これでは、変数重要度は分かりそうな気がしてもモデルとしてはちょっとよろしくない感じです。

そこで、久保先生のアドバイスに従って負の二項分布によるGLMを導入してみましょう。負の二項分布については群馬大の青木先生による解説や、Matlabヘルプサイトの解説が分かりやすいかと思います。

僕も勉強中なので理解不足なのですが、このモデルではポアソン分布に比べて平均値と分散との乖離が大きいようなカウントデータまで含めて広汎に表現できる負の二項分布を使い、さらに正確にモデリングすることを目指している、とのことです。その回帰式は以下の通り。

これをRのglm(){stats}でやろうとしてもオプションがないので、新しくglm.nb(){MASS}でやることになります。

> 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 

これで"Residual deviance" / "degree of freedom"の値が1.26ぐらいと妥当な数字に収まり、モデルとしての当てはまりも悪くないという結果になりました。

偏回帰係数で統計的に有意なのはd1とd3のみで、（ポアソン回帰モデルの拡張だから当たり前なんですが）だいたいポアソン回帰モデルの時と近い値になっています。

ランダムフォレスト回帰でやってみる

実は、この手の回帰モデルもランダムフォレストである程度計算することができます。ただしランダムフォレストでは偏回帰係数の「符号」までは出せませんので、基本的には変数重要度を見るだけということになります。

> dat1.rf<-randomForest(click~.,dat1,ntree=2000)
> print(dat1.rf)

Call:
 randomForest(formula = click ~ ., data = dat1, ntree = 2000) 
               Type of random forest: regression # 分類問題ではなく回帰問題
                     Number of trees: 2000
No. of variables tried at each split: 4

          Mean of squared residuals: 181438.4
                    % Var explained: 40.54
> importance(dat1.rf)
    IncNodePurity
d1      5556729.2
d2      1951166.3
d3       979182.8
d4       679803.4
d5       591530.5
d6       891261.9
d7       678428.5
d8      1087023.2
d9      3094248.6
d10      661204.0
d11      557049.7
d12      486588.2

ということで順位付けをしていくと、d1 -> d9 -> d2 -> d8 -> d3 -> d6 -> …の順に影響が大きいということが分かりました。これ、上のポアソン分布 / 負の二項分布によるGLMで算出した偏回帰係数の「符号（プラスorマイナス）」と見比べてみると、影響の大きさで上位に来る6変数のうち、d2とd8がマイナス（＝クリック数減につながる）ということになるので、結構嫌なデータだと思います。

参考資料など

今回のテーマで記事を書くに当たって、以下のページを参考にしました。

• Rとカテゴリカルデータのモデリング(2)
• GLM入門編1 - Murakami's Memorandum

本当は多項ロジットモデルとかまで行きたかったんですが、今回はこの辺で。。。そもそも今回作ったサンプルデータはマルチコとか何も一切考慮してなかったりして厳密なことを言おうとすると面倒なので（笑）。

そうそう、今回も僕にとっては勉強中のテーマです。間違っている点・あやふやな点などあれば、ぜひぜひご指摘下さい。

追記（SVRでやってみた＋予測してみた）

ここまでやっておいてSVR（サポートベクター回帰）しないのは何だかなぁと思ったので、{e1071}パッケージのSVM*4でやってみました。

> library("e1071", lib.loc="C:/Program Files/R/R-3.0.1/library")
Loading required package: class
> dat1.svm<-svm(click~.,dat1)
> print(dat1.svm)

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


Parameters:
   SVM-Type:  eps-regression 
 SVM-Kernel:  radial 
       cost:  1 
      gamma:  0.08333333 
    epsilon:  0.1 


Number of Support Vectors:  66

これだけだとどうモデルが違うのかよく分からないと言われそうなのでpredict()で予測しようと思ったんですが、折角なので他の回帰モデルも併せて予測してみることにしました。ちなみに予測のために入力値として選んだのは、d1-d12まで全て1で表現される「全部盛り」ベクトルxinです。

> predict(dat1.glm,newdata=xin,type="response")
       1 
1509.956 
> predict(dat1.glmnb,newdata=xin,type="response")
       1 
1373.509 
> predict(dat1.rf,xin)
       1 
1424.997 
> predict(dat1.svm,xin)
       1 
1287.155 

おや、結構モデルごとに予測結果が違いますね。ちなみにd1-d12まで全て0で表現される「全部なし」ベクトルxin0を入力して予測すると、

> predict(dat1.glm,newdata=xin0,type="response")
       1 
608.7171 
> predict(dat1.glmnb,newdata=xin0,type="response")
      1 
658.246 
> predict(dat1.rf,xin0)
       1 
865.2397 
> predict(dat1.svm,xin0)
       1 
805.4882 

これまた結構違いますね。最後に、ポアソン回帰モデルで得られた偏回帰係数の「符号」に合わせて1 / 0を振り分けた「多分クリック数が最大値になるはず」ベクトルxinfを入力して予測すると、

> predict(dat1.glm,newdata=xinf,type="response")
       1 
2624.216 
> predict(dat1.glmnb,newdata=xinf,type="response")
       1 
2589.834 
> predict(dat1.rf,xinf)
       1 
1679.118 
> predict(dat1.svm,xinf)
       1 
1692.253 

生データの最大値が2325なんですが、ポアソン回帰モデルと負の二項分布モデルではこれを上回る予測値になっています。一方で、ランダムフォレストとSVRによる予測値はこれを下回ってます。えーと、これ何でだったっけ。。。勉強しときます（汗）。多分、今回のサンプルデータのうち予測に効く変数が少なくて、ランダムフォレストでもSVRでも全体の分散に対する説明力が弱い（実際ランダムフォレストでも40%ぐらいしか説明できてない）からなんだと思いますが。。。

*1:今回も適当に僕がrnorm()を組み合わせて作ったシミュレーションデータです

*2:ご指摘を頂いたので訂正しました。 http://b.hatena.ne.jp/uncorrelated/20130919#bookmark-161101866 および http://d.hatena.ne.jp/uncorrelated/20130919/1379599682

*3:この比率が1前後でないとモデルとしては妥当ではないらしい

*4:つまりあの有名なマルチリンガルパッケージLIBSVMのR版