ちょっと前の記事でこんなネタをやってみたわけですが。
世の中には色々な「データ分析」のやり方があるなぁと思った時に、この同じ2013年のテニス四大大会のデータからそれぞれのやり方をしている人たちがどんな異なるアプローチを取るのかなぁとふと想像したもので、半分ネタ的に書いてみました。便宜的に以下のようにステージを分けてあります。
なお、データは以前の記事と同じこちらのものをお使い下さい。
その上で、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も問われる妙に実践的なデータセットとも言えますねｗ
ということで、ここでは「プレイヤーの勝敗に関連する要因は何か」を追究するというのが目的であるとして、ステージ別に色々やってみることにします。
1. 集計ステージ
例えばこのデータの分析をExcelでやれと言われたら、僕が思い付くのは「Resultが0 / 1それぞれの場合で23個ある変数の平均値を取ってみる」というものでしょうか。これは以前一緒に働いたことのあるマーケッターの人たちもやっていたやり方です。
で、男子のデータだけそれでやってみました。単に"men.txt"をコピペしてResultの結果でソートし、36列それぞれの平均をAVERAGE関数で計算するだけです。ついでに棒グラフも描いてみることに。
うん、何か分かる気がする。でも何も分からない気もする。。。ちなみに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
ところどころ「勝者の方が大きい値」とか「敗者の方が大きい値」とかが見えますが、これだけではちょっと分からないことばかりだなぁという印象がありますね。ただし大まかな傾向はこれで分かります。
2. 検定ステージ
データ集計が当たり前になった現場では一般に「検定」が好まれる傾向があると個人的には感じています。特にt検定とかカイ二乗検定とかが好まれるようです。このデータセットに対しても同じように（例えば）t検定を使ってみる、ということを考える人は結構な数いると思います。
例えば、男子のACE.1のデータなんかは明らかに差がありそうなので、検定で大小を確かめてみる価値はあるんじゃないかと考えられます。ではこれをやってみようと思うわけですが。。。
うあー、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. 重回帰分析ステージ
そんなわけで検定では飽き足らなくなってきた人が興味を持ってちらほらやり始めるのが「重回帰分析」。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とそこそこ高いし、なかなか悪くないんじゃないでしょうか。ということで調子に乗って男子の結果・女子の結果それぞれに当てはめてみると。。。
> 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じゃない結果だった。。。これではちょっと使えないですね。二値の場合に適した手法を選択しなければならないようです。
4. 機械学習を含めたモデリングステージ
基本的に「テニスの勝敗」というのは二値データです。ということは、機械学習における分類器を用いるのがおそらく相性が良くて、その上で例えばどの説明変数が分類のために重要だったかを調べる、というのが妥当だろうと考えられますね。
ということで、ここまで来てようやく以前の記事同様の流れにたどり着きます。おさらいの意も込めてもう一度やってみましょう。まず、二値カテゴリ型であることを明示するために目的変数を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が得たブレーク機会数）が、ロジスティック回帰の偏回帰係数＆ランダムフォレストの変数重要度のいずれから見ても、最も勝敗を左右するらしい、ということが分かりました。
5. 厳密性に拘るステージ
ここから先はやってもやらなくても同じとか、厳密性を期すならこういうことをやってもいいのではないかとかそういう話です。例えば今回最も気にしなければならないのは多重共線性（マルチコ：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)
これ、データの性質を考えれば当たり前なんですがプレイヤー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入れるべきと思いながらどっちも入ってませんが。。。