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

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

同じデータセットに対するアプローチの違いから見る「データ分析のステージ」

追記 (2015/02/21)


いくつか抜けてるところがあったなぁと思ったので、後から追記や加筆修正してみました。最初のオリジナル版から少し内容が変わっているところがありますがご了承ください。


ちょっと前の記事でこんなネタをやってみたわけですが。



世の中には色々な「データ分析」のやり方があるなぁと思った時に、この同じ2013年のテニス四大大会のデータからそれぞれのやり方をしている人たちがどんな異なるアプローチを取るのかなぁとふと想像したもので、半分ネタ的に書いてみました。便宜的に以下のようにステージを分けてあります。

  1. 集計ステージ
  2. 検定ステージ
  3. 相関ステージ
  4. 重回帰分析ステージ
  5. 機械学習を含めたモデリングステージ
  6. 厳密性に拘るステージ


なお、データは以前の記事と同じこちらのものをお使い下さい。

その上で、Rで分析する際は以下のように前処理しておきます。単にプレイヤー名・獲得ゲーム数・総獲得ポイント数など、以前の記事で不要と判断されたものを削除するだけですが。

> dm<-read.delim("men.txt")
> dw<-read.delim("women.txt")
> dm<-dm[,-c(1,2,16,17,18,19,20,21,34,35,36,37,38,39)]
> dw<-dw[,-c(1,2,16,17,18,19,20,21,34,35,36,37,38,39)]

ここまで終わったら、適当にやってみましょう。ちなみに「集計ステージ」は多分Excelでもできると思うので、Excelでやってみようかと思います。


そうそう、各変数の説明をもう一度下に貼っておきます。

Result Result of the match (0/1) - Referenced on Player 1 is Result = 1 if Player 1 wins (FNL.1>FNL.2)
FSP.1 First Serve Percentage for player 1 (Real Number)
FSW.1 First Serve Won by player 1 (Real Number)
SSP.1 Second Serve Percentage for player 1 (Real Number)
SSW.1 Second Serve Won by player 1 (Real Number)
ACE.1 Aces won by player 1 (Numeric-Integer)
DBF.1 Double Faults committed by player 1 (Numeric-Integer)
WNR.1 Winners earned by player 1 (Numeric)
UFE.1 Unforced Errors committed by player 1 (Numeric)
BPC.1 Break Points Created by player 1 (Numeric)
BPW.1 Break Points Won by player 1 (Numeric)
NPA.1 Net Points Attempted by player 1 (Numeric)
NPW.1 Net Points Won by player 1 (Numeric)
FSP.2 First Serve Percentage for player 2 (Real Number)
FSW.2 First Serve Won by player 2 (Real Number)
SSP.2 Second Serve Percentage for player 2 (Real Number)
SSW.2 Second Serve Won by player 2 (Real Number)
ACE.2 Aces won by player 2 (Numeric-Integer)
DBF.2 Double Faults committed by player 2 (Numeric-Integer)
WNR.2 Winners earned by player 2 (Numeric)
UFE.2 Unforced Errors committed by player 2 (Numeric)
BPC.2 Break Points Created by player 2 (Numeric)
BPW.2 Break Points Won by player 2 (Numeric)
NPA.2 Net Points Attempted by player 2 (Numeric)
NPW.2 Net Points Won by player 2 (Numeric)


これは前の記事でも書いたように、テニス経験者ならすぐそれぞれがどういう意味を持つ量かが分かるかと思います。そういう意味ではドメイン知識の有無*1も問われる妙に実践的なデータセットとも言えますねw


ということで、ここでは「プレイヤーの勝敗に関連する要因は何か」を追究するというのが目的であるとして、ステージ別に色々やってみることにします。


1. 集計ステージ


例えばこのデータの分析をExcelでやれと言われたら、僕が思い付くのは「Resultが0 / 1それぞれの場合で23個ある変数の平均値を取ってみる」というものでしょうか。これは以前一緒に働いたことのあるマーケッターの人たちもやっていたやり方です。


で、男子のデータだけそれでやってみました。単に"men.txt"をコピペしてResultの結果でソートし、36列それぞれの平均をAVERAGE関数で計算するだけです。ついでに棒グラフも描いてみることに。

f:id:TJO:20150220114629p:plain

うん、何か分かる気がする。でも何も分からない気もする。。。ちなみにRの方が多分手間がかからないです。

> colMeans(dm[dm$Result==0,-1])
    FSP.1     FSW.1     SSP.1     SSW.1     ACE.1     DBF.1     WNR.1     UFE.1     BPC.1     BPW.1     NPA.1     NPW.1     FSP.2 
60.860558 47.860558 39.139442 21.015936  7.701195  4.374502 24.366534 24.924303  3.278884  5.840637 17.043825 20.625498 62.147410 
    FSW.2     SSP.2     SSW.2     ACE.2     DBF.2     WNR.2     UFE.2     BPC.2     BPW.2     NPA.2     NPW.2 
51.956175 37.852590 23.043825 10.306773  3.557769 30.227092 21.414343  6.621514 10.243028 18.310757 20.685259 
> colMeans(dm[dm$Result==1,-1])
    FSP.1     FSW.1     SSP.1     SSW.1     ACE.1     DBF.1     WNR.1     UFE.1     BPC.1     BPW.1     NPA.1     NPW.1     FSP.2 
