G formula
今ひとつNested Structual Modelが理解できていない気がして、Stata Journalから以下をまとめる。
http://fmwww.bc.edu/repec/bocode/g/gformula_Daniel_DeStavola_Cousens.pdf
要約
この記事では、gformulaコマンドについて説明する。このコマンドはG-formula手法の実装であり、時間変動する曝露が時間変動する交絡因子（曝露自体の影響を受ける交絡因子を含む）を伴うアウトカムに与える因果効果を推定するために使用される。この手法は、曝露が中間変数を介してアウトカムに因果効果を及ぼす場合に、制御された直接効果および自然直接/間接効果を推定する関連問題にも対応できる。特に、中間変数とアウトカムの関係を交絡する因子が曝露の影響を受ける場合に有効である。理論の簡単な概要、コマンドおよびそのオプションの説明、そして2つのシミュレーション例を用いた実例を示す。
キーワード: gformula, 因果推論, G-計算式, 時間変動する交絡, 媒介, 直接効果と間接効果
導入
時間変化する交絡因子
設定
縦断研究では、複数の時点でデータが収集されることが一般的であり、疫学、臨床試験、生態学、社会学、計量経済学など、多くの研究分野で広く行われている。特に本記事では、関心のある説明変数（または複数の変数）が時間とともに変化し、複数の固定された時点で各ユニット（または被験者）について測定される状況を扱う。このような状況において、以下のいずれかに対する時間変動する説明変数の因果効果に関心が向けられる。すなわち、(1) 研究終了時に測定される興味のあるアウトカム、または (2) フォローアップ期間中の任意の時点で発生する可能性がある興味のある事象の発生時間。ただし、この事象は離散的な時間（すなわち各訪問時）で測定され、前回の訪問以降に事象が発生したかどうかが評価される。
この因果効果を測定する際には、交絡因子の役割を考慮することが重要である。交絡因子とは、説明変数とアウトカムの両方に影響を与える変数のことであり（非公式な定義であり、詳細は後述）、これを考慮しない場合、通常は因果効果の推定値にバイアスが生じる。
交絡に関する一般的なテーマについては、数多くの文献で議論されている（Pearl, 2009; Rothman et al., 2008; Morgan and Winship, 2007; Angrist and Pischke, 2009）。本記事では、時間変動する交絡因子の特定の問題に焦点を当てる。すなわち、時間変動する説明変数とアウトカムとの因果関係を交絡する可能性があり、時間とともに変化し、研究期間中に繰り返し測定される要因である。特に、時間変動する交絡因子が関心のある時間変動する説明変数自体の影響を受ける場合、交絡を処理するための従来の方法（すなわち回帰調整）はもはや適用できない（Robins, 1986; Robins and Hernán, 2009）。
この状況は、因果ダイアグラム（Greenland et al., 1999）を用いると最もよく説明できる。図1に示されたように、矢印は因果的影響の仮定された方向を表している。\(A_0–A_T\)は、時間点\(0, 1, \cdots, T\)で測定される関心のある説明変数を表す。\(L_0–L_T\)は、時間点\(0, 1, \cdots, T\)で測定される潜在的な交絡因子を表し、\(L_t\)は\(A_t\)の直前に発生するものと仮定される。この図では、\(Y\)は研究終了時（訪問\(T+1\)）に測定される興味のあるアウトカムであり、\(U\)は\(L_0–L_T\)および\(Y\)に影響を与える未測定の要因の集合である。
注意すべき点は、\(U\)から\(A_0–A_T\)への矢印や、\(A_0–A_T\)と\(Y\)の共通原因（例えば\(V\)）からの矢印がないことである。これは、\(\{L_0, L_1, \cdots, L_t\}\)および\(\{A_0, A_1, \cdots, A_{t-1}\}\)を条件とした場合、\(A_t\)が\(Y\)に因果効果を持たない場合には、\(A_t\)が\(Y\)と独立であるという仮定を表している。この仮定は「未測定の交絡因子が存在しない」という仮定として知られており、各訪問時点\(t\)で、\(A_t\)と\(Y\)の関係を交絡する要因の十分なセットが測定されていることを意味する。
また、\(A_t\)から\(L_{t+1}\)への矢印も注意が必要である。これらの矢印は、ある時間点での交絡因子が、直前の時間点での説明変数の影響を受ける可能性があることを表している（なお、図を読みやすくするために、\(A_0\)から\(L_2\)への矢印や\(L_0\)から\(A_1\)への矢印などは省略されているが、\(U\)から\(A_0–A_T\)への矢印の省略は重要である）。
一方、アウトカムが時間-事象データ（タイム・トゥ・イベント）であり、離散的な時間点\(1, 2, \cdots, T+1\)で測定される場合、適切な因果ダイアグラムは図2のようになる。この場合、各\(Y_t\)は、その時間間隔\((t-1, t]\)で事象が発生したかどうかを示す二値変数である（ここでも、\(A_0\)から\(L_2\)、\(L_0\)から\(A_1\)、\(L_0\)から\(Y_2\)への矢印などを含めることも可能であるが、\(U\)から\(A_0–A_T\)への矢印は含められていない）。
従来のモデルの限界
標準的な方法では、\(L\)による交絡を調整するために、回帰分析で\(L_0\)から\(L_t\)までを条件付けする。もし\(A_t\)から\(L_{t+1}\)への矢印がなければ、この戦略は\(A\)から\(Y\)へのいわゆるバックドアパス（Greenlandら、1999）をすべて遮断することに成功し、\(A_0\)から\(A_t\)までの\(Y\)（またはイベント発生までの時間アウトカムの場合は\(Y_1, \cdots, Y_{t+1}\)）への共同因果効果を一貫して推定することができる。しかし、図1と図2に示されるように、交絡因子が過去の説明変数の値に影響される状況では、\(L_0\)から\(L_t\)までを条件付けすることは\(2\)つの理由で有効ではない。例えば、\(A_0\)から\(Y\)への因果効果を考えてみる（図1）。\(L_0\)を条件付けすることで、\(A_0 \leftarrow L_0 \leftarrow U \rightarrow Y\)というバックドア（非因果）パスを遮断することに成功した。しかし、\(L_1\)（および将来の\(L_t\)）を条件付けすることで、\(A_0\rightarrow L_1 \rightarrow L_2 \rightarrow \cdots \rightarrow Y\)（および他の多くのパス）を遮断してしまうが、これは\(A_0\)から\(Y\)への因果効果の一部を表している。さらに、\(U\)と\(A_0\)がともに\(L_1\)に影響を与えるため、\(L_1\)を条件付けすることで\(A_0\)と\(U\)の間に非因果的な関連が生じ（Greenlandら、1999を参照）、これにより\(A_0\)から\(U\)を通って\(Y\)への新たなバックドアパスが開かれてしまう。これは「コライダー条件付けバイアス」と呼ばれ、Hernánら（2004）によって分類された選択バイアスの一形態である。この問題は図2にも同様に適用される。
例1
HIV研究における抗レトロウイルス療法（ART）の縦断的研究では、\(A_t\)は時点\(t\)で被験者がARTを処方されているかどうかを示す二値変数、\(L_t\)は時点\(t\)で測定されたCD4カウント、\(Y_t\)は区間\((t-1, t]\)で被験者がAIDSを発症するかどうかを示す二値変数である。観察研究では、特定の時点でARTを治療するかどうかの決定は、患者の現在のCD4カウントに影響されると予想される。また、ARTは患者のCD4カウントを上昇させることで効果を発揮するため、前述の理由から標準的な回帰分析でCD4を調整することは賢明ではない。この例については、4.1節で再度取り上げる。
媒介
設定
内容的には異なるものの、方法論的には密接に関連した問題として、曝露\(X\)がアウトカム\(Y\)に与える因果効果を、媒介変数\(M\)を通じて作用する間接効果と、\(M\)によって媒介されない直接効果に分解したい場合に生じるものがある。
図3：曝露の影響を受ける媒介変数-アウトカム交絡因子を伴う媒介を示す因果ダイアグラム。
従来の手法の限界
この場合の標準的なアプローチは、(i) \(X\)（および\(X – Y\)間の交絡因子\(C\)）を条件とする\(Y\)の回帰モデルを適合させ、次に(ii) 同じモデルに\(M\)を追加することである。モデル(i)と(ii)の間で\(X\)の係数がどのように変化するかを調べることで、\(X\)が\(Y\)に与える効果のどの程度が\(M\)によって媒介されているかを解釈することがある。より形式的には、モデル(i)における\(X\)の係数は、\(X\)が\(Y\)に与える総効果を表し、モデル(ii)における係数はしばしば\(M\)によって媒介されない\(X\)から\(Y\)への直接効果を表すと解釈される。
このアプローチは、図3に示されるように、\(M – Y\)関係の交絡因子\(L\)が存在する場合には無効である。なぜなら、第2のモデルは\(X\)から\(Y\)への直接効果を一貫して推定することができないからである。\(M\)を条件付けることで、\(X\)と\(L\)の間に関連が生じ、\(X\)から\(L\)を経由して\(Y\)に至るバックドアパスが開いてしまう。このバックドアパスを遮断するには\(L\)を条件付けする必要がある。しかし、図3に示されるように\(L\)が\(X\)の影響を受ける場合、\(L\)を条件付けることで\(X \rightarrow L \rightarrow Y\)というパスも遮断される。このパスは\(M\)によって媒介されない\(X\)から\(Y\)への直接効果の一部であるため、\(L\)を条件付けても\(M\)を条件付けたことで生じる問題を解決することはできない。
さらに、標準的な回帰アプローチでは、\(M\)が\(Y\)に与える効果が\(X\)によって修飾されないことが要求される。
時間変化する交絡因子との関係
2つの設定の関連性を理解するために、図3が\(T = 1\)の場合の図1と同じであることに注目する。このとき、\(L_0 = C, A_0 = X, A_1 = M, L_1 = L\)となる。したがって、時間依存の交絡の例では、\(A_0\)から\(A_T\)までの\(Y\)に対する因果効果は\(T+1\)の直接効果から構成される。つまり、\(A_0\)から\(Y\)への\(A_1\)から\(A_T\)までによって媒介されない直接効果、\(A_1\)から\(Y\)への\(A_2\)から\(A_T\)までによって媒介されない直接効果、そしてその後も同様に続く。
直接効果と間接効果についての詳細
直接効果と間接効果が何を意味するのかを正確に議論するために、反事実記法を用いる（Robins and Greenland, 1992; Pearl, 2001）。\(Y(x, m)\)を、反事実的である可能性があるが、介入によって\(X\)が\(x\)に設定され、\(M\)が\(m\)に設定された場合の潜在アウトカムとする。制御された直接効果\(\text{CDE}_m\)は、\(m\)を固定したまま、\(x\)の異なる値に対する\(\text{E}\{Y(x, m)\}\)の比較である。例えば、\(X\)が単変量で二値の場合、制御された直接効果 (\(m\)での) を次のように具体的に表現できる。$$\text{CDE}_m = \text{E}\{Y(1, m)\} – \text{E}\{Y(0, m)\}$$
次に、\(M(x)\)を、反事実的である可能性があるが、\(X\)が\(x\)に設定された場合の媒介変数の潜在値とする。総因果効果\(\text{TCE}\)は、\(x\)の異なる値に対する\(\text{E}[Y\{x, M(x)\}]\)の比較である。同様に、\(X\)が二値の場合、次のように表される。$$\text{TCE} = \text{E}[Y\{1, M(1)\}] – \text{E}[Y\{0, M(0)\}]$$
これらの量を使用して、総効果と直接効果の差として間接効果を推論することが望ましい。しかし、制御された直接効果が\(m\)の関数であるため、これが難しくなる。\(\text{CDE}_m\)は\(m\)の値ごとに異なる可能性がある。
この理由から、自然直接効果 (\(\text{NDE}_{x_0}\) は、\(x_0\)を固定したまま、\(x\)の異なる値に対する\(\text{E[Y\{x, M(x_0)\}]}\)の比較として定義される（通常、自然な選択肢が存在する場合、\(X\)の「基準値」に設定される）。言い換えれば、基準介入下で\(M\)がその自然な値を取る場合における\(X\)が\(Y\)に与える効果である。二値の\(X\)において、次のように表される。$$\text{NDE}_0 = \text{E}[Y\{1, M(0)\}] – E[Y\{0, M(0)\}]$$
次に、自然間接効果 (\(\text{NIE}_{x_1}\)は、総因果効果と自然直接効果の差として定義される。したがって、これは\(x_1\)を固定したまま、\(x\)の異なる値に対する\(\text{E}[Y\{x_1, M(x)\}]\)の比較である。これを二値の\(X\)に適用すると、自然間接効果は次のようになる。$$\text{NIE}_1 = \text{E}[Y\{1, M(1)\}] – \text{E}[Y\{1, M(0)\}]$$
自然直接効果と自然間接効果が適切に定義されているかどうかについては議論がある（Robins, 2003; Didelez, 2006; Hafeman, 2009）。これらの（興味深く重要な）哲学的な議論は本記事の範囲を超えていまうが、我々はこれらの効果が取り上げる状況では適切に定義されているという前提で進める。
例2
アルコール摂取が収縮期血圧（SBP）に因果効果を持つことは広く信じられているが、この因果効果がどのようなメカニズムを通じて作用するかは十分に理解されていない。一つの仮説として、アルコール摂取が肝酵素であるGGTのレベルに影響を与え、その結果、GGTがSBPに影響を与えるというものがある。したがって、アルコール摂取がSBPに及ぼす因果効果のうち、どの程度がGGTによって媒介されているかを知ることは興味深い。BMIはGGTとSBPの両方に影響を与えると考えられており、社会経済的地位（SEP）はアルコール摂取、BMI、SBPに影響を与えると考えられている。さらに、アルコール摂取はBMIに因果効果を持つことが知られている。この状況は図3で表されており、\(X = \) アルコール摂取、\(M = \text{GGT}\)、\(Y = \)収縮期血圧、\(L = \text{BMI}\)、\(C = \text{SEP}\)となる。この例については4.2節で再度取り上げる。
今後の方向性
前述の2つの状況での交絡を扱うには、標準的な回帰調整に代わる方法が必要である。その一つの方法として、g-computationがある。この方法はRobins（1986）によって最初に提案され、その後Robins and Hernán（2009）やTaubman et al.（2009）によってさらに議論された。次節では、この手法の概要を簡単に説明し、その後、第3節でStataの新しいコマンド（gformula）を使用した本研究の実装について述べる。第4節では、上述の2つの例を用いた具体的な例を示し、最後に第5節でまとめの考察を行う。
g-computationの手順
時間変化する交絡因子
基本的な考え方
g-computationは、まず観察データに見られる変数間の関係をモデル化することから始まる。これらのモデルを使用して、観察データのように自然に進展した場合ではなく、変数 \(A_0\)から\(A_T\)が介入によって決定された場合に、研究対象者に何が起こったかをシミュレーションする。このモデル化とシミュレーションは時間に沿って「順方向」に行われる。具体的には、まず時点\(0\)のデータを基に時点\(1\)のデータをモデル化する。これにより、時点\(0\)の曝露に対するさまざまな仮定された介入下での時点\(1\)のデータをシミュレーションして比較することが可能になる。次に、時点\(0\)と時点\(1\)のデータを基に時点\(2\)のデータをモデル化し、時点\(0\)および時点\(1\)の曝露に対するさまざまな介入下での時点\(2\)のデータをシミュレーションする。このようにしてシミュレーションを繰り返す。各介入下で、基準時以降の交絡因子やアウトカムもすべてシミュレーションされる。その後、異なる介入下でのアウトカムを比較することで、これらが無作為化実験から生成されたかのように因果推論を行うことが可能になる。
モデルの当てはめ
まず、\(L_0\)と\(A_0\)を条件とした\(L_1\)の条件付き分布に対してパラメトリックモデルを指定する。（時間固定の交絡因子がある場合は、これらを\(L_0\)に含める。）\(L_1\)が連続変数の場合、このモデルから得られる\(f_{L_1| L_0, A_0}(l_1| l_0, a_0;\alpha_1)\)は確率密度関数となる。\(L_1\)が二値変数の場合、これは条件付き確率を表す。
このモデルを\(L_1, L_0, A_0\)の観察データに適合させることで、パラメータ\({\alpha_1}\)の推定値\(\hat{\alpha_1}\)を得る。同様に各\(t \in [2, T]\)に対して、 \(L_0, A_0, \cdots, L_{t-1}, A_{t-1}\)を条件とした\(L_t\)の条件付き分布のモデルを指定し、観察データから\(f_{L_t| L_0, A_0, \cdots, L_{t-1}, A_{t-1}}(l_1 | l_0, a_0, \cdots, l_{t-1}, a_{t-1};\alpha_t)\)のパラメータ\(\alpha_t\)の推定値\(\hat{\alpha_t}\)を得る。
追跡終了時に測定される1つのアウトカム\(Y\)の場合、\(L_0, A_0, \cdots, L_T, A_T\)を条件とした\(Y\)の条件付き分布のモデルを指定し、観察データから\(f_{Y|L_0, A_0, \cdots, L_T, A_T}(y | l_0, a_0, \cdots, l_T, a_T; \beta)\)のパラメータ\(\beta\)の推定値\(\hat{\beta}\)を得る。
アウトカムがイベント発生までの時間、つまり\(Y_1, Y_2, \cdots, Y_{T+1}\)という一連の二値変数で記述される場合は、各\(t \in [1, T+1]\)に対して、\(L_0, A_0, \cdots, L_{t-1}, A_{t-1}\)を条件とし、 \(Y_{t-1} = 0\)を条件とした\(Y_t = 1\)の条件付き確率のモデルを指定する。そして、観察データのうち\(Y_{t-1} = 0\)（リスク集合にまだ属するもの）のサブセットから、\(f_{Y_t | L_0, A_0, \cdots, L_{t-1}, A_{t-1}}(y_t | l_0, a_0, \cdots, l_{t-1}, a_{t-1}; \beta_t)\)のパラメータ\(\beta_t\)の推定値\(\hat{\beta_t}\)を得る。
現在、gformulaではこれらのモデルの適合に使用できるのは、regress（線形回帰）とlogit（ロジスティック回帰）のみである。これらのモデルを各時点ごとに個別に適合させるか、すべての時点のデータを統合してパラメータを推定するかを選択できる。ただし、個別に適合させる場合でも、各時点でモデルの形式は同じである必要がある。詳細は第3節で説明する。
仮定された介入下でのシミュレーション：追跡終了時に測定された単一のアウトカムの場合
治療がすべての被験者において常に行われなかった場合、すなわち介入\(A_0 = 0, A_1 = 0, \cdots, A_T = 0\)の下で、研究対象者に何が起こったかをシミュレーションすると仮定する。
この介入下での\(L_0, L_1, \cdots, L_T\)の（シミュレーションされた）値をそれぞれ\(L_0*, L_1*, \cdots, L_T*\)と表す。\(L_0\)は\(A_0\)に先行するため、\(A_0\)から\(A_T\)に関する介入の影響を受けない。このため、\(L_0* = L_0\)となる。
\(L_1*\)は、分布 \(f_{L_1|L_1, A_0}(l_1 | L_0*, 0;\hat{\alpha_1})\)に基づいてシミュレーションされる（モデルの適合に関する前節を参照）。具体的には、観察データから推定された\(L_1\)の条件付き分布（\(L_0\)と\(A_0\)を条件とする）を使用し、\(A_0\)を\(0\)に、\(L_0\)を介入下での値\(L_0*\)に置き換えた後、この分布から\(L_1*\)をシミュレーションする。
• \(L_1\)が連続変数の場合、\(L_1*\)は密度\(f_{L_1| L_0, A_0}(l_1 | L_0*. 0;\hat{\alpha_1})\)に基づく確率的な抽出値である。
• \(L_1\)が二値変数の場合、\(L_1*\)は確率 \(f_{L_1|L_0, A_0}(1 | L_0*, 0; \hat{\alpha_1})\)を持つベルヌーイ分布からの確率的な抽出値である。
同様に、各\(t\in [2, T]\)に対して、\(L_t*\)は分布\(f_{L_0, A_0, \cdots, L_T, A_{T-1}}(l_1 | L_0*, 0, \cdots, L_{t-1}*, 0; \hat{\alpha_t})\)に基づいてシミュレーションされる。
最後に、\(Y*\)は分布\(f_{Y| L_0, A_0, \cdots, L_T, A_T}(y | L_0*, 0, \cdots, L_T*, 0; \hat{\beta})\)からシミュレーションされる。 \(Y*\)は仮定された介入下での潜在アウトカムとして知られている。これは、指定された介入下でアウトカムがどうなっていたかを表す。
このようにして、基準時以降のすべての変数（潜在アウトカムを含む）が、観測されていない交絡因子が存在しないという仮定およびモデル適合段階でのモデリング仮定の下で、介入「治療がすべての被験者において常に行われない」のシナリオに基づいてシミュレーションされる。注意すべき点は、シミュレーションの各段階で使用される条件付き密度が、曝露と交絡因子の過去の値にのみ条件付けされることである。この点が、標準的な回帰調整法とこのアプローチの本質的な違いである。後者では未来と過去の両方に条件付けるため、導入部で説明した問題が生じる。
なお、\(U\)は観測されていないため、シミュレーションは\(U\)の未観測分布にわたって周辺化されるが、\(U\)が\(A\)と\(Y\)の関係の交絡因子ではないため、これによりバイアスは生じない（Daniel et al., 2010）。
イベント発生までの時間アウトカムのシミュレーション
イベント発生までの時間をアウトカムとする場合、 \(Y_1*\)は確率\(f_{Y_1 | L_0, A_0}(1 | L_0*, 0; \hat{\beta_1})\)を持つベルヌーイ分布からシミュレーションされる。\(Y_1* = 0\)の対象者に対しては、\(Y_2 * \)が確率 \(f_{Y_2| L_0, A_0, L_1, A_1}(1 | L_0*, 0, L_1*, 0; \hat{\beta_2})\)を持つベルヌーイ分布からシミュレーションされる。同様に、最後に\(Y_T* = 0\)の対象者に対しては、\(Y_{T+1}*\)が確率 \(f_{Y_{T+1}| L_0, A_0, \cdots, L_T, A_T}(1| L_0*, 0, \cdots, L_T*, 0; \hat{\beta}_{T+1})\)を持つベルヌーイ分布からシミュレーションされる。これらの\(\{Y_1*, \cdots, Y_{T+1}*\}\)は、介入がすべての対象者において常に治療を行わないという仮定の下での潜在的なイベント発生までの時間アウトカムを表す。
複数の仮定された介入の比較
前述の「治療を行わない」という介入を「常に治療を行う」という介入に変更し、同様のシミュレーションを繰り返すことができる。この場合、\(A_0\)から\(A_T\)の各値を\(0\)ではなく\(1\)に置き換える。同様の原則を用いて、さらに多くの仮定された介入を比較することが可能である。たとえば、「偶数時点のみ治療を行う」や「時点\(3\)以降に治療を行う」などの介入をシミュレーションすることができる。
\(A_0\)から\(A_T\)が連続変数（または多変量）の場合、\(A_0\)を\(a_0\)、\(A_T\)を\(a_T\)に設定するようなさまざまな値の組み合わせを用いて異なる仮定された介入を比較することができる。
追跡終了時に測定された単一のアウトカムの場合、各介入について全被験者の平均的な潜在アウトカムを計算することで仮定された介入を比較できる。この平均は全被験者について計算されるため、すべての背景変数について周辺化されたものとなる。この意味で、g-computationは時間依存の曝露に対して有効な標準化の一形式と見なすべきである。
イベント発生までの時間アウトカムの場合、異なる仮定された介入下での平均発生率や累積発生率を比較し、それぞれの介入についてKaplan-Meier曲線をプロットすることができる。同様に、これらはすべての被験者を異なる介入下で比較することに基づいているため、他のすべての変数について周辺化されている。
未測定の交絡因子が存在しないという仮定およびモデル適合段階で用いられたパラメトリックモデリングの仮定の下で、平均的な潜在アウトカム、平均発生率、累積発生率、またはKaplan-Meier曲線の差異（有限サンプル誤差やモンテカルロシミュレーション誤差を超えるもの）は、曝露の因果効果に帰することができる。
シミュレーションされる被験者数
計算を高速化するため、シミュレーションされる被験者を元のデータセットのランダムなサブセットにすることができる。選ばれたサブセットが完全にランダムであれば、バイアスは生じないが、モンテカルロシミュレーション誤差が増加し、精度が低下する。この方法は、元のデータセットが非常に大きく、計算時間が受け入れ難いほど長い場合を除いて推奨されない。
ブートストラップを用いた標準誤差と推論
標準誤差と信頼区間はブートストラップを用いて得られる。ブートストラップサンプルは元のデータセットから被験者単位で抽出される。モンテカルロシミュレーションの回数が元のサンプルサイズより少ない場合、モンテカルロ用のサブセットはブートストラップサンプルから選ばれる。
動的介入の比較
これまでに考慮した介入はすべて「静的」と呼ばれるもので、（これらの介入が実施される仮想的な世界では）治療の経路が研究開始時点から完全に既知であると仮定している。これらの介入下での研究参加者の仮定された行動は、（曝露が交絡因子の過去の値に依存する）観察データを用いて構築されているが、我々の目的はこの依存性がないデータをシミュレーションすることにある。
別の介入カテゴリとして「動的介入」と呼ばれるものがある。この場合、治療の経路が交絡因子の経路に事前に指定された方法で依存することが許される。HIV研究における動的介入の例として、「CD4カウントが\(x\)未満になった場合に治療を開始する」というものがある。g-computationは、前述の通りの手順で使用することで、異なる値の\(x\)の下で研究参加者に何が起こるかをシミュレーションできる。\(x\)の値をさまざまに試すことで、このクラスの中で最適な介入（たとえば、HIV研究では期待されるAIDSフリーの生存期間を最大化する介入）が探求できる。この場合、介入値\(A_t*\)は初めから固定されるのではなく、\(L_0*, A_0*, \cdots, A_{t-1}*, L_{t}*\)に依存する。たとえば、特定の被験者について\(L_0* > x, \cdots, L_{t-1}* > x\)であるが\(L_t* < x\)である場合、\(A_0* = A_1* = \cdots, A_{t-1}* = 0\)となり、 \(A_t*\)は\(1\)に設定される。これは「決定論的動的介入」の例である。この場合、過去の交絡因子および治療の値が与えられたとき、動的介入を定義するルールは曝露に対して確率\(1\)で値を割り当てる。確率論的動的介入の例については、以下の「観察された介入下でのシミュレーション」に関する節を参照すること。
マージナル構造モデル (MSM) のパラメータ推定
これまで、g-computationを、異なる仮定された介入下での\(Y\)（または時間依存アウトカムの場合は\(Y_1, \cdots, Y_{T+1}\)の分布をシミュレーションする方法として説明してきた。そして、実際にそれがg-computationの本質である。しかし、比較をより簡潔に要約する方法が必要になる場合がある。この目的のために、マージナル構造モデル (Marginal Structural Model, MSM) を使用することができる（Robins et al., 2000）。
マージナル構造モデル (MSM) は、潜在アウトカムの分布のある特性を仮定された介入変数の関数として表現する。たとえば、追跡終了時に測定されたアウトカムの場合、\(Y(a_0, a_1, \cdots, a_T)\)を、二値変数\(A_0\)を\(a_0\)、\(A_1\)を\(a_1\)などに設定した仮定された介入下での潜在アウトカムとした場合、MSMは次のように表現される可能性がある。$$E\{Y(a_0, a_1, \cdots, a_T)\} = \gamma + \phi \sum_{t = 0}^{T}{a_t}$$
このように、このMSMが前提とする仮定の下では、 \(A_0, \cdots, A_T\)が\(Y\)に与える因果効果を、すべての潜在的なアウトカム間の（多数の）ペア比較ではなく、1つのパラメータ（\(\phi\)）で要約することができる。多数の異なる仮定された介入の下でシミュレーションを行い、\(Y*\)を\(\displaystyle \sum_{t=0}^{T}{A_t*}\)に回帰させることで、各シミュレーションされた介入データセットを連結した合計データセットから、\(\gamma\)および\(\phi\)の推定値を得ることができる。現在、gformulaではregress（線形回帰）とlogit（ロジスティック回帰）を用いたMSMの適合のみがサポートされている。
時間依存アウトカムの場合、\(Y_1–Y_{T+1}\)のモデルのパラメータがすべての時点にわたるプールされたロジスティック回帰から推定される場合、MSMとして自然な選択はマージナル構造Coxモデル（D’Agostino et al., 1990）である。gformulaはこのようなMSMをstcoxを使用して適合させることをサポートしている。
標準誤差と信頼区間は再びブートストラップによって得られる。
MSMの定義によれば、MSMは必然的に静的介入の比較に使用されるため、動的介入を比較する際にはgformulaでこのオプションを指定することはできない。動的MSMは開発されている（Cain et al., 2010）が、現時点ではgformulaではサポートされていない。
追跡脱落への対処
本節で検討されているような縦断研究では、追跡終了前に一部の被験者が脱落する可能性が常にある。追跡脱落がランダムに発生する（Little and Rubin, 2002）、すなわち、脱落が未観測データと条件的に独立しており、観測データ（脱落前に観測されたデータ）にのみ依存するという仮定の下では、このような追跡脱落をg-computationに非常に簡単に組み込むことができる。脱落は潜在的な治療経路の1つと見なすことができ、たとえば「常に治療を受け、脱落しない」といった経路のシミュレーションが行われる。この場合、観測不能交絡因子が存在しないという仮定に、欠測がランダムであるという仮定が暗黙的に含まれる。脱落を明示的にモデル化する必要がないことは、MAR（Little and Rubin, 2002）の下での尤度分析における無視可能性（ignorability）の一例である。
死亡によるセンサリングへの対処
上記で述べたことの例外は、センサリングが死亡によるものである場合である。その被験者が特定の介入下で死亡する時点以降のデータを仮定された介入下でシミュレーションすることは、不自然であり（実際には潜在的に誤解を招く可能性がある）、適切ではない。このため、サバイバル（生存）は、前述のAIDSフリー生存期間と同様にシミュレーションすべき追加のアウトカムプロセスとして見なすことができる。まず、「この被験者は区間 \((t-1, t]\)の間に生存していたか？」という質問が設定され、その回答が「はい」とシミュレーションされた場合に限り、「この被験者は区間\((t-1, t]\)の間にAIDSを発症したか？」という2番目の質問が設定され、その回答がシミュレーションされる。
時間固定変数の欠測値とアウトカムおよび時間変動変数の断続的な欠測への単一確率的補完を用いた対処
縦断研究では、被験者が研究から離脱して戻らない「脱落」に加え、断続的（または非単調な）欠測がよく見られる。これは、特定の訪問または複数の訪問で欠測が生じた後、次の訪問で被験者が戻るようなケースである。加えて、時点\(t\)において一部の観測値が欠測しているが、他の観測値は観測されている場合や、ベースライン（時間固定）変数の一部が欠測している場合もある。
欠測がランダムである（MAR）という仮定の下では（非単調な欠測パターンに対するこの仮定はやや議論の余地がある（Robins and Gill, 1997））、このような非単調な欠測パターンには、多重補完法（MICE: Multiple Imputation using Chained Equations）を用いて対処できる（van Buuren et al., 1999）。van Buurenらが説明した方法では、Little and Rubin（2002）が述べた意味での「適切な補完」を複数回行う。「適切な補完」とは、欠測データが観測データに基づく分布（この分布のパラメータの最尤推定値を用いる場合は「不適切」とされます）からではなく、この分布のパラメータの事後分布からのベイジアンな抽出値に基づいて行われることを意味する。この事後分布からの抽出は、複数の補完ごとに独立して行われる。
複数の適切な補完を行うことで、Rubinのルール（Little and Rubin, 2002）を使用して、最終的なパラメータ推定値の標準誤差を推定できる。
gformulaでは標準誤差を解析的に推定するのではなく、ブートストラップを用いて推定するため、gformulaに含まれる補完方法は、MICEを用いた単一の確率的補完である。この方法は、van Buurenら（1999）によって説明されたものと同一ですあるが、各欠測値に対して1回だけ補完を行い、その補完は「不適切な」補完である。Rubinの分散推定量を使用しない場合、このアプローチが有効であることが示されている（Tsiatis, 2006, 第14章）。
現在、gformulaでサポートされている補完コマンドは regress（線形回帰）、logit（ロジスティック回帰）、および mlogit（多項ロジスティック回帰）のみである。
観察データのレジーム下でのシミュレーション
複数の仮定された介入下で研究参加者に何が起こったかをシミュレーションすることに加え、介入が行われなかった場合、すなわち、観察データと同じメカニズムで被験者が治療を選択した場合に何が起こったかをシミュレーションすることも可能である。このためには、治療と交絡因子の履歴を条件とした治療割り当てのパラメトリックモデルを指定する必要がある。このレジーム下での\(A_1*\)から\(A_T*\)の値は、前述の\(L_1*\)から\(L_T*\)の場合と同様にシミュレーションされる。
追跡脱落がほとんどない場合、実際の観察データと観察データのレジーム下でシミュレーションされたデータを比較することで、この手法の成功度をある程度確認できる。この2つが大きく異なる場合、少なくともいくつかのパラメトリックモデリング仮定や、未観測の交絡因子が存在しないという仮定が成立していないことを示唆する。しかし、実際の観察データと観察レジーム下でのシミュレーションデータがよく一致していても、それらの仮定が成立していることを保証するものではない。
観察レジームは動的レジーム（前述参照）の一例である。実際、各\(A_t*\)は分布から抽出されており（観察設定で見られる変動を模倣するため）、これは確率論的動的レジームの一例である。
観察レジームと特定の介入との比較は、研究対象の集団にその介入を実施した場合の影響を評価する際にしばしば関心を引く（Taubman et al., 2009）。
媒介
媒介の例におけるg-computationは、基本的には同様の方法で機能する。ただし、自然直接効果と間接効果を推定する場合、異なる仮定された介入下でのシミュレーションを組み合わせる必要がある。たとえば、\(X\)が二値の場合、\(M\)は\(X =1\)および\(X = 0\)の下でそれぞれシミュレーションされ、それにより\(M(1)\)および\(M(0)\)が得られる。自然直接効果を推定するために必要な\(Y\{1, M(0)\}\)をシミュレーションするには、\(X\)を\(1\)に設定すると同時に、\(M\)を介入\(X = 0\)の下でシミュレーションされた値（つまり\(M(0)\)）に設定する。
\(X\)が二値でない場合、または多変量の場合、総因果効果、制御された直接効果、自然直接効果/間接効果を計算するための自然な比較（たとえば、\(1\)対\(0\)）が存在しないことがある。この場合、1.2節の式は以下に置き換えられる。$$\begin{eqnarray}\text{CDE}_m = \text{E}\{Y(X, m)\} – \text{E}\{Y(0, m)\} \\ \text{TCE} = \text{E}\{Y(X, M(X)\} – \text{E}\{Y(0, M(0))\} \\ \text{NDE}_0 = \text{E}\{Y(X, M(0)\} – \text{E}\{Y(0, M(0))\} \\ \text{NIE}_X = \text{E}\{Y(X, M(X)\} – \text{E}\{Y(X, M(0))\}\end{eqnarray}$$
ここで、\(0\)は依然として\(X\)の「基準値」を表すが、観察データで自然に生じる\(X\)の分布と比較される。このような比較は、しばしば関心のある因果的な質問に対応する。
いずれかの変数に欠測値がある場合は、前述の単一確率的補完（連鎖方程式を使用）を使用して対処できる。
gformulaのコマンド
文法
gformula mainvarlist [ if ] [ in ], outcome(varname) commands(string) equations(string)
[ idvar(varname) tvar(varname) varyingcovariates(varlist) intvars(varlist) interventions(string)
dynamic eofu pooled death(varname) derived(varlist) derrules(string) fixedcovariates(varlist)
laggedvars(varlist) lagrules(string) msm(string) mediation exposure(varlist) mediator(varlist)
control(string) baseline(string) base_confs(varlist) post_confs(varlist) impute(varlist)
imp_cmd(string) imp_eq(string) imp_cycles(#) simulations(#) samples(#) seed(#) obe all graph
saving(string) replace ]mainvarlist には、コマンドで使用するすべての変数を含める。ただし、変数名の省略形や、x1-x3 のように x1 x2 x3 を表す変数リストの使用はサポートされていないことに注意する。カテゴリカル変数は、その名前（例: agecat）のみを使用し、i.（例: i.agecat）のプレフィックスを使用しないで記載する必要がある。
データ構造
時間依存交絡因子オプションの場合（媒介オプションとは異なる—詳細は後述）、データはロング形式でなければならない（[D] reshape を参照）。つまり、各被験者の各時点に対して個別のレコードが必要である。アウトカムが時間依存のイベントの場合、各被験者のアウトカムデータは、図2で示されているように、各時点で測定された一連の二値変数として与える必要がある。
死亡や追跡脱落、または（時間依存のイベントアウトカムの場合）その時点以前にイベントを経験した被験者については、その時点以降のレコードをデータセットに含めるべきではない。
補完される値（断続的な欠測訪問の値を含む。これらの訪問についてはレコードを含める必要があります）は、Stataの規約に従って「.」で示す必要がある。
媒介オプションの場合、各被験者に正確に1つのレコードが必要である。同様に、補完される欠測値は「.」で示す必要がある。
各状況におけるデータの構造例は、第4節に記載されている。
オプション
時間依存交絡因子オプション
• outcome(varname)
varname をアウトカム変数として指定する。
• commands(string)
各パラメトリックモデルの適合に使用するコマンド（regress または logit）を指定する。変数名の後にコロン（\(:\)）、その後にコマンド名を記載し、異なる変数はコンマ（\(,\)）で区切る（例は第4節参照）。
アウトカム変数、時間変動交絡因子、時間変動曝露のモデルに対してコマンドを指定する必要がある。死亡によるセンサリングがある場合は、死亡のモデルに使用するコマンドも指定する。
時間依存アウトカム（および死亡の場合）では、アウトカム（または死亡）が一連の二値変数として与えられるため、logit が適切なコマンドである。
• equations(string)
上記で指定されたモデルの右辺に使用する方程式を指定する。従属変数の名前の後にコロン（\(:\)）、その後に独立変数のリストを記載し、異なる従属変数の方程式はコンマ（\(,\)）で区切る（例は第4節参照）。
データはロング形式で保存されているため、過去の訪問のデータに依存性を組み込むには遅延変数を使用する必要がある（後述参照）。
時間変動交絡因子\(L\)のような特定の変数の方程式は、各訪問で同じでなければならない。
方程式の右辺にカテゴリカル変数を使用する場合、変数の前に\(i.\)を付ける必要がある。
• idvar(varname)
varname を被験者を識別する数値変数として指定する。
• tvar(varname)
varname を時点を識別する数値変数として指定する。
• varyingcovariates(varlist)
時間変動共変量を varlist で指定する。遅延バージョンを使用する場合でも、非遅延バージョンのみをこのリストに含める。
• intvars(varlist)
介入が指定される変数を varlist で指定する。遅延バージョンを使用する場合でも、非遅延バージョンのみをこのリストに含める。
• interventions(string)
比較する具体的な介入を指定する。異なる介入はコンマ（\(,\)）で区切り、1つの介入内の異なるコマンドはバックスラッシュ（\(\backslash\)）で区切る（例は第4節参照）。
• dynamic
比較するレジームが動的であることを指定する。このオプションを指定しない場合、比較されるレジーム（観察レジームを除く）はすべて静的と見なされる。
• eofu
アウトカムが追跡終了時のみ測定されることを指定する。このオプションを指定しない場合、アウトカムは時間依存イベントであると見なされる。
• pooled
上記のcommandおよびequationオプション（該当する場合はimp_cmdおよびimp_eqオプションを含む）で定義されたモデルをすべての訪問データに対して一括して適合させることを指定する。このオプションを指定しない場合、各訪問ごとに個別にモデルを適合させる。
• death(varname)
各時点での二値変数（その時点で被験者が生存している場合は\(0\)、前の時点と現在の時点の間に死亡した場合は\(1\)）を指定する。死亡後のレコードは元のデータセットに含めるべきではない。
死亡\(=0\)である場合のセンサリングはすべて追跡脱落によるものと見なされる。この場合、死亡はあるが追跡脱落がない状況を模倣するシミュレーションが行われる。
deathオプションが指定されない場合、すべてのセンサリングは追跡脱落によるものと見なされ、追跡脱落がない状況を模倣するシミュレーションが行われる。
• derived(varlist)
他の変数から派生するすべての変数をリストする（例: 交互作用）。遅延変数そのものは含めないが、1つ以上の遅延変数を使用して派生した変数は含める。派生変数は元のデータセット内に存在する必要がある。
• derrules(string)
他の変数から派生変数を生成する方法を指定する。例: \(a\)と\(l\)の積として変数\(al\)を作成する場合、コードは\(\text{derrules}(al: a*l)\)とし、\(al\)を derived(varlist) に含める。複数の派生変数を生成するルールはコンマで区切る。
• fixedcovariates(varlist)
時間固定の共変量をリストする。これらは時間変動曝露に依存せず、シミュレーションされない。
• laggedvars(varlist)
遅延変数をリストする。遅延変数は元のデータセット内に存在する必要がある。
• lagrules(string)
遅延変数の詳細を指定する。例: 変数\(a_lag\)が\(a\)の遅延バージョンであり、\(a_lag2\)が二重遅延バージョンである場合、\(\text{lagrules}(a_lag: a 1, a_lag2: a 2)\)と記載する。
• msm(string)
マージナル構造モデルの形式を指定する。例: \(msm(regress\ y\ \text{a_lag}\ \text{a_lag2})\)または\(msm(stcox\ \text{a_lag}\ \text{a_lag2}\)\)。現在サポートされているのはregress、logit、およびstcoxのみである。このオプションは dynamic オプションと併用できない。
• impute(varlist)
単一確率的補完を使用して欠測値を補完する変数をリストする。
• imp_cmd(string)
各補完モデルを適合させる際に使用するコマンド（regress、logit、またはmlogit）を指定する。構文はcommandsオプションと同様である。
• imp_eq(string)
各補完モデルを適合させる際に使用する方程式の右辺を指定する。構文はequationsオプションと同様である。
• imp_cycles(#)
補完手続きで使用する連鎖方程式のサイクル数を指定する。デフォルトは10。
• simulations(#)
モンテカルロシミュレーションされたデータセットのサイズを指定する。デフォルトは観察データセットと同じサイズだが、計算上の理由でこれを小さくすることができる。
• samples(#)
ブートストラップサンプルの数を指定する。デフォルトは1,000。
• seed(#)
ランダムナンバーのシードを指定する。
• all
すべてのブートストラップ信頼区間（正規分布、パーセンタイル、バイアス補正、バイアス補正加速型）を表示することを指定する。デフォルトは正規分布に基づくブートストラップ信頼区間のみを表示する（[R] bootstrap を参照）。
• graph
各介入下での生存曲線のKaplan-Meierプロットを表示する。このオプションは時間依存交絡因子解析で時間依存アウトカムの場合にのみ関連する。
• saving(string)
元の観察データおよびすべてのモンテカルロシミュレーションを含むデータセットを、stringという名前のStataデータセットとして保存する。このデータセットには、変数_intが含まれている。この変数は、観察データの場合は値\(0\)、介入\(1\)に対応するシミュレーションの場合は値\(1\)、以下同様に指定された各介入について値を持つ。最後に、「介入なし」レジームの下でのモンテカルロシミュレーションがデータセットの末尾に追加され、その際の_intの値は\(m+1\)となる（\(m\)は指定された介入の数）。
• replace
saving(string) オプションで指定された.dtaファイルがすでに存在する場合、そのファイルを上書きすることを指定する。
媒介オプション
• mediation
分析が媒介分析であることを指定する。このオプションを指定しない場合、時間依存交絡因子解析が想定される。
• outcome(varname)
varname をアウトカム変数として指定する。
• commands(string)
シミュレーションの基礎となるパラメトリックモデルを適合させる際に使用するコマンド（regress または logit）を指定する。変数名の後にコロン（\(:\)）、その後にコマンド名を記載し、異なる変数はコンマ（\(,\)）で区切る（例は第4節参照）。
媒介変数、アウトカム、および曝露の影響を受ける媒介変数とアウトカムの関係の追跡後交絡因子のモデルを指定する必要がある。
• equations(string)
上記で指定されたモデルの右辺に使用する方程式を指定する。従属変数の名前の後にコロン（\(:\)）、その後に独立変数のリストを記載し、異なる従属変数の方程式はコンマ（\(,\)）で区切る（例は第4節参照）。
方程式の右辺にカテゴリカル変数を使用する場合、変数の前に i. を付ける必要がある。
• derived(varlist)
他の変数から派生するすべての変数（例: 交互作用）をリストする。
• derrules(string)
他の変数から派生変数を生成する方法を指定する。例: 変数\(al\)を\(a\)と\(l\)の積として作成する場合、コードは\(\text{derrules}(al: a*l)\)とし、\(al\)を\(**\text{derived(varlist)}**\)に含める。複数の派生変数を生成するルールはコンマで区切る。
• exposure(varlist)
曝露変数をリストする。
• mediator(varlist)
媒介変数をリストする。
• control(string)
制御された直接効果を計算する際に媒介変数を固定すべき値を指定する（例は第4節参照）。このオプションを指定しない場合、自然直接/間接効果のみが推定される。
• baseline(string)
曝露の基準値として使用する値を指定する（例は第4節参照）。
• base_confs(varlist)
曝露とアウトカムの関係の交絡因子をリストする。
• post_confs(varlist)
媒介変数とアウトカムの関係の交絡因子をリストする。
• impute(varlist)
単一確率的補完を使用して欠測値を補完する変数をリストする。
• imp_cmd(string)
各補完モデルを適合させる際に使用するコマンド（regress、logit、または mlogit）を指定する。構文はcommandsオプションと同様。
• imp_eq(string)
各補完モデルを適合させる際に使用する方程式の右辺を指定する。構文はequationsオプションと同様。
• imp_cycles(#)
補完手続きで使用する連鎖方程式のサイクル数を指定する。デフォルトは10。
• simulations(#)
モンテカルロシミュレーションされたデータセットのサイズを指定する。デフォルトは観察データセットと同じサイズである。
• samples(#)
ブートストラップサンプルの数を指定する。デフォルトは1,000。
• seed(#)
ランダムナンバーのシードを指定する。
• obe
曝露が\(1\)つの二値変数であることを指定する。この場合、比較は\(X = 1\)と\(X = 0\)の間で行われる。このオプションを指定しない場合、比較は観察データでの自然な\(X\)の分布と基準値の間で行われる。
• all
すべてのブートストラップ信頼区間（正規分布、パーセンタイル、バイアス補正、バイアス補正加速型）を表示することを指定する。デフォルトは正規分布に基づくブートストラップ信頼区間のみを表示する（[R] bootstrap を参照）。
• saving(string)
元の観察データおよびすべてのモンテカルロシミュレーションを含むデータセットを、stringという名前のStataデータセットとして保存する。
• replace
saving(string) オプションで指定された.dtaファイルがすでに存在する場合、そのファイルを上書きすることを指定する。
2つの例を用いた説明
例1：時間変化する共変量
データ
2つのデータセットが、1.1.3節で説明された方法に従って\(T = 9\)の条件でシミュレーションされている。1つ目のデータセットでは、\(A_0\)から\(A_9\)は二値の治療変数であり、訪問時点\(t\)において被験者が抗レトロウイルス療法（ART）を処方された場合は\(A_t = 1\)、そうでない場合は\(A_t = 0\)である。\(L_0\)から\(L_9\)は各訪問時点でのCD4カウントの対数値である。\(Y_1\)から\(Y_{10}\)は二値変数であり、時点\(t\)間隔\((t-1, t]\)において被験者がAIDSを発症した場合は\(Y_t = 1\)、そうでない場合は\(Y_t = 0\)である。すべての被験者はベースライン時点でAIDSフリーであり（したがって\(Y_1\)は\(Y\)の最初の記録値となる）、もし\(Y_t = 1\)であれば、その時点以降\(t+1\)からの記録はその被験者に対して含まれない。以下は、1つ目のデータセットにおける最初の3名の被験者のデータである。
被験者\(1\)および\(3\)は追跡終了時までAIDSを発症しなかった。被験者\(2\)は時点\(6\)と\(7\)の間にAIDSを発症した。\(\text{a_log}\)は変数\(a\)の遅延版であり、時点\(0\)を除いて\(a\)の前回の値を含む。時点\(0\)では、すべての被験者において\(a\)は0である。同様に、\(\text{l_lag}\)は変数\(l\)の遅延版である。\(\text{cuma}\)は時点\(t\)におけるその被験者の値\(a\)の時点\(t\)までの累積和を表す。また、\(\text{cuma_log}\)はその遅延版である
データ（1,000人の被験者で構成）は以下のように生成された。
• \(U\)は平均\(0\)、分散\(0.25\)の正規乱数変数である。
• \(L_0\)は平均\(5.5 + U\)、分散\(0.04\)の正規乱数変数である。
• \(A_0\)は確率$$\frac{\exp{(5-L_0)}}{1 + \exp{(5-L_0)}}$$に基づいてベルヌーイ分布から生成される。
• その後、各\(t \in [1, 9]\)において、かつ\(Y_{t-1} = 0\)の被験者について、\(Y_t, L_t\)および\(A_t\)は以下のように生成される。\(Y_t\)は次の確率に基づいてベルヌーイ分布から生成される。$$\frac{\exp{(-8 + L_{t-1}-0.3\sum_{s=0}^{t-1}{A_s}-U)}}{1+\exp{(-8+L_{t-1}-0.3\sum_{s=0}^{t-1}{A_s}-U)}}$$\(L_t\)は次の平均と分散を持つ正規分布から生成される。平均\(0.9L_{t-1} + A_{t-1} + 0.1 U\)、分散\(0.01\)。\(A_t\)は次の確率に基づいてベルヌーイ分布から生成される。$$\frac{\exp{(A_{t-1}+4.5L_t)}}{1+\exp{(A_{t-1}+4.5L_t)}}$$
• 最後に、\(Y_9 = 0\)の被験者について、\(Y_{10}\)は次の確率に基づいてベルヌーイ分布から生成される。$$\frac{\exp{(-8 + L_{9}-0.3\sum_{s=0}^{9}{A_s}-U)}}{1+\exp{(-8+L_{9}-0.3\sum_{s=0}^{9}{A_s}-U)}}$$
2つ目のデータセットは、死亡および追跡脱落の両方によるセンサリングがある点を除いて、まったく同じ方法で生成される。すべての被験者が時点\(0\)では観察されている。その後、時点\(t\)における追跡脱落は、次の平均を持つベルヌーイ乱数変数として生成される。$$\frac{\exp{(-10 + L_{t-1}-0.3\sum_{s=0}^{t-1}{A_s}-U)}}{1+\exp{(-10+L_{t-1}-0.3\sum_{s=0}^{t-1}{A_s}-U)}}$$
死亡または追跡脱落が発生していない場合、 \(Y_t, L_t\)、および\(A_t\)は前述の方法で生成される。以下は、この\(2\)つ目のデータセットからの\(4\)人の被験者のデータである（\(d\)は死亡を示す変数である）。
被験者\(1\)は訪問\(1\)と\(2\)の間に死亡した。被験者\(2\)は追跡終了時までAIDSを発症しなかった。被験者\(9\)は訪問\(5\)と\(6\)の間にAIDSを発症した。被験者\(30\)は訪問\(7\)の後に追跡脱落した。
コマンド
g-computationは、以下のコマンドを使用して、1つ目のデータセットに適応される。
gformula y a l a_lag l_lag cuma cuma_lag id t, out(y) com(y:logit, ///
l:regress, a:logit) eq(y:l_lag cuma_lag, l:l_lag a_lag, a:l a_lag) ///
id(id) t(t) var(l) intvars(a) interventions(a=1 if t<10, ///
a=0 if t<=1 \ a=1 if t>1 & t<10, a=0 if t<=3 \ a=1 if t>3 & t<10, ///
a=0 if t<=5 \ a=1 if t>5 & t<10, a=0 if t<=7 \ a=1 if t>7 & t<10, ///
a=0 if t<=9) pooled lag(l_lag a_lag cuma_lag) lagrules(l_lag: l 1, ///
a_lag: a 1, cuma_lag: cuma 1) msm(stcox cuma_lag) derived(cuma) ///
derrules(cuma:cuma_lag+a) seed(79)6つの静的レジームが比較される。
- 1. \((A_0, A_1, A_2, A_3, A_4, A_5, A_6, A_7, A_8, A_9) = (1, 1, 1, 1, 1, 1, 1, 1, 1, 1)\)
- 2. \((A_0, A_1, A_2, A_3, A_4, A_5, A_6, A_7, A_8, A_9) = (0, 0, 1, 1, 1, 1, 1, 1, 1, 1)\)
- 3. \((A_0, A_1, A_2, A_3, A_4, A_5, A_6, A_7, A_8, A_9) = (0, 0, 0, 0, 1, 1, 1, 1, 1, 1)\)
- 4. \((A_0, A_1, A_2, A_3, A_4, A_5, A_6, A_7, A_8, A_9) = (0, 0, 0, 0, 0, 0, 1, 1, 1, 1)\)
- 5. \((A_0, A_1, A_2, A_3, A_4, A_5, A_6, A_7, A_8, A_9) = (0, 0, 0, 0, 0, 0, 0, 0, 1, 1)\)
- 6. \((A_0, A_1, A_2, A_3, A_4, A_5, A_6, A_7, A_8, A_9) = (0, 0, 0, 0, 0, 0, 0, 0, 0, 0)\)
2回目の分析では、同じデータセットを使用するが、動的レジームを比較する。以下がそのコードである。
gformula y a l a_lag l_lag cuma cuma_lag id t, out(y) com(y:logit, ///
25l:regress, a:logit) eq(y:l_lag cuma_lag, l:l_lag a_lag, a:l a_lag) ///
id(id) t(t) var(l) intvars(a) interventions(a=0 if t<10 & l>6.9 ///
\ a=1 if t<10 & l<=6.9, a=0 if t<10 & l>6.55 \ a=1 if t<10 & l<=6.55, ///
a=0 if t<10 & l>6.2 \ a=1 if t<10 & l<=6.2, a=0 if t<10 & l>5.3 \ ///
a=1 if t<10 & l<=5.3, a=0 if t<10 & l>4.6 \ a=1 if t<10 & l<=4.6) ///
dynamic pooled lag(l_lag a_lag cuma_lag) lagrules(l_lag: l 1, ///
a_lag: a 1, cuma_lag: cuma 1) derived(cuma) derrules(cuma:cuma_lag+a) ///
seed(801)比較される動的レジームはすべて「時点\(t\)において\(L_t < x\)の場合にのみ治療を行う」という形式であり、比較される\(5\)つの異なるレジームにおいて\(x\)の値はそれぞれ\(6.9, 6.55, 6.2, 5.3, 4.6\)となっている。最後に、追跡脱落と死亡によるセンサリングを含む\(2\)つ目のデータセットを分析し、上記に記載された同じ\(6\)つの静的レジームを以下のコードを使用して比較する。
gformula y a l a_lag l_lag d cuma cuma_lag id t, out(y) com(y:logit, ///
l:regress, a:logit, d:logit) eq(y:l_lag cuma_lag, l:l_lag a_lag, ///
a:l a_lag, d:l_lag cuma_lag) id(id) t(t) var(l) intvars(a) ///
interventions(a=1 if t<10, a=0 if t<=1 \ a=1 if t>1 & t<10, ///
a=0 if t<=3 \ a=1 if t>3 & t<10, a=0 if t<=5 \ a=1 if t>5 & t<10, ///
a=0 if t<=7 \ a=1 if t>7 & t<10, a=0 if t<=9) pooled ///
lag(l_lag a_lag cuma_lag) lagrules(l_lag: l 1, a_lag: a 1, ///
cuma_lag: cuma 1) msm(stcox cuma_lag) derived(cuma) ///
derrules(cuma:cuma_lag+a) death(d) seed(79)出力
以下は、最初の分析（追跡脱落も死亡もない場合の静的レジームの比較）の（省略された）出力である。
\(3\)つの表はすべて治療の有益な効果を示している。すなわち、被験者が治療を受ける回数が多いほど、AIDSを発症せずに生存する期間が長くなる。これは、MSMの結果から得られる累積治療に関連する負のログハザード比（HR \(0.805\)、\(95\%\) CI \([0.745, 0.870]\)に対応）や、他の\(2\)つの表で下に進むにつれて増加する平均ログ発生率および累積発生率から確認できる。治療が行われなかった場合（介入\(6\)）には、研究参加者の\(43\%\)が仮想研究中にAIDSを発症するとシミュレーションされたが、治療が常時処方された場合にはわずか\(17\%\)しかAIDSを発症しないとシミュレーションされた。
観察レジームの下で、シミュレーションデータと観察データの間にはわずかな差がある（累積発生率は\(31\%\)対\(28\%\)、平均ログ発生率は\(-3.24\)対\(-3.37\)）。しかし、これらの差は標準誤差に比べて小さいため、このチェックに起因する懸念はない。以下は、\(2\)回目の分析（動的レジームの比較）の（省略された）出力である。
上述の通り、この分析からはマージナル構造モデルのパラメータを推定することはできない。しかし、平均ログ発生率および累積発生率の結果は治療が有益であることを確認しており、ART（抗レトロウイルス療法）が投与される閾値\(x\)が最も高い動的レジームの下で、AIDSフリーの生存期間がより長いことが示されている。観察レジームの下でのシミュレーションデータと観察データの一致は非常に良好であり、観察レジームはAIDSフリー生存期間に関してレジーム\(4\)とレジーム\(5\)の間に位置していると示唆されている。最後に、追跡脱落および死亡によるセンサリングを含む\(3\)回目の分析の（省略された）出力を以下に示す。
この分析から得られる結論は似ているが、死亡が競合するイベントと見なされるため、解釈はより難しくなっている。また、観察レジームの下でシミュレーションデータと観察データを比較することも難しくなっている。というのも、シミュレーションデータには追跡脱落が含まれていないのに対し、観察データには追跡脱落が含まれているためである。
標準的な解析との比較
以下に、累積治療を考慮したAIDSフリー生存期間に対する標準Cox回帰分析を示す。分析は時間変動交絡因子である\(\text{log(CD4)}\)を調整した場合と調整しない場合で行われている。これらは、最初のシミュレーションデータセット（死亡や追跡脱落によるセンサリングなし）の結果である。
両方の分析は、治療が有害であることを示唆している（ただし、調整された分析ではこの効果に対する証拠は非常に弱い）。データをシミュレーションした方法から、治療が有益であることは分かっており、このことはg-computationによって確認されている。上記の未調整の分析は、治療を受けた被験者が任意の訪問時点で治療を受けていない被験者より健康状態が悪い（治療を行うかどうかの決定がその訪問時点での\(\text{CD4}\)カウントに依存するため）ことを考慮していないため、バイアスがかかっている。
しかし、\(\text{log(CD4)}\)を調整すると、治療の多くの有益な効果が見えなくなるため、治療が有害であると誤って示唆されている（将来の\(\text{CD4}\)カウントを条件付けることで、治療の効果が隠れてしまうためである）。
同様の結果（実際にはより極端なもの）が、\(2\)つ目のデータセット（追跡脱落および死亡によるセンサリングを含む）に対して標準的な分析を行った場合にも見られる。
例2：媒介
データ
セクション1.2.5で説明された方法に従い、\(10,000\)人の被験者を含むデータセットが以下のようにシミュレーションされた。
• 社会経済的地位（SEP） は、\(30\%\)の被験者に対して\(1\)（低）、\(50\%\)に対して\(2\)（中）、残り\(20\%\)に対して\(3\)（高）として生成される。
• アルコール摂取量（A）（1日あたりの単位数）は、ゼロ膨張した歪んだ分布として生成される。確率\(0.8I(SEP = 1) + 0.7I(SEP = 2) + 0.9I(SEP = 3)\)に基づき、ベルヌーイ乱数変数を生成する。この二値変数が\(0\)の場合、\(A = 0\)とする。それ以外の場合、\(\text{log(A)}\)は平均\(I (SEP = 1) + 0.7I(SEP = 2) + 1.2 I(SEP = 3)\)、分散\(0.25\)の正規分布から生成される。
• 体格指数（BMI）の対数（log(BMI)） は、平均\(23 + I(SEP = 1) + 0.4A\)、分散\(4\)の正規分布から生成される。
• GGTの対数（log(GGT）、単位はグラム/リットル）は、平均\(2.5+0.02BMI + 0.1A\)、分散\(1\)の正規分布から生成される。
• 最後に、収縮期血圧（SBP、単位はmmHg） は、平均\(80 + 0.5BMI + 6A + 7\text{log(GGT)} -\text{log(GGT)A}- 5(SEP -3)\)、分散\(100\)の正規分布から生成される。
独立して完全にランダムに、被験者の\(5\%\)はアルコール変数が欠測、\(5\%\)はBMI変数が欠測、\(5\%\)はGGT変数が欠測している。その結果、\(8,620\)人の被験者が完全なデータを持ち、\(442\)人はGGTのみ欠測、\(424\)人はBMIのみ欠測、\(427\)人はアルコールのみ欠測している。さらに、\(26\)人はアルコールとBMIの両方が欠測（GGTは観測済み）、\(26\)人はアルコールとGGTの両方が欠測（BMIは観測済み）、\(34\)人はBMIとGGTの両方が欠測（アルコールは観測済み）である。最後に、\(1\)人の被験者は\(3\)つの変数すべてが欠測している。
データセット内の最初の\(3\)人の被験者（すべて\(SEP = 1\)）のデータが以下に示されている。\(\text{log_ggt~c}\)は\(\text{log_ggt_alc}\)（\(\text{log_ggt}\)と\(\text{alc}\)の積）の略であり、\(\text{alc_sbp}\)は\(\text{alc}\)と\(\text{sbp}\)の積、\(\text{log_ggt~p}\)は\(\tect{log_ggt_sbp}\)（\(\text{log_ggt}\)と\(\text{sbp}\)の積）の略である。
コマンド
g-computationは、以下のコマンドを使用して適用された。
gformula sep alc bmi log_ggt sbp log_ggt_alc alc_sbp log_ggt_sbp, ///
mediation out(sbp) eq(bmi:i.sep alc, log_ggt:bmi alc, sbp:bmi alc ///
log_ggt log_ggt_alc i.sep) com(bmi:regress, log_ggt:regress, sbp:regress) ///
ex(alc) mediator(log_ggt) control(log_ggt:3) baseline(alc:0) ///
post_confs(bmi) base_confs(sep) derived(log_ggt_alc alc_sbp log_ggt_sbp) ///
derrules(log_ggt_alc:log_ggt*alc, alc_sbp:alc*sbp, log_ggt_sbp:log_ggt*sbp) ///
impute(alc bmi log_ggt) imp_cmd(alc:regress, bmi:regress, log_ggt:regress) ///
imp_eq(alc:i.sep bmi log_ggt sbp log_ggt_sbp, bmi:i.sep alc log_ggt sbp ///
log_ggt_alc, log_ggt:i.sep alc bmi sbp alc_sbp) seed(79)出力
以下が（省略された）出力になる。
ここでの結論は、アルコール摂取が収縮期血圧に因果的な影響を与えるということである。もし全員が飲酒をやめた場合、平均収縮期血圧（SBP）は\(7.51\)単位低下する（\(95\%\)信頼区間\([7.08, 7.95]\)）と推定される。この低下のごく一部（\(1.13\)単位）はGGTを介した媒介効果によるものである。影響の大部分は直接的なものであり、BMIや他の経路を通じて作用している。
標準的な解析との比較
以下は、導入部で説明したように、このデータに対して使用される可能性がある標準的な分析である。欠測データに対処するために、上記と同じ補完モデルを使用して連鎖方程式による多重補完を行う。各欠測値に対して\(5\)回の適切な補完を生成し（Stataのiceを使用）、結果をmimコマンドで分析および統合する。今回は、標準誤差をブートストラップではなく解析的に推定するため、適切な補完を複数回行う必要がある。
これらの推定値は、g-computationを用いて得られた推定値と直接比較することはできない。上記の分析でのアルコールの係数は、\(1\)日あたりの摂取量の単位変化に対応しているためである。シミュレーションデータセットにおける\(1\)日あたりの平均摂取量は\(2.22\)単位であるため、標準回帰で推定される総効果の相当値はおおよそ\(3.27 \times 2.22 = 7.25\)となり、予想されるように、先に得られた\(7.51\)と類似している。しかし、\(2\)回目の回帰分析でのアルコールの係数を単純に解釈すると、これはGGTを介さない直接効果を表すと考えられ、間接効果はおおよそ\(7.25-(2.54\times 2.22) = 1.62\)となり、この分析ではg-computationよりも大きい値が得られるように見える。言い換えれば、標準的な分析では、効果のより大部分がGGTを介して媒介されていると誤って結論付ける可能性がある。これは予想されることであり、アルコールがSBPに及ぼす直接効果（すなわち、GGTを介さない効果）の一部がBMIを介して作用しているためである。この部分の効果は、上記の標準分析では正しく割り当てられておらず、直接効果の過小評価につながる。
計算時間に関する注意
gformulaコマンドは計算負荷が非常に高く、時点数が増加するにつれて計算時間が指数関数的に増加する。上記の時間依存交絡因子の例では、\(T = 9\)の場合、パラメトリックモデルを適合させ、各介入下でデータをシミュレーションし、シミュレーションされた各データセットを分析するのに約\(30\)秒が必要である。したがって、\(1,000\)回のブートストラップサンプルが必要な場合、全体の分析には\(8\)時間以上かかる。しかし、ブートストラップはタスクの分散処理に非常に適しており、高性能コンピュータクラスタを使用することで処理時間を大幅に短縮できる。
最終的な考察
時間変動交絡因子や媒介に関する問題では、曝露によって交絡因子が影響を受ける場合、標準的な回帰分析が無効であることを繰り返し強調してきた。一方で、g-computationは、過去の曝露によって交絡因子が影響を受けることを許容する、より弱い仮定の下で有効である。この手法が有効であるために必要な構造的仮定は、十分な交絡因子セットが測定されていることである。さらに、この手法では、観察データにおける追跡開始後の変数について正しいパラメトリックモデルを仮定する必要がある。
代替的なセミパラメトリックモデルおよび推定方法も提案されています（Robins et al., 1992, 2000）。これらは、構造ネストモデルのg推定法およびマージナル構造モデルの逆確率重み推定法である。これらの代替手法は、パラメトリックモデリングの仮定が少なく、モデルの誤指定バイアスに対してより強靭である。また、これらのセミパラメトリックアプローチではモンテカルロシミュレーションを必要とせず、計算負荷が軽減される。これらの方法のStataでの実装も報告されている（Sterne and Tilling, 2002; Fewell et al., 2004）。
しかしながら、g-computationにはこれらの方法に対する明確な利点がある。例えば、複雑な多変量（または結合的）介入をより簡単に扱うことができる点や、仮定された介入を観察レジームと比較することが可能である点である。これらは政策決定の指針となる場合に重要である（Taubman et al., 2009）。これらの利点に加え、g-computationは、より強いモデリング仮定の代償として、統計的効率性を向上させるという利点もある（Daniel et al., 2010）。
我々は、g-computationが多くの状況において価値のあるツールであると考えている。この手法は1986年にRobinsによって最初に提案されたが、その見た目の複雑さやソフトウェアルーチンの不足により、近年のSASにおけるGFORMULAマクロ（Taubman et al., 2009）まで広く利用されることはなかった。このStataルーチンが、この方法をより多くの応用研究者に利用しやすくすることを願っている。
コメント