[連載] フリーソフトによるデータ解析・マイニング 第19

Rと樹木モデル()

 


1.樹木モデル関数rpartmvpart

関数rpart(recursive partitioning and regression trees)1997Terry TherneauBeth Atkinsonにより実装され、関数mvpartrpartのモデルの当てはめとプロットのための包括関数で、オーストラリア海洋研究所のGlenn De’athより実装された。

関数rpartの書き式を次に示す。より多くの引数やデフォルトの設定などはhelp(rpart)で確認することができる。

rpart(formula, method,・・・)

formulaの書き式は、treeと同じである。引数methodは、表1の中からひとつを選択することができる。引数methodが指定されていない場合は、formulaに指定されている応答変数(被説明変数)のデータ型から判断する。

 

表1 引数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)

 

 この結果の樹木のグラフは、次のように関数plottextを用いて作成することができる。

 

> 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は図と4辺との間の余白(マージン)を調整する。引数marginの値が大きいほど図が小さくなり、余白が大きい。

関数textにも複数の引数が用意されている。上記のコマンドで用いた引数use.nは、各ノードに含まれる個体の数を表示するか、しないかを指定する。デフォルトにはuse.n=FALSEになっている。引数がuse.n=TRUEである場合は、ノードに含まれる数を表示する。

 

(2) 樹木の剪定

関数rpartは樹木を成長させると同時に、交差確認法の結果も計算する。関数printcpでその結果を返す。デフォルトのn重交差確認のnxval=10に設定されている。

関数printcpは、樹木の剪定のための複雑さのパラメータ(cp)を返す。この結果は、用いたデータセットをランダムに分割して交差確認を行った結果であるので、同じデータを用いて同じの条件のもとでrpartを繰り返しても、まったく同じの結果になるとは限らない。

 

>printcp(iris.cp)

 

通常、樹木のサイズはxerrorの最小値からその標準偏差1倍の範囲内の最大のxerror値を選ぶ。これをMin+1SE方法と呼ぶことにする。

上記に返された最終行のxerroxstdはそれぞれ0.070.0258であり、Min+1SE0.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になっている。

関数plotcpは、樹木の選定に必要となる交差確認に関する情報のプロットを返す。

 

>plotcp(iris.rp)

2 iris.rpplotcpのプロット

 

グラフの下部の横軸がcpで、上部の横軸が樹木の葉の数で示すサイズである。図の中の水平直線Min+1SEは、xerrorの最小値に1倍の標準偏差を足した値である。その直線の下方に最も近い点をオレンジ色で、xerrorの最小値の点を赤い色で示している。オレンジ色の点が樹木を剪定する目安である。

2ではcp=0.094、樹木のサイズが3に対応するところにオレンジ色の点が位置している。

この結果とplotcpでの剪定の目安cp=0.02が一致していないことに疑問があるに違いない。これは、printcpでは分岐の数とcpの対応関係で、plotcpは樹木のサイズとcpの対応関係であるからである。樹木の剪定にはどの結果を用いても良い。

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)

        iris.rp2r

iris.lab C  S  V

       C 24  0  1

       S  0 25  0

       V  3  0 22

 

(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は分岐をする際に、ノードに含まれる最少の個体の数をコントロールする。デフォルトは5に設定されている。値が小さいほど大きな木が作成される。

引数minbucketは葉における最少の個体数をコントロールする。デフォルトはminsplitの三分の一になっている。よって、minsplitminbucketの中でひとつのみが指定されているときには、自動的に3*minbucket=minsplitによって計算される。

引数cp(complexity parameter)はすでに説明したように、複雑さのパラメータである。明らかに価値のない木が成長し過ぎることを制御するためのパラメータである。デフォルトは0.01になっている。値が小さいほど樹木が茂る。

maxcompetemaxsurrogateusesurrogatesurrogatestyleは分岐の候補、欠損値に関して設定を行う引数である。これらの結果は、summaryから確認できる。

引数xvaln重交差確認のnを設定する。デフォルトは10になっている。通常310の値を指定する。

引数maxdepthは、樹木の最大の深さをコントロールする。根の深さは0とカウントされる。デフォルトは30になっているが、32ビットマシンではこれ以上設定することは意味がないであろう。

 関数rpart.controlの引数は、関数rpartで直接引数として用いることができる。

 

3.回帰木

回帰木で用いる分岐基準は、実測値とこのセルの平均値との差の2乗の和である。

ここでも実データを用いて回帰木モデルの作成などを説明する。関数treeの説明と同じくデータcarsを用いることにする。

データcars50個体、2変数のデータフレームである。変数は、車の速度(speed)とブレーキを掛けた後車が停車するまでの距離(dist)の計測値である。ここでは速度(speed)を説明変数、距離(dist)を被説明変数とする。

 

(1) 回帰木の作成

 関数rpartによる回帰木の作成は、関数treeの書き式と同じである。

 

>(cars.rp<-rpart(dist~speed,data=cars))

<紙面の都合で返される結果は省略する>

 

返される結果から、回帰木は6の葉を持っていることが分かる。

 

(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を用いて多変量回帰木を説明する。

データspider12種類の蜘蛛狩りの分布と砂丘周辺の環境の特徴に関する2818列のデータフレームである。

データの初めの12列は12種類の蜘蛛狩りに関する数値データの分布で、それに続く6列は砂丘周辺の環境の特徴に関するデータである。

ここでは砂丘周辺の環境の特徴に関する6列のデータを説明変数とし、12種類の蜘蛛狩りの多変量回帰木を作成する例を示す。関数rpartを用いて多変量回帰木を作成することもできるが、関数mvpartを用いた方がもっと便利である。

関数rpartmvpartを用いて多変量回帰木を成長させる際には、用いた多変量の応答変数(被説明変数)をマトリックス形式にする必要がある。関数mvpartにはrpartより豊富な引数が用意されている。

その全てを紹介する紙面がないので、主な引数の使用方法を例示する。

 関数mvpartでは、直接交差確認法の情報に基づいて剪定を行った樹木のグラフを返すように設定されている。もちろん、その設定は後術の引数を用いて替えることができる。

 次にデータspiderを用いた多変量回帰木を作成する。

 

>data(spider)

>spider[1,13:18]

  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などが用意され、これらを用いた多変量回帰分析を行うこともできる。