62.979167 50.687500 37.020833 21.837500 10.429167  3.458333 30.379167 20.945833  6.908333 10.454167 18.183333 20.516667 60.462500 
    FSW.2     SSP.2     SSW.2     ACE.2     DBF.2     WNR.2     UFE.2     BPC.2     BPW.2     NPA.2     NPW.2 
45.554167 39.537500 20.545833  7.575000  4.554167 23.858333 25.958333  2.804167  5.575000 18.066667 22.037500 


ところどころ「勝者の方が大きい値」とか「敗者の方が大きい値」とかが見えますが、これだけではちょっと分からないことばかりだなぁという印象がありますね。ただし大まかな傾向はこれで分かります。


一般に、集計のみである程度分析していこうと思ったらかなり緻密に絞り込んだ「仮説」が必要になります。このデータで言えばやはりテニスそのものについてのドメイン知識に基づき、例えば「UFE.1/2(アンフォーストエラー=自発的なミスの数)が少ない方が勝ちやすい」とか「WNR.1/2(ウィナー=ストロークエースの数)が多い方が勝ちやすい」とか、「ACE.1/2(サービスエースの数)が多い方が勝ちやすい」とか、そういうところをチェックするわけです。


ちなみに今回のデータを見ると、

Player1の勝敗 UFE.1の平均 UFE.2の平均 WNR.1の平均 WNR.2の平均 ACE.1の平均 ACE.2の平均
24.9 21.4 24.4 30.2 7.7 10.3
20.9 25.9 30.4 23.9 10.4 7.6


ということで、確かに上記の仮説通りになっています。なので、これらのデータが勝敗の予測に重要そうだということは集計するだけでも言えそうです。とは言え、ここで結論付けられることは限定的かなと。何故かというと、その分散を計算すると

> apply(dm[dm$Result==0,],2,var)
    Result      FSP.1      FSW.1      SSP.1      SSW.1      ACE.1      DBF.1      WNR.1 
  0.000000  54.968478 287.112478  54.968478  78.863745  39.554359   9.459187 436.161116 
     UFE.1      BPC.1      BPW.1      NPA.1      NPW.1      FSP.2      FSW.2      SSP.2 
429.102247  10.481912  26.854502 177.114072 191.587187  49.830183 247.810072  49.830183 
     SSW.2      ACE.2      DBF.2      WNR.2      UFE.2      BPC.2      BPW.2      NPA.2 
 63.714072  42.717514   7.695649 489.808223 360.187633  12.868175  29.320701 191.791044 
     NPW.2 
208.808542 
> apply(dm[dm$Result==1,],2,var)
    Result      FSP.1      FSW.1      SSP.1      SSW.1      ACE.1      DBF.1      WNR.1 
  0.000000  56.514209 232.901935  56.514209  63.986036  41.467765   6.341353 525.868183 
     UFE.1      BPC.1      BPW.1      NPA.1      NPW.1      FSP.2      FSW.2      SSP.2 
327.800401  12.744700  26.483246 216.510181 169.623152  54.082270 287.084920  54.082270 
     SSW.2      ACE.2      DBF.2      WNR.2      UFE.2      BPC.2      BPW.2      NPA.2 
 77.487430  38.705649   8.432200 389.963110 450.776499   7.948937  21.726569 210.221478 
     NPW.2 
234.973483 


となっており、先ほどの比較表と同じように並べてみると

Player1の勝敗 UFE.1の分散 UFE.2の分散 WNR.1の分散 WNR.2の分散 ACE.1の分散 ACE.2の分散
429.1 360.2 436.2 489.8 39.6 42.7
327.8 450.8 525.9 390.0 41.5 38.7


というように意外と分散が大きい印象があります。そこでExcelで棒グラフを描く際によく用いられる「平均±標準偏差」というやり方でプロットしてみると、以下のような感じになります。


f:id:TJO:20150221162616p:plain


ものすごーーーく標準偏差が大きい気がします(汗)。ということで、単にこのプロットを見ているだけではちょっと結論を出せる気がしないですね。そこで次のステージに進む必要がありそうだと思われます。


2. 検定ステージ


データ集計が当たり前になった現場では一般に「検定」が好まれる傾向があると個人的には感じています。特にt検定とかカイ二乗検定とかが好まれるようです。このデータセットに対しても同じように(例えば)t検定を使ってみる、ということを考える人は結構な数いると思います。


実際、上記の集計結果のようにエラーバー(標準偏差)が思ったよりも大きく、棒グラフにしてもエラーバー同士が重なるようなケースでは、なかなか差があるとは確信が持てないこともあるかと。


例えば、男子のACE.1のデータなんかは明らかに差がありそうなので、検定で大小を確かめてみる価値はあるんじゃないかと考えられます。ではこれをやってみようと思うわけですが。。。

f:id:TJO:20150220115533p:plain

うあー、Excelの「分析ツール」だとそもそもどれを使ったらいいんだか分かりにくいですね。でも結果としては「勝者と敗者とでACE.1(サービスエースの本数)には有意差がある」ということのようです。ちなみにRだともっとサクッと簡単にできます。

