読了：Therneau, Atkinson & Mayo Foundation (2022) Rのrpartパッケージによる決定木・回帰木の舞台裏

Therneau, T.M., Atkinson, E.J., Mayo Foundation (2022) An Introduction to Recursive Partitioning Using the RPART Routines. Octover 21, 2022.

　決定木・回帰木を提供するRの標準パッケージrpartのビネット。たまにこのパッケージを使うとき、パラメータの意味を忘れちゃってて混乱するので、このたび目を通した。

1. イントロダクション
[略]

2. 表記
　分類問題について考える。
　クラス数を\(C\)、クラス\(i\)の事前確率を\(\pi_i\)とする。ある事例の予測子ベクトルを\(x\)として、その事例の真のクラスを\(\tau(x)\)とする。
　クラス\(i\)の事例を\(j\)に分類してしまったときの損失を\(L(i,j)\)とし、\(L(i,i)=0\)とする。

　モデルの終端ノード数を\(k\)とする。モデルのノード\(A\)が終端ノードだったとき、そのノードにいずれかのクラスが与えられるわけだが、そのクラスのことを\(\tau(A)\)と書く。
　観察事例数を\(n\), うちクラス\(i\)の事例数を\(n_i\)、ノード\(A\)の事例数を\(n_A\)とする。[明記されてないけど、クラス\(i\)の事例のうちノード\(A\)に落ちた事例数を\(n_{iA}\)としている]

　ノード\(A\)の確率を$$ P(A) = \sum_{i=1}^C \pi_i Prob(x \in A | \tau(x) = i) $$とする。その推定量は\(\sum_{i=1}^C \pi_i (n_{iA}/n_i)\) である。
　ノード\(A\)におけるクラスの条件付き確率を$$ p(i|A) = Prob(\tau(x) = i | x \in A) = \pi_i Prob(x \in A| \tau(x) = i) / Prob(x \in A)$$ とする。その推定量は\(\pi(n_{iA}/n_i) / \sum \pi_i(n_{iA}/n_i) \)である。
　ノード\(A\)のリスクを$$ R(A) = \sum_{i=1}^C p(i|A) L(i, \tau(A))$$ とする。ツリー\(T\)のリスクを$$ R(T) = \sum_{j=1}^k P(A_j) R(A_j)$$とする(\(A_j\)は終端ノード)。

3. ツリーの構築
3.1 分割基準
　ノード\(A\)をノード\(A_L, A_R\)に分割したとき、$$ P(A_L)R(A_L) + P(A_R)R(A_R) \leq P(A)R(A)$$ となることが証明できる。だからリスク減少を最大化するような分割を求めればよい。と思うでしょ？ そうは問屋が卸さない。
　たとえば、損失はすべて等しく、データの80%がクラス1だとしよう。いま、\(A_L\)の54%がクラス1, \(A_R\)の100%がクラス1だとする。リスクを最小にする予測は\(\tau(A_L) = \tau(A_R) = 1\)である。ってことは、分割してもリスクは減少しない。でも分割しないってのはおかしい。
　そこでrpartパッケージはどうしているかというと、ノードのimpurity指標を使っている。ノード\(A\)のimpurityを$$ I(A) = \sum_{i=1}^C f(p(i|A))$$とする。\(f(p)\)はなんらかの関数で、\(A\)が純粋なときは\(I(A)=0\)であってほしいから、つまり\(f(0) = f(1) = 0\)であるようなconcaveな関数にする必要がある。とりあえず、情報指標\(f(p) = -p \log (p)\)を使うという手と、Gini指標\(f(p) = p(1-p)\)を使うという手がある。どっちを使っても大差ない。
　まあとにかく、ノード\(A\)のimpurity \(I(A)\)が定義できたとして、impurityの減少$$ \Delta I = p(A)I(A) – p(A_L) I(A_L) – p(A_R)I(A_R)$$ を最大化するような分割を探せばいいわけだ。
　[ここで、\(C\)個のクラスを2個の上位クラスに分けるにはどうするかという話。rpartには実装されていないので省略]

