4 단계 g-formula 알고리즘과 IPW vs 표준화

1 정의

정의: Plug-in g-formula 추정량

표준화 식별식 \(\mathrm{E}[Y^a] = \sum_l \mathrm{E}[Y|A=a, L=l] \Pr(L=l)\) 의 우변을

1. \(\Pr(L=l)\) 를 표본 경험 분포로 대체.
2. \(\mathrm{E}[Y|A=a, L=l]\) 을 모수 모형 적합 후 예측치로 대체.

하여 추정한 양:

\[\widehat{\mathrm{E}}[Y^a]_\text{plug-in} = \frac{1}{n} \sum_{i=1}^{n} \widehat{\mathrm{E}}[Y | A=a, L_i]\]

각 환자 \(i\) 의 \(L_i\) 에 가상 처치 \(a\) 를 대입한 예측치를 평균.

직관 — Plug-in 의 의미: “식별식의 미지수를 데이터 추정값으로 끼워 넣는다(plug in)” 는 절차. 비모수 부분 (\(\Pr(L)\)) 은 표본 분포로 그대로 사용하고, 모수 부분 (\(\mathrm{E}[Y|A,L]\)) 만 모형으로 추정 — semi-parametric 추정량.

2 13.3 4 단계 표준화 알고리즘

2.1 알고리즘의 4 단계

데이터 확장 → 적합 → 예측 → 평균화

Step 1 — 데이터 확장: 원 데이터 \(n\) 행을 3 블록으로 복제 (\(3n\) 행). - 블록 1: 원 데이터 (관측된 \(A\), \(Y\)). - 블록 2: 모든 행에 \(A := 0\), \(Y := \text{NA}\). - 블록 3: 모든 행에 \(A := 1\), \(Y := \text{NA}\).

Step 2 — 결과 모형 적합: 블록 1 만으로 \(\mathrm{E}[Y|A,L]\) 회귀 (NA 제외).

Step 3 — 예측: 적합된 모형으로 블록 2·3 의 결측 \(Y\) 를 예측.

Step 4 — 평균화: 블록 2 의 예측 평균 = \(\widehat{\mathrm{E}}[Y^{a=0}]\), 블록 3 의 평균 = \(\widehat{\mathrm{E}}[Y^{a=1}]\).

차이 \(\widehat{\mathrm{E}}[Y^{a=1}] - \widehat{\mathrm{E}}[Y^{a=0}]\) 가 ATE 추정량.

직관 — 데이터 확장의 trick: 블록 2·3 은 “모두가 \(A=0\) 받았다면” 또는 “모두가 \(A=1\) 받았다면” 의 가상 시나리오. 결측 \(Y\) 라 모형 적합에는 영향 안 주지만 예측에는 사용된다. 반사실을 데이터 행으로 표현하는 코딩 트릭 — 직관적이고 자동화 쉬움.

직관 — 왜 같은 모형 한 번만 적합: 블록 1 의 적합은 “관측된 데이터에서 \(Y\) 가 \(A\) 와 \(L\) 에 어떻게 의존하는지” 를 학습. 그 학습된 함수를 블록 2·3 의 가상 데이터에 적용해 예측 — 학습된 관계가 가상 시나리오에도 그대로 통한다는 가정이 표준화의 핵심.

2.2 NHEFS 사례 결과

NHEFS 표준화 결과 (Hernan, Program 13.3)

• 블록 2 (모두 비금연) 의 평균 예측: \(\widehat{\mathrm{E}}[Y^{a=0, c=0}] = 1.66\) kg
• 블록 3 (모두 금연) 의 평균 예측: \(\widehat{\mathrm{E}}[Y^{a=1, c=0}] = 5.18\) kg
• ATE 추정값: \(5.18 - 1.66 = 3.52\) kg
• Bootstrap 95% CI: \((2.6, 4.5)\)

직관 — 두 평균의 의미: 1.66kg = “이 표본의 baseline 분포를 가진 사람들이 모두 비금연으로 살아간 가상 평행 우주의 평균 체중 변화”. 5.18kg = “같은 사람들이 모두 금연한 가상 우주의 평균”. 두 우주의 차이가 인과 효과 — 반사실의 정량화.

2.3 작은 데이터에서의 동등성

표 2.2 의 20 명 가상 데이터에서 4 단계 알고리즘과 직접 표준화 (Section 2.3) 의 결과가 정확히 일치 — 모두 비모수이므로 (saturated 모형 + 표본 분포).

직관 — 알고리즘이 이론과 같은 답을 내는 이유: 비모수 한도에서 식별식의 우변이 정확히 계산된다. 모수 모형을 도입해야 차이가 발생 — 그 차이는 모형 specification 의 흔적이다.

2.4 변수 종류별 처리

