21  G-Methods for Time-Varying Treatments

이전 챕터에서 시간에 따라 변화하는 처치와 그 처치와 교란변수 간에 일종의 피드백(feedback)이 존재할 경우에 기존에 시점이 고정된 상태에서 교란변수들을 조정해주는 것만으로는 충분하지 않다는 것을 확인하였다. 이러한 피드백이 존재할 경우 기존의 방식으로 조정하더라도 인과효과가 없어야 하는(null) 조건에서도 인과 효과가 존재하는 것으로 추정된다는 것이다. 이 챕터에서 저자들은 g-formula, IP 가중치, g-estimation 등과 이중-강건 일반화를 통해 이 세 가지 g-methods가 인과 효과가 없을 경우를 제대로 잡아낸다는 것을 보여주고 있다.

21.1 The g-formula for time-varying treatments

이 섹션은 시간에 따라 변하는 처치를 평가하기 위한 g-formula에 대해 설명하고 있다. 예를 들어, 시간에 따라 변하는 치료를 받는 환자들의 데이터를 가지고 있을 때, 전통적인 방법으로는 교란변수의 영향력을 제대로 조정할 수 없다는 것을 보여준다. 반면, g-formula를 사용하면, 시간에 따라 변하는 처치가 있는 경우, 처치와 교란변수 사이의 피드백이 있다 하더라도 올바른(무효) 효과에 대한 추정치를 얻을 수 있다.

시점이 고정된 처치의 효과에 관심이 있다고 하자. 이 경우, 우리가 비교하고 싶은 것은 평균적인 반사실적 결과 \(E[Y^{a_1=1}], E[Y^{a_1=0}]\)라고 할 수 있다. 시계열적인 무작위 실험이 수행된 경우, 처치인 \(\bar{A} = (A_0, A_1)\)은 시간에 따라 변화하며, 이 경우 처치-교란변수 간 피드백이 존재하게 된다. 시계열적으로 수행된 무작위 실험이 시계열적 교환가능성의 조건을 충족한다고 하더라도, 이전 처치가 현재 처치에 미치는 영향은 기존의 방식으로는 조정될 수 없다는 점에서 편향된 추정치를 산출하게 된다. 하지만 g-formula는 시계열적인 무작위 실험에서 반사실적 평균 \(E[Y^{a_0, a_1}]\)를 계산하는 데 사용될 수 있다.

• 간단히 말해서, g-formula는 시간이 지남에 따라 변할 수 있는 처치와 공변량(교란변수)을 고려하여 결과의 조건부 평균을 가중 평균으로 계산하고, 이를 공변량 분포에 표준화하여 조정한다.

  • g-formula에서 사용하는 가중치는 과거의 처치에 대한 정보를 바탕으로 현재의 공변량 \(L_1\)의 분포를 나타낸다. 여기서 ’과거의 처치값’이란 이전 시점에서의 처치 상태를 의미한다.

• ’Sequential exchangeability’는 연구 대상이 시간에 따라 여러 번 처치를 받을 때, 각 시점에서 처치받은 그룹과 처치받지 않은 그룹이 비교 가능하다는 가정이다. 즉, 처치의 순서와 결과 사이에 다른 영향을 주는 교란변수가 없다는 것을 의미한다.

• 저자들은 또한 가상의 모집단을 시뮬레이션 하는 방법과 g-formula의 값이 어떻게 결정되는지를 설명하고 있다. ’가상의 모집단(counterfactual population)’이란 특정 처치 전략이 전체 모집단에 적용되었을 때의 가상 시나리오를 일컫는다.

  가상 모집단을 시뮬레이션하기 위해 다음과 같은 두 단계가 요구된다:

  1. 각 시점 \(k=0, k=1\)에서 처치 \(a_0 = 1\), \(a_1 = 1\)를 받을 확률을 1이라 하자. 이는 처치 전략이 모든 시점에서 “항상 처치가 이루어지는 경우”를 나타낸다.

  2. 원래 연구 인구에서 사용된 것과 동일한 확률 \(\Pr[L_1=l_1|A_0=a_0]\)과 동일한 평균 \(E[Y|A_0=a_0, A_1=a_1, L_1=l_1]\)을 가상 모집단에 할당한다.

  이때 g-formula의 값은 \(L\)에 무엇이 포함되었는지, 또는 포함되었다면 무엇인지에 따라 달라진다. 처치 \(A_1\)이 \(L_1\) 에 의해 영향을 받는 경우, \(L_1\)을 포함하지 않는 g-formula는 더 이상 인과적으로 해석할 수 없다.

  또한, g-formula가 인과 해석을 가질 수 있더라도 그 구성 요소들 각각이 인과적으로 해석되지 않을 수도 있다.

