==== 第2週 ====

今週の{{:mselab:2013:stat:week3:note-statistics-3-20131029.pdf|配付資料}}。グラフと表のみによるデータ解析と比較的単純な学習機械によるデータ解析に続いて、今回は次の2つの課題に取り組んでもらう。

PPDACの5つの要素のうち、Pが変わって、Dにデータが追加されて、Aも変更になる。

  - P(roblem)は「V86の契約の有無と他の変数との関係の分析」から「V86の契約の有無の他の変数からの予測モデルの構築」に変わる
    - 先週は「学習機械の当てはめ」によるデータの分析、今週は「学習機械による予測結果」を評価してどれを用いるかを検討
    - 各学習機械には様々な「パラメータ」があり、それを変えることで学習結果も予測精度も変わる
    - 線形判別分析(lda)、二次判別分析(qda)、決定木(rpart)だけでなく、より柔軟な学習機械も用いる候補に含める
  - D(ata)に先週に解析してもらったtic.learnに加えて、tic.evalを追加する
    - tic.learnは学習用、tic.evalは検証用
  - A(nalysis)に「予測精度の評価」が加わる
    - 前回は学習しただけ、今回は学習したモデルを予測に用いる(仮想的な現場投入)
    - tic.learnで学習したモデルを、tic.evalの予測に用いて、予測精度を評価する
  - C(onclusion)は「どのモデルが予測精度が良いか」とその考察
  - キャンペーンの提案書の作成
    - 一番上手に予測できた「制御パラメータを調整済み」の学習機械に基づいて、予測精度を自慢しつつ、効率が良いと思われるキャンペーンを提案する提案書を起草せよ。

１回目と２回目とで、各手法に対する印象が少し異なってくると嬉しい。

=== 学習と予測 ===
== 学習 ==

例えばあるデータに回帰分析を行う。

<code>
tic.lm <- lm(V86~.,data=tic.learn)
</code>

学習理論の言葉では、これは「回帰分析」という学習機械を「最小二乗法」という学習アルゴリズムで学習させたことに相当する。
そして学習の精度は、「実際の個々の観測値(tic.learn$V86)」と、モデルを推定した結果による個別の平均の推定値である「当てはめ値(fitted(tic.lm))」とを比較して検討する。

まずは横軸を観測値、縦軸を当てはめ値にした散布図。
<code>
plot(tic.learn$V86, fitted(tic.lm),
     xlab="V86",
     ylab="Fitted Value")
</code>

あまりよくわからないので、箱ひげ図に直してみる。
<code>
plot(fitted(tic.lm)~as.factor(tic.learn$V86),
     xlab="V86",
     ylab="Fitted Value")
</code>

ヒストグラムを描いてもいい。
<code>
par(mfrow=c(2,1))
hist(fitted(tic.lm)[tic.learn$V86==0],
     main="V86=0",
     xlab="Fitted Value")
hist(fitted(tic.lm)[tic.learn$V86==1],
     main="V86=1",
     xlab="Fitted Value")
</code>

あまり変わらないように見える。
平均の差のt検定(Welchの検定)ぐらいをやってみると、この差に統計的な有意差はあることが分かる。

<code>
t.test(fitted(tic.lm)[tic.learn$V86==0],fitted(tic.lm)[tic.learn$V86==1],var.equal=F)
</code>
さりとて、当てはめ値が0付近なのはあまりよい回帰分析とは言えない。

当てはまりの良さ(悪さ)は、残差という量でも検討できる。残差の定義は、「実際の個々の観測値(tic.learn$V86)-当てはめ値(fitted(tic.lm))」であり、データの個々の観測値の変動のうちのモデルで説明しきれなかった部分である。
関数residuals()を用いて計算させることができる。
<code>
plot(residuals(tic.lm))
</code>
単に、レコード番号順に残差を計算させるとこうなる。

残差の大きさや傾向を検討をするときには、観測値と残差の傾向をみる散布図
<code>
plot(tic.learn$V86,residuals(tic.lm),
     xlab="V86",
     ylab="Residuals")
