Blogs

Blogs

経時的に変化する治療 (Time-varying treatments) に対する因果推論

0
By Topics | Analytics Students & Educators

はじめに

統計的因果推論に関する多くの書籍や文献では、その効果に興味のある特定の治療※１が1度だけ行われる場合の因果推論について紹介がされています。しかし、現実的には治療を複数回行い、それらの治療効果に興味がある場面も存在します。治療（介入）時点が複数ある例としては、新型コロナワクチンの接種がその最たるものかと思います。それ以外にもある集団（個人）に対して商品・サービスのレコメンド（e.g., DM、メール）を行い、その後の行動 (e.g., 購入) につながったかどうかを評価するなど、医学に限らず適用分野は多岐にわたります。

しかし、治療が複数回行われ各時点での治療レベルが異なりうる（時間の経過とともに治療が変化する）場合の因果推論は、その統計学的理論が比較的高度になってしまうこともあり、それほど認知されていないかと思います。そこで本コラムでは「治療が複数回行われた場合には因果推論を行うことがなぜ難しいのか」、「因果効果の推定はどのようにすればよいか」について可能な限り簡単に解説をし、SASでの実装方法にも触れていきたいと思います。

なお本コラムに関しては因果推論の基本的な理解があることを前提としています。必要がありましたら適宜関連書籍や因果推論に関する連載コラムを参考にしていただければ幸いです。

※１　介入 (intervention) 、曝露 (exposure) であっても同様の議論

 

治療の分類（"time-fixed" or "time-varying"）

治療時点が複数ある場合の因果推論を理解するために重要なことは、まず治療自体を十分に理解することです。今回は簡潔な議論のために、治療は2時点 (t = 0, 1) で行われる2値変数 （e.g., 治療を受ける or 受けない） であり、アウトカム（目的変数）は研究の最終時点でのみ測定されるものとします※2。また解析の対象となる研究集団はfixed study populationであるとし※3、研究開始以降の治療による因果効果のみを興味の対象とします（研究開始以前の治療による効果は議論しない）。

ここで複数回行われる治療は、各時点の治療の取り方から以下の2つのいずれかに分類がされます。

  • Time-fixed treatment
  • Time-varying treatment

まず前者からその分類の定義を説明します。

治療がtime-fixedであるとは、研究対象集団におけるすべての被験者において、最初に行われる治療（ベースライン治療）が研究を通して変わらず継続的に実施されることを意味します。異なる言い方をするのであれば、一番初めの治療レベルがその後全ての治療レベルを決定する（同一なものとする）ということになります。

Time-fixed treatmentsに対する因果推論は比較的容易です。なぜならば、一番初めの治療がその後すべての治療を完全に決定するため、このコラムの後半で説明する時間依存性交絡 (time-dependent confounding) の問題は生じず、治療が一度だけ行われる場合と同様に、研究開始時点の共変量（ベースライン共変量）の調整を行うという議論に帰着するためです。すなわち、その調整を行う手法（回帰、マッチング、傾向スコアを用いた手法など）を適用することが可能です※4。

これに対し、time-fixedではない治療すべて、つまり各時点の治療が一番初めの治療レベルと異なりうる治療のことをtime-varying treatmentsと呼びます。例えば、治療を月に1回行われる運動指導とし各時点で取りうる行動が「自宅で何もしない（指導を受けない）」もしくは「運動指導を受けにいく」であるとしましょう。そして各時点での治療の選択が個人の意思に委ねられているとすると、気が向けば指導を受けますし、そうでなければ指導受けないといった判断が発生します。すなわち、各時点での治療レベルが異なる可能性があるわけです。

ここで、上記の状況における治療とアウトカムを以下のような記号で表すとします※5。

  • At：時点tにおける二値治療 (t = 0, 1)
    • A0：t = 0（研究開始時点）での治療
    • A1：t = 1での治療
    • Āt = (A0, A1)：治療のパターン
  • Y：研究の最終時点において観測されるアウトカム（生活習慣病の発生の有無）