여러 시점 \(k\) 를 가진 고차원(high-dimensional) 조건에 g-formula를 일반화하는 공식은 다음과 같다:

\[ \sum_l E[Y|\bar{A}=\bar{a}, \bar{L}=l] \prod^{K}_{k=0} f(l_k|\bar{a}_{k-1}, \bar{l}_{k-1}) \]

• \(E[Y∣A=a,L=l]\)는 주어진 처치 \(A\)와 공변량 \(L\)의 조건 하에 결과 \(Y\)의 기대값을 나타낸다.

• \(f(l_k​∣a_{k−1}​,l_{k−1}​)\)는 시간 \(k\)에서 공변량 \(L\)의 조건부 분포를 나타내며, 이전 시간 \(k−1\)에서의 처치와 공변량에 기초한다.

• 곱셈 부분(production)은 시간에 따른 공변량의 조건부 분포를 모델링하며, 이는 처치가 시간에 따라 어떻게 공변량에 영향을 미치는지를 보여준다.

그런데 이때 자료가 지나치게 고차원일 경우(공변량이 많을 경우 등), g-formula 의 각 요소들을 직접적으로 컴퓨팅하기는 어렵다. 따라서 관측된 데이터의 분포에 특정한 함수적 형태를 가정하는 모델(선형회귀나 로지스틱회귀)을 사용하여 g-formula 의 구성 요소를 추정하는 선택을 할 수 있다. 이를 plug-in-g-formula 라고 하며, 모수적 모델로부터 얻은 추정치를 g-formula에 대입하는 방식을 일컫는다.

처치 전략 \(f^{int}\) 하에서 시점 \(k\) 에서의 처치 \(a_k\) 에 대한 조건부 확률을 반영하여 일반화한 g-formula는 다음과 같다:

\[ \sum_{a,l} E[Y|\bar{A}=\bar{a}, \bar{L} = \bar{l}]\prod^K_{k=0}f(l_k|\bar{a}_{k-1}, \bar{l}_{k-1})\prod^K_{k=0}f^{int}(a_k|\bar{a}_{k-1},\bar{l}_k). \]

• \(\sum_{a,l}\) : 모든 가능한 처치 조합 \(a\) 와 공변량 \(l\) 에 대한 합산으로 가능한 모든 시나리오를 고려한다는 것을 의미한다.

• \(E[Y|\bar{A}=\bar{a} , \bar{L} = \bar{l}]\): 주어진 처치 \(a\) 와 공변량 \(l\) 의 조건 하에서 결과 \(Y\) 의 기대값으로 결과의 평균을 나타낸다.

• \(\prod^K_{k=0}\): 시점 \(k\) 전체에 걸친 변화를 나타낸다.

• \(f(l_k|\bar{a}_{k-1}, \bar{l}_{k-1})\): 시간 \(k\)에서의 공변량 \(L\)의 분포로, 이전 시간 \(k−1\)의 처치 \(A\)와 공변량 \(L\)에 조건부로 정의된다.

• \(f^{int}(a_k|\bar{a}_{k-1},\bar{l}_k)\): 가상의 처리 전략 \(f^{int}\) 하에서 시간 \(k\)에서의 처치 \(A\)의 조건부 분포를 나타낸다 이는 주어진 공변량 \(L\)과 이전 시간의 처치 \(A\)에 기초한다.

21.2 IP weighting for time-varying treatments

