異質性治療効果を推定するための統計手法・機械学習手法の現況と医療分野への応用について
下川敏雄
和歌山県立医科大学 医学部 ／ 附属病院 臨床研究センター

今回のお話における問題設定と統計的課題◼ はじめに◼ 樹木モデルを基調とした治療効果モデル◼ Meta-learner◼ 変換アウトカム/修正共変量法◼ 共有基底条件付き平均回帰◼ おわりに

はじめに

今回のお話における問題設定と統計的課題
問題設定

Target Therapy (TT)OutcomeStandard Therapy (ST)

Outcome

治療効果 (TE; Treatment Effect)

患者

いずれか一方の治療しか施行できないため，TEの観測は不可能(Neyman-Rubinの反事実的事象)個別化医療に向けての統計的課題
個別化医療を実現するには，患者の医療情
Target Therapy (TT)

患者

医療情報
(背景データ)報(背景データ)から，最適な治療法を推論する必要がある．

Decision

そのため，個々の患者のTEを何らかの方法でStandard Therapy (ST)

推定しなければならない．

今回のお話における問題設定と統計的課題
Heterogeneous Treatment Effect

引用：Bica et al. (2021). From Real-world patient data to individualized treatment effects using machine learningcurrent and future methods to address underling challenges, Clinical Pharmacology & Theraputics, 109(1),87-100.

一般的に利用されるATE(Average Treatment Effect)は，すべての患者に対して等しい治療効果が得られると仮定する．

他方，被験者の背景因子等(Patient Characteristic)によって治療効果は異なるはずであり，共変量によって条件づけられた治療効果，すなわち，異質性治療効果(Heterogeneous Treatment Effect)が重要である．

RWD等を活用した個別化医療に対するアプローチ
Subgroup Identification (SI)
患者情報に基づいてサブグループを構成し，サブグルー
プ内の平均値の差を評価

治療効果モデル

サブグループ同定法
e.g. SIDES (Lipkovich et al., 2011)
QUINT (Dusseldorp et al., 2014)

推定HTEを応答としたSI(Treatment effect model)
患者情報xにおける(条件づけられた)もとでの治療効
果の条件付き期待値

SIは治療法から患者を診るためのアプローチ(解釈優先)

Unconfoundedness Assumption

は異質性治療効果(HTE; Heterogeneous Treatment
Effect)と呼ばれる．治療効果モデルは，HTEを推定するための統計モデルである．

TEMは患者から治療効果を予測するためのアプローチ(予測確度優先)

個(患者)

治療

治療

個(患者)

・TEが顕著なレスポンダーのプロ
ファイリング・治療効果によるClassification

治療効果モデルの諸型 (主なものだけに限定しています)
樹木モデルを基調とした方法
因果樹木法

変換アウトカム/修正共変量手法(Athey & Imbens, 2016)

因果Forest法

因果Boosting法(Wager & Athey, 2018)(Powers et al., 2018)

一般化RandomForest法

F-learner

修正共変量法(Künzel et al.2019)(Tian et al, 2018)

U-learner

A-learner(Künzel et al.2019)(Chen et al, 2017)(Athey et al., 2019)

Orthgonal RandomForest法(Oprescu et al., 2019)

Meta-learner
S-learner(Künzel et al.2019)

Virtural TWINS法 (Foster et al., 2011), BART法 (Hill, 2011)

共有基底条件付き平均回帰

T-learner

DR-learner

PTO-Foreset(Künzel et al.2019)(Kennedy, 2022)(Powers et al, 2018)

X-learner

DA-learner

因果MARS法(Künzel et al.2019)(EconML)(Powers et al., 2018)

R-learner

因果RuleFit法(Nie & Wager, 2019)(Wan et al., 2023)

樹木モデルを基調とした治療効果モデル

因果樹木法 (Athey & Imbens, 2016)
Athey & Inbens(2016)は，CART法(Breiman et al., 1984)では，ふし内平均(回帰パラメータの推定)が樹木構築(終結ふしの推定)の影響を受けることから，パラメータの独立性が担保できないことを問題視した．
このことを解決するための方策として，Honest推定を提案している．
データ集合
学習標本

Honest型

推定標本

樹木モデルの構成
終結ふしでの予測

テスト標本
評価

データ集合
樹木モデルの構成

Adaptive型

学習標本

終結ふしでの予測

テスト標本評価
Athey, S. & Imbens, G.W. (2016). Recursive partitioning for heterogeneous treatment effect. Proceedings of the National Academy of Science, 113(27), 735–7360.

因果樹木法 (Athey & Imbens, 2016)
Athey & Imbens (2016)は，樹木モデルを構築するための損失関数として，修正MSEを提案するとともに，その期待値(EMSE; 期待MSE)の負値を最大にすることを提案している．そのもとで，HTE推定に修正している．
因果樹木の樹木構築で用いられる損失関数(EMSE):学習標本の標本サイズ，:推定標本の標本サイズ，
：ふし

：ふし

：ふし

:ふし

における

の割合

における推定HTE(応答平均の差)における群kの応答平均における群kの不偏分散

EMSEの負値の最大化を行うことが，因果樹木法の樹木構築の目的である．

Athey, S. & Imbens, G.W. (2016). Recursive partitioning for heterogeneous treatment effect. Proceedings of the National Academy of Science, 113(27), 735–7360.

因果樹木法の構成アルゴリズム

樹木モデルの構築
学習標本

を用いる

学習標本においてEMSEが最小になるようにモデルを構築する推定HTEの計算
推定標本

を用いる

推定標本を樹木モデルに当てはめることで，各終結ふし(leaf)での推定標本における推定HTE(終結ふし内での応答平均の差)を計算し，因果樹木での推定値とする．

最適樹木の選定
テスト標本

を用いる

樹木系列のそれぞれに対して，テスト標本のもとで
EMSEを評価し，最小になる
樹木を最適サイズ樹木(right
sized tree)とする．

Athey, S. & Imbens, G.W. (2016). Recursive partitioning for heterogeneous treatment effect. Proceedings of the National Academy of Science, 113(27), 735–7360.

因果樹木法の医学研究への応用：Look AHEAD Study
これは，米国16施設において2型糖尿病を伴う肥満患者5145人を対象にカロリー制限および運動といった生活習慣に対する介入群とコントロール群に無作為に割り付けた無作為化比較試験である．主要評価項目は，Cause specific survival time(心血管系が原因の死亡，致死的心筋梗塞，非致死的脳卒中，狭心症)である．
体重

腰囲

体力(MET)

糖化ヘモグロビン

主論文(The Look AHEAD research group, 2013)では，介入によるダイエット効果が認められたが，それがイベント発現に差異を及ぼすには至らなかったと結論付けている
The Look AHEAD Research Group. (2013). Cardiovascular effect of intensive lifestyle intervention in type 2 diabetes, NEJM, 369, 145-154.

因果樹木法の医学研究への応用：Look AHEAD Study
疑問：どのような患者には積極的に栄養・運動介入を行うべきか？

Intervention

Control

ControlIntervention

つまり，HbA1cが4.8(腎機能がある程度保たれている)場合において，SF-36(患者の健康状態
QOLの指標)が48点以上あった場合には有効だ
が，そうでない場合には悪化させる．
Baum, A. et al. (2017). Targeting weight loss interventions to reduce cardiovascular complications of type 2 diabetes: a machine learning-based post-hoc analysis ofheterogeneous treatment effects in the Look AHEAD trial, Lancet Diabetes Endocrinol., 5(10), 808-815.

因果樹木法を基本学習器とするアンサンブル学習法：Causal Forest法
Wager & Athey (2018)は，因果樹木に基づくRandomForest法を提案している．そこでは，樹木モデル(base learner)として，
Double-sample trees及び，Propensity treesを提案している：

Double-sample trees
Double-sample treesでは，観測値 L = { X n , Yn , Z n }n =1 をランダムに半分に分け( L = {I , J } )，次のようなアルゴリズムによってモデルを構築している．
N

STEP.1：樹木を成長させる．このとき，標本 J を用いる分岐候補の選定には，Causal treeのEMSEを用いる．
STEP.2：各終結ふしにおける応答の予測値 ˆlの計算には，Iのみを用いる．

→ したがって，Athey & Imbens (2016)のHonest型のCausal treesと同じ

Propensity trees

Propensity treesでは，群 Z n を応答と見做したもとで樹木を構築したもとで，終結ふしに対する治療効果の予測値 ˆlを応答に基づいて推定する．
STEP.1：樹木を成長させる．このとき，( X n , Z n )に対して，分類樹木(Breiman et al., 1984)によって樹木を構築する．
STEP.2：各終結ふしにおける応答の予測値 ˆlを計算する．

なお，オリジナルでは，サブサンプルを用いるような記法になっている(後述のCausal forestに拡張するため)．
Causal forestは，double-sample treeあるいはpropensity treeをbase learnerとしたもとでrandom forestのアルゴリズムによってでアンサンブルさせた方法である．ただし，ブートストラップ標本は，重複なしとする．
さらに，Wager & Athey (2018)では，Causal Forest法での漸近正規性および，分散推定量を提案している．
Wager, S. & Athey, S. (2018) Estimation and inference of heterogeneous treatment effects using random forests, J. Amer. Statist. Assoc., 113(523), 1228-1242.

因果樹木法を基本学習器とするアンサンブル学習法：因果Boosting法
Powers et al. (2018)は，因果樹木を基本学習器(base learner)に用いるBoostingとして，因果Boosting法を提案している．アルゴリズムは，通常の勾配ブースティング(Friedman,2001)と同じ
因果Boosting法のアルゴリズム
Require: ,：基底関数の重み，

Initialize:

,

For

do

in
データ集合

：Boosting回数：疑似残差，に因果樹木：因果Boostingによる推定治療効果モデルを当てはめる

疑似残差の更新
推定治療効果モデルの更新

End for
ReturnPowers et al. (2018)は，最適な樹木数を計算するための方法として，hold-out validationを提案している．

Friedman, J.H. (2001). Greedy function approximation: a gradiet boosting machine. Ann. Statist., 29(5), 1189–1232.
Powers, S., et al. (2018). Some methods for heterogeneous treatment effect estimation in high dimensions. Statistics in Medicine, 37(11), 1767-–1787.