변수 종류	표준화 처리
이산 (성별, 인종)	더미 변수
다범주 (교육)	하나-핫 인코딩
연속 (나이, 체중)	직접 사용 + quadratic
연속 + 비선형	spline 또는 다항식
결측	imputation 또는 명시적 missing 카테고리

직관 — Continuous L 의 적분이 표본 평균이 되는 이유: 식별식의 적분 \(\int \mathrm{E}[Y|A=a,L=l] f_L(l) dl\) 는 표본의 경험 분포로 추정 시 단순 표본 평균 \(\frac{1}{n} \sum_i \widehat{\mathrm{E}}[Y|A=a, L_i]\) 가 된다. 적분이 합으로 변환되는 표준 estimator 의 형태.

3 13.4 IPW vs 표준화

3.1 비모수 동등성

Hernan, Technical Point 2.3 에서 증명:

\[\sum_l \mathrm{E}[Y|A=a, L=l] \Pr(L=l) = \mathrm{E}\left[\frac{\mathbb{1}\{A=a\} Y}{f(A|L)}\right]\]

좌변(표준화) 와 우변(IPW) 가 정확히 같은 양. 비모수 추정량으로 같은 값으로 수렴.

직관 — 같은 인과량의 두 표현: \(A\) 와 \(L\) 의 결합 분포를 두 가지 다른 분해로 표현. 표준화는 \(f(Y|A,L) f(L)\), IPW 는 \(f(Y, A | L) / f(A|L)\). 베이즈 정리의 두 측면. 어느 쪽으로 풀어도 같은 답.

3.2 모수 모형에서의 차이

측면	IPW	표준화
모형 대상	처치 \(\Pr(A|L)\)	결과 \(\mathrm{E}[Y|A,L]\)
Misspecification 영향	가중치 부정확 → 점추정 편향	외삽 부정확 → 점추정 편향
Positivity 위반	가중치 폭발	외삽 (가정 의존)
NHEFS 결과	3.4 kg	3.5 kg

직관 — 0.1kg 차이의 의미: 두 모형의 misspecification 이 약간 다른 방향으로 영향을 준 결과. 진짜 효과는 두 추정값 사이 어딘가, 또는 둘 다 같은 방향으로 약간 어긋났을 수도. 차이 자체는 작지만 측정 변동의 한계 안에 있다.

직관 — 큰 차이가 나면 무엇을 의심: NHEFS 의 0.1kg 차이는 안정적이지만, 다른 데이터에서 IPW 와 표준화가 1kg 이상 차이가 나면 (a) 처치 모형 미명세, (b) 결과 모형 미명세, (c) 양의 확률 근접 위반 중 하나. 진단 도구로 가중치 분포·결과 모형 잔차를 점검.

3.3 Doubly Robust 의 미리보기

Doubly Robust Plug-in 추정량 (Hernan, Fine Point 13.2)

처치 모형 \(\widehat{f}(A|L)\) 와 결과 모형 \(\widehat{m}(A,L)\) 모두 적합한 후, 둘 중 한쪽만 옳아도 일치 추정인 추정량.

대표 형태 — Bang & Robins (2005) 의 plug-in DR:

1. IPW 가중치 \(W^A = 1/\widehat{f}(A|L)\) 계산.
2. 결과 모형에 추가 변수 \(R\) 포함: \(R = W^A\) if \(A=1\), \(R = -W^A\) if \(A=0\).
3. 확장된 결과 모형으로 표준화 (4 단계 알고리즘).

직관 — 두 번의 기회: 처치 모형이 옳으면 IPW 보정이 작동. 결과 모형이 옳으면 표준화가 작동. 둘을 결합한 DR 은 어느 한쪽만 옳아도 답이 나오므로 모형 misspecification 위험이 효과적으로 줄어든다.

직관 — Augmented IPW (AIPW) 도 같은 아이디어: \(\widehat{\mathrm{E}}[Y^a]_\text{AIPW} = \widehat{m}(a, L) + \mathbb{1}\{A=a\} (Y - \widehat{m}(A,L)) / \widehat{f}(A|L)\) 의 표본 평균. 첫 항이 표준화, 둘째 항이 잔차에 대한 IPW 보정. 형태가 다르지만 같은 doubly robust 성질.

3.4 실무 권장 — 두 방법 모두 사용

Hernan 의 권장 (Ch.13.4)

• 둘 중 하나를 선택하지 말고 둘 다 사용.
• 결과가 비슷하면 specification 안정성의 증거.
• 결과가 다르면 모형 미명세의 경보 — 도메인 가설 재점검.
• 가능하면 doubly robust 도 함께 보고.

직관 — “한 도구만” 의 실무적 위험: 어떤 도구 한 개만 사용하면 그 도구의 약점이 결과의 약점이 된다. 두 도구를 함께 사용하면 각자의 약점이 보완 — IPW 의 양의 확률 위반 약점은 표준화로, 표준화의 외삽 위험은 IPW 로.

