統計解析環境Rによる樹木構造接近法の実践

統計解析環境Rによる
樹木構造接近法の実践
下川敏雄
山梨大学大学院医学工学総合研究部
分類樹木の適用例
４種類の酵素が，４種類の肝臓疾患患者に対して取られている
(Albert & Harris, 1987)．
グループ変数: 肝臓病の種類
・ D1: 急性ウィルス性肝炎(57例)
・ D3: 活動性慢性肝炎 (40例)
・ D2: 持続性慢性肝炎 (44例)
・ D4: 破死後性肝硬変(77例)
説明変数：4種類の酵素
・アスパラギン酸アミノ基転移酵素(AST)
・グルタミン酸デヒドロゲナーゼ（GLDH）
・アラニン酸アミノ基転移酵素 (ALT)
・オルニチン・カルバモイル・トランスフェラ
ーゼ (OCT)
いずれの疾患もウィルス性である．アルコール性と異なり，処置を行わなければ，
肝臓がんを併発する惧れがある．
活動性慢性肝炎
持続性慢性肝炎
重症度が高い
破死後性肝硬変
急性ウィルス性肝炎
AST, ALT, GLDH, OCT:健康診断などの血液検査で一般的にとられている検査値
データの入力と内容
データの入力(本データはパッケージpmlrのなかのデータ集合enzymesにある)
> install.packages("pmlr")
> data(enzymes, package="pmlr")
> summary(enzymes)
Patient
Min.
: 1.00
1st Qu.:14.00
Median :28.00
Mean
:29.66
3rd Qu.:41.75
Max.
:77.00
Group
1:57
2:44
3:40
4:77
Group
1:57
2:44
3:40
4:77
AST
Min.
: 20.00
1st Qu.: 54.25
Median : 94.00
Mean
: 159.53
3rd Qu.: 165.75
Max.
:2298.00
Groupの番号
# パッケージpmlrのインストール
# データの入力
# データの要約の出力
ALT
Min.
: 18.0
1st Qu.: 52.0
Median : 88.0
Mean
: 244.9
3rd Qu.: 329.2
Max.
:1595.0
GLDH
Min.
: 2.00
1st Qu.: 7.00
Median : 11.00
Mean
: 14.65
3rd Qu.: 18.00
Max.
:160.00
OCT
Min.
: 18.0
1st Qu.: 116.2
Median : 228.5
Mean
: 399.6
3rd Qu.: 424.5
Max.
:5067.0
1：急性ウィルス性肝炎(acute viral hepatitis)
2：持続性慢性肝炎(persistent chronic hepatitis)
3：活動性慢性肝炎(aggressive chronic hepatitis)
4：破死後性肝硬変(post-necrotic cirrhosis)
rpartパッケージとその留意点
> library(rpart)
> mdl.rpart <- rpart(Group~AST+ALT+GLDH+OCT, data=enzymes)
> plot(mdl.rpart, margin=0.1)
> text(mdl.rpart)
ALT>=227.5
|
注意!!
関数rpart()のplotではそのまま描くと
テキストが切れてしまう．そのため，
marginを入れて，テキストが切れない
ようにしなければならない．
GLDH< 24.5
1
AST< 40.5
3
ALT>=85.5
2
GLDH< 14.5
4
2
3
rpartの引数とデフォルトの関係
関数rpart()では，下記の形式で成長する樹木のパラメータを設定できる．
rpart(formula, control=rpart.control(optons) )
ここで，optionに記載できるパラメータとその意味は下表のとおり．
形式
デフォルト
説明
minsplit
minsplit = 20
分岐ふしにおける最小個体数
minbucket
minbucket=round(minsplit/3)
終結ふしにおける最小個体数
ALT>=227.5
|
maxcompete maxcompete
=4
分岐候補として競合した代替変数の保持する数
xval
xval = 10
交差確認回数
maxdepth
maxdepth = 30
樹木の深さの最大値
cp
cp = 0.01
樹木の複雑度パラメータの最小値
GLDH< 24.5
1
3
TM.1
TM.2
AST< 40.5
ALT>=85.5
2
TM.3
GLDH< 14.5
4
2
3
TM.4
TM.5
TM.6
さきほどの事例では，TM.1(n=59), TM.2(n=7),
TM.3(n=28), TM.4(n=17), TM.5(n=27),
TM.6(n=80)だった．すなわち，今回の樹木
は，TM.2がminbucketのパラメータに引っか
かっただけである．
rpart()のオプションを変更した場合
> mdl.rpart2 <- rpart(Group~AST+ALT+GLDH+OCT, data=enzymes,
+
control=rpart.control(minbucket=2)) # 終結ふしにおける最小個体数を2に変更
> plot(mdl.rpart2, margin=0.1)
> text(mdl.rpart2, cex=0.9) # cexは文字数を変更するオプション(デフォルト1)
ALT>=227.5
注意!!
|
関数rpart()は分岐過程が終了しただけ
であり，交差確認法による最適樹木
の選定までを行っているわけではな
い．
GLDH< 24.5
1
AST< 40.5
3
ALT>=85.5
2
増加
GLDH< 14.5
AST< 135.5
GLDH< 6.5
3
4
1
2
4
交差確認法による最適樹木の選定
理論編でも述べたが，CART法では交差確認法を用いて，樹木系列のなかから最
適なサイズの樹木を選ぶ．このとき，全データでの樹木系列と交差確認法で得
られた樹木系列は，複雑度パラメータで結び付けられる．
> library(rpart)
> mdl.rpart3 <- rpart(Group~AST+ALT+GLDH+OCT, data=enzymes,
+
control=rpart.control(cp=0)) # 可能な限り大きな樹木を構成
> j0 <- which.min(mdl.rpart$cptable[,"xerror"]) # 誤分類率の交差確認推定値が最小
> rtree_0se <- prune(mdl.rpart3,mdl.rpart$cptable[j0,"CP"])
> plot(rtree_0se, margin=0.1)
> text(rtree_0se)
ALT>=227.5
|
GLDH< 24.5
1
結果的には，rpart()のデフォルトと同様の結果にな
ったが，関数rpart()をCARTアルゴリズムに合致させ
るためには，上記の記述に沿うことが必要．
AST< 40.5
3
ALT>=85.5
2
GLDH< 14.5
4
2
3
交差確認法による損失関数の交差確認推定値
> plotcp(mdl.rpart3)
size of tree
2
3
4
5
6
10
終結ふしの個数
0.6
0.8
ここが，最小交差確認推定値に
なる．
0.4
X-val Relative Error
1.0
1
Inf
0.26
0.11
0.063
cp
0.045
0.011
0
複雑度パラメータ
グラフィカル出力の工夫：partykitの利用
関数rpart()の樹木の出力は貧相である．そのため，パッケージparty
のフォーマットに変換するほうが綺麗なグラフ表現が可能である．
ALT>=227.5
|
GLDH< 24.5
1
> install.packages("party")
> library(party)
> party_0se <- as.party(rtree_0se)
> plot(party_0se)
AST< 40.5
3
ALT>=85.5
2
GLDH< 14.5
4
2
# partyのインストール
# partyの読み込み
# partyフォーマットに変換
# party樹木の出力
3
1
ALT
>= 227.5
< 227.5
2
5
GLDH
AST
< 40.5
>= 40.5
7
ALT
< 24.5
>= 24.5
>= 85.5
8
< 85.5
終結ふしのグラフィク
スが見やすくなる(左ほ
ど1，右ほど4)ように娘
ふしの置き方が工夫さ
れている．
GLDH
< 14.5
Node 3 (n = 59)
1
0.8
0.6
0.4
0.2
0
Node 4 (n = 7)
1
0.8
0.6
0.4
0.2
0
1 2 3 4
Node 6 (n = 28)
1
0.8
0.6
0.4
0.2
0
1 2 3 4
Node 9 (n = 17)
1
0.8
0.6
0.4
0.2
0
1 2 3 4
1 2 3 4
>= 14.5
Node 10 (n = 27)
1
0.8
0.6
0.4
0.2
0
1 2 3 4
Node 11 (n = 80)
1
0.8
0.6
0.4
0.2
0
1 2 3 4
交差確認法による最適樹木サイズの選定：1SEルールの利用
CART法では，最小損失関数の不安定性を考えて．保守的な樹木を構築することが考
えられる．この場合，最小損失+損失関数以下となる最小サイズの樹木を選択する．
このことを1SEルールと呼ぶ．
1SEルールによる最適樹木サイズの選定(one_se.R)
one.se.rule <- function(tree){
j0 <- which.min(tree$cptable[,"xerror"])
Rcv.min <- tree$cptable[j0,"xerror"]
one.se <- tree$cptable[j0,"xstd"]
j1 <- 1
while(tree$cptable[j1,"xerror"] > Rcv.min + one.se)
{j1 <- j1+1}
return(tree$cptable[j1,])
}
oe_se.Rが"C:/lecture/"にあると想定
> plotcp(mdl.rpart3)
> out_1se <- one_se_rule(mdl.rpart3 )
> rtree_1se <- prune(mdl.rpart3, cp=out_1se["CP"])
> party_1se <- as.party(rtree_1se)
> plot(party_1se)
1SEルールによって選定された部分樹木
1
ALT
>= 227.5
< 227.5
3
AST
< 40.5
>= 40.5
5
ALT
>= 85.5
< 85.5
6
GLDH
< 14.5
Node 2 (n = 66)
1
0.8
0.6
0.4
0.2
0
Node 4 (n = 28)
1
0.8
0.6
0.4
0.2
0
1
2 3
4
Node 7 (n = 17)
1
0.8
0.6
0.4
0.2
0
1
2 3
4
>= 14.5
Node 8 (n = 27)
1
0.8
0.6
0.4
0.2
0
1 2
3
4
Node 9 (n = 80)
1
0.8
0.6
0.4
0.2
0
1
2
3 4
1
2 3
4
急性ウィルス性肝炎
(ALT ≧ 227.5)
持続性慢性肝炎：
(ALT < 228.5)∩(AST<40.5)，
(85.5≦ALT < 228.5)∩(AST≧40.5)∩(GLDH<14.5)
活動性慢性肝炎
(85.5≦ALT < 228.5)∩(AST≧40.5)∩(GLDH≧14.5)
破死後性肝硬変
つまり，ALTが低く，ASTが高い状況が「破死後
(ALT < 85.5)∩(AST≧40.5)
性肝硬変」の惧れ
変数重要度のプロット
rpartのバージョン4.03以降より，オブジェクトのなかに変数重要度が追加されるよ
うになった．
60
40
20
AST, ALTの変数重要度が同程度に高
かった．
0
variable importance
80
100
> vimp <- mdl.rpart3$variable.importance
> rimp <- vimp/max(vimp)*100
> barplot(sort(rimp, decreasing=TRUE), ylab="variable importance", col="blue")
AST
ALT
OCT
GLDH
分類樹木における基準の変更
自動者価格データ(パッケージrpartに入っているファイル(Chambers & Hastie, 1992))
- 応答(品質)：非常に悪い，悪い，平均的，良い，非常に良い(本来は順序尺度)
data(cu.summary)
## Gini係数による分類樹木(default)
fit1 <- rpart(Reliability ~ Price + Country + Mileage + Type,
data = cu.summary, parms = list(split = 'gini'))
## エントロピーによる樹木
fit2 <- rpart(Reliability ~ Price + Country + Mileage + Type,
data = cu.summary, parms = list(split = 'information'))
## 樹木の描写
par(mfrow = c(1,2), mar = rep(0.1,4))
plot(fit1, margin = 0.1); text(fit1, use.n = TRUE, cex = 0.7)
plot(fit2, margin = 0.1); text(fit2, use.n = TRUE, cex = 0.7)
出力結果
エントロピーによる
分類樹木
Gini係数による
分類樹木
Country=dghij
Country=dghij
|
|
Type=e
Much better
0/0/3/3/21
Type=bcef
Much better
0/0/3/3/21
Type=e
Type=bcf
Country=gj
Much w orse
7/0/2/0/0
w orse
4/6/1/3/0
Country=gj
average
7/4/16/0/0
w orse
4/6/1/3/0
average
0/2/4/2/0
Much w orse
7/0/2/0/0
average
0/2/4/2/0
average
7/4/16/0/0
不純性に対して，エントロピー＞ Gini係数＞誤分類率の順でふし内不均一性が高く
なることが反映している(最適な指標に関する報告は少ない)．
回帰樹木の適用例
インドの小児の栄養失調データ (Fenske et al, 2009)．
データの入力(本データはパッケージexpectregのなかのデータ集合indiaにある)
> install.packages("expectreg")
> data(india, package="expectreg")
> summary(india)
stunting
Min.
:-599.0
1st Qu.:-287.0
Median :-176.0
Mean
:-175.4
3rd Qu.: -65.0
Max.
: 564.0
distH
Min.
: 5.0
1st Qu.:118.0
Median :232.0
Mean
:235.4
3rd Qu.:351.0
Max.
:482.0
cbmi
Min.
:10.03
1st Qu.:14.23
Median :15.36
Mean
:15.52
3rd Qu.:16.60
Max.
:25.95
# パッケージpmlrのインストール
# データの入力
# データの要約の出力
cage
Min.
: 0.00
1st Qu.: 8.00
Median :17.00
Mean
:17.23
3rd Qu.:26.00
Max.
:35.00
mbmi
Min.
:13.14
1st Qu.:17.85
Median :19.36
Mean
:19.81
3rd Qu.:21.21
Max.
:39.81
mage
Min.
:13.00
1st Qu.:21.00
Median :24.00
Mean
:24.41
3rd Qu.:27.00
Max.
:49.00
stunting：栄養不良スコア
cbmi：小児のBMI
cage：小児の月年齢
mtmi：母親のBMI
mage：母親の年齢
distH：インドの地区のコード番号
rpartの一連の過程
> library(rpart)
> mdl.regress <- rpart(stunting~cbmi+cage+mbmi+mage, data=india,
+
control=rpart.control(cp=0))
# 可能な限り大きな樹木を構成
> out_1se <- one_se_rule(mdl.regress)
# 誤分類率の交差確認推定値が最小
> regree_1se <- prune(mdl.regress, cp=out_1se["CP"])
#1SEルールによるモデル選択
> partreg_1se <- as.party(regree_1se)
# partyオブジェクトに変換
> plot(partreg_1se)
# partyの樹木で表示
1
cage
>= 11.5
母親のBMI
小児の月年齢
< 11.5
2
7
mbmi
cage
< 20.025
小児の月年齢
>= 20.025
4
cbmi
子供のBMI
>= 4.5
< 4.5
>= 18.01768< 18.01768
Node 3 (n = 1590)
Node 5 (n = 130)
600
600
400
400
200
200
0
0
-200
-200
-400
-400
-600
-600
Node 6 (n = 945)
600
400
200
0
-200
-400
-600
Node 8 (n = 831)
600
400
200
0
-200
-400
-600
Node 9 (n = 504)
600
400
200
0
-200
-400
-600
変数重要度
100
> vimp <- mdl.regress$variable.importance
> rimp <- vimp/max(vimp)*100
> barplot(sort(rimp, decreasing=TRUE), ylab="variable importance", col="blue")
子供のつき年齢の影響が最も強く，
80
次いで，子供のBIM，母親のBMIの
20
40
60
は非常に低かった．
0
variable importance
影響が強かった．母親の年齢の影響
cage
cbmi
mbmi
mage
応答とCART樹木形式との関係
関数rpart()では，回帰，分類，計数，生存時間(指数分布)といった応答に対応して
いる．オプションで，metod="手法名"で指定できるが，変数の形式で省くことが
できる．
意味
method="○○"
短縮できる応答の型
"anova"
回帰樹木
応答が数値型(numeric)
"class"
分類樹木
応答が因子型(factor)
"poisson"
Poisson回帰樹木
応答が2列(成功数，個体数)
"exp"
生存時間(指数回帰)樹木
応答がSurvによる生存時間型
partyオブジェクトへの変換後のグラフィクスはそれぞれに対応している．
意味
method="○○"
終結ふしの形式
"anova"
回帰樹木
ボックスプロット
"class"
分類樹木
帯グラフ(2値)，バーチャート(多値)
"poisson"
Poisson回帰樹木
バーチャート
"exp"
生存時間(指数回帰)樹木 Kaplan-Meierプロット
パッケージmvpartのなかのrpart関数
■ mvpartパッケージにも関数rpart()が存在する．ただし，デフォルトが異なるもの
があるので注意．
形式
rpartパッケージのデフォルト
mvpartパッケージのデフォルト
minsplit
minsplit = 20
minsplit = 5
minbucket
minbucket=round(minsplit/3)
minbucket=round(minsplit/3)
maxcompete
maxcompete = 4
maxcompete = 4
xval
xval = 10
xval = 10
maxdepth
maxdepth = 30
maxdepth = 30
cp
cp = 0.01
cp = 0.01
■ methodが追加される
意味
method="○○"
短縮できる応答の型
"mrt"
多変量回帰樹木
応答が行列型 (変数×個体)
"dist"
分類樹木
応答が行列型 (個体×個体)
mrtの個体間距離 dissim="euc"：Euclid距離(L2ノルム距離)
dissim="man"：Manhattan距離(L1ノルム距離)
関数mvpartの概要
■ mvpartパッケージの関数mvpartはrpartのラッパーであり，オプションが若干異なる．
形式
説明
xv
交差確認法による最適モデル推定
("1se"：1SEルール，"min"最小値，"none"：なし)
plot.add, text.add
コマンド実行時に樹木/説明テキストの表示(default: TRUE)
xval
交差確認回数
xvmult
多重交差確認回数
多重交差確認法
 損失関数に対する交差確認推定値をL回繰り返し，その平均値をとるのが多重交