因果樹木法を基本学習器とするアンサンブル学習法：因果Boosting法
hold-out validation

を用いてBoosting樹木を構築

Causal Forest法と異なり，Causal BoostingではBase learner (因果樹木)が増えるほど過学習に陥る可能性がある．
ただし，通常の交差確認法を用いることができないため，Powers et al.(2018)では，新たに，hold-out validation法を提案している．

Step.2

Step.1を用いてBoosting樹木を構築

Data set

k=1

k=2

ｖ

K

・・・

Test
sampleStep.2で構築した樹木の終結ふしの推定
HTEをテスト標本で再推定し，モデルを推定する．

・・・kR

Training
sampleStep.3k個の樹木から構成される因果Boostingモ
デルに対するvalidation errorを計算する．

k=1での推定HTE：
k=2での推定HTE：
・・・

・・・
・・・

k=K での推定HTE：

が最小となる樹木数Powers, S., et al. (2018). Some methods for heterogeneous treatment effect estimation in high dimensions.
Statistics in Medicine, 37(11), 1767–1787.

を選定する．

meta-learner

Meta-learnerの諸型
S (Single)-learner(Künzel et al., 2019)
いま，

を共変量

および治療指標 (:新規治療，:既存治療)による応答に対する回帰モデルとするとき，S-learnerによるHTEはで与えられる．
つまり，S-learnerとは，共変量に治療指標 を加えた回帰モデルを推定し，そのもとで，HTEを推定するだけの卑近な方法である．
S-learnerに基づく方法としては，■ BART (Bayesian Additive Regression Trees; Chipman et al., 2007, 2010)を用いた方法 (Hill, 2011)■ Random Forest (Breiman, 2001)を用いた方法：Virtual Twins(Foster et al., 2011)などがある．
Künzel, S.R., et al. (2019) Metalearners for estimating heterogeneous treatment effects using machine learning, PNAS, 116(19), 4156-4165.
Hill, J. L. (2011) Bayesian Nonparametric Modeling for Causal Inference, J. Compt. Graph. Statist., 20(1), 217-240.
Foster, J.C. et al. (2011) Subgroup identification from randomized clinical trial data, Statist. Med., 30(24), 2867-80.

Meta-learnerの諸型
T (Two) -learner(Künzel et al., 2019)
いま，新規治療群に対する回帰モデル，既存治療群に対する回帰モデルが与えられたとき，T-learnerによるHTEは，で与えられる．
つまり，T-learnerとは，治療指標 によりデータを2群に分けたうえで，各群において，回帰モデル，を推定したうえで差をとるだけの卑近な方法である．

Powers et al.(2018)は，T-learnerを条件付き平均(conditional mean regression)と称している．

S/T-learnerの問題点■ S/T-learnerはHTEを推定するのみであり，(推定)HTEに対する要因構造が提示できない．
■ S/T-learnerでは，治療法に対する選択バイアスを調整する手段は考慮されていない．したがって，いわゆる観察研究と同様に取り扱われるRWD(Real World Data)では不適切な結果をもたらす恐れがある．
Powers et al. (2018) Some methods for heterogeneous treatment effect estimation in high-dimensions, Statist. Med., 37(11), 1767-1787.

Meta-learnerの諸型
X-learner(Künzel et al., 2019)
新規治療群に対する回帰モデル：，既存治療群

新規治療群のアウトカム・共変量：,実施された治療
既存治療

治療群

新規治療

に対する回帰モデル：，既存治療群のアウトカム・共変量：,

疑似アウトカム(pseudo outcom)

を応答，

を説明変数とした回帰モデルを当てはめる：

を応答，

を説明変数とした回帰モデルを当てはめる：

傾向スコア：

治療効果モデル：Künzel, S.R., et al. (2019) Metalearners for estimating heterogeneous treatment effects using machine learning, PNAS, 116(19), 4156-4165.

Meta-learnerの諸型
DA-learner(EconMLに実装)
X-learnerの問題点■ X-learnerではおよびの推定において選択バイアスが考慮されていない．
■ X-learnerの治療効果モデルは，(推定)HTEに対する影響要因を評価できない．
DA-learnerは，上記の問題点を解消する形式で開発されたX-learnerの派生である．
重みの計算

各群の応答に対する回帰モデルの重み付き推定

傾向スコア

疑似アウトカム

治療効果モデル：： に対して何らかの回帰モデルを当てはめ
https://econml.azurewebsites.net/

Meta-learnerの諸型
DR-learner(Kennedy, 2022)
T-learnerの問題点■ T-learnerではおよびの推定において選択バイアスが考慮されていない．
■ T-learnerの治療効果モデルは，(推定)HTEに対する影響要因を評価できない．
DR-learnerは，2重ロバスト推定の枠組みで拡張された方法である．
新規治療群に対する回帰モデル：，既存治療群に対する回帰モデル：

傾向スコアとするとき，HTE(疑似アウトカム)における2重ロバスト推定はで与えられる．DR-learnerにおいても，DA-learnerと同様に疑似アウトカム に対して何らかの回帰モデルを当てはめ()，それを治療効果モデルとしている．
なお，DR-learnerの2重ロバスト性等の数理的な議論については，Kennedy (2022)を参照
Kennedy, E.H. (2022). Towards optimal doubly robust estimation of heterogeneous causal effects, arXiv:2004.14497v3

Meta-learnerの諸型
R-learner
Nie & Wager (2021)は，Robinson変換に基づく新たなmeta-learnerとしてR-learnerを提案している．

モデル：

ただし，

(unconfounderedness assumption)
このとき，上記モデルの条件付き平均関数(両辺に期待値をとる)とき，である．ここに，は傾向スコアである．この2つの式の差分を取ることで，となる(Robinson変換)．したがって，R-learnerによる治療効果モデルは，である．Nie & Wager (2020)は，正則化項を付与したもとで，経験的推定値をのように定義している．
Nie, X. & Wager, S. (2021). Quasioracle estimation of heterogeneous treatment effects. Biometrika, 108(2), 299-319.

Meta-learnerにおける過学習に対する対応：Cross-fitting approach
DA- /DR- /R- learnerでは，■ アウトカムを推定するための回帰モデル (アウトカム・モデル)■ 傾向スコアを推定するための回帰モデル (傾向スコアモデル)■ HTE (疑似アウトカム)を推定するための回帰モデル (治療効果モデル)の3種類の回帰モデルを推定している．このとき，過学習によるバイアス(overfitting bias)を解決するための方法
として，cross-fitting approachが存在する．
TEM.1:

HTEの推定には各モデルから得られた推定HTEの平均値を用いる

TEM.2:

Data
Set

TEM.3:
TEM.4:
・・・

・・・
ランダムに
M個に分割

Jacob(2020)の数値事例

アウトカム・モデル
，傾向ス
コアモデル
，治療効果モデルを異なるサブサンプルを用いて推定する．

Jacob, D. (2020) Cross-fitting and averaging for machine learning estimation of heterogeneous
treatment effects, ITRG 1792 Discussion Paper, 2020-014.

Meta-learnerに対するシミュレーション評価 (Okasa, 2022)
Okasa, G 2022). Meta-learners for estimation of causal effect: finite
sample cross-fit performance, arXiv:2201.12692

データ生成
共変量：

P：共変量の数(100個),

傾向スコア：

治療指標：:バランスを調整するパラメータ(モデルで異なる)誤差：
既存治療での真のアウトカム：
新規治療での真のアウトカム：
Simulation 2：

Simulation 5：

HTE・治療指標のバランスのパターン

新規治療群の割合：50%

Simulation 1：
新規治療群の割合：50%

Simulation 4：
新規治療群の割合：25%

Simulation 3：
新規治療群の割合：15%

Simulation 6：

新規治療群の割合：25%

新規治療群の割合：25%

評価尺度

C：シミュレーション回数 (10,000回)

推定HTEの標準偏差

RMSE

バイアス

Jarque-Bera統計量

尖度，歪度に関する統計量

歪度

手法

尖度

・S-learner (S)・T-learner (T)・X-learner (X-F)
・DR-learner (DR-F)

・標本を分割してのX-learner (X-S)
・標本を分割してのDR-learner (DR-S)
・標本を分割してのR-learner (DR-S)
・Cross fittingを伴うX-learner (X-C)
・Cross fittingを伴うDR-learner (DR-C)

・R-learner (DR-F)

・Cross fittingを伴うR-learner (DR-C)・初期分岐に治療指標を用いるS-learner (S-W)

Note：すべての回帰モデルにはRandomForest法(樹木数100)を用いる．
Okasa, G 2022). Meta-learners for estimation of causal effect: finite sample cross-fit performance,
arXiv:2201.12692

Meta-learnerに対するシミュレーション評価 (Okasa, 2022)
Model.6 での平均
標本サイズ

その他のモデルはSupplementary

・標本サイズが小さい場合には， 通常のX-learner(X-F)での性能が最良であるが，標本サイズが大きい
場合には，Cross-Fittingを伴うDR-learner(DR-C)が最良だった．
・S-learnerにおいて初期分岐に治療指標を含んだ場合，T-learnerと性能がほぼ変わらなかった．
・標本サイズが小さい場合にサブサンプルに分ける(因果樹木の形式)では，性能が著しく低下する．
Okasa, G 2022). Meta-learners for estimation of causal effect: finite sample cross-fit performance, arXiv:2201.12692

変換アウトカム／修正共変量法

変換アウトカム法

Tian, L., Alizadeh, A.A., Gentles, A.J., Tibshirani, R.: A simple method for estimating interactions between treatment
and a large number of covariates, Journal of the American Statistical Association, 109, 1517-1532, 2014.

ここでは議論を簡略化するため，線形モデルで考える．
いま，治療指標Zを以下のように表す：,(1：新規治療，-1：既存治療)

交互作用を含んだ次のような線形モデルを考える

共変量の影響

治療の影響

治療×共変量の影響