> t.test(dm$ACE.1[dm$Result==0],dm$ACE.1[dm$Result==1])

	Welch Two Sample t-test

data:  dm$ACE.1[dm$Result == 0] and dm$ACE.1[dm$Result == 1]
t = -4.7461, df = 486.716, p-value = 2.73e-06
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -3.857323 -1.598619
sample estimates:
mean of x mean of y 
 7.701195 10.429167


なお、このケースではそれぞれの分散は

> var(dm$ACE.1[dm$Result==0])
[1] 39.55436
> var(dm$ACE.1[dm$Result==1])
[1] 41.46776


ということで、等しい可能性もありますが厳密には異なるようです。このように迷う状況であれば、Welchの検定を用いるのが正しいとみるべきでしょう*2


ただ、ここで個別にACE.1がどうのとか、SSW.1がどうのとか検定しながらあれこれ言っても、どの値が勝敗に関連するかは見えてこないように思われます。というより、これ以外の山のような変数で全て対比較するように検定するのはあまりにも非効率なやり方だと考えられますね*3


3. 相関ステージ


一般に全て対比較するのは恐ろしくつらいので、何かもっと一括して全体の関係性を見られる方法があった方が良いと思う人は多いです。で、結構よく使われるのが相関係数。このケースだとResultフラグとの相関を計算することである程度自動的に「勝敗と最も強く関連しているのはどのデータか」を占うことができると期待されます。


Excelでやっても良いのですが(CORREL関数でベタベタやれば良い)、何となくRで同じことをやってみます。

> cor(dm)[1,-1]
       FSP.1        FSW.1        SSP.1        SSW.1        ACE.1        DBF.1        WNR.1 
 0.140744830  0.087377313 -0.140744830  0.048578332  0.209955001 -0.160787529  0.136183284 
       UFE.1        BPC.1        BPW.1        NPA.1        NPW.1        FSP.2        FSW.2 
-0.101752613  0.471084156  0.408427895  0.040697570 -0.004053553 -0.116344509 -0.192573843 
       SSP.2        SSW.2        ACE.2        DBF.2        WNR.2        UFE.2        BPC.2 
 0.116344509 -0.147447547 -0.209574829  0.173187767 -0.150184323  0.112455740 -0.508852459 
       BPW.2        NPA.2        NPW.2 
-0.419437013 -0.008627817  0.045453847 


cor関数で計算した1行目だけを見れば、Resultフラグとの相関を計算した結果が得られます。ついでなので棒グラフにしてみましょう。Excelで同じことをやったものも下に並べておきます。

f:id:TJO:20150221163659p:plain

f:id:TJO:20150221164033p:plain

これでBPC.1/2(ブレーク機会数), BPW.1/2(ブレークポイント獲得数)が強く寄与するということはよく分かりました。けれども、ここで計算した単相関は実は変数間の関係性によっては誤った結果になっていることもある(なぜ項目ごとに単純な集計をするより、多変量解析(重回帰分析)をした方が正確な結果を返すのか - 銀座で働くデータサイエンティストのブログ)ので、本当は偏相関を計算しなければいけない場面でして。。。となるとExcelでは難しいかも。


その辺の問題は重回帰分析を含む重回帰モデルである程度解決可能であり、ゆえに重回帰分析がこの手の多変量解析で多用される所以でもあります。


4. 重回帰分析ステージ


そんなわけで検定や相関では飽き足らなくなってきた人が興味を持ってちらほらやり始めるのが「重回帰分析」。Excelでもやりようによっては出来るし、Rならお手軽に出来る。というわけでさすがにここからはRでやってみましょう。もちろん皆さんご存知lm関数で簡単に実行できます。

> dm.lm<-lm(Result~.,dm)
> summary(dm.lm)

Call:
lm(formula = Result ~ ., data = dm)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.83231 -0.16572 -0.00958  0.15599  0.91011 

Coefficients: (2 not defined because of singularities)
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.5622390  0.2790697   2.015  0.04451 *  
FSP.1        0.0005964  0.0031624   0.189  0.85050    
FSW.1        0.0148113  0.0021471   6.898 1.71e-11 ***
SSP.1               NA         NA      NA       NA    
SSW.1        0.0146962  0.0033140   4.435 1.15e-05 ***
ACE.1        0.0008919  0.0028465   0.313  0.75417    
DBF.1       -0.0129911  0.0053800  -2.415  0.01613 *  
WNR.1        0.0045199  0.0015284   2.957  0.00326 ** 
UFE.1       -0.0015501  0.0015712  -0.987  0.32437    
BPC.1        0.0330029  0.0049656   6.646 8.38e-11 ***
BPW.1        0.0227510  0.0034899   6.519 1.83e-10 ***
NPA.1        0.0036844  0.0025294   1.457  0.14589    
NPW.1       -0.0037674  0.0023192  -1.624  0.10495    
FSP.2       -0.0007562  0.0032574  -0.232  0.81653    
FSW.2       -0.0191987  0.0022646  -8.478 3.02e-16 ***
SSP.2               NA         NA      NA       NA    
SSW.2       -0.0126902  0.0031146  -4.074 5.42e-05 ***
ACE.2        0.0017076  0.0027408   0.623  0.53357    
DBF.2        0.0017989  0.0055318   0.325  0.74517    
WNR.2       -0.0031550  0.0016141  -1.955  0.05121 .  
UFE.2        0.0007894  0.0015282   0.517  0.60573    
BPC.2       -0.0324298  0.0053581  -6.052 2.93e-09 ***
BPW.2       -0.0174061  0.0035668  -4.880 1.46e-06 ***
NPA.2       -0.0017918  0.0025950  -0.690  0.49023    
NPW.2        0.0043482  0.0021816   1.993  0.04683 *  
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2708 on 468 degrees of freedom
Multiple R-squared:  0.7203,	Adjusted R-squared:  0.7071 
F-statistic: 54.78 on 22 and 468 DF,  p-value: < 2.2e-16