</code>
や、当てはめ値と残差の散布図
<code>
plot(fitted(tic.lm),residuals(tic.lm),
     xlab="Fitted Values",
     ylab="Residuals")
</code>
を見て、何か説明できていない傾向が残っていないかどうかを検討する。
他に、てこ比を見たり、半正規プロットなどを検討することもある。

重回帰分析に限って、上で紹介したグラフの多くは
<code>
plot(tic.lm)
</code>
で自動的に描画される。

他にも要約の出力
<code>
summary(tic.lm)
</code>
を見て、学習精度を検討する。

== 予測 ==

予測とは、学習した結果(学習済みの学習機械)を別のデータに基づく予測に用いることを言う。
そのためには、学習用データと同じ構造を持つ、別のデータを用意する。
予測精度を評価するので、ここでは評価用データと呼ぶ。
今年の実験ではtic.learnが学習用データ、tic.evalが評価用データであり、前回のデータ読み込みの手順で両方ともダウンロードされる。

予測値の算出には関数predictを用いる。学習済み機械tic.lmに、tic.evalの中の各レコードのV1からV85に基づいてV86を予測させるにはpredict()を
<code>
predict(tic.lm, newdata=tic.eval)
</code>
と用いる。学習にtic.evalを用いていないことから、この結果は予測値として扱える。

lmを用いて予測したので、
<code>
v86.lm <- predict(tic.lm, newdata=tic.eval)
</code>
と名付けておく。

<code>
plot(tic.eval$V86, v86.lm,
     xlab="V86",
     ylab="Predicted Value")
</code>

ヒストグラムを描いてもいい。
<code>
par(mfrow=c(2,1))
hist(v86.lm[tic.eval$V86==0],
     main="V86=0",
     xlab="Predicted Value")
hist(v86.lm[tic.eval$V86==1],
     main="V86=1",
     xlab="Predicted Value")
</code>
なんかダメそうに見える。

さて残差は「当てはめ値-観測値」だったが、予測の良さを測るには予測誤差を用いる。
予測誤差の定義は「予測値-正解」である。
<code>
plot(tic.eval$V86, tic.eval$V86-v86.lm,
     xlab="V86",
     ylab="Prediction Error")
</code>
予測誤差にも傾向があり、とてもではないが、重回帰分析はこのままではよい予測を与えてくれる気がしない。
箱ひげ図を描くまでもないが、念のため。

<code>
plot((tic.eval$V86-v86.lm)~as.factor(tic.eval$V86),
     xlab="V86",
     ylab="Prediction Error")
</code>

あまり変わらないように見える。
平均の差のt検定(Welchの検定)ぐらいをやってみると、この差に統計的な有意差はあることが分かる。

<code>
t.test((tic.eval$V86-v86.lm)[tic.eval$V86==0],
       (tic.eval$V86-v86.lm)[tic.eval$V86==1],
       var.equal=F)
</code>

まあだめそう。。。

== 線形判別分析(lda)の場合 ==

学習機械は予測精度が高くて初めて、有用性が認められる。そのためには「学習に用いなかったデータ」に関して予測を行い、精度を評価する必要がある。

例えば線形判別分析(lda)に学習用データ(tic.learn)を学習させるとする。

<code>
lda(V86~.,data=tic.learn)
</code>

その学習結果は次のようにオブジェクトに収めることができる。

<code>
tic.lda <- lda(V86~.,data=tic.learn)
</code>

これを用いて、評価用データ(lda.eval)に関するCARAVAN(V86)の予測を行うには、

<code>
predict(tic.lda,newdata=tic.eval)
</code>

を実行すれば良い。この結果もまたオブジェクトの中に収めることができる。

<code>
v86.lda <- predict(tic.lda,newdata=tic.eval)
</code>

v86.ldaは「リスト」と呼ばれるオブジェクトであり、

<code>
is.list(v86.lda)
</code>

を実行して、リストかどうかを尋ねるとTRUEが返ってくる。

リストの中身を見るにはnames関数を用いる。

<code>
names(v86.lda)
</code>

