1．樹木モデル関数rpartとmvpart
関数rpart(recursive partitioning and regression trees)は1997年Terry TherneauとBeth Atkinsonにより実装され、関数mvpartはrpartのモデルの当てはめとプロットのための包括関数で、オーストラリア海洋研究所のGlenn De’athより実装された。
関数rpartの書き式を次に示す。より多くの引数やデフォルトの設定などはhelp(rpart)で確認することができる。
rpart(formula, method,･･･)
formulaの書き式は、treeと同じである。引数methodは、表1の中からひとつを選択することができる。引数methodが指定されていない場合は、formulaに指定されている応答変数(被説明変数)のデータ型から判断する。
表１ 引数methodのオプション
|
応答変数 |
引数の書き式 |
|
|
method="anova" |
|
|
method="poisson" |
|
|
method="class" |
|
|
method="mrt" |
|
|
method="dist" |
|
|
method="exp" |
関数rpartを使用するためにはrpart、あるいはmvpartを読み込まなければならない(library(rpart)、library(mvpart))。
関数オブジェクトrpartのクラスに属する主なメソッド関数を次の表に示す。
表2 rpartクラスの主な関数
|
関数名 |
主な機能 |
|
plot.rpart |
樹木の図示 |
|
text.rpart |
文字列の操作 |
|
plotcp |
交差確認の結果の図示 |
|
predict |
モデルによる当てはめ |
|
print.rpart |
樹木の出力 |
|
printcp |
複雑さのパラメータの出力 |
|
prune.part |
樹木の剪定を行う |
|
residuals.rpart |
残差を返す |
|
rpart.control |
樹木をコントロール |
|
summary.rpart |
要約を返す |
2．分類木
関数rpartでは、分岐の基準としてはGini係数
Gini Index=
あるは情報量エントロピー
Entrory=
が用いられている。式の中の
はノード
内のクラス
のデータである。
デフォルトではGini係数が設定されている。分岐基準エントロピーの指定は、引数split=”information”を用いる。
(1) 樹木の作成
関数rpartを用いた樹木の作成について例を用いて説明する。ここでもirisのデータを用いることにする。
>library(mvpart) #あるいはlibrary(rpart)
> iris.rp<-rpart(Species~.,data=iris)
上記のコマンドの実行で、デフォルトに設定された条件の下でrpartの分類木が作成される。作成された樹木の結果は、次のように関数printで返すことができる。関数のprintに用いたdigitは返す値の小数点以下の桁数を指定する引数である。返される結果はdigitの数値より一桁多い。
>print(iris.rp,digit=2)
この結果の樹木のグラフは、次のように関数plotとtextを用いて作成することができる。
> plot(iris.rp,uniform=T,branch=0.6,margin=0.05)
> text(iris.rp,all=T,use.n=T)
図1 rpartによるirisの分類木
関数plot.rpartには樹木のデザインのための多くの引数が用意されている。
上記のコマンドに用いたuniformはノード間の間隔に関する引数である。デフォルトはuniform=FALSEで、木のノードの間の間隔は分類のエラーの数に比例するように作成する。引数をuniform=TRULにすると、ノートの間の間隔は等間隔となる。
引数branchは枝の角度を調整する。指定する値は0から1までである。0の場合の角度が最も大きく、1の場合は垂直となる。デフォルトはbranch=1になっている。
引数marginは図と４辺との間の余白(マージン)を調整する。引数marginの値が大きいほど図が小さくなり、余白が大きい。
関数textにも複数の引数が用意されている。上記のコマンドで用いた引数use.nは、各ノードに含まれる個体の数を表示するか、しないかを指定する。デフォルトにはuse.n=FALSEになっている。引数がuse.n=TRUEである場合は、ノードに含まれる数を表示する。
(2) 樹木の剪定
関数rpartは樹木を成長させると同時に、交差確認法の結果も計算する。関数printcpでその結果を返す。デフォルトのn重交差確認のnはxval=10に設定されている。
関数printcpは、樹木の剪定のための複雑さのパラメータ(cp)を返す。この結果は、用いたデータセットをランダムに分割して交差確認を行った結果であるので、同じデータを用いて同じの条件のもとでrpartを繰り返しても、まったく同じの結果になるとは限らない。
>printcp(iris.cp)
通常、樹木のサイズはxerrorの最小値からその標準偏差1倍の範囲内の最大のxerror値を選ぶ。これをMin+1SE方法と呼ぶことにする。
上記に返された最終行のxerro、xstdはそれぞれ0.07、0.0258であり、Min+1SE＝0.07 +0.258=0.0958となる。この値は、第2行のxerror=0.7より小さく、第3行のxerror=0.07より大きいので、第3行のcp(complexity parameter)=0.02がひとつの目安となる。
関数rpartでは、複雑度cpで樹木をコントロールすることができる。cpの値が小さいほど樹木が複雑になる。デフォルトはcp=0.01になっている。
図2 iris.rpのplotcpのプロット
図3に関数pruneを用いてcp=0.02で剪定したrpartの樹木を示す。
>iris.rp1<-prune(iris.rp,cp=0.02)
> plot(iris.rp1,uniform=T,branch=0.6)
> text(iris.rp1,use.n=T)
図3 cp=0.02で剪定したirisの分類木
(3) 樹木モデルによる予測
作成した樹木のモデルを用いて未知のデータについて予測・判別を行うためには、モデルの作成に用いていないデータが必要である。そこでここでも関数treeを説明するときと同じくirisデータを奇数行と偶数行に分けて学習用と予測のテスト用にする。
>even.n<-2*(1:75)-1
>train.data<-iris[even.n,]
>test.data<-iris[-even.n,-5]
>iris.lab<-factor(c(rep("S",25),rep("C",25),rep("V",25)))
>train.data[,5]<-iris.lab
>iris.rp2<-rpart(Species~.,train.data)
>plotcp(iris.rp2)
図4 irisの奇数行のprintcpプロット
上記の結果から、さらに剪定する必要がないと判断し、作成したモデルを用いて予測・判別を行う。
>iris.rp2r<-predict(Iris.rp2,test.data,type="class")
> table(Iris.lab,Iris.rp2r)
(4) 樹木の剪定とコントロール
関数rpartが返す樹木は関数rpart.controlによりコントロールされている。その引数の説明は、help(rpart.control)で確認できる。その主な引数の設定を説明するため、rpart.controlのデフォルトの設定を次に示す。
rpart.control(minsplit=5, minbucket= round(minsplit/3), cp=0.01, maxcompete=4, maxsurrogate=5, usesurrogate=2, xval=10,
surrogatestyle=0, maxdepth=30, ...)
引数minsplitは分岐をする際に、ノードに含まれる最少の個体の数をコントロールする。デフォルトは５に設定されている。値が小さいほど大きな木が作成される。
引数minbucketは葉における最少の個体数をコントロールする。デフォルトはminsplitの三分の一になっている。よって、minsplitとminbucketの中でひとつのみが指定されているときには、自動的に3*minbucket=minsplitによって計算される。
引数cp(complexity parameter)はすでに説明したように、複雑さのパラメータである。明らかに価値のない木が成長し過ぎることを制御するためのパラメータである。デフォルトは0.01になっている。値が小さいほど樹木が茂る。
引数maxdepthは、樹木の最大の深さをコントロールする。根の深さは０とカウントされる。デフォルトは30になっているが、32ビットマシンではこれ以上設定することは意味がないであろう。
回帰木で用いる分岐基準は、実測値
とこのセルの平均値
との差の2乗の和である。
ここでも実データを用いて回帰木モデルの作成などを説明する。関数treeの説明と同じくデータcarsを用いることにする。
データcarsは50個体、2変数のデータフレームである。変数は、車の速度(speed)とブレーキを掛けた後車が停車するまでの距離(dist)の計測値である。ここでは速度(speed)を説明変数、距離(dist)を被説明変数とする。
(1) 回帰木の作成
関数rpartによる回帰木の作成は、関数treeの書き式と同じである。
>(cars.rp<-rpart(dist~speed,data=cars))
＜紙面の都合で返される結果は省略する＞
返される結果から、回帰木は６の葉を持っていることが分かる。
(2) 回帰木の剪定
まず関数plotcpを用いて成長させた回帰木について剪定の必要性について考察する。
>
plotcp(cars.rp)
図5 回帰木のplotcpプロット
図5は成長させた回帰木を4つの葉に剪定する必要があることを示唆している。
次のように関数pruneを用いることにより4つの葉に剪定される。
> cars.rp1<-prune(cars.rp,cp=0.044)
>plot(cars.rp1,uniform=T,margin=0.05)
>text(cars.rp2,use.n=T)
＜紙面の都合で返される結果は省略する＞
この回帰木は関数treeを用いて作成した結果と同じである。
通常、作成したモデルを用いて予測を行うときには、モデル作成に使用していない説明変数のデータを用いて予測値を求める。関数rpartでは、predictを用いる。その書き式は分類木の場合と同じであるのでここでは省略する。残差は関数residualsで返される。
4. 多変量回帰木
多変量回帰(multivariate regression ) は、被説明変数が複数である回帰である。
ここでは、パッケージmvpartに同梱されているデータspiderを用いて多変量回帰木を説明する。
データspiderは12種類の蜘蛛狩りの分布と砂丘周辺の環境の特徴に関する28行18列のデータフレームである。
データの初めの12列は12種類の蜘蛛狩りに関する数値データの分布で、それに続く6列は砂丘周辺の環境の特徴に関するデータである。
ここでは砂丘周辺の環境の特徴に関する6列のデータを説明変数とし、12種類の蜘蛛狩りの多変量回帰木を作成する例を示す。関数rpartを用いて多変量回帰木を作成することもできるが、関数mvpartを用いた方がもっと便利である。
関数rpart、mvpartを用いて多変量回帰木を成長させる際には、用いた多変量の応答変数(被説明変数)をマトリックス形式にする必要がある。関数mvpartにはrpartより豊富な引数が用意されている。
その全てを紹介する紙面がないので、主な引数の使用方法を例示する。
関数mvpartでは、直接交差確認法の情報に基づいて剪定を行った樹木のグラフを返すように設定されている。もちろん、その設定は後術の引数を用いて替えることができる。
次にデータspiderを用いた多変量回帰木を作成する。
>data(spider)
water sand moss reft
twigs herbs
1 9 0 1 1 9 5
>fit.mv<-mvpart(as.matrix(spider[,1:12])~water+twigs+reft+herbs+moss+sand,spider)
図6 spiderの多変量回帰木
xv = c("1se",
"min", "pick", "none")
引数xv = "1se"はMin+1SEの方法による木を返し、xv="min"は最も当てはめの良い木を返す。xv="pick"はrpartにおけるprintcpと同じのグラフと木を返し、xv="none"は交差確認を行わない木を返す。デフォルトではxv
= "1se"になっている。
回帰木のルールを返すには関数printを用いる。
引数pcaは、樹木モデルと被説明変数の主成分分析の結果をリンクをコントロールする。デフォルトはpca = FALSになっている。次にpca = TRUEにした結果を示す。木のグラフが返されたら、グラフの画面でマウスを右クリックすると被説明変数の主成分分析のグラフが返される。
>mvpart(as.matrix(spider[,1:12])~water+twigs+reft+herbs+moss+sand,spider,pca=T)
(a)
(b)
図7 データspaiderの多変量回帰木と主成分の散布図とのリンク
主成分分析のグラフは、関数rpart.pacを用いて作成することもできる。また、パッケージmvpartの中には9種類の距離を求める関数gdist、古典的多次元尺度法と回帰木とのリンク関数cmds.dissなどが用意され、これらを用いた多変量回帰分析を行うこともできる。