とりあえず決定木モデルを実行のためのメモ

Rでの決定木分析（分類木、回帰木）の実行に関して、こうではないかとのメモを記す。

CARTアルゴリズムによる決定木分析を行うパッケージはrpartとmvpartがある。

mvpartはrpartに機能拡張を加えた上位パッケージとのこと。

2重に読み込んだ場合、以下のメッセージが出力される。
library(rpart)
library(mvpart)

次のパッケージを付け加えます: 'mvpart'

The following object(s) are masked from 'package:rpart':

    meanvar, na.rpart, path.rpart, plotcp, post, printcp, prune, prune.rpart,
    rpart, rpart.control, rpconvert, rsq.rpart, snip.rpart, xpred.rpart

決定木モデルの作成にはrpart()関数を用い、plot()やplotcp()関数などで結果を確認する。mvpartパッケージを使用した場合（上書きでも、単独でも）、rpartパッケージの出力（主に見た目）が異なる箇所がある。
上のメッセージはrpartパッケージで定義されている（表示されている）各関数が、後で読み込んだmvpartパッケージ中にも含まれています、ということだと思う。
それで、後から読み込んだものが適用されるのだと思われる。

が、mvpartのマニュアルとrpartのマニュアルの中のそれらの関数の説明は内容が違うように見えない…。

決定木モデルは、分岐数の多さによってモデルの複雑さが決定される。

回帰木（目的変数が連続値）では；
グループAの平方和とグループBの平方和の和が最小になるように、分岐に使用する変数と分割点を決定する。

分類木（目的変数がカテゴリカル）では；
ジニ係数（Gini Index）、もしくはエントロピー（Entropy, Cross-entropy）を最大化するように、分岐に使用する変数と分割点を決定する。
rpart()関数の引数にparms=list(split="gini")と指定するとジニ係数が用いられ、parms=list(split="information")とするとエントロピーが用いられるようである。
また、Loss行列なども指定できるらしい。

どこまで細かく分割するかはrpart() 関数の引数で指定できる。

CARTは、指定した条件まで木を一旦分割し、その後適当な大きさに剪定し、最終的なモデルを決定する。

データの作成。
y <- rnorm(100)
x1 <- y + 0.5*rnorm(100)
x2 <- y + 0.75*rnorm(100)
x3 <- y + rnorm(100)

Data1 <- data.frame(y,x1,x2,x3)
cor(Data1)

　
           y        x1        x2        x3
y  1.0000000 0.8810408 0.8086425 0.6582551
x1 0.8810408 1.0000000 0.6947170 0.5983025
x2 0.8086425 0.6947170 1.0000000 0.5468263
x3 0.6582551 0.5983025 0.5468263 1.0000000

yとの相関がx1>x2>x3となるデータを意図的に作成した。

決定木（回帰木）の実行。
result1 <- rpart( y~., data=Data1, control=rpart.control(minsplit=5, cp=0.01) )

control=rpart.control(minsplit=5, cp=0.01)において、どこまで細かく分岐するかを決定する。

minsplitはマニュアルによると、その枝にあるデータ数が指定した数以下になると、分岐が行われない。
cp(complexity parameter)は、木のサイズとあてはまりの良さをトレードオフする指標で、要するに、木がある大きさのときのモデルの複雑さの指標のようなものと考えられる。小さく指定すると細かく分岐が行われる。詳細はよく分かっていない。

result1

　
n= 100 

node), split, n, deviance, yval
      * denotes terminal node

 1) root 100 82.74123000 -0.08323216  
   2) x1< -0.118119 50 20.52873000 -0.76585870  
     4) x2< -1.121383 18  5.31604000 -1.36086100  
       8) x1< -1.620781 5  0.99465230 -2.05751900  
        16) x1< -2.258496 2  0.05018492 -2.55654600 *
        17) x1>=-2.258496 3  0.11437500 -1.72483500 *
       9) x1>=-1.620781 13  0.96139660 -1.09291600 *
     5) x2>=-1.121383 32  5.25564900 -0.43116960  
      10) x2< 0.3988502 27  3.16579100 -0.53866160  
        20) x1< -0.5092024 12  1.31139200 -0.73501100 *
        21) x1>=-0.5092024 15  1.02165200 -0.38158210 *
      11) x2>=0.3988502 5  0.09323544  0.14928700 *
   3) x1>=-0.118119 50 15.61460000  0.59939430  
     6) x1< 1.336195 40  8.11672700  0.42480680  
      12) x2< 1.20844 34  4.89382400  0.32000810  
        24) x3< 0.4031585 21  1.99932900  0.16314720 *
        25) x3>=0.4031585 13  1.54309500  0.57339870 *
      13) x2>=1.20844 6  0.73347440  1.01866600 *
     7) x1>=1.336195 10  1.40171400  1.29774500 *

詳細な結果はsummary(result1)で確認できる。

10のグループに分けられた。

木のプロットを書くと理解しやすい。
plot(result1)
text(result1, use.n=TRUE, cex=0.5)
引数cexで、ラベルの大きさを調整できる。

f:id:High_School_Student:20130703193216j:plain

ここで、path.rpart(result1)を実行すると、上プロット上にマウスを持っていくと＋マークになる。ノードの箇所でクリックすると分岐履歴が確認できるので面白い。

また、分岐にどの変数が多く使われているかを見るとやはりx1である、これはyとの相関をx1>x2>x3としたデータであるためである。

変数重要度のような指標を出力する関数があるかは、不明。

モデルの複雑性を確認し、剪定する箇所を決定し、最終的なモデルを作成する。

cp統計量を確認。
printcp(result1)

　
Regression tree:
rpart(formula = y ~ ., data = Data1, control = rpart.control(minsplit = 5, 
    cp = 0.01))

Variables actually used in tree construction:
[1] x1 x2 x3

Root node error: 82.741/100 = 0.82741

n= 100 

         CP nsplit rel error  xerror     xstd
1  0.563176      0   1.00000 1.01980 0.140235
2  0.120340      1   0.43682 0.46874 0.069506
3  0.073677      2   0.31648 0.40655 0.048833
4  0.040608      3   0.24281 0.37062 0.046885
5  0.030087      4   0.20220 0.36497 0.046381
6  0.024131      5   0.17211 0.34970 0.045379
7  0.016333      6   0.14798 0.30970 0.039702
8  0.010064      7   0.13165 0.28018 0.038337
9  0.010032      8   0.12158 0.27143 0.036761
10 0.010000      9   0.11155 0.27143 0.036761

プロットも行う。
plotcp(result1)

f:id:High_School_Student:20130703193327j:plain
k分割交差検証法による結果である様子。ちなみに、kの数は、control=rpart.control(xval=k)で指定できるそう。

結果の出力はいまいち不明。XerrorはRMSE (Root Mean Square Error：誤差の標準偏差)かな、と思うのだが、xstdがなにかも分からない。rel errorもいまいち不明..。

プロットは、黒い下のラインが結果の出力のrel errorと一致する。上の青いラインが検証用データによるものだと推測される。

縦軸（X-Val Relative Error）が0.3の箇所に引かれている横線がMin+1SEとのことで、最大分岐モデルのxerror+1*xstdで、0.27143 + 0.036761 = 0.308191とのこと。この横線の下側に最も近い検証用データによる結果の分岐数にモデルを決定するのが基準とのこと。詳細はますます不明..。

では、cp=0.010032、分岐数が8にモデルを剪定する。

prune() 関数を用いるとのこと。
result2 <- prune(result1, cp=0.01005)

plot(result2)
text(result2, use.n=TRUE, cex=0.5)

f:id:High_School_Student:20130703193401j:plain