統計解析環境Rによる 樹木構造接近法の実践 下川敏雄 山梨大学 大学院医学工学総合研究部 分類樹木の適用例 4種類の酵素が,4種類の肝臓疾患患者に対して取られている (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 | tT 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.621113-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.621113-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.621113-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) Rにおける極限回帰法のパッケージ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 最適基底関数サイズの選定 結びに代えて:樹木構造接近法をRで実行するための留意点 ■ rpartパッケージ 関数rpartだけでは最適樹木サイズを選定しない.pruneを用いて刈り込み を行うことが必要である. 樹木図の終結ふしをグラフィカルに表示するにはpartykitが便利である. 1SEルールによる樹木選定を行う関数がないので,自分で用意しなければ ならない. ■ partyパッケージ 並べ替え検定のp値を用いるので,大規模データの場合には,有意水準α の調整などを行わなければ複雑な樹木を構成する惧れがある(樹木サイズ に対する評価がないため). ただし,分岐の根拠がp値に依拠するため,樹木の分離度がわかり易い. ■ earth, mdaパッケージ MARS法を適用する場合にはearthパッケージを利用することが推奨され る. 最大基底関数の個数,基底関数の次数の選定はMARS法の一つの問題点で ある.例えば,残差平方和に対するGCV推定値を用いて諸種の間隔で評 価することが必要になる.
© Copyright 2025 ExpyDoc