差確認法である．
 標本サイズが小規模の場合に，交差確認推定値が乱数の影響を受けやすい．多
重交差確認法では，このことを緩和することができる．
計数応答に対する樹木(Poisson回帰樹木)
ガラパコス諸島の植物生態データ(Johnson & Raven, 1973)
これは，ガラパコス諸島の各島の生物多様性を評価するために，それぞれの島に
繁殖する植物の種数を調査したデータである．
> install.packages("faraway")
> data(gala, package="faraway")
> summary(gala)
Species
Min.
: 2.00
1st Qu.: 13.00
Median : 42.00
Mean
: 85.23
3rd Qu.: 96.00
Max.
:444.00
Scruz
Min.
: 0.00
1st Qu.: 11.03
Median : 46.65
Mean
: 56.98
3rd Qu.: 81.08
Max.
:290.20
Endemics
Min.
: 0.00
1st Qu.: 7.25
Median :18.00
Mean
:26.10
3rd Qu.:32.25
Max.
:95.00
Adjacent
Min.
:
0.03
1st Qu.:
0.52
Median :
2.59
Mean
: 261.10
3rd Qu.: 59.24
Max.
:4669.32
Area
Min.
:
0.010
1st Qu.:
0.258
Median :
2.590
Mean
: 261.709
3rd Qu.: 59.237
Max.
:4669.320
Elevation
Min.
: 25.00
1st Qu.: 97.75
Median : 192.00
Mean
: 368.03
3rd Qu.: 435.25
Max.
:1707.00
Nearest
Min.
: 0.20
1st Qu.: 0.80
Median : 3.05
Mean
:10.06
3rd Qu.:10.03
Max.
:47.40
Species：植物の種数
Endemics：固有種の数
Area：島の面積(km2)
Elevation：島の最大対角(m)
Nearest：一番近い島までの距離(km)
Scruz：Santa Cruz 島からの距離(km)
Adjacent：付近の島の面積(Adjacent)
樹木の構成(1)：損失関数の交差確認推定値の省察
> ptree0 <- rpart(Species~., data=gala, method="poisson"
+
,control=rpart.control(cp=0,minsplit=5))
> plotcp(ptree0)
size of tree
2
3
4
5
6
7
8
9
10
11
1.0
0.5
殆ど変化がない．樹木の安定性と
簡単さを考えれれば，1SEルール
がふさわしい．
0.0
X-val Relative Error
1.5
1
Inf
0.35 0.06 0.02
0.0094
cp
0.0028
0.0013
0
樹木の構成(2)：1SEルールによる樹木表現
> out_1se <- one_se_rule(ptree0)
> ptree1<- prune(ptree0,cp=out_1se["CP"])
> plot(as.party(ptree1))
1
固有種の数
Endemics
< 51
>= 51
2
Endemics
固有種の数
固有種の数の影響のみが存在している．
< 14
>= 14
Node 3 (n = 13)
Node 4 (n = 12)
Node 5 (n = 5)
400
400
400
300
300
300
200
200
200
100
100
100
0
0
0
順序カテゴリカル応答の場合
順序カテゴリカル応答に対する樹木は，パッケージrpartScoreに存在する．
> install.packages("rpartScore")
> library(rpartScore)
rpartScore(formula, Options, control=rpart.control(optons) )
関数rpartScore()では，2種類のふし内不均一性測度が提案されている．いずれも
Piccarreta(2009) が提案した順序付き一般化Gini係数による方法の拡張版(Galimberti et
al., 2012)である．分岐基準は，J個の順序付きスコア{s1 < s2 <･･･< sJ}に対して
J
分岐過程に対す
るオプション
誤分類による損失
rGG (t )  C ( j | j ') Pr( j | t ) Pr( j ' | t )
j 1 j   j
C ( j | j ') | s j  s j |
Options
split="abs"
C ( j | j ')  ( s j  s j ) 2 split="quad"
なお，刈り込み過程における複雑度の定義は下記のとおり：
刈り込み過程に対
するオプション
R (T )   R (t )   | T |
tT
Options
split="mr"
誤分類率
応答のスコアと推定 prune="mc"
スコアの差の絶対値
マンモグラフィ検査に対するアンケートデー (Hosmer & Lemeshow, 2000)
・応答：マンモグラフィ検査の受診経験
1：受診したことがない，2：一年以内に受診，3：受診中
説明変数
質問「あなたは乳がんの兆候がなければマンモグラフィ検査を受ける必
要がない」：SYMPT
0：強く思う，1：思う，2：思わない，3：強く思わない
質問「母親あるいは姉妹に乳がん経験があるか否か」：HIST
0：いいえ，1：はい
質問「自らの触診によるしこりの発見法を指導してもらった経験がある
か」：BSE
0：いいえ，1：はい
質問「マンモグラフィ検査により乳がんを発見できると思うか」：DECT
0：強く思う，1：思う，2：思わない，3：強く思わない
マンモグラフィ検査の利点に対する質問に対する点数(高いほど正解が多
い：PB
順序カテゴリカル応答に対するCART法の実行(1/5)
分岐基準のコストがスコア差の絶対値の場合(複雑度はコストに基づく)
> data(mammoexp, package="TH.data") # TH.dataに入っているデータmammoexpを入力
> ord.rpart.abs <- rpartScore(as.numeric(ME)~., data=mammoexp, split="abs", prune="mc",
+
control=rpart.control(cp=0,minsplit=5))
> plotcp(ord.rpart.abs)
rpartScore関数では，応答が順序付き因子型でなくnumeric型でないといけない(それぞ
れのカテゴリに対するスコアに対応させるため)
size of tree
4
5
6
7
11
13
15
21
0.8
0.9
1.0
1.1
2
最小損失による選択
1SEによる選択
0.7
X-val Relative Error
1
Inf
0.063
0.034
0.019
0.014
0.0084
cp
0.0028
0
順序カテゴリカル応答に対するCART法の実行(2/5)
> out_1se <- one_se_rule(ord.rpart.abs)
# 誤分類率の交差確認推定値が最小
> ord.rpart.1se <- prune(ord.rpart.abs, cp=out_1se["CP"]) #1SEルールによるモデル選択
> plot(ord.rpart.1se, margin=0.05)
> text(ord.rpart.1se, cex=0.9)
SYMPT< 2.5
|
HIST< 1.5
PB< 8.5
1
2
マンモグラフィー検査
1：受診したことがない
2：一年以内に受診
3：受診中
注意!!
BSE< 1.5
1
1
2
rpartScoreでは応答がnumeric
型になるので，partyの樹木
グラフを描写することがで
きない． 1SE を探す関数は
rpartと同様に利用できる．
27
順序カテゴリカル応答に対するCART法の実行(3/5)
分岐基準のコストがスコア差の2乗値の場合(複雑度はコストに基づく)
> ord.rpart.quad <- rpartScore(as.numeric(ME)~., data=mammoexp, split="quad",prune="mc",
+
control=rpart.control(cp=0,minsplit=5))
> plotcp(ord.rpart.quad)
size of tree
4
5
7
8
13
19
25
0.0019
0
0.8
0.9
1.0
1.1
2
1SEによる選択
最小損失による選択
0.7
X-val Relative Error
1
Inf
0.063
0.034
0.015
0.0089
cp
0.0056
0.0032
順序カテゴリカル応答に対するCART法の実行(4/5)
> out_1se <- one_se_rule(ord.rpart.quad)
# 誤分類率の交差確認推定値が最小
> ord.rpart.1se <- prune(ord.rpart.quad, cp=out_1se["CP"]) #1SEルールによるモデル選択
> plot(ord.rpart.1se, margin=0.05)
> text(ord.rpart.1se, cex=0.9)
SYMPT< 2.5
|
HIST< 1.5
PB>=8.5
1
2
マンモグラフィー検査
1：受診したことがない
2：一年以内に受診
3：受診中
BSE< 1.5
1
1
2
いずれの場合にも，最終的に得られた樹木は変わらなかった．
順序カテゴリカル応答に対するCART法の実行(5/5)
> vimp.abs <- ord.rpart.abs$variable.importance
> vimp.quad <- ord.rpart.quad$variable.importance
> rimp.abs <-sort(vimp.abs/max(vimp.abs), decreasing=TRUE)*100
> rimp.quad <-sort(vimp.quad/max(vimp.quad), decreasing=TRUE)*100
> par(mfrow = c(1,2))
> barplot(rimp.abs, ylab="variable importance", col="blue", main="abs")
> barplot(rimp.quad, ylab="variable importance", col="blue", main="quad")
80
60
40
0
20
variable importance
60
40
20
0
variable importance
80
100
quad
100
abs
SYMPT
PB
DECT
BSE
HIST
SYMPT
PB
BSE
刈り込まれた末端付近の樹木による違いが反映されたようである．
HIST
DECT
ふし間分離度に基づく方法：条件付き推測樹木
条件付き推測樹木は，変数選択とふし点選択を個別に
行い，変数選択では，帰無仮説
親ふし
H 0p : F ( y )  F ( y | x p )
娘ふし
VS
娘ふし
を並べ替え検定により行い，最小のp値をもつ変数に
対して，ふし点選択を行う．任意のふし点を用いて，
2値化を行ったうえで，同様の検定を行い，最小のp
値をもつふし点で分岐する．
partyパッケージでは，上記の検定を条件付き推測パッケージcoinを用いて行う．
coinパッケージでは，検定に用いる変数(応答)の形式により，回帰樹木，分類樹
木，順序カテゴリカル応答に対する樹木，および生存時間樹木に適用可能であ
る．
検定統計量に基づく方法では，分岐過程のみで終了する．そして，停止条件に
はそれぞれのふしに対して行う条件付き推測に対して設定された有意水準αを用
いる．
関数ctreeの概要
ctree
回帰樹木
応答が数値型(numeric)
分類樹木
応答が因子型(factor)
順序カテゴリカル応答
応答が順序付き型(ordered)
に対する樹木
生存時間樹木
応答がSurvによる生存時間型
形式
説明 (緑色はデフォルト)
teststat
分岐における分岐点選択基準*1("quad", "max")
testtype
検定統計量の計算方式*2 ("Bonferroni", "MonteCarlo", "Univariate", "Teststatistic")
mincriterion
1-有意水準α (α=0.95)
minsplit
分岐ふしにおける最小個体数 (20)
minbucket
終結ふしにおける最小個体数 (7)
nresample
Monte-Carlo シミュレーション回数 (9999)
maxdepth
樹木の深さの最大値(0：指定無し)
*1：検定統計量の形式(詳しくは，Strasser & Weber (1999)参照)
*2：条件付き推測樹木では分岐変数選択，分岐点選択において多重比較を行う．
このときの調整形式を選択する．
条件付き推測樹木の例示(1)：分類樹木
> data(enzymes, package="pmlr") # 肝臓病データ
> class.party <- ctree(Group~AST+ALT+GLDH+OCT, data=enzymes)
> plot(class.party)
1
ALT
p < 0.001
224
224
2
ALT
p < 0.001
85
7
GLDH
p < 0.001
85
3
AST
p < 0.001
41
Node 4 (n = 29)
1
0.8
0.6
0.4
0.2
0
1 23 4
24
24
41
Node 5 (n = 78)
1
0.8
0.6
0.4
0.2
0
1 2 34
Node 6 (n = 45)
1
0.8
0.6
0.4
0.2
0
1234
Node 8 (n = 59)
1
0.8
0.6
0.4
0.2
0
1 23 4
Node 9 (n = 7)
1
0.8
0.6
0.4
0.2
0
1 2 34
多重比較調整を変更すると...
> class.party1 <- ctree(Group~AST+ALT+GLDH+OCT, data=enzymes,
+
controls = ctree_control(testtype="Bonferroni"))
> class.party2 <- ctree(Group~AST+ALT+GLDH+OCT, data=enzymes,
+
controls = ctree_control(testtype="MonteCarlo"))
> class.party3 <- ctree(Group~AST+ALT+GLDH+OCT, data=enzymes,
+
controls = ctree_control(testtype="Univariate"))
> class.party4 <- ctree(Group~AST+ALT+GLDH+OCT, data=enzymes,
+
controls = ctree_control(testtype="Teststatistic"))
> par(mfrow = c(2,2))
> plot(class.party1)
> plot(class.party2)
> plot(class.party3)
> plot(class.party4)
多重比較調整を変更すると...
1
ALT
p < 0.001
Bonferroni
224
MonteCarlo
224
2
ALT
p < 0.001
85
7
GLDH
p < 0.001
41
24
24
41
4
AST
p < 0.001
41
35
Node 5 (n = 78)
1
0.8
0.6
0.4
0.2
0
1
0.8
0.6
0.4
0.2
0
1 2 3 4
1 2 3 4
35
Node 5 (n = 22)
1
0.8
0.6
0.4
0.2
0
1234
1 2 3 4
224
85
85
Node 5 (n = 22)
1
0.8
0.6
0.4
0.2
0
35
224
41
24
1234
1234
Node 9 (n = 17)
1
0.8
0.6
0.4
0.2
0
1234
1234
1234
Node 12 (n = 59)
1
0.8
0.6
0.4
0.2
0
1234
Node 9 (n = 65)
1
0.8
0.6
0.4
0.2
0
1234
3
AST
p < 0.001
4
41
41 7
AST
OCT
p < 0.001
p < 0.001
63
63
85
9
AST
p < 0.001
48
24
Node 13 (n = 7)
1
0.8
0.6
0.4
0.2
0
1234
Node 10 (n = 45)
1
0.8
0.6
0.4
0.2
0
1234
2
ALT
p < 0.001
24
Node 12 (n = 59)
1
0.8
0.6
0.4
0.2
0
1234
Node 13 (n = 7)
1
0.8
0.6
0.4
0.2
0
1234
1234
1
ALT
p < 0.001
224
224
85 18
GLDH
p < 0.001
10
10 20
AST
p < 0.001
11
ALT
p < 0.001
12
61
61
AST
p < 0.001
84
84
24
GLDH
p < 0.001
25
20
20
GLDH
p < 0.001
23
GLDH
p < 0.001
24
24
48
10
Node 10 (n = 28)
1
0.8
0.6
0.4
0.2
0
24
63
Node 8 (n = 13)
1
0.8
0.6
0.4
0.2
0
35
35
Node 7 (n = 78)
1
0.8
0.6
0.4
0.2
0
11
GLDH
p = 0.015
85
Teststatistic
11
GLDH
p < 0.001
35
Node 6 (n = 7)
1
0.8
0.6
0.4
0.2
0
224
7
OCT
p = 0.002
63
Node 6 (n = 7)
1
0.8
0.6
0.4
0.2
0
1234
1 2 3 4
8
GLDH
p = 0.043
10
41
224
Node 9 (n = 7)
1
0.8
0.6
0.4
0.2
0
1 2 3 4
2
ALT
p < 0.001
3
AST
p < 0.001
41
Node 8 (n = 59)
1
0.8
0.6
0.4
0.2
0
1
ALT
p < 0.001
Univariate
4
AST
p = 0.028
Node 6 (n = 45)
1
0.8
0.6
0.4
0.2
0
2
ALT
p = 0.009
85
3
AST
p < 0.001
85
3
AST
p < 0.001
Node 4 (n = 29)
1
ALT
p < 0.001
13
ALT
p < 0.001
50
50
159
159
9
9
Node 5 (n = 22)
Node 6 (n = Node
7)
8 (n = 13)
Node 10 (n =Node
7) 14 (n =Node
15) 15 (n =Node
11) 16 (n =Node
11) 17 (n =Node
21) 19 (n =Node
17) 21 (n =Node
16) 22 (n =Node
12) 26 (n =Node
11) 27 (n = Node
39) 28 (n =Node
9) 29 (n = 7)
0.8
0.8
0.8
0.8
0.8
0.8
0.8
0.8
0.8
0.8
0.8
0.8
0.8
0.8
0.8
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1
1
1
1
1
1
1
1
1
1
1
1
1
1
Bonfferoni < MonteCarlo < Univariate < Teststatisticの順で樹木サイズが大きくなって
いる．つまり，Bonfferoniが最も保守的な樹木である．
1
条件付き推測樹木の例示(2)：回帰樹木
> data(india, package="expectreg")
> regress.party <- ctree(stunting~cbmi+cage+mbmi+mage, data=india,
+
controls = ctree_control(testtype="Bonferroni", mincriterion = 0.99))
> plot(regress.party)
注意!!
1
cage
p < 0.001
11
2
cage
p < 0.001
3
mbmi
p = 0.008
4
4
6
mbmi
p < 0.001
7
cage
p = 0.009
19
11
22.34
14
cbmi
p < 0.001
13
mbmi
p < 0.001
20.02
10
cbmi
p = 0.002
13.83813.838
17
cbmi
p < 0.001
18
mbmi
p < 0.001
18.14218.142
7
20.02
22.34
19
7
標本サイズが多い場合に，
条件付き推測樹木の適用(と
くに回帰)には注意が必要
24.94
19
cbmi
p < 0.001
18.009 18.009
24.94
15.854 15.854
21
cage
p = 0.009
27
27
Node 4 (n = 180)
Node 5 (n = 324)
Node 8 (n = 348)
Node 9 (n = 384)
Node 11 (n = 10)
Node 12 (n =Node
89) 15 (n = 1481)
Node 16 (n = 109)
Node 20 (n = 506)
Node 22 (n = 200)
Node 23 (n =Node
98) 24 (n = 141)
Node 25 (n = 130)
600
600
600
600
600
600
600
600
600
600
600
600
600
200
200
200
200
200
200
200
200
200
200
200
200
200
-200
-200
-200
-200
-200
-200
-200
-200
-200
-200
-200
-200
-200
-600
-600
-600
-600
-600
-600
-600
-600
-600
-600
-600
-600
-600
36
条件付き推測樹木の例示(3)：順序カテゴリカル応答に対する樹木
> data(mammoexp, package="TH.data") # TH.dataに入っているデータmammoexpを入力
> ord.ctree <- ctree(ME~., data=mammoexp)
> plot(ord.ctree)
1
SYMPT
p < 0.001
Agree
Agree
3
PB
p = 0.012
8
Node 2 (n = 113)
8
Node 4 (n = 208)
Node 5 (n = 91)
1
1
1
0.8
0.8
0.8
0.6
0.6
0.6
0.4
0.4
0.4
0.2
0.2
0.2
0
0
Never
Within a Year
Over a Year
0
Never
Within a Year
Over a Year
Never
Within a Year
Over a Year
多変量適応型回帰スプライン
mars.to.earth
でオブジェ
クトを変換
可能
mdaパッケージ
mars関数
mars(x, y, options)
mars関数は柔軟判別分析(Hastie et al., 1994)に応用された簡略版である．多
変量応答に拡張されている．
earthパッケージ
earth関数
earth(x, y, options) OR earth(formula, options)
Friedman (1991)を拡張している．ただし，glm関数を用いているため，glmで
解析可能な応答はすべて用いることができる．
Optionsの説明
mars
earth
degree
degree
nk
nk
一般化交差確認法におけるペナルティ
penalty
penalty
前進ステップワイズ過程において，隣接する区分点のあ
いだに必要とする最小個体数
thresh
thresh
最小スパン
-
minspan
一般化線形モデルの関数glm()での記述
-
glm
前進ステップワイズ法において，新たな変数を追加する
ときのペナルティ
-
newvar.penalty
交差確認法による評価
-
nfold
基底関数の次数
前進ステップワイズ過程での最大値基底関数個数
MARS法の問題点とその解決
MARS法では前進ステップワイズ法で構成される基底関数の最大個数および基底関数
の次数を任意に与えないといけない．
これらを評価するための関数を用意する．
seq_earth <- function(X, Y, degrees, Mmaxs){
ndeg <- length(degrees)
nMax <- length(Mmaxs)
ncol <- ndeg*nMax
GRres <- matrix(numeric(ncol), ncol=ndeg)
colnames(GRres) <- degrees
rownames(GRres) <- Mmaxs
Mres <- GRres
k <- 1
for (i in 1:ndeg){
for (j in 1:nMax){
md <- earth(X, Y, nk=Mmaxs[j], degree=degrees[i])
GRres[j,i] = md$grsq
Mres[j,i] = length(md$selected.terms)
k <- k+1
}
}
続く
MARS法の問題点とその解決
続く
layout(matrix(1:2,ncol=2))
matplot(Mmaxs, GRres, pch=1:ndeg, type="o", col=rainbow(ndeg)
, xlab=expression(M[max]),ylab="GCV R-squared")
legend(x=min(Mmaxs), y=max(GRres) , legend=degrees
, col=rainbow(ndeg), title="Degree", pch=1:ndeg)
matplot(Mmaxs, Mres, pch=1:ndeg, type="o", col=rainbow(ndeg)
, xlab=expression(M[max])
, ylab="Number of selected base function")
legend(x=min(Mmaxs), y=max(Mres) , legend=degrees
, col=rainbow(ndeg)
,title="Degree", pch=1:ndeg)
}
ここで紹介している関数は，
http://www.ccn.yamanashi.ac.jp/~shimokawa/R-tree/
に公表します．
MARS法の事例：earthでの出力
事例検討：体脂肪測定データ(Garcia et al, 2005)
75 名のドイツ人女性に対して，2 重エネルギーX 線吸収測定法(DXA) によって測定さ
れた体脂肪値が体系に関する特性値とともに記録されている．
変数名
age
DEXfat
waistcirc
hipcirc
elbowbreadth
kneebreadth
説明
年齢
DXA
腹囲(ウエスト)
腰囲(ヒップ)
肘幅
膝幅
このデータはパッケージmboostのデータbodyfatに入っている．
> data(bodyfat, package="mboost")
> y <- bodyfat$DEXfat
> x <- bodyfat[,c(1,3:6)]
seq_earthの実行
6.0
> y <- bodyfat$DEXfat
> x <- bodyfat[,c(1,3:6)]
> seq_earth(x,y,1:3, seq(10,50,by=10))
Degree
5.4
5.6
5.8
1
2
3
5.0
5.2
0.864
0.865
Number of selected base function
1
2
3
0.863
GCV R-squared
0.866
Degree
10
20
30
Mmax
40
50
10
20
30
40
50
Mmax
基底関数の次数2のときには，最大基底関数に依らずGCV寄与率，基底関数の個数と
もに変化がなく，最も高かった．
MARS法の実行と推定結果
> mars.mdl <- earth(x,y,degree=2)
# 基底関数の次数=2でのMARSモデルの構成
> summary(mars.mdl,style="pmax",digits=3) # 基底関数の出力
Call: earth(x=x, y=y, degree=2)
y =
38.1
- 0.422
- 0.621
+ 0.0503
+ 0.0164
+ 0.483
0.97
*
*
*
*
*
*
pmax(0,
pmax(0,
pmax(0,
pmax(0,
pmax(0,
pmax(0,
89.8
113
age
58
waistcirc
hipcirc
-
waistcirc)
hipcirc)
58)
age)
89.8)
113)
*
*
*
*
pmax(0,
waistcirc
pmax(0,
waistcirc
pmax(0, kneebreadth
pmax(0, elbowbreadth
-
89.8)
89.8)
9.8)
6.5)
Selected 7 of 21 terms, and 5 of 5 predictors
Importance: hipcirc, waistcirc, kneebreadth, elbowbreadth, age
Number of terms at each degree of interaction: 1 2 4
GCV 17
RSS 723
GRSq 0.863
RSq 0.915
■ 推定されたMARSモデルの結果
y=38.100-0.422 89.8-waistcirc  -0.621113-hipcirc   0.050  age-58  waistcirc-89.8
+0.0164 58-age   waistcirc-89.8 +0.483  waistcirc-89.8  kneebreadth-9.8
-0.970  hipcirc-113  elbowbreadth-6.5
MARS法の結果とグラフィカル表示：手書き
■ 推定されたMARSモデルの結果
y=38.100-0.422 89.8-waistcirc -0.621113-hipcirc   0.050  age-58  waistcirc-89.8
+0.0164 58-age  waistcirc-89.8 +0.483  waistcirc-89.8  kneebreadth-9.8
-0.970  hipcirc-113  elbowbreadth-6.5
■ 推定されたMARSモデルの樹木表現
38.100
waistcirc<89.8
hipcirc<113.0
-0.422
<58
0.0164
age
>58
0.050
-0.621
kneebreadth > 9.8
0.483
kneebreadth > 6.5
-0.970
earthオブジェクトの表示形式の違い (1/2)
> summary(mars.mdl,style="max",digits=3) # 基底関数の出力(maxで表示)
y =
38.1
- 0.422
- 0.621
+ 0.0503
+ 0.0164
+ 0.483
0.97
*
*
*
*
*
*
max(0,
max(0,
max(0,
max(0,
max(0,
max(0,
89.8
113
age
58
waistcirc
hipcirc
-
waistcirc)
hipcirc)
58)
age)
89.8)
113)
*
*
*
*
max(0,
waistcirc
max(0,
waistcirc
max(0, kneebreadth
max(0, elbowbreadth
-
> summary(mars.mdl,style="h",digits=3) # 基底関数の出力(marsと同じ表示)
coefficients
(Intercept)
38.05
h(89.8-waistcirc)
-0.42
h(113-hipcirc)
-0.62
h(age-58) * h(waistcirc-89.8)
0.05
h(58-age) * h(waistcirc-89.8)
0.02
h(waistcirc-89.8) * h(kneebreadth-9.8)
0.48
h(hipcirc-113) * h(elbowbreadth-6.5)
-0.97
89.8)
89.8)
9.8)
6.5)
earthオブジェクトの表示形式の違い (2/2)
> summary(mars.mdl,style="C",digits=3) # 基底関数の出力(maxで表示，変数が行列)
y =
38.053
- 0.42238
0.6206
+ 0.050304
+ 0.016364
+ 0.48276
- 0.97006
*
*
*
*
*
*
max(0,
max(0,
max(0,
max(0,
max(0,
max(0,
89.8
113
x[0]
58
x[1]
x[2]
-
x[1])
x[2])
58)
x[0])
89.8)
113)
*
*
*
*
max(0,
max(0,
max(0,
max(0,
x[1]
x[1]
x[4]
x[3]
- 89.8)
- 89.8)
- 9.8)
- 6.5)
> summary(mars.mdl,style="bf",digits=3) # 基底関数の出力
y =
38.1
- 0.422
- 0.621
+ 0.0503
+ 0.0164
+ 0.483
0.97
続き
*
*
*
*
*
*
bf1
bf2
bf3
bf5
bf4
bf7
*
*
*
*
bf4
bf4
bf6
bf8
続く
bf1
bf2
bf3
bf4
bf5
bf6
bf7
bf8
h(89.8-waistcirc)
h(113-hipcirc)
h(age-58)
h(waistcirc-89.8)
h(58-age)
h(kneebreadth-9.8)
h(hipcirc-113)
h(elbowbreadth-6.5)
MARSモデルにおける変数重要度
earthパッケージでは変数重要度(ここでは相対測度)を表示する関数としてevimpがある．
evimp(model)
関数evimpではオプションとしてtrimがある(DefaultはTRUE)．これは基底関数に含まれ
ない変数(変数重要度が0)を予め削除するか否かを規定する．
> ev <- evimp(mars.mdl)
> par(mfrow = c(1,2))
> barplot(ev[,4], ylab="Variable importance", main="GCV", col="blue") #GCVによる変数重要度
> barplot(ev[,6], ylab="Variable importance", main="RSS", col="blue") # RSSによる変数重要度
100
RSS
100
GCV
60
80
データにおける変数重要度
と見做すことができる．
0
20
40
Variable importance
60
40
20
0
Variable importance
80
予測の観点から見たときの
変数重要度と見做すことが
できる．
hipcirc
waistcirc
kneebreadth
elbowbreadth
age
hipcirc
waistcirc
kneebreadth
elbowbreadth
age
それぞれの基底関数の挙動
earthパッケージでは，個々の基底関数の挙動を省察するための方法と
してplotmo関数がある．ただし，変数対の場合は同じグラフとなる．
≧c
xa
＜c
> plotmo(mars.mdl)
同じグラフ
-0.422 89.8-waistcirc
80
90
100
110
2 waistcirc: kneebreadth
90
100
110
1 age: waistcirc
120
3 hipcirc: elbowbreadth
e
ag
70
2 hipcirc
20 30 40 50 60
20 30 40 50 60
1 waistcirc
-0.621113-hipcirc
y earth(x=x,y=y,degree=2)
130
c
ci r
ist
a
w
0.050  age-58  waistcirc-89.8
+0.0164 58-age  waistcirc-89.8
kn
dth
c
cir
rea
0.483  waistcirc-89.8  kneebreadth-9.8
hip
b
ee
c
cir
ist
a
w
elb
a
b re
w
o
dt h
-0.970  hipcirc-113  elbowbreadth-6.5
余談：PRIM法，XR法，AIM法のパッケージとその意味(1/3)
■ PRIM法 (パッケージ：prim)
prim.box(x, y, threshold, threshold.type, option)
x: 説明変数， y：応答， threshold：閾値
threshold.type：応答の方向 (0：全方向, -1：負方向，1：正方向)
体脂肪データにおいて，箱内応答平均が80パーセント点以上になる箱を計算
> prim.res <- prim.box(x, y, threshold=quantile(y,0.8), threshold.type=1)
> summary(prim.res)
box-fun box-mass threshold.type
box1
46.27846 0.1830986
1
box2*
27.30966 0.8169014
NA
overall 30.78282 1.0000000
NA
prim.boxでは，箱のルールを表示するための関数がない．下川・杉本・後藤(2013)で
は，この問題を解決するための関数を用意している．
余談：PRIM法，XR法，AIM法のパッケージとその意味(2/3)
■ XR法 (パッケージ：XReg)
Ｒにおける極限回帰法のパッケージXRegは，バージョン2.10.1以降では動作
しない．そのため，CRANにも存在しない．利用を希望する場合には，
R2.10.1をインストール後，XRegをアーカイブとして保存しているサイト
からダウンロードしなければならない．
> data(bodyfat, package="mboost")
> y <- bodyfat$DEXfat
> x <- bodyfat[,c(1,3:6)]
> outstepxr1=stepXR(y,x,ypred=y,xpred=x,niterfull=10,
+ niter=5,step.size=.3,kstep=8,kbest=4,penalty=3,survival=FALSE)
> varlist=outstepxr1$steplist[[outstepxr1$bestgcv]]
> outxrs=XReg(y=y,x=x,varlist=varlist,step.adapt=TRUE,
+ niter=20,limsize=.05,step.size=.5,survival=FALSE)
> etaxr=predictfun(xpred=x,varlist=varlist,outxrs$beta0list)$predicted
> levelout=level.set(quantile(etaxr,.80),outxrs,pred=TRUE,xpred=x)
> print(levelout$cutpoints)
余談：PRIM法，XR法，AIM法のパッケージとその意味(3/3)
■ AIM法 (パッケージ：aim)
応答：連続変数(単回帰型)
主効果のみ
関数
説明
交互作用
関数
説明
lm.main
AIMの実行
lm.interaction
AIMの実行
cv.lm.main
最適基底関数サイズの選定
cv.lm.interaction
最適基底関数サイズの選定
応答：2値(0,1) (ロジスティック回帰型)
主効果のみ
関数
説明
交互作用
関数
説明
logistic.main
AIMの実行
logistic.interaction
AIMの実行
cv.logistic.main
最適基底関数サイズの選定
cv. logistic.interaction
最適基底関数サイズの選定
応答：生存時間(Surv(time,event))(比例ハザード・モデル型)
主効果のみ
交互作用
関数
説明
関数
説明
cox.main
AIMの実行
cox.interaction
AIMの実行
cv.cox.main
最適基底関数サイズの選定
cv. ox.interaction
最適基底関数サイズの選定
結びに代えて：樹木構造接近法をＲで実行するための留意点
■ rpartパッケージ
 関数rpartだけでは最適樹木サイズを選定しない．pruneを用いて刈り込み
を行うことが必要である．
 樹木図の終結ふしをグラフィカルに表示するにはpartykitが便利である．
 1SEルールによる樹木選定を行う関数がないので，自分で用意しなければ
ならない．
■ partyパッケージ
 並べ替え検定のp値を用いるので，大規模データの場合には，有意水準α
の調整などを行わなければ複雑な樹木を構成する惧れがある(樹木サイズ
に対する評価がないため)．
 ただし，分岐の根拠がp値に依拠するため，樹木の分離度がわかり易い．
■ earth, mdaパッケージ
 MARS法を適用する場合にはearthパッケージを利用することが推奨され
る．
 最大基底関数の個数，基底関数の次数の選定はMARS法の一つの問題点で
ある．例えば，残差平方和に対するGCV推定値を用いて諸種の間隔で評
価することが必要になる．

Download Report