すると、被験者が受ける治療パターンĀt = (A0, A1)としては次の4つが考えられます。

  • (A0, A1) = (0, 0) ：ともに指導を受けない
  • (A0, A1) = (0, 1) ：t = 0では指導を受けず、t = 1では指導を受ける
  • (A0, A1) = (1, 0) ：t = 0では指導を受け、t = 1では指導を受けない
  • (A0, A1) = (1, 1) ：ともに指導を受ける

この治療のパターンはtreatment regime（治療レジメン）と呼ばれます※6。また、それぞれのレジメンに対応した潜在アウトカム（の期待値）が、一時点のみ治療が行われる場合と同様に想定されます。

  • E[Ya0=0, a1=0]：ともに指導を受けない場合の潜在アウトカムの期待値
  • E[Ya0=0, a1=1]：t = 0では指導を受けず、t = 1では指導を受ける場合の潜在アウトカムの期待値
  • E[Ya0=1, a1=0]：t = 0では指導を受け、t = 1では指導を受けない場合の潜在アウトカムの期待値
  • E[Ya0=1, a1=1]：ともに指導を受ける場合の潜在アウトカムの期待値

一般に、研究を行う際に我々の関心の対象となるのは上記の潜在アウトカムの期待値の比較によって定義がされる平均因果効果です。例えば、E[Ya0=1, a1=1] - E[Ya0=0, a1=0] という比較は、何もしなかった（二時点とも指導を受けなかった）と比べて二回とも指導を受けると、どの程度アウトカムに因果効果があったのかということを意味しています。つまり、time-varying treatmentsに対しての因果効果を推定するためには、ある治療レジメンの潜在アウトカムの期待値E[YA]を推定できれば良いことになります。

Time-varyingに対する因果推論が難しい第一のポイントは、単一の時点で治療が行われる場合と同様に十分な治療の定義を行うことに加え、治療がどのような治療の経路（治療レジメン）をとるのかということも考慮する必要があることです。なお治療レジメンは、次のセクションで紹介するようにいくつかの観点から分類することが可能です。

※2　時点数や治療の取りうるレベルが今回のケース（ともに2）よりも増えた場合や、アウトカムが最終時点だけでなく各時点での治療の前に繰り返し測定される場合、アウトカムがtime-to-event（生存時間）の場合にも今回の議論は適用可能
※3　各被験者の追跡開始日が明確に定義されたclosed cohort（研究開始以後に追加される被験者は存在しない集団）
※4　各推定手法が要求する仮定は満たされる必要がある
※5　大文字は確率変数、小文字はその実数値、太文字はベクトルを意味する
※6　regime以外にもstrategy, plan, protocol, policyなどとも呼ばれますが、既にデータが存在している場合の研究（観察研究）においては治療歴 (treatment history) と呼ばれることが多い

 

治療レジメンの分類（"dynamic" or "static" / "deterministic" or "stochastic"）

治療レジメンについて注目するにあたり、必要となる追加の記号を以下のように定義します。

  • Xt： 時点tにおける共変量ベクトル
    • e.g.) X0：ベースライン共変量
    • X‾t: ある時点tまでの共変量歴
  • Ht = (X‾t, At-1) ： 時点tにおける治療までの共変量歴および治療歴
    • e.g.) H0 = X0：ベースライン共変量

このように記号を置くと、各被験者が受ける治療のパターンを意味するtreatment regimeは以下の2つの観点からその分類を行うことが可能です。

  1. "static" or "dynamic"
    •  static regime（静的レジメン）または non-dynamic regime（非動的レジメン）
      • 任意の時点tでの治療atが時間のみに依存して決定される
    • dynamic regime（動的レジメン）
      • 任意の時点tでの治療atが時間と他の変数htに依存して決定される
  2. "deterministic" or "stochastic"
    • deterministic regime（決定論的レジメン）
      • 任意の時点tにおいて、いずれかの治療を受ける確率が100％（決定的に定まる）
    • stochastic (random)  regime（確率論的レジメン）
      • 任意の時点tにおいて、いずれかの治療を受ける確率が変動する（確率的に定まる）