見た感じちらほら有意なパラメータもあるみたいだし、R^2 = 0.7071とそこそこ高いし、なかなか悪くないんじゃないでしょうか。


ちなみに単相関を計算した時にはFSP.1(1stサーブパーセンテージ) > FSW.1(1stサーブ時獲得ポイント数)だったのですが、重回帰分析で得られた偏回帰係数(coefficients)を見るとFSP.1 < FSW.1となっています。これ、よくよく考えたら当たり前のことで、FSP.1が高くなるほどFSW.1の母集団サンプルが増えるんですね。意外と見落としがちな変数間相関がこんなところにあるわけで、その辺も重回帰であれば解決できるというわけです。


ということで、回帰モデルのご利益のもう一つが「予測(or当てはめ)」。調子に乗って男子の結果・女子の結果それぞれに当てはめてみると。。。

> cor(dm$Result,predict(dm.lm,newdata=dm[,-1]))
[1] 0.8487034
Warning message:
In predict.lm(dm.lm, newdata = dm[, -1]) :
   ランク落ちした状態で当てはめを行っているため、誤った予測をしている可能性があります 
> cor(dw$Result,predict(dm.lm,newdata=dw[,-1]))
[1] 0.8540651
Warning message:
In predict.lm(dm.lm, newdata = dw[, -1]) :
   ランク落ちした状態で当てはめを行っているため、誤った予測をしている可能性があります 


学習データである男子のデータセットへの自己当てはめは良いのですが、女子のデータセットへの当てはまりが良くないのかも? ということで女子だけ別にやってみると、

> dw.lm<-lm(Result~.,dw)
> cor(dw$Result,predict(dw.lm,newdata=dw[,-1]))
[1] 0.8807626


自己当てはめに改めてみたら、結果が良くなりました。ということで個々の勝敗の結果を男子のデータセットからランダムに6サンプル抽出してきて、男子の重回帰分析結果を使って予測してみようとすると、

> idx<-sample(491,6,replace=F)
> predict(dm.lm,newdata=dm[idx,-1])
        29        471        379        371        193        416 
0.23144971 0.50944218 0.54150387 0.14365001 1.31922667 0.03788432 
Warning message:
In predict.lm(dm.lm, newdata = dm[idx, -1]) :
   ランク落ちした状態で当てはめを行っているため、誤った予測をしている可能性があります 


うわ、そもそも0/1じゃない結果だった。。。これではちょっと使えないですね。二値の場合に適した手法を選択しなければならないようです。


5. 機械学習を含めたモデリングステージ


基本的に「テニスの勝敗」というのは二値データです。ということは、機械学習における分類器を用いるのがおそらく相性が良くて、その上で例えばどの説明変数が分類のために重要だったかを調べる、というのが妥当だろうと考えられますね。


ということで、ここまで来てようやく以前の記事同様の流れにたどり着きます。おさらいの意も込めてもう一度やってみましょう。まず、二値カテゴリ型であることを明示するために目的変数をfactor型に直します。

> dm$Result<-as.factor(dm$Result)
> dw$Result<-as.factor(dw$Result)


その上で、基本形としてロジスティック回帰SVMランダムフォレストをやってみます。以前の記事同様、男子のデータセットで学習させ、女子のデータセットでテストします。まずロジスティック回帰から。

> dm.glm<-glm(Result~.,dm,family=binomial)
> summary(dm.glm)

Call:
glm(formula = Result ~ ., family = binomial, data = dm)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-2.91243  -0.14289  -0.00641   0.08835   3.04139  