任意の共変量X=xが与えられたときのHTEΔ(x)は次のように表すことができる．
( x ) = E(Y (1) − Y ( −1) | X = x )( x ) = E(Y (1) − Y ( −1) | X = x ) = E(Y | X = x , Z = 1) − E(Y | X = x , Z = −1)= E(Y | X = x , Z = 1) − E(Y | X = x , Z = −1)となり，治療効果を予測するための線形モデルになる．また，Zが確率1/2で1,あるいは-1をとることから，で表すことができる．このことは，アウトカムを2YZに変換したうえで，通常の重回帰モデルの取り扱いで計算できることを意味する．そのため，この方法は，修正アウトカムモデルと呼ばれる．

変換アウトカム法
観測値が与えられたとき，Tian et al. (2014)による変換アウトカム法では，を応答として回帰モデルf(x)を推定する：

ただし，Tian et al. (2014)の方法は，確率1/2で割り付けられたRCTのみに焦点が当てられている．
Künzel et al. (2019) は，傾向スコアを用いることで，観察研究に対する変換アウトカム法として，二つの方法を紹介している．
F-learner：

U-learner：,,:推定傾向スコア:治療指標を伴わないxのみでの回帰モデル

なお，U-learnerは，R-learnerの縮小推定を伴わない場合と同様である．
Künzel, S.R., et al. (2019) Metalearners for estimating heterogeneous treatment effects using machine learning, PNAS, 116(19), 4156-4165.

修正共変量法
Tian et al.(2014)は，修正アウトカムモデルの拡張系として，2値データ，生存時間データに利用できる修正共変量手法を提案している．
いま，次のような回帰モデルを考える：
観測値が与えられたとき，最小2乗損失Lossは以下のように定義される．

緑色の部分に注目すると，応答
り，治療効果は，に対して，±1/2を掛け合わせて最小2乗法を行ったとみることができる．つまActiveControl
Treatment effect

つまり，修正共変量法では修正アウトカムモデルの半分で考えていることと同じである．

なので，最小2乗損失Lossに基づくパラメータ推定は，治療効果を推定する修正アウトカムモデルと同じになる．

アウトカムを円関する代わりに，共変量に重みづけする方法を修正共変量法という．

修正共変量法
修正共変量法の特徴は，2値応答あるいは生存時間データへの拡張が可能な点である．基底関数をW(x)とするとき，各アウトカムに対する修正共変量法のモデルは下記の通り：
HTE：

連続応答：
2値応答：
HTE：

Cox比例ハザード・モデルに基づく生存時間データへの拡張は，Tian et al.(2014)を参照されたい．
Chen et al. (2017)は，傾向スコアを用いることで，修正共変量法を観察データに拡張した方法として，
A-learnerを提案している．

ここに，は傾向スコアである．

Chen, S., et al. (2017). A general statistical framework for subgroup identification and comparative treatment scoring, Biometerics, 73(4), 1199-1209.

実際例(2)：局所進行性上咽頭癌患者に対する化学放射線療法に対する導入化学療法の上乗せ効果のレスポンダー探索
局所進行性上咽頭癌患者に対する導入化学療法の効果には，個人差があることが無作為化比較試験で明らかになった．この研究では，コホート
データを用いるとともに，修正共変量手法を用いて導入化学療法スコアを構築するとともに，その妥当性について検討している．

学習標本

導入化学療法推奨群

導入化学療法非推奨群

修正共変量手法の治療効果スコアに基づいて3 群に分けたときの上下2 グループでのKaplan-Meierプロット

推定されたモデルを用いることで，放射線化学療法に対する導入化学療法の推奨度を測ることができる．

テスト標本

メモ
本研究のもとのデータは，コホート研究に基づいており，本来は，傾向スコアの逆確率重み付け(IPW)などを用いて事前に加重化などの調整を必要があったと思われる(検証標本にRCTを利用)．
Dong, D. et al.: Development and validation of a novel MR
imaging predictor of response to induction chemotherapy inlocoregionally advanced nasopharyngeal cancer: a randomized
controlled trial substudy (NCT01245959) , BMC Medicine, 17:190.
2019.

共有基底条件付き平均回帰治療効果モデルに対する新たなアプローチ：共有基底条件付き平均回帰■ 条件付き平均回帰はHTEを推定するのみであり，(推定)HTEに対する要因構造が提示できない．
■ 条件付き平均回帰では，治療法に対する選択バイアスを調整する手段は考慮されていない．
Powers et al. (2018)は，条件付き平均回帰(T-learner)の欠点を解決するためのフレームワークとして，共有基底条件付き平均回帰を提案している．
新規治療に対する回帰モデル

条件付き平均回帰(T-learner)

既存治療に対する回帰モデル,
共有基底条件付き平均回帰では，同じ基底関数をもつことを仮定している．

HTE：

アウトカムに対する回帰モデル

共有基底
条件付き平均回帰HTE：

共有基底条件付き平均回帰では，(1) 基底関数の推定， (2) パラメータ推定，2段階での推定を行われる．
Powers et al. (2018) Some methods for heterogeneous treatment effect estimation in high-dimensions, Statist. Med., 37(11), 1767-1787.

共有基底条件付き平均回帰の諸型：PTO Forest法 (Powers, et al., 2018)
基底関数(複数の樹木)の推定
傾向スコアの推定

変換アウトカムの計算

基底関数の推定回帰パラメータ(leafにおけるHTE)の推定
b番目の樹木
推定値は下式で与えられる
m番目のleaf
Leafに属する部分集合：
Leafに属する個体数：()

Powers et al. (2018)は，このプロセスを受粉(Pollinate)と称している．

Leafのルール；
PTO Forest法の推定モデル：

共有基底条件付き平均回帰の諸型：因果MARS法 (Powers, et al., 2018)
MARS(Multivariate Adaptive Regression Spline; Friedman, 1996)法とは，打ち切りベキ乗基底関数に基づく非線形回帰手法である．

■MARS法のモデル： f MARS ( x ) =JmM   H (sm =1m

( j ,m )( x p ( j ,m ) − c p ( j ,m ) ))j =1

c p ( j ,m )詳しくは付録参照

打ち切りべき乗関数(括弧内が負値ならば0,正値ならばそのまま)

[ x p ( j ,m ) − c p ( j ,m ) ]+H ( s( j ,m ) ( x p ( j ,m ) − c p ( j ,m ) )) = [c p ( j ,m ) − x p ( j ,m ) ]+x p ( j ,m ), s( j , m ) = 1, s( j ,m ) = −1

MARS法は，(1) 前進ステップワイズ法，(2)後退ステップワイズ法(モデル系列の構成)，(3) 最適モデルの選択(通常はGCV(Generalized Cross-Validation))，で構成される．
前進ステップワイズ法におけるモデル構築課程では，以下の制約がある．
MJm  H (sf MARS ( x ) =  m(1) 交互作用を伴うため，基底関数の積
m =1( j ,m )( x p ( j ,m ) − c p ( j ,m ) )) のなかに，同じ変数は含まない．

j =1(2) 基底関数は [ x p ( j ,m ) − c p ( j ,m ) ]+ , [c p ( j ,m ) − x p ( j ,m ) ]+ の2個を追加する．

(3) 各ステップにおける基底関数の選定および回帰パラメータの推定の損失関数には残差平方和が用いられる．

Friedman, J.H.: Multivariate adaptive regression splines, The Annals of Statistics, 19(1), 1-67, 1991.

共有基底条件付き平均回帰の諸型：因果MARS法 (Powers, et al., 2018)
因果MARS法のモデル

既存治療群

新規治療群

因果MARS法においても，(1) 前進ステップワイズ法，(2) 後退ステップワイズ法，の順番でモデル推定を行う．ただし，前進ステップワイズ法において，基底関数の推定と回帰パラメータの推定は，個別に行われる(2段階推定)．

因果MARS法のアルゴリズムの略図
Step 1：基底関数の選定
基底関数の追加ペア候補

基底関数の
追加ペア

M回繰り返す直前のモデルでの残差を応答として，ステイジワイズの形式で当てはめ
主効果のペアのモデル

治療指標を加えたモデル最大となるペアを選定

Step 2：パラメータ推定

Step 1で選定された基底関数を加えたもとで，通常最小2乗法により回帰パラメー
タを推定する

残差推定されたモデルのもとで残差を計算する．

共有基底条件付き平均回帰の諸型：因果MARS法 (Powers, et al., 2018)
因果MARS法では，後退ステップワイズ法および交差確認法(cross-validation)によるモデル選択を行うが，Powers
et al. (2018)は，代替として，Baggingによるアンサンブルを行うことを推奨している．
数値検証のモデル (シミュレート回数：100回)

CM(1)：因果MARS (加法)CM(2)：因果MARS (2次交互作用)
BCM(1)：Bagging因果MARS (加法)BCM(2)：Bagging因果MARS (2次交互作用)

P=20

CM(1)：因果MARS (加法)CM(2)：因果MARS (2次交互作用)
BCM(1)：Bagging因果MARS (加法)BCM(2)：Bagging因果MARS (2次交互作用)

P=50

シミュレーションによる評価(Powers et al., 2018より引用)
Algorithm
1: 共変量の生成：X pn ~ N(0,1)，治療群の生成：Z n ~ Bernoulli( ( x))2: 主効果(main effect)・局所治療効果(local treatment effect)による応答の生成関数は表を参照- main effect function ( x)
- local treatment effect function ( x)
- Generate： N (  ( X n ) + ( Z n − 1/ 2) ( X n ),  Y2 )各群の平均アウトカムは 1 ( x) =  ( x) +  ( x) / 2, 1 ( x) =  ( x) −  ( x) / 2
分散 Y (ノイズ)は表を参照
2

Simulation Senario
Functionsf1 ( x) = 0,
f 2 ( x) = 5  I ( x1  1) − 5,
f 3 ( x) = 2 x1 − 4,f 4 ( x) = x2 x4 x6 + 2 x2 x4 (1 − x6 ) + 3 x2 (1 − x4 ) x6 + 4 x2 (1 − x4 )(1 − x6 )1, 9
2, 10 3, 11+ 5(1 − x2 ) x4 x6 + 6(1 − x2 ) x4 (1 − x6 ) + 7(1 − x4 )(1 − x4 ) x6n+ 8(1 − x2 )(1 − x4 )(1 − x6 ),
200
200
300
f5 ( x) = x1 + x3 + x5 + x7 + x8 + x9 − 2,
400
400
300
p
f 6 ( x) = 4 I ( x1  1) I ( x3  0) + 4 I ( x5  1) I ( x7  0) + 2 x8 x9 ,
 ( x) f8(x) f5(x) f4(x)
