【Rによるデータサイエンス】樹木モデル（続き） 〜 回帰木
「【Rによるデータサイエンス】樹木モデル」の続き。
回帰木のケーススタディについて勉強する。
◎使用するデータ
車のスピードと停止距離を表すcarsデータを使用する。
speed dist 1 4 2 2 4 10 3 7 4 4 7 22 5 8 16
１．回帰木の作成
目的変数を距離(dist)とし、説明変数をスピード(speed)として関数rpartを使用して回帰木を作成する。
> (cars.rp<-rpart(dist~speed,data=cars)) n= 50 node), split, n, deviance, yval * denotes terminal node 1) root 50 32538.9800 42.98000 2) speed< 17.5 31 8306.7740 29.32258 4) speed< 12.5 15 1176.4000 18.20000 8) speed< 9.5 6 277.3333 10.66667 * 9) speed>=9.5 9 331.5556 23.22222 * 5) speed>=12.5 16 3535.0000 39.75000 * 3) speed>=17.5 19 9015.6840 65.26316 6) speed< 23.5 14 2846.8570 55.71429 12) speed>=18.5 10 1323.6000 52.20000 * 13) speed< 18.5 4 1091.0000 64.50000 * 7) speed>=23.5 5 1318.0000 92.00000 *
これから回帰木は６つの葉を持っていることが分かる。らしいが、私にはどこを見れば葉の数が分かるのか分からない！(追記)※がついているのが葉。コメントでr-de-rさんに教えてもらいました。
なので、次のplotcpの実行によって得られるグラフの上の横軸(Size of tree)を見て葉の数を知ることにする。
２．回帰木の剪定
剪定の必要性についての考察。
> plotcp(cars.rp)
成長させた回帰木を４つの葉に剪定する必要があることを示唆している。
そういうわけで、回帰木を４つの葉に剪定する。そのためには、関数pruneを使用する！
> cars.rp1<-prune(cars.rp,cp=0.044) > plot(cars.rp1,uniform=T,margin=0.05) > text(cars.rp1,use.n=T)
相変わらず期待する図を描くことができない・・・。
【Rによるデータサイエンス】樹木モデル（続き） 〜 多変量回帰木
「Rによるデータサイエンス】樹木モデル（続き） 〜 回帰木」では回帰木のケーススタディを実施した。
今回は多変量回帰木のケーススタディを行う。
なお、多変量回帰とは、目的変数が複数である回帰分析のことを呼ぶ。
◎使用するデータ
パッケージmvpartに同梱されているデータspiderを使用する。
> nrow(spider);ncol(spider) [1] 28 [1] 18
28行、18列のデータである。
18列のうち、左の12列は異なる蜘蛛の種類で、残りの6列は蜘蛛狩りの環境に関するデータ。
１．データをマトリックス形式に変換する
２．mvpart関数を使用して多変量回帰木を作成する
> spider[1,13:18] water sand moss reft twigs herbs 1 9 0 1 1 9 5 > spider.mv<-mvpart(as.matrix(spider[,1:12])~water+sand+moss+reft+twigs+herbs,data=spider)
実行すると回帰木が作成される。
次に、目的変数の主成分分析の結果も表示させるために、引数pcaにTRUEを指定する。
> spider.mv<-mvpart(as.matrix(spider[,1:12])~water+sand+moss+reft+twigs+herbs,data=spider,pca=T)
グラフを一度クリックすると、目的変数の主成分分析のグラフが表示される。
さて、同じことを関数rpartを使用して実施する。
１．rpart関数を使用して多変量回帰木を作成する
> spider.rp<-rpart(as.matrix(spider[,1:12])~water+sand+moss+reft+twigs+herbs,data=spider)
可視化する。
> plot(spider.rp,uniform=T,margin=0.05) > text(spider.rp,use.n=T)
とりあえず樹木を描くことができたが、全ての葉が必要かを確認する（剪定の必要性を確認する）ために関数plotcpを実行する。
> plotcp(spider.rp)
> printcp(spider.rp) rpart(formula = as.matrix(spider[, 1:12]) ~ water + sand + moss + reft + twigs + herbs, data = spider) Variables actually used in tree construction: [1] herbs moss reft sand water Root node error: 1418/28 = 50.642 n= 28 CP nsplit rel error xerror xstd 1 0.518641 0 1.00000 1.06149 0.121164 2 0.144890 1 0.48136 0.56700 0.072972 3 0.075375 2 0.33647 0.43982 0.074040 4 0.046781 3 0.26109 0.41701 0.075281 5 0.035301 4 0.21431 0.41637 0.071690 6 0.034792 5 0.17901 0.43751 0.072407 7 0.018806 6 0.14422 0.42978 0.077079 8 0.010000 7 0.12541 0.45108 0.077910
ここまでの結果から、葉（Size of tree)は3つで十分だということが分かるので、No3剪定して可視化する。
> spider.rp1<-prune(spider.rp,cp=0.075375) > plot(spider.rp1,uniform=T,branch=0.6) > text(spider.rp1,use.n=T)
期待通りに剪定された。ただし、ラベルをつけてくれていないなど、mvpartの方が一見便利そう。可視化方法は別途学習が必要。。。
【Rによるデータサイエンス】樹木モデル
【定義】
樹木モデルとは。。。
非線形回帰分析、非線形判別分析の1つの方法である。
回帰問題では回帰木（regression tree)、分類問題では分類木(classification tree)あるいは決定木(decision tree)と呼ばれている。
これだけでは抽象的すぎて、分からない。
何をインプットとして？どのようにして？どういう形で出力してくれるのか。
「何を？」という部分は、ざっくり言って「データ」なので、「どのようにして？」という部分が気になる。
◎樹木モデルの基礎
樹木モデルでは、CHAID、C4.5/C5.0/See5、CARTをベースとしたアルゴリズムが広く用いられている。
| アルゴリズム名 | 説明 | 分岐基準 |
| CHAID | Hartiganが1975年に、Morgannらによって提案されたAIDを発展させたもの。 | カイ2乗統計量、F統計量 |
| C4.5/C5.0/See5 | オーストラリアのJ.ROSS Quanlanが機械学習のアプローチで1986年に発表したID3を改良・発展させたもの。 | 利得比 |
| CART | 2進再帰分割を行い説明変数を2進木に分岐させる。 | ジニ係数、情報利得 |
◎ケーススタディ 〜 分類木
１．木の作成
関数rpartを使用して分類木を作成する。
rpartでは、分岐点の計算はジニ係数とエントロピーを用いる。デフォルトにはジニ係数が使用される。エントロピーを使用する場合は、引数split="information"を用いる。
> library(mvpart) > set.seed(20) > iris.rp<-rpart(Species~.,data=iris) > iris.rp n= 150 node), split, n, loss, yval, (yprob) * denotes terminal node 1) root 150 100 setosa (0.33333333 0.33333333 0.33333333) 2) Petal.Length< 2.45 50 0 setosa (1.00000000 0.00000000 0.00000000) * 3) Petal.Length>=2.45 100 50 versicolor (0.00000000 0.50000000 0.50000000) 6) Petal.Width< 1.75 54 5 versicolor (0.00000000 0.90740741 0.09259259) 12) Petal.Length< 4.95 48 1 versicolor (0.00000000 0.97916667 0.02083333) * 13) Petal.Length>=4.95 6 2 virginica (0.00000000 0.33333333 0.66666667) * 7) Petal.Width>=1.75 46 1 virginica (0.00000000 0.02173913 0.97826087)
可視化する。
plot関数の引数は次のとおりとする。
- uniform引数にTを設定し、ノードが等間隔になるようにする。デフォルト（uniform=FALSE)では、木のノードの間の間隔が分類のエラーの数に比例するようになる。
- branch引数にて、枝の角度を調整する。指定する値は0から1の間の数で、0の場合が最も角度が大きい。1の場合は垂直となる。
- margin引数は、図の外側の余白を調整する。
> plot(iris.rp,uniform=T,branch=0.6,margin=0.05) > text(iris.rp,use.n=T,all=T)
どういうわけか図が期待通りに書けない・・・。
２．木の剪定
次に、交差確認法の結果を確認する。関数rpartでは、樹木を成長させると同時に交差確認法の結果も計算している。
関数printcpは、樹木の剪定のための複雑さのパラメータ(cp)を返す。
> printcp(iris.rp) Classification tree: rpart(formula = Species ~ ., data = iris) Variables actually used in tree construction: [1] Petal.Length Petal.Width Root node error: 100/150 = 0.66667 n= 150 CP nsplit rel error xerror xstd 1 0.50 0 1.00 1.19 0.049592 2 0.44 1 0.50 0.68 0.060970 3 0.02 2 0.06 0.11 0.031927 4 0.01 3 0.04 0.10 0.030551
通常、木のサイズはxerrorの最小値を中心としたその標準偏差1倍の範囲内の最大のxerror値を選ぶ。これをMin+1SE法と呼ぶことにする。
上で得られたxerrorの結果に対してMin+1SEの値を計算すると次のようになる。
| No | xerror | xstd | Min+1SE |
| 1 | 1.19 | 0.049592 | 1.239592 |
| 2 | 0.68 | 0.060970 | 0.740970 |
| 3 | 0.11 | 0.031927 | 0.141927 |
| 4 | 0.10 | 0.030551 | 0.130551 |
この結果を見ると、最終行のMin+1SEの値はNo3のxerrorの値より大きくて、No2のxerrorの値より小さい。
そのため、第3行のcp(complexity parameter)=0.02が1つの目安となる。
関数rpartでは、複雑度cpで木をコントロールすることができる。cpの値が小さいほど木が複雑になる。
デフォルトはcp=0.01となっているが、ここではrpartで生成した木を剪定する関数pruneを使用してcp=0.02で樹木を剪定する。
> iris.rp1<-prune(iris.rp,cp=0.02) > plot(iris.rp1,uniform=T,branch=0.6) > text(iris.rp1,use.n=T)
相変わらず期待通りにグラフを作成できない・・・。絵心はまた別途勉強することとする・・・。
ところで、実はrpartにはcpを視覚化する関数plotcpがある。
これを用いて剪定に必要となる交差確認に関する情報を可視化できる。
オレンジ色の点が樹木を剪定する目安である。xerrorの最小値は赤色部分。
３．判別
作成した樹木モデルを用いて未知のデータについて予測・判別を行うためには、モデルの作成に使用していないデータが必要である。
そこで、ここでもirisデータを奇数行と偶数行に分けて学習用とテスト用にする。
> plotcp(iris.rp) > plotcp(iris.rp) > even.n<-2*(1:75)-1 > iris.train<-iris[even.n,] > iris.test<-iris[-even.n,] > set.seed(20) > iris.rp2<-rpart(Species~.,iris.train) > plotcp(iris.rp2)
作成したモデルを用いて予測・判別を行う。
> iris.rp3<-predict(iris.rp2,iris.test[,-5],type="class") > table(iris.test[,5],iris.rp3) iris.rp3 setosa versicolor virginica setosa 25 0 0 versicolor 0 24 1 virginica 0 3 22