Coefficients: (2 not defined because of singularities)
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -1.378486   5.064233  -0.272 0.785468    
FSP.1        0.116125   0.065947   1.761 0.078258 .  
FSW.1        0.198896   0.043094   4.615 3.92e-06 ***
SSP.1              NA         NA      NA       NA    
SSW.1        0.290083   0.069708   4.161 3.16e-05 ***
ACE.1       -0.040074   0.048374  -0.828 0.407430    
DBF.1       -0.136684   0.082271  -1.661 0.096636 .  
WNR.1        0.074055   0.029644   2.498 0.012485 *  
UFE.1       -0.007523   0.030277  -0.248 0.803779    
BPC.1        0.393523   0.102120   3.854 0.000116 ***
BPW.1        0.313412   0.069018   4.541 5.60e-06 ***
NPA.1        0.012960   0.042315   0.306 0.759398    
NPW.1       -0.014493   0.036409  -0.398 0.690578    
FSP.2       -0.065368   0.066044  -0.990 0.322288    
FSW.2       -0.263955   0.047725  -5.531 3.19e-08 ***
SSP.2              NA         NA      NA       NA    
SSW.2       -0.247359   0.063257  -3.910 9.22e-05 ***
ACE.2        0.026067   0.038823   0.671 0.501937    
DBF.2        0.117291   0.084828   1.383 0.166759    
WNR.2       -0.032345   0.028113  -1.151 0.249925    
UFE.2       -0.024302   0.028906  -0.841 0.400490    
BPC.2       -0.329055   0.111769  -2.944 0.003239 ** 
BPW.2       -0.351452   0.087206  -4.030 5.57e-05 ***
NPA.2       -0.036414   0.059266  -0.614 0.538935    
NPW.2        0.047514   0.042770   1.111 0.266603    
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 680.42  on 490  degrees of freedom
Residual deviance: 150.87  on 468  degrees of freedom
AIC: 196.87

Number of Fisher Scoring iterations: 8

> table(dw$Result,round(predict(dm.glm,newdata=dw[,-1],type='response'),0))
# confusion matrixはtable関数で計算する
# round関数を入れているのはロジスティック回帰の出力が確率値であるため
   
      0   1
  0 211  16
  1  17 208