画面には次のように表示される。

<code>
> names(v86.lda)
[1] "class"     "posterior" "x"      
</code>

これはv86.ldaの中には次の3つが入っていることを意味する。

|オブジェクト名|中身|
|class|クラス|
|posterior|事後確率|
|x|判別関数の値|

この内訳は、学習に使用したアルゴリズム毎に異なる。

ここからが少し、問題。学習機械に基づく予測結果は
<code>
v86.lda$class
</code>
に収められている。ただしこれは因子(factor)オブジェクトであり、加減乗除ができない。
<code>
v86.lda$class+1
</code>
はエラーになる。そこで、正答率を出すには
<code>
table(v86.lda$class, tic.eval$V86)
</code>
のようにクロス集計を実行して、予測結果(v86.lda$class)と実際(tic.eval$V86)を比較する。

== 二次判別分析(qda)の場合 ==

他の学習機械についても同様に
<code>
tic.qda <- qda(V86~.,data=tic.learn)
v86.qda <- predict(tic.qda,newdata=tic.eval)
table(v86.qda$class, tic.eval$V86)
</code>
などと、予測性能の評価が可能になる。ただし二次判別分析(qda)はエラーを表示して実行できない。この時は、エラーの内容から
<code>
tic.qda <- qda(V86~V1,data=tic.learn)
tic.qda <- qda(V86~V1+V2,data=tic.learn)
tic.qda <- qda(V86~V1+V2+V3,data=tic.learn)
tic.qda <- qda(V86~V1+V2+V3+V4,data=tic.learn)
tic.qda <- qda(V86~V1+V2+V3+V4+V5,data=tic.learn)
tic.qda <- qda(V86~V1+V2+V3+V4+V5+V6,data=tic.learn)
</code>
などと、ひとつひとつ増やしていってみるか、例えば
<code>
tic.qda <- qda(V86~.,data=tic.learn[,c(1:30,86)])
</code>
のように最初の30変数を使って86番目を予測する、とするなどの工夫が必要になる。

== 決定木(rpart)の場合 ==

再帰的分割(rpart)を用いる場合は、次のようなデータの変換が必要になる。
<code>
tic.learn.2 <- tic.learn
tic.learn.2$V86 <- as.factor(tic.learn$V86)
</code>
V86を因子オブジェクト(factor)に変換すると、rpartが分類木を学習する。この変換をしないと、回帰木を学習してしまう。
<code>
tic.rpart <- rpart(V86~.,data=tic.learn.2)
v86.rpart <- predict(tic.rpart,newdata=tic.eval,type="class")
table(v86.rpart, tic.eval$V86)
</code>

rpartの制御はたとえば
<code>
tic.rpart <- rpart(V86~.,data=tic.learn.2, control=rpart.control(cp=0.00001, minsplit=5))
</code>
のように指定する。

誤判別率(誤った予測をした割合)などを算出しつつ、各学習機械の制御パラメータを調整してみよ。

=== 柔軟な学習機械を使ってみる ===