3.2 損失の統合
　上の考え方には損失関数が含まれていない。CARTの世界では、impurity基準に損失を含める方法として、一般化Gini指標と修正事前分布という方法がある。rpartは後者を採用している。
　[一般化Gini指標の説明。省略]
　修正事前分布とは何か。リスクの定義は$$ R(A) = \sum_{i=1}^C p(i|A) L(i, \tau(A)) = \sum_{i=1}^C \pi_i L(i, \tau(A)) (n_{iA}/n_i)(n/n_A) $$ であった[なんだか真値と推定量がごっちゃになっていると思うけど、まあそれはいいとして]。ってことは、いま新しい事前確率\(\tilde{\pi}\)と新しい損失関数\(\tilde{L}(i,j)\)があっても、任意の二つのクラス\(i,j\)について\(\tilde{\pi}_i \tilde{L}(i,j) = \pi_i L(i,j)\) が成り立つならば、新しい事前確率と損失関数に差し替えても\(R(A)\)は変わらないわけだ。
　そこで、2クラスのときは、\(i\neq j\)のとき\(L_i = L(i,j)\)として、$$ \tilde{\pi_i} = \frac{\pi_i L_i}{\sum_j \pi_j L_j} $$ とする。3クラス以上のときは、\(L_i = \sum_{k \neq i} L(i,k)\)として同じ式を使う。[…中略…]

3.3 事例: ステージC前立腺がん
[略]

3.4 変数重要性
[いまちょっと関心ないのでパス]

4. 木の剪定
4.1 定義
　複雑性パラメータを\(\alpha (> 0)\)とする。終端ノード数を\(|T|\)として$$ R_\alpha(T) = R(T) – \alpha|T|$$ とする。[ああなるほど。終端ノード数に対するペナルティね]
　\(T_1, T_2\)がともに\(T\)のサブツリーで、\(R_\alpha(T_1) = R_\alpha(T_2)\)であるならば、いっぽうが他方のサブツリーであるということが示せる。つまり、\(\alpha\)が決まれば、\(R_\alpha(T)\)を最小化するツリーのうち最小のツリーがひとつ決まるわけである。以下ではこれを\(T_\alpha\)と呼ぶ。
　\(\alpha > \beta\)ならば、\(T_\alpha\)は\(T_\beta\)と同じか、\(T_\beta\)の厳密なサブツリーである。[…中略…]

4.2 交差妥当化
[CVで\(\alpha\)を決めるという話。パス]

4.3 事例: 確率的な数の再認問題
[パス]

5. 欠損データ
5.1 分割の選択
　rpartは、目的変数に値があり、少なくともひとつの独立変数に値がある事例ならば、モデリングに使用する。
　最大化したいのは$$ \Delta I = p(A)I(A) – p(A_L) I(A_L) – p(A_R)I(A_R)$$ である。第1項はいいとして、第2,3項をどうするか。
　仮に、欠損値があるせいで\(A_L, A_R\)のどっちに落としたらいいかわかんない事例を取り除いたとしよう。極端な話、どっちかに落とせる事例がもし2事例しかなかったら、\(A_L, A_R\)はどちらもpureだということになるので、そういう分割が選ばれてしまうことになる。これは困る。