すなわち、2×2の4通りにtreatment regimeは分類されます。この説明だけですと分かりにくいので、それぞれの場合の具体例を紹介しておきます。

  • static かつ deterministic な regime
    • 治療開始から3ヶ月経ったら治療を行う
  • static かつ stochastic な regime
    • 治療開始から3ヶ月までは50%の確率で治療群に割り当て、それ以降は80%の確率で治療群に割り当てる
  • dynamic かつ deterministic な regime
    • 血圧が150を上回ったら治療を行う（治療時点に加えて、血圧という共変量によって治療の実施を決定）
  • dynamic かつ stochastic な regime
    • 血圧が150を上回ったら80%の確率で治療群に割り当てる

 

因果効果の推定に必要な条件

Time-varying treatmentsに対して因果推論を行う際には、こちらのコラムでも紹介している識別可能条件 (identifiablity assumptions) と呼ばれる3つの条件を、以下のように複数時点に拡張する必要があります※7。それぞれが意図することは、治療が単一の時点でのみ行われる場合と大きく変わりはありませんので、もし概念的に詳しく知りたい方は上記のコラムもご覧いただけますと幸いです。今回は重要なポイントのみの説明です。

  • Consistency（一致性）
    • If Ā = ā for a given subject, Yā = Y for that subject
      • 治療レジメンāを取る場合の潜在アウトカム、すなわちāで表される一連の治療を受ける場合に観測されるであろう結果（ifの結果）は、実際にその治療レジメンを受けた場合に観測される結果と一致する
  • Sequential exchangeability（逐次交換可能性）
    • Yā  ∐ At | ht for all regime ā and all X‾t
      • 全ての時点において、各治療レベル間の交換可能性が保たれており群間比較が妥当な状況である
      • 各時点での治療の割り当てがランダムに行われる研究デザイン （逐次ランダム化実験、Sequential randomized experiments）において期待されるが、それ以外の研究デザインでは交絡の調整に十分な共変量が測定されていることと同義
  • Positivity（正値性）
    • P( 0 < P(At =1 | Ht = ht) < 1) = 1 for all ā
      • 各時点において特定の治療を受ける確率が0 or 1になるような被験者は存在しない。
      • 全ての時点において、任意の治療を絶対に受ける（もしくは受けない）といった状況は存在しない
※7　厳密には治療レジメンのタイプによって置かれる仮定がやや異なり、今回のコラムではstatic regimeを想定した

 

時間依存性交絡の問題

逐次ランダム化実験以外の研究デザインではランダム化が行われていないため、任意の時点までにある一連の治療を行った集団において実際に治療を受けた集団と受けない集団の共変量のバランスは一般的には異なります（交絡が発生）。すなわち各時点での治療群と対照群の交換可能性が保たれていないため、結果として各治療レジメンをとる集団は比較可能ではありません。これは例えば、被験者の予後が共変量であるとすると、t = 0, 1でともに治療を受けなかった (a0=0, a1=0) 群と受けた (a0=1, a1=1) 群でそれを比較すると下図のようにその分布がアンバランスになっているという状況です。

つまり、交絡を十分に制御（共変量の調整）をしなければ、各集団間のアウトカムの差は、治療による因果効果によるものか、それとも集団の特性によるものであるかが判断できません。Time-varying treatmentsに対する因果推論が難しい最大の理由は、こういった時間に依存する交絡、時間依存性交絡 (time-varying confounding) の問題です※8。

時間依存性交絡を理解するために以下のようなDAGで表される変数間の関係を考えてみます（DAGの構造は正しいとする）。

ここでは冒頭の例と同じようにA0とA1をそれぞれt=0とt=1での運動指導の有無、そしてX1を体脂肪率、Yを生活習慣病の発生の有無とします。また、各時点での運動指導の有無はアウトカムである生活習慣病に対して影響を与え（A0→Y, A1→Yの存在）、初めに受けた運動指導はその後の体脂肪率にも影響を与えている（A0→X1の存在）ものとします。さらに体脂肪率は次の時点で運動指導を受けるかどうかという意思決定にも影響があり（X1→A1の存在）、加えて最終的な生活習慣病にも因果的な効果をもたらしている（X1→Yの存在）とします。