1 2
2
2
2
2
f 7 ( x) = ( x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 − 11) ,2
 ( x) f1(x) f2(x) f3(x)1
f8 ( x) =( f 4 ( x) + f5 ( x))
1
1/4
1
 Y2
2

Scenarios

4, 12

5, 13

6, 14

7, 15

8, 16

300

400

400

1000

1000

300

200

200

100

100

f7(x)

f3(x)

f1(x)

f2(x)

f6(x)

f4(x)

f5(x)

f6(x)

f7(x)

f8(x)

1/4

1

1

4

4

Powers et al.:Some methods for heterogeneous treatment effect estimation in high dimensions, Statistics in Medicine, 37, 1767-1787, 2018.

Note：詳しい結果は付録を参照

シミュレーションの考察■ RCTを想定した場合• Scenario 4, 5, 7, 8ではBoosted因果MARSが最良の結果を示した(Scenario 4においてもPTO Forestと同程度)• Scenario 2, 3, 4, 6ではPTO Forestが最良の結果を示した．
Null：Null model，TO：Transformed outcome forest，DB：Different-basis forest
CF：Causal forest，PTO0：PTO forest ，CB0：因果Boosting，BCM0：Boosted因果MARS

Senario

主効果

治療効果

1位

順位
2位

3位

1

線形

NULL

CB0

NULL

CF

2

線形

指示関数

PTO0

CB0

BCM0

3

非線形

線形

PTO0

TO

BCM0/CF

4

非線形

非線形

PTO0/BCM0

-

CB0

5

線形

線形

BCM0

CB0

PTO0

6

NULL

指示＋交互作用

PTO0

CB0

BCM0

7

指示関数

非線形

BCM0

PTO0

TO

8

非線形

指示＋交互作用BCM0

PTO0/CB0-Note：詳しい結果は付録を参照

シミュレーションの考察■ 観察研究を想定した場合• 因果Boosting(CB1)が8個の数値検証のうち，6個において最良の結果を示した．
Null：Null model，CF：Causal forest，PTO0：PTO forest ，CB0：因果Boosting(PS調整なし)，
CB1：因果Boosting(PS調整あり)， BCM0：Boosted因果MARS(PS調整なし)，Boosted BCM1(PS調整あり)

主効果

Senario

治療効果

1位

順位
2位

3位

9

線形

NULL

CB1

BCM1

CF

10

線形

指示関数

CB1

BCM1

CB0

11

非線形

線形

CB1

BCM1

CF

12

非線形

非線形

CB1

CB0

PTO/CF

13

線形

線形

BCM0

BCM1

CB1

14

NULL

指示＋交互作用

PTO

CB1

CB0

15

指示関数

非線形

CF/PTO-

BCM1/BCM016

非線形

指示＋交互作用CB1

CF/BCM1-

いずれのScenarioにおいても因果Forest(Causal forest)が惨憺たる結果になっているが，Generalized Forest法
において，その性能が大幅に向上しているので，そちらとの比較になると，その限りではないと思われる．

共有基底条件付き平均回帰の諸型：因果RuleFit法 (Wan, et al., 2023)
RuleFit法(Friedman & Popescu, 2008)
回帰パラメータ

事後重み付けアンサンブル・モデル

・・・

RuleFit法

ルールに対してlasso等を用い・・・ て重みづけを行う．

各樹木にlasso等を用いて重みづけを行う

ふしを含めてルールに分割

RuleFit法のモデル：

Boosting樹木を用いて得られたルール

修正線形項：しきい値(qおよび(1-q)分位点)

回帰パラメータαk，βpの推定にはlassoあるいはElastic Netなどの縮小推定を利用
Friedman, J.H. & Popescu, B.E. (2008). Predictive learning via rule ensemble, Ann. Appl. Stat. 2(3), 916-954.

共有基底条件付き平均回帰の諸型：因果RuleFit法 (Wan, et al., 2023)
因果RuleFit法のモデル

Step 1：ルール項の探索
ルール項の探索には，変換アウトカムを応答としたGradient Boosting Machine (Friedman, 2001)を用いる．
変換アウトカム：

Step 2：回帰パラメータの推定

Group

Group

同一基底関数のもとでの回帰パラメータをグループとしたう
えで，Group lasso (Yuan & Lin, 2007)を用いて回帰パラメー
タを推定する．
これにより，基底関数自体が刈り込める．
Wan, K., Tanioka, K., Shimokawa, T. (2023) Rule ensemble method with adaptive group lasso for heterogeneous treatment effect estimation, Statist. Med. in press.

因果RuleFit法の利点：基本学習器の評価

因果RuleFit法では，基本学習器(ルール or 修正線形項)に基づいて，「どのようなルールがHTEに影響を及ぼしているのか」を評価できるため，HTEの予測とサブグループの検討を同時に行うことができる．
Wan, K., Tanioka, K., Shimokawa, T. (2023) Rule ensemble method with adaptive group lasso for heterogeneous treatment effect estimation, Statist. Med. in press.

さいごに治療効果モデルに対する最近の展開：time-to-eventデータへの拡張
生存時間解析では，生存時間の推定に関心があるとは限らない(むしろ少ない)．

そのため，様々なHTEの定義があり，それらに対する治療効果モデルが提案されている．
HTEの
定 義

生存時間の差Hu et al. (2021)
Xu et al. (2022)

生存率の差
Tabib & Larocque(2017)

MSTの差
Zhang et al. (2017)hRMSTの差
Cui et al. (2023)

Zhu & Gallego. (2020)
Xu et al. (2022)• Hu et al. (2021)は，加速モデルのもとで拡張されたBART法(Henderson, 2020)，Random Survival Forest法(Ishiwaran et al., 2008)，
Deep NN(Katzman et al., 2018)をS-learnerのもとで，治療効果モデルに拡張している．
• Xu et al.は，Meta-learnerにおいて，S-/T-は生存率の差，X-/R-生存時間の差として，生存時間に拡張する方法を提示している．
なお，前者ではRandom Survival Forest，後者ではIPCW(Inverse Probability Censoring Weighting)を伴うRF法を用いている．
• Tabib & Larocque(2017)は，平均生存時間(KMで推定したうえで右側打ち切りは指数分布で補間)の差としてHTEを定義したうえでRF法によって治療効果モデルを構築している．

• Zhu & Gallego(2020)は，T-learnerの枠組みでSuper Learner(van der Laan et al., 2007)を当てはめている．
• Cui(2022)は，一般化Random Forest法の損失関数にRMSTを用いるとともに，IPCWによる重み付けを伴う形式で生存時間データへ拡張している．

治療効果モデル(treatment effect model)とサブグループ同定法(subgroup identification)― 個人的な意見 ―サブグループ
同定法

治療効果モデル• サブグループ内の群間差に基づいて分類する．
• Responderを抽出する．

臨床研究への応用が存在する．

探索的サブグループ解析(Lipkovich et al, 2016)

HTEを推定する