4 G-formula 의 일반화

정의: G-formula (일반)

조건부 교환가능성 + 양의 확률 + 일관성 아래

\[\mathrm{E}[Y^a] = \int \mathrm{E}[Y|A=a, L=l] f_L(l) \, dl\]

이 식이 시간변동 처치 (Part III) 로 확장되면

\[\mathrm{E}[Y^{\bar{a}}] = \int \cdots \int \mathrm{E}[Y|\bar{A}=\bar{a}, \bar{L}=\bar{l}] \prod_{t=0}^{T} f(L_t | \bar{L}_{t-1}, \bar{A}_{t-1}) \, dl_0 \cdots dl_T\]

여기서 \(\bar{a} = (a_0, a_1, \ldots, a_T)\) 가 처치 시퀀스. 이 일반 형태가 g-formula 라 불리며, Robins (1986) 의 핵심 기여.

직관 — “g” 의 의미: “general” 또는 “g-computation algorithm formula”. 단일 시점 표준화의 시간변동 일반화. Part III 의 모든 인과 분석의 기반 식이며, 시간변동 처치-교란 피드백 같은 어려운 문제도 같은 형태로 표현 가능.

5 외삽과 비외삽의 구분

추정 영역	데이터 존재	모형 의존도
비외삽 (interpolation)	데이터 점들 사이	낮음
외삽 (extrapolation)	데이터 범위 밖	매우 높음
약한 양의 확률	데이터 거의 없음	높음 (사실상 외삽)

직관 — 결과 보고 시 외삽 명시: NHEFS 사례에서 baseline 흡연량 25 개비/일 이하 부분 표본을 사용했다. 이 부분 표본 외부의 흡연량 변화 추정값은 외삽 — 결과 보고 시 “외삽 영역에서는 모형 가정에 의존” 명시 권장.

6 응용 분야

• 임상 코호트의 잠재 결과 분석: NHEFS 같은 사례의 표준 분석 패턴
• 공중보건 정책 평가: 정책 시나리오의 비교
• HEOR: 비용·QALY 의 잠재 결과 추정
• 마케팅 incremental analysis: 캠페인 적용 vs 미적용 시나리오
• 온라인 실험의 보정 분석: 무작위성 깨졌을 때

7 코드 — 4 단계 알고리즘 + 비교

import pandas as pd
import numpy as np
import statsmodels.api as sm
import statsmodels.formula.api as smf

nhefs = pd.read_csv("nhefs.csv").dropna(subset=["wt82_71"]).reset_index(drop=True)

formula = (
    "wt82_71 ~ qsmk + I(qsmk * smokeintensity) "
    "+ sex + race + C(education) + age + I(age**2) "
    "+ smokeintensity + I(smokeintensity**2) "
    "+ smokeyrs + I(smokeyrs**2) + C(exercise) + C(active) "
    "+ wt71 + I(wt71**2)"
)

# Step 1-2: 결과 모형 적합 (블록 1 만)
out_model = smf.ols(formula, data=nhefs).fit()

# Step 3: 가상 시나리오 데이터 (블록 2, 3) 예측
treated_block = nhefs.copy(); treated_block["qsmk"] = 1
untreated_block = nhefs.copy(); untreated_block["qsmk"] = 0

mean_y1 = out_model.predict(treated_block).mean()    # ~5.18
mean_y0 = out_model.predict(untreated_block).mean()  # ~1.66

# Step 4: 차이
ate_std = mean_y1 - mean_y0   # ~3.52
print(f"Standardization ATE: {ate_std:.2f} kg")

# Bootstrap CI
def bootstrap_std(data, n_boot=1000, seed=42):
    rng = np.random.default_rng(seed)
    estimates = []
    for _ in range(n_boot):
        idx = rng.integers(0, len(data), len(data))
        boot = data.iloc[idx].reset_index(drop=True)
        m = smf.ols(formula, data=boot).fit()
        t = boot.copy(); t["qsmk"] = 1
        u = boot.copy(); u["qsmk"] = 0
        estimates.append(m.predict(t).mean() - m.predict(u).mean())
    return np.array(estimates)

boot_estimates = bootstrap_std(nhefs, n_boot=1000)
print(f"95% CI: ({np.percentile(boot_estimates, 2.5):.2f}, "
      f"{np.percentile(boot_estimates, 97.5):.2f})")

# IPW 비교 — Ch.12 sequence
ps_model = smf.logit(
    formula.replace("wt82_71 ~", "qsmk ~ ").replace("qsmk + I(qsmk * smokeintensity) + ", ""),
    data=nhefs,
).fit(disp=False)
nhefs["ps"] = ps_model.predict()
nhefs["w"] = np.where(nhefs["qsmk"] == 1, 1/nhefs["ps"], 1/(1-nhefs["ps"]))

