読者です 読者をやめる 読者になる 読者になる

東京に棲む日々

データ分析、統計、ITを勉強中。未だ世に出ず。

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

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

予測は、predict(result2, 予測データ) として行う。


と、不明な点は多いのだが、とりあえずお手本に従った使い方をしてみた。

今回はCARTを使ったが、決定木モデルはCHAIDなど含めいろいろと派生型がある。
個人的には、JMPソフトウェアの「分岐」、「剪定」ボタンで実行できる対話的な方式も使い勝手が良いと思う。

決定木であるが、長所は、モデルが単純な分岐の繰り返しなので解釈が容易、と言われる。

こういった、グループ分けの解釈を決定木に求めるのであればJMPの対話形式は非常に良い。

ただし、あてはまりの良いレベルまで木を分岐させると、たいていの場合、解釈が簡単でないまで木が成長してしまう..。

それに、決定木モデルは不安定とも言われる。あてはめを行うそのデータ(学習データ)にとても依存するモデルで、外れ値などの影響で、ある変数がたまたま上位で分岐されると、その下位の分岐は、その分岐のもとに成り立ってしまうことになる。

この不安定性をバギングと分岐毎の分岐の候補変数の無作為抽出によって克服し、予測精度を高めた決定木モデルの応用がランダムフォレストである。

ランダムフォレストによる予測精度の向上だけではなく、前回(「重回帰と変数選択」)でも触れたが、沢山変数がある場合の重要な変数の抽出に関しても次回以降触れたい。

 

Rによるデータサイエンス - データ解析の基礎から最新手法まで

Rによるデータサイエンス - データ解析の基礎から最新手法まで