시간에 따라 변하는 처치와 공변량을 가진 상황에서 역확률 가중치(Inverse Probability Weights, IPW)를 어떻게 일반화할 수 있을까? 역확률 가중치는 관찰 연구에서 처치 효과를 추정할 때, 처치를 받을 확률에 따라 각 개체의 데이터를 가중하는 방법으로 두 가지 종류로 나누어 살펴볼 수 있다.

1. 비안정화(nonstabilized) IP 가중치 \((W^{\bar{A}})\): 이는 각 시점에서 처치를 받을 조건부 확률의 역수로 계산된다. 이 가중치는 개인이 각 시점에서 처치를 받을 확률을 기반으로 하여, 처치를 받지 않았을 경우에 개인이 받았을 결과를 가중한다: \(W^{\bar{A}} = \frac{1}{f(A_0|L_0) \times f(A_1|A_0, L_0, L_1)} = \prod_{k=0}^{1} \frac{1}{f(A_k | A_{k-1}, L_k)}\)
2. 안정화(stabilized) IP 가중치 \((SW^{\bar{A}})\): 이는 처치를 받을 한계확률과 처치를 받을 조건부 확률의 비율로 계산된다. 안정화 가중치는 표본에서 처치 받은 사람들과 받지 않은 사람들의 분포를 보정하는 데 사용된다: \(SW^{\bar{A}} = \frac{f(A_0)}{f(A_0|L_0)} \times \frac{f(A_1|A_0)}{f(A_1|A_0, L_0, L_1)} = \prod_{k=0}^{1} \frac{f(A_k | A_{k-1})}{f(A_k | A_{k-1}, L_k)}\)

여기서 \(A_0, A_1\) 은 각각 시점 0과 1에서의 처치를, \(L_0, L_1\)​은 해당 시점에서의 공변량을 나타낸다. \(f\) 는 조건부 확률을 의미한다. \(A_{-1}\) ​은 정의에 따라 0이다. 276쪽에서 예제 데이터인 Table 21.1에 역확률 가중치를 적용하여 살펴보면, g-formula 와 마찬가지로 인과 효과를 타당하게 추정하는 것을 확인할 수 있다.

단, 역확률 가중치의 계산과정에서 양의 확률(positivity) 가정이 위반되거나 계산된 weight의 극단 값으로 인한 추정치의 불확실성이 클 수 있다.

또한, 문제는 실제 연구, 혹은 관찰 연구 자료에서는 우리가 보정해야할 것들을 제대로 보정했는지 확신할 수 없다는 것에 있다. 즉, 모델을 잘못 특정할 경우(misspecification)에는 결과적으로 편향된 추정치를 얻게 된다는 것이다. 따라서 과연 우리가 모델을 제대로 특정하였는지를 판단하기 위한 방법으로 저자들은 이중-강건 추정량(doubly-robust estimators)을 살펴볼 것을 권하고 있다.

21.3 A doubly robust estimator for time-varying treatments

처치 효과를 모델링하는 모델과 결과에 대한 모델 중 적어도 하나가 정확하다면, 이 이중-강건 추정량을 사용하여 올바른 인과 효과를 추정할 수 있다는 가능성에 기초하고 있다. 시간에 따라 변하는 처치를 가진 자료를 다룰 때, 처치의 이력(history) \(\bar{a}\)가 평균 결과에 미치는 효과는 처치의 누적량 \(\mathrm{cum}(\bar{a})=\sum^K_{k=0} a_k\) ​​에 대해 선형적으로 증가한다고 가정할 수 있다.

\[ E[Y|\bar{a}] = \beta_0 + \beta_1\mathrm{cum}(\bar{a}) \]

• 위의 일반화된 한계 구조적 평균 모델에서 \(\beta_1\) 은 시간에 따라 변화하는 처치 \(A\) 의 평균 인과효과를 보여주는 모수이다.

• 이때, 평균 인과효과는 \(E[Y|\bar{a}] - E[Y|\bar{a} = 0]\) 는 \(\beta_1\times\mathrm{cum}(\bar{a})\) 와 같다.

이중-강건 추정량은 결과 모델의 추정량 \(E[Y|A=a, L=l]\) 과 처치 모델의 추정량 \(\Pr[A=1|L]\) 을 결합하여 산출한다. 이는 세 가지 단계를 따르는데:

1. 전체 관측 결과에 대해서 모델을 적합한다(pooled).

2. 각 시점 \(m\) 에서의 별개의 결과 모델 \(b_m (\bar{L}_m; \beta_m)\) 을 적합한다.

3. 값 \(\hat{E}[Y^{\bar{a} = \bar{1}}]\) 에 대한 추정치를 모든 개인에 대해 컴퓨팅하고 표본평균으로 계산한다.

   • 만약 (1) 모든 \(m\) 에 대해 결과 모델이 올바르게 특정되었거나, (2) 모든 \(m\) 에 대해 처치 모델이 올바르게 특정되었다면, \(\hat{E}[Y^{\bar{a} = \bar{1}}]\) 는 \(E[Y|\bar{a} = 1]\) 에 대해 편향되지 않을 것이다.

이러한 이중-강건 추정량은 종종 Targeted minimum loss-based estimator (TMLE)라고도 한다.

21.4 G-estimation for time-varying treatments

여기서는 시간에 따라 변하는 처치의 인과 효과를 추정하기 위한 G-estimation 방법을 다룬다. 처치의 시간 변동성을 고려하여, 처치와 결과 간의 관계를 모델링하는 데 필요한 복잡성을 해결하는 다양한 방법들이 제시된다.

첫 번째 단계에서는 구조적 중첩 평균 모델(Structured Nested Mean Model)을 사용하여 각 시간 단계에서의 처치가 결과에 미치는 평균적인 효과를 모델링한다. 이 모델은 처치가 시간에 따라 변할 수 있음을 고려하여, 처치의 각 시간 단계가 결과에 미치는 효과를 모델링에 고려한다. 예를 들어, 시간에 따라 변하는 처치 \(A\) 와 결과 \(Y\) 사이의 관계를 모델링하는 공식은 다음과 같다:

\[ Y_{a_0, a_1} = Y_{a_0, 0} + \psi_{11}a_1 + \psi_{12}a_1L_1 + \psi_{13}a_1a_0 + \psi_{14}a_1a_0L_1 \]

• 여기에서 \(Y_{a_0, a_1}\) 은 처치 \(a_0\) 과 \(a_1\) 을 받았을 때의 결ㅏ를 나타내며, \(\psi\) 모수들은 처치의 다양한 조합이 결과에 미치는 효과를 크기를 나타낸다.

두 번째 단계에서는 G-estimation 기법을 사용하여 모델의 모수(parameters)를추정한다. G-estimation은 관찰된 데이터를 바탕으로 이론적 모델의 모수를 추정하는 과정으로, 시간에 따라 변하는 처치와 결과 사이의 관계를 정량화하기 위해 사용된다.

마지막으로, G-estimation을 통해 얻어진 결과를 해석하고, 이를 바탕으로 시간에 따라 변하는 처치의 인과 효과에 대한 추정치를 제공한다.

예를 들어, 289쪽의 내용을 살펴보자. 시간에 따라 변하는 처치를 가진 구조적 중첩 평균 모델의 G-estimation 을 산출하고자 한다고 하자. 이제까지는 단순화를 위해, 처치의 효과가 시간 \(k\) 또는 처치 및 공변량의 이력에 따라 변하지 않는다고 가정하는 단일 파라미터 \(\beta_1\)​을 가진 모델을 고려해보자. 그러나 이제 \(\beta_1\) 를 벡터라고 가정하고, 모델이 5차원이라고 가정할 때, 다음과 같은 모델을 고려할 수 있다:

\[ \gamma_k (\bar{a}_{k-1}, l_k, \beta) = \beta_0 + \beta_1 k + \beta_2 a_{k-1} + \beta_3 l_k + \beta_4 l_k a_{k-1} \]

• 여기서 \(\bar{a}_{k-1}\) 은시간 \(k-1\)까지의 처치 역사, \(l_k\)​는 시간 \(k\)에서의 공변량을 의미한다.

• 처치 \(a_k\)​가 결과에 미치는 효과를 시간, 이전의 처치, 그리고 공변량에 따라 모델링한다.

• 5개의 파라미터 \(\beta\)를 추정하기 위해, 다음과 같은 로지스틱 회귀 모델을 적합시킬 수 있다:

\[ \text{logit} \Pr[A_k = 1 | H_k (\psi^{\dagger}), L_k, \bar{A}_{k-1}] = \alpha_0 + H_k (\psi^{\dagger})(\alpha_1 + \alpha_2 k + \alpha_3 A_{k-1} + \alpha_4 L_k + \alpha_5 L_k A_{k-1}) + \alpha_6 W_k \]

여기서 \(H_k(\psi^{\dagger})\) 의 임의의 추정치 \(\psi^{\dagger}\) 를 사용하여 계산된 가정적인 반사실적 \(A_k\) 는 시간 \(k\) 에서의 처치, \(L_k\) 는 공변량, \(\bar{A}_{k-1}\)은 이전 시간까지의 처치 역사를 나타내다. \(W_k\) 는 가중치로, 특정 모델에서 추정된 추가 공변량을 나타낼 수 있다.

이 모델은 처치가 주어진 경우와 주어지지 않은 경우 간의 로그값을 취한 오즈를 모델링한다. 각 \(a\) 는 로지스틱 회귀 모델에서 추정되어야 하는 파라미터이다. 이 모델은 \(\beta\) 의 점 추정치의 일관성에는 영향을 미치지 않지만, 신뢰 구간의 너비만을 결정한다.

• 모델이 올바르게 지정되고 가정이 유지된다면, 추정치는 샘플 크기가 증가함에 따라 실제 매개변수 값에 수렴하게 된다.

• 그러나 매개변수 추정치의 신뢰 구간의 폭은 공변량의 선택에 의해 영향을 받게 된다. 이는 모델에 각 공변량이 도입할 수 있는 변동성 때문이다.

• 모델 지정의 유연성(flexibility of model specification)을 가지고 있지만, 개념 자체 복잡하여 이해 및 적용이 쉽지 않다.

21.5 Censoring is a time-varying treatment

센서링은 시간에 따라 변하는 처치로 간주하고 인과 효과를 추정할 때 사용되는 중요한 공식은 IP 가중치를 사용하여 센서링으로 인한 선택 편향을 조정하는 과정에서 나타난다. 예를 들어, 센서링으로 인해 발생할 수 있는 선택 편향을 조정하기 위한 IP 가중치 \(W_C\)는 다음과 같이 계산된다:

\[ W_C = \prod^{K+1}_{k=1}\frac{1}{\Pr (C_k = 0|C_{k-1} = 0, A_{k-1, L_k-1})} \]

• 여기서 \(C_k\)​는 시간 \(k\)에서의 센서링상태를 나타내며, \(A_{k-1}\)과 \(L_{k-1}\)​은 이전 시간의 처치와 공변량을 나타낸다. 이 가중치를 사용하여, 센서링 없이 관찰되었을 모든 데이터에 대한 분석을 수행함으로써, 센서링으로 인한 선택 편향을 제거할 수 있게 된다다.

g-estimation 절차는 다음과 같은 단계를 포함한다:

1. 비안정화된 IP 가중치 \(W_{\bar{C}}\)를 추정하여 모든 참가자가 연구를 완료한 상태를 모사하는 가상의 모집단을 생성한다.

2. 생성된 가상의 모집단에 g-estimation을 적용하여 구조적 중첩 모델의 파라미터를 추정한다.

3. 추정된 모델 파라미터를 사용하여 처치의 인과 효과에 대한 추정치를 얻는다.

21.6 The big g-formula

관찰 연구에서는 인과관계 다이어그램이라고 할 수 있는 그래프가 실제로 인과관계 다이어그램인지 확실하게 알 수 없다.

g-formula는 모수적 모형 (parametric model)을 사용하기 때문에 다양한 고차원 항 또는 교호작용 항을 모형에 반영하는 것에는 한계가 있다.

또한 R 패키지인 gfoRmula는 생존 효과로 인한 편향을 제어하기 위하여 g-formula를 적용하는 것이 더 적절한 방법이나, 현재 g-formulaR package는 반복측정 되는 시간 간격이 일정하지 않은 경우에 대한 분석은 제공하고 있지 않다.