これはディスプレイ広告に限った話ではないと思うんですが、あるPC / スマホ上の何かしらのクリエイティブに対するクリック数がそのデザインの良し悪しによって左右されるということは、web業界ではよく知られているかと思います。
そういう場合「どんなデザインがクリック数を増やすのに有効か？」というのは、厳密にはきちんと条件統制をかけて実験計画法に基づいてデザインしたA/Bテストなどで調べるべきなんでしょうが、そこまで綿密にやっている余裕のない現場も結構多いはず。
そこで、今回は既に計測済みの各広告のクリック数（CTRでしか得られていないようであれば実クリック数に直すものと想定する）データが得られているものと仮定して、それを各広告のデザイン要素を表すインデックス（二値orカテゴリカルデータ）のデータと組み合わせて、「どんなデザインをすればクリック数が増えるか？」を推定するというケースを想定してRでどうすれば良いかを考えてみようと思います。
ちなみに、今回用いるポアソン回帰モデルを初めとするGLM周りの基礎知識についてはこちらの本が最も分かりやすいと思います。書評を途中まで書いて、恐れ多くて結局書くのをやめてしまったんですが、分かりやすさでは他の追随を許さない名著です。
データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)
- 作者: 久保拓弥
- 出版社/メーカー: 岩波書店
- 発売日: 2012/05/19
- メディア: 単行本
- 購入: 16人 クリック: 163回
- この商品を含むブログ (17件) を見る
また、同じ久保先生の手による解説ページ（生態学データ解析 - 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がマイナス（＝クリック数減につながる）ということになるので、結構嫌なデータだと思います。
参考資料など
今回のテーマで記事を書くに当たって、以下のページを参考にしました。
本当は多項ロジットモデルとかまで行きたかったんですが、今回はこの辺で。。。そもそも今回作ったサンプルデータはマルチコとか何も一切考慮してなかったりして厳密なことを言おうとすると面倒なので（笑）。
そうそう、今回も僕にとっては勉強中のテーマです。間違っている点・あやふやな点などあれば、ぜひぜひご指摘下さい。
追記（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%ぐらいしか説明できてない）からなんだと思いますが。。。