HTEに対する予測確度等を重要視するため，アンサンブル樹木あるいはStackingなど，モデル
をブラックボックス化する手法が研究の主題になっている．
要因構造の探索を行う場合には，結果の解釈が困難なため，部分従属プロット，変数重要度による解釈に終始する (解釈が主観的になる)．
臨床への応用が難しい場面が少なくないサブグループ同定法(subgroup identification)の臨床応用例[臨床応用例(1)]「重度の制御不能な喘息患者に対するTralokinumabの有効性・安全性に関する無作為化比較第III相試験」に対するレスポンダー探索
Gotthlow et al.(2019)は，STRATOS1およびSTRATOS2試験における5つのバイオマーカー(Periostin(喘息と関連がある遺伝子)，
DPP-4(消化管ホルモンの不活性を行う酵素)，Eosinophilis(好酸球，アレルギー反応の制御を行う白血球の一種)，FeNO(呼気中一酸化窒素濃度)，IgE(免疫グロブリンE)について，TralokinumabのレスポンダーをSIDES法を用いて検討している．
Lipkovich I, Dmitrienko A, Denne J, EnasG. Subgroup identification based on differential effect search (SIDES): A recursive partitioning method for establishing response totreatment in patient subpopulations. Statistics in Medicine 2011, 30:2601-2621.
Gottlow, M. et al: Application of structured statistical analyses to identify a biomarker predictive of enhanced tralokinumab efficacy in phase III clinical trials for severe,
uncontrolled asthma, BMC Pulmonary Medicine, 19(129), 2019.

[臨床応用例(2)]「小細胞肺癌に対するNavitoclaxの単群第II相試験」のデータを用いたバイオマーカーの最適カットオフ値の探索
Navitoclaxは，小細胞肺癌に対するBcl-2およびBcl-xLを選択的に阻害する分子標的薬である．この試験では，Navitoclax投与による予後に影響を及ぼすバイオマーカーを評価することにある．評価対象となるバイオマーカーは，CYFRA21-1，NSE，CTCである．このとき，BATTing法を利用して，これらの最適なカットオフ値(応答は奏効割合)を探索している．
Devanarayan, V., Cummins, D.J., Tanzer, L.R.,, et al.: Application of GAM and tree models for assessing the role of drug resistance proteins in leukemia chemotherapy.
Proceedings of the American Statistical Association Joint Statistical Meetings, 1999.
Rudin, C.M. et al.: Phase II study of single-agent Navitoclax (ABT-263) and biomarker correlates in patients with relapsed small cell lunge cancer, Clinical Cancer Research,
18(11), 3163-3169, 2012.

治療効果モデルの活用事例

Zhu, J. Gallego, B. (2020). Targeted estimation of heterogeneous treatment effect in observational survival analysis, Journal oBiomedical Informatics, 103474.

Zhu & Gallego (2020)は，サブグループ同定のためにT-learnerにより推定されたHTEに対して，BART法を当てはめることで変数重要度(VI)を推定している．

HTEの推定VIを計算

Subgroupの摘出

結局はサブグループ同定が目的であり，治療効果モデルは，共変量のスクリーニングに応用されるに過ぎない．

事例検討：膵体尾部猟奇疾患に対する切除術に対する無作為化比較第III相試験
Yamada et al. (2021)は，膵体尾部領域疾患に対する切除術として，膵静脈剥離個別切除(Active：個別切除群)と膵静脈同地切断(Control：同時切除群)の多機関共同無作為化比較第III相試験を実施している．
出血量の比較

ここでは，因果Boosting法での結果について解釈する．

2500
2000
1500
1000
0

500

Amount of bleeding

共変量：年齢，性別，BMI，PS(Performance Status)，
糖尿病の有無，ステロイド服用の有無，抗凝固薬
服薬の有無，心血管障害の有無，放射線療法の有無，化学療法の有無，膵切除予定部位の厚み，主膵管経，膵切除ラインの位置，血液・生化学検査(WBC，Hb，PLT，T-Bil，Cre，Alb，Ami，HbA1c)

3000

アウトカム：出血量

Control

ActiveYamada, S. et al. (2021). Safety of combined division vs separate divisionof the splenic vein in patients undergoing distal pancreatectomy:A noninferiority randomized clinical
trial, JAMA Surgery, 156(5), 418–428.

lower risk

higher risk

Amount of bleeding

500

500

1000

2500
2000
1000

1500

Amount of bleeding

500
0higher risk

lower risk

0

0

-500

Estimated heterogeneous treatment effect

1500

3000

推定されたHTEの推移

Control

Active
tx

0

50

100

150
pts

200

250

300

Control

Activetx

推 定 HTE が 0 よ り も 大 き な グ ル ー プ (higher risk) で は ， 個 別 切 除 群(Active)のほうが出血量が多いのに対して，推定HTEが0以下のグループ(lower risk group)では，同時切断群(Control)のほうが出血量が多い．

化学療法の
有無

20

放射線療法
の有無

ステロイド
服用の有無

15

抗凝固薬服
薬の有無

40

10

心血管障害
の有無

性別

20

5

PS

尿病の有無

T.Bil

年齢

Alb

WBC

PLT

主膵管経

膵切除ライ
ンの位置

Hb

Ami

BMI

HbA1c

Cre

膵切除予定
部位の厚み

0

-200

-100

同時切除risk

60

0

100

個別切除risk

100

Partial Dependence

80

200

• 変数重要度より，膵切除予定部位の厚みが切除術による出血量の違いに影響を及ぼすことがわかった．
• 部分従属度より，水切除予定部位の厚みが10mm未満の場合には，個別切除での出血量が多く，20mm以上では，同時切除での出血量が多いことが示唆された．

部分従属度

25
30

膵切除予定部位の厚み

変数重要度

HTEの臨床での適用場面の一つが探索的サブぐr－
ぷ解析であることに留意すると，Virtual TWINS法(Foster, et al., 2011)のアルゴリズムを用いることができるのではないか？

治療効果モデルの多くがアンサンブル学習法を用いるため，モデルがブラックボックス化される．
そのため，結果の解釈が変数重要度や部分従属度
で「示唆」を与えるだけになる．

higher risk
0.765 subgroup

1

Node 4 (n = 22)10Node 5 (n = 37)

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

1500

Amount of bleeding

1000

800

0

Control

Active

0

1

0

1

1

0

Node 2 (n = 245)

0

0

0.765

600

Cre
Cre

400

3

200

10.35

Amount of bleeding

10.35

1000

膵切除予定部位の厚み
Atsumi

500

1

2000

2500

1200

(個別切除群のほうが出血リスクが高いリスク集合を考える)

3000

1400

推定HTEが正値の場合に1，負値の場合に0としたときの分類樹木の結果

0

Control

ActiveFoster, J.C., Taylor, J.M.G., Ruberg, S.J.: Subgroup identification fromrandomized clinical trial data, Statistics in Medicine, 30:2867–2880,
2011.

治療効果モデルに対する適応上の留意点：再現性の問題
Raghavan et al.(2022)は，2型糖尿病に対する血糖コントロールの有効性に関する2つの無作為化比較試験(ACCORD
study, VADT Study)の主要エンドポイント(全死因死亡率)に対して因果Forest法を適用した．
因果Forest法の2試験に対する変数重要度

2試験のKendallの順位相関係数は0.717であり，決して高い数値ではないことから，同じ試験デザインであっても，HTEを再現することは困難であることを指摘している．
(限定的な)解決策として，Raghavan et al.(2022)は，複数の試験データを併合することで，より再現性の高い結果をもたらすことができると結論付けている．

Raghavan, S. et al. (2022) Generalizability of heterogeneous treatment effects based on causal forests applied to two randomized clinical trials ofintensive glycemic control, Ann. Epidemiol., 65, 101-108.

おわりに• 本講演では，異質性治療効果(HTE)に対する統計モデルとして，治療効果モデルを紹介した．
• 治療効果モデルは，統計学・機械学習において多数の提案がなされているが，ここでは，(1) 樹木に基づくモデル(2) Meta-learner(3) 変換アウトカム/修正共変量法(4) 共有基底条件付き平均回帰に分けて紹介した．
• 治療効果モデルの多くがアンサンブル学習法を用いることからモデルがブラックボックス化される．このことは，治療効果モデルを臨床応用するうえにおいて，多大な欠点となり得ることを指摘した．

• HTEを用いる一つの適応場面が探索的なサブグループ解析であることから，Virtual TWINSの流儀に倣い，事後的にCART樹木を用いて解釈する方法を紹介した．
• 結局のところ，HTEが推定できる統計モデルを今後，どのように活用するのかについては，今後の課題である．

Thank you for your kind attention
toshibow2000@gmail.com

shimokaw@wakayama-med.ac.jp

Supplementary

付録１：多変量適応型回帰スプライン法

多変量適応的回帰スプライン法 (MARS; Multivariate Adaptive Regression Splines)
Friedman (1991)◼ MARS法では，ステップ関数による基底関数を打ち切りベキ乗基底関数に変更することで，CART法の問題を解決している．
◼ MARS法では，データの層別と回帰分析を同時に遂行することが可能であり，かつ非線形構造や交互作用も捉えることができる．
◼ MARS法のモデル構築アルゴリズムは，前進ステップワイズ過程，後退ステップワイズ過程，モデル系列からの最適モデルの選択によって構成される．

c p ( j ,m )■MARS法のモデル：
M

Jm

m =1

j =1

f MARS ( x ) =   m  H ( s( j ,m ) ( x p ( j ,m ) − c p ( j ,m ) ))x p ( j ,m )[−( x p ( j ,m ) − c p ( j ,m ) )]+[+( x p ( j ,m ) − c p ( j ,m ) )]+打ち切りべき乗関数(括弧内が負値ならば0,正値ならばそのまま)

[ x p ( j ,m ) − c p ( j ,m ) ]+H ( s( j ,m ) ( x p ( j ,m ) − c p ( j ,m ) )) = [c p ( j ,m ) − x p ( j ,m ) ]+, s( j , m ) = 1, s( j ,m ) = −1

Friedman, J.H.: Multivariate adaptive regression splines, The Annals of Statistics, 19(1), 1-67, 1991.

1.0

MARS法のイメージ

CARTモデル

0.0

MARSモデル-0.5

・ふし点の数=ステップの数なので，複雑な樹木になっている（樹木が複雑になると解釈が困難）．
・一般に予測確度に劣る（関数の頂点付近では当てはまっていない）．

MARSモデル

-1.0

y

0.5

CARTモデル

真のモデル

0.0

0.5

x

1.0

・区分点の数+1=基底関数の数なので，
CARTに比べて簡素な結果．
1.5・非線形構造，交互作用項構造をCARTに比べて高い予測確度で捉えることができる．

前進ステップワイズ法

前進

後退

ﾓﾃﾞﾙ
選択

予め設定された基底関数の最大数Mmax ，および交互作用の最大次数Jmaxまで，過度に大きなモデルを生成する．
このとき，以下の点に留意する：(1) CARTと異なり，分割(交互作用項の追加)を行った場合でも，親ふし(追加前の基底関数)を削除し
ない．
CART法
xp

MARS法
初期分岐にxpが選定された
場合，その後は，xpとの交互作用のみが探索される．

複数の主効果が評価できる．

(2) 基底関数のなかに同一の変数を含めない．
(3) [ x p ( j ,m ) − c p ( j ,m ) ]+ , [c p ( j ,m ) − x p ( j ,m ) ]+ を各ステップで一度に追加する．
■ パラメータ推定(最小2乗法)

N( M .Jm )β = arg min β ( fˆMARS( x) − y)2n =1■ 基底関数の追加 → 残差平方和が最小になるように基底関数を追加する．

前進ステップワイズ法のイメージ

前進

後退

̂ 0
c p (1,1) J =1

x p (1,1)

̂1

̂ 3
c p ( 3, 4 ) J =3x p ( 2, 2 )̂ 7

x p ( 3, 4 )x p (1,3)

̂ 5

̂ 2
c p ( 2, 2 ) J =2c p (1,3) c p (1,1) c p (1,3) ̂ 6c p ( 2, 2 ) ̂ 4
c p ( 3, 4 ) ̂ 8[Step.0]yˆ = ˆ0[Step.1]yˆ = [Step.0] + ˆ1[−( x p (1,1) − c p (1,1) )]+ + ˆ2 [+( x p (1,1) − c p (1,1) )]+[Step.2]yˆ = [Step.1] + ˆ3 [+ ( x p (1,1) − c p (1,1) )]+ [−( x p (2,2) − c p (2,2) )]+ + ˆ4 [+ ( x p (1,1) − c p (1,1) )]+ [+ ( x p (2,2) − c p (2,2) )]+[Step.3]yˆ = [Step.2] + ˆ5 [−( x p (1,3) − c p (1,1) )]+ + ˆ6 [+( x p (1,3) − c p (1,1) )]+[Step.4]yˆ = [Step.3]+ ˆ7 [+( x p (1,1) − c p (1,1) )]+ [−( x p ( 2, 2 ) − c p ( 2, 2 ) )]+ [−( x p (3, 4 ) − c p (3, 4 ) )]++ ˆ8 [+( x p (1,1) − c p (1,1) )]+ [−( x p ( 2, 2 ) − c p ( 2, 2 ) )]+ [+( x p (3, 4 ) − c p (3, 4 ) )]+

ﾓﾃﾞﾙ
選択

後退ステップワイズ法

前進

後退

ﾓﾃﾞﾙ
選択

前進ステップワイズ法で得られた，大きなモデルを後退すてっぷワイズ法により，切片項まで逐次に削除する．

これにより，最大モデルから切片のみ(Nullモデル)までのモデル系列が得られる．
(1) CARTの刈り込み過程と異なり，親ふし(基底関数)が削除されても，子ふしは削除しない．

(2) 前進ステップワイズ法では，2個の打ち切りベキ乗基底関数を追加したが，後退ステップワイズ法では，1個のみが削除される．

CART法

MARS法
例1:

CART法では分岐が刈り込まれる

例2:

MARS法では基底関数が刈り込まれる．

後退ステップワイズ法のイメージ

前進

̂ 0
̂1
(6)

̂ 3
(5)

̂ 7
(2)

̂ 6
(1)

(4)
̂ 4

(7)
̂ 8

ﾓﾃﾞﾙ
選択

注：アニメーションで見ること
̂ 5
(8)

̂ 2
(3)

後退

(a)

削除された基底関数(括弧内は順番)yˆ = ˆ0+ ˆ1 [−( x p (1,1) − c p (1,1) )]+ (6)+ ˆ 2 [+( x p (1,1) − c p (1,1) )]+ (3)+ ˆ3 [+( x p (1,1) − c p (1,1) )]+ [−( x p ( 2, 2 ) − c p ( 2, 2 ) )]+ (5)+ ˆ 4 [+( x p (1,1) − c p (1,1) )]+ [+( x p ( 2, 2 ) − c p ( 2, 2 ) )]+ (4)+ ˆ5 [−( x p (1,3) − c p (1,1) )]+ (8)+ ˆ6 [+( x p (1,3) − c p (1,1) )]+ (1)+ ˆ7 [+( x p (1,1) − c p (1,1) )]+ [−( x p ( 2, 2 ) − c p ( 2, 2 ) )]+ [−( x p (3, 4 ) − c p (3, 4 ) )]+ (2)+ ˆ8 [+( x p (1,1) − c p (1,1) )]+ [−( x p ( 2, 2 ) − c p ( 2, 2 ) )]+ [+( x p (3, 4 ) − c p (3, 4 ) )]+ (7)最適モデルの選定

前進

後退

ﾓﾃﾞﾙ
選択

■ 予測確度に焦点を当てたモデル選択 (OriginalのMARSにおいて採用)
一般化交差確認法(GCV: Generalized Cross Validation)の修正版(GCV基準) (y − yˆ )Nn(M ) 2nN

GCV ( M ) = n =1
1 − M ( + 1) / N 2

M ：モデルの番号yˆ n( M )：M 番目のモデルに当てはめたときの予測値

：基底関数の最適化をはかるパラメータ(大きいほど小さなモデルが選択される(経験的に2    4 )■ 情報量規準に基づくモデル選択 (Survival MARS (LeBlanc & Crowley, 1999)において採用)
赤池の情報量規準(AIC) 等 (重回帰のモデル選択と同様)

AIC ( M ) = −2l ( M ) + 2 m

M ：モデルの番号，
m ：基底関数の個数l ( M ) ：M番目のモデルの対数尤度

：罰則項(通常はα=3)■ その他のアプローチ
・縮小推定の利用 (後退ステップワイズ法以降の代替に縮小推定を使用)：Taylan & Weber (2018)
・Bagging-MARS法 (前進ステップワイズ法のみでBaggingによるアンサンブル)

LeBlanc, M. and Crowley, J.: Adaptive regression splines in the Cox model, Biometrics, 55, 204-213, 1999.
Taylor, P. and Weber, G.W.: CG-Lasso Estimator for Multivariate Adaptive Regression Spline, Mathematical Methods in Engineering, 24, 121-136, 2018.

GCV基準による最適モデルの選定のイメージ
M =8

GCV (8)M =7GCV (7)M =6

GCV (6)最適モデルとして選択

M =5

GCV (5)M =4

GCV (4)M =3

GCV (3)yˆ = ˆ0+ ˆ [−( x1p (1,1)− c p (1,1) )]++ ˆ2 [+( x p (1,1) − c p (1,1) )]++ ˆ3 [+( x p (1,1) − c p (1,1) )]+ [−( x p (2,2) − c p (2,2) )]++ ˆ4 [+( x p (1,1) − c p (1,1) )]+ [+( x p (2,2) − c p (2,2) )]++ ˆ5 [−( x p (1,3) − c p (1,1) )]++ ˆ6 [+( x p (1,3) − c p (1,1) )]++ ˆ7 [+( x p (1,1) − c p (1,1) )]+ [−( x p (2,2) − c p (2,2) )]+ [−( x p (3,4) − c p (3,4 ) )]+M =2

GCV (2)M =1

GCV (1)M =0

GCV (0)

すべてのMに対してGCV(M)を計算し，最小のGCVをもつモデルを最適モデルとして選択する．

+ ˆ8 [+( x p (1,1) − c p (1,1) )]+ [−( x p (2,2) − c p (2,2) )]+ [+( x p (3,4) − c p (3,4) )]+Supplementary

付録２：アンサンブル樹木

樹木に基づくアンサンブル学習法
基本学習器には，単独では予測確度に乏しいモデルが用いられる．
例) 線形モデル，樹木，スタンプ(2分岐)など
基本学習器1

基本学習器2

基本学習器3(base learner)(base learner)(base learner)

f1 ( x )f2 ( x)f3 ( x )・・・

基本学習器B(base learner)f B ( x)・・・

統合(Committee)
Bootstrap型(Bagging, RandomForest)

Boosting型(Adaboost, MART)B

1 B
Fbts ( x ) =  f b ( x )B b =1

Fbst ( x ) =  b f b ( x ) b：重みパラメータ

基本学習器はBootstrap標本により生成される(個々の基本学習器は独立)．

b番目の基本学習器はb-1個の基本学習器の予測確度を補完するように生成される(ステイジワイズ戦略)．

b =1

樹木に基づくアンサンブル学習法の基本的アルゴリズム(Friedman & Popescu, 2008)
N
Input：学習標本：{ yi xi }1
アンサンブル回数：B

sb ( ) ：η≦Nの部分標本
N

初期モデル： F0 ( x ) = arg min  L( yi ,  )i =1

For b=1 to B do
pb = arg min p  L{ yi , Fb −1 ( xi ) + f ( xi , p)}iSb ( )Fb ( x ) = Fb −1 + · fb ( x )ensumble

Output：ensumble= { fb ( x )}B1＋

pb：基本学習器(樹木)のパラメータ(RandomForestでは，樹木の分岐選択変数集合をランダム化に選択，刈り込みなし)

 ：結合パラメータ(MART：任意, Bagging/RandomForest：0)

fb ( x ) = f ( xi , pb )End

(MART：非復元抽出, Bagging/RandomForest：復元抽出)B

MART：Fbst ( x ) =  b f b ( x )
b =1

1 B
Bagging/RandomForest： Fbts ( x ) =  f b ( x )B b =1＋ ・・・すなわち，樹木に基づくアンサンブル学習法は，樹木モデルの加法形式でモデルが構築される．

Friedman, J.H.and Popescu, B.E.: Predictive learning via rule ensemble, Annals of Applied Statistics, 2(3), 916-954, 2008.

Boosting接近法に基づくアンサンブル学習法

樹木の当てはめ

樹木の当てはめTree.1

Tree.2モデルの更新モデルの更新モデル
Tree.1

モデル
Tree.1 + Tree.2

・・・

修正標本
樹木の当てはめ

当てはめ結果に基づく更新

修正標本

当てはめ結果に基づく更新

学習標本・・・・・・Tree.B
モデルの更新・・・・・・

最終モデル

反復更新によりモデルを精緻化：多重化法型回帰樹木法(Friedman, 2001)などが含まれる

勾配Boosting樹木のアルゴリズム
N

F
(
x
)
=
argminInitialize： 0   L ( yn ,  )  n =1

For m = 1 to M do(a) n = 1, 2, , N に対して以下を計算
 L( yn , f ( xn )) 
rnm = −  f ( xn )
 f = fm−1

疑似残差の計算(b)樹木モデルを疑似残差 rnm に当てはめて終結ふしの領域
R jm ( j = 1, 2, , J M ) を計算する．
(c)領域 R jm ( j = 1, 2, , J M ) に対して回帰パラメータを計算
する．





 xn R jm

 jm = argmin    L( yn , Fm −1 +  ) 

領域の計算

各領域における回帰パラメータの推定．このとき，これまでの回帰パラメータについては再推定していないことに注意(ステイジワイズ過程)(d)モデル更新
Jm

Fm ( x ) = Fm −1 ( x ) +   jm I ( x  R jm )End For

j =1

M

Output： Fˆm ( x) =  f m ( x)m =1

このアルゴリズムは，Mを過度に大きくすると過学習を引き起こす恐れがある．Friedman(2001)は，モデル更新(赤色の数式)に縮小パラメータν(経験的に0.1未満)を付与することで過学習を回避できることを指摘している(ニューラルネットワークにおける慣性項と同様の役割)．すなわち，
この更新式は以下に変更される．

Jm

Fm ( x ) = Fm −1 ( x ) +    jm I ( x  R jm )j =1

勾配Boosting樹木の例示

Iteration

インドの小児の栄養失調に関する人口保険調査データ(Fenske et al., 2009)

0

200

400

600

800

1000

00

200
200

400
400

600
600

アンサンブル回数B

アンサンブル回数B

Iteration

Iteration
Iteration

800
800

8.5

7.5

6.5

5.5

7.5
7.5
6.5
6.5

v=1.00 (縮小推定なし)
v=0.10v=0.01

5.5
5.5

損失関数

6.5

7.5

v=1.00 (縮小推定なし)
v=0.10v=0.01CV Mean Squared Error

樹木の深さ = 3の場合

8.5
8.5

樹木の深さ = 2 (stump)の場合

CV Mean Squared Error
CV Mean Squared Error

Sample=0.5,Depth=3.0
Depth=1.0
Sample=1.0,

8.5

Sample=1.0, Depth=1.0

5.5

損失関数

CV Mean Squared Error

応答：栄養不良スコア， 説明変数：小児のBMI,小児の月年齢，母親のBMI

1000
1000

8.5
.5

ared Error

8.5.5

ared Error◼ 縮小パラメータν=1.00 (縮小推定なし)の場合，アンサンブル回数の増加とともに損失関数が減少している(過学習を起こしている)．
Sample=1.0, Depth=3.0
Sample=0.5, Depth=3.0◼ 縮小パラメータν=0.01の場合，早期に樹木が上昇傾向を示している．つまり，縮小パラメータを小さくしす
ぎると，殆どアンサンブルが終了しないままになる．

Fenske, N., Kneib, T. and Hothorn, T.：Identifying risk factors for severe childhood malnutrition by boosting additive quantile regression, Technical Report 052, University of Munich.

確率勾配Boosting樹木のアルゴリズム
Friedman(1999)は，部分標本を用いることで，勾配Boostingに比べて計算負荷が小さく，かつ，予測確度に優れたモデルを構築できることを指摘している．この方法は，確率勾配Boosting樹木と呼ばれる．

Jm

Fm ( x ) = Fm −1 ( x ) +    jm I ( x  R jm )End For

0

j200
=1

400

600

M

Output： Fˆm ( x) =  f m ( x)m =1

Iteration

800

1000

8.5

200

400
600
800
アンサンブル回数B
Iteration

8.5
7.5
6.5

400

ほぼ同じだが，より早期に収束してIte
いる

Sample=1.0, Depth=3.0

0

200

1000

8.5

Sample=0

7.5

5.5(d)モデル更新

0

6.5



1000

5.5

 xn R jm

8.5 7.5
7.5 6.5
6.5 5.5 jm = argmin    L( yn , Fm −1 +  ) 

400
600
800
アンサンブル回数B
Sample=0.5,
Depth=1.0
Iteration

η =0.5 (サブサンプル=1/2)

5.5

7.5



6.5

(c)領域 R jm ( j = 1, 2, , J M ) に対して回帰パラメータを計算
する．

Error
CV Mean Squared
Squared Error
CV Mean
損失関数

8.5

CV Mean Squared Error(c)樹木モデルを疑似残差 rn*mに当てはめて終結ふしの領域
Sample=1.0, Depth=1.0
R jm ( j = 1, 2, , J M )を計算する．

200

5.5

8.5
0

CV Mean Squared Error L( yn* , f ( xn* )) 
rn*m = − f ( xn* ) f = fm−1

CV Mean Squared Error(b)抽出された標本 n* = 1, 2, , N に対して以下を計算

v=1.00 (縮小推定なし)
v=0.10v=0.01

7.5

For m = 1 to M do(a)抽出割合  (0.0    1.0) に対して N = N  個の部分標本を非復元抽出する．

η =1.0 (サブサンプルなし)

6.5

 n =1

Initialize：F0 ( x ) = argmin   L( yn ,  ) Sample=05.5

Sample=1.0, Depth=1.0

CV Mean Squared Error
損失関数



N

Hastie et al. (2014)は経験的に1/2で十分であることを指摘
0
200
400
600
800
1000
0
200
400
している．

数値例での結果

8000
6000
4000
2000

テスト標本の残差平方和

10000

Boston住宅価格(Harrinson and Rubinfeld(1978), MASSパッケージ)データにおいて，テスト標本(N=200)をランダムに抽出し(非復元抽出)，残りのデータを学習標本に用いて，1,000回シミュレーションしたときの結果

Depth=1

Depth=3

Depth=5

縮小パラメータ ν=1.0
サブサンプル = 1.0 (なし)

Depth=1

Depth=3

Depth=5

縮小パラメータ ν=0.1
サブサンプル = 1.0 (なし)

Depth=1

Depth=3

Depth=5

縮小パラメータ ν=1.0
サブサンプル = 0.5

Depth=1

Depth=3

Depth=5

縮小パラメータ ν=0.5
サブサンプル = 0.5

Bootstrap接近法に基づくアンサンブル学習法
学習標本
0

リサンプリング

Bootstrap標本1

Bootstrap標本2

Bootstrap標本3

樹木の当てはめ

樹木の当てはめ

樹木の当てはめ・・・・・・

Bootstrap標本B
樹木の当てはめ・・・・・・Tree.1

Tree.2

Tree.3

最終モデル

Tree.B

モデルの統合(Committee)

Bootstrap標本を用いて，個別に樹木を構成したうえで統合する．Bagging(Breiman, 1996)およびRandom
Forest(Breiman, 2001)がこれに該当する．

Bagging法のアルゴリズム (回帰の場合)

Breiman, L.: Bagging predictors, Machine Learning, 24, 123-140, 1996.

Require ブートストラップ回数 B
For b = 1 to B do(b )N(a)ブートストラップ標本を生成 ( Lboot = { yu ( n ) , xu ( n ) }n =1 , b = 1,  B )b)(b)ブートストラップ標本L(bootにCART樹木を当てはめる(Tb ( x))End Forb1
Ensemble： Fˆbag ( x ) =  Tb ( x )B b =1

分類の場合には，B個の樹木における多数決ルールに基づいてクラス予測値が決定される．

単一樹木より，なぜブートストラップ標本を用いた複数樹木のほうが予測確度に優れるのか？
いま，データ集合 L = { yn , xn } によるCART樹木をT ( x) ，L が従う確率分布をP ，Lの経験分布をPN とする．また，Lのブートストラップ標本を Lboot ，ブートストラップ分布を PN とするとき，分布収束 N (PN − PN ) ~(b )大きなN(N → ∞)に対して，

1 bb)
Tb ( x) = Eˆ T ( xn ; L(boot)  → E T ( xn ; L ) 
B b =1

N (PN − P ) より，十分(１)である．また，E L T ( xn ; L ) 2   E L T ( xn ; L )  であることから，
2

E L  yn − T ( xn ; L )    yn − E L T ( xn ; L ) 


2

2

(２)

b)
Eˆ T ( xn ; L(boot)  で予測すれば，単一樹木でのとなる．式(1)，式(2)より，ブートストラップ法を用いて多数の樹木を生成し，残差平方和に比べて小さくなることがわかる．詳しい証明はBreiman(1996)を参照．

数値例での結果

6000
5000
4000
3000
2000

テスト標本の残差平方和

7000

Boston住宅価格(Harrinson and Rubinfeld(1978), MASSパッケージ)データにおいて，テスト標本(N=200)
をランダムに抽出し(非復元抽出)，残りのデータを学習標本に用いて，1,000回シミュレーションしたときのCART樹木とBagging樹木の結果

CART

CART

B=10

B=20

B=50

B=100

B=200

Bagging樹木

Harrison, D. and Rubinfeld, D.L.：Hedonic prices and the demand for clean air. J. Environ. Economics and Management 5, 81–102, 1978.

RandomForest法
Bagging法とRandomForest法の違い
Bagging法
Bootstrap標本

・・・

CART樹木の当てはめ

RandomForest法

・・・

Bootstrap標本

Bootstrap標本

CART樹木の当てはめ

CART樹木の当てはめ(分岐過程のみ)

Bootstrap標本
CART樹木の当てはめ(分岐過程のみ)・・・

・・・
分岐に用いる変数についてもランダム化

分岐に用いる変数についてもランダム化

RandomForestでは，・各樹木(基本学習器)の生成は，ブートストラップ標本が用いられる (Baggingと同じ)．
・CARTの分岐過程のみが行われる(過剰適合するほど大きな樹木を構築する)・各樹木の分岐選択では，P 個の変数のなかからランダムに選択された P ( P  P ) が用いられる．

RandomForest法のアルゴリズム (回帰の場合)
Require ブートストラップ回数 B
For b = 1 to B do(b )N(a)ブートストラップ標本を生成 ( Lboot = { yu ( n ) , xu ( n ) }n =1 , b = 1,  B )b)(b)ブートストラップ標本L(bootにCART樹木を当てはめる(Tb ( x))．このとき，各樹木の構成は以下のルールに基づく．

・CART法の分岐過程のみを行い，終結ふしの最小個体数 nminまで分岐させる．
・各CART樹木の各分岐では，P 個の説明変数のなかからランダムに選択された P 個を用いる．
End Forb1
Ensemble： Fˆbag ( x ) =  Tb ( x )B b =1

分類の場合には，Bagging法と同様に多数決ルールに基づいて選定される．

なお，上記のアルゴリズムは，Breiman(2001)では，RandomForest-RIと呼ばれている．その他に，RandomForest-RCという方法がある．この方法では，P 個の変数の線形結合に基づいてカットオフ値が選定される．

P の選定において，Hastie et al.(2014)は，回帰問題の場合には，P = P ，分類問題の場合には P = P / 3 を推奨している．
過剰適合するほど樹木を成長させること，並びに分岐変数選択にランダム性を導入することで予測確度が向上することの証明は，Breiman(2001)，杉本他(2005)を参照されたい．

杉本知之・下川敏雄・後藤昌司：樹木構造接近法と最新の発展，計算機統計学，123-164，2005．
Breiman, L.: Random Forests, Machine Learning, 45(1), 5-32, 2001.

数値例での結果

8000
6000
4000
2000

テスト標本の残差平方和

Boston住宅価格(Harrinson and Rubinfeld(1978), MASSパッケージ)データにおいて，テスト標本(N=200)をランダムに抽出し(非復元抽出)，残りのデータを学習標本に用いて，1,000回シミュレーション
したときのCART樹木とBagging，RandomForestの結果

CART

CART

B=10

B=20

B=50

B=100

Bagging

B=200

B=500

B=10

B=20

B=50

B=100

RandomForest

B=200

B=500

樹木に対する縮小推定
アンサンブル・モデル： F ( x ) =BB T ( x) =      I ( x  t )
b =1

b

b =1

t

tT

(b)

b番目の樹木の終結ふしの集合 T(b)

Boosting樹木： =  , t =ステイジワイズ過程において推定された回帰パラメータ1
Bagging, RandomForest樹木：  = ,  t = y (t ) (CART法と同様に推定された終結ふしの予測値)
B
青色の部分のみに注目すると，アンサンブル・モデルは，B個の変数の線形モデルと見做すことができる．
そのため，アンサンブル学習を行った後に，アンサンブル樹木を再評価し，不要な樹木を削除する(あるいは各樹木に重み付けを行う)ことが考えられる．
事後重み付けを伴うアンサンブル・モデル： F ( x ) =B  T ( x)
b =1

b b

最も直截的なアプローチは，lassoなどの縮小推定の方法を用いてパラメータ  bを推定することである．すなわ
ち，
2
N
B
B
である．

 


ˆ
{b } = argmin b   yn −  bTb ( xn )  +   | b |
b =1
b =1

 n =1 前進ステイジワイズ線形回帰(forward stagewise linear regression)
一方で，Hasite et al.(2007)は，前進ステイジワイズ線形回帰が，係数軌道が単調に増加するlasso(monotone
lasso)と一致することを示している．
Initialize ˆ b = 0 (b = 1, 2, , B)
For m = 1 to M do
2
N
B




(a){ * , b*} = argmin  ,b   yn −   ˆ k Tk ( xn )  − Tb ( xn ) *(b)ˆ b*  ˆ b* + ・  sign( )n =1End For k =1bOutput： FˆFS ( x) =  ˆ bTb ( x)b =1

前立腺がんデータに対する線形回帰に対するlassoと前進ス
テイジワイズ線形回帰の回帰パラメータの推移を表している(Hastie et al.(2009)より引用)
軌道が非常に似ていることがわかる(前進ステイジワイズ
過程においてステップ関数なのは，εによって上昇するためである)．
Hastie, T. et al. (2009) The Element of Statistical Learning (2nd edition),
Springer [杉山将他 (2014). 統計的機械学習の基礎，共立出版].

損失関数(残差平方和)に対して，最も影響のある樹木およびパラメータを推定する．

(a)において選ばれた変数をεだけ増加あるいは減少させる．

縮小推定の有無による予測確度の違い

テスト標本(N=200)をランダムに抽出し(非復元抽出)，残りのデー
タを学習標本に用いて，200回シミュレーションしたときの結果

Boston住宅価格データ(説明変数 = 13個)
7000
6000
5000
3000

4000

テスト標本の残差平方和

6e+06
5e+06

2000

4e+06

テスト標本の残差平方和

7e+06

インドの小児栄養データ(説明変数 = 4個 )LinearCART

Bagging

RandomForest

Boosted Trees

Stagewise Trees

LinearCART

Bagging

RandomForest

Boosted Trees

Stagewise TreesSupplementary
付録３：Okasa(2022)のシミュレーション結果

Model 1(Null,
等標本サイズ)Model 2
(無相関,
等標本サイズ)Model 3(一定,
不等標本サイズ)

Model 4(ステップ関数,
不等標本サイズ)Model 5(線形関数, 不等標本サイズ)

Supplementary
付録４：Powers et al.(2018)のシミュレーション結果

シミュレーションの結果：π(x)=1/2 (RCTの場合) 1/2
Null：Null model，TO：Transformed outcome forest，DB：different-bassis forest，CF：Causal forest
PTO0：PTO forest ，CB0：因果Boosting，BCM0(因果Boosting)

主効果：線形
治療効果：なし

主効果：線形，治療効果：指示関数

主効果：非線形治療効果：線形

主効果：非線形治療効果：非線形

シミュレーションの結果：π(x)=1/2 (RCTの場合) 2/2
Null：Null model，TO：Transformed outcome forest，DB：different-bassis forest，CF：Causal forest
PTO0：PTO forest ，CB0：因果Boosting，BCM0：因果MARS

主効果：なし治療効果：指示＋交互作用

主効果：線形
治療効果：線形

主効果：指示関数
治療効果：非線形

主効果：非線形治療効果：指示＋交互作用シミュレーションの結果：π(x)=exp(μ(x)-τ(x)/2)/(1+ μ(x)-τ(x)/2) (観察研究の場合) 1/2
Null：Null model，CF：Causal forest，PTO0：PTO forest ，CB0：因果Boosting(PS調整なし)，
CB1：因果Boosting(PS調整あり)， BCM0：因果MARS(PS調整なし)，BCM1(PS調整あり)
主効果：線形
治療効果：指示関数

主効果：線形
治療効果：なし

主効果：非線形治療効果：線形

主効果：非線形治療効果：非線形

シミュレーションの結果：π(x)=exp(μ(x)-τ(x)/2)/(1+ μ(x)-τ(x)/2) (観察研究の場合) 2/2
Null：Null model，CF：Causal forest，PTO0：PTO forest ，CB0：因果Boosting(PS調整なし)，
CB1：因果Boosting(PS調整あり)， BCM0：因果MARS(PS調整なし)，BCM1(PS調整あり)
主効果：なし治療効果：指示＋交互作用

主効果：線形
治療効果：線形

主効果：非線形治療効果：指示＋交互作用

主効果：指示関数
治療効果：非線形Supplementary
付録5：Virtual TWINS法大局アウトカムモデリング

Virtual Twins法 Foster et al.(2011)

無作為化比較試験では，新治療あるいは既存治療のいずれかに割り付けられて，治療が行われる，新治療が評価既存(標準)治療が評価

新治療群

データが存在データなし

既存(標準)治療群

データなし

データが存在

つまり，共変量Xiをもつi番目の被験者の治療効果(新治療と既存治療のアウトカムの差)は，いずれかのデータがないため，評価できない．

新治療が評価
新治療群

既存(標準)治療群

既存(標準)治療が評価
治療効果の推定値

治療指標を加えてモデル化： f ( x , z )( x ) = f ( x, z = 1) − f ( x, z = 2)

さらに，Virtual Twins法では，サブグループ抽出のために  ( x ) を応答にしたもとでCART法を用いている．
Foster, J.C., Taylor, J.M.G., Ruberg, S.J.: Subgroup identification from randomized clinical trial data, Statistics in Medicine, 30:2867–2880, 2011.

Virtual Twins法
いま，N人の被験者に対して，治療指標zn(1：薬A，2：薬B)に対して説明変数xn，応答ynが与えられているとする．

[STEP.1(Global Modeling過程)] ：治療指標をダミー変数として加える
RandomForest(RF)による被験者毎の治療効果の推定
新治療が評価

新治療群
既存(標準)治療群

既存(標準)治療が評価

RandomForestモデル fˆRF ( x , z )Foster et al.(2011)では，RandomForest法を用いているが，実際には，交互作用を含む非線形モデルであれば，何でもよい．

STEP.1によって推定されたモデルを用いて治療効果を推定する：ˆ n ( x ) = fˆRF ( x , z = 1) − fˆRF ( x, z = 2)

ˆ ( x ) = Pˆ − Pˆ ，Pˆ = Pr( y = 1| x , z = k )2値データの場合には，モデルにより推定される帰属確率を計算：
n
1n
2n
kn
n
n
n

[STEP.2(Subgroup Identification過程)] ：治療群のダミー変数は除外する
推定治療効果に対するCART樹木を当てはめ(アプローチによって2種類が存在)■ VT(R)：ˆ n を応答とする場合(回帰樹木の利用)．

■ VT(C)：ˆ n にしきい値を設定する場合(分類樹木の利用)．

Virtual twinsの例示
肺血症に関する無作為化比較試験を想定したシミュレート・データ(aVirtualTwinパッケージのテストデータ)• Outcome：生存の有無• 共変量：THERAPY(1：新治療，0：既存治療)，最初の治療がfailureしてからの期間(TIMFIRST) ，年齢(AGE)，ベースラインでの血小板数(BILLPLAT)，ベースラインSOFAスコア(心血管，血液，肝腎，呼吸器機能のスコア；blSOFA)，ベースラインでのクレアチニン(BILLCREAT)，ベースラインでの臓器不全の個数(ORGANUM)，輸血前のAPACHE IIスコア(PRACHE)，ベースライン
でのGlasgowスケール(BLGCS)，ベースラインでの血清IL-6濃度(BLIL6)，ベースラインでの活動性スコア(BLADL)，ベースラ
インでの局所ビリルビン(BLLILBILI)
0.057
100%
yes

PRAPACHE < 27 no

輸血前のAPACHE IIスコア-0.077
67%

0.32
33%

AGE < 62

AGE < 55

年齢

年齢

 APACHE IIは重症患者の重症度を測る指標として広く用いられている指標である．
 APACHE IIが27点未満(重症度が高くない)被験者に対して，新治療は有効であり，とくに，55歳以上の被験者においては，その傾向は顕著である．
メモ
Virtual twinsは，RのパッケージaVirtualTwinsで実行できる．ただし，
このパッケージには樹木図を描くための関数が用意されていないため，パッケージrpart.plotを入れる必要がある．記載例は次の通り，-0.009425%> data(sepsis)> vt.o <- vt.data(sepsis, "survival", "THERAPY", T)> vt.f <- vt.forest("one", vt.o)> vt.tr <- vt.tree("reg", vt.f, threshold = 0.01)> rpart.plot(vt.tr$tree)

BLGCS >= 14

Glasgowスケール

-0.12
41%

-0.051
14%

0.047
11%

0.13
10%

0.41
24%

新治療(active)が既存治療(Control)に比べて良好なサブグループ

Supplementary
付録6：Rのパッケージ

Rの関数一覧
パッケージ

場所

因果樹木法

causalTree

github(susanathey/causalTree)

因果Forest法

grf

CRAN

因果BoostingcausalLearning github(saberpowers/causalLearning)

PTO Forest

causalLearning github(saberpowers/causalLearning)

因果MARS

causalLearning github(saberpowers/causalLearning)

Python〇〇(github)Rの関数一覧
Meta-learnerは，pythonのライブラリ(causalml，EconML)に多数収録

causalml

S-learner, T-learner, X-learner,
R-learner, DR-learner

パッケージ
R-learner
DR-learner

S-learner, T-learner, X-learner,
EconML DA-learner, DR-learner

場所

Pythonrlearner

github(xnie/rlearner)〇(causalml)

drlearner

github(christophergandrud/drlearner)〇(causalml)(EconML)

Rでは，包括的なmeta-learnerのパッケージが存在しないようである．
ただし，自作してもさほど大変ではない．