ここで、t=0, 1での治療の有無によるアウトカムへの因果効果に興味があるとします。すなわち、以下の2つの因果効果が推定したい対象になります。

  1. A0からYへの因果効果 (直接効果：A0→Y、間接効果：A0→X1→Y, A0→X1→A1→Y)
  2. A1からYへの因果効果 (直接効果：A1→Y)

しかし変数間の関係がこのDAGで表されるとき因果効果の推定にあたりある問題が発生します。それはX1の存在です。X1に着目して再度上記のDAGを見てみると、次の2つを満たしていることがわかるかと思います。

  • A1とYの交絡因子 (A1←X1→Y)※9
  • A0とYの中間因子 (A0→X1→A1→Y)

まずA1とYの交絡因子ですので、A1の治療効果をバイアスなく推定するためにはバックドアパスと呼ばれるバイアス経路 (A1←X1→Y) を閉じることが必要であり、X1を条件付ける必要があります。しかし、同時にX1はA0とYの中間因子ですので、条件付けをしてしまうとA0からYへの因果効果のうちX1を介す間接効果の推定結果にバイアスが含まれてしまいます。つまりX1を条件付けても付けなくても因果効果推定にバイアスが含まれてしまうという二項対立の状況が発生してしまっているわけです。よって、条件付けに基づく推定方法である回帰分析やマッチングなどの標準的な手法では、時間依存性交絡が存在する場合に推定結果にバイアスが含まれてしまいます※10。また、上記の2つの要件を満たす共変量X1のことを時間依存性交絡因子 (time-varying confounders / time-dependent confounders) ないしは時間依存性共変量 (time-varying covariates / time-dependent covariates) と呼びます。

※8　単一の時点でのみ治療が行われる場合と同様に、交絡という事象ありその調整に十分な変数として交絡因子が定義される
※9　もしくはA0と共通の原因を持つということが時間依存性共変量である条件になる
※10　No feed back (X1→A1がない)という仮定（ある時点の治療は以前の共変量の影響を受けない）があれば可能だが非現実的

 

Time-varying treatmentに対する因果効果の推定手法