ひとつ目の課題は、より柔軟な学習機械を適用してもらい、先週までの方法との違いを検討してもらう。こちらについては今回も、同志社大学の金(じん)先生が公開されてらっしゃる

  * [[http://mjin.doshisha.ac.jp/R/31/31.html|Ｒとカーネル法・サポートベクターマシン]]
  * [[http://mjin.doshisha.ac.jp/R/32/32.html|Ｒと集団学習]]

を参考に、進めて欲しい。

サポートベクターマシン(ksvm)はkernlabパッケージに、バギング(bagging)はipredパッケージに、ブースト(ada)はadaパッケージに、ランダムフォレスト(randomForest)はrandomForestパッケージに、それぞれ含まれているので、まずは必要なパッケージをインストールする。
<code>
Sys.setenv("http_proxy"="http://130.153.8.19:8080/")
install.packages("kernlab", dependencies=TRUE)
library(kernlab)
install.packages("ipred", dependencies=TRUE)
library(ipred)
install.packages("ada", dependencies=TRUE)
library(ada)
install.packages("randomForest", dependencies=TRUE)
library(randomForest)
</code>

サポートベクターマシンに学習させる実行例：
<code>
tic.svm <- ksvm(V86~., data=tic.learn.2, 
                 kernel="rbfdot",
                 kpar=list(sigma=0.01))
v86.svm <- predict(tic.svm, newdata=tic.eval) 
table(v86.svm,tic.eval$V86)
</code>
kparの中身を変更しないと、全く予測精度が稼げない感じ。

バギングに学習させる実行例：
<code>
tic.bagging <-bagging(V86~., data=tic.learn.2,
                  nbagg=40) 
v86.bagging <- predict(tic.bagging,newdata=tic.eval,type= "class") 
table(v86.bagging, tic.eval$V86)
</code>
これもバッグの数が40でいいか分からない微妙な感じ。

Adaブーストに学習させる実行例：
<code>
tic.ada <- ada(V86~., data=tic.learn.2,
                iter=20) 
v86.ada <- predict(tic.ada, newdata=tic.eval, type= "vector") 
table(v86.ada, tic.eval$V86)
</code>
これも繰り返し数(iter)が20ではぜんぜんだめな感じ。

ランダムフォレストに学習させる実行例：
<code>
tic.randomForest <- randomForest(V86~., data=tic.learn.2) 
v86.randomForest <- predict(tic.randomForest, newdata=tic.eval) 
table(v86.randomForest, tic.eval$V86)
</code>

以上は「弱い学習機械」の指定と、個々の学習機械の設定パラメータが、デフォルトのままであることとは注意しておく。それぞれの予測結果を見れば、パラメータのチューニングが必要なことは明らか。

なお、この課題では「過学習」について言及していないが、過学習は大事な問題であるので、各自、少し調べて気にすることを進める。

=== 提案書 ===

2つ目の課題は、前回までのレポートの考察として、どのような顧客層に重点的にキャンペーンを展開するのがよいか、提案書を起草してもらう。「提案書」という書類の形式については、成書を参照するのがよいが、参考までに幾つかのウェブサイトへのリンクを掲げておく。

  * [[http://www.fujixerox.co.jp/support/xdirect/magazine/rp0811/08111a.html|ロジカル＆ビジュアルなドキュメンテーション（提案書編）]](富士ゼロックス)
  * [[http://www.bizocean.jp/doc/category/260/|提案書のテンプレート]]
  * [[http://itpro.nikkeibp.co.jp/article/COLUMN/20130123/451421/|提案書作成スキルを身につけよう]](日経BP)
  * [[http://japan.zdnet.com/blog/netcommerce/2013/07/13/entry_30022576/|提案書の最初を見れば営業の力量が見えてくる]](ZDnet)
  * [[http://www.microsoft.com/japan/office/previous/xp/suminaka/powerpoint/teian/teian_con2_1.htm|顧客への提案活動の実施]](マイクロソフト)

なおレポートには、訪問リストを生成するための手順を記述しておけばよく、訪問リストそのものの添付は不要である。また、最終的な提案に機械学習の結果を用いる場合は、以下の点に留意されたい。

  * 決定木や判別機械を学習させた結果を用いて訪問リストを作るとき、それらの学習に用いた設定や、学習結果は説明してほしい。
  * サポートベクターマシンを学習させた結果を用いて訪問リストを作るとき、その訪問リストの傾向は検討してほしい。アンサンブル学習を用いた場合も同様。
    * 例えば「サポートベクターマシンによって選ばれた顧客リスト1000名に営業をかけるのが効果的です」と宣言したら、「サポートベクターマシンって何かね」「柔軟な学習機械です」「どう柔軟なんだね」「非常に柔軟です」「・・・」というバッドエンドなコミュニケーションは可能な限り回避するのがいい。
    * 「サポートベクターマシンによって、例えば以下のような傾向を持つ顧客が対象に選ばれました。それぞれの傾向は、契約獲得の可能性を高めることが予想されています。」として選ばれた顧客と選ばれなかった顧客を、集団同士として比較するような分析を加えた方が、説得力は増す。