X = sm.add_constant(nhefs["qsmk"])
gee = sm.WLS(nhefs["wt82_71"], X, weights=nhefs["w"]).fit(cov_type="HC0")
ate_ipw = gee.params["qsmk"]
print(f"IPW ATE: {ate_ipw:.2f} kg")

# 두 추정값 비교
print(f"Standardization vs IPW: {ate_std:.2f} vs {ate_ipw:.2f}, diff = {abs(ate_std-ate_ipw):.2f}")

8 한 줄 요약

4 단계 plug-in g-formula 알고리즘은 데이터를 3 블록으로 확장 → 결과 모형 적합 → 가상 처치 시나리오의 예측치 평균 → 차이 계산. NHEFS 에서 5.18 - 1.66 = 3.5kg, IPW 의 3.4kg 와 거의 일치. 두 방법은 비모수 한도에서 동등하며 모수 차이는 specification 흔적. Doubly robust 추정량은 두 모형 중 한쪽만 옳아도 일치 — 모형 위험을 효과적으로 분산.

9 관련 주제

선행 지식

• Ch.13 표준화·g-formula 개관
• 표준화 절차 + 결과 모형 추정
• Ch.12 IPW와 MSM

후속 주제

• 추정값 신뢰성 — Ch.13.5
• G-estimation 과 SNMM — Ch.14
• 성향점수 — Ch.15
• Doubly Robust ML — Ch.18
• 시간변동 g-방법 — Ch.21

다른 카테고리 연결

• Bootstrap — 분산 추정의 비모수 도구
• 회귀 진단

10 G-formula 의 단일 시점 vs 시간변동

단일 시점 vs 시간변동 (Part III 미리보기)

단일 시점 (Ch.13): \[\mathrm{E}[Y^a] = \int \mathrm{E}[Y|A=a, L=l] f_L(l) dl\]

시간변동 (Ch.21): \[\mathrm{E}[Y^{\bar{a}}] = \int \cdots \int \mathrm{E}[Y|\bar{A}=\bar{a}, \bar{L}=\bar{l}] \prod_{t=0}^{T} f(L_t | \bar{L}_{t-1}, \bar{A}_{t-1}) dl_0 \cdots dl_T\]

시간변동은 매 시점의 \(L\) 분포가 과거 \(A, L\) 에 의존 — 식이 길어짐.

직관 — 시간 차원의 폭발: 시점 1 개에서 적분 1 개. 시점 10 개면 적분 10 개. 각 적분이 표본 평균으로 추정 가능하지만 분산이 시점 곱으로 누적. 단일 시점 g-formula 의 깔끔함이 시간변동에서는 복잡성과 분산 폭발의 위험으로 변환.

직관 — Monte Carlo 시뮬레이션 g-formula: 시간변동 g-formula 의 다중 적분은 해석적으로 풀기 어려워, 시뮬레이션으로 가상 코호트의 결과를 생성한다 — 각 시점에 처치를 강제 적용하고 \(L\) 의 진화를 표본 분포에서 추출. Robins (1986) 의 핵심 알고리즘.

11 Doubly Robust 의 두 형태 비교

DR plug-in vs AIPW

DR plug-in (Bang & Robins 2005): 결과 모형에 IP 가중치 \(W^A\) 를 추가 변수로 포함 후 표준화.

AIPW (Robins, Rotnitzky, Zhao 1994): \[\widehat{\mathrm{E}}[Y^a]_\text{AIPW} = \frac{1}{n} \sum_i \left[ \widehat{m}(a, L_i) + \frac{\mathbb{1}\{A_i=a\}}{\widehat{f}(A_i|L_i)} (Y_i - \widehat{m}(A_i, L_i)) \right]\]

두 형태가 점근적으로 동등하지만 유한 표본에서 미세하게 다른 행동.

직관 — AIPW 의 두 항: 첫 항이 표준화 (모형 옳음 시 일치), 둘째 항이 잔차에 대한 IPW 보정 (처치 모형 옳음 시 잔차의 평균이 0). 결과 모형 옳음 + 처치 모형 틀림 → 첫 항이 답, 둘째 항의 평균이 0 으로 수렴. 처치 모형 옳음 + 결과 모형 틀림 → 첫 항이 편향, 둘째 항이 그 편향을 정확히 상쇄. 두 시나리오 모두에서 답이 진짜 효과로 수렴하는 마법.

직관 — DR 의 정규화 편향 보호: 결과 모형이 ML (lasso, random forest 등) 로 적합되면 정규화 편향이 발생. AIPW 의 둘째 항이 이 편향을 첫 항에 끼치는 것을 자동 보정. Ch.18 의 Targeted Maximum Likelihood Estimation (TMLE) 가 이 원리의 정점.