　そこでrpartは代替変数surrogate variablesというアイデアを使う。
　たとえば、どうやって決めたのかは置いといて、「40歳未満と40歳以上に分割する」ということが決まったとする。するとrpartは次の2通りのやり方を考える。
　その1, 「多数派に流れる」ルール。年齢が欠損である事例は、「40歳未満」と「40歳以上」のうち人数の多いほうに送る。
　その2、年齢が欠損である事例について、それが40歳未満か40歳以上かを、他の独立変数を用いた決定木(再帰なし)によって予測する。事前確率や損失はなし、リスクは単に(誤分類数)/(事例数)である。こうして得た予測のことを代理変数という。
　代理変数はたくさんできるわけだが、分割先のいずれに対しても2事例以上を送るような奴しか採用しない。当然ながら、代理変数にも欠損がありうる。
　代理変数群と「多数派に流れる」ルールを、誤分類率が小さい順に並べる。で、まずひとつめの奴をつかって、年齢が欠損である事例をどんどんどっちかに送っていく。ひとつめの奴に欠損があって送れなかった事例については、ふたつめの奴を使ってどっちかに送る。さらに…という風に続けていく。願わくば、「多数派に流れる」ルールが使われる前に、すべての事例がどっちかに流されますように、という話である。

5.2 事例: ステージC前立腺がん
[略]

6. そのほかのオプション
6.1 プログラム・オプション
　rpart()の主な引数について解説しよう。

  • formula: lm()とかと同じ。右辺には連続変数と因子が使える。左辺が2水準以上のときは時間がかかる。[rpartは左辺がベクトルでないといけないみたいね。mvpart::mvpart()は行列でもいけるのだが]
  • data, weights, subset
  • method: classification, anova, poisson, exponentialのいずれか。
  • parms: オプション・パラメータのリスト。分類の場合、prior(事前分布のベクトル)、loss(損失行列)、split(分割指標。giniとinformation)を指定できる。
  • na.action: デフォルトはna.rpartで、反応変数が欠損、ないしすべての予測子が欠損の行が除外される。
  • control: コントロールパラメータのリスト。rpart.control()を使って指定する。主な引数は:
    • minsplit: 事例数がそれ未満のノードは分割しない
    • minbucket: 終端ノードの最小数
    • maxcompete: 採用しなかった分割も画面に出力する
    • xval: CVの数
    • usesurrogate: {0:欠損を持つ事例はその下のノードに送られない, 1:「多数派に流れる」ルールぬき, 2:上記の通り}
    • maxsurrogate: ノードあたりに保持される代理変数の最大数
    • cp: 複雑性パラメータ。分割のない木を\(T_1\)として \(\alpha = cp * R(T_1)\)となる。cp=1とすると分割は起きなくなる。回帰木の場合は、「全体の\(R^2\)の上昇がcp以上であるような分割がみつからなかったらストップしろ」という意味になる。[疑問: rpartの知ったこっちゃないだろうけど、もし多変量回帰木だったら\(R^2\)はどう定義するんだろうか？]

6.2 事例: 自動車Consumer Report
[略]

6.2 事例: 脊柱後彎症
[略]

7. 回帰
7.1 定義
　手続きを以下のように拡張する。

  • 分割基準。分類課題の場合はGiniか対数尤度だった。anovaの場合は、あるノードについて\(SS_T = \sum (y_i – \bar{y})^2\)とし、その子ノードについて\(SS_L, SS_R\)を求めて、\(SS_T – (SS_L + SS_R)\)を分割基準とする。
  • ノードの記述。分類課題の場合は、クラスの予測、そしてクラス確率のベクトルである。anovaの場合はノードの平均となる。
  • ノードの誤差。分類課題の場合は損失の予測である。anovaの場合は\(y\)の分散である。
  • 新事例についての予測誤差。anovaの場合は\(y_{new} – \bar{y}\)となる。
  • 初期化。適宜変更する。

7.2 事例: 自動車Consumer Report

7.3 事例: ステージC前立腺がん

8. ポアソン回帰
[いま関心ないのでパスするけど、いずれ読まないと…]

9. 描画オプション
[略]

10. その他のオプション
[交差妥当化のための関数 xpred.rpart() の説明。略]

11. テスト事例
[超簡単な事例について数値を手計算で示している。すごい]
——–
　わかりやすい説明で、ありがたいです。表記が統一されていなくてちょっとイラっとしたけど、