上記で説明したようにtime-varying treatmentsに対しては、ベースライン共変量の調整を行う一般的な手法では推定結果にバイアスが含まれてしまいます。そこで、この問題に対応するものとしてハーバード大学のJames M. Robins教授により一般化 (generalized) されたのが、g-methodsと総称される以下の3つの手法です (Robins' g-methods) 。

  1. G-computaion algorithm formula ("g-formula")
  2. Inverse probability of treatment weighting (IPTW) of marginal structural models (MSMs)
  3. G-estimation of structural nested models (SNMs)

これらは、線形回帰のような手法よりもより弱い制約条件の下で興味のある因果効果を推定することが可能です。一部簡単にそれぞれの特徴をまとめたものが下表です。

気になることとしては一体どの方法が良いのかということかと思いますが、残念ながら明確な答えはありません。実際に解析を行う場合には、上記でまとめたような各手法のメリットやデメリット、仮定への頑健性などから総合的にどの手法を選択すべきかの判断をしていく必要があるかと思います。またg-methods以外の推定手法についても近年研究が進んでおり、いくつか学術論文も出版されています※11。

※11 Edward H. Kennedyによる提案など

 

SASでの実装方法 (IPTW of MSMs)

今回はg-methodsの中で最も理解のしやすいIPTW of MSMsについて、そのSASでの実装方法を紹介していきます。まず行うことはmarginal structural model（周辺構造モデル）の特定になります。これは潜在アウトカムの周辺期待値E[Ya0, a1]と、治療計画āの関係を示したものになります。今回は周辺構造モデルをE[Ya0, a1]  = β0 + ψ0a0 + ψ1a1として考えます。つまり、各時点での治療効果は加法的に（足し算の形で）潜在アウトカムの期待値に因果効果を及ぼしており、治療間や治療と共変量との交互作用はないものとします。

余談ですが、なぜ”周辺”という言葉を用いるかといいますと、個人レベルの潜在アウトカムYa0, a1を期待値E[Ya0, a1]という集団ベルにまで周辺化 (marginalize) させるためであり、これは我々が興味のある平均因果効果 (average causal effect) を、しばしば周辺因果効果 (marginal causal effect) とも呼ぶ理由と同じです。

さて、上記の周辺構造モデルを想定したときt = 0, 1でともに治療を受けた場合に、ともに治療を受けなかった場合と比べてどの程度因果効果があるのかは、下記のように式展開を行っていくと2つの因果効果に分解することが可能であり、ψ0, ψ1という2つの因果パラメータを推定できれば因果効果をその組み合わせによって興味のある因果効果を算出することができます。

  • E[Ya0=1, a1=1] - E[Ya0=0, a1=0]  = (β0 + ψ0 + ψ1) - β0  = ψ0 + ψ1
    • ψ0: t = 0（ベースライン）での治療による因果効果
    • ψ1: t = 1での治療による因果効果

そして、このパラメータ推定を行う際に利用されるのがIPTW法です。IPTW法を併用する場合には、まず各個人のアウトカムをそれぞれが受けた実際の治療 Ā を受ける確率の逆数で重み付けし、その重み付けした後の集団（疑似母集団、pseudo population）に対して想定した周辺構造モデルを当てはめ、通常の回帰分析と同様にパラメータ推定を行います。 これによってなぜ興味のある因果パラメータを推定することができるのかについては、理論的な話になってしまうため今回は省略しますが簡単に理由を述べますと、傾向スコアが適切に推定されそれを用いて作成された疑似母集団における条件付きアウトカムの期待値は無条件での潜在アウトカムの期待値と一致するためです (E[Ya0, a1] = Eps[Y | A0=a0, A1=a1) 。

 

ここからが具体的なSASでの実装方法に移ります。まず以下のDAGで表される想定する変数間の構造を考えます。

ここで各記号は以下を示しているものとします。

    • X0：t = 0（ベースライン）での共変量
    • X1：t = 1での共変量（時間依存性共変量）
    • A0：ベースラインでの治療
    • A1：t = 1での治療（２回目の治療）
    • Y：観測されるアウトカム

すなわち、ある時点の治療はそれ以前の共変量と治療の影響を受けており、かつアウトカムに対しても因果効果を持つものとし、これは共変量に関しても同様です。この構造となり、ψ0 = -10, ψ1 = -4となるようなデータを生成します（この値を算出することが目標）。

/* データのシミュレーション */
data df;
	call streaminit(12345); /* シード値の指定 */
	keep ID Gender Age Press TG;
	N=5000; /* サンプルサイズの指定 */
	do i=1 to N;
/* 変数の生成 */
		ID = i;
		Gender = rand('BINORMINAL',0.5,1);
		Age = round(20+45*rand('UNIFORM'),1);
			if Gender=0 then Press = round(rand('NORMAL',100,20)+0.6*Age,1);
			if Gender=1 then Press = round(rand('NORMAL',110,20)+0.6*Age,1);
		TG = round(rand('NORMAL',90,15) +0.6*Age + 10*Gender,1);
		output;
	end;
run;
/* t=0での治療変数の追加 */
data df;
	set df;
	call streaminit(12345); /* シード値の指定 */
	keep ID Trt0 Gender Age Press TG;
	PS0_True=logistic(0.5+0.8*Gender-0.03*Age);/* t=0での真の傾向スコアモデル */
	Trt0 = rand('BINOMINAL',PS0_True,1);
run;
/* t=1での治療変数、アウトカムの追加 */
data df;
	set df;
	call streaminit(12345); /* シード値の指定 */
	keep ID Trt0 Trt1 Outcome Gender Age Press TG;
　	PS1_True=logistic(0.6-0.4*Gender-0.01*Age+0.5*Trt0+0.01*Press-0.015*TG);/* t=0での真の傾向スコアモデル */
	Trt1 = rand('BINOMINAL',PS1_True,1);
　	epsilon = rand('NORMAL',0,10);
	Outcome = round(130 - 10*Trt0 - 4*Trt1 +Gender +Age +Press+ TG+ epsilon,1);
run;
/* 変数の順序変更 */
data df;
	format ID Trt0 Trt1 Outcome Gender Age Press TG;
	set df;
run;

 

因果パラメータの推定にあたって、我々が行うステップは大きく3つです。

  1. 各時点での傾向スコアの算出
    • 今回推定しなければいけない傾向スコアは2つ※12
      • t = 0での治療を受ける確率
      • t = 1での治療を受ける確率
    • それぞれロジスティック回帰モデルによって推定を行う
  2. 各個人への重みの計算
    • 通常の治療レジメンを受ける確率の逆数であるweightではなく、stabilized weightを採用するため各時点の集団においてある治療を受ける割合も算出
  3. 因果パラメータの推定
    • 周辺構造モデルを重み付けした後の集団に適用し、モデルを当てはめ (GENMOD Procedure)

これらのプログラミング・コードが以下になります。

/*各時点の傾向スコアの算出*/
/*t=0*/
proc logistic data=df descending;
	class Gender(ref="0") /param=ref ref=first;
	model Trt0=Gender Age;
	output out=df predicted=PS0;
run;
/*t=1*/
proc logistic data=df descending;
	class Gender(ref="0")/param=ref ref=first;
	model Trt1=Gender Age Trt0 Press TG;
	output out=df predicted=PS1;
run;
/*SWの分子の計算*/
proc logistic data=df descending;
	　model Trt0=;
　	output out=df predicted=Pr0;
run;
proc logistic data=df descending;
　	class Trt0/param=ref ref=first;
	model Trt1=Trt0;
　	output out=df predicted=Pr1;
run;
/*各個人のweightの計算*/
data df_wt;
	set df;
　	if Trt0=1 then wt0=Pr0/PS0; else wt0=(1-Pr0)/(1-PS0);
　	if Trt1=1 then wt1=Pr1/PS1; else wt1=(1-Pr1)/(1-PS1);
	wt=wt0*wt1;
run;
/*因果パラメータの推定*/
proc genmod data=df_wt;
	class ID Trt0 Trt1/param=ref ref=first;
　	model Outcome = Trt0 Trt1/ error=normal link=id;
　	weight wt;
	repeated subject=ID/ type=ind;
run;

上記を実行する次のような結果となり、赤枠で囲った部分からも分かるように当初にシミュレーションした因果パラメータをバイアスなく推定できていることが分かります。

※12 傾向スコアモデルが正しく特定される必要がある

 

補足

Time-varying treatmentに対する因果推論について今回紹介を行いましたが、簡潔な説明のために下記の要素も仮定（条件）しています。

  • 共変量と治療は一定の等間隔で収集される
    • e.g.) 週ごとや月ごとに測定されるデータ
  • 打ち切りや脱落、変数の測定誤差は存在しない（完全なデータである）
  • アウトカムは研究最終時点においてのみ測定される単一の確率変数
  • 未測定交絡は存在しない（No unmeasured confounders）

これらの仮定が満たされない場合には、今回紹介した各手法を対応するよう拡張する必要があります。この点についても興味がある方がいらっしゃいましたら、本コラムにコメントいただけますと幸いです。また今回紹介したtime-varying treatmentsに対する因果推論の概要や、その手法の1つである周辺構造モデルにおけるIPTW法については今年度のSASユーザー総会2022において口頭発表を行っています。以下がその際の講演資料を修正したものになりますのでご自由にご活用ください。

Tags JPN Academic 統計的因果推論
Share

About Author


Academic Staff

公共ソリューショングループ / 東京医科大学大学院医学研究科 生物統計学（Biostatistics）と統計的因果推論（Causal Inference）を専攻しています。お気軽にご連絡ください。

Leave A Reply

Back to Top