Warning message:
In predict.lm(object, newdata, se.fit, scale = 1, type = ifelse(type ==  :
  prediction from a rank-deficient fit may be misleading
> sum(diag(table(dw$Result,round(predict(dm.glm,newdata=dw[,-1],type='response'),0))))/nrow(dw)
[1] 0.9269912 # ロジスティック回帰のテスト正答率は92.7%
Warning message:
In predict.lm(object, newdata, se.fit, scale = 1, type = ifelse(type ==  :
  prediction from a rank-deficient fit may be misleading


それなりの結果になりました。今度はSVM

> library(e1071)
 要求されたパッケージ class をロード中です 
> dm.tune<-tune.svm(Result~.,data=dm)
> dm.tune$best.model

Call:
best.svm(x = Result ~ ., data = dm)


Parameters:
   SVM-Type:  C-classification 
 SVM-Kernel:  radial 
       cost:  1 
      gamma:  0.04166667 

Number of Support Vectors:  225

> dm.svm<-svm(Result~.,dm,cost=1,gamma=dm.tune$best.model$gamma)
> table(dw$Result,predict(dm.svm,newdata=dw[,-1]))
   
      0   1
  0 210  17
  1  19 206
> sum(diag(table(dw$Result,predict(dm.svm,newdata=dw[,-1]))))/nrow(dw)
[1] 0.920354 # 正答率92.0%


お、意外とSVMのパフォーマンスがロジスティック回帰を下回る結果になりました。最後はランダムフォレスト。

> tuneRF(dm[,-1],dm[,1],doBest=T)
mtry = 4  OOB error = 12.83% 
Searching left ...
mtry = 2 	OOB error = 13.03% 
-0.01587302 0.05 
Searching right ...
mtry = 8 	OOB error = 11.61% 
0.0952381 0.05 
mtry = 16 	OOB error = 11.41% 
0.01754386 0.05 

Call:
 randomForest(x = x, y = y, mtry = res[which.min(res[, 2]), 1]) 
               Type of random forest: classification
                     Number of trees: 500
No. of variables tried at each split: 16

        OOB estimate of  error rate: 10.59%
Confusion matrix:
    0   1 class.error
0 228  23  0.09163347
1  29 211  0.12083333
> dm.rf<-randomForest(Result~.,dm,mtry=16)
> importance(dm.rf)
      MeanDecreaseGini
FSP.1         2.067423
FSW.1         3.936941
SSP.1         1.711354
SSW.1         2.888828
ACE.1         3.007018
DBF.1         2.894268
WNR.1         2.918830
UFE.1         2.007570
BPC.1        70.118502
BPW.1        13.768396
NPA.1         2.255394
NPW.1         2.306841
FSP.2         2.303675
FSW.2         8.031076
SSP.2         2.138432
SSW.2         4.207970
ACE.2         3.311075
DBF.2         3.347946
WNR.2         3.580537
UFE.2         2.624922
BPC.2        85.534696
BPW.2        14.580636
NPA.2         2.474571
NPW.2         2.958757
> table(dw$Result,predict(dm.rf,newdata=dw[,-1]))
   
      0   1
  0 219   8
  1  28 197
> sum(diag(table(dw$Result,predict(dm.rf,newdata=dw[,-1]))))/nrow(dw)
[1] 0.920354 # 正答率92.0%


以前も見たようにランダムフォレストとSVMとでパフォーマンスが同じになりました。ということで、ロジスティック回帰がベストパフォーマンスの分類器(正答率92.7%)となると同時に、算出された変数重要度からBPC.1(プレイヤー1が得たブレーク機会数)が、ロジスティック回帰の偏回帰係数&ランダムフォレストの変数重要度のいずれから見ても、最も勝敗を左右するらしい、ということが分かりました。


6. 厳密性に拘るステージ


ここから先はやってもやらなくても同じとか、厳密性を期すならこういうことをやってもいいのではないかとかそういう話です。例えば今回最も気にしなければならないのは多重共線性(マルチコ:multi-colinearityの略)で、ロジスティック回帰ならステップワイズ法で変数を絞りたいところ。この場合は{MASS}パッケージのstepAIC関数で、AICに基づいて変数を絞り込むのが定番かと。

> library(MASS)
> dm.glm.step<-stepAIC(dm.glm)
Start:  AIC=196.87
Result ~ FSP.1 + FSW.1 + SSP.1 + SSW.1 + ACE.1 + DBF.1 + WNR.1 + 
    UFE.1 + BPC.1 + BPW.1 + NPA.1 + NPW.1 + FSP.2 + FSW.2 + SSP.2 + 
    SSW.2 + ACE.2 + DBF.2 + WNR.2 + UFE.2 + BPC.2 + BPW.2 + NPA.2 + 
    NPW.2

# 中略

Step:  AIC=185.86
Result ~ FSP.1 + FSW.1 + SSW.1 + DBF.1 + WNR.1 + BPC.1 + BPW.1 + 
    FSP.2 + FSW.2 + SSW.2 + DBF.2 + UFE.2 + BPC.2 + BPW.2 + NPW.2

        Df Deviance    AIC
<none>       153.86 185.86
- DBF.2  1   155.88 185.88
- FSP.2  1   156.11 186.11
- NPW.2  1   156.16 186.16
- DBF.1  1   156.25 186.25
- FSP.1  1   158.54 188.54
- UFE.2  1   161.87 191.87
- WNR.1  1   165.91 195.91
- SSW.2  1   176.82 206.82
- SSW.1  1   178.59 208.59
- BPW.2  1   180.71 210.71
- FSW.1  1   181.31 211.31
- BPC.2  1   182.61 212.61
- BPC.1  1   187.76 217.76
- BPW.1  1   189.80 219.80
- FSW.2  1   200.18 230.18
> summary(dm.glm.step)

Call:
glm(formula = Result ~ FSP.1 + FSW.1 + SSW.1 + DBF.1 + WNR.1 + 
    BPC.1 + BPW.1 + FSP.2 + FSW.2 + SSW.2 + DBF.2 + UFE.2 + BPC.2 + 
    BPW.2 + NPW.2, family = binomial, data = dm)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.8678  -0.1372  -0.0069   0.1008   3.0421  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.67573    5.02335  -0.135  0.89299    
FSP.1        0.12962    0.06182   2.097  0.03602 *  
FSW.1        0.18332    0.03890   4.713 2.44e-06 ***
SSW.1        0.29546    0.06685   4.420 9.87e-06 ***
DBF.1       -0.11832    0.07914  -1.495  0.13490    
WNR.1        0.05049    0.01536   3.287  0.00101 ** 
BPC.1        0.38807    0.08934   4.344 1.40e-05 ***
BPW.1        0.32564    0.06178   5.271 1.36e-07 ***
FSP.2       -0.09027    0.06117  -1.476  0.14003    
FSW.2       -0.25356    0.04424  -5.732 9.94e-09 ***
SSW.2       -0.26049    0.06132  -4.248 2.16e-05 ***
DBF.2        0.11584    0.08149   1.422  0.15516    
UFE.2       -0.04035    0.01475  -2.735  0.00624 ** 
BPC.2       -0.39068    0.09954  -3.925 8.67e-05 ***
BPW.2       -0.31849    0.06928  -4.598 4.28e-06 ***
NPW.2        0.02078    0.01396   1.488  0.13665    
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 680.42  on 490  degrees of freedom
Residual deviance: 153.86  on 475  degrees of freedom
AIC: 185.86

Number of Fisher Scoring iterations: 8

> table(dw$Result,round(predict(dm.glm.step,newdata=dw[,-1],type='response'),0))
   
      0   1
  0 213  14
  1  17 208
> sum(diag(table(dw$Result,round(predict(dm.glm.step,newdata=dw[,-1],type='response'),0))))/nrow(dw)
[1] 0.9314159 # 正答率93.1%


不要な説明変数を削除したことで、確かにロジスティック回帰の分類性能が向上しました。


一方、機械学習分類器であれば次元削減したいので、例えばPCAでまとめられる説明変数を探すという手があります。で、これをやってみると。。。

> dm.pc<-princomp(scale(dm[,-1]))
> summary(dm.pc)
Importance of components:
                          Comp.1    Comp.2    Comp.3     Comp.4     Comp.5     Comp.6     Comp.7     Comp.8     Comp.9    Comp.10
Standard deviation     2.4711108 1.9169025 1.5689050 1.38949782 1.28753087 1.24644013 1.17840200 0.99883671 0.97008773 0.80688314
Proportion of Variance 0.2549521 0.1534173 0.1027703 0.08061018 0.06921329 0.06486598 0.05797772 0.04165462 0.03929128 0.02718288
Cumulative Proportion  0.2549521 0.4083694 0.5111396 0.59174980 0.66096309 0.72582907 0.78380679 0.82546141 0.86475269 0.89193557
                          Comp.11    Comp.12    Comp.13    Comp.14     Comp.15     Comp.16     Comp.17     Comp.18     Comp.19     Comp.20
Standard deviation     0.78618564 0.69539057 0.59254087 0.54542966 0.424759371 0.396639175 0.362198898 0.323533789 0.285738461 0.269436122
Proportion of Variance 0.02580622 0.02018979 0.01465922 0.01242086 0.007532864 0.006568488 0.005477324 0.004370322 0.003408879 0.003030999
Cumulative Proportion  0.91774179 0.93793157 0.95259079 0.96501165 0.972544516 0.979113004 0.984590327 0.988960650 0.992369529 0.995400528
                           Comp.21     Comp.22      Comp.23 Comp.24
Standard deviation     0.246627290 0.222120451 1.449068e-08       0
Proportion of Variance 0.002539548 0.002059924 8.767008e-18       0
Cumulative Proportion  0.997940076 1.000000000 1.000000e+00       1
> biplot(dm.pc)

f:id:TJO:20150220125123p:plain


これ、データの性質を考えれば当たり前なんですがプレイヤー1/2それぞれの説明変数って結構マルチコが強いはずなんですよね(鏡像的な関係にあるので)。累積寄与度を見ると、10個の説明変数で全体の90%弱、12個の説明変数で全体の95%弱が説明できてしまうことが分かります。ということで、biplotの結果を見て試しにFSP.1, FSP.2, WNR.1, NPW.1, NPW.2, FSW.2, ACE.2, SSW.1, DBF.1, BPW.1, SSP.1, SSP.2の12変数に絞ってみます。

> dm2<-data.frame(Result=dm$Result,FSP.1=dm$FSP.1, FSP.2=dm$FSP.2, WNR.1=dm$WNR.1, NPW.1=dm$NPW.1, NPW.2=dm$NPW.2, FSW.2=dm$FSW.2, ACE.2=dm$ACE.2, SSW.1=dm$SSW.1, DBF.1=dm$DBF.1, BPW.1=dm$BPW.1, SSP.1=dm$SSP.1, SSP.2=dm$SSP.2)
> dm2$Result<-as.factor(dm2$Result)
> tuneRF(dm2[,-1],dm2[,1],doBest=T)
mtry = 3  OOB error = 23.42% 
Searching left ...
mtry = 2 	OOB error = 24.85% 
-0.06086957 0.05 
Searching right ...
mtry = 6 	OOB error = 24.64% 
-0.05217391 0.05 

Call:
 randomForest(x = x, y = y, mtry = res[which.min(res[, 2]), 1]) 
               Type of random forest: classification
                     Number of trees: 500
No. of variables tried at each split: 3

        OOB estimate of  error rate: 24.85%
Confusion matrix:
    0   1 class.error
0 189  62    0.247012
1  60 180    0.250000
> dm2.rf<-randomForest(Result~.,dm2,mtry=3)
> dw2<-data.frame(Result=dw$Result,FSP.1=dw$FSP.1, FSP.2=dw$FSP.2, WNR.1=dw$WNR.1, NPW.1=dw$NPW.1, NPW.2=dw$NPW.2, FSW.2=dw$FSW.2, ACE.2=dw$ACE.2, SSW.1=dw$SSW.1, DBF.1=dw$DBF.1, BPW.1=dw$BPW.1, SSP.1=dw$SSP.1, SSP.2=dw$SSP.2)
> table(dw2$Result,predict(dm2.rf,newdata=dw2[,-1]))
   
      0   1
  0 169  58
  1  48 177
> sum(diag(table(dw2$Result,predict(dm2.rf,newdata=dw2[,-1]))))/nrow(dw2)
[1] 0.7676991 # え?76.8%???


うわ、超絶パフォーマンスが悪化したorz biplotの結果から目視で変数選択するのは良くないようです。そこで諦めてstepAICのベストモデルに従ってみることにします。

> dm3<-data.frame(Result=dm$Result,FSP.1=dm$FSP.1, FSW.1=dm$FSW.1, SSW.1=dm$SSW.1, DBF.1=dm$DBF.1, WNR.1=dm$WNR.1, BPC.1=dm$BPC.1, BPW.1=dm$BPW.1, FSP.2=dm$FSP.2, FSW.2=dm$FSW.2, SSW.2=dm$SSW.2, DBF.2=dm$DBF.2, UFE.2=dm$UFE.2, BPC.2=dm$BPC.2, BPW.2=dm$BPW.2, NPW.2=dm$NPW.2)
> dw3<-data.frame(Result=dw$Result,FSP.1=dw$FSP.1, FSW.1=dw$FSW.1, SSW.1=dw$SSW.1, DBF.1=dw$DBF.1, WNR.1=dw$WNR.1, BPC.1=dw$BPC.1, BPW.1=dw$BPW.1, FSP.2=dw$FSP.2, FSW.2=dw$FSW.2, SSW.2=dw$SSW.2, DBF.2=dw$DBF.2, UFE.2=dw$UFE.2, BPC.2=dw$BPC.2, BPW.2=dw$BPW.2, NPW.2=dw$NPW.2)
> dm3$Result<-as.factor(dm3$Result)
> dw3$Result<-as.factor(dw3$Result)

# SVMの場合

> dm3.tune<-tune.svm(Result~.,data=dm3)
> dm3.tune$best.model

Call:
best.svm(x = Result ~ ., data = dm3)


Parameters:
   SVM-Type:  C-classification 
 SVM-Kernel:  radial 
       cost:  1 
      gamma:  0.06666667 

Number of Support Vectors:  205

> dm3.svm<-svm(Result~.,dm3,cost=1,gamma=dm3.tune$best.model$gamma)
> table(dw3$Result,predict(dm3.svm,dw3[,-1]))
   
      0   1
  0 210  17
  1  18 207
> sum(diag(table(dw3$Result,predict(dm3.svm,dw3[,-1]))))/nrow(dw3)
[1] 0.9225664 # SVMの正答率は92.3%に向上した

# ランダムフォレストの場合

> tuneRF(dm3[,-1],dm[,1],doBest=T)
mtry = 3  OOB error = 12.63% 
Searching left ...
mtry = 2 	OOB error = 13.24% 
-0.0483871 0.05 
Searching right ...
mtry = 6 	OOB error = 9.16% 
0.2741935 0.05 
mtry = 12 	OOB error = 9.98% 
-0.08888889 0.05 

Call:
 randomForest(x = x, y = y, mtry = res[which.min(res[, 2]), 1]) 
               Type of random forest: classification
                     Number of trees: 500
No. of variables tried at each split: 6

        OOB estimate of  error rate: 10.39%
Confusion matrix:
    0   1 class.error
0 227  24  0.09561753
1  27 213  0.11250000
> dm3.rf<-randomForest(Result~.,dm3,mtry=6)
> table(dw3$Result,predict(dm3.rf,dw3[,-1]))
   
      0   1
  0 222   5
  1  26 199
> sum(diag(table(dw3$Result,predict(dm3.rf,dw3[,-1]))))/nrow(dw3)
[1] 0.9314159 # ランダムフォレストも正答率93.1%になってロジスティック回帰に追い付いた


最終的にAICに基づいて変数選択を行ったロジスティック回帰のテスト正答率と、同じ変数選択を採用したランダムフォレストのテスト正答率が揃って93.1%と、ひとまずのベストパフォーマンスに到達しました。


ということで、こういう試行錯誤をしながらより良いモデルを探していきましょうというのがこのステージです。なお、テニスのデータに限って言えばこれにさらにL1/L2正則化を入れるとか、TrueSkill Modelを共変量を伴うものに発展させてモデリングするとか、学習データが少ないので微妙ですがDeep Learningにぶち込むとか、工夫は色々出来るかと思います。


最後に


とりあえず、僕がこれまで色々な現場で見てきた「データ分析」の「アプローチ」の違いを、テニス四大大会のデータに基づいて再現してみました。こんな感じの使い分けが現場ではRによるアドホック分析に限らず、例えばPython, Javaなどでバックエンドで実装しているような場面でも行われていると思っていただければと*4


おまけ


言った手前やってみたくなったので、{h2o}でDeep Belief Netにぶち込んでみました。

> library(h2o)
> write.table(dm3,file='dm3.csv',quote=F,col.names=T,row.names=F,sep=',')
> write.table(dw3,file='dw3.csv',quote=F,col.names=T,row.names=F,sep=',')
> localH2O <- h2o.init(ip = "localhost", port = 54321, startH2O = TRUE, nthreads=-1)
> dmData<-h2o.importFile(localH2O,path='dm3.csv')
> dwData<-h2o.importFile(localH2O,path='dw3.csv')

# 中略(山のようなパラメータチューニングの試行錯誤)

> dm.dl<-h2o.deeplearning(x=2:16,y=1,data=dmData,activation="Rectifier",hidden=c(100,100,50),epochs=100)
> dm.dl.prd<-h2o.predict(object=dm.dl,newdata=dwData)
> dm.dl.prd.df<-as.data.frame(dm.dl.prd)
> table(dw3$Result,dm.dl.prd.df[,1])
   
      0   1
  0 217  10
  1  20 205
> sum(diag(table(dw3$Result,dm.dl.prd.df[,1])))/nrow(dw3)
[1] 0.9336283 # テスト正答率93.4%


一応*5、ベストパフォーマンスが出ました。Deep Learningの売りでもあるDropoutは入れてないし、Deep Belief NetならRBMかAutoencoder入れるべきと思いながらどっちも入ってませんが。。。


追記



L1/L2正則化を{glmnet}でやってみた話を書いたんですが、その最後の方がこの記事の追記になっているのでよろしければこちらもお読みくださいー。

*1:1stサーブが重要とか、ブレークポイントの意味合いとか

*2:そしてRは勝手に最初からWelchを選択する

*3:そもそも全ての対で対比較したまま有意水準を0.05に固定するのは多重比較の問題も抱える

*4:ただしバックエンド側のチューニングはもっとしんどい印象が。。。

*5:散々パラメータチューニングした後なので本当